On understanding mountainous carbonate basins of the Mediterranean using parsimonious modeling solutions

The study aims to demonstrate that an effective solution can be implemented for modeling complex carbonate basins, in the situation of limited data availability. Considering the alternative modeling approaches under circumstances of data shortage is more significant knowing the vulnerability and effectiveness of these kinds of basins to drought and climate change conditions. In this regard, a hybrid approach that combines time series analysis and reservoir modeling is proposed to describe behavior in carbonate basins. Time series analysis estimates the contributing area and response time of the fractured carbonate system beyond the catchment's hydrographic boundaries. The results obtained align with previous literature-based field surveys. This information is then used to develop a conceptual reservoir system using the GEOframe modeling system. The model is validated using in situ discharge observations and Earth observations (EO) data on evapotranspiration and snow. Model reliability is assessed using traditional goodness of fit indicators, hydrological signatures, and a novel statistical method based on empirical conditional probability. This approach enables detailed analysis and investigation of water budget components in Mediterranean carbonate catchments, highlighting their response to significant precipitation deficits.

Overall, our results demonstrate that flows from carbonate rock areas outside the hydrographic boundaries significantly impact the water budget of the upper Nera River. The storage capacity of the carbonate basin plays a crucial role in sustaining river discharge during drought years. In a single dry year, meteorological drought is considerably attenuated, while in subsequent dry years, it is slightly intensified. Multi-year droughts result in slower recovery due to the time required for precipitation to replenish the depleted storage that supported river discharge in previous dry years. This unique behavior makes these basins particularly vulnerable to the more severe and frequent drought episodes expected under future climate change.

1 Introduction

Carbonate/karst landscapes represent approximately 7 %–12 % of the Earth's continental area, and they provide a significant challenge for hydrologists (Hartmann et al., 2014). Due to the capability of these landscapes to retain water for a longer period (i.e., long-term hydrological memory catchments), their storage plays an important role in the control of drought propagation and delayed hydrological recovery (Alvarez-Garreton et al., 2021).

Generally, a carbonate/karst landscape forms when the percolated precipitation dissolves the subterranean carbonate bedrock and creates extensive fissures, open fractures, conduits, and caves. This can result in a complex network of groundwater flowpaths occurring within the same or adjacent aquifers (Kiraly et al., 1995). To model these types of systems one powerful solution is to use distributed, process-based models (PB) (e.g., Hartmann et al., 2014), which are based on solvers for the groundwater partial differential equation. Yet, the main challenge of this kind of distributed model is that they require a large amount of hydrogeological data and extensive field analysis to set appropriate physical parameter values and correct boundary conditions. On top of that, high computational power is needed to run these models (Li et al., 2022).

An alternative to PB is black-box models based on machine learning (MLM) in which not all details about the structure of the aquifer and the hydrodynamics parameters are needed, e.g., Tapoglou et al. (2014) and Castilla-Rho et al. (2015). Although the implementation of MLM is easy, the model parameters do not have a physical meaning and are only indirectly related to the characteristics of the carbonate system (Zhou et al., 2019). Furthermore, MLM does not explicitly solve the water budget, and, thus, it is not possible to have information about the dynamics of all water budget components.

Hydrological dynamical systems (lumped models, HDSys) represent another type of model, based on a set of ordinary differential equations (ODEs) that conceptualize the entire carbonate system as a series of reservoirs (e.g., Bancheri et al., 2019; Hartmann et al., 2014; Rimmer and Hartmann, 2012; Butscher and Huggenberger, 2008; Tritz et al., 2011; Jukic and Denić-Jukić, 2009; Dubois et al., 2020a). Instead of explicitly considering spatial variables, HDSys specify the interconnection of fluxes between different reservoirs, which leads to reducing the computational complexity. However, HDSys still require the definition of model parameters, which typically rely on calibration and inverse modeling using monitored discharge data or other relevant data sources (Hartmann et al., 2014). Several studies have also explored modeling the fast and slow drainage from carbonate systems using tracer information (e.g., Rimmer and Hartmann, 2012; Dubois et al., 2020a). This approach involves introducing an artificial tracer into a sinkhole and then tracking the tracer's movement in the surrounding areas at different times (Hartmann et al., 2014; Zhang et al., 2021; Nanni et al., 2020). While this technique can be useful, it is time consuming and may not always be feasible due to accessibility issues.

These HDSys can be conjugated by techniques that rely on the correlation between precipitation, and discharge can provide valuable insights about the behavior of carbonate systems, particularly in situations where field information about water circulation is limited. It could also provide useful information in the situation of missing tracer test analysis. For example, Fiorillo and Doglioni (2010) used cross-correlation analysis to estimate the time that water requires to flow through fissured aquifers. Another useful method, borrowed from applied economics (Kristoufek, 2014, 2015), was employed by Giani et al. (2021) to estimate the basin response time of hydrographs to precipitation, with successful results. However, according to the authors' best knowledge, to date, this data analysis technique has not been applied to complex carbonate systems to determine their hydrological response to precipitation.

This study aims to address the following five research questions (RQs):

1. Can the complex response of carbonate catchments to precipitation be modeled with HDSys relying only upon streamflow and precipitation time series, aided by cross-correlation analysis?

2. What type of modeling solution is suitable for this task and is a parsimonious modeling approach appropriate?

3. Are the classic goodness of fit scores enough to evaluate the reliability of the models?

4. What is the impact of external contributing areas on streamflow in catchments with fractured carbonate rocks? To what extent does this contributing area affect the total streamflow from small headwater catchments to the main outlet?

5. What is the role of storage in sustaining streamflow during the years with significant precipitation deficit in these types of catchments?

