# Water cycle modelling strengthened by probabilistic integration of field data for groundwater management of a quite unknown tropical volcanic hydrosystem

Marc Dumont, Valérie Plagnes, Patrick Lachassagne, Roger Guérin, Bayu

Nugraha, Febriwan Mohamad, Ludovic Oudin, Arif Fadillah, Danièle

Valdés-Lao, Gilles Brocard, et al.

► To cite this version:

Marc Dumont, Valérie Plagnes, Patrick Lachassagne, Roger Guérin, Bayu Nugraha, et al. Water cycle modelling strengthened by probabilistic integration of field data for groundwater management of a quite unknown tropical volcanic hydrosystem. Comptes Rendus. Géoscience, 2023, 355 (S1), pp.207-229. 10.5802/crgeos.192. hal-04888430

HAL Id: hal-04888430

https://hal.science/hal-04888430v1

Submitted on 15 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Comptes Rendus
Géoscience

Sciences de la Planète

Marc Dumont, Valérie Plagnes, Patrick Lachassagne, Roger Guérin, Bayu Nugraha, Febriwan Mohamad, Ludovic Oudin, Arif Fadillah, Danièle Valdès, Gilles Brocard, Jean-Luc Bonjour, Mohamed Saadi, Anne-Sophie Esneu, Aswar Muhammad, Hendarmawan and Nathalie Dörfliger

Water cycle modelling strengthened by probabilistic integration of field data for groundwater management of a quite unknown tropical volcanic hydrosystem

Volume 355, Special Issue S1 (2023), p. 207-229

Online since: 26 January 2023

**Part of Special Issue:** Geo-hydrological Data & Models

**Guest editors:** Vazken Andréassian (INRAE, France), Valérie Plagnes (Sorbonne Université, France), Craig Simmons (Flinders University, Australia) and Pierre Ribstein (Sorbonne Université, France)

https://doi.org/10.5802/crgeos.192

<figure><img src="image_2.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/

<figure><img src="image_3.png" /></figure>
Research article

Geo-hydrological Data & Models

Water cycle modelling strengthened by probabilistic
integration of field data for groundwater
management of a quite unknown tropical volcanic
hydrosystem

Marc Dumont<sup>a, *, a, b</sup>, Valérie Plagnes<sup>a, *</sup>, Patrick Lachassagne<sup>a, c</sup>, Roger Guérin<sup>a, *</sup>, Bayu Nugraha<sup>d</sup>, Febriwan Mohamad<sup>b, d</sup>, Ludovic Oudin<sup>a, *</sup>, Arif Fadillah<sup>e</sup>, Danièle Valdès<sup>a, *</sup>, Gilles Brocard<sup>f, *</sup>, Jean-Luc Bonjour<sup>g</sup>, Mohamed Saadi<sup>a, h</sup>, Anne-Sophie Esneu<sup>a</sup>, Aswar Muhammad<sup>e</sup>, Hendarmawan<sup>d</sup> and Nathalie Dörfliger<sup>g</sup>

<sup>a</sup> Sorbonne Université, CNRS, EPHE, UMR 7619 METIS, F-75005 Paris, France

<sup>b</sup> Urban and Environmental Engineering, Applied Geophysics, University of Liege,
B-4000 Liege, Belgium

<sup>c</sup> HSM, Univ. Montpellier, CNRS, IRD, IMT Mines Alès, Montpellier, France

<sup>d</sup> Faculty of Geological Engineering, Universitas Padjadjaran, Jatinangor, Sumedang
45363, Indonesia

<sup>e</sup> Danone Aqua group, Department of Water Resources, Jakarta, Indonesia

<sup>f</sup> Archéorient, Maison de l'Orient et de la Méditerranée, Université Louis
Lumière-Lyon 2, France

<sup>g</sup> Water Institute by Evian, Water Resources and Sustainability Division, Danone
Waters, Evian-les-Bains, France

<sup>h</sup> Institut de Mécanique des Fluides de Toulouse (IMFT)—Université de Toulouse,
CNRS-INPT-UPS, Toulouse, France

E-mails: mdumont@uliege.be (M. Dumont), valerie.plagnes@upmc.fr (V. Plagnes), patrick.lachassagne@umontpellier.fr (P. Lachassagne), roger.guerin@sorbonne-universite.fr (R. Guérin), geo.bayunugraha@gmail.com (B. Nugraha), febriwan@unpad.ac.id (F. Mohamad), ludovic.oudin@sorbonne-universite.fr (L. Oudin), arif.fadillah@danone.com (A. Fadillah), daniele.valdes_lao@upmc.fr (D. Valdès), gilles.brocard@mom.fr (G. Brocard), jeanluc.bonjour@danone.com (J.-L. Bonjour), mohamed.saadi@toulouse-inp.fr (M. Saadi), anne-sophie.esneu@ifp.fr (A.-S. Esneu), azwar.muhammad@danone.com (A. Muhammad), hendarmawan@unpad.ac.id (Hendarmawan), nathalie.dorfliger@danone.com (N. Dörfliger)

<footnote>* Corresponding author.</footnote>
**Abstract.** Andesitic volcanic hydrosystems in Indonesia are mostly hydrogeologically unknown despite their socio-economic importance. The development of robust and easy-to-implement methodologies to conceptualize and quantify the water cycle components becomes a prerequisite for their sustainable management.

We developed a lumped hydrological model to mimic the structure and functioning of a previously unknown hydrosystem located on the flanks of the Salak volcano (West Java). The structure of the aquifers was revealed with electrical resistivity tomography. The distinction between springs fed by the extensive artesian aquifer and others fed by shallow perched aquifers was obtained mostly using hydrochemistry. The elevation of the recharge area was identified using isotopic analysis of spring water.

After designing the hydrological model structure, we carried out a probabilistic parameters exploration using the multiple-try differential evolution adaptive Metropolis algorithm to calibrate aquifer discharge. Multiple Markov chains allow a better exploration of the parameter values. The Bayesian approach provides the best water cycle simulation with a parameter uncertainty analysis, improving the accuracy and representation of the water cycle appropriate for previously unknown hydrosystems.

**Keywords.** Groundwater, Groundwater management, Multidisciplinary approach, Lumped model, Probabilistic, Andesitic setting.

Manuscript received 6 April 2022, revised 13 July 2022 and 17 October 2022, accepted 28 November 2022.

# 1. Introduction

Freshwater is an indispensable component of ecosystems and human activities' essential needs. The sustainability of human societies depends on the constant supply of clean fresh water. Currently, only 10% of the Earth's renewable freshwater is used by humans, while the other 90% sustains natural ecosystems and biodiversity [de Marsily, 2020]. To improve the human management of renewable water and sustain societal development, we need to increase our understanding of water interaction with the other components of the water cycle. Nevertheless, the consideration of water fluxes in all natural compartments is hampered by the availability of relevant data and/or conceptual models. This problem is all the more important in countries where observation networks of the water cycle have only recently been developed, and often sparsely. In such conditions, the first challenge is to accurately monitor and quantify the water fluxes from the atmosphere to the subsurface. While the dichotomy of the water cycle between the water reservoirs may be difficult to grasp in complex geological settings, our study presents a methodology combining multidisciplinary field measurements and probabilistic simulations to (i) quantify the water cycle within the different compartments of a previously unknown watershed, and (ii) define field monitoring needs to increase the accuracy of the water cycle representation.

Our study focuses on a watershed located along the flank of Mount Salak, an andesitic volcano at the western part of Java island, Indonesia (Figure 1A).

In this tropical region, the humid tropical climate generates significant rainfall (from 2 to 5 m per year) year-round accompanied by high evapotranspiration (from 1 to 1.4 m per year). Java Island climate is highly impacted by Asian and Australian monsoon leading to intense rainfall during the wet season [from November to March; Aldrian and Djamil, 2008]. Conversely, the dry season is characterized by low rainfall inducing notably a decrease in surface water and spring discharge. Population growth and urban sprawl lead to a constant increase in anthropogenic pressure on water resources. Upstream, surface water is largely used for the irrigation of rice fields [Khasanah et al., 2021]. Downstream, industrial and urban wastewater discharges to the rivers have led to a degradation of surface water quality [Buwono et al., 2021, Suriadikusumah et al., 2021]. This expansion of human activities has led to a significant increase in anthropogenic pressure on groundwater resources throughout the watersheds [Khasanah et al., 2021]. This situation could lead to water supply shortages during the dry season and degradation of groundwater quality.

In this context, it becomes crucial to understand the distribution of renewable water flows through the different compartments of the watershed. However, such quantification is hampered by the heterogeneity of the subsurface, which is composed of volcanic rocks and/or detrital formations that host spatially variable heterogeneous aquifers [Lachassagne et al., 2014, Selles et al., 2015]. The summit of the volcanoes, characterized by steep slopes above 20° and up to 40°, is usually mainly com-
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (A) Location map of the study area in Indonesia. (B) Map of the investigated surface watershed (delineated in white) with location of surveyed groundwater springs (solid circles), and of the twelve numbered ERT profiles. The springs are colored as a function of the water electrical conductivity and labeled with spring discharge (in L/s).</figcaption></figure>
posed of andesitic lava flows and pyroclastic projections. At the foot of these steep slopes, lava flows transition to alternations of detrital and distal ash fallout layers that form peripheral detrital aprons. Lower down, the aprons become alluvial floodplains, with alternating laharic formations and fine-grained soils, which can prevent the upwelling of underlying groundwater, creating artesian aquifers [Toulier et al., 2022]. This layering of geological formations related to the distribution of volcanic formations spatially disconnects the surface watersheds from underground ones [Charlier et al., 2011, Dumont et al., 2021, Vittecoq et al., 2014]. At the same time, geological discontinuities along slopes may allow the release of groundwater at spatially-restricted emergence points, which may substantially contribute to surface water flows [Vittecoq et al., 2019]. These andesitic hydrosystems are complex due to their lateral and vertical heterogeneity and are poorly or not instrumented.

Hydrogeological models are useful tools to better assess water fluxes between the compartments of the water cycle. The existing modelling strategies cover a large spectrum organized along a two-dimensional continuum according to spatial resolutions and process complexities [Hrachowitz and Clark, 2017], from lumped, conceptual models to process-based models describing the system from a microscale perspective. Both approaches were applied in andesitic settings. Lumped or semi-distributed holistic reservoir models were developed to accurately capture water exchanges between the atmosphere, the soil, runoff and the subsurface [Cerbelaud et al., 2022, Charlier et al., 2008, Gómez-Delgado et al., 2011, Mottes et al., 2015]. Generally calibrated on river flows, they simplify the complexity of the aquifers that is induced by the andesitic volcano structure. Process-based models solving Richards' equation and Darcy's law have also been developed to represent physical processes [Chang et al., 2020, He et al., 2008, Lo et al., 2021, Sapapona et al., 2018]. In this case, the finite element mesh is constructed by incorporating geological and hydrogeological features to simulate the water table with a calibration on monitored boreholes. Nevertheless, the surface/underground exchanges are significantly simplified limiting the representativeness of the aquifers in the watershed water cycle. Both approaches rely heavily on data, either to calibrate model parameters or to prescribe the geometry and

hydrodynamic properties of geological and hydrogeological features. To calibrate and/or to validate the models, continuous river discharge or piezometric data are required over long record periods to represent annual and interannual variability. This traditional calibration framework does not allow for the integration of non-continuous data or expert knowledge that could improve the performance of the simulations [Gharari et al., 2014].

To model andesitic hydrosystems, we argue that representing the surface/underground exchanges is mandatory while at the same time keeping the representativeness of the aquifers consistent with the heterogeneity of the environment. In addition, the management of weakly instrumented watersheds would benefit from the development of an iterative approach between field measurements and/or monitoring and modelling approach to incrementally improve water cycle quantification. Such an approach may provide early-stage estimations of water fluxes and highlight which data are missing for more accurate estimation. The development of such water cycle simulation is strongly needed in order to sustainably manage water resources and to quantify the impact of climate on water availability.

To achieve this, our study develops a methodology to structure the progressive characterization of a minimally instrumented hydrosystem. The objective is to obtain a first simulation usable for water use management while guiding future surveys and monitoring. This approach was designed for a previously hydrogeologically unknown mountainous watershed in West Java (Figure 1A). Located on the volcano chain crossing the island, this mountainous area is strategic for a wide variety of water use, such as for population supply, irrigation, and industry. However, these hydrosystems are considered hydrogeologically unknown because their geological structures, extensions and dynamics have not yet been determined, which is necessary to improve their management. The only hydrological data available in the watershed is the discharge of the artesian springs used to supply the local population and/or industrial uses. Consequently, a better understanding of the hydrosystems and their functioning is necessary to define sustainable water management. Our study is divided into three steps: (i) hydrosystem characterization through a multidisciplinary study to set up a conceptual model; (ii) the transformation of the
conceptual model into a numerical one; and (iii) the
calibration of this latter using field data.

**2. Material and method**

**2.1. Methodological approach**

To set up the conceptual model, we need to de-
fine the structure, characterize the compartiments
of the water cycle and estimate their behaviors to
construct a meaningful model. First, the volcanic
structure and aquifer(s)/impervious bodie(s) geome-
tries have to be determined. For decades, geophys-
ical methods demonstrated their strength to image
volcanic structures and identify aquifer geometries
[Descloitres et al., 1997, Dumont et al., 2019, Gresse
et al., 2021, Vittecoq et al., 2019]. We thus imaged
the hydrogeological structure of the watershed us-
ing 2D electrical resistivity tomography (ERT) com-
bined with a geological review of surface outcrops
and of artesian borehole lithological data. Second,
the different aquifers need to be identified and char-
acterized between shallow, perched, and deep reser-
voirs [Join et al., 2016]. Geochemistry has proven
its efficiency in identifying the recharge areas and
the main characteristics of the aquifer [Bénard et al.,
2019, Join et al., 2016, Toulier et al., 2019]. To do so,
we used geochemical analyses and discharge mea-
surements to identify both perched and deep arte-
sian aquifers. This further allowed defining the ex-
tent of each identified aquifer type. In addition, the
elevation of aquifer recharge was assessed using wa-
ter stable isotope of rainfall within the watershed and
deep artesian aquifers. Lastly, hydrological monitor-
ing of river or spring discharges are an important data
to compute water balance between the different wa-
ter reservoirs of the watershed [Charlier et al., 2011,
Vittecoq et al., 2019]. To overcome the lack of such
dataset, a rough surface water balance was obtained
from point measurements of river flow, which are
used to estimate surface water/groundwater flows.

As the model relies on limited information while
integrating the whole water cycle, we used a lumped
model. The production function adopts the equa-
tions of Thornthwaite [1948], while the model
storage discharge equation is derived from Perrin
et al. [2003]. We applied a Bayesian approach, the
multiple-try differential evolution adaptive Metrop-
olis (MT-DREAM<sub>ZS</sub>) to account for the uncertainty in

model parameters calibration [Laloy et al., 2015]. The
starting probability distribution of the parameters
is set as uniform within a threshold estimated from
field measurements. Multiple runs of Markov chains
then allow transforming the a priori distribution
functions into a posteriori probability distribution
functions [de Marsily et al., 1999]. This probabilistic
exploration of the model parameters space, framed
by field interpretations, allows (i) to improve the
representability of the water cycle simulation, and
(ii) to guide future field characterization by estimat-
ing the influence of each parameter on the simula-
tion uncertainty. Thus future field surveys or instru-
mentation can be optimized using the simulation
guidelines.

2.2. *Study area*

The study site is located on the southern flanks of
the Salak volcano, which culminates at 2900 m a.s.l.
(Figure 1). This volcano is part of the mountain belt
that stretches along the center of the island of Java
from west to east. The study area straddles the sur-
face drainage limit between the northern and south-
ern sides of the island. Water resources of the water-
shed supply local populations as well as agricultural
and industrial activities. While the main agricultural
use consists of irrigating rice fields by river diversion,
other activities are supplied by spring captation or
pumping into wells. It is therefore essential to char-
acterize renewable water volumes between the dif-
ferent components of this hydrosystem to maintain
a base flow sufficient to sustain local biodiversity and
local communities.

2.3. *Field measurements*

2.3.1. *Geomorphology*

The geological description of the Salak volcano is
still limited. Therefore, we applied a geomorpholog-
ical analysis to define its main geological units. In
particular, slopes and land roughness allow to distin-
guish areas underlained by lava flows from those un-
derlained by detrital deposits, as described by Selles
et al. [2015]. Surface topography characterization
was performed on a digital elevation model (DEM)
calculated from the Advanced Land Observing Satel-
lite dataset [ALOS—Tadono et al., 2014]. With a res-
olution of 5 m, the ALOS DEM provides an accurate
depiction of the slopes morphology of Salak volcano. The DEM geodesic reference is EPSG 32748—WGS 84/UTM zone 48S. The DEM was used to extract slope maps and delineate surface watersheds using the QGIS software (v3.10.11—https://www.qgis.org/). A pit-filling algorithm (r.fill.dir from the GRASS v7.8.4 package) was used to remove local depressions and other artifacts. The filled DEM was then used to define the watersheds and the hydrographic network using the “r.watershed” function from the GRASS v7.8.4 package. The slope map was calculated using the slope function of the GDAL v3.1.4 package.

### 2.3.2. Geophysics

To confirm the geomorphological classification of land surfaces into geological units and to provide some insights on the underlying geology, geophysical and geological campaigns were conducted during two field surveys. The first, conducted in December 2020, focused initially on the immediate surroundings of the artesian springs (Figure 1), before proceeding upstream in order to explore the artesian structure extension. Artesian spring water is collected into shallow boreholes, which logging data provide deep geological information. This initial survey was meant to understand the geological parameters controlling the location of the artesian springs, and to collect the geological information at depth needed to assist geophysical imaging interpretations. A second survey was conducted in June 2021 in order to complete geophysical imaging, such as to obtain a nearly continuous longitudinal geophysical section from the artesian springs area up to the break-in-slope marking the transition from the core of the volcano to the surrounding detrital accumulation zone.

During these two surveys, a total of 12 ERT profiles were acquired using an IRIS Instruments Syscal Pro resistivity meter (http://www.iris-instruments.com/syscal-pro.html). All 2D lines were acquired using 64 electrodes with a 10 m spacing and some roll-along specified in Table 1. Acquisition sequences were defined with Wenner-Schlumberger set up, which provides the best trade-off between shallow resolution and depth of investigation [Dahlin and Zhou, 2004]. Each measurement was stacked 3 to 6 times with imposed electrical potential of 800 mV. The measurement time window was 500 ms.

After field acquisition, ERT and roll-along data were merged and processed using Prosys III software (v1.05, https://www.iris-instruments.com/fr/download.html#processing). In order to delete gap-fillers measurements, as well as measurements with negative apparent resistivity, low intensity and/or low electrical potential data were filtered out. Finally, data with poor standard deviation (above 10%) and insufficient stacking were also removed.

Data were then inverted using 11 norm in Res2driv software version 3.59.121. The result is presented after 3 iterations. Similar parameters were used to invert all ERT profiles. The width of model grid cells was set to half the electrode spacing. The thickness of the top layer corresponds to 120% of the cell width, and the layer thickness increases by a factor of 1%. The depth of the model is defined with a vertical extension of 5%. The damping factors used are: (i) 0.15 for the initial, (ii) 0.02 for the minimum, and (iii) 5 for the top layer. The flatness ratio is set at 1 to ensure the imaging of both vertical and horizontal structures. Finally, the depth of investigation is estimated using a sensitivity estimation of the inverted resistivity model [Marescot et al., 2003, Oldenburg and Li, 1999]. This approach uses two inversions with different initial models and calculates the difference between the two inverted results. Here, we considered that if that difference exceeds 10%, the inverted resistivity model is not valid. Inversions yielding differences between 5 and 10% are regarded as poorly reliable but still informative.

2D-inverted resistivity profiles were groundproofed near the artesian springs using available geological logs defined with the drilling cores. They allowed us to identify the geological formations corresponding to the imaged geoelectrical layers. The ERT profiles obtained during the second survey allow us extending these interpretations up the slope, combined with surface geological observations.

### 2.3.3. Hydrology and geochemistry

The study site is located in Java Island characterized by a tropical climate. The year is divided in two seasons: the wet season characterized by strong moonson rainfall (October–March); the dry season with less or no rainfall at all [April–September; Dumont et al., 2022]. The annual rainfall varies from 2 to 5 m/year, while the daily average temperature ranges from 18 to 25 °C.
<table><caption>**Table 1.** Summary of ERT acquisition parameters, which provides for each of the 12 profiles the total length, number of electrodes and the roll-along. No roll-along was implemented along profiles #8, #11 and #12. However, profile #8 is composed of two profiles with an overlap of 310 m, whereas profiles #11 and #12 do not overlap. The two processing columns present the number of data obtained after (without gap-fillers) and before processing, with indication of the percentage of data retained for the inversion. The last column summarizes the absolute error of the inverted model for each profiles. All the profiles are located in Figure 1</caption><thead><tr><th>ERT profile</th><th>Length (m)</th><th>Electrodes</th><th>Roll-along</th><th>After processing</th><th>Before processing</th><th>Absolute error (%)</th></tr></thead><tbody><tr><td>#1</td><td>630</td><td>64</td><td>—</td><td>1100</td><td>966 (88%)</td><td>6.20</td></tr><tr><td>#2</td><td>630</td><td>64</td><td>—</td><td>1100</td><td>1023 (93%)</td><td>10.45</td></tr><tr><td>#3</td><td>630</td><td>64</td><td>—</td><td>1010</td><td>821 (81%)</td><td>8.99</td></tr><tr><td>#4</td><td>950</td><td>64</td><td>32</td><td>1944</td><td>1592 (82%)</td><td>6.68</td></tr><tr><td>#5</td><td>950</td><td>64</td><td>32</td><td>1920</td><td>1720 (90%)</td><td>8.25</td></tr><tr><td>#6</td><td>950</td><td>64</td><td>32</td><td>1944</td><td>1716 (88%)</td><td>9.72</td></tr><tr><td>#7</td><td>950</td><td>64</td><td>32</td><td>1944</td><td>1698 (87%)</td><td>15.09</td></tr><tr><td>#8</td><td>950</td><td>2 × 64</td><td>—</td><td>2200</td><td>1700 (77%)</td><td>14.36</td></tr><tr><td>#9</td><td>950</td><td>64</td><td>32</td><td>1920</td><td>1749 (91%)</td><td>16.75</td></tr><tr><td>#10</td><td>1190</td><td>64</td><td>32 + 24</td><td>2520</td><td>2434 (97%)</td><td>13.25</td></tr><tr><td>#11</td><td>630 × 2</td><td>2 × 64</td><td>—</td><td>2110</td><td>2015 (95%)</td><td>7.96 &amp; 3.15</td></tr><tr><td>#12</td><td>630 × 2</td><td>2 × 64</td><td>—</td><td>2200</td><td>1944 (88%)</td><td>4.03 &amp; 4.78</td></tr></tbody></table>

Alongside the geophysical surveys, springs were inventoried. Their discharges and physico-chemical properties were measured in the field and water samples were taken for further analyses, in particular laboratory measurements of major elements and isotopic ratios. Major elements analysis has been proceeded by UnPad ITB laboratory using an ion chromatography instrument. The water isotopic ratios have been estimated by the BATAN laboratory (National Nuclear Energy Agency of Indonesia) using a liquid-water stable isotope analyzer Los Gatos Research DLT-100. Springs were classified according to discharge, major element composition, and ground-water electrical conductivity [Bénard et al., 2019, Charlier et al., 2011, Join et al., 1997, 2016]. Water stable isotopic ratios were used to define the average elevation of the recharge area using the local meteoric line, which we defined from analysis of weighted rainfall collected within the watershed as described in Toulier et al. [2019]. For the watershed, rainfall was collected at four different elevations during the hydrological cycle 2006–2007, both the amount of rainfall and its mean monthly isotopic composition were measured. To characterize the artesian aquifer

recharge, three campaigns of water sampling for water stable isotope analysis were carried out from 2006 to 2007. The variations in water stable isotope ratios over the period allowed us to assess the seasonal variability of the recharge processes and to obtain an estimate of the average elevation of the aquifer recharge. Thus the spread of isotopic results informs on changes in the dynamics of the recharge from dry to wet seasons.

Once the hydrogeological structure of the subsurface was defined, it was necessary to understand the distribution of water flow and its dynamics. This step was hampered by the absence of long-term monitoring of river discharges within the surface watershed. Besides, artesian aquifers extension is generally large and disconnected from the surface drainage network. On the one hand, part of the subsurface flow could not be quantified with hydrological monitoring. On the other hand, the model must represent the functioning of the artesian springs and not the whole system in order to be calibrated on the available spring discharges. The river flow measurements undertaken during the 2019 dry season and the 2020 wet season were used as a constraint for sur-
face discharge evaluation in the model. During these surveys, the discharge was repeatedly measured at downstream locations whereas the upstream rivers were measured only during the wet season. The specific discharge rates (as a function of surface area) derived from these site-specific measurements were used as proxies of the specific runoff of all upstream sectors (central and proximal) in order to compute the total runoff. The difference in flow between the dry and wet seasons provides some preliminary insights on annual variations in the dynamics of the rivers. Nevertheless, this information is not comprehensive enough to be used as input in the models and was therefore only used to constrain the simulation results.

Finally, to constrain the water cycle in the system, we use the discharge of the artesian springs tapped for the drinking water supply of the population (Figure 2). These springs have been gradually instrumented since 2015. It is thus possible to have a partial measurement of the artesian aquifer outlet. Indeed, these measurements do not integrate all the artesian springs identified on the field, nor the underground flows that may not be drained by the springs. Nevertheless, we summed the discharge of all the inventoried springs issuing from the artesian aquifer to obtain one discharge time series of more than 5 years, from 2015 to mid-2021, that recorded the seasonal fluctuations of the system. The climate forcing during this period evolved from one year to another with notably an ENSO (El-Niño Southern Oscillation) event between 2014 and 2016. These data allow us to understand how the system reacts to variable forcings.

## 2.4. Water cycle simulation

### 2.4.1. Lumped model structure

Since the final structure of the hydrological model depends directly on the characterization of the hydrosystems, this section will focus on the description of the equations used. The model was chosen to run at the daily timestep in order to account for the strong rainfall events that occur in such a tropical context.

The production function of the model consists of a soil production storage (SW) defined to estimate

the balance between rainfall, evaporation, and vegetation transpiration. The Thornthwaite equation is used to compute the effective rainfall (ER) as:

$$ER = R - AET + \frac{\Delta SW}{\Delta t} \quad (1)$$

where R is rainfall, AET is actual evapotranspiration, and SW is soil water content (SW₀ ≤ W ≤ SWmax), where SWmax represents the soil water holding capacity and SW₀ is the lower limit of evaporable soil water. The AET is calculated as follows:

$$AET_n = \begin{cases} PET_n & \text{for } P_n + SW_{n-1} \ge PET_n \\ P_n + SW_{n-1} & \text{for } PET_n - P_n \ge SW_{n-1} \\ 0 & \text{else.} \end{cases} \quad (2)$$

Where n is the daily timestep and PET is the potential evapotranspiration estimated using the Penman-Monteith equations [Monteith, 1965], which explicitly accounts for climatic and vegetation characteristics:

$$PET = \frac{0.408 \times \Delta(R_n - G) + \frac{\rho_a * C_p}{r_a} * (e_s - e_a)}{\Delta + \gamma * \left(1 + \frac{r_s}{r_a}\right)} \quad (3)$$

with Δ the slope of the saturation vapor pressure curve as a function of the air temperature (kPa·°C⁻¹); Rₙ the net radiation (MJ·m⁻²·day⁻¹); G the soil heat flux density (MJ·m⁻²·day⁻¹); ρₐ the air density; rₐ the aerodynamic resistance (s·m⁻¹); rₛ the stomatal resistance (s·m⁻¹); eₛ the saturation vapor pressure at the mean air temperature (kPa); eₐ the actual vapor pressure of the air (kPa); γ the psychometric constant equal to 0.066 kPa·°C⁻¹; Cₚ the specific heat of the air (MJ·kg⁻¹·°C⁻¹).

Weather data used for the water cycle simulation come from the Stamet Citeko station of BMKG network (Badan Meteorologi, Klimatologi, dan Geofisika—http://www.bmkg.go.id) located at 14 km north of the watershed. While the station is not within the watershed, the rainfall analysis of Dumont et al. [2022] demonstrated the representativeness and accuracy of this weather station monitored from 2000 to 2021 for the studied area (Figure 2). For the estimation of potential evapotranspiration, the Stamet Citeko station is the only station providing long-term climate data (temperature, radiation, wind, sun exposure time per day and relative air humidity) without long gaps. However, due to punctual gaps in the climate data, they have been averaged over the time series to obtain the climatological potential evapotranspiration for each day of the year
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Daily rainfall (R), potential evapotranspiration (PET) and artesian springs discharge (Q) time series from 2015 to mid-2021. Note that the higher noise of the data in 2015 and 2017 is a consequence of the pressure device; the noise is not representative of the real discharge of the spring.</figcaption></figure>

(Figure 2). This long-term climatic time series allows us to assign a warm-up period to the model from the beginning of 2000 to the end of 2014, the latter corresponding to the beginning of spring discharge measurements. This warm-up period ensures that the model during the studied period is not impacted by the choices of model initialization.

Below the SR reservoir, the ER feeds runoff (RO) and infiltration. Flux allocation is defined by the Baseflow Index (BFI) parameter, used to calculate the infiltration/effective rainfall ratio. Deeper in the ground, the model structure is not defined at this stage, as it depends on spatially varying hydrogeological characteristics. In the model, each reservoir is characterized by a recharge and by the infiltration or the discharge from a shallower reservoir. Its discharge is calculated as a function of the reservoir state using the GR4J model equations [Perrin et al., 2003]. For a reservoir level R (mm), the outflow $Q_R$ thus depends on its current filling level R and its maximum capacity $R_{max}$ (mm):

$$Q_R = R \left\{ 1 - \left[ 1 + \left( \frac{R}{R_{max}} \right)^4 \right]^{-\frac{1}{4}} \right\}. \quad (4)$$

### 2.4.2. Probabilistic parameter exploration

The objective of our methodology is to define the range of possible values for input model parameters, based on the interpretation of field data. This will allow a better linkage between the numerical simulations and field reality while improving the future characterization of the studied hydrosystems. For this reason, the parameterization of the artesian springs flow simulation is based on a probabilistic approach. To do so, we used the multiple-chains differential evolution adaptive Metropolis algorithm [MT-DREAM<sub>ZS</sub>—Laloy et al., 2015] to infer the posterior probability density functions (PDF) of model parameters. At the start, an a priori probability distribution is defined for each parameter. At each iteration, the MT-DREAM<sub>ZS</sub> algorithm uses a latin hypercube sampling to sample parameter values following the statistical methodology presented in Laloy et al. [2015]. Once the parameter values are accepted, the water cycle simulation is run in order to estimate a log-likelihood function used to compare simulation results and measured discharges. The use of multiple Markov chains allows us improving parameter space exploration, limiting the impact of local maximum likelihood on the search for best-fit parameters. In this study, five Markov chains are used with 50,000
iterations per chain. At each iteration, a candidate parameter set is tested via the Metropolis ratio. If it is accepted, the parameter set is validated and integrated into the calibration process.

**3. Results and discussion**

This section is divided into three parts. First, the conceptual model of the hydrological system is defined using field measurements and observations. Second, this conceptual model is used as the basis to build and set up the numerical lumped model which is used for the probabilistic parameters exploration. Third, the optimal parameter set is used to simulate the spring discharge.

**3.1. From observations to a conceptual model**

To understand the geological structure of the watershed, our methodology relies on a combination of geological, geomorphological and geophysical information (Figure 3). In addition, hydrogeological knowledge such as the location and properties of the springs are of prime importance to decipher the structure of the conceptual model.

**3.1.1. Geomorphological and geophysical results**

The slopes map allows dividing the watershed into three zones: the core of the volcano is characterized by steep slopes (average: 23.9°) and, on the opposite, the distal zone by lower slopes (average: 5.5°); the intermediate proximal zone shows variable slopes between 5° and 20° (Figures 3A and 4). In the latter, field geological observations showed that several morphological forms correspond to lava flow fronts, indicating that it is the meeting zone between lava flows tracking from the core of the volcano and the detrital deposits of the distal zone. Downstream of these massive flows, the valley is particularly eroded with slopes of up to 12°, which contrasts with the low erosion of the distal area. The outcrop of lava flow may have disconnected this area from the upper part of the watershed, limiting alluvium deposit or mark the upper limit of regressive erosion coming from downstream. This particular feature allows the underlying formations to outcrop.

The ERT profiles provide an imagery of the structures at depth (Figure 3B-D). The lowest profile

stretches from the artesian springs zone to the lower limit of the proximal zone. These profiles were geologically and hydrogeologically calibrated from the geological data of springs and boreholes located near the springs area (Figure 3B). Thus, these ERT profiles reveal a long, from 10 up to 50 m thick, superficial conductive layer (<60 Ω·m), which extends to about 570 m a.s.l. and thickens over a distance of 2400 m. This conductive layer corresponds to a clayey laharic formation identified in the wells. There it constitutes an aquiclude that seals the underlying resistive (>60 Ω·m) artesian aquifer. The artesian aquifer corresponds in the wells to mildly reworked volcanic breccias, pyroclastic rocks, and sandstones. This main deep reservoir feeds several high discharge (>10 L/s) springs observed in the springs area. Groundwater flows out from the aquifer through the locally thinner (notably in incised valleys) and/or locally more permeable laharic formations (Figure 4). Additional ERT profiles performed immediately downstream the spring area show that the aquifer is bounded there to the south-east by clayey conductive formations; this is an important feature as it surely prevents any groundwater to flow downstream the artesian springs. The western and eastern lateral ERT profiles confirm that both the aquifer and its confining layer extend laterally West and East. The deep resistive aquifer layer is less resistive in the western part, while it deepens towards the eastern limit. Unfortunately, these lateral limits of the aquifer were not reached by the current #5 and 8 ERT profiles (Figure 3C). However, given the absence of other known artesian springs laterally (see Figure 7 below), the aquifer is not expected to extend much farther than the chosen surface watershed boundaries.

At the transition between the distal and the proximal zones, a few tens of meters in elevation above the upper limit of the aquifer confining layer, the ERT profiles reveal a thick resistive layer (80–150 Ω·m) that reaches the surface on a steep slope (8–15°). Both the morphology and outcrops confirm a lava flow front issuing from the North. From this elevation, a conductive layer consistent with laharic outcrops is present in the areas with the shallow slopes (<8°). At the limit with the core of the volcano, the conductive layer deepens and becomes thinner (as seen in #9). This conductive layer is thus comprised between two resistive bodies (Figure 3D). The proximal zone appears as a transition area between a
<figure><img src="image_6.png" /></figure>

**Figure 3.** (A) Slopes map of the watershed extracted from the ALOS DEM with the location of the ERT profiles, (B) the downstream profiles starting from the artesian springs area (#4, #11 and #12—green), (C) lateral profiles oriented west–east (#5 and #8—red), and (D) upstream profiles (#9 and #10—blue). The distance along the profiles is presented in m (X axis), as well as the elevation (Y axis). Profiles not presented here are in orange in the map (A). In the slopes map, three geomorphological zones are delineated as: core of the volcano (slope > 15°), proximal (slope from 8° to 15°) and distal (slope < 8°) zones. The three field pictures illustrate the different types of vegetation cover and geological outcrops encountered in the catchment, with clayey rice fields (P1), laharic deposits of the proximal zone (P2), and lava formations of the core of the volcano (P3). In the lowest profile, two artesian wells are located with laharic formations in grey and pyroclastic ones in red. The two nearest springs are located with blue triangles.

reworked-pyroclastics dominant area (downstream) and a lava-flow dominant area (upstream). There, groundwater may flow from lava-flow to sandstone aquifers. There might also be surface-groundwater interactions, with both springs issuing from aquifers and, locally, streams infiltrating to aquifers. In the core of the volcano, only lava flows and breccias deposits crop out, which is consistent with alternating thinner conductive (lahars) and thicker resistive (lava flows) layers (Figure 3D).

### 3.1.2. *Geochemical and isotopic results*

Location and characterization of springs (discharge, geochemical and isotopic analyses) combined with hydrological measurements on

streams complete this geological/geomorphological/geophysical analysis. In the distal zone, there are no springs outside the spring zone (blue points in Figure 1). There, artesian springs are characterized by high discharges (between 10 and 40 L/s) and the highest measured electrical conductivities (130 to 150 µS/cm—Figure 5A). During the wet season, additional artesian springs were identified a little further upstream in the field. This highlights that the artesian aquifer discharge varies according to the hydraulic head in the aquifer which is controlled by seasonal recharge. When this latter increases, temporary springs appear at a higher elevation. Several of these springs only flow during the wet season, accommodating fluctuations in discharge of the
<figure><img src="image_7.png" /><figcaption>**Figure 4.** 2D conceptual hydrogeological model of the volcano slope. The water flows are represented by blue arrows. A zoom of the artesian springs area provides information on the artesian aquifer water level and the overflow process.</figcaption></figure>

aquifer throughout the year (Figure 4). Upstream, in the proximal zone, all springs have a lower discharge rate (≤5 L/s) and lower electrical conductivities (<120 µS/cm—Figure 5A) which suggests the existence of only local perched aquifers of small spatial extent (Figure 4).

Globally, the groundwater electrical conductivity is well correlated to the discharge of the springs as well as to the carbonate mass concentration, which may be related to the contact time between groundwater and the aquifer rocks (Figure 5A and B). The sodium/chloride (Na/Cl) ratio, defined by Join et al. [1997], tends to be well correlated with flow discharge (Figure 5C). Higher ratios (longer water/rock interactions) are usually reached in the deep artesian aquifer while the shallow aquifers exhibit lower Na/Cl ratios (Na/Cl < 1.6). The Na/Cl ratios support previous insights indicating that the proximal springs are fed by small extention aquifers while the artesian springs are supported by an extensive deep one (Figure 4).

Then, oxygen and hydrogen stable isotope ratios allowed to estimate the elevation of the recharge area of the artesian aquifer (Figure 6). The mea-

surement of monthly weigthed rainfall isotopic ratios allow to estimate the water local meteoric line of Salak volcano (SMWL; $\delta^2\text{H}_{\%0} = 7.83 \times \delta^{18}\text{O} + 12.35$). This meteoric line is consistent with Toulier et al. [2019] Bromo-Tenger meteoric line (BMWL; $\delta^2\text{H}_{\%0} = 7.71 \times \delta^{18}\text{O} + 12.37$). As for Toulier et al. [2019], the SMWL appears above the global and local meteoric line, GMWL and LMWL respectively, demonstrat- ing the local topographic effect on rainfall compo- sition. The oxygen and deuterium ratios measured in the artesian aquifer have a low range of variation throughout the sampling periods compared to the variation observed in other springs and rivers sam- pled in the same time in the watershed (Figure 6A). This stability confirms the important inertia of the artesian aquifer.

In addition, the monitoring of stable isotopes of rainwater sampled at four different elevations (475, 628, 760 and 923 m) during the 2006–2007 hydrologi- cal year was used to calibrate the local elevation me- teoric line: $\delta^{18}\text{O} = -0.0026 \times \text{elevation} - 4.9672$. We can thus estimate that the aquifer recharge occurs within the elevation range between 710 and 920 m a.s.l. according to Figure 6B, indicating a median
<figure><img src="image_8.png" /></figure>

**Figure 5.** (A) Spring discharge as a function of the groundwater (GW) electrical conductivity; (B) carbonate mass concentration as a function of the groundwater (GW) electrical conductivity; and (C) spring discharge as a function of sodium/chloride ratio (concentrations are transformed in meq/L to estimate this ratio). The dot colors correspond to groundwater electrical conductivity clusters defined in Figure 1. Blue and green dots correspond to artesian springs.

recharge area of about 800 m a.s.l. Considering that the artesian aquifer recharge area can hardly extend below about 600 m a.s.l. (see section above on the geophysical interpretations; 200 m below the 800 m a.s.l. median recharge area), we can conclude that rainfall precipitated up to about 1000 m a.s.l. may recharge the artesian aquifer (200 m above the 800 m a.s.l. median recharge area), provided effective rain-

fall does not vary too much within this elevation range [Toulier et al., 2019]. However, as recharge may also result from runoff after rainfall events, the upper limit of the recharge area is probably lower than 1000 m a.s.l. (Figure 4) as runoff may come from higher elevations than the recharge upper limit. More detailed research in that area should be performed to precise this upper limit.

### 3.1.3. Hydrological results

The final step in the characterization of the hydrosystem is the surface water dynamics. Concerning surface water, river discharge measurements showed limited flows in the upstream sectors with specific flows of 17 to 20 L/s/km² (river discharges between 120 and 200 L/s for a subcatchment surface of 6.5 to 11 km²). Downstream, the river flow increases up to 4140 L/s for a contributive surface of 40 km² (specific discharge of 104 L/s/km²). The overall watershed discharge evolves during the year with a decrease of the discharge by a factor of 8 during the dry season (a discharge of 510 L/s and a specific one of 12.75 L/s/km²). In fact, during the dry season, rainfall decreases drastically resulting in the drying up of perched springs and upstream rivers. The map of specific discharges is presented in Supplementary Material A.

The analysis of the artesian springs discharge time series allows a better understanding of the seasonal variations of hydraulic heads in the aquifer. These variations are buffered in the deep aquifer, with low water periods between October and March depending on the year. The annual variations range between 150 and 190 L/s. The suspected overflow effect described previously from the temporary discharge of the springs located above the springs area zone appears to be supported by a flow threshold around 190 L/s.

In this section, the combination of geological observations, geophysical imagery, slope analysis, hydrological and hydrogeological observations as well as geochemical and isotopic analyses allowed to characterize the hydrosystem's structure and limits, and provided a first robust hydrological conceptual model (Figure 4). The numerical model will then be based on this conceptual model. Nevertheless, the numerical model will also enable challenge and, if possible, confirm the coherence of this conceptual model.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (A) Oxygen and deuterium isotopic ratios of weighted rainfall (in light blue), artesian springs (in yellow, blue and green according to the sampling periods) and other springs and rivers of the watershed (in black). Meteoric lines are plotted with dashed lines: dark blue for the global meteoric water line [Craig, 1961]; light blue for the local meteoric line (GNIP station of Jakarta); orange for Bromo line [Toulier et al., 2019]; and the Salak line (inferred in this study). (B) Analysis of water stable isotopes to estimate the mean elevation of the recharge of the artesian aquifer.</figcaption></figure>

### 3.1.4. Water cycle characterization summary

In order to assist in understanding the transition from the conceptual to the numerical model, the characteristics of the hydrosystem are summarized below. In the studied area, a thick clayey layer extends over the distal area. It thins out at the center of the distal zone of the watershed, at places where the slope increases locally. At depth, the underlying aquifer is confined below this impervious layer over most of the area. Locally, thanks to the thinnings of this clayey layer, it outflows through several artesian springs (Figure 7). Outside the artesian spring area (blue circle in Figure 7), no significant springs were identified in this distal zone. The layering imaged at low elevation in the distal zone appears to change in the proximal zone farther upslope. One of the reasons is the presence of lava flow fronts in the steep areas of the proximal zone, which could act as local, perched aquifers. These perched aquifers feed springs with different geochemical signatures (Figure 5), characterized by a lower mineralization, and shorter residence time than the deep aquifer, typical of local, perched aquifers [Join et al., 2016]. Detritic and laharic accumulations probably produce the conductive clayey layers imaged near the surface. Farther upstream, slopes are directly underlying by high permeability pyroclastic flows with no interstratified low-permeability detrital formations.

There, the increase in slope, and progressive disappearance of clayey deposits blanketing the land surface, agree with the average recharge elevation estimated from stable isotopes (Figure 7). However, this result does not prevent the existence of a recharge from the core volcano slopes and lower proximal zone.

## 3.2. Probabilistic lumped model

### 3.2.1. Lumped model structure

The structure of the hydrological model is then set up based on hydrosystem architecture described in Figure 8. The lumped model relies on the equations previously presented in Section 2.4.

At the top of the model, seasonal fluctuations in rainfall (R) and actual evapotranspiration (AET) forcings control the water level in the soil water (SW) reservoir (Equations (1)–(3)—Figure 8). When the SW reservoir is full (>SWmax), excess water is transformed into effective rainfall (ER). The ER is then divided between runoff (RO), which supplies stream water, and infiltration to the shallow perched aquifer reservoir (P). The parting of ER fluxes is defined by the BFI parameter. The discharge of P (Qp) is calculated using Equation (4) and divided into perched springs discharge (PSD) and the recharge of the artesian aquifer (A), according to the perched flow index
<figure><img src="image_10.png" /></figure>

**Figure 7.** Slopes of the watershed extracted from ALOS DEM with the delimitation of the core of the volcano, proximal and medial zones. The artesian springs located in the artesian spring area (blue circle) are characterized by high electrical conductivity (130–250 µS/cm) and perched springs by a lower electrical conductivity (60–120 µS/cm). The mean recharge area is delimited by the lower (710 m, orange line) and higher (920 m, green line) limits estimated from water stable isotopes.

ratio (PFI). This latter corresponds to the ratio be-
tween water feeding perched springs discharge, and

water reaching the artesian aquifer. The discharge of
reservoir A (Q_A) represents the artesian spring's out-
<figure><img src="image_11.png" /></figure>

**Figure 8.** Structure (A) and spatial organization (B) of the lumped model. In the structure, the parameters are represented in red, climatic forcing in green, model fluxes in light blue and model states in dark blue. In the map, the surface of the recharge area (S) of the model is represented in green. This surface may evolve according to the results of the isotopic ratios (±ΔS). The artesian (blue circle) and the overflow springs areas (yellow extension) are delineated according to the eroded area in the slopes map (Figure 7).

flow, which is simulated and compared to spring discharge measurements during the calibration process (blue circle—Figure 8B). However, considering that some springs only discharge during the wet season, we considered instead that, above a specific water level in reservoir A (threshold $A_{OvF}$), the discharge of the springs no longer increases. This parameter consists of a filling ratio of reservoir A. This can be explained by the drainage of groundwater towards

non-perennial springs (not instrumented) located at slightly higher elevations (yellow extension—Figure 8B). To represent that shortcut in the model, an overflow function ($OvF$) was implemented into the artesian reservoir A. The excess water is thus transformed into a non-perennial springs discharge ($Q_{OvF}$).

Consequently, the numerical model resulting from the proposed conceptual model is thus com-
posed of three reservoirs: *SW*, *P* and *A* (Figure 8A). Each of them possesses a maximum filling level that constitutes 3 calibrated parameters: *SW*<sub>max</sub>, *P*<sub>max</sub>, *A*<sub>max</sub>. The other four parameters of the model are indexes *BFI* and *PFI*, overflow threshold level *A*<sub>OuF</sub>, and the extension of the recharge area *S* that supplies the artesian springs (Figure 8A). As the discharge series we want to simulate correspond to only a part of the water flows in the watershed, the model extension should include only a part of the watershed. We delineate this feeding catchment *S* as the area enclosed within the elevation range of the recharge area, inferred from the analysis of the isotopic ratios of the artesian springs (Figure 6B). This surface can evolve towards a higher or lower altitude according to the hydrological cycles and its extension can also change.

### 3.2.2. Probabilistic exploration of lumped model parameters

The calibration of the model is conducted using the probabilistic approach MT-DREAM<sub>ZS</sub>. For this purpose, uniform a priori probability distributions are defined for each parameter. Field data are used to define the range of parameter values to improve the convergence of the simulations, prevent convergence toward local solutions, and improve the representativeness of spring discharge simulations.

The first parameter $SW_{max}$ represents the thickness of the soil reservoir that serves as an interface between the atmosphere, the vegetation and the deep reservoirs. In the watershed, soils are thin on lava flows but thicken substantially on detrital formations. The a priori distribution of soil thickness $SW$ was allowed to vary between 50 and 500 mm (of water). The ranges of values of reservoir filling parameters of the perched and deep reservoirs, $P_{max}$ and $A_{max}$, were set based on geophysical imaging. In the central and proximal zones, resistant horizons, interpreted as perched aquifers, have thicknesses ranging from 1 to 20 m. At depth, the artesian aquifer displays thicknesses ranging from 10 m up to 60 m. The size of the feeding catchment $S$ is about 7.2 km² from the area enclosed within the elevation of the recharge area, inferred from the analysis of the isotopic ratios. However, the recharge may take place over a larger area as demonstrated by Toulier et al. [2019]. Therefore, the contributing catchment can

expand up to 30 km² according to the surface watershed (Figure 8B). An in-situ observable range of values for the remaining three parameters can not be constrained with field data. We therefore allowed these parameters to vary over a very broad range. Indexes *BFI* and *PFI* thus vary from 0 to 1, and overflow $A_{OuF}$ varies between 10 and 50% of the *A* reservoir thickness.

Once the size of the parameters space was set, 50,000 MT-DREAM<sub>ZS</sub> iterations were run on 5 Markov chains (see Supplementary Material B). Over all these iterations, 1401 sets of parameters were validated by the Metropolis ratio. Each parameter quickly converged to a best-fit value. Parameters appear uncorrelated with each other, highlighting their independence (see the correlation matrix in Supplementary Materials C). Out of the 1401 simulations, the last 25% are considered as a posteriori simulation and are used for the analysis of the results. The posterior PDFs are thus reported on Figure 9.

The parameter exploration converges to a restricted range of posterior distributions. The best-fit values of $SW_{max}$, $P_{max}$ and $A_{max}$ are 471 mm, 1.01 and 22.83 m respectively, each one is close to the probability peak (Figure 9—top panel). Their posterior distribution at the 95% confidence level is restricted between 447 and 495 mm, 1 and 1.11 m, and 22.2 and 23.66 m, respectively. These results demonstrate fairly good confidence in the parameters estimates despite the low constraint on the deep reservoir. The area is also fairly well constrained, between 17.1 and 17.6 km², with an optimum at 17.3 km². The two indexes, *BFI* and *PFI*, are also focused around their optimal values 0.4 and 0.43 respectively (range from 0.38 to 0.42 and 0.41 to 0.44). Finally, the optimum of $A_{OuF}$ is located at 0.174 (range from 0.173 to 0.175). The good convergence of the Bayesian approach is remarkable on the simulated flow curve since the 95% prediction interval is narrow (Figure 9—bottom panel).

With a KGE score of 0.87, the simulation represents fairly well the dynamics of the artesian springs, especially at the beginning (2015–2017) and for the hydrological cycle (2019–2020). In between, the quality of the simulation is slightly altered. For example, in 2018, the simulation overestimates the discharge during the wet season and early 2019, the simulated discharge keeps decreasing while the observations increased. Also, the simulation underestimates the
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Probabilistic exploration of lumped model parameters. The red lines represent the seven probability density functions (PDFs) of the parameters. In each line, the posterior distribution of each parameter is defined by a dark square. The optimum parameter value is pointed with a blue diamond. The bottom part represents the artesian springs observed discharge with yellow dots, the optimum simulated discharge in blue with the 95% range parameters uncertainty. The 95% range parameter uncertainty represents the variability in simulated spring discharge as a function of the posterior distribution of parameters. It does not represent a comprehensive assessment of the model uncertainties. The parameter uncertainty is small since the posterior distribution of the parameters is narrow.</figcaption></figure>

recharge of the aquifer at the end of 2020 and the start of 2021. Although the simulation is less efficient during these periods, the model error does not exceed 6% of the observed discharge.

The last model validation is based on the use of river specific flows measured in the field. To do this, our model allows computing the discharge of the different water flows over the watershed (Figure 10). The simulated watershed runoff (RO) varies between 0 and 7 m³/s on average with larger floods up to 16 m³/s. The perched springs discharge (PSD) varies between 0.2 and 1.5 m³/s. The artesian springs discharge ($Q_A$) is completely buffered between 150 and 180 L/s.

From these simulations, we extracted the average values corresponding to the months of field measurements. In December 2020, the runoff mea-

sured in the field was 18 to 20 L/s/km². In comparison, the median specific flow of the model during December is 13 L/s/km². However, the simulated floods occurring during this period have a strong impact on the average specific flow (50 L/s/km²). Considering the dry/wet period comparison, while in the field the river flow is 8 times higher during the wet season, the model simulates an increase of a factor of 5.5. Nevertheless, the high-frequency variability of the river flow simulated during the wet season is directly related to rainfall dynamics. This high variability complicates the analysis of punctual river measurements. It seems therefore impossible to constrain the model without continuous monitoring of river flows. Herein lies the most important margin of progress to improve the simulations. The implementation of continuous monitoring of river flow, even limited in time,
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Simulation inputs and outputs simulated at daily basis: rainfall (R), actual evapotranspiration (AET), runoff (RO), perched spring discharge (PSD), and the artesian spring discharge (Qₐ). The climatic forcings are in mm/day while the surface and groundwater flows are in m³/s.</figcaption></figure>

will improve the robustness of the simulations. This monitoring should be made at the outlet of the surface watershed as well as upstream in the proximal zone to calibrate perched and the artesian aquifer water balance. To do so, the upstream river should be instrumented to calibrate BFI parameter, as well as the downstream one for PFI.

### 3.2.3. Interest of probabilist approach for andesitic hydrosystem simulation

In this study, we developed a hydrological model aiming to represent the water cycle on the flanks of an andesitic volcano. To do so, we followed a two-step framework including (i) the proposal of a model structure as simple as possible to represent the complex and heterogenous surface- and ground-water behaviors; (ii) a probabilistic approach aiming at determining model parameters using a wider diversity of punctual measurements data.

Regarding the first issue, we focused on lumped reservoir models to balance the number of parameters with available data. The flexible structure of such models can be easily adapted to the complex

hydrogeological characteristics of the environment [Dubois et al., 2020, Mazzilli et al., 2019]. Such a modelling approach was already implemented in a French West Indies watershed by Charlier et al. [2008] who set up a conceptual hydrogeological model considering both surface and groundwater flows. Our study goes further by integrating multidisciplinary analyses combining geophysical, geochemical, and isotopic measurements to further constrain the model structure. Both studies show the ability of parsimonious reservoir models to represent complex subsurface flows of minimally instrumented watersheds. In these conditions, we assume that this parsimonious approach of lumped parameters models is the most consistent approach to represent the complexity of the volcanic environment. Process-based models would provide a finer representation of physical processes but would require extensive additional data to set up and validate. Given the amount of available data and the large heterogeneities of the studied areas, the feasibility of such a modelling approach is still to be proven.

Regarding the second aspect, we developed a
probabilistic approach aiming at integrating punc-
tual measurements data from a wide variety of disci-
plines. Prior distributions of parameter values were
derived using geophysical imagery, geochemical
and isotopic analyses as well as expert knowledge,
which allows strengthening the confidence in model
simulations. In the future, other datasets could be
integrated to improve the calibration/validation
process as piezometric data [Charlier et al., 2008],
geophysical monitoring [Lesparre et al., 2020], or
groundwater dating [Kolbe et al., 2016]. The in-
tegration of a wider diversity of dataset in water
simulation increase the confidence in the pro-
cesses captured by the model while limiting equi-
finality issues during the calibration of the model
parameters.

Finally, the probabilistic approach provides pos-
terior distribution of parameter values and an esti-
mation of parameter uncertainty instead of a unique
set of parameters. This allows representing fairly the
parametric uncertainties of the model that are in-
herently large in poorly gauged environments. How-
ever, these parametric results have to be carefully in-
terpreted as (i) DREAM algorithm try to disentan-
gle error source but Bayes' law suffers from the in-
teraction between individual sources of error (input,
output, parameter, and model structure error); and
(ii) the DREAM algorithm tends to provide a tighter
posterior distribution of parameters than other in-
formal Bayesian approach [Vrugt et al., 2009]. Be-
sides, the analysis of posterior distribution is a first
step towards the identification of ways to improve the
model structure and parametrization by further field
measurement.

**4. Conclusion**

The objective of this work was to integrate the
multidisciplinary characterization of a previously
unknown watershed to build a representative hy-
drological model of the water cycle in an andesitic
volcanic aquifer. This characterization was needed
for a better management of the wide variety of water
uses for human needs within the watershed, while
supporting biodiversity. Unfortunately, this char-
acterization was hampered by a significant lack of
measurements. To overcome this difficulty, we de-
fined a parsimonious and realistic conceptual model

in order to understand the hydrological behavior in
function of the climate forcing.

The multidisciplinary field characterization pro-
vides a characterization of hydrosystem structure as
well as the identification of the aquifers distributed
over the watershed. Our study revealed an exten-
sive artesian aquifer at the foot of the andesite vol-
cano confined by detritic clayey deposits. When
the thickness of this latter decreases, due to erosion
and/or thinner deposits, the aquifer supports signif-
icant artesian springs (150–190 L/s). On the slopes
of the volcano, the alternation of lava flows and de-
tritic formations allows the development of perched
aquifers with small extension. While the recharge of
both types of aquifers occurs on the volcano slopes,
extension and outlet location of the aquifers induce
spatially different behaviors. This characterization
based on data analysis and field observations allowed
to set up a conceptual lumped model characterized
by various reservoirs representing closely the differ-
ent functions observed. This model is used to simu-
late the artesian springs discharge.

Due to data scarcity, we improved the hydrologi-
cal simulation by implementing a probabilistic esti-
mation of the model parameters. The field charac-
terization allowed to define a uniform a priori distri-
bution for some of the parameters. From these distri-
butions, the use of several Markov chains allowed to
test a wide range of parameter sets while limiting the
influence of local optima. The Bayesian simulation
allows providing a PDF to analyze the uncertainty of
parameter values. Finally, the uncertainty could be
represented on the simulated water flows to estimate
its impact.

This combination of a field-based approach and
the probabilistic integration of field measurements
allows reinforcing the simulation of ungauged wa-
tersheds. Despite a limited constraint on surface
river flow, this simulation allows today to better man-
age the groundwater resource in this strategic sec-
tor. In the future, new field acquisitions could be
guided by the simulation results. To go further, it
would be possible to integrate into the model other
reservoirs to represent the main anthropogenic activ-
ities that develop rapidly in the watershed. First, an
exchange function with urban areas could be added.
This function would withdraw groundwater from the
deep aquifer and discharge water to the surface. A
second reservoir representing the rice terraces would
allow the model to be refined by integrating the delayed effect of surface flows, the increase in evaporation, and the influence on the recharge of agricultural practices.

## Conflicts of interest

Authors have no conflict of interest to declare.

## Acknowledgements

This project has been carried out within the collaboration between Sorbonne University (METIS laboratory), Water Institute by Evian, Danone Aqua group and the Universitas Padjajaran (UNPAD). This research was funded by Sorbonne University and Danone Aqua, in the framework of Danone Aqua Waterstewardship acceleration plan on watershed preservation. The authors are grateful to Damien Jougnot (Sorbonne University—METIS laboratory) for his help in the implementation of the MT-DREAM<sub>ZS</sub> algorithm.

## Supplementary data

Supporting information for this article is available on the journal's website under https://doi.org/10.5802/crgeos.192 or from the author.

## References

Aldrian, E. and Djamil, Y. S. (2008). Spatio-temporal climatic change of rainfall in East Java Indonesia. *Int. J. Climatol.*, 28(4), 435–448.

Bénard, B., Famin, V., Agrinier, P., Aunay, B., Lebeau, G., Sanjuan, B., Vimeux, F., Bardoux, G., and Deza-yes, C. (2019). Origin and fate of hydrothermal fluids at Piton des Neiges volcano (Réunion Island): A geochemical and isotopic (O, H, C, Sr, Li, Cl) study of thermal springs. *J. Volcanol. Geotherm. Res.*, 392, article no. 106682.

Buwono, N. R., Risjani, Y., and Soegianto, A. (2021). Distribution of microplastic in relation to water quality parameters in the Brantas River, East Java, Indonesia. *Environ. Technol. Innov.*, 24, article no. 101915.

Cerbelaud, A., Lefèvre, J., Genthon, P., and Menkes, C. (2022). Assessment of the WRF-hydro uncoupled hydro-meteorological model on flashy watersheds of the Grande Terre tropical island of New Caledonia (South-West Pacific). *J. Hydrol. Reg. Stud.*, 40, article no. 101003.

Chang, S. W., Chung, I.-M., Kim, M.-G., and Yifru, B. A. (2020). Vulnerability assessment considering impact of future groundwater exploitation on coastal groundwater resources in northeastern Jeju Island, South Korea. *Environ. Earth Sci.*, 79(22), article no. 498.

Charlier, J.-B., Cattan, P., Moussa, R., and Voltz, M. (2008). Hydrological behaviour and modelling of a volcanic tropical cultivated catchment. *Hydrol. Process.*, 22(22), 4355–4370.

Charlier, J.-B., Lachassagne, P., Ladouche, B., Cattan, P., Moussa, R., and Voltz, M. (2011). Structure and hydrogeological functioning of an insular tropical humid andesitic volcanic watershed: A multi-disciplinary experimental approach. *J. Hydrol.*, 398(3–4), 155–170.

Craig, H. (1961). Isotopic variations in meteoric waters. *Science*, 133(3465), 1702–1703.

Dahlin, T. and Zhou, B. (2004). A numerical comparison of 2D resistivity imaging with 10 electrode arrays. *Geophys. Prospect.*, 52(5), 379–398.

de Marsily, G. (2020). Will we soon run out of water? *Ann. Nutrition Metabolism*, 76(1), 10–16.

de Marsily, G., Delhomme, J.-P., Delay, F., and Buoro, A. (1999). Regards sur 40 ans de problèmes inverses en hydrogéologie. *C. R. Acad. Sci. Ser. IIA - Earth Planet. Sci.*, 329(2), 73–87.

Descloitre, M., Ritz, M., Robineau, B., and Courteaud, M. (1997). Electrical structure beneath the eastern collapsed flank of Piton de la Fournaise volcano, Reunion Island: Implications for the quest for groundwater. *Water Resour. Res.*, 33(1), 13–19.

Dubois, E., Doummar, J., Pistre, S., and Larocque, M. (2020). Calibration of a semi-distributed lumped karst system model and analysis of its sensitivity to climate conditions: The example of the Qachquouch karst spring (Lebanon). *Hydrol. Earth Syst. Sci. Discussions*, 24, 1–25.

Dumont, M., Peltier, A., Robin, E., Reninger, P.-A., Barde-Cabusson, S., Finizola, A., and Ferrazzini, V. (2019). Imagery of internal structure and destabilization features of active volcano by 3D high res-
olution airborne electromagnetism. *Sci. Rep.*, 9(1),
1-11.

Dumont, M., Reninger, P. A., Aunay, B., Pryet, A., Jougnot, D., Join, J. L., Michon, L., and Martelet, G. (2021). Hydrogeophysical characterization in a volcanic context from local to regional scales combining airborne electromagnetism and magnetism. *Geophys. Res. Lett.*, 48(12), article no. e2020GL092000.

Dumont, M., Saadi, M., Oudin, L., Lachassagne, P., Nugraha, B., Fadillah, A., Bonjour, J.-L., Muhammad, A., Hendarmawan, Dörfliger, N., and Plagnes, V. (2022). Assessing rainfall global products reliability for water resource management in a tropical volcanic mountainous catchment. *J. Hydrol. Reg. Stud.*, 40, article no. 101037.

Gharari, S., Hrachowitz, M., Fenicia, F., Gao, H., and Savenije, H. H. G. (2014). Using expert knowledge to increase realism in environmental system models can dramatically reduce the need for calibration. *Hydrol. Earth Syst. Sci.*, 18(12), 4839–4859.

Gómez-Delgado, F., Roupsard, O., le Maire, G., Taugourdeau, S., Pérez, A., van Oijen, M., Vaast, P., Rapidel, B., Harmand, J. M., Voltz, M., Bonnefond, J. M., Imbach, P., and Moussa, R. (2011). Modelling the hydrological behaviour of a coffee agroforestry basin in Costa Rica. *Hydrol. Earth Syst. Sci.*, 15(1), 369–392.

Gresse, M., Uyeshima, M., Koyama, T., Hase, H., Aizawa, K., Yamaya, Y., Morita, Y., Weller, D., Rung-Arunwan, T., Kaneko, T., Sasai, Y., Zlotnicki, J., Ishido, T., Ueda, H., and Hata, M. (2021). Hydrothermal and magmatic system of a volcanic island inferred from magnetotellurics, seismicity, self-potential, and thermal image: An example of Miyakejima (Japan). *J. Geophys. Res. Solid Earth*, 126(6), article no. e2021JB022034.

He, B., Takase, K., and Wang, Y. (2008). Numerical simulation of groundwater flow for a coastal plain in Japan: Data collection and model calibration. *Environ. Geol.*, 55(8), 1745-1753.

Hrachowitz, M. and Clark, M. P. (2017). HESS Opinions: The complementary merits of competing modelling philosophies in hydrology. *Hydrol. Earth Syst. Sci.*, 21(8), 3953–3973.

Join, J. L., Coudray, J., and Longworth, K. (1997). Using principal components analysis and Na/Cl ratios to trace groundwater circulation in a volcanic

island: The example of Reunion. *J. Hydrol.*, 190(1-2), 1-18.

Join, J.-L., Folio, J.-L., Bourhane, A., and Comte, J.-C. (2016). Groundwater resources on active Basaltic Volcanoes: Conceptual models from La Réunion Island and Grande Comore. In Bachelery, P., Lenat, J.-F., Muro, A. D., and Michon, L., editors, *Active Volcanoes of the Southwest Indian Ocean*, pages 61–70. Springer, Berlin, Heidelberg.

Khasanah, N., Tanika, L., Pratama, L. D. Y., Leimona, B., Prasetiyo, E., Marulani, F., Hendriatna, A., Zulkarnain, M. T., Toulier, A., and van Noordwijk, M. (2021). Groundwater-extracting rice production in the Rejoso Water-Shed (Indonesia) reducing urban water availability: Characterisation and intervention priorities. *Land*, 10(6), article no. 586.

Kolbe, T., Marçais, J., Thomas, Z., Abbott, B. W., de Dreuzy, J.-R., Rousseau-Gueutin, P., Aquilina, L., Labasque, T., and Pinay, G. (2016). Coupling 3D groundwater modeling with CFC-based age dating to classify local groundwater circulation in an unconfined crystalline aquifer. *J. Hydrol.*, 543, 31–46.

Lachassagne, P., Aunay, B., Frissant, N., Guilbert, M., and Malard, A. (2014). High-resolution conceptual hydrogeological model of complex basaltic volcanic islands: A Mayotte, Comoros, case study. *Terra Nova*, 26(4), 307–321.

Laloy, E., Linde, N., Jacques, D., and Vrugt, J. A. (2015). Probabilistic inference of multi-Gaussian fields from indirect hydrological data using circulant embedding and dimensionality reduction. *Water Resour. Res.*, 51, 4224–4243.

Lesparre, N., Girard, J.-F., Jeannot, B., Weill, S., Dumont, M., Boucher, M., Viville, D., Pierret, M.-C., Legchenko, A., and Delay, F. (2020). Magnetic resonance sounding measurements as posterior information to condition hydrological model parameters: Application to a hard-rock headwater catchment. *J. Hydrol.*, 587, article no. 124941.

Lo, W., Purnomo, S. N., Sarah, D., Aghnia, S., and Haridini, P. (2021). Groundwater modelling in urban development to achieve sustainability of groundwater resources: A case study of Semarang City, Indonesia. *Water*, 13(10), article no. 1395.

Marescot, L., Loke, M. H., Chapellier, D., Delaloye, R., Lambiel, C., and Reynard, E. (2003). Assessing reliability of 2D resistivity imaging in mountain permafrost studies using the depth of investigation in-
dex method. *Near Surf. Geophys.*, 1(2), 57–67.

Mazzilli, N., Guinot, V., Jourde, H., Lecoq, N., Labat, D., Arfib, B., Baudement, C., Danquigny, C., Dal Soglio, L., and Bertin, D. (2019). KarstMod: A modelling platform for rainfall–discharge analysis and modelling dedicated to karst systems. *Environ. Model. Softw.*, 122, article no. 103927.

Monteith, J. L. (1965). Evaporation and environment. *Symp. Soc. Experiment. Biol.*, 19, 205–234.

Mottes, C., Lesueur-Jannoyer, M., Charlier, J.-B., Charles, C., Guéné, M., Le Bail, M., and Malézieux, E. (2015). Hydrological and pesticide transfer modeling in a tropical volcanic watershed with the WATP-PASS model. *J. Hydrol.*, 529, 909–927.

Oldenburg, D. W. and Li, Y. (1999). Estimating depth of investigation in dc resistivity and IP surveys. *Geophysics*, 64(2), 403–416.

Perrin, C., Michel, C., and Andréassian, V. (2003). Improvement of a parsimonious model for streamflow simulation. *J. Hydrol.*, 279(1), 275–289.

Satapona, A., Putra, D. P. E., and Hendrayana, H. (2018). Groundwater flow modeling in the Malioboro, Yogyakarta, Indonesia. *J. Appl. Geol.*, 3(1), 11–22.

Selles, A., Deffontaines, B., Hendrayana, H., and Violette, S. (2015). The eastern flank of the Merapi volcano (Central Java, Indonesia): Architecture and implications of volcaniclastic deposits. *J. Asian Earth Sci.*, 108, 33–47.

Suriadikusumah, A., Mulyani, O., Sudirja, R., Sofyan, E. T., Maulana, M. H. R., and Mulyono, A. (2021). Analysis of the water quality at Cipeusing river, Indonesia using the pollution index method. *Acta Ecol. Sin.*, 41(3), 177–182.

Tadono, T., Ishida, H., Oda, F., Naito, S., Minakawa, K.,

and Iwamoto, H. (2014). Precise global DEM generation by ALOS PRISM. *ISPRS Ann. Photogramm. Remote Sens. Spatial Inform. Sci.*, II4, 71–76.

Thornthwaite, C. W. (1948). An approach toward a rational classification of climate. *Geogr. Rev.*, 38(1), 55–94. JSTOR.

Toulier, A., Baud, B., de Montety, V., Lachassagne, P., Leonardi, V., Pistre, S., Dautria, J.-M., Hendrayana, H., Miftakhul Fajar, M. H., Satrya Muhammad, A., Beon, O., and Jourde, H. (2019). Multidisciplinary study with quantitative analysis of isotopic data for the assessment of recharge and functioning of volcanic aquifers: Case of Bromo-Tengger volcano, Indonesia. *J. Hydrol. Reg. Stud.*, 26, article no. 100634.

Toulier, A., Lachassagne, P., Hendrayana, H., Fadilah, A., and Jourde, H. (2022). A cost-effective device and methodology to compute aquifer transmissivity and piezometry from free-flowing artesian wells. *Hydrogeol. J.*, 30, 1917–1931.

Vittecoq, B., Deparis, J., Violette, S., Jaouën, T., and Lacquement, F. (2014). Influence of successive phases of volcanic construction and erosion on Mayotte Island’s hydrogeological functioning as determined from a helicopter-borne resistivity survey correlated with borehole geological and permeability data. *J. Hydrol.*, 509, 519–538.

Vittecoq, B., Reninger, P.-A., Lacquement, F., Martelet, G., and Violette, S. (2019). Hydrogeological conceptual model of andesitic watersheds revealed by high-resolution heliborne geophysics. *Hydrol. Earth Syst. Sci.*, 23(5), 2321–2338.

Vrugt, J. A., ter Braak, C. J. F., Gupta, H. V., and Robinson, B. A. (2009). Equifinality of Formal (DREAM) and Informal (GLUE) Bayesian Approaches in Hydrologic Modeling? *Stoch. Environ. Res. Risk Assess.*, 23(7), 1011–1026.