Elsevier

Diario de Hidrología

Volumen 600, Septiembre de 2021, 126508
Journal of Hydrology

Documentos de investigación
Desafía de modelización de Karst 1: Resultados de la modelización hidrológica

https://doi.org/10.1016/j.jhydrol.2021.126508Obtenga derechos y contenido
Bajo licencia Creative Commonslicense
acceso abierto

Aspectos destacados

    -
  • -

    , los modelos de depósito, se comparaban directamente los modelos semi- y completamente distribuidos.

  • -

    Cómo se pudieron evaluar los resultados? Los criterios de evaluación de Kling-Gupta fueron reconocidos como los mejores criterios, aunque no perfectos.

  • -

    Los pasos de tiempo por hora tanto para los datos como para la modelización son más adecuados que los diarios para la mayoría de los sistemas kárstico.

  • -

    La mayoría de los modelos son razonablemente buenos, pero mal predecir los bajos caudales de agua.

Resumen

La complejidad de la modelización del flujo de aguas subterráneas kársticas se refleja en la cantidad de enfoques de simulación. El objetivo del Karst Modelling Challenge (KMC) es comparar diferentes enfoques en un solo sistema utilizando el mismo conjunto de datos. Trece equipos con diferentes modelos computacionales para simular variaciones de descarga en los muelles kársmicos han aplicado sus respectivos modelos en un solo conjunto de datos procedentes del Sistema Hidrogeológico Milandre Karst (MKHS). Los enfoques incluyen redes neuronales, modelos de embalse, modelos semidistribuidos y modelos de aguas subterráneas totalmente distribuidos. Se proporcionaron cuatro años y medio de entrada meteorológica por hora o diario y datos de descarga por hora para la calibración modelo. La validación comunífice de un año de alta, sin los datos de descarga observados. El rendimiento del modelo se evaluó utilizando la conservación del volumen, la eficiencia de Nash-Sutcliffe (NSE) y la eficiencia de Kling-Gupta (KGE) aplicadas en el vertido total y componentes de flujo individuales. Como resultado, la comparación de los rendimientos de los modelos es una tarea difícil debido a las diferencias en la arquitectura del modelo, pero también requiere pasos de tiempo: algunos de los modelos requieren pasos diarios agregados mientras que otros podrían ejecutarse utilizando datos por hora, lo que proporcionó algunas diferencias interesantes dependiendo de cómo se transformaron los datos. El uso de datos instantáneos (por ejemplo, valor al mediodía) produce menos sesgo que promediando los datos por hora a lo largo de un día. La transformación de los datos por hora en diarios produce una disminución de Nash y KGE de 0,05 a 0,05 a 0,058 (es decir, de 1 a 0,73). Las simulaciones resultantes (valores predecidos para el año 2016) produjeron KGE que oscilaban entre 083 y 0,37 (0,83 a 0,24 euros para NSE). Aunque las simulaciones coinían con los flujos monitoreados razonablemente bien, la mayoría de los modelos lucharon para simular las condiciones del flujo de base con precisión. En general, los modelos que mejor realizaron para este ejercicio fueron los globales (Jardeia y Varkarst), con un número limitado de parámetros, que se pueden calibrar mediante procedimientos de calibración automática. Los modelos de red neuronal también mostraron un potencial justo, con uno proporcionando resultados razonables a pesar del conjunto de datos relativamente corto disponible para el calentamiento (4,5 años). Los modelos semi-y completamente distribuidos también sugirieron que con algún esfuerzo más podrían funcionar bien. La precisión de las predicciones de los modelos no parece aumentar mediante el uso de modelos con más de 9-12 parámetros de calibración. Una evaluación de los errores relativos entre los valores previstos y los observados reveló que para la mayoría de los modelos, el 50% de los valores previstos contenían más del 50% de la diferencia con la tasa de descarga observada, con un 25% con una diferencia superior al 100%. Una parte significativa de los valores mal previstados correspondió al flujo de base, lo que fue sorprendente dado que como flujo de base es generalmente mucho más fácil de predecir que el flujo máximo. Por lo tanto, esto demuestra que los enfoques de modelización y los criterios para la calibración están demasiado orientados a las secciones de flujo máximo de los hidrógrafos, y que las mejoras podrían obtenerse más centradas en el flujo de base.

Palabras clave

Modelización
Karst
Comparación
Criterios de eficiencia
Recargo
Flujo de base

1. Introducción

1.1. Idea

El objetivo del reto de la modelización era invitar a equipos de investigación de todo el mundo a comparar sus enfoques y herramientas de modelización aplicándolos en el mismo conjunto de datos. Había que definir un método de evaluación común y los resultados de cada equipo han sido analizados con los mismos criterios. De esta manera todos los resultados podrían compararse directamente. Se han probado varias funciones objetivas (varias formas de calcular la diferencia entre las series temporales previstas y las series temporales observadas), con tres posteriormente utilizadas para la comparación final.

1.2. Definición de los problemas

Para comenzar con una pregunta bastante simple (Pase 1) se decidió centrarse en el comportamiento hidrológico de los sistemas hidrogeológicos kárstico (KHS), es decir, la relación entre los parámetros que controlan la entrada de agua en el kárstico (principalmente precipitación y temperatura) y la tasa de descarga de los muelles kárstico a la salida del sistema (Fig. 1). La cuestión que deben resolver los respectivos equipos es pronosticar con la mayor precisión posible las tasas de alta en el jarlet kárstico de los datos de entrada meteorológica.

  1. Descargar Descargar imagen de alta resolución (346KB)
  2. Descargar Descargar imagen a tamaño completo

Fig. 1. Modelo conceptual de un sistema hidrológico kárstico (izgato). El agua de precipitación (alógena y autógena) fluye a través del macizo kárstico a lo largo de conductos (flejas rojas) y fisuras (flejas azules) con el pico de descarga correspondiente en la primavera después de unas pocas horas o días. Para KMC step 1 sólo consideramos la relación entre precipitación y descarga. Los subsistemas pueden o no ser considerados en los enfoques de modelización aplicados. (Para la interpretación de las referencias al color en esta leyenda figura, el lector es referido a la versión web de este artículo.)

Cabe señalar que se prevén nuevas medidas, con la simulación del flujo y la distribución de la cabeza en el espacio y el tiempo (Pase 2), y la simulación de los procesos de transporte (Pase 3). Sin embargo, la evaluación de este paso 1 del reto no tiene en cuenta si los modelos podrán abordar o no los objetivos de los nuevos pasos.

En este trabajo, la palabra "Modelo" se refiere a los enfoques y las herramientas numéricas correspondientes utilizadas por los equipos respectivos para su ejercicio de modelado. La simulación de la palabra se utiliza para nombrar los resultados del trabajo de modelado.

1.3. Fluen a través de un macizo kárstico

Como se describe en la mayoría de los libros de texto (por ejemplo, Ford y Williams, 1989), los macizos kársticos se caracterizan por una geomorfología e hidrología específicas (Fig. 1). Las formas de tierra kratra se deben a la disolución de la roca en el agua de precipitación. Las formas terrestres de superficie como karrenfields, dolines y cuevas se producen por infiltración de en la roca. Los ejes verticales, las cuevas (vainas o fratic) y los grandes muelles están relacionados con el flujo y la salida de agua.

Desde un punto de vista hidrológico la precipitación (llugar y derretimiento de la nieve) se infiltra a través de la cubierta vegetal y los suelos en la capa superior de la roca, que es más o menos climatizado (epikarst). Una parte de las aguas de precipitación regresa de vuelta a la atmósfera a través de evapotranspiración. En algunos casos, las aguas superficiales pueden formar un pequeño arroyo o incluso un río, que se traga en una ubicación concentrada (aguño de traga). Si el arroyo drena una zona no kárstica esta parte de la cuenca se llamará "allogenic".

Epikarst distribuye el agua parcialmente hacia ejes verticales que conducen el líquido directamente a la zona freática. Otra parte del agua está atrapada en el fondo del epikarst y se filtra lentamente hacia abajo a través de grietas y articulaciones. En algunos puntos el agua es drenada de nuevo por la red de canales kársticos, que drena eficientemente la masa de roca, y conduce el agua en la salida del sistema kárstico: el manantial kárstico (Mangin, 1975, Williams, 1983, Smart y Friederich, 1986). El volumen de agua almacenada en fisuras es significativo en comparación con el almacenado en los canales kársticos (Király, 1975, Drogue, 1980, Smart y Hobbs, 1986).

La relación entre la recarga y la descarga en los muelles parece bastante directa en la mayoría de los sistemas kárstico, ya que los eventos de lluvia son seguidos por picos claros de las tasas de descarga en los muelles kárstico generalmente en unas pocas horas a días. Por lo tanto, la respuesta es bastante similar a la que se suele observar en la hidrología superficial.

El grado de flujo difuso en la matriz rocosa en comparación con el flujo concentrado en conductos, así como la cantidad de recarga concentrada en agujeros de golondrina vs difamar difuso a través de suelos y epikarst puede ser significativamente diferente dependiendo de las características de la roca y de contextos geomorfológicos y climatológicos.

Se han desarrollado varios modelos para simular la relación funcional entre los picos de precipitación y . Algunos de ellos incluyen procesos físicos que tienen lugar en cada subsistema atravesado por el agua, otros simplifican la naturaleza en subsistemas lineales, y otros sólo consideran la relación matemática entre una señal de entrada y una salida. El objetivo del paso 1 del Karst Modelling Challenge es comparar estos enfoques, y su capacidad para simular un hidrografo de una serie de tiempos de precipitación.

1.4. Perspectivas cortas sobre el modelado del flujo de aguas subterráneas en kárstico

El análisis de los hidrografos comenzó principalmente después de la publicación de Maillet (1905). Sin embargo, la monitoreo continuo de los muelles kársticos realmente comenzó sólo desde la década de 1950. Para el análisis de los manantiales kársticos, los primeros modelos ocurrieron a lo largo de la década de 1960 (Schoeller, 1962, Forkasiewicz y Paloc, 1967, Castany, 1968, Brown, 1970, Brown, 1973, Mangin, 1970, Mangin, 1975, Chemin, 1974, Bezes, 1976, Milanovic, 1976, Atkinson, 1977, Dreiss, 1982, Dreiss, 1983, Dodge, 1983, Bonacci, 1987, Ford y Williams, 1989). Al principio se reconocieron dos componentes, uno explicando la rápida e intensa respuesta de las tasas de descarga de primavera a las precipitaciones, y el otro explicando la lenta disminución de las tasas de descarga (recesión) durante largos períodos de tiempo sin eventos de recarga.

Varias ideas conceptuales se desarrollaron a través de varios tipos de modelos matemáticos, que pueden agruparse en tres familias numerosas (Teutsch y Sauter, 1991, Kovács y Sauter, 2007, Ghasemizadeh, 2012, Fiorillo, 2014):

  • 1)

    Modelos basados en datos o de caja negra, incluidas funciones de transferencia de reservorio y no paramétricas

  • 2)

    Modelos hidrogeológicos (flujo de matriz de cierre con o sin conductos)

  • 3)

    Modelos de flujo de tuberías (basado en la hidráulica de flujo en conductos)

1.4.1. Modelos de caja de datos o de caja negra

En estos modelos, el sistema sólo se considera como una caja negra o una combinación de cajas negras, cada caja transformando una señal de entrada (influjo) en una salida (salida).

La transformación se basa en una función puramente matemática, sin considerar los procesos físicos que controlan el flujo en la clandestinidad o simplificándolos fuertemente (Dewandel et al., 2003, Fiorillo, 2014). El modelo más sencillo de este tipo sería un modelo de una sola caja, tomando la precipitación total como entrada y transformándolo en una serie de tiempo de descarga. El modelo clásico de hidrografía .unirtario, suponiendo una sola función de transferencia lineal (núcleo), introducido por primera vez por Sherman (1932), fue probado en sistemas kárstico. Rápidamente resultó que KHS no es lineal ni estable y que un modelo tan simple no puede describir los acuíferos kárstico de una manera significativa. Así, los esquemas más sofisticados, acoplando varios núcleos, eran necesarios para adaptarse a los datos observados de alguna manera (Dreiss, 1982). Las redes neuronales son los desarrollos más recientes de esta categoría de modelos.

Otra categoría de modelos agrupados fue intentado por varios autores, utilizando depósitos paralelos en cascada y paralelos (Forkasiewicz y Paloc, 1967, Bezes, 1976, Halihan y Wicks, 1998). Por lo general, un reservorio simula el componente de flujo rápido, y los demás el componente más lento de KHS. A continuación se produjo un nuevo parámetro, es decir, la distribución de la recarga en los componentes respectivos. Los mejores resultados se obtuvieron con una combinación de tres a cinco embalses, algunos de ellos parcialmente vaciados por evapotranspiración. Este enfoque condujo a resultados de simulación bastante precisos (por ejemplo, Largo, 2009). En la mayoría de los casos, los parámetros de entrada no se distribuyen, es decir, cada embalse se supone que se expande a través de todo el tamaño de la cuenca. De hecho, este enfoque también puede distribuirse en el espacio, y puede tener en cuenta algunas diferencias espaciales de las características en las partes respectivas de la zona de captación (por ejemplo. Bittner et al., 2018). Con tres a cinco embalses distribuidos en decenas o cientos de zonas el número de parámetros se vuelve muy alto, lo que dificulta la calibración de los modelos. Gracias a los procedimientos de calibración automática, esta tarea se puede lograr ahora con un esfuerzo razonable (Pianosi et al., 2016, Mazzilli et al., 2017).

1.4.2. Modelos hidrogeológicos (modelos distribuidos)

Darcy (1856) describió empíricamente la relación matemática entre los gradientes de agua subterránea (diferencias de cabeza) y la cantidad de flujo de aguas subterráneas. Este fue de hecho el primer modelo de flujo de agua subterránea, que era unidimensional y de estado estacionario. Se desarrollaron soluciones analíticas para algunas geometrías y condiciones de límites típicos. En la década de 1970 el desarrollo de computadoras permitió resolver las ecuaciones numéricamente, haciendo posible calcular el flujo a través de una gama mucho más amplia de situaciones (por ejemplo. Bredehoeft y Pinder, 1970). La aplicación de este enfoque del kárstico se intentó por primera vez a mediados de la década de 1970 (Király y Morel, 1976), pero rápidamente se hizo evidente que KHS no puede ser fácilmente modelado usando este enfoque, ya que la simulación de eventos rápidos e intensos de inundación, así como un flujo de base sostenido, como se observó en la mayoría de los resortes karst resultó difícil. Los modelistas decidieron introducir una red de conductos dentro de una matriz con una conductividad hidráulica baja. Con contrastes de 10 -4 a 10 -5 m/s entre matriz y conductos, los hidrografías kársticos podrían ser simulados. Sin embargo, el proceso de recarga también tuvo que mejorarse, lo que se logró añadiendo una capa con una alta conductividad hidráulica en la parte superior del modelo (epikarst), capaz de absorber la precipitación y distribuirla a los conductos y la matriz (Klyirá et al., 1995). Dado que el alto contraste de las conductividades hidráulicas en los modelos espaciales puede ser bastante desafiante desde el punto de vista numérico, se han desarrollado varias soluciones, como doble continuidad o modelos medianos dobles (Teutsch y Sauter, 1991).

Dos mejoras más tarde ocurrieron: flujo turbulento en conduits y flujo parcialmente saturado (Thierrien y Sudicky, 1996, Annable y Sudicky, 1998).

En estos modelos completamente distribuidos el número de parámetros es muy alto (con varios valores para cada célula) haciendo manual de calibración muy lento y prácticamente imposible. La zona de parámetros y los procedimientos de calibración automática mejoraron este aspecto (Doherty et al., 1994, Borghi et al., 2016). Sin embargo, dado que primero hay que definir la estructura básica de la red, su geometría y topología no suele ajustarse durante el proceso de calibración, aunque puede desempeñar un papel significativo en algunas situaciones (Kovács, 2003) especialmente en la zona epifreática (Jeannin, 2001). Suponiendo que la física considerada en el modelo sea correcta, los modelos pueden revelar sin embargo estructuras poco realistas y, por lo tanto, pueden actuar como una herramienta muy avanzada que proporciona más información (Enemark et al., 2019, Gill et al., 2020, Duran y Gill, 2021).

Estos modelos pertenecen a modelos de flujo de agua subterráneos basados en procesos de acuerdo con la clasificación propuesta por Anderson et al. (2015).

1.4.3. Modelos de flujo de tuberías (modelos seemidistribuidos)

Los modelos de flujo de tuberías sólo consideran una red de conductos y flujo de negligencia a través de la matriz, o agregarlo como una entidad no espacializada. La idea no es nueva, ya que muchos autores originales (Martel, 1921, Livre, 1915, Livre, 1940, Trombe, 1948), explorando el kárstico con una perspectiva hidrogeológica y espeleológica, sugirió que el flujo en kárstico fuera bien descrito por las ecuaciones de la tubería hidráulica. Estos modelos también pertenecen a modelos basados en procesos según Anderson et al., (2015), pero el proceso principal considerado es la hidráulica de flujo a través de conductos en las zonas fréticas y epihreáticas (turbulenta y parcialmente saturada). El flujo de Darcys no es o sólo marginalmente considerado.

Unos pocos autores utilizaron entonces modelos de tuberías para analizar sus datos (blanco y blanco, 1970, Atkinson, 1977, Boegli, 1980, Lauritzen et al., 1985, Smart, 1988). Todos ellos demostraron que el flujo es turbulento en la mayoría de los conductos naturales (cuevas), al menos durante las condiciones medias y altas del agua. El flujo puede ser transitorio o incluso laminar durante condiciones de caudal muy bajas. Ford y Williams (1989) y Worthington (1991)

Jeannin y Maréchal (1995) y Jeannin (2001) evaluaron las cabezas de cabeza debido a las curvas y los cambios en las secciones transversales en comparación con los relacionados con la fricción a lo largo de las paredes. Resulta que estos últimos suelen estar dominando, y los valores típicos podrían ser asignados para conductos naturales. También calcularon la hidráulica de una red de conductos naturales para el flujo de modelado en la parte aguas abajo de la cueva de Hálloch (Suiza). El papel de los conductos epifráticos en la hidráulica y la velocidad del flujo se encontró que era significativo para muchos KHS. Este comportamiento específico no se puede simular con los paquetes de software hidrogeológicos habituales. Por esta razón, los modelos de flujo de tuberías (por ejemplo, el software SWMM desarrollado por US EPA y otros programas de drenaje urbano) se han utilizado cada vez más en el modelado de hidrología kárstica: Campbell y Sullivan, 2002, Peterson y Wicks, 2006, Wu et al., 2008, Gill et al., 2013, Chen y Goldscheider, 2014, Jeannin et al., 2015, Kaufmann et al., 2016, Malard, 2018, Vuilleumier et al., 2019, Morrissey et al., 2020, ,2020. Estos modelos permiten simular el flujo en una compleja red de tuberías con diversos diámetros, rudicidad y formas transversales.

Si este enfoque es eficiente para vincular las cabezas observadas y las tasas de descarga, no es adecuado vincular la precipitación para descargar hidrografía. De hecho, siempre es necesario hacer un modelo de recarga para hacerlo, con el fin de transformar los datos meteorológicos en recarga de aguas subterráneas.

1.4.4 Recarga de modelos de sistemas hidrogeológicos kárstico (KHS)

La recarga se define aquí como la cantidad de agua infiltrando y fluyendo a través del sistema kárstico. Corresponde aproximadamente a la precipitación total menos evapotranspiración y la escorridad terrestre eludiendo el kárstico en algunas circunstancias. También se considerará el almacenamiento de nieve y hielo en los procesos de recarga.

La recarga depende en gran medida de la evapotranspiración. Se han desarrollado varias fórmulas para evaluar la Evapotranspiración Potencial Thornthwaite, 1948, Penman, 1948, Turc, 1961). Estos modelos proporcionan una evaluación de la cantidad de agua que se evaporaría o transpiraría por vegetación si hubiera agua suficiente disponible en los suelos.

En realidad, durante períodos de bajo agua, los suelos se vuelven fuertemente insaturados y las plantas pueden sufrir un déficit hídrico, reduciendo su capacidad de traer agua del suelo a la atmósfera. La Evapotranspiración Real (RET) es por lo tanto menor que la PET. En los climas húmedos la diferencia puede ser insignificante, pero se vuelve significativa o incluso extrema en los climas más secos.

Muchos modelos se pueden utilizar para la evaluación de PET y RET, desde muy simples a muy sofisticados. Las mediciones reales de RET son difíciles, especialmente a la escala de captación.

1.5. Objetivos del desafío de la modelización

El objetivo general del desafío de modelización, del que presenta este trabajo, es explorar modelos para su capacidad de reproducir los siguientes aspectos: recarga subterránea, velocidad de flujo de aguas subterráneas en el espacio y el tiempo, cabezales de agua subterránea en el espacio y el tiempo, hidrógrafos de descarga de primavera.

El primer paso del desafío se limita a la relación entre los parámetros meteorológicos (parámetros de recarga de las aguas subterráneas) y los hidrógrafos de primavera (tasa de descarga a la producción del sistema), es decir, se centra en aspectos de la recarga de las aguas subterráneas y los hidrógrafos de primavera.

Cualquier grupo o persona interesada en comparar los resultados de su modelo con los demás fue invitada abiertamente a participar 11.

2. Enfoques de modelado en comparación en KMC

Todos los tipos de modelos descritos anteriormente se aplicaron para el desafío: Las dos categorías de modelos de parámetros agrupados (redes neuronales y depósitos), así como modelos semidistribuidos y completamente distribuidos.

2.1. 2.1. Modelos a base de datos 1: Redes neuronales

Dos equipos de investigación aplicaron redes neuronales para la simulación de hidrográficas kársticas (KIT-Karlsruhe con tres modelos y IMT Mines Als con uno). Este último equipo utiliza Perceptor MLP) el antiguo equipo aplica el modelo NARX, Convolutional Neural Networks (CNN), así como Long Short-Termry Memory Networks (LSTM).

Las redes neuronales artificiales (ANN) forman una rama de inteligencia artificial. Se basan en operaciones de una sola neurona, que se arreglan y conectan en una arquitectura específica utilizando parámetros (también llamados pesos). El comportamiento del modelo está programado por los valores asignados a los parámetros establecidos. Durante la etapa de entrenamiento, los datos se presentan a la red y los algoritmos de entrenamiento especiales se ajustan a los datos objetivo modificando los parámetros entre las neuronas. Una neurona calcula dos valores: (i) la suma ponderada de su vector de entrada con sus parámetros, y (ii) su salida, que transforma la suma ponderada en un valor escalar utilizando funciones predefinidas (por ejemplo, lineales o sigmoid), conocidas como funciones de activación. La respuesta no lineal de las funciones de activación sigmoide de la calificación permite a la ANN capturar relaciones no lineales dentro de los datos de entrenamiento.

En esta etapa, es importante notar que los modelos ANN no tienen función predefinida antes de entrenar. El objetivo del entrenamiento es construir tanto la función como para calcular sus parámetros. Para ello, es necesaria una base de datos que represente todo el tipo de conductas. Por lo general, para representar el comportamiento con buena calidad, la longitud necesaria de la base de datos depende de los comportamientos a considerar. Cuando la base de datos es continua con un pequeño paso de tiempo (1 h), como en este estudio, se llevan a cabo varios comportamientos en diferentes escalas de tiempo; por lo tanto, la base de datos debe ser lo suficientemente larga para permitir que el modelo capte los fenómenos subyacentes: 15 años deben ser necesarios (Artigue et al., 2012, Coutouis et al., 2016). En este estudio sólo se disponía de 4,5 años en dos períodos diferentes de tiempo. Por lo tanto, es importante notar que la base de datos utilizada en esta labor no es suficiente para los enfoques de ANN. Los resultados mostrados a lo largo de la siguiente A continuación pueden considerarse como una prueba de viabilidad.

A continuación se resumen los diversos enfoques (modelos) utilizados en este estudio. Como estas técnicas son diversas y no se utilizan comúnmente en la modelización de kárstico se dan más detalles como material suplementario.

Todos los modelos de red neuronal diseñados en este trabajo utilizaron precipitación y temperatura por hora, datos diarios de evapotranspiración remodelados a paso hora por hora, y descarga de datos a un paso de hora. Los rangos de datos utilizados para la formación o validación difieren ligeramente dependiendo del modelo.

2.1.1. Perceptor multicapa profunda recurrente (ANN/rec-MLP)

El perceptor multicapa elegido es una red neuronal recurrente (RNN) con una capa oculta de hneuronas ocultas y una neurona de salida. Las neuronas de capa oculta no son lineales y aplican una función sigmoide. La neurona de salida es lineal; su salida es igual a la suma ponderada de sus entradas. Esta arquitectura específica muestra las dos propiedades de aproximación universal y parsimonia, esta última debido a la nolinealidad relativa a los insumos y parámetros del modelo (Barrón, 1993). La aproximación Universal (Hornik et al., 1989) significa que el modelo es capaz de identificar cualquier función diferente, siempre que la existencia de una base de datos que represente el comportamiento para identificar. Para más detalles sobre esta arquitectura, el lector interesado se refiere a Dreyfus (2005). En este estudio se utiliza un perceptor multicapa profundo añadiendo una capa profunda, conectada para implementar un preesqueque de datos. La neurona de salida es sigmoide. El modelo contiene 381 parámetros.

2.1.2. Modelo NARX

El modelo NARX (red autorista lineal con entradas exógenas) utilizada por KIT (Karlsruhe) es bastante similar al modelo anterior utilizado por la escuela minera en Als (Francia). Las diferencias son, por ejemplo, el número de capas ocultas, los datos de entrada que se utilizan y la arquitectura.

Los modelos NARX también pueden pertenecer al grupo de redes neuronales recurrentes (RNN). NARX en particular tiene una conexión de retroalimentación global entre su salida y capa de entrada, permitiendo el flujo de información en diferentes direcciones de la red, y haciéndolos especialmente adecuados para modelar sistemas dinámicos no lineales. Esto también significa que cada valor de salida no sólo se retrocede en las señales de entrada independientes, sino también en las señales de salida anteriores. Dependiendo del modelo elegido, los RNN pueden tener dificultades para capturar dependencias a largo plazo mayores a 10 pasos de tiempo debido al problema de la desaparición de los gradientes (Bengio et al., 1994); sin embargo, NARX puede mantener la información hasta tres veces más tiempo que los RNN simples (Lin et al., 1996, Lin et al., 1996b). Para construir y aplicar los modelos NARX, utilizamos Matlab 2019a (Mathworks Inc.) y su Deep Learning Toolbox. El modelo contiene 371 parámetros.

2.1.3. Red de memoria a corto plazo Largo

Enfoques de aprendizaje profundo como Long Short-Term Memory Networks (LSTM) recibieron mucha atención en los últimos años, debido a éxitos significativos en varias disciplinas. Los LSTM también pertenecen al grupo de RNNs y se aplican ampliamente al modelo de datos secuenciales como series de tiempo o lenguaje natural. A diferencia de otros RNN, los LSTM han sido diseñados explícitamente para superar el problema de los gradientes desaparecidos. Además del estado oculto inherente a todo tipo de RNNs, los LSTM tienen una memoria celular (o estado celular) para almacenar información y tres puertas para controlar el flujo de información. Esto permite que la información permanezca en la memoria celular, razón por la cual los LSTM pueden manejar dependencias a largo plazo (Hochreiter y Schmidhuber, 1997). El modelo LSTM que simula la salida de la cueva de Milandre tiene 11.111 parámetros entrenables.

2.1. 2.1.4. Redes neuronales convolucionales

Convolutional Neural Networks (CNN) son un tipo de red neuronal de aprendizaje profundo multicapa que fue diseñada originalmente para manejar eficientemente los datos de imagen, pero que también se ha aplicado con éxito a la predicción de series temporales, tratando una secuencia de observaciones como una imagen unidimensional. La CNN que simula la salida de la cueva de Milandre tiene 950.921 parámetros entrenables.

2.2. Modelos basados en datos 2: embalse

Los modelos Gardenia (BRGM), RCD-Seasonal (TU-Freiberg), CHLEM (Uni-Zuerich) y KarstMod (SNO KARST) son herramientas muy similares diseñadas para simular los principales procesos del equilibrio hídrico a escala de cuenca utilizando leyes físicas simplificadas que representan flujos a través de una sucesión de embalses. Utilizando series temporales climáticas (precipitación, potencial evapotranspiración, temperatura del aire) en una zona de recarga, estos modelos son capaces de calcular el flujo en la salida de un KHS (primación).

2.2.1. CCD-seasonal

El sistema kárstico está representado por tres almacenes, que están conectados por enlaces: un almacenamiento para representar el sistema de carga de carga (suelo y epikarst), y otros dos para el cachorro c y el .

2.2.2. KarstMod

KarstMod es una plataforma modular para modelar la relación a nivel de lluvia en las cuencas kársticas desarrolladas por el francés SNO Karst. En su forma más completa (Mazzilli et al., 2019) ofrece 4 compartimentos organizados como una estructura de dos niveles. En este estudio, se utilizó para probar el rendimiento de una función de transferencia de dos parámetros con un tiempo característico infinito (Guinot et al., 2015). El modelo propuesto tiene una estructura de dos niveles. Un compartimento superior está destinado a explicar el flujo dentro del suelo y epikarst. Sobre un cierto nivel de agua se desborda a un compartimento inferior, lo que significa la infiltración y las zonas saturadas. La respuesta de unidad del compartimiento inferior es una función tiempo, que permite tiempos característicos infinitos. Tal respuesta puede ser muy adecuada para reproducir los efectos de memoria a largo plazo de KHS.

2.2.3 Gardenia

El código numérico GARDENIA (Thiéry, 2014, Thiéry, 2015, http://gardenia.brgm.fr/) tiene 4 embalses: i) un reservorio superficial S que representa la capacidad de retención de agua del suelo (ley cuadrática), ii) un embalse lineal intermedio que representa aproximadamente la zona insaturada, y dos depósitos lineales subterráneos G1 y G2 que representan compartimentos de la zona saturada. La descarga total es la suma de los flujos rápidos, lentos y más profundos, que se enruta a la salida.

2.2.4. CHLEM (Uni-Zuerich)

Modelo de Elemento Lumped de la Cueva (CHLEM) proporciona diferentes tipos de elementos hidráulicos que se pueden unir libremente. Para KMC, el sistema hidráulico se simuló con los conductos paralelos que producen la descarga del modelo, que se compara con las mediciones. Cada una de las cuatro tuberías consiste en un elemento de almacenamiento y un elemento de retardo y recibe agua de un elemento de entrada con o sin evapotranspiración. En los elementos de almacenamiento, la descarga depende del volumen de almacenamiento después de una ley de alimentación, y está limitada por un parámetro ajustable. Cada uno de los elementos tiene de uno a tres parámetros ajustables que fueron ajustados con un algoritmo de optimización (COBYLA de la biblioteca de optimización de nlopt).

2.3. Modelos semidistribuidos

KRM_1 (SISKA), Varkarst (Uni-Freiburg), and InfoWorks (TCD Dublin) models are all semi-distributed models but with significant differences. The number of parameters for these models is moderate.

KRM_1 is very similar to lumped parameter models with reservoirs but input parameters can be semi or fully distributed. It has been developed as an alternative to models (daily rainfall-runoff models based on two reservoirs and three parameters (Edijatno and Michel, 1989), in order to further refine interception and soil infiltration processes which are assumed to be the main recharge controls for lowland and vegetated karst aquifers. Storage capacities of interception reservoirs vary according to the type and respective proportion of land-uses (forests, cultures, urban areas, denudated areas, etc.) and according to seasons in order to reproduce the dynamic of the vegetation. Interception and ETP parameters are not uniform over the whole catchment. Two slow reservoirs are introduced to mimic processes in the soil/epikarst and in the deep vadose zone. They allow underflow and overflow depending on recharge intensity. Parameters of soil/epikarst reservoirs vary according to a drought index while those of the deeper reservoirs are unique. In addition, an “exchanger” makes the recharged water pass from the low permeability volumes (LPV) to the karst conduits and vice-et-versa depending on pseudo hydraulic relations. The complete description of KRM_1 is available in Malard [2018]. The spatial distribution of rainfall over the catchment area can be taken into account, but this was not applied in the present application (rain is spatially uniform, but vegetation cover is semi-distributed). The number of parameters for the calibration is 18.

Varkarst (Hartmann et al., 2013) is a semi-distributed model that considers the spatial heterogeneity using a distribution function. The meteorological forcing is assumed homogeneous over the catchment area. With shape parameters describing the distributions of soil and epikarst storages, vertical hydraulic conductivities, the separation into diffuse & concentrated recharge, and the groundwater hydraulic conductivities, the model consists of 15 compartments reproducing the spatial variability of recharge and storage dynamics. This value was tested in Hartmann et al., (2012), and shown to be enough to obtain stable integrated discharge and its variability. It has also been confirmed by several following applications (Hartmann et al., 2013, Hartmann et al., 2021, Liu et al., 2021). Each compartment includes three superposed reservoirs (soil, epikarst and groundwater). The first 14 compartments represent the dynamics of diffuse recharge and matrix flow, while the last compartment represents the dynamics of concentrated recharge and conduit flow (Hartmann et al., 2014). The discharge of the main spring is determined by the addition of flow from the groundwater reservoir of all model compartments, which are controlled by the variable groundwater storage coefficients and conduit storage coefficient, respectively. The number of parameters for the calibration is 8.

TCD-Dublin developed a semi-distributed 1D model of the catchment based on the drainage software InfoWorks ICM (Innovyse), which models the hydraulics of the karst conduit network in both open channel and pressurized pipe flow (Gill et al., 2013). The governing model equations are the Saint-Venant equations of conservation of mass and momentum. Diffuse recharge from rainfall is modelled per sub-catchment via a series of reservoirs: rainfall runoff, soil and groundwater stores in series, all yielding different delayed discharges in parallel into the pipe network. Flows can also discharge into permeable pipes, one connected for each sub-catchment to represent the primary and secondary permeability expressed by Darcy’s Law (Morrissey et al., 2020, Schuler et al., 2020). The higher the density of these pipes included in this type of model moves towards becoming a fully distributed.

In the InfoWorks model, the Milandre catchment area is explicitly subdivided into 30 sub-catchments of 0.45 to 0.475 km2. The sizes of the sub-catchments are the result of modelling experience with respect to balancing the numerical stability of the model along with reasonable diameters of the draining pipes and heads. Every sub-catchment generates recharge and storage through a hierarchical sequence of processes related to a) runoff (quick recharge); b) soil store contribution (intermediate recharge); and c) groundwater store contribution (slow recharge). 19 parameters control these dynamics, and in this model, all sub-catchments have the same parameter settings.

2.4. Fully distributed models

In KarstFLOW (Pardo-Igúzquiza et al., 2018) recharge is split between fast recharge and slow recharge. Fast recharge is a percentage of rainfall that reaches the water level through the conduit network instantaneously (for practical purposes). The slow (diffuse) recharge takes place over the whole recharge domain, while fast (concentrated) recharge is limited to specified regions (Fig. 2). Diffuse recharge is estimated by the method of Thornthwaite applied to the percentage of rainfall that does not undergo quick recharge. The model domain is discretized into voxels with absolute altitudes using a digital elevation model. The diffuse flow along each column of voxels is described by the Richards equation (assuming a porous medium). When the non-saturated flow reaches the phreatic zone, the one-dimensional vertical flow couples, as recharge, with a two-dimensional domain where flow through conduits and fractures is simulated by the equations of a two-dimensional Darcian flow in a porous medium. The karst spring discharge is simulated by using a drain. The parameters needed by the model are described in Fig. 2.

  1. Download : Download high-res image (228KB)
  2. Download : Download full-size image

Fig. 2. Parameterization used in KARSTFLOW methodology.

MODFLOW-CFPv2, in short CFPv2 (Fig. 3), is a flow and transport discrete-continuum model which combines CAVE heat and solute transport routines (Liedl et al., 2003) with MODFLOW2005-CFP Mode1 flow code (Shoemaker et al., 2008), considering some additional improvements (Reimann et al., 2018). In CFPv2 the aquifer encloses a network of conduits with turbulent/laminar flow embedded within a matrix of lower permeability (porous medium). Recharge can be split between diffuse recharge through the low permeability porous medium, infiltration through the epikarst reservoir (e.g., Chang et al., 2019), and direct infiltration into conduits and their associated fast-response reservoirs, namely CADS (see Kavousi et al., 2020). Accordingly, several recharge conditions were conceptualized and initiated as numerical model variants of the catchment. In the presented simulations, recharge is calculated by subtracting potential evapotranspiration (ETP) to total precipitation (P), considering a simple degree day snow-cover-snowmelt approach for estimation of hourly snowmelt, as well (for the methodology of degree-day approach, see Rango and Martinec (1995) among many others). It should be pointed out that, the distributed recharge was given here as a single value all across the modeling domain. Moreover, only six input parameters were homogenously assigned to the model, keeping in mind that the overall aim of karst modeling challenge as a proof of concept, not calibration statistics. It should be noted that spatial variation of recharge and model input parameters (which were not also provided with current state catchment data), could improve the performance of CFPv2 process-based model. The inverse problem was solved using the PEST parameter estimation code (Doherty, 2015).

  1. Download : Download high-res image (259KB)
  2. Download : Download full-size image

Fig. 3. Conceptualization on karst aquifers by MODFLOW CFPv2 (See Reimann et al., 2014 for detailed explanation); note that sinkholes and losing streams are not present within the Milandre catchment.

2.5. Comment on the various approaches

The challenge is focused on the recharge-discharge relationship only (mainly the uncertainty of forecasted hydrographs at the system output). Table 1 provides an outlook of the main characteristics of the respective models. Lumped parameter models are thus basically sufficient to address the question. Although semi- and fully-distributed models have been applied, they neither considered the spatial variability of infiltration nor of the aquifer parameters (parameters were assumed to be the same all over the catchment area). Only the spatial variability of the interception of precipitation and evapotranspiration was considered in some of the models.

Table 1. Overview of the main characteristics of the 13 models applied for the Karst Modelling Challenge step 1.

Model summaryApproachMain characteristicsParameter #Calibration typeObjective functionComputational cost (for 1 run)Remark
BRGM, FranceGardeniaLumped, reservoir4 reservoirs9AutomNashSeconds
Uni-Freiburg, GermanyVarkarstSemi-distributed15 compartments8Autommin. residualsSecondsDistribution functions are used to reflect spatial variability
IMT Mines Alès, FranceANN/rec_MLPLumped, neuralRecurrent Deep Multilayer Perceptron381Autom.mean quadratic errorSecondsNot enough hydrologic situations in the database (too short) to design a robust neural-network model.
SISKA-SwitzerlandKRM_1Semi-distributed3 reservoirs + 1 exchanger22ManualVolume conservation, similar shape of peaks and base flowSecondsThe upper reservoir (soil) is semi-spatially distributed
IGME Madrid, SpainKarstFLOWFully distributed2 submodels (recharge + aquifer)12Autom.minimization of residuals15 minAquifer simulated with Darcy's law
TCD Dublin, IrelandInfoWorksSemi-distributedDouble recharge reservoirs (diffuse/concentrated) + pipe flowmin. 19ManualVolume conservation, similar shape of peaks and base flow5-10 min.sSubcatchments linked by a pipe network --> spatial distribution of recharge
KIT-Karlsruhe, GermanyCNNLumped, neuralConvolutional Neural Networks950921Autom.MSESeconds
KIT-Karlsruhe, GermanyNARXLumped, neuralNonlinear AutoRegressive networks with eXogenous inputs371Autom.MSESeconds
SNO KARST, FranceKarstModLumped, reservoir2 reservoirs6 for KMC
(up to 12)
Autom.MSESecondsInfinite characteristic time configuration with fixed ra, emin, kec and calibrated hmax, tau0, alpha
TU-Dresden, GermanyCFP-modifiedFully distributedDiscrete-continuum modelmin. 6Autom.sum of squared errors, SSEMinutes
TU BA-Freiberg, GermanyRCD-SeasonalLumped, reservoir3 reservoirs9Autom.peak max and flow minSecondsCalibration needs 40 seconds (10 runs)
Uni-Zürich, SwitzerlandCHLEMLumped, reservoir9 Elements (flexible)min. 30Autom.RMS difference0.01 secTracks flux, temperature and dye, flexible chaining of elements
KIT-Karlsruhe, GermanyLSTMLumped, neuralLong Short-Term Memory networks11311Autom.MSESeconds

Most models include some parameters which are more or less physically-based, and others which aren’t. InfoWorks, KarstFLOW and MODFLOW CFPv2 are mostly based on physical descriptions of flow with most of their parameters having a physical meaning. However, the physical concepts used in these models are considerably different and simplified compared to processes really taking place in nature.

3. Study site and data

The Milandre karst hydrogeological system (MKHS) is located in the Tabular Jura, in the front part of the Jura Mountains in Northern Switzerland. The site is known as being the karst laboratory used by SISKA for conducting surveys and experiments since early 1990s. The system outputs are formed by two perennial springs: Saivu (20–200 L·s−1, 373 m.a.s.L) and Font (12–600 L·s−1, 369 m.a.s.L). During medium to high water periods water overflows at Bâme spring (0–3000 L s−1, 375 m.a.s.L.). All springs are located on the left (West) side of the Allaine river (Fig. 4). MKHS is fed by a recharge area of ~ 13 km2 (Grasso and Jeannin, 1994, Jeannin, 1996, Perrin et al., 2003) which consists in a limestone plateau at an altitude of ~ 550 m.a.s.L. This area is occupied by forests (~30% of the catchment area), pastures (~30%), cultivated lands (~30%) and urban regions (~5%). The limestone is almost completely covered (95%) by soil with a thickness in the order of 0.3 to 0.5 m in forested parts and up to ~ 2 m in cultivated land. GW recharge is purely autogenic and diffuse (no surface stream and swallow-holes).

  1. Download : Download high-res image (596KB)
  2. Download : Download full-size image

Fig. 4. Map of the Milandre catchment with its main characteristics. The catchment area is about 13 km2 large and encloses the Milandre underground stream.

The mean annual precipitation is 1070 mm as measured by the MeteoSwiss weather station at Fahy, located 7 km away from the center of the catchment area. The mean annual effective precipitation (precipitation minus real evapotranspiration) is 520 mm. An overview of the hydrogeological settings of the area is provided by Kovacs and Jeannin (2003). The aquifer is hosted by Upper Jurassic limestone and is underlain by the Oxfordian marls, which act as a regional aquiclude. As the active conduits lie almost directly above this impervious formation, there is no deep phreatic zone and the system is qualified as a shallow karst system (Perrin, 2003).

The downstream part of the catchment contains a 10.5 km speleological network, the Milandre cave, which hosts a 4 km long perennial cave stream, the Milandrine. This is the main drainage axis of the catchment. The stream flows into a sump around ~ 500 m upstream from the Saivu spring. Two main underground tributaries feed the Milandrine and each of them contributes to 25 to 35% to the cumulated discharge of the springs (Grasso and Jeannin, 1994).

The underground stream is monitored along its course at three places for discharge, temperature and electrical conductivity. Several drillholes located around the main karst conduits, have been monitored showing the groundwater behavior around the conduit network. This data will be important for the next steps of the challenge.

4. Methods

4.1. Workflow of the karst modelling challenge (KMC)

For step 1 of KMC, dedicated to the relationship between meteorological parameters and spring hydrograph, the discharge rates of the three spring of the system (Saivu, Font and Bâme) were simply summed in order to provide one single discharge time series for the whole KHS.

4.1.1. Calibration data

In February 2017, each team received 4 years of meteorological data, from 1.1.1992 to 31.12.1995. For the main meteorological station (Fahy), located 7 km south-west of center of the catchment area, hourly precipitation, hourly temperature and daily PET were given. Hourly precipitation data were also given for the Maira station, which is located in the middle of the catchment area. Daily precipitation were also given for a third station (Mormont) located ~ 2 km south-east of the catchment. Each team also received ~ 2.5 years of hourly discharge data (sum of Saivu, Font and Bâme springs), from 24.9.1992 to 28.3.1995.

This data set was used for a first calibration of the models, as well as for adjusting criteria for the comparison. All teams returned their results, and a first comparison was made in 2018. Based on this first exercise, it was then decided to provide a second dataset for a further evaluation of model performances.

4.1.2. Model evaluation

In June 2019 a new dataset was sent to all teams with two years of data (2014–2015) for the meteorological station of Fahy and for the total combined discharge rates of Saivu, Font and Bâme springs. Data from Maira and Mormont were not provided because they were declared as not useful by the teams after the analysis of the first dataset.

Teams were asked to apply their model to the second dataset with the parameters calibrated with the first dataset, to run the quality criteria, and to write a comment on it. If they wished, they had the possibility to improve the calibration using years 2014–2015.

4.1.3. Test data

Finally, all teams were invited to forecast discharge time series for the year 2016, for which they only received input data. That meant that they could use their calibration data based on the 1992 to 1995 time series data for the simulation (forecast) of 2016, or based on the 2014–2015 time series, or a combination of both. They sent their best model to the first author, who applied the same evaluation criteria to all received time series.

Among the 15 teams originally interested in the challenge, 8 of them sent their results and 5 additional teams joined afterwards, making a total of 13 modelling teams comparing their approaches and results.

4.2. Uncertainty of the data

We briefly discuss here two different types of uncertainties, which are attached to measured data used for the challenge:

4.2.1. Measurement uncertainty

Air temperature is measured with a high precision, better than + -0.1 °C. Precipitation is measured with a high-quality rain gage (calibrated tipping bucket of MeteoSwiss), with a precision in the order of 0.2 mm. Both are measured every 10 min and averaged (T) or cumulated (P) over one hour.

For the spring discharges the water level is measured every 15 or 30 min and is transformed into discharge using a rating curve established empirically for the measurement section. The obtained discharge rates are averaged over one hour. For this step 1 of the challenge, the output discharge rate is obtained by summing the three outlet points (springs) of the MKHS. Rating curves were derived by taking discharge rate measurements at the springs across a range of various rates (or water levels). Most measurements were made by salt dilution method. Measurements are repeated every few years in order to look at potential changes in the rating curves.

Uncertainties on discharge rate data are in the order of ±20% in absolute values. However, this uncertainty must be considered as global and not as a noise related to every single value.

4.2.2. Uncertainty related to spatial heterogeneity

Meteorological parameters are measured at an official meteorological station of the Swiss Meteorological Agency (SMA). The station is located 7 km away from the center of the catchment area, at a similar elevation. Summer showers locally can be quite different from what is measured at the rain gage. However, for most major events, the rain distribution is rather uniform, as given by a comparison with other meteorological stations in the area, including the one located in the middle of the catchment area, which is rather small (10–15 km2) and rather flat. This latter rain gage was compared but not used for the challenge, as its measurements appeared to be less liable than those of the official SMA station.

Land-use is not uniform. As detailed previously the catchment area is covered by forest (30%), pastures and cultivated land (60%), and about 10% of urban area. The effect on recharge of the differences in RET models for these respective regions is not directly measurable, inducing another source of spatial uncertainty.

4.3. Evaluation of model performances

The only aim in the present step of the challenge was to simulate the observed discharge rate at the system output “as well as possible”. The quality criteria (or efficiency criteria) must therefore result from the comparison between the observed and the simulated (forecasted) time series. Six methods or objective functions were considered, providing an evaluation of model performances as a function of different objectives (Madsen et al., 2002, Moriasi et al., 2007, Biondi et al., 2012): Mean Squared Error (MSE), Normalized Mean Squared Error (NMSE), Variance (VAR), Nash criteria (NASH), Kling-Gupta Efficiency criteria (KGE), and Volume Conservations Coefficient (VCC). Each method summarizes the difference between the forecasted and the observed curves into one single value, which is obtained with a different calculation.

During the first exercise using data from the 1990 s, discussions between participants of the challenge showed that only NASH, KGE and VCC were applicable for the comparison. The other criteria are not normalized, and can therefore hardly be compared from one simulation to another. Hence, although they have been calculated, they will not be further discussed in this paper.

VCC is the ratio between the forecasted volume of flow and the observed one. Ideally the value must be 1. Quality classes are given in Table 2 as well as the range between very good to low.

Table 2. Quality classes for VCC, NASH and KGE criteria, and for the necessary effort.

NASH takes into account the ratio between the mean squared error (MSE) and the variance (VAR) as defined by Nash and Sutcliffe (1970). The NASH coefficient ranges from -∝ to 1 (1 being the ideal value). For NASH less than 0, the mean of the observed values is a better indicator than the calibrated model. For NASH = 0, the results provided by the model are as accurate as the mean of the observed values. Usually, the model may be considered as reliable when NASH coefficient is higher than ~0.75.

The Nash criterion (least squares) is well known for giving a lot of weight to flood periods, but in return being less sensitive with regards to the simulation of low water levels. This motivated Gupta et al. (2009) to introduce an improved criteria (KGE) taking into account a normalized distance between the observed and forecasted curves. KGE ranges between 1 and 0.

These three criteria were applied to the whole time series. However, in order to evaluate more closely the apparent strengths and weaknesses of the respective models, we also split the observed curve into four components: peak rising limb, peak recession, base flow, and undetermined. We calculated the quality criteria for the overall time series, as well as for each of its four components.

In addition, a fourth criteria was used to assess the applicability of the respective models, i.e. the effort necessary to set up, calibrate and run the model. Classes are given in Table 2.

A “final note” was arbitrarily calculated using a weighted mean between KGE (weight 2/5), VCC (weight 2/5) and Nash (weight 1/5).

Relative errors are also used for characterizing the difference between the forecasted and the observed values. For each value (hour) the difference is calculated and divided by the larger value of observed or forecasted discharge rate. Doing so, a forecasted value 3 times lower than the observed one makes a relative error of − 300%, and a forecasted value 3 times larger than the observed one a relative error of + 300%.

4.4. Time step

Hourly Discharge Measurements (HMs) of the Milandre springs show fast fluctuations on a time-scale which is shorter than a day which would suggest that model simulations should be computed at an hourly time-step.

However, some models are designed for daily time-steps and provide daily values (DVs). These may be compared to the reference (observed) time series (HMs) either by taking daily instantaneous values (DIVs), e.g. the measured value at noon, or by averaging values over 24 h (DAVs).

In the framework of KMC, results of the simulations had to be compared to HMs for computing the efficiency criteria. Therefore, simulated DVs were retransformed into hourly values (HVs). A linear decomposition of DVs into HVs was thus necessary.

Conversions from hourly values to daily values and then to hourly values again may significantly affect the quality of the simulation. Two biases were tested: the effect of resampling/decomposing, and the way to resample (from instantaneous or from average values).

One question raised by participants was if it was better to use average or instantaneous values. Therefore, we made a comparison between reference HMs and hourly values issued from the resampling of DIVs and DAVs. The new time-series are called HDIVs and HDAVs respectively (Fig. 5).

  1. Download : Download high-res image (102KB)
  2. Download : Download full-size image

Fig. 5. Workflow applied for testing the effect of resampling and of using daily average values or instantaneous values (e.g. the measure value at noon).

Fig. 6 shows that significant differences do exist between HMs, HDIVs and HDAVs. Most of the time, instantaneous values (HDIVs) show a better fit with HMs than average values (HDAVs): HDAVs systematically underestimate peak-flows and often overestimate the rising limb of the floods. HDIVs usually show a better fit with peak-flow except when variations of peak-flow are clearly shorter than one day (e.g. the flood on Dec. 5th, 1992). For this case, the averaged value may provide a better fit. HDIVs better reproduce low-flows than HDAVs.

  1. Download : Download high-res image (234KB)
  2. Download : Download full-size image

Fig. 6. Zoom on the three time-series HMS, HDIVs and HDAVs (period 1/11/1992 to 31/12/1992).

In summary HDAVs show a more systematic bias than HDIVs.

HDIVs and HDAVs are thus compared to HMs using the provided efficiency criteria (see Table 3). In both cases it appears that VCC (volume conservation coefficient) is well conserved. NASH and KGE are about ~0.9 to ~0.95. meaning that conversion from hourly to daily and then to hourly values again degrades the time series in the same order of magnitude of what is usually targeted (and obtained) with (the best) simulations.

Table 3. Comparison of the efficiency criteria for HDIVs and HDAVs.

Empty CellFlow componentsEmpty CellEmpty CellEmpty Cell
HDIVsRising limbBase flowFlood recessionUndet. FlowTotal
VCC0.951.011.011.041.002
NASH0.8160.9850.9170.9270.92
KGE0.860.960.950.930.95
Empty CellFlow componentsEmpty CellEmpty CellEmpty Cell
HDAVsRising limbBase flowFlood recessionUndet. FlowTotal
VCC0.971.010.991.060.995
NASH0.8180.9720.9180.9340.91
KGE0.830.960.900.910.92

It should be observed that NASH and KGE for HDIVs are slightly better than those computed for HDAVs. VCC is a bit lower for HDAVs, which makes sense as HDVAs systematically underestimate peak flows.

Our advice for KMC was then to work at hourly time-step as far as possible for simulations at the Milandre test site. But if this was not possible, then teams were encouraged to use DIVs from HMs to compare daily values obtained with their models.

Finally, as all teams were required to compute their final efficiency criteria at an hourly time step, their simulated results, if calculated at daily time-step, needed to be resampled into hourly data using a linear decomposition.

5. Results

5.1. Overall performance

As shown on Fig. 7 the different model forecasts generally seem to reproduce the correct order of magnitude of discharge rates for most of the prediction period. Models tend to underestimate discharge during recessions occurring after significant high-water conditions (e.g. March and June on Fig. 7).

  1. Download : Download high-res image (195KB)
  2. Download : Download full-size image

Fig. 7. Forecasted discharge rates in red (envelope of results coming from all models) compared to the observed rates (Black line). Most observed peaks are forecasted. Models tend to underestimate discharge during recessions occurring after significant high-water conditions (e.g. March and June). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Given that most models provided reasonable results, quality criteria were then applied in order to identify strengths and weaknesses of the respective models. The following comparison takes into account only the best prediction for the year 2016 provided by the participating teams. Most are obtained with calibration on years 2014–2015, but some used a mix with 1992–1995 as well.

In Table 4 model results have been classified according to their “final score” as defined in chapter 3. Effort is given as indication. Concerning KGE and VCC all models are acceptable to good, but none is very good.

Table 4. Comparison of the model performances. All models produced acceptable results. Color ratings are given in Table 2.

Forecasted curves for the year 2016 are given for all individual models in Fig. 8. The log scale provides a better view of the relative errors between the forecasted and the observed curve than the linear scale.

  1. Download : Download high-res image (592KB)
  2. Download : Download full-size image

Fig. 8. Forecasted hydrographs of the thirteen models applied for KMC-step 1 compared to the observed curve (pale red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

All models reproduce the overall shape of the curve of discharge variations, covering periods of time of low flow and other periods with peaks. The amplitude of the peaks is in the right order of a magnitude. Low flow conditions are correctly identified. However, a closer look shows that low flow conditions are rather poorly forecasted by most models whereby often the forecasted discharge rate is too low and too flat.

It appears that curves predicted by reservoir models and semi-distributed models reproduce somehow better the observed behavior with many peaks from January to June, and less peaks from July to October. InfoWorks and KRM1 are semi-distributed models with their recharge processes based on reservoir models; their curves are thus quite similar to those of reservoir models. VarKarst is based on a different approach (it explicitly considers spatial variability), and it can be seen that the peaks are pretty well predicted, but the simulation does not reflect the observed behavior well during low flow conditions. RCD-seasonal oversimplifies recession during summer with just a single long recession. The distributed models fail at forecasting low flow conditions as too many little peaks are produced, and the subsequent recessions are much too quick. Peaks are underestimated in winter and overestimated in Summer. CHLEM, which is a reservoir model, produces curves which are similar to those of distributed models. Lumped models are not bad for the first half of the year, but have a significant drift in the second part of the year, with, ANN-2 performing better overall than the three others.

Relative errors are given in Fig. 9. Errors on peaks are of short duration, but may be of considerable intensity (usually larger than 400%, positive or negative, for some events along the year). The flow-rate during the recession in March is underestimated by most models. CNN is the only model providing a prediction close to observations for this event. VarKarst, KarstMOD and CHLEM are particularly low. The main recession in July tends also to be underestimated by most models with Gardenia, Karstflow and KRM1 performing better than the others.

  1. Download : Download high-res image (1MB)
  2. Download : Download full-size image

Fig. 9. Relative errors of the respective models for year 2016. All predictions show periods of time with errors larger than 200%. The green band shows the domain with + -50% of relative error. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

It is noticeable that most models using neural-derived approaches display a significant drift in the last months of the simulations except the MLP model.

By classifying the relative errors, we can display the percentage of the forecasted data (hours/year) for which the error of the model is lower than a value (Fig. 10). For instance, the LSTM model (intense pale blue) has 22% of the forecasted values with an error lower than 50% (e.g. between 67 and 150 L/s instead of 100 L/s). Gardenia has 72%: the higher is the curve, the better the performance of the model.

  1. Download : Download high-res image (366KB)
  2. Download : Download full-size image

Fig. 10. Percentage of forecasted data (of time) as a function of relative errors. The best model has 72% of forecasted data within a relative error of + -50%. The worse model has only 22% of forecasted data with a relative error lower than 50%. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

5.2. Performance by flow components

The observed times series was divided into four components: rising limb of the peaks, recession, base flow and undetermined flow. Base-flow was empirically determined when the discharge rate is decreasing over five successive hours with a rate lower than 1 L/s per hour. The latter component includes points which are not clearly related to one of the other three classes.

For the year 2016 the significance of the four components is given in Table 5 (in hours). Variations from one year to the other are in the order of ±5%:

Table 5. Part of the total hydrograph in the respective components.

Flow componentPercentage in hoursPercentage in volumeMean discharge rate
Rising limb1121608 L/s
Recession3357555 L/s
Base flow99106 L/s
Undetermined2813156 L/s

Each value of the observed curve is attributed to one single flow component. For the forecasted curves, the components of the observed curve are used for the separation, meaning that all models are compared on the basis of the same components. As can be seen on Fig. 11, most values belonging to the undetermined flow correspond in fact to base-flow conditions. They could not be attributed to base-flow because they include slight increases of the discharge rate, which may be related to some technical problems (e.g. temperature variations on the electronics, or changes of the atmospheric pressure which are not perfectly corrected).

  1. Download : Download high-res image (153KB)
  2. Download : Download full-size image

Fig. 11. Decomposition of the time series of observed discharge-rates. The undetermined component mainly corresponds to base flow. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The respective models have specific results for each component. Gardenia has among the best performances for most criteria and components. This is less the case for Varkarst, especially concerning volume conservation (VCC). Most models have one component with KGE close to or below 0.5. Nash criteria is low for base flow in most models. CHLEM, CFP-modified and RCD_seasonal, three models with the lowest scores, are better than most other models for base-flow according to VCC and KGE.

The Nash criteria is quite poor for most models as this criterium cannot be good for a trending time series (as the mean value is taken as reference in the denominator of the Nash-Sutcliff formula).

Other quality criteria have also been calculated for the respective models, mainly Mean squared error and the variance test. However, these criteria do not provide more information than those presented in Table 6.

Table 6. Table showing the three criteria applied to the four components of the 13 models compared for the challenge. Models are sorted according to their global performance (see Table 4).

5.3. Comparison of calibrated parameters

Given that the models are based on different concepts, their parameters are all different and cannot be compared directly. Some aspects, however, can be roughly compared for several models, i.e. catchment size, storage capacity, conduit/diffuse flow components.

5.3.1. Catchment size

In karst regions catchment size is never known exactly, and so this was considered to be a calibration parameter by some of the models. The calibrated size area varies between 11.6 and 16.24 km2, which ranges between-9 and + 28% of the assessed surface area (12.7 km2). The main reason for these variations is related to the model applied for assessing evapotranspiration: the size was calibrated to adjust the water balance on the calibration period.

Neural networks models do not provide any catchment size. The role of this parameter is approached by a series of unidentified functions and parameters.

5.3.2. Storage capacity

The attenuation between recharge and discharge is related to the characteristics of the storage in all models. However, the way it is generated is different in the respective models. Most models have quick and slow storages included somehow, which reflects what is commonly accepted in karst aquifers since the 1970s as summarized e.g. by Bakalowicz (2005). Their number and characteristics are, however, different and hardly comparable. It may be only mathematical, or based on the infilling of a reservoir, or on the routing of water, etc. Somehow, all models generate peaks within 0.2–5 days after a rain event, and a base-flow usually corresponding to the emptying of some kind of reservoir (or mathematical function) with slow recession keeping water flowing out of the system over more than 6 months without recharge. The storage volume of this slow component is in the order of 100–250 mm.

5.3.3. Conduit/diffuse flow components

Another long-lasting debate among karst hydrogeologist (see e.g., Bakalowicz, 2005) is about the ratio between the quick-flow component (assumed to be related to conduits) and the slow-flow component (assumed to be related to matrix or diffuse flow). From our results we were hoping to be able to compare the amount of recharge flowing through both components of the respective models. However, as every model has a different architecture and concept, this could hardly be done. The only thing we can compare is the amount of peak flow (rising limb and recession) and base-flow (including undetermined flow) for each model. The comparison shows, however, that differences between models are not very significant, the amount of peak flow ranging between 64 and 78%, and the volume of base flow ranging between 22 and 36%.

6. Discussion and model comparison

The comparison shows that forecasted curves with similar values for their quality criteria may have different shapes. This illustrates the fact that quality criteria for model comparison could still be improved. Our first assessment of the results was that all models produced fair to good simulations. However, this assessment is probably rather qualitative than quantitative. Indeed, when looking at relative errors, we realized that most models produce values with errors larger than 50% for about 50% of the time or more. Results are clearly worse during low-flow conditions. Absolute differences are low during these periods of time, but relative ones are significant.

By looking at the respective curves (Fig. 8) and at the relative errors (Fig. 9), one can say that the three criteria we selected are not sufficient to characterize the overall quality of the results. For instance, the shape of the curve forecasted by the VarKarst model reproduces very well the peaks, be very poorly the values during the dry-season. Gardenia is better for base-flow, but still include significant differences in March and October of the simulated year. The fact that most models obviously tend to underestimate or poorly reproduce low-flow (e.g., LSTM or VarKarst) could be related to the way calibration is implemented. It is usually mainly targeted on minimising the squared errors metric, which gives much weight to the peaks and not much to low-flow values. A criterium minimizing relative errors instead of absolute errors could possibly improve the weighting of base-flow conditions in the overall results. Another option could also be to calibrate the parameters which control base-flow first as they not dependent much on the recharge model, and then to calibrate the other parameters on the residuals of this first calibration step. In the case of RCD_seasonal and possibly VarKarst, the main reason for the poor forecasting of low-flow conditions seems to be related to the recharge model, which completely fails to produce any recharge between June and November or December. The water deficit in the recharge reservoir is obviously too high over this period compared to the reality.

Interestingly the spatial distribution of precipitation appeared not to be significant and was abandoned along the process in agreement with all participants. Obviously, the number of recharge event for which differences are significant (e.g. summer storms) is low. This may be the case for the Milandrine KHS, which is small and located in a temperate and hilly region. This may be different for a larger catchment with more relief and contrasting climate.

KRM_1 model is the only model for which land-use was taken into consideration. However, this does not appear to improve the results in a way to produce better results than the other models.

The number of parameters in the respective models ranges between less than 10 and several thousands. Despite automatic calibration/optimization procedures, models with less parameters tend to be better than the models with many parameters, i.e. distributed models did not provide better results than global ones. Not considering catchment characteristics in the respective models could be a reason for this. Therefore, it seems that the transformation between meteorological parameters and discharge rates is mainly controlled by a few general characteristics of the system (such as size, recharge process and one or two storages) rather than by the detailed spatial distribution of parameters controlling these). Modellers using neural networks, with many parameters, claimed that data delivered for calibration were not sufficient. As no basic structure is provided in their models, compared to the models based upon a more physical function, they have to build this structure using data. They require therefore at least 15 years of data without any gap, which is challenging!

These observations raise the question of uncertainty and overfitting associated with numerical models. In fact, we can probably conclude that the best models are as good as the data we provided.

MODFLOW-CFPv2 and KarstFLOW are fully distributed models. These models are based on the principle that the transformation of the input signal into the output one is mainly related to the characteristics of the aquifer. This approach is more deductive than for the other models because these models are constructed from the supposed characteristics of the catchment and of the aquifer. A consequence is that the recharge part of the models is simplified. This is the main reason why peaks in the summer season are strongly exaggerated with both models (see Fig. 8). For the given exercise, this approach is also not the most efficient as it requires more work for digitizing a lot of data, but at the end the results include more bias than many other models because recharge was not considered properly. Would a spatial discretization of the model parameters improve the results? Whilst this was not explicitly investigated in this research, we can argue that the main discrepancy of these models is related to (too) large peaks in summer, which seem to be more related to an oversimplification of the recharge model, rather than to the spatialization of the model parameters. Recharge seems therefore really to be the key factor in the present case.

On the contrary, lumped parameter models are more inductive because they include less or no constraining hypotheses. The disadvantage is that none of their parameters can be compared somehow with natural characteristics (not even the size of the catchment area!). They provide purely functional results. Another problem is that they require long and continuous data sets for learning the input–output relationship of a natural flow system.

In semi and fully distributed models, the spatial distribution of recharge (land-use, vegetation, precipitation) can be taken into account. However, as the performance of these models turned out to be no better than results of similar models (reservoir) which did not take it into account, one can infer that, at least in the present case-study, the role of spatial distribution of the recharge is low. The effect of spatial distribution was, however, not analysed in any depth and so could be further investigated (e.g., Bittner et al., 2018, Gill et al., 2020). As the Milandre site is rather flat, small and without allogenic recharge, the use of a complimentary catchment area with more contrasted conditions would probably be better for this analysis. It should be noted that the relative significance of the spatial distribution of recharge function in any models will vary according to the characteristics of the karst catchment being simulated: for example, whether there is significant allogenic recharge or not, the size of the catchment area, and/or the steepness of the topography may all be important factors.

An evaluation of the number of parameters hydrologically constrained in a rainfall-runoff model in surface hydrology was investigated by Jakeman and Hornberger (1993). They show that most hydrological systems are well described by models with 4 to 6 parameters. The contribution of any supplementary parameter seems to be within the data uncertainty related to monitoring data. This explains somehow the fact that, in our study, models with many parameters (greater than20) do not provide better results than those with less parameters. In karst, it is difficult to obtain reasonable results with less than 6 parameters: three for the routing and three for the transformation of precipitation to recharge. Because of some visible threshold processes in karst systems, 9 parameters seems to be a more reasonable minimum. Table 1 indicates that some of the models had less than 9 parameters, however they include implicit assumptions fixing some of the parameters (e.g. distribution functions, fixed topology or geometry, etc.). Jakeman and Hornberger (1993) also suggest that a standardized simple model could be applied systematically to any hydrological system. This would help comparing systems’ characteristics and identifying drifts in data related to technical or natural reasons (e.g. climate change). The model comparison presented in our paper could therefore be a good basis for defining such a standardized karst model.

7. Conclusion and outlook

7.1. Conclusion

Before making a conclusion, it is important to have in mind that the comparison between the models was only carried out a one-year simulation period. Due to the non-stationary nature of hydrologic phenomena, some models could be very good for one kind of situation and not so good for others. The presented interpretations must thus be considered as “snapshot” interpretations.

Compared to a similar exercise conducted for surface catchment (Holländer et al., 2009) all models performed reasonably well! Although they are based on strongly different modelling approaches all models lead to reasonable and somehow similar results for the proposed exercise (input–output hydrological modelling). Despite considerable simplifications with respect to reality, Gardenia provides excellent results in the present case. It fits the observed data almost within the range of uncertainty on the data themselves. However, despite very good scores, the forecasted hydrograph displays significant differences with the observed one, especially in low water conditions. This is even more obvious for results issued from VarKarst model.

Apart from the different aspects of the modelling approach it can also be concluded from the modelling exercise that recharge is crucial in the modelling of any precipitation-discharge relationship. More specifically it can be observed that the real evapotranspiration must be carefully approached, otherwise models will fail at reproducing discharge rates from precipitation data, at least in temperate regions. It can even be inferred that hydrogeological processes taking place within the aquifer are subsidiary compared to recharge processes taking place in the soil and epikarst.

Lumped parameters models are interesting options if they can include non-linear functions and adequate parameters for taking evapotranspiration into account. Regarding ANN models, their limitation is that they require long and complete datasets for calibration and validation, and that they usually do not provide any physical indication about the real system, even if this “knowledge extraction” has been successfully attempted (as suggested by Kong-A-Siou et al., 2013, Darras et al., 2015) these models remain mainly functional. Taking into account the limited duration covered by the available data, several ANNs show an acceptable behavior.

From the comparison, as well as from tests of some models on other sites, it also seems that the order and priority in which the calibration of the models is carried out has a strong effect on the result. Calibrating base-flow first and peak-flow afterwards appears to improve the results. This point was not properly investigated in this study and should be further examined. Some physical reasons could support this idea as storage characteristics of karst aquifer are probably relevant concerning base-flow. As this part of the hydrograph is more independent on recharge processes, it may be meaningful to calibrate the storage components of the models first using one or two storage reservoirs. This is also supported by ideas derived from Kovacs et al. (2005). Then the recharge model must focus on modelling the difference between precipitation and the assessed base-flow.

It should also be highlighted here that the scoring scheme used to evaluate of strengths and weaknesses needs to be adapted for the purpose for the model is being developed for. Models used to provide flood assessments may not be assessed in the same way as those models focussing on the implications of low flow/drought conditions. The scores for the respective flow components (Table 5) shows the differences between the respective models.

It should be noted here that some of the modelling approaches could be combined in order to use their respective strengths. For instance, semi-distributed models could be used in combination with neural methods in order to improve their recharge assessment component. Such optimizations will need to be adjusted to specific conditions of the site, available data and question to be addressed.

7.2. Outlook

In the case of the Milandre Catchment area, only one model (KRM1) within the challenge attempted to take the spatial distribution of evapotranspiration (land-use) into account. This did not improve the model results in a measurable way. However, as real evapotranspiration is very significant, and as we know that this differentiation is necessary in catchments with stronger altitudinal gradients, it would be meaningful to further investigate this aspect.

This report covers the very first step of a model comparison. It was focused only on the temporal discharge response of the spring to recharge events. We would like to expand the challenge to further aspects of the modelling of karst hydrogeology. The next step is to take into account the spatial distribution of flow and heads within the aquifer. A first question with this respect is the fact that the discharge rate given at the spring for the challenge in reality corresponds to the flow-rate coming out of three different springs: Font, Saivu and Bâme, which are located several hundred meters apart. We also have head measurements at various locations within the aquifer, and it would be interesting to see how various modelling approaches could reproduce those data. For this new exercise, it seems that lumped models will be hardly usable, and that distributed models will be more adequate. It will be therefore interesting to see how the respective teams will adapt their approach for this new challenge.

Once heads and discharge rates will be approached, we will try to simulate flow velocities and some transport processes for which we also acquired data. This may also result in better performance of the models to the overall flows simulated in step 1 as some parameters will be better constrained.

CRediT authorship contribution statement

Pierre-Yves Jeannin: Conceptualization, Methodology, Formal analysis, Investigation, Resources, Data curation, Writing - original draft, Visualization, Supervision, Project administration. Guillaume Artigue: Software, Validation, Formal analysis, Investigation. Christoph Butscher: Software, Validation, Formal analysis, Writing - review & editing. Yong Chang: Software, Validation, Formal analysis. Jean-Baptiste Charlier: Methodology, Software, Validation, Formal analysis, Writing - review & editing. Lea Duran: Software, Validation, Formal analysis. Laurence Gill: Methodology, Software, Validation, Formal analysis, Writing - review & editing. Andreas Hartmann: Methodology, Software, Validation, Formal analysis, Writing - review & editing. Anne Johannet: Methodology, Software, Validation, Formal analysis, Writing - review & editing. Hervé Jourde: Software, Validation, Formal analysis. Alireza Kavousi: Software, Validation, Formal analysis, Writing - review & editing. Tanja Liesch: Software, Validation, Formal analysis. Yan Liu: Software, Validation, Formal analysis. Martin Lüthi: Software, Validation, Formal analysis. Arnauld Malard: Methodology, Software, Validation, Formal analysis, Writing - review & editing. Naomi Mazzilli: Software, Validation, Formal analysis, Writing - review & editing. Eulogio Pardo-Igúzquiza: Methodology, Software, Validation, Formal analysis, Writing - review & editing. Dominique Thiéry: Software, Validation, Formal analysis. Thomas Reimann: Methodology, Software, Validation, Formal analysis, Writing - review & editing. Philip Schuler: Software, Validation, Formal analysis. Thomas Wöhling: Software, Validation, Formal analysis. Andreas Wunsch: Methodology, Software, Validation, Formal analysis, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was supported by the Swiss National Research Foundation (SNF) for Swisskarst and Thermokarst projects (Grant Numbers 406140-125962 and 200021_188636), the Deutsche Forschungsgemeinschaft, DFG for the iKarst project (Grant Numbers: LI 727/31-1 and RE 4001/2-1). The work of E. Pardo-Igúzquiza was supported by research project PID2019-106435 GB-I00 of the Ministerio de Ciencia e Innovación of Spain. This TCD-Dublin was conducted within the Irish Centre for Research in Applied Geosciences (ICRAG), supported in part by a research grant from Science Foundation Ireland (SFI) under Grant Number 13/RC/2092. The authors would like to thank the French Karst National Observatory Service (SNO KARST) initiative at the INSU/CNRS for their diffusion of KarstMod platform. It strengthens dissemination of knowledge and promote cross-disciplinarily research on karst systems at the French national scale.

The data used in the project were acquired along several projects supported by the Federal Roads Office (FEDRO) and by the University of Neuchâtel. Marc Hessenauer (MFR SA, Delémont) and Pierre-Xavier Meury (Géo & Environment, Delémont) also participated to the data collection. The caving club (Spéléo-Club Jura was also very supportive). We also would like to thank Prof. Carol Wicks and an anonymous reviewer as well as Dr. Junbing Pu (associate editor) for their constructive comments.

Appendix. Supplementary data

What’s this?

The following are the Supplementary data to this article:

Download : Download Acrobat PDF file (3MB)

Supplementary Data 1.

Download : Download zip file (7MB)

Supplementary Data 2.

Referencias

Citado por (30)

Ver todos los artículos de citas en Scopus
1

El reto está sucediendo: cualquier persona interesada puede utilizar los datos disponibles en el material suplementario para este artículo, y puede enviarla/su resultado a la autora de los corrs para obtener los valores de los criterios.