We have examined the water budget of the Nera River basin, which is a significant tributary of the Tiber River, the second largest river in Italy. The Nera River basin contributes nearly 50 % of the total discharge of the Tiber River and is characterized by a significant portion of fissured and fractured carbonate rocks feeding the river discharge by releasing a large amount of groundwater into the river bed from streambed springs. Thus, this catchment is a good representative of the carbonate catchments for answering the RQs. Additionally, groundwater data shortage is a problem that is not unique to the Upper Nera River area, and the findings of this study could help inform water management and policy decisions in other carbonate basins as well. By providing a comprehensive analysis of the water cycle in this area, this study could also help identify potential sources of water stress and suggest strategies to mitigate them.

2 Study area and datasets

2.1 Study area

The Nera River is the largest tributary of the Tiber River and its sources are in the Sibylline Mountains in central Italy. It is 116 km long and flows almost entirely in a deep valley called Valnerina, through limestone formations that constitute huge aquifers that are drained by the river. The landscape is mainly hilly and mountainous and is almost totally covered by forests, with pastures at higher elevations (from 1200 to 2200 m a.s.l). The Upper Nera River basin up to the Visso River station (our study area) covers an area of around 110 km² and the elevation ranges between 570 and 2200 m a.s.l., with a mean basin slope of 25 %. The basin is characterized by a Mediterranean climate, with precipitation concentrated mostly in the autumn–spring period, when floods generally occur. The annual precipitation and the average temperature of the basin are around 1100 mm and 10 ^∘C, respectively.

The discharge of the Nera River is contributed to by a set of permanent streambed springs fed by large limestone aquifers, already studied in Boni et al. (1986), that give rise to complex groundwater–surface water interactions. Nanni et al. (2020) and Mastrorillo et al. (2019) showed how these aquifers extend beyond the limits of the river basin into the wide and complex hydrogeological boundary of the Sibylline Mountains. Mastrorillo et al. (2019) estimated that the total contributing area of the fractured carbonate system outside the hydrographic boundaries of the basin (our study area) is around 97 km². Fronzi et al. (2021) did several tracer tests and showed that the river is fed (from the southeast of the basin) by carbonate aquifers with an area almost four times larger than that enclosed by the river station of Castelsantangelo, located upstream of Visso (Fig. 1a). Similar findings were found for the Ussita River, the main tributary of the Nera River at Visso, which is characterized by a real contributing area almost twice that of its hydrographic boundaries (Mastrorillo et al., 2019).

2.2 Terrain data and ground meteorological network

The terrain data for the geomorphological analysis of the basin were provided by the Marche Region Authority. The Horton Machine toolbox Abera et al. (2017) was used to define the basin and the hydrographic boundaries shown in (Fig. 1a); the basin was then further subdivided into 47 hydrologic response units (HRUs) (Fig. 1b).

Figure 1 (a) The study area with the DEM of the basin, additional external area, and the locations of the hydrometric stations on the study area. (b) Delineation of the hydrologic response units (HRUs) in the hydrological basin, based on the generic method (shape and morphology of the basin). The background map was retrieved from Google Satellite Hybrid (© OpenStreetMap contributors 2022. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.).

In the study, we used the meteorological network of the area provided by the Marche Regional Authority from which we selected 32 precipitation gauges, 21 thermometers, and 3 hydrometric stations that are distributed throughout the basin (these data are provided in the supplemental material). The monitoring network provides 15 min data for which a quality check to remove anomalous values and a re-sampling to the hourly resolution were performed. Streamflow data were calculated by transforming water levels measured at the hydrometric stations via rating curves updated bi-yearly (Fig. 2). The stations used in the study are Visso and Castelsantangelo (CSA) on the Nera River, and Madonna dell'Uccelletto (MU) on the Ussita River. To prevent any confusion, from this point forward in the text, “MU” will be used to refer to the station, while “Ussita” will be used to denote the river basin.

Figure 2 The discharge time series at the outlet of (a) Castelsantangelo (CSA), (b) Madonna dell’Uccelletto (MU), and (c) Visso (outlet) stations. The discharge values are not available for the gray color periods and have been eliminated from the hydrological analysis. All the time series are hourly.

The Visso station is the main outlet section of the study area and its data, covering the period 2007–2021, were used for the hydrological analysis. CSA station, with data available in the period of 2010–2016, is located upstream of the basin and is affected by a significant proportion of groundwater discharge coming from the external carbonate area (Fronzi et al., 2021). It should be mentioned that the continental deposits preserve this carbonate aquifer from the direct dissolution processes limiting the mature karst development in the saturated zones (Petitta et al., 2022). So the basin is not considered a fully karst system. The two stations were affected by the seismic sequences of 2016–2018, and thus these data have been excluded from the analyses. The earthquake altered the groundwater contribution of the fractured system, determining an abrupt and sustained change in the river and spring discharges in several parts of the basin (Di Matteo et al., 2020, 2021a). The MU hydrometric station on the Ussita River is characterized by about 3 years of hourly data since November 2018.

Table 1 summarizes the total contributing area (from the basin and from outside the hydrographic boundaries of the basin) of the three hydrometric stations (based on the literature Sect. 2.1=). These basin areas were derived from the terrain analysis and literature and can be estimated to be equal to 207 km² for Visso, 85 km² for MU, and 87 km² for CSA (see Fig. 1a and Table 1).

Table 1 Areas obtained from the classical hydrographic catchment delineation and from the hydrogeological survey at the outlet Visso, Castelsantangelo, and MU hydrometric stations.

2.3 Remote sensing data

Remote sensing data of evapotranspiration (ET) and snow depth were also used to complement the validation. In particular, we used MODIS actual ET (Mu et al., 2013) and the Sentinel-1 snow depth (Lievens et al., 2019). The global MODIS ET dataset (MOD16A2/MYD16A2) provides evaporation from wet and moist soil, evaporation from rainwater intercepted by the canopy before it reaches the ground, and transpiration through stomata on plant leaves and stems with 500 m spatial and 8 d temporal resolutions (Mu et al., 2013). The dataset can be downloaded from (https://modis.gsfc.nasa.gov/data/dataprod/mod16.php, last access: 12 December 2023).

The Sentinel-1 snow depth retrieval algorithm is based on an empirical change-detection method applied to the measurements of the cross-polarization ratio (${\mathit{\sigma }}_{vh}^{\mathrm{0}}/{\mathit{\sigma }}_{vv}^{\mathrm{0}}$; in dB) (Lievens et al., 2019). The Sentinel-1 snow depth retrievals are available online at https://ees.kuleuven.be/project/c-snow (last access: 12 December 2023).

3 Characterization of the basin

Data availability plays an important role in selecting a suitable modeling approach for carbonate basins. Prior to determining the groundwater contribution to river discharge with respect to the other components of the water balance, we analyzed the precipitation and streamflow time series by focusing on the area derived from the terrain analysis (i.e., without external area contribution, Table 1). Figure 3a1 shows a comparison of the cumulative precipitation against cumulative river discharge observed at CSA. It can be seen that the cumulative river discharge is much higher than the precipitation. In particular, Fig. 3a2 shows that the runoff coefficient of the basin at CSA – obtained by dividing the discharge at CSA by the precipitation time series – ranges between 4 and 5. This means that, assuming as the null hypothesis that the mean precipitation is constant in the red-shaded area of Fig. 1, the external groundwater contributing to the CSA outlet is at least four or five times larger than the extension of the basin delineated by terrain analysis, which is in line with Fronzi et al. (2021). Figure 3b1 and b2 also show evidence of external groundwater contributions for MU, determining a runoff coefficient higher than 1. The relative contribution of external groundwater to river discharge tends to reduce by moving downstream, as manifested by the smaller runoff coefficient observed at Visso (Fig. 3c1 and c2). The area of CSA, Ussita, and Visso increases in order, with CSA being the smallest and Visso being the largest, encompassing both CSA and Ussita. Figure 3 illustrates that the contribution from the carbonate (red) catchment decreases with increasing catchment size, which is a reasonable expectation. For further clarification about the lower runoff coefficient of Visso, the reader is referred to Sect. 5.2 and Fig. 12. Note that the panels of Fig. 3a1, b1, and c1 are obtained using the rain gauges closest to the river stations; the plots related to other rain gauges located inside and outside the analyzed basins are available in the Supplement (Figs. S1 and S2).

Figure 3 The relationship between cumulative observed discharge and cumulative precipitation at three different locations: Castelsantangelo (CSA) (a), Madonna dell’Uccelletto (MU) (b), and Visso (c). For each panel, the two plots are (1) cumulative observed discharge vs. cumulative precipitation recorded at the closest station, the green line representing a 1:1 relationship between discharge and precipitation, and (2) the runoff coefficient time series computed by dividing the discharge by the precipitation time series recorded at different stations. The runoff coefficient varies considerably in the basin. For Castelsantangelo (CSA) (a) and Madonna dell’Uccelletto (MU) (b) the runoff coefficient is approximately 4 and 1.5, respectively, while at the outlet of the basin (Visso, c), the coefficient is around 1.

After understanding that the external contribution to the basin is significantly greater than that provided by terrain analysis, it becomes necessary to determine the time delay between precipitation input and groundwater release into the basin.

To accomplish this, the approach proposed by Giani et al. (2021), which characterizes the time lag at which precipitation and streamflow are better correlated, is employed. To the best of the authors' knowledge, this method has not previously been used in the literature to estimate the travel time of groundwater in reservoir-based models.

Figure 4 shows the application of the method of Giani et al. (2020) for streamflow at the CSA and MU stations and the precipitation over the basin to understand the time lag at which these two variables are most strongly correlated. Figure 4a shows that ∼ 30 d (700 h) is the time window at which the precipitation and discharge of CSA are most correlated, which is in line with the results of Nanni et al. (2020), who showed that the mean tracer transit time (days) is around 26($+/-\mathrm{3}$) d for CSA. We attribute the high correlation of the first time window (3 d/70 h) to the subsurface river discharge response and that of the second to the groundwater contribution to the total river discharge. Tracer tests conducted later by Fronzi et al. (2021) demonstrated that springs emerge along the Nera River Valley and feed the river directly, as well as some streambed springs that emerge a few kilometers downstream. For MU, Fig. 4b shows that 167 d (4000 h) is the window length in which the discharge and the precipitation are most correlated. Evidence of a delayed external groundwater contribution in the Ussita River discharge is also present in Mastrorillo et al. (2019), who showed that a fraction of the total external area contribution to the Nera River feeds the Ussita River basin before MU (see Table 1).

Figure 4 The correlation between precipitation time series and discharge for (a) Castelsantangelo (CSA) and (b) Madonna dell’Uccelletto (MU), using different time window lengths. For Castelsantangelo (CSA), the strongest correlation is related to the 3 d (around 70 h) and 30 d (700 h) time lags. The earlier correlation is likely linked to the surface basin's response to the precipitation, while the latter is associated with water recharge from fissured rocks. For Madonna dell’Uccelletto (MU), the highest correlation between precipitation and discharge is at 167 d (around 4000 h) time lags. The higher correlation is likely associated with the slow response of the aquifer which responds slowly to precipitation events.

4 Methods for modeling basins with external groundwater contribution

Generally speaking, the topography and the derived shape of the basin give a guide on how to separate the pathways of water and obtain the precipitation volumes for each sub-basin. This general approach is evidently not applicable in the case of the Upper Nera basin, as also evidenced by the values of the runoff coefficients in Sect. 3. However, according to the field survey and measurements conducted by Mastrorillo et al. (2009) and Nanni et al. (2020), the extension of the contributing basin area has been determined, but the groundwater pathways are still unknown, and this poses important challenges for models based on fully coupled 3D or 2D surface–groundwater contributions.

4.1 Model structure

To model the fissured systems we split the basin area into two main parts: the surface catchment (SC) and the “external aquifer” catchment (AC). As mentioned earlier, the SC has been partitioned into the 47 HRUs, while the ACs are considered as unique sub-basins where the water flows into the SC from the CSA and MU.

The snow contribution is taken into account for both SC and AC, as shown by the extended Petri net (EPN) in Fig. 5. In particular, the EPN shows how precipitation flux is partitioned into snowfall and rainfall, according to the air temperature (Formetta et al., 2014). Snow then melts, increasing the liquid water in snow, and thereafter it can either refreeze or flow and be transferred to the vegetation reservoir. Afterward, the input flux into the canopy reservoir either evaporates from the canopy or falls through to the soil (Fig. 5b). The part that falls to the soil is partitioned into three reservoirs, presented in Fig. 5c. For more details on the functionality of the reservoirs, the reader is referred to the main references such as Formetta et al. (2014); Bancheri et al. (2019). To compute actual evapotranspiration (ET), we rely upon the Prospero model (Bottazzi et al., 2021), which uses a sun-shade canopy model (de Pury and Farquhar, 1997) that solves the energy balance for both sun-lit and shaded vegetation, extending the one recently developed by Schymanski and Or (2017) to canopy level.

Figure 5 The expressions and extended Petri net (EPN) of different components associated with snow (a), canopy (b), runoff (c), and external groundwater runoff (d). For more details about the components of GEOframe-Newage, the reader is referred to Formetta et al. (2014) and Bancheri et al. (2019).

The idea behind creating a specific reservoir arrangement for the karst external catchment (AC) is that the fate of water in the soil is different in the surface catchment (SC) and the AC one. Therefore, a different arrangement of reservoirs has been included for CSA and MU (see Fig. 5d). Following the analysis done in Sect. 3, the AC water flowing into the CSA is conceived as a single reservoir with a travel-time parameter set to 30 d, while MU is modeled using a reservoir with a travel time equal to 167 d. In this study, the primary focus is on the temporal variation of the precipitation-recharge-discharge behavior of the AC water flowing from CSA and Ussita rather than the spatial variability of the carbonate system's behavior. This allows us to specifically investigate the impact of single or multiple-year drought events on the basin storage, as discussed in Sect. 5.2. So, a lumped modeling approach is more appropriate for providing more insights into the temporal system's response to drought conditions and its implications for basin storage. Furthermore, incorporating more spatial variability for the carbonate areas would result in an increased number of model parameters. This introduces additional uncertainty into the model. Given the limited availability of data for calibrating these parameters, using two separate lumped systems has been considered as an efficient strategy for the modeling. The results will also demonstrate that the temporal behavior of the AC water and its response to drought events could be investigated properly by this modeling approach.

4.2 Experiment setup, calibration, and validation

The model was calibrated with hourly river discharge, observed at the hydrometric stations of CSA, MU, and Visso, and cross-compared with EO products of evaporation and snow cover, as well as validated with river discharge observations. The comparison with EO data provides insights into the model's robustness in reproducing spatial patterns of ET and snow across the study area.

The model was spun up at the CSA sub-basin in the period January 2005–December 2009, while the period from January 2010 to December 2015 was used to calibrate the model. Also for Visso, the warm-up period involved the first 5 years of data (January 2005–December 2009), while the following years (January 2010 to December 2017) were used for calibration. Finally, the period from January 2019 to December 2021 was applied for the calibration of model parameters for MU. The model was validated using only data from the Visso station between January 2019 and December 2021 because CSA data were not available in this period, and MU data is quite short for calibrating the model. Note that the warm-up period was chosen based on the features of the sub-basins.

Based on the results of a sensitivity analysis (not shown), 18 parameters were chosen for the calibration process. The model calibration process was performed based on LUCA (Hay et al., 2006), using the available observed discharge at the mentioned stations. The calibration algorithm was used to maximize the Kling–Gupta Efficiency (KGE) (Gupta et al., 2009) value between the observed and simulated discharge time series.

In addition, different scores and hydrological signatures (Addor et al., 2017) were used to evaluate the robustness of the optimized model (Table 2), including mean daily discharge, high flow, low flow, flow duration curve slope, and low-flow duration frequency. In particular, high flow and low flow represent the 95th and 5th percentile of the discharge, respectively. The low-flow duration frequency signature is defined as the frequency of the days with discharge lower than 20 % of the mean discharge value ($\mathrm{0.2}\stackrel{\mathrm{‾}}{Q}$). Since the Nera River discharge does not show significant variations, this signature is considered as the frequency of days with discharges lower than the mean discharge value ($\stackrel{\mathrm{‾}}{Q}$) instead of 20 % of the mean discharge ($\mathrm{0.2}\stackrel{\mathrm{‾}}{Q}$). On top of the classic goodness scores (i.e., NS and KGE) and the mentioned signatures, the goodness of discharge simulations was evaluated in a more statistically elaborated way. This was done to understand the predictive uncertainty of the model better and, thus, the information brought by the model based on the observed river discharge. In particular, we extracted the most likely expected discharge and the range of probable discharge values by using an empirical-based conditional probability approach as described in the next section.

4.3 Assessing the reliability of the model with empirical conditional probability (ECP)

According to the literature, predictive uncertainty is defined as the probability of real values of any variable of interest (i.e., the predictand) conditional on all available information up to the present (Todini, 2008; Krzysztofowicz, 1999). This information is provided by any deterministic model (i.e., the predictor). According to this definition, we propose a method that identifies the conditional probability (i.e., the probability density of the predictand conditional on the model simulation, the predictor) based on a non-parametric parsimonious method to estimate the posterior distribution of a predictand. This approach overcomes the challenges of most predictive uncertainty-based statistical methods that have to deal with Gaussianity assumption (which can be far from reality in many cases) and problems such as extrapolation to extreme values.

The process of computing the empirical conditional probability (ECP) involves the following steps:

• Combining the observed discharge and the corresponding simulated values into a single dataset.

• Grouping the dataset into n classes (bins) according to the simulated discharge values. The quantile-based discretization method has been applied for binning data into different classes.

• Computing the empirical cumulative distribution function (ECDF) for each class j using the formula:

  $$\begin{array}{}\text{(1)}& {\mathrm{ECDF}}_j\left(Q\right)={\displaystyle \frac{\mathrm{1}}{m_j}}\sum _{i=\mathrm{1}}^{m_j}{I}_{X_i<Q}.\end{array}$$

  Here, ECDF represents the empirical cumulative distribution function of the jth class, m_j is the number of measures in the group, X_i denotes the ith measure in the group, and

  $$\begin{array}{}\text{(2)}& I_{X_i<Q}=\left\{\begin{array}{cc}\mathrm{1}& \text{if }{X}_i<Q\\ \mathrm{0}& \text{otherwise}\end{array}\right..\end{array}$$

• Computing the features of empirical distribution function (i.e., mode, maximum, minimum, and mean of the discharge) for each class.

From the ECDF, we can derive the empirical distribution functions for different classes (e.g., the one shown in Fig. 6b and c), which are assigned to each time step. Figure 6 shows an example of implementing this method at two time steps where the observed discharge values are equal to 3.6 and 4.51 m³ s⁻¹. In Fig. 6, the discharge corresponding to the highest probability and the observation discharge are indicated by a green bar and a red line, respectively. The x axis shows the range of probable discharge at each time step. For instance, as shown in Fig. 6b, the variation range of probable discharge values is between 3.25 to 5.1 m³ s⁻¹. The difference between the observed and the highest probable discharge value can be considered as an estimate of the predictive model error. Figure 6b also shows that the observed discharge is well matched with the highest probable discharge, so the model performance is considered reliable. On the other hand, for time step t₂, as Fig. 6c presents, the observed discharge does not match with the highest probable discharge value, and the difference between observed discharge and highest probable discharge is a measure of model reliability.

The histograms obtained for different bins (e.g., Fig. 6b and c) are dedicated to the time steps visualized as green dots and gray-shaded areas in Figs. 7, 8, and 9. The green dots and gray area illustrated in the figures provide an indicator of the reliability of the simulations, according to previous simulated and observed data. In particular, the disparity between the measured and the mode values (green dots) can be considered as a measure of this reliability. The complete estimation procedure is thoroughly documented in a specific Notebook, accessible in the supplementary material. It is important to highlight that the number of classes (bins) has been carefully chosen to ensure a meaningful histogram for each discharge class. Even for the shortest available dataset (at the MU station, which encompasses approximately 26 000 hourly data points), a reasonable number of samples are available for each discharge class.

Figure 6 (a) Two empirical probability distribution functions (EPDF), represented as histograms, which are conditional on two different values of simulated discharge. (b) The examples of computing the predictive error at time step t₁ with the observed discharge Q=3.61 m³ s⁻¹ and (c) at time step t₂ with the observed discharge Q=4.51 m³ s⁻¹ in Visso basin. The green bar in each histogram represents the mode of the EPDF, while the observed discharge is represented by a red line. The difference between the observed and the most probable discharge is a proxy for the predictive model error.

Figure 7 (a) Simulated discharge at the outlet of Castelsantangelo (CSA) during the calibration period. The part of discharge with low quality was ignored in the analysis (light gray-shaded area). The variation range of probable discharge values obtained by ECP is shown in the gray-shaded area. The green dots show the value of the discharge with the highest probability in each time step. The closer to the observed discharge, the more reliable model performance. The interval is not wide, and this is evidence of the capability of the model in the Castelsantangelo (CSA) discharge simulation. (b) The flow duration curve for the observed and modeled Castelsantangelo (CSA) discharge are in blue and red, respectively.

Table 2 Classic signatures applied to evaluate the performance of the model

5 Results

In this section, the evidence needed to answer the research questions (RQ) presented in the Introduction is provided. In Sect. 5.1, it is investigated whether the parsimonious modeling solution is efficient to model the water budget of complex systems (e.g., the Nera river basin with a huge external groundwater contribution). Moreover, different approaches to evaluate the model's reliability are investigated and, eventually, a new statistical approach is proposed to understand the model's uncertainty and reliability. Section 5.2 is concerned with the variability of different water budget components and the capability of the basin to sustain streamflow during single and multi-year drought episodes.

5.1 Model suitability to simulate river discharge (RQ1, RQ2, RQ3)

Figure 7a shows the simulated hourly discharge at the outlet of CSA, considering the external reservoir in the model with its travel time parameter set to 30 d, as found in Sect. 3 and described in Sect. 4.1. During the calibration period, the KGE (see Table 2) and R² were obtained equal to 0.51 and 0.62, respectively. Furthermore, looking at Fig. 7b and the value of the slope of the flow duration curve for CSA in Table 3, it seems that there is no considerable discrepancy between the high (95th percentile) and low (5th percentile) streamflow in CSA. Figure 7a also shows predictive uncertainty results obtained with the ECP method. In this figure, red dots and the blue line represent simulated and observed discharge, respectively. Due to the rapid frequency of oscillation, the discharge appears as clouds in the plot. Green dots represent the highest probable discharge based on ECP analysis. Higher differences between the green dots and blue line reveal the lower capacity of the model to predict discharge. Generally, the relatively small difference between the observed and the highest probable discharge values in each time step suggests the adequacy of the model to reproduce river discharge for this river section. In both Autumn 2010 and Autumn 2012, a notable disparity was observed between the simulated discharge (represented in red) and the actual measured discharge (depicted in blue). The pattern of the modal discharge (green dots) is also similar to the simulated one. This suggests that these two particular seasons deviated from the norm in comparison to other periods. Therefore, it would likely be beneficial to separately conduct more detailed studies for these anomalous seasons.

Table 3 Different signatures, after Addor et al. (2017), of simulated and observed flows during the calibration period at CSA, Visso, and MU hydrometric stations. It should be noted that the calibration periods are different for the three hydrometric stations and that the average daily discharges cannot be compared with each other.

Further information is derived from Table 3. It shows that for the majority of the days of the year, the discharge values are lower than the average river flow, indicating that the distribution of discharges is left-skewed. All the signatures in the Table are well reproduced, even though the calibration was not targeted at them directly. High-flow statistics show a discrepancy of −5 %, while low-flow statistics bias is greater with a difference of around −10 %. The duration curve of the simulated behavior is 30 % steeper, meaning that the actual discharges are, on average, greater than the simulated ones, even though the latter have higher extremes. The discrepancy between simulated and observed mean discharge is, however, limited to less than 3 %, which can arguably be considered below the heuristically expected uncertainty of the forecast. It should be noted that due to different periods of calibration at the three hydrometric stations, the average daily discharges cannot be compared to each other.

Figure 8a shows the simulated discharge at the outlet of the basin (the Visso River station). The model simulation here yields R² and KGE values of 0.77 and 0.83 for the calibration period and 0.77 and 0.87 for the validation period, respectively, which is a sensible improvement with respect to CSA. Considering Table 3, the frequency of low-flow (LFD) score for Visso, which is 205 d for the observed discharges, is similar to the observed value at CSA, demonstrating that low flows dominate the river even at the Visso outlet (except for the November 2013 event). However, at Visso, the model estimates more days with low flow than observed (211 vs. 205), while at CSA the opposite happens (180 vs. 208). In Table 3, a comparison of the slopes of the flow duration curves for Visso and CSA shows a lower oscillation of flow in CSA (see also Fig. 8b). Figure 8a also shows how the model performs statistically with regard to historically observed discharge. The model accurately reproduces the jump in discharge observed in November 2014, although the average discharge is slightly overestimated. Overall, the green, red, and blue lines exhibit close agreement.

Figure 8 (a) Simulated discharge at the Visso outlet during the calibration (2008–2017) and validation periods (2019–2021), depicted in the same plot. The discharge affected by the earthquake was ignored for analysis. The variation range of the probable observed discharge values in each time step is shown in the gray-shaded area. The green dots show the value of the discharge with the highest probability in each time step. (b) The flow duration curves for observed and modeled Visso discharge are shown in blue and red, respectively.

Considering the travel time parameter of the external reservoir equal to 167 d (see Sects. 3 and 4.1), the R² and KGE values of the model for simulating MU river discharge during the calibration period were equal to 0.71 and 0.68, respectively. Figure 9a–b also shows the simulated discharge at MU during the calibration period. Table 3 presents different scores related to the signatures, which show relatively good agreement. Although the long-term average value of discharge is almost the same at CSA and MU, the difference between low and high flows at MU is higher (Table 3). This confirms the dissimilar river regime behavior in these two parts of the Nera basin. Furthermore, the flow duration slope values in Table 3 demonstrate a higher variation of discharge at MU than CSA. Additionally, Fig. 9a shows the variation range of probable discharge at each time step (gray-shaded area) for MU, which seems to be generally narrower compared to those of Visso, suggesting that the model has more uncertainty to give information about the river discharge variability at Visso station. It is interesting to note that this comes with a calibration score for simulating Visso discharge better than those of MU, providing evidence about the information that ECP can give on top of classical metrics.

Figure 9 (a) Simulated discharges at the Madonna dell’Uccelletto (MU) outlet of the Ussita subbasin for the calibration period (2019–2021). The uncertainty analysis results obtained by the ECP method for the simulated discharge are shown by the gray area (b). The flow duration curve for observed and modeled discharge are shown in blue and red, respectively.

5.1.1 Evapotranspiration

As was mentioned in Sect. 4.2, ET was assessed by comparing GEOframe-Prospero ET (GET) and actual ET (MET) from the MODIS product (Mu et al., 2013), based on the Penman–Monteith relation and various measured optical quantities. Figure 10a1–a2 and a3–a4 compare the GET and MET for two CSA sub-basins characterized by different elevations of 696 and 1615 m above sea level, respectively. The two ET time series show good correlation, with correlation coefficients of 0.86 and 0.78, respectively. Furthermore, the figure shows MET with much higher values than GET (i.e., it is less sensitive to stresses than Prospero).

Likewise, the comparison of GET and MET for two sub-basins located at 778 and 824 m a.s.l is shown in Fig. 10b1–b2 and b3–b4, with the correlation values of 0.76 and 0.77, respectively. In terms of bias, some positive systematic differences were expected as the MET is partially based on large-scale forcings that embrace a much larger area at lower elevations than the study area. Comparing the plot of MET and GET reveals that higher ET values are more frequent in MET than in GET, resulting in a fatter tail for the MET than for the GET. The absolute difference between the two estimates can be used as an indicator of the error in these simulations.

The dynamics of ET is similar for both products but the peculiarity of the upper Nera catchment, located in a very mountainous basin, along with the uncertainty of the EO dataset makes the evaluation of the bias very uncertain. Overall, we are more confident in the estimates provided by Prospero, given that in a complex topography region like the one of the study area, the radiation component of the ET can play an important role in evaporative fluxes, and this is not explicitly considered in the MET product.

5.1.2 Snow

The current version of the GEOframe snow component provides the snow water equivalent (SWE) and not the snow depth. Thus, a comparison between the model results and the Sentinel-1 EO product was done only in terms of the spatial correlation between the snow cover obtained by the model and Sentinel-1. Figure 11a, and the associated box plot, shows the spatial correlation of Sentinel-1/GEOframe snow cover over time for CSA. To avoid uncertainties related to snow compaction (SWE vs. snow depth), only the period from December to March was considered in the analysis, which corresponds to the snow accumulation period. The figure shows the 75th and 25th percentiles of correlation are in the range of 0.79 and 0.5 for the CSA basin.

Figure 10 Comparison of actual evapotranspiration from MODIS and GEOframe-Prospero, and the associated scatter plots, for two Castelsantangelo sub-basins located at 696 m a.s.l (a1 and a2) and 1615 m a.s.l (a3 and a4), and two Nera sub-basins located at 778 m a.s.l (b1 and b2) and 827 m a.s.l (b3 and b4). The green and red lines in the scatter plots show the regression and 1:1 lines, respectively. At higher elevations, a larger discrepancy between the MODIS and GEOframe-Prospero actual ET is observed.

Figure 11 The spatial correlation of Sentinel-1 and GEOframe snow cover over time and the corresponding box plots for the (a) Castelsantangelo (CSA), (b) Nera, and (c) Ussita basins.

Figure 11b presents the spatial correlation between snow cover obtained by Sentinel-1 and GEOframe for the HRUs of the Nera River basin at Visso in the period of December to March varying quite widely. The 75th and 25th percentiles of correlation are in the range of 0.75 and 0.25, respectively. This agreement is relatively less for the Ussita sub-basin, as can be clearly seen in the same figure (Fig. 11c), but we do not have sufficient elements to explain the reason for this. Overall, although the agreement of Sentinel-1 and in situ snow depth has been investigated (Fig. S3); further investigations on the ground are necessary for snow modeling.

5.2 Water budget analysis and the effect of groundwater discharge on the river regime (RQ4 and RQ5)

In this section, we present the interannual variability of the different water budget components of the catchment to understand the response of Mediterranean carbonate systems to climate variability. In particular, we focus on the ability of these basins to sustain streamflow during periods of significant precipitation deficit, such as the one experienced in the region in 2012 (Di Matteo et al., 2021b). For each delineated sub-catchment, the budget can be expressed as

$$\begin{array}{}\text{(3)}& {{\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}}^{\mathrm{•}}=P+Q_{AC}^{\mathrm{•}}-E{T}^{\mathrm{•}}-Q^{\mathrm{•}},\end{array}$$

where • denotes quantities estimated by GEOframe, $\mathrm{d}S/\mathrm{d}t$ is the variation in water storage in the sub-catchment, P is the total precipitation in the sub-catchment, Q_AC is the external discharge supplied by fissured rocks, ET is the evapotranspiration, and Q is the discharge at the outlet. Note that the water budget is resolved over a hydrological year (from September of the previous year to October of the current year). In this section, all the simulated values of the different components of the budget during the period 2010–2021 are considered, even the years impacted by the earthquake which showed different behavior.

Figure 12 displays the different annual water budget components from 2010 to 2021, closed at CSA and Visso (outlet) stations. The carbonate areas and fissured rocks located upstream of CSA and MU are now called “external CSA” and “external Ussita”, respectively. In the figure, for each year the left-hand side bars are related to the input fluxes, including precipitation and the additional discharge supplied by the external CSA and external Ussita, the middle bars are associated with the output fluxes, which contain actual ET and river discharge at the outlet, and the right-hand side bar is the variation of the basin storage. Figure 12a (CSA) shows that the external CSA dominates the input fluxes and, consequently, the discharge at the outlet of CSA is significantly generated by this external area. Similarly, the external water from external Ussita and external CSA combined make a notable contribution to the input fluxes of the whole catchment (i.e., Visso station, Fig. 12b), almost equaling the average precipitation falling within the hydrographic boundaries of the catchment.

Regarding the dynamics of the external CSA (Fig. 12a), it can be seen that it is not constant, but rather it follows the interannual precipitation variability. For instance, during the dry year of 2012, the external CSA was significantly reduced with respect to the other years, while in 2014 the basins experienced the maximum annual precipitation for the observation period, resulting in the greatest increase of discharge from external CSA. Because of the low precipitation in 2012 CSA shows a slightly larger output flux than input flux, with a negative storage change, meaning that the basin was able to sustain river discharge during periods of significant precipitation deficit. For Visso (Fig. 12b), the storage differences remain positive for the period of interest, indicating a potentially infinite stored water accumulation over multiple years. Possible explanations are either that the basin feeds the groundwater system (not simulated by the model) between the CSA and Visso hydrometric stations (see also the observation discharge at CSA and Visso stations) or, considering the long-term memory of the basin, that 10 years (2010–2021) is a relatively short period for observing storage changes. While both assumptions remain unanswered, the complexity of the system makes the first assumption plausible; however, further investigations are needed to provide compelling evidence of this. Moreover, the first assumption could also justify the lower runoff coefficient at Visso station.

Overall, ET plays a minor role at CSA, but its contribution to Visso increases significantly. In a future drier and warmer Mediterranean, this could become a frequent situation, making this basin more vulnerable to drought episodes than others.

To better understand the catchment dynamics at Visso, in Fig. 13 we get rid of absolute values of the different water balance components and instead show the anomalies of precipitation, discharge, storage variation, and actual ET with respect to the average. Based on Fig. 13, it is clear that the pattern of anomalies for precipitation and storage are well in agreement. It can be seen that during dry years the storage anomalies are negative, meaning that the storage is able to sustain river discharge in the catchment during these years. In particular, during single dry years like 2017 and 2019, the meteorological drought was significantly attenuated (the precipitation deficit was larger than the river discharge deficit), while for subsequent dry years, like 2011 and 2012, the meteorological drought was slightly exacerbated (the precipitation deficit was smaller than the river discharge deficit). These results confirm the findings of Bruno et al. (2022), where a dataset of catchments in Italy indicated that carbonate basins are capable of attenuating meteorological droughts during single dry years. Also, the study of Alvarez-Garreton et al. (2021) shows that for basins characterized by large hydrological memory (i.e., large residence time of water within different stores, e.g., snow, GW, soil moisture) multiple dry years can result in an exacerbation of the meteorological drought. This is explained by the recovery from drought for these basins, which is slower given the time that precipitation takes to replenish the depleted storages that sustained discharge during the previous drought year. This can be observed by comparing the anomalies of discharge in 2013 and 2018 in Fig. 13. Regarding evapotranspiration, a relatively constant pattern is observed, with a tendency to increase during 2012-2013, likely due to higher evaporative demand sustained by water stored in the unsaturated zone in weathered bedrock (i.e., rock moisture, Rempe and Dietrich, 2018).

Figure 12 The annual water budget components during 2010–2021 over (a) Castelsantangelo (CSA), and (b) Visso (the whole catchment). The left-hand side bars are related to the input fluxes, i.e., precipitation and the additional discharge coming from the external Castelsantangelo (CSA) and external Ussita (purple bars), the middle bars are associated with the output fluxes, i.e., actual ET and river discharge, and the right-hand side bar is the variation of storage over each year.

Figure 13 The anomalies pattern of annual precipitation, discharge, actual ET, and storage variation for the catchment during the water years of 2010–2021 at Visso station.

6 Conclusions

Carbonate catchments supply a significant fraction of water for domestic water supply, energy production, agriculture, and industry, and they are strategic in dry climate areas like the Mediterranean region. Yet, the modeling of these complex systems is challenging because of the differences between the hydrogeological and hydrographic basin delineation. Regarding the malfunctionality of the conventional approach to delineate river basins, our findings demonstrated that it is still feasible to model complex carbonate rock catchments effectively by utilizing streamflow time series and additional information about the contributing area, specifically in case of lack of data. In particular, we leveraged the time series analysis to assist the hydrological modeling of the Upper Nera River basin – a complex fissured rock catchment located in the Apennines in central Italy. It is necessary to conduct a preliminary, yet straightforward, check on the water balance closure whenever doubts arise about the presence of external groundwater inputs in the river basin (see Sect. 3). Specifically, for this study area the runoff coefficients in the upper part of the basin range from 4 to 5, indicating a substantial contribution of runoff in that particular area (Figs. 3 and 4). A time series analysis has been applied to determine the response times of the external catchments which are aligned with the estimates derived from field surveys incorporating tracer tests. Thereafter, by incorporating additional groundwater reservoirs into the modeling solution, with their average response time obtained through rainfall-discharge time series analysis, a significant improvement in the model's performance was observed. This highlights the importance of considering the presence of groundwater reservoirs and their response times for accurately simulating the hydrological behavior of the catchments. The configurable structure of GEOframe played a crucial role in the investigations, as it enabled us to customize the modeling solutions according to the specific characteristics of the hydrological system in a straightforward manner. This flexibility allowed us to adapt the model to the unique features and complexities of the study area, ensuring a more accurate representation of the hydrological processes. Although difficulties still remain in river discharge estimation, we do not have sufficient elements to say whether these difficulties can be overcome with more complex and data-demanding modeling solutions. Overall, having more data with a longer period of overlapping records would probably be beneficial to improve the simulation of such a complex basin behavior. Although one of the limitations of our study is the limited number of stations with overlapping records, employing a physically based (albeit lumped) modeling approach together with a robust correlation analysis could mitigate the data shortage issue. The results also revealed that external groundwater discharge, primarily originating from fissured rocks, has a significant impact on the water budget of the basin, particularly in the upstream areas (CSA). Although this influence diminishes downstream, it remains a substantial component of the water budget, nearly equivalent to the average precipitation within the hydrographic boundaries of the basin.

To evaluate the reliability of simulated discharges, besides more traditional indicators, we employed a method based on the empirical probability of the observed discharge, conditioned on the simulated discharge. This methodology effectively assessed the model's performance. The relatively small variation range of probable discharge for CSA (the gray-shaded area) generally suggests the adequacy of the model to reproduce river discharge for this river section. The ECP analysis has provided compelling evidence that the discharge levels in Autumn 2010 and 2012 were anomalously low compared to the anticipated average. This finding strongly suggests that these seasons need further detailed investigation. Concurrently, these discrepancies underscore the uncertain inherent of parameter calibration. This remains true even for modeling procedures intended to be “physically based”. Despite the better classic scores for simulating Visso discharge than those of MU, the variation range of probable discharge values at each time step for Visso seems to be generally wider than that of MU, suggesting that the model has more uncertainty to give information about the river discharge variability at Visso station.

Additionally, the model performance is cross-validated using Earth observation ET and snow products. The results consistently showed a slight underestimation of ET obtained by the GEOframe compared to the MODIS ET product. However, assuming MODIS-based ET as a reliable reference would imply lower discharges or reduced water accumulation in the groundwater system, which is inconsistent with the budget analysis of the CSA. While this discrepancy cannot be dismissed for the entire Nera basin, the CSA budget suggests opposition. Moreover, the model snow simulations are not in reliable agreement with EO snow information, which warrants further investigation and analysis.

The water budget analysis of CSA shows that ET is not a significant component of the budget due to the substantial groundwater contribution from the carbonate area. The behavior of soil moisture/groundwater storage at the CSA station oscillates around zero, indicating a balance between recharge and discharge. The water budget analysis at the outlet of Visso shows a consistently positive accumulation of groundwater, which suggests the presence of a groundwater flow feeding by the river that is not adequately captured by the modeling solution. Therefore, further investigation is needed to better understand and incorporate these factors into the model to improve the representation of groundwater dynamics at the Visso closure.

We also examined the role of storage in sustaining river discharge during periods of significant precipitation deficit. The results revealed that during single dry years, such as 2017 or 2019, the anomaly of river discharge remained positive, indicating that groundwater storage played a crucial role in maintaining streamflow. However, during multi-year droughts, such as 2011–2012, slight drought conditions were observed.

Additionally, our research determined that periods of consecutive dry years led to a proportionally slower recovery of streamflow compared to single-year droughts as derived from the analysis of the discharge anomalies in 2013 and 2018. This discovery bears significant implications, especially in light of the Mediterranean region's trend towards increased aridity and warming, as underscored by preceding studies (Giorgi, 2006).