**Coupling SWAT Model and CMB Method for Modeling of High-Permeability Bedrock Basins Receiving Interbasin Groundwater Flow**

#### **Javier Senent-Aparicio 1,\*, Francisco J. Alcalá 2,3, Sitian Liu <sup>1</sup> and Patricia Jimeno-Sáez <sup>1</sup>**


Received: 10 February 2020; Accepted: 27 February 2020; Published: 29 February 2020

**Abstract:** This paper couples the Soil and Water Assessment Tool (SWAT) model and the chloride mass balance (CMB) method to improve the modeling of streamflow in high-permeability bedrock basins receiving interbasin groundwater flow (IGF). IGF refers to the naturally occurring groundwater flow beneath a topographic divide, which indicates that baseflow simulated by standard hydrological models may be substantially less than its actual magnitude. Identification and quantification of IGF is so difficult that most hydrological models use convenient simplifications to ignore it, leaving us with minimal knowledge of strategies to quantify it. The Castril River basin (CRB) was chosen to show this problematic and to propose the CMB method to assess the magnitude of the IGF contribution to baseflow. In this headwater area, which has null groundwater exploitation, the CMB method shows that yearly IGF hardly varies and represents about 51% of mean yearly baseflow. Based on this external IGF appraisal, simulated streamflow was corrected to obtain a reduction in the percent bias of the SWAT model, from 52.29 to 22.40. Corrected simulated streamflow was used during the SWAT model calibration and validation phases. The Nash–Sutcliffe Efficiency (NSE) coefficient and the logarithmic values of NSE (lnNSE) were used for overall SWAT model performance. For calibration and validation, monthly NSE was 0.77 and 0.80, respectively, whereas daily lnNSE was 0.81 and 0.64, respectively. This methodological framework, which includes initial system conceptualization and a new formulation, provides a reproducible way to deal with similar basins, the baseflow component of which is strongly determined by IGF.

**Keywords:** SWAT model; CMB method; interbasin groundwater flow; Castril River; baseflow filter

#### **1. Introduction**

Hydrological models have become essential tools for water management issues due to their ability to simulate the hydrological cycle through integrated and multidisciplinary approaches, along with their skills to simulate climate change scenarios, land use, and water management [1]. Reliability of such models depends on the spatial and temporal scales covered, as well as the capacity to conceptualize the system functioning [2–4]. Those hydrological models operating at a basin scale are powerful decision support tools because they can provide insights into water resource management [5]. Among them, the SWAT (Soil and Water Assessment Tool) model, a physically-based and semi-distributed eco-hydrological open access model [6] stands out. SWAT can simulate the quality and quantity of surface water and groundwater balance components at different catchment scales to predict the impact

of climate change on the water balance of large watersheds [7,8] and deduce the effect of human-induced actions on water resources, such as irrigation practices and land-use changes [9]. A recent review of water quality and erosion models reveals that SWAT is, by far, the most used model [10]. However, a downside of SWAT is related to the simplified groundwater concept [11]. The simplified representation of the groundwater discharge and aquifer storage processes has been highlighted by several authors as something that may lead to a misunderstanding of the hydrological processes that occur in groundwater dominated watersheds [12].

One of these limitations is SWAT's inability to consider interbasin groundwater flow (IGF). IGF can be defined as the naturally occurring groundwater flow beneath the topographic divide that defines the basin boundary introduced in the SWAT model and in other hydrological models. It contributes to the baseflow of another basin different from that from which it was generated [13]. The magnitude of IGF may be especially relevant in high-permeability bedrock areas, such as steep karst areas. IGF maintains permanent streamflow in dry seasons, thus significantly altering the water balance of a region [14]. Despite the fact that IGF is a common hydrological process in large karst areas, it is often difficult to estimate, even tentatively [15]. Methodologies to identify and quantify IGF have traditionally relied on physical techniques, including the soil–water budget and water fluctuation, when there are sufficient data [16,17], tracer techniques measuring mostly environmental chemicals and stable isotope contents of precipitation and stream water [18,19], and groundwater modeling tools for indirect evaluations [20,21]. Some studies have aimed to assess IGF using the SWAT model. More specifically, Palanisamy and Workman (2014) [22] developed the KarstSWAT model to simulate IGF in watersheds dominated by typical karst features, determining input (recharge in sinkholes) and output (discharge in springs) water component dynamics. Malagó et al. (2016) [23] developed the KSWAT model, which was based on a combination of two previous SWAT applications: (1) a SWAT model adaptation to consider fast infiltration through caves and sinkholes up to the deep aquifers developed by Baffaut and Benson (2009) [24] and (2) a karst flow model in Excel to simulate spring flow discharge developed by Nikolaidis et al. (2013) [25]. More recently, Nguyen et al. (2020) [14] proposed a two linear reservoir model to represent the duality of aquifer recharge and discharge processes in a karst-dominated area in Germany. However, this interest is ongoing and, to our knowledge, the SWAT model has yet to be combined with the CMB method to improve hydrological cycle simulation in those basins where there is a difference between groundwater flow divides and surface topographic divides.

The evaluation of IGF is a complex, uncertain task when groundwater system functioning is partially unknown and the spatiotemporal coverage of data is too low to implement suitable evaluation techniques. In general, spatiotemporal coverage of environmental variables decreases in mountainous areas, thereby limiting the range of suitable techniques to assess IGF and other water balance components. In ungauged areas, IGF can be indirectly assessed when enough is known about groundwater system functioning to assert that IGF equals net aquifer recharge, which is the typical circumstance in most mountainous karst areas in a natural regime. In steep basins with gaining streams under a natural (undisturbed) regime, long-term net aquifer recharge (R) and discharge can be equated when groundwater abstraction, direct evapotranspiration from shallow aquifers, and underflow to deep aquifers are virtually null [26,27]. In such undisturbed hydrological functioning, net aquifer discharge equals the baseflow component of streamflow [28–31], and the problem shrinks to a matter of implementing suitable and viable techniques to determine R. Note that R is the infiltration amount that effectively contributes to the aquifer storage after some delay, smoothing out the variability inherent to precipitation events [32,33]. To assess renewable groundwater resources that finally reach streambeds, R is the governing variable [4,34].

Different methods can be used for R [35,36]. An independent, well-known method to determine R is the atmospheric Chloride Mass Balance (CMB) method [37–41]. The CMB method has been widely applied in different orographic, climatic, and geological contexts to yield mostly long-term (steady) R estimates when recharge water salinity can be attributed to the atmospheric salinity that reaches the water table. The CMB method was recently used to assess distributed mean R from precipitation and its uncertainty over continental Spain by verifying that the CMB variables were steady long-term [34,42]. This data availability was the reason the CMB method was chosen to assess IGF in areas with no human activities. In other territories, other techniques strictly intended for regional R can be selected for IGF evaluations when there are available data sets of similarly sufficient confidence.

This paper aims to evaluate the reliability of coupling the SWAT model and the CMB method to improve the modeling of streamflow in high-permeability bedrock basins receiving IGF. To that end, two main methodological steps are introduced. The first conceptualizes the hydrogeological functioning to confidently estimate IGF from existing CMB datasets for a control period and introduces a new formulation to generate a long-term baseflow series. The second step integrates corrected baseflow series into the SWAT model to improve the streamflow simulation. This methodology has been applied to the Castril River basin (CRB), which is an undisturbed, high-permeability bedrock area, characterized by the strong contribution of IGF to streamflow, as evidenced by a preliminary surface runoff coefficient greater than one.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The Castril River is an aquifer-fed mountain stream located 37◦47'–37◦59' north and 2◦40'–2◦50' west at the headwater of the Guadalquivir River watershed (GRW) in the province of Granada in southern Spain, adjacent to the Segura River watershed (SRW) (Figure 1a). The Castril River basin (CRB) headwater extends from the GRW-SRW divide (peak elevation is 2130 m a.s.l. in the north) to the Portillo Reservoir (outlet is 837 m a.s.l. in the south), covers an area of about 120 km2, and flows southward among the Sierra de Castril (west) and Sierra Seca (east) Mountains [43] (Figure 1b).

The climate is comparable to the continental Mediterranean, according to the Köppen classification [44]. Average annual precipitation is about 770 mm, with a coefficient of variation of 0.31 over the period 1951–2017. Precipitation is generated by Atlantic weather fronts coming in from the west and by short, intense Mediterranean convective storms. Most precipitation occurs during the autumn and spring. In winter, wet westerly and cold northerly winds predominate, whilst in summer and autumn, wet easterly and warm southerly winds blow [45]. Based on the period 1951–2016, the average annual temperature is about 8 ◦C, with the lowest temperatures in January and the highest in August. Average annual potential evapotranspiration is about 800 mm [46].

Geologically, the area belongs to the inner Prebetic domain of the external zone of the Alpine Betic Chain, which includes the following synthetic succession from bottom to top [47,48]: (1) Triassic gypsum-rich marls and clays (Keuper facies) with occasional limestones (Muschelkalk facies); (2) Lower and Middle Jurassic dolomites and oolitic limestones; (3) Upper Jurassic nodule limestones, calcareous marls, and marls; (4) Lower Cretaceous dolomites, dolomitic limestones, and marls; (5) Upper Cretaceous calcareous marls, marls, and dolomitic limestones; and (6) Paleocene to Middle Miocene limestones, marls, and calcarenites [48,49]. Small Late Quaternary alluvial deposits intermittently fill the valleys (Figure 1c).

From a hydrogeological point of view, geological materials can be classified into four groups, attending to the permeability type and storage capacity reported by the Geological Survey of Spain (IGME) (1988, 1995, 2001) [50–52]: (1) Triassic marls and clays are low-permeability materials that form the impervious boundary of local aquifers; (2) Jurassic and Cretaceous carbonate materials form highly permeable aquifers as thick as 300 m and have manifest karst features; (3) Jurassic and Cretaceous marls and calcareous marls are low-permeability materials, often confining the above Jurassic and Cretaceous carbonate materials; and (4) Late Quaternary alluvial deposits form temporary unconfined aquifers (Figure 1c).

**Figure 1.** (**a**) Location of the Castril River basin (CRB) within the Guadalquivir River watershed (GRW) in southern Spain, adjacent to the Segura River watershed (SRW). (**b**) Discretization of the CRB and 29 sub-basins using the 25 m resolution Digital Elevation Model (DEM) from the Spanish National Geographic Institute, showing other features cited in the text. (**c**) After the Geological Survey of Spain (IGME) (1988, 1995, 2001) [50–52] and direct field observations, a hydrogeological map of the CRB (scale 1:200,000), the schematic hydrogeological functioning of the CRB and the hydraulically connected upstream GRW and SRW contributing areas, a hydrogeological cross-section A–A showing aquifer dimensions, CBR location, groundwater divides and flow paths, and the 10 km × 10 km cells for distributed net aquifer recharge (R) in the part of continental Spain [34,42] covered in the study area was developed.

Hydrogeological functioning of the area was defined by IGME (1988, 1995, 2001) [50–52]. Aquifer boundaries, groundwater divides, and groundwater flow paths were established from hydrogeological maps, piezometry, and chemical and isotopic data. These hydrogeological criteria enabled experts to identify preferred areas for aquifer recharge in summits and for aquifer discharge, at the precise place where the incisive valley topography intersects the piezometry of Cretaceous carbonate aquifers to generate intermittent (upstream) and permanent (downstream) springs (Figure 1c). Downstream, outside the study area, Pliocene and Quaternary alluvial formations fill the Castril River valley and form an unconfined aquifer that hydraulically connects to the stream [43].

The study area is within the Sierra de Castril Natural Park, which has been an environmentally protected space since 1989, and is catalogued as a zone of special conservation for wildlife by the European Natura 2020 network. With respect to land use, forest, grassland, woodland and shrubland, sparsely vegetated areas, bare areas, and marginal rain-fed crops occupy most of the basin surface (Figure 2). Other marginal uses are seasonal livestock (sheep and goats), irrigated traditional crops, and riverine forest of *Pinus nigra* and *Pinus Halepensis* [53]. Neither permanent human settlements nor relevant water uses exist.

**Figure 2.** Land-use map (scale 1:25,000) from the Andalusian Environmental Information Network (REDIAM).

#### *2.2. Overall Model Description*

A coupled application of the SWAT model and the CMB method to improve streamflow simulation by considering IGF is introduced. This application includes four steps, shown as bulleted lists (Figure 3). The first step uses the CMB method to assess the magnitude of IGF contributing to the CRB baseflow from another upstream GRW and SRW areas, as described in Section 2.3. The second step uses the automated digital filter program BFLOW to split daily streamflow records into baseflow and surface runoff components, as a prerequisite to correct streamflow records by adding IGF to the baseline component, as described in Section 2.4. The third step uses the SWAT model to compare simulated streamflow with and without IGF, as described in Section 2.5. The fourth step uses the SWAT model for standard calibration and validation of simulated streamflow by considering the IGF correction, as described in Section 2.5.

**Figure 3.** Flow diagram for the coupled Soil and Water Assessment Tool (SWAT) model and chloride mass balance (CMB) method application to model streamflow of hydrological basins subjected to interbasin groundwater flow (IGF).

#### *2.3. CMB Method*

#### 2.3.1. CMB Method Application for Aquifer Recharge over Continental Spain

The atmospheric chloride mass balance (CMB) is one of the most widely used methods to estimate net aquifer recharge (R) from precipitation in different orographic, climatic, and geological contexts [35,36]. The CMB is a global method based on the principle of mass conservation of a conservative tracer, in this case the chloride ion, atmospherically contributing to the land surface. This technique yields mostly long-term (steady) R estimations when recharge water salinity can be attributed to the atmospheric salinity that reaches the water table [37–41].

The CMB method was recently applied to estimate distributed spatial mean R and its natural uncertainty (standard deviation) over continental Spain. For a confident application, the long-term steady condition of the CMB variables: atmospheric chloride bulk deposition, chloride export flux by surface runoff, and recharge water chloride content was verified [34,42,54]. This evaluation examined the influence of hydraulic properties (mostly permeability and storability) of different aquifer lithologies on R estimates, as well as the potential contribution of non-atmospheric sources of chloride [55]. For local usage, the reliability and hydrological meaning of distributed R were evaluated by comparing them with local, presumably trustworthy R estimates; one of these local cases was the CRB. Ordinary kriging was used to regionalize the CMB variables at the same 4976 nodes of a 10 km × 10 km grid. In each grid node a mean R value was estimated. Nodal R values were affected by two main types of uncertainty, the natural variability of the CMB variables and the error from its mapping. These uncertainties were identified and estimated [34,42].

The evaluation covered a 10-year period, which represented the critical balance period for the CMB variables to reach comparable steady means and standard deviations. This 10-year period matches the decadal global climatic cycles acting on the Iberian Peninsula, with irregular ~5-year positive and negative phases that follow the North Atlantic Oscillation trend [45,56]. Considering that (1) at least a 10-year balance period is required for reliably steady R evaluations in continental Spain and (2) the CMB datasets preferably spanned the period 1994–2007, the control period (1996–2005), which span a full 10-year long NAO climatic cycle, was chosen to estimate R in this work. Other authors have also implemented these CMB datasets for reliable local R evaluations in different climatic and geological settings, such as Andreu et al. (2011) [27] in Sierra de Gádor karst Massif in Southern Spain, Raposo et al. (2013) [57] in varied geological contexts in Galicia, on the coast of Northern Spain, and Barberá et al. (2018) [58] in the high-mountain, weathered-bedrock Bérchules basin in Southern Spain.

#### 2.3.2. IGF Series Generation

As introduced in Section 1, long-term steady R and IGF can be made equivalent in steep basins with gaining stream under natural (undisturbed) regime when groundwater abstractions, direct evapotranspiration from shallow aquifers, and underflow to deep aquifers are virtually null, and the hydrogeological functioning is well-defined [28–31]. This is the case of the CRB, as described in Section 2.1, because human water use is virtually null and there is enough hydrogeological information. The fraction of R produced in upstream contributing areas can be used as a reliable proxy for the additional baseflow fraction contributing to the CRB baseline [26,27]. To assess the IGF, R is the significant factor [4,34].

As described above, several cells, each one yielding a nodal mean R value and its standard deviation for the control period (1996–2005) instead of a yearly R series, cover the CRB and upstream contributing areas. Adapted from Pulido-Velázquez et al. (2018) [59], a procedure is introduced to obtain the yearly R series by adopting the temporal structure of the yearly P series for the control period (1996–2005). This model uses a correction function that forces the control R series to have the same relative deviation as the control P series, while maintaining the magnitude of its initial mean and standard deviation. The calibrated function is applied to obtain a yearly R series, presuming the correction function does not change. The calculation includes the following steps:

Average change in mean and standard deviation of P and R for the same control period (1996–2005):

$$\Delta \mathbf{m} = \frac{\mathbf{m}(\boldsymbol{\mathcal{R}}) - \mathbf{m}(\boldsymbol{P})}{\mathbf{m}(\boldsymbol{P})} \tag{1a}$$

$$
\Delta \sigma = \frac{\sigma(\mathbf{R}) - \sigma(P)}{\sigma(P)} \tag{1b}
$$

where Δm is the change in mean and Δσ is the change in standard deviation.

Normalization of the yearly *P* series:

$$P n\_i = \frac{P\_i - \overline{P}}{\sigma\_{\overline{R}}} \tag{2}$$

where *Pi* is *i*–year *P* and *Psi* is its normalized value, *P* is mean *P*, and σ*<sup>R</sup>* is standard deviation of mean *R*. Generation of yearly *R* series from yearly *P* series:

$$R\_{\bar{l}} = \mathbf{m}\_{\bar{\mathsf{C}}} + \sigma\_{\mathsf{C}^\*} \mathbf{P} \mathbf{s}\_{\bar{l}} \tag{3}$$

where *Ri* is *i*–year *R*, and m*<sup>C</sup>* and σ*<sup>C</sup>* are expressed as:

$$\mathbf{m}\_{\mathbb{C}} = \mathbf{m}(P) \cdot (1 + \Delta \mathbf{m}) \tag{4a}$$

$$
\sigma\_{\mathbb{C}} = \sigma(P) \cdot (1 + \Delta \sigma) \tag{4b}
$$

When Equation (4) is applied to the control yearly P series, the generated yearly R series adopts the same mean and standard deviation as the control R series from Alcalá and Custodio (2014, 2015) [34,42]. When this procedure is applied to the full yearly P series, the equivalent yearly R series is obtained by assuming that the bias correction remains invariant over the full observation period.

#### *2.4. BFLOW Program*

The automated digital filter program (BFLOW) to split daily streamflow records into the baseflow and surface runoff components was used. Nathan and McMahon (1990) [60] were the first to implement this recursive digital filter technique for baseflow analysis. The hypothesis of BFLOW is that low-frequency signals represent the baseflow component while high-frequency signals represent the runoff component [61]. This technique gives results similar to those obtained using other automated

models or manual techniques despite having no physical basis. BFLOW has been used in many studies related to the SWAT model [30,62]; see Arnold and Allen (1999) [63] for more details about this technique. The baseflow obtained by BFLOW was increased by adding IGF, calculated in previous section.

#### *2.5. SWAT Model*

#### 2.5.1. Description of the SWAT Model

The Soil and Water Assessment Tool (SWAT) is a long-term watershed hydrological model with strong physical mechanisms, developed jointly in 1994 by the Agricultural Research Service of the United States Department of Agriculture (USDA-ARS) and Texas A&M AgriLife Research, part of the Texas A&M University System. The SWAT model simulation includes atmospheric precipitation, surface runoff, subsurface flow, evapotranspiration, groundwater flow, river network flow concentration, and other intermediate water balance components subjected to variable delays [6]. The SWAT model first divides the study area into several hydrological response units (HRUs) based on the digital elevation model (DEM), land use, soil type, and meteorological data. Then, SWAT establishes a hydrophysical conceptual model of each HRU, calculates runoff in each HRU, and finally connects the entire set of the HRU runoff responses through the river network of the study area toward the basin outlet. The hydrological processes simulated by the SWAT model are based on the following water balance equation:

$$\text{SW}\_{\text{t}} = \text{SW}\_{0} + \sum\_{i=1}^{t} \left( \mathbf{P}\_{\text{day}} - \mathbf{E}\_{\text{day}} - \mathbf{Q}\_{\text{surf}} - \mathbf{W}\_{\text{secp}} - \mathbf{Q}\_{\text{gw}} \right)\_{\text{i}} \tag{5}$$

where SWt is final soil–water content (mm H2O), SW0 is initial soil–water content (mm H2O), *t* is time (days), Pday is precipitation on day *i* (mm H2O), ETday is evapotranspiration on day *i* (mm H2O), Qsurf is surface runoff on day *i* (mm H2O), Wseep is water amount that enters the vadose zone from the soil profile on day *i* (mm H2O), and Qgw is groundwater return flow on day *i* (mm H2O).

#### 2.5.2. Data, Model Set-Up, Calibration, and Validation

Datasets used to implement the SWAT model were: (1) the 25-m resolution DEM from the Spanish National Geographic Institute (IGN); (2) the land-use map (scale 1:25,000) from the Andalusian Environmental Information Network (REDIAM); (3) the 1-km resolution georeferenced soil data from the World Soil Coordination Map; (4) the 5-km resolution nodal daily precipitation series in Spain from the Spanish National Weather Service (AEMET) grid version 1.0, which cover the period 1951–2017; (5) the 10-km resolution nodal daily temperature series in Spain from the fifth version of the high-resolution SPAIN02 grid, which cover the period 1951–2016; and (6) the 24-h streamflow records downloaded from the Spanish Centre for Public Works Studies and Experimentation (CEDEX) website. The open source QGIS interface for SWAT (QSWAT 1.8) was used to set up the SWAT model.

The SUFI-2 algorithm of SWAT-CUP (Calibration and Uncertainty Programs) to calibrate and validate the SWAT model was used. Based on our previous modeling experiences [64,65], twenty-one widely used flow calibration parameters and their ranges were initially selected. Aimed at reaching an acceptable calibration, two iterations (representing 500 simulations each) were performed; the first included 13 parameters on a monthly scale, the latter included 8 parameters on a daily scale. To mitigate the effect of initial soil–water condition, a five-year warm-up period was imposed [66]. The periods 1995–1997 and 1982–1984 were, respectively, selected for the calibration and validation phases. As the downloaded daily streamflow (discharge) series was discontinuous, time intervals for calibration and validation were carefully selected to minimize the effect of existing data gaps.

As the CRB is a singular aquifer-fed mountain stream, some quantitative information to crossvalidate the SWAT model results were used. For model efficiency criteria, Nash-Sutcliffe efficiency coefficient (NSE), logarithmic form of the NSE (lnNSE), coefficient of determination (R2), percent bias (PBIAS), Root Mean Square Error (RMSE), and RMSE relative to standard deviation of the observed data (RSR) were used (Table 1).



1 n is the total number of observations, *Qobs i* and *Qsim i* are observed and simulated streamflow at observation *i*, *Q* is the mean of the observed data over the simulation period, and *Qsim i* isthemeanofthesimulateddataoverthesimulationperiod.

#### **3. Results and Discussion**

#### *3.1. Using the CMB Datasets to Estimate IGF*

As shown in Figure 1c, the entire hydrogeological system that contributes to streamflow at the CRB outlet covers the CRB surface itself and some hydraulically connected adjacent areas from GRW and SRW. The methodology described in Section 2.3 was applied to the nodal R values gathered from those 10 km × 10 km cells covering the CRB and those contributing to upstream GRW and SRW areas. Attending to the hydrogeological functioning (Figure 1c) and existing land and water uses (Figure 2), nodal mean values and standard deviations of R and baseflow can be assumed to be equal.

In this area, for the control period (1996–2005), nodal mean R varied within the range 143–332 mm year–1, which means recharge–precipitation ratios were in the 0.29–0.37 range; the standard deviation of mean R varied within the 39–90 mm year–1 range, which placed the given coefficients of variation of mean annual R (mean value-to-standard deviation ratio) in the 0.27–0.30 range (Table 2). For the control period (1996–2005), fitting parameters were calculated to generate the yearly R data series in the CRB and upstream GRW and SRW contributing areas, which are in Table 3, whereas the generated surface-weighted yearly P and R series are in Table 4. In each area, yearly R and P series for the control period (1996–2005) were compared. The resulting parametric functions allowed for the extension of the calculated yearly R series to cover the yearly P full record (1951–2016) (Figure 4). Figure 5 shows the full yearly baseflow series generated within the CBR, as well as the yearly surface-weighted IGF series contributed by upstream GRW and SRW areas. As observed, IGF is somewhat higher than baseflow, generated within the CRB. IGF is about 51% of total CRB baseflow.


**Table 2.** For the 10 km × 10 km cells covering the CRB and upstream GRW and SGW contributing areas, the CMB datasets gathered from Alcalá and Custodio (2014, 2015) [34,42].

<sup>1</sup> Cell ID as in Figure 1c. <sup>2</sup> S is surface in km2, P and R are, respectively, mean precipitation and mean net aquifer recharge over the control period (1996–2005) in mm year–1; and CVP and CVR are the dimensionless coefficients of variation of mean P and R over the control period (1996–2005) as fractions. <sup>3</sup> SWA is surface-weighted average.

**Table 3.** Fitting parameters for the CRB and upstream GRW and SGW areas.


<sup>1</sup> Δm and Δσ are dimensionless, and m*<sup>C</sup>* and σ*<sup>C</sup>* are in mm year–1.


**Table 4.** For the control period (1996–2005), surface-weighted yearly series of (i) P and R in the CRB and in upstream GRW and SRW areas, and (ii) IGF from the GRW and SRW area contributing to CRB.

<sup>1</sup> P and R are, respectively, annual precipitation and net aquifer recharge in mm year−1, and Psi is dimensionless normalized yearly P. <sup>2</sup> IGF is interbasin groundwater flow in mm year−1. <sup>3</sup> Mean and SD are mean and standard deviation over the control period (1996–2005) in mm year−1, and CV is dimensionless coefficient of variation as a fraction.

**Figure 4.** For the control period (1996–2005), parameterization of yearly P–R functions in the CRB and upstream GRW and SRW contributing areas; yearly R equals yearly baseflow. Yearly IGF series refers to the surface-weighed sum of upstream R = baseflow from GRW and SRW areas contributing to the CRB streamflow. In all cases, the Pearson coefficient of correlation is 1.

#### *3.2. Comparison of SWAT Model Results with and without IGF*

Based on DEM analysis and after SWAT model implementation, the CRB was discretized into 29 sub-basins. Based on the combination of land uses, soil types, and slope ranges (<2%, 2%–8%, >8%), 149 HRUs were defined. The thresholds for defining HRUs were set to 5% to optimize model processing. The Hargreaves non-global method was used to simulate potential evapotranspiration [68]. As a result, only precipitation and temperature data to run the SWAT model were needed.

As described in Section 3.1, the IGF from upstream GRW and SRW areas greatly contributes to the Castril River streamflow. For the period 1995−1997, the SWAT model was doubly implemented on a monthly scale with and without IGF. The result was a large difference between observed and initial simulated streamflow when IGF was omitted (Figure 6). When IGF was included as an additional baseflow fraction, the difference between observed and corrected simulated streamflow clearly narrowed. This preliminary trial at model performance showed that the statistics NSE and PBIAS improve when IGF was included (Figure 6). Overall model performance increased about 80% in absolute terms.

**Figure 5.** For the full period (1951–2016), (**a**) surface-weighted yearly P series in the area compiled from the Spanish National Weather Service (AEMET) grid version 1.0 and cumulative deviation (CD) from mean yearly P in mm year<sup>−</sup>1; and (**b**) generated yearly baseflow series in the CRB and yearly surface-weighed IGF series from upstream GRW and SRW contributing areas in mm year<sup>−</sup>1, and IGF fraction relative to total CRB baseflow (IGF–CRB) dimensionless ratio. The control period (1996–2005) is grey shallowed (CP). Vertical dotted lines indicate selected time intervals for the SWAT model warm-up (W), calibration (C), and validation (V) phases.

**Figure 6.** For the selected calibration period (1995−1997) and on a monthly scale, observed streamflow compared to (i) initial simulated streamflow without IGF and (ii) corrected simulated streamflow with IGF. The statistics NSE and PBIAS show the model performance achieved in each simulation.

#### *3.3. Calibration and Validation of SWAT Model Including IGF*

A total of 21 SWAT parameters were optimized using the SUFI-2 algorithm from SWAT-CUP. As described in Section 2.5.2, parameter selection was based on our previous research experiences in similar basins in Southern Spain [64,65]. The SWAT calibration phase covered a 3-year period (1995–1997). The final ranges used and the final fitted values of these parameters are given in Table 5.

The magnitude of calibrated GW\_REVAP, ESCO, LAT\_TTIME, GWQMN, and ALPHA\_BF parameters is quite similar to those obtained with similar orography, geology, climate, and land use [64,65]. The ESCO is also similar to that fitted in other Mediterranean karst areas, where yearly actual evapotranspiration is typically 0.7–0.9 fold yearly precipitation [26,27,46]. The low value of ALPHA\_BF indicates a slow aquifer response [69]. This is corroborated by the long-delayed responses to recharge events in similar karst aquifers in the region, reported by Moral et al. (2008) [70].



(r\_) change, i.e., parameter multiplied by (1 calibration), (v\_) existing parameter replaced by the value obtained in calibration, and (a\_) refers to absolute change, i.e., the fitted value must be added to the existing value of the parameter.

Corrected streamflow records with IGF were used for model calibration (1995–1997) and validation (1982–1984) phases. Observed streamflow was compared to corrected simulated streamflow on monthly (Figure 7) and daily (Figure 8) scales during the calibration and validation periods. In the CRB, the fitted SWAT model replicated, almost identically, the trend of the streamflow hydrograph. The higher fluctuations in the simulated peaks and the lower ones in low flows were found, both in monthly and daily streamflow simulations.

**Figure 7.** On a monthly scale, observed streamflow compared to corrected simulated streamflow with SWAT model for the (**a**) calibration and (**b**) validation phases.

**Figure 8.** On a daily scale, observed streamflow compared to corrected simulated streamflow with SWAT model for the (**a**) calibration and (**b**) validation phases.

As observed in Figure 8, low flows predominate in the CRB daily streamflow record. As many other SWAT models reported for other similar aquifer-fed karst areas, the SWAT model performance for high and normal flows decreased in the face of predominant low flows [65,71]. Therefore, following the

suggestion of Krause et al. (2005) [72], NSE and InNSE were used to measure, respectively, the high and the low flows, to reduce the problem of the squared differences and the resulting sensitivity to extreme values of NSE. Calibration and validation of monthly corrected streamflow showed good agreement of simulated and observed data, as indicated by the model performance statistics for monthly and daily simulations given in Table 6.


**Table 6.** SWAT model performance statistics for corrected simulated monthly and daily streamflow during calibration and validation phases.

As finally deduced, the SWAT model performs well and can be used for further analysis in the CRB and in other similar high-permeability bedrock basins, the baseflow of which is strongly determined by IGF. For this, there must be a confident evaluation of IGF or, at minimum, a reliable external evaluation available. In this paper, the CMB method and the available CMB datasets in continental Spain [34,42] were used for this purpose. However, as described for other groundwater and surface water coupling models, the SWAT-CMB application presented here reveals the following overall problems: (1) the spatial (size and volume) and temporal (renovation rate) scales of groundwater and surface water bodies differ, whereas the coupling model can only simulate the same spatial and temporal scale in both types of bodies; (2) the coupling models had defects in the coupling mechanism processing, which demanded substantial simplification of the coupling process, thereby causing model distortion; and (3) this coupling model was established for a certain region or a specific problem and, although good results have been achieved, there is no general adaptability, so additional hydrogeological knowledge of local applications is needed to consider the changes in scale effects and actual flow conditions.

#### **4. Conclusions**

This paper presents the combined application of the SWAT model and the CMB method to model streamflow in the CRB, a representative high-permeability bedrock basin where the streamflow is significantly determined by IGF from upstream contributing areas. The CMB method and available CMB datasets in continental Spain were used for the IGF that adds to the baseflow generated within the CRB. The SWAT model performance improved noticeably when simulated streamflow with IGF was used. Using the CMB datasets for streamflow correction, the SWAT model showed good performance both in daily and monthly simulations. Some overall remarks from this research are highlighted below.

The influence of IGF on basins like the CRB is remarkable. Therefore, IGF must be considered to improve water resource evaluation and management in this kind of basin located in the headwater of large river watersheds. In the CRB, IGF means about 51% of total baseflow. We do not suggest using the SWAT model alone for modeling of these aquifer-fed mountain basins. It must be coupled with other specific methods to accurately assess IGF. The CMB was revealed to be a suitable method for IGF, because of the favorable hydrogeological setting and the negligible groundwater abstraction, which allowed for equivalent net aquifer recharge to the baseflow contributing to streamflow in the headwater of large river watersheds. In other areas that reflect different patterns of groundwater use and hydrogeological features, assessment of IGF must rely on other techniques coupled with the SWAT model.

**Author Contributions:** J.S.-A. and F.J.A. conceived and designed the research; P.J.-S. and S.L. implemented the coupled application; S.L. compiled data. All authors contributed to writing the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors acknowledge the AEMET, IGN, IGME, CEDEX, and REDIAM Public Institutions for the data provided through their information service websites. We also acknowledge Papercheck Proofreading and Editing Services. We also wish to express our gratitude to two anonymous reviewers for their valuable advice and comments.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **The Ecosystem Resilience Concept Applied to Hydrogeological Systems: A General Approach**

**África de la Hera-Portillo 1,\*, Julio López-Gutiérrez 1, Pedro Zorrilla-Miras 2, Beatriz Mayor <sup>2</sup> and Elena López-Gunn <sup>2</sup>**


Received: 30 May 2020; Accepted: 22 June 2020; Published: 25 June 2020

**Abstract:** We have witnessed the great changes that hydrogeological systems are facing in the last decades: rivers that have dried up; wetlands that have disappeared, leaving their buckets converted into farmland; and aquifers that have been intensively exploited for years, among others. Humans have caused the most part of these results that can be worsened by climate change, with delayed effects on groundwater quantity and quality. The consequences are negatively impacting ecosystems and dependent societies. The concept of resilience has not been extensively used in the hydrogeological research, and it can be a very useful concept that can improve the understanding and management of these systems. The aim of this work is to briefly discuss the role of resilience in the context of freshwater systems affected by either climate or anthropic actions as a way to increase our understanding of how anticipating negative changes (transitions) may contribute to improving the management of the system and preserving the services that it provides. First, the article presents the basic concepts applied to hydrogeological systems from the ecosystem's resilience approach. Second, the factors controlling for hydrogeological systems' responses to different impacts are commented upon. Third, a case study is analyzed and discussed. Finally, the useful implications of the concept are discussed.

**Keywords:** ecosystems; hydrogeological system; sustainability; significant damage; resilience

#### **1. Introduction**

Groundwater and surface water resources are heavily exploited in many parts of the world, and freshwater demands are increasing globally [1]. Any alteration of the baseline conditions of the system may lead to an undesirable state (degradation of quality or quantity). Anthropogenic effects disturb the natural processes of aquifers, and the equilibrium within the unsaturated and saturated zones, and may also increase the contents of undesired substances in groundwater. There are different types of impact affecting hydrogeological systems at different scales. The impacts or disturbances may be either natural or anthropic. Natural disturbances include any type of catastrophe that can affect a hydrogeological system, such as earthquakes, climate (extreme events but also climate variability), or fires. Anthropic disturbances include pumping and various polluting activities such as discharges; agricultural, industrial and nuclear activities; the filtration of substances stored underground; injection into wells; and urban solid waste deposits, among others. Groundwater quality and quantity degradation owing to intensive aquifer exploitation is recorded in many countries [2–8].

This article aims to contribute to develop a better understanding of the concept of resilience when applied to hydrogeological systems, which, in turn, will help develop a better understanding of the buffering capacity of hydrogeological systems. This represents a step to be able to anticipate the potential impacts on the system of specific changes and the system's response. This requires focusing attention on the internal variables which are more sensitive to impacts. The concept of resilience can be helpful to avoid these problems. For example, based on the established groundwater baseline patterns, changes can be identified from the very beginning and can provide an early warning signal to make decisions on sustainable groundwater management. Resilience is also related to water security as groundwater is regarded as one of the most reliable yet also vulnerable sources of drinking water in many countries. This is important since the substantial decline of groundwater levels may affect the water security of a growing economy [9].

There exists an extensive literary record dealing with the resilience of natural systems to different impacts, although most of this literature does not deal with resilience explicitly. Until very recently, the term "resilience" did not appear in hydrogeology glossaries. The idea of the resilience of a "hydrogeological space" or "hydrogeological medium" was developed by LV Demidyuk, NI Lebedeva and GA Golodkovskaya in the 1970s and 1980s, but unfortunately these reports were published in Russian only, e.g., Golodkovskaya and Eliseyev (1989) [10]. There are barely 20 publications in the Web of Science (WOS) returned by the search terms "resilience" and "hydrogeology", and most of these are articles do not treat resilience as a central topic [11–20]. In the field of hydrogeology, the most frequent works dealing with resilience are specific and local, and the concept of resilience is not approached from a generic point of view. Our work aims to contribute to these conceptual reflections.

There is often some confusion in the literature regarding the application of the concept of resilience. Some works apply the concept to groundwater or water resources (liquid phase), and others to aquifers (physical environment). This is an important difference as the intrinsic properties of the system vary in each case. The confusion stems from the fact that the descriptor variable most frequently selected is the same: the piezometric level.

Most literature uses aquifer resilience (AR) as a conclusion derived from the research (Cuthbert et al. (2019) [11] Maurice et al. (2019) [12], Mazi et al. (2014) [14], Chinnasamy et al. (2018) [16], De Eyto et al. (2016) [17], Hejazian et al. (2017) [21]), with the largest number of works focused on the analysis of drought as a disturbing element (Lorenzo-Lacruz et al. (2017) [13], McDonald et al. (2017) [22]. Two references present a deeper study of the resilience of aquifers. First, Bouska et al. (2019) [23] apply the concept of general resilience to the restoration of large river ecosystems in the Upper Mississippi. However, their approach is more biological [24,25] than hydrogeological, i.e., they develop indicators for three principles of general resilience: diversity and redundancy, connectivity, and controlling variables. The latter includes historical water level fluctuations, water clarity, nutrient concentration, and invasive aquatic species. Wurl et al. (2018) [26] adopt an approach that is closer to the resilience concept from a hydrogeological perspective. The authors designed and used a set of indicators as outcomes for combined human–water systems to predict water trajectories under different human impacts.

This literature review identifies a number of tools that have been used indirectly to build the application of resilience into hydrogeology: the use of tracers or isotopes to determine groundwater age [27], the analysis of piezometric evolution trends in aquifers [28], and the quantification of water consumption for agriculture in restricted aquifers [29]. Quantification also relies on qualitative indicators as an auxiliary tool [16]. However, the term "resilience" in the hydrogeological literature is presented as an attribute that is not analyzed in itself, but which is frequently cited within the framework of processes (recharge, precipitation) associated mostly with climatic variability, and, to a much lesser degree, with other processes (extreme events).

The aim of this work is to increase our understanding of how hydrogeological systems deal with disturbances as a way to anticipate transitions (changes in the system affecting their functioning) and to propose a conceptual model for the analysis of the resilience of aquifers. In turn, this knowledge can support the more sustainable management of the system. The paper compiles the key knowledge around this concept that can be applied to hydrogeological systems and discusses a relevant case for illustrative purposes.

The article is structured as follows: first, the article presents the basic concepts applied to hydrogeological systems, and the conceptual model proposed, including factors controlling the responses of aquifers to different types of impacts. Second, a series of considerations are made regarding its scope. It concludes with a discussion and considerations of how to deal with the resilience property in the framework of long-term data series to obtain analytical and useful results.

#### **2. Hydrogeological Systems and the Resilience Concept**

#### *2.1. Hydrogeological Systems as Complex Systems*

Hydrogeological water systems operate with a certain behavior for a certain period at the human time scale [30–35]. For example, depending on their water regime, rivers can be classified as permanent, seasonal, temporary or intermittent; wetlands can be permanent, temporary, seasonal or erratic based on their hydroperiod; and aquifers can be classified as unconfined, confined, semi-confined or perched according to their operation.

In many areas of the planet, hydrogeological systems (water systems) are subject to different types of stressors, such as the extraction of water for irrigated agriculture, or in urban coastal areas contamination from localized or extended sources and changes in recharge regimes due to climate change, which may cause the systems to exceed the limits of sustainability [36]. If this occurs, the systems become unbalanced, leading to a tipping point [37] in their behavior, bringing them to a new state of equilibrium [38].

#### *2.2. The Resilience Concept and Theory*

Resilience is fundamentally directed to the way a system responds to a disturbance [39,40] that can be punctual or a long-term process, for example one implying gradual alterations (slow onset changes) [41]. This concept originated in metallurgy and has subsequently been applied in many other disciplines. In the field of ecology, the concept is based on the theory that systems are in a natural state of flux rather than an equilibrium [42]. In this work, we define resilience as a system's ability to recover a situation of equilibrium or metastability (known as state, see the definition below), characterized by a known behavior. We do not see it as a return to its pristine conditions for two reasons: (a) in many areas, there is no information available to define the pristine conditions, i.e., it refers to a period prior to registration for which there is no data; (b) in general terms, the systems are altered (by humans or by other natural processes) in one way or another, so that their return to an pristine state may be a goal that is significantly outside the realms of possibility.

The literature on the resilience of complex systems is highly fragmented [43,44]. One of the main problems that emerges is the lack of terminological consensus among the authors [44,45]. Different authors use a range of terms to refer generically to the same concept of "complex systems", including "complex adaptive ecosystems" [46,47], "complex adaptive systems", "complex, coupled Socio-Ecological Systems (SESs)", and "complex, multi-scalar SESs" [48]. In some cases, they include references to impacts, events, shocks, pulses, threats or stress, leading to a terminological confusion that must be avoided in a scoping study such as the one presented here.

From a biological perspective, the concept of resilience is especially applicable to natural systems that adapt to different degrees of disturbance while maintaining the same processes and structures that reinforce each other [49,50], and whose connection is known as a "regime" [51] (Table 1). The transition from one regime to another (regime shift) occurs through thresholds (Figure 1); and the new regime is characterized by a different set of processes and structures (behavior) [52]. Regime changes are typically associated with significant consequences in processes or structures (e.g., a change in water composition that leads to a loss of water quality), and do not always occur in sudden leaps or at turning points in their trajectory, but may be the result of long system periods [53] or slow and progressive changes [41]. Not all regime changes entail a tipping point. Indeed, Bertalanffy (1968) [54] identifies six mechanisms that can trigger a regime change (slow–fast cyclic transition, stochastic resonance, noise-induced transition, long transient upon extreme events, big stepwise changes in drivers), but only one of them involves a tipping point: slowly changing driver to tipping point. Four out of these six transition indicators can be used as early-warning signals [37,55,56]. In the field of hydrogeological ecosystem research, the aspects of the dynamics of changes between the states of equilibrium are relatively unexplored.

**Figure 1.** Scheme of the aquifer resilience (AR) through the temporal evolution of a descriptive variable. The black lines indicate the time series observed of the aquifer system, and the red lines represent the time series of the underlying environmental conditions. (**a**) In a theoretical example, the pumping of groundwater from an aquifer (disturbance) prolonged over time will cause a water table depletion (parameter) which could reflect the aquifer response to the impact. (**b**). If the pumping ceases (disturbance stops), the aquifer may recover a new equilibrium (State-3). Both behaviors (**a**,**b**) are useful to understand aquifer resilience and must be taken into account in the interpretation of the AR (Source: authors' own).


**Table 1.** Extended definition of terms associated with resilience in a biological environment. Text in italics indicates the literal definition of the authors. Text in brackets are comments from the authors. When the source is not indicated, the reference is proper.


**Table 1.** *Cont*.

The scope of the concept of ecosystems resilience is broader than initially considered. When discussing ecosystems alone, resilience is closely related to sustainability ([60,61]). Scheffer et al. (2001) [62] report that a loss of ecosystem resilience generally paves the way for a change to an alternative state and suggest that sustainable management should be directed towards maintaining ecosystem resilience [62].

Resilience also refers to the system's adaptive capacity [57] and vulnerability, given that it offers another approach to the changes produced by a disturbance; however, this idea is controversial, as some authors consider that the concept of adaptive capacity is muddled with multiple meanings in current use often being indistinguishable from resilience [41]. Resilience is an intrinsic property of the system that emerges from certain changes. Some authors define the adaptive capacity of ecosystems as a latent potential quality to alter resilience in response to change [57].

Resilience has implications in the socio-economic and political spheres, since knowing the dimensions of this attribute enables managers to intervene in the natural environment before a change of regime or favoring one state of equilibrium over another.

Although the term "resilience" is increasingly used by political and environmental managers, it remains vague, variable and difficult to quantify [39]. This work clarifies what it means from the hydrogeological perspective, without attempting to review the state of the art of the concept or to list an inventory of works in which the concept is applied to geological and hydrogeological studies.

#### *2.3. Resilience from Ecology to Hydrogeology: A Conceptual Framework for Its Analysis*

The analysis of resilience focuses on the dynamics of the system, particularly looking at two areas: the cause effect (disturbance), and its consequence in the system (system response to the disturbance). The disturbance leads the system to alter certain internal variables, which define the new equilibrium state (regime). Based on the concept of resilience and on the literature dealing with the resilience of hydrogeological systems, we propose a conceptual model that can use aquifer resilience to support its management (Figure 2).

**Figure 2.** Proposed conceptual model for using the resilience concept in the management of aquifers (Source: authors' own).

Step 1. The first step must be the system definition, including its geological boundaries and its main characteristics, flows and functions. An important issue to be addressed in any resilience analysis is the scale of work (space–time dimensions).

To understand hydrogeological system resilience, cause–effect relationships, and impacts, it is necessary to have information about its functioning. Groundwater flow systems depend on both the hydrogeological characteristics of the soil/rock material and on the landscape position [3]. These factors control the permeability conditions and the hydraulic gradient differences which regulate the groundwater flow movement. Not every impact affecting the Groundwater Flow System (GWFS) affects its hierarchical structure and functioning; this will depend on the nature, magnitude and duration of the impact, and the factors controlling the response of GWFS to those pressures.

Step 2. The second step is to describe the process triggered by the disturbance (some authors refer to them as descriptive variables of the change and the interactions between them) that is, to describe its magnitude, duration and scale. This involves monitoring the variables, describing the change before, during and after the change. The second step is therefore to define and describe the disturbance that acts on the system (state 1) causing a series of internal changes that lead to a situation of instability for a certain or indefinite time (Figure 1). The usual process is for systems to tend to a new state of equilibrium (state 2) through internal changes and interactions between the variables that describe the change.

In describing the disturbance, there are a set of elements that need to be considered when analyzing resilience. One of them is the time scale of the system and of the disturbing forces. In the natural environment, some internal changes occur over short time periods and are visible on a human time scale (for example, change in the eutrophication conditions of a lake, reduction in the population of a certain insect, etc.). In other cases, the effect of a disturbing agent may not become evident for years, millennia or millions of years, and thus be difficult to determine, for example given the different temporal scales of geological processes.

Another element to consider is the possible overlap of effects due to the vast dimensions of a system and the difference in the periodicity and breadth of the various antagonistic processes. For example, in a large detrital aquifer with an immense storage capacity, a short extremely dry period can be obliterated by the hyper annual natural recharge of average and humid years, i.e., the overlapping of the previous and subsequent average recharge would have cushioned these effects.

Step 3. The third step is to describe the new state of equilibrium (regime) [45].

The analysis of resilience focuses on the disturbance, the processes of change, and the system's recovery into a new state of equilibrium. If the concept of resilience is to be made operational, we need to find ways to measure it. Therefore, an element to consider is the fact that the transition from one state (or "alternative state" according to some authors) to another occurs through thresholds [63]. Some authors argue that these thresholds are crucial for measuring resilience as these offer a way to quantify how much disturbance can be absorbed by a system before switching to another regime. The identification of these thresholds requires experimental or observational data on the changes between regimes in a certain system and, if possible, on its recovery trajectory [39]. Although there is considerable work done on transition indicators in biological systems [56], this is not the case for geological systems. A long time series of data is not always possible. In these cases, there are other resources to carry out this analytical study, as we show in the next section.

#### **3. Conceptual Model Applied to a Real Case: The Upper Guadiana Basin in Central-West Spain**

The case study presented here is for illustrative purposes to show through an example why the lens of resilience is valuable to gain a better, more anticipative knowledge of the system. It is also a pertinent case because it is applied to the functioning of a very large aquifer, its relation with an important groundwater-dependent wetland (the Tablas de Daimiel National Park, TDNP), and the dynamic of the system from a starting point, to highly degraded systems (both aquifer and wetland), and to a current significant recovery.

In the following sections, a description of the main steps mentioned above is presented. Figure 3 summarizes the main factors considered.

#### *3.1. First Step: Description of the System, Flows and Functions*

This case provides, on the one hand, the aquifer perspective through the Mancha Occidental aquifer, and, on the other hand, the wetland perspective represented by a groundwater-dependent ecosystem affected by multiple impacts with enough data to allow its analysis. The equilibrium in wetland ecosystems is very fragile, showing high sensibility and vulnerability [64–67].

Tablas de Daimiel is a groundwater-dependent ecosystem subject to different types of impacts, both climatic and anthropogenic. This wetland is the main discharge outlet of the Upper Guadiana basin's aquifers (Figure 4), in such a way that it can be considered the "thermometer" of the 16,000 km2 groundwater system [4]. The Tablas de Daimiel wetland has existed for over 250,000 years, evolving

from a deep lake to a fluctuating shallow system, with different reversible intermediate phases depending on hydroclimatic conditions [8]. In 1960, the system water inflows combined brackish surface water from the Cigüela River with freshwater inputs from the Guadiana River and the underlying aquifer.

**Figure 4.** (**A**) General geographical setting. (**B**) Upper Guadiana basin. (**C**) The Mancha Occidental Aquifer (currently divided into three groundwater bodies: Mancha Occidental I, Mancha Occidental II and Rus-Valdelobos) (Source: Authors' own).

The analysis of the flood data series for the period 1944–1974 (Figure 5), prior to the overexploitation of the aquifer, reveals that in that period changes in rainfall (Figure 6) determined changes in water variability: changes in rainfall determined changes in the wetland surface covered by water [68].

**Figure 5.** Flooded area evolution of Tablas de Daimiel National Park. Data from Sánchez-Carrillo et al. 2016 [68]. At the end of May 2020, the flooded area is less than 80 ha, and a new water transfer has been requested by the Tablas de Daimiel National Park (TDNP) Director to the Government (Source: Authors' own).

In resilience terms, 1944–1970 corresponds to the pre-pumping stage of the system and can be used as a reference stage for our purpose.

**Figure 6.** Long-term rainfall patterns in Mancha Occidental Aquifer 1904–2014. Rainfall data from weather stations 4121 and 4121C (data from [8]) (Authors' own).

#### *3.2. Water Quality in Tablas de Daimiel National Park: Baseline Stage (September–October 1974)*

Electrical conductivity ranged between 800 and 1600 μS/cm, with a minimum of 300 and a maximum of 5400 μS/cm. The dominant anions were sulphate and bicarbonate, whereas the dominant cations were calcium and magnesium [69]. Calcium bicarbonate waters predominated to the northwest and northeast of the national park, as well as in the vicinity of the Ojos del Guadiana springs, with the lowest conductivity sampled. Meanwhile, calcium sulphate waters predominated around the left bank of the wetland [4].

#### *3.3. Second Step: Description of the Perturbation*

Since the mid-1970s, the intensive exploitation of the aquifer for agricultural irrigation caused the desiccation of the wetland and neighboring springs. The cause of the hydrological situation of the Tablas de Daimiel in the 1980s can be well explained due to the length of the time series covering period 1975–2008. Extensive descriptive publications exist about this period examining its origin, reasons and social-ecological consequences [5–7].

Based on data provided by Aguilera and Moreno (2018) [69], the impact of drought caused a decrease in groundwater levels, which, at the same time, produced the burning of peat underlying the wetland, causing a smoldering peat fire in 2009 [69]. In this case, the descriptive variables observed were soil moisture, temperature and organic matter content for the period 2006–2010 [70]. Continuous soil moisture and temperature monitoring is recommended as an indicator of potential combustion and auto-ignition fire risk but does not work as alert system for an already active fire. In fact, the presence of active smoke columns is a late warning. A new fire means to arrive late. Fire modifies irreversibly the physical structure of affected soils, which implies a damage to the ecosystem.

In resilience terms, this analysis shows how the system faces the burning of peat impact. This provides the added value of the reaction to the change in the hydrological conditions of the soil.

#### *3.4. Third Step: Description of the New State of Equilibrium*

The longest water table records show the evolution of the system through the different described impacts (Figure 7). With depleted piezometric levels, the TDNP wetland operates as a recharge system for a local shallow perched multi-layer aquifer disconnected from the deeper regional groundwater flow towards main irrigation areas [71]. Water-table records show the tendency of the system to behave in a roughly similar manner across its entire extension.

**Figure 7.** Water table evolution in a representative monitoring well (1930.4.0040) located near the Ojos del Guadiana springs and Tablas de Daimiel National Park. It ranks among the longest available piezometric records. The figure also shows the main historical events occurred in the (**a**) Tablas de Daimiel National Park: 01/1986–87: peat fires; 02/1987: Tablas de Daimiel was dry up for the first time; 03/1988: the first water transfer from Tagus-Segura Aqueduct; 04/1994: peat fires; 50/1995: flooded area, 30 out of 2000 ha; 06/1996–98: heavy rains; 07/2009: peat fires; 08/2010: water transfer from Tagus-Segura Aqueduct; 09/2012: groundwater feed Tablas de Daimiel. (**b**) Mancha Occidental Aquifer: 01/1983: Ojos del Guadiana springs dried; 02/1987: provisional declaration of aquifer overexploitation; 03/1988: maximum pumping peak (570 million m3); 04/1994: declaration of aquifer overexploitation (pumping regulation); 05/1995: water level depth 47 m in Ojos del Guadiana springs; 06/1996–98: heavy rains; 07/2000: Water Framework Directive (2000/60/EC); 08/2009: water level depth 35 m in Ojos del Guadiana springs; 09/2010: heavy rains; 10/2014: Declaration of the Upper Guadiana basin's groundwater bodies (GWBs) at risk of not achieving the good quantitative and chemical status. (Data from [7]) (Authors' own).

For three decades, the wetland remained in precarious hydrological conditions, with the only exception of rapid floods due to extreme rainfall events and sporadic water transfers from the Tagus river basin. The water transfers from the Tajo-Segura Aqueduct started in 1988 and have often been carried out during spring or summer when evaporation and infiltration rates are highest and when increased demands of water for irrigation promote illegal extractions. Flooding the wetland with

pumped groundwater was a management tool used constantly during dry periods to keep a minimum flooded area. Both management tools induce quantitative and qualitative impacts on the system [69].

There are several stressors acting simultaneously on the Tablas de Daimiel National Park: droughts (drying) and flooding situations, pumping, fires, treated wastewater, water transfer from other basins, and land use changes, among others. It is not possible to isolate just one cause and one effect for the analysis.

#### *3.5. Water Quality in Tablas de Daimiel National Park: Pumping Stage*

In 2007, the number of analyses is smaller than in the other two stages due to the water table dropping below the bottom of some piezometers and all springs drying up [7]. Electrical conductivity ranged between 2000 and 11,000 μS/cm. In the vicinity of the wetland, the main cations evolved into sodium and magnesium, while the dominant anions leaned towards sulphate and chloride. The surface water exerts little influence on groundwater chemistry across most of the system.

In the case study, an unusual situation has recently occurred: in 2011, a decrease in groundwater abstraction and an extraordinary wet period reversed the trend. Following a wet period (2006–2009) capped by an exceptionally humid year (2010), the aquifer experienced an unexpected recovery of groundwater levels (almost 20 m in some areas), restoring groundwater discharge to springs and wetlands, which came back to life for the first time since the early 1980s (Figure 7). For the sake of brevity, this history may be found in [8,68]; for a more recent perspective and hydrogeological functioning, details are presented in Castaño et al. (2018) [7] and Martínez-Santos et al. (2018) [8]. Here, we just mention those aspects relevant to the aim of the paper.

#### *3.6. Water Quality in Tablas de Daimiel National Park: Restoring Stage*

Data from 2014 show that the hydrological recovery has not yet been mirrored by a similar recovery in water quality. In fact, the groundwater is more saline than it used to be before the 1970s, and the predominant hydrochemical facies have shifted with meaningful spatial gradients [7]. The descriptive variable is the water quality of the groundwater around the wetland considering three stages: (i) prior to degradation of the wetland (1974, baseline); (ii) during a period of major degradation (2007); and (iii) after the most important recovery on record (2014). The pictures of Stiff diagrams obtained for each one of these milestones respond to a different state of the wetland, in this case, without data of transition among them, which is key to identify the beginning of the changings and provides the most valuable information to make responsible management decisions. In spite of this, the second picture (2007) shows an important change in water quality which indicates a remarkable internal change in the ecosystem after 33 years of pumping. The 2014 hydrochemical data indicate that the hydrological recovery of the system refers exclusively to the water balance and not to the water quality. This means that the lag time for water quality to reach a new equilibrium is shorter than the lag time of water levels to reach the position close to that of the 1970s. In fact, the reaction of groundwater levels to rainfall and decrease of pumping is fast (a short period of a few months), whereas the evolution of groundwater quality is much slower. Moreover, there are no studies to predict whether or not it will be reached or when. Considering the high level of hydrogeological knowledge existing in this area, including hydrogeological flow models, some research could be done along these lines in order not only to predict whether a new equilibrium will be reached, but also whether the system is able to keep such a state for a long time. Perhaps new actions will be needed in the context of sustainable management decisions.

#### *3.7. Surface Flaming Fires and Smoldering Peat Fires*

Both surface flaming fires and smoldering peat fires have been relatively frequent in the TDNP surroundings (1977, 1987 and 1991) [7,72]. In fact, most natural peatlands outside the park limits have disappeared. Smoldering peat fires have even been reported inside the park in 1986, 1987 and 1994 [67] (see Figure 6). But they occurred under relatively wet soil conditions, with a shallow water table located less than 1 m below the surface, and only affected small areas. In the 2009 fires, on the contrary, the soil moisture was much lower, and the water table was located deep below the surface, so the fires represented a much bigger problem [70]. This means that this fire posed an enormous risk for both the physical structure supporting the ecosystem and the quality of groundwater beneath it. The analysis of key parameters monitored in several locations of the TDNP at different depths shows that there was enough previous evidence to foresee the peat self-combustion and the risk that any surface fire could be transmitted to the subsoil. Data were taken in the vadose zone of the TDNP up to a depth of 2 meters at 12 points.

In resilience terms, this analysis shows that the soil's organic carbon content and moisture are two key variables in smoldering fires. The first is related to the amount of fuel available in the soil. The second represents a threshold condition, as below a certain moisture content peat can burn. This means that the peat combustion could be predictable allowing for pro-active management.

In resilience terms, this analysis shows that the system is more resilient to water quantity than to groundwater quality. This provides an added value for the assessment of conservation strategies.

#### **4. Discussion—The Need for Good Quality Long-Term Data**

From the analysis undertaken above, on understanding a complex hydrogeological system like the Upper Guadiana aquifer, through the lens of resilience we can gather new insights on the functioning of the system. For example, this wetland has been exposed to impacts that are not always evident and reversible [73]. The first signals of a change in water quality could have been detected if an adequate system of monitoring and data interpretation had been performed since 1970. An earlier intervention in the system could have avoided the degradation of water quality suffered currently in the wetland. Most of the processes triggered by global changes were not detectable in the short term; instead, it is necessary to adopt a longer decadal scale to understand their dynamics and evaluate their consequences [68].

The most frequent controlled variables in flow river systems are discharge data. In large rivers, these records are normally well registered. Nevertheless, in many areas, gauging stations do not have continuous records, making it difficult to undertake long-term series analysis. It is important to note that over recent decades, hydrological regimes have been changing at a very fast pace. Some progress has been made in extracting long-term signals of change from hydroclimatic data. However, further studies investigating long-term changes in river runoff, and focusing on the detection of underlying mechanisms and the disentanglement of their effects are needed [73].

The most frequent variable observed in hydrogeological studies are groundwater level fluctuations and periodical groundwater samples analysis. These are the variables controlled in most groundwater monitoring networks. It is important, therefore, to have good quality data records on groundwater abstractions to investigate the links between groundwater abstractions and their potentiometric surfaces to better understand future aquifer responses to climatic and anthropogenic stresses [74].

It is not easy to find continuous flow records from springs, since only a reduced number of springs have this information available. However, in order to identify flow patterns, other variables are needed, such as electrical conductivity and temperature, that would need to be monitored simultaneously and synchronically to the flow record. In wetlands, as in any other surface water body, monitoring the height of water in the wetland flooding control [68], the groundwater level in some close wells, and water samples from the wetland's water and groundwater [75] should also be recorded continuously [67].

While many groundwater and surface water flow systems have long-term operation histories, they do not have long-term series of data to assess such operation in depth. Very frequently, this is due to a lack of budget or a lack of staff and time to interpret these records. The hidden information behind those long-term series data should be extracted through long-term trends analyses from which some processes can be identified. Galassi et al. (2014) [15] studied the results of the effects of a 6.3 Mw earthquake on 6 April 2009 on the Gran Sasso karst aquifer in L'Aquila (Italy) by comparing biotic and abiotic data from two years prior to the event (1997 and 2005) and another post-event (year 2012), although not with contiguous hydrological years. This highlights the lack of data available to conduct this type of analysis.

A sufficient length of the time series is vital to be able to distinguish between different impacts, for example, natural climate variability and signals of climate change. When adopting strict criteria regarding data length and data quality, the available information probably decreases.

Different time horizons of observations and measurements can lead to different conclusions, and therefore time extrapolations are always risky. Hierarchy theory states that there is a control mechanism in the temporal order of systems and suggests that long-term processes (that operate mainly on wide spatial scales) restrict fast processes (that work on small spatial scales), which limits their degrees of freedom [76,77]. Although these concepts have theoretical strength, their empirical evidence has not been widely demonstrated so far due to the lack of data sets [68].

It is thus argued that better management decisions could be made if they were based on managing the resilience of systems rather than maintaining them as if these systems were inherently static and thus aim to return them to the statistics of, e.g., 50 years ago. In reality, natural systems are dynamic, and even more so when combined with anthropic influences, with the impacts of multiple factors of global change not present or, at least, not having the same intensity [78].

It is essential to identify the descriptive variables of the changes in order to monitor the analysis. Without adequate monitoring of these variables, it is impossible to understand the change dynamics and their scope, duration and characterization. This monitoring should be ex ante, taking a good design of the spatial observation network into account as well as an adequate periodicity of reading or sampling. Moreover, this monitoring should be permanent in order to allow data from the pre-disturbance phase to be available during and after the disturbance. Precisely, one of the reasons for the scarcity of these studies in hydrogeology is undoubtedly the absence of monitored information on geological processes. Long-term data records are required, combined with an observation network with good spatial coverage [79,80], to facilitate the analysis of the system's resilience. Furthermore, data limitations and the lack of information on mechanisms and processes pose significant limitations to research in many systems [19,75]. Some anthropogenic interventions may imply permanent or long-term durable changes, like the construction of buildings or the start of a new groundwater competing sector, such as agriculture. Permanent or long-term durable changes (in a human scale) imply that it is often difficult to return a system to its initial conditions in the temporal and socio-economic spheres. It makes more sense to talk about resilience today in terms of considering that the system attains a new state of equilibrium under the new conditions of change caused by a disturbance [79,80]. After the cessation of a disturbance, it can be possible to return to the initial state if no new disturbances occur, although this depends on the recovery capacity of the system. Nevertheless, it is important to point out that initial does not mean pristine. In cases where an initial state is known, the management could try to return the system to this initial state as opposed to the pristine state.

To date, transition indicators can be defined in certain natural systems, including volcanic (terrain changes prior to eruption) and hydrological systems (surface water and groundwater quality changes). However, many system state changes, let alone state change thresholds, can only be roughly recognized [79,80]. Also, early-warning signals are an open field of work and will be the next step once the system transitions can be identified.

The study of the resilience of natural systems requires multidisciplinary research including teams of experts in different fields. For the case of wetlands, it is necessary to know not only their characteristic functions but also the interrelation between hydrogeological and biological processes, and particularly the dynamic of governing and socio-ecological variables. The compartmentalization of the natural environment into separate disciplinary fields merely reduces the visual scope and skews the dynamics of the natural processes that develop between both spheres conceptually. Collaboration between the different disciplines enriches this vision, guarantees a more effective approach to reality, and is the safeguard of true scientific advancement.

#### **5. Conclusions**

The analysis of the resilience of any natural system must be defined based on its twofold nature, that is, from what to what, and how. The work must focus on analyzing the type of disturbance and identifying the changes produced by this disturbance in the internal dynamics of the system (the operation of the system), both within and between states.

It is important to consider the different scales of analysis and precisely define the trigger for the changes (impact) as the external agent acting from a higher scale, and the space–time scope of what is considered the "system" under study. Within this system, the variables that describe the change must be identified, along with their evolution, the interaction between them, and their potential recovery once the disturbance ceases or the shock is cushioned.

Our exemplification through the case of the case of Tablas de Daimiel describes how a set of changes have been caused by a series of impacts. These impacts act simultaneously making difficult any correlation of cause–effect binomial. Since the system is changing in a complex way, as a consequence of global changes, the only option is to obtain good quality data with a long enough time series to be able to discriminate a system's responses to different causes. The changes occurred in the Tablas de Daimiel are a response to multiple disturbances, and their interpretation varies with the time scale considered [68]. The response of many processes triggered by different changes is reflected in a time lag that is impossible to detect with short observation periods. The Tablas de Daimiel is currently a system in a new state. In this paper, some examples have been shown covering how and at what speed ecosystems have moved their structure to this new state, using flooding, rainfall and groundwater records jointly with additional short periods of specific data (water quality, and soil moisture, temperature and organic matter content). The Tablas de Daimiel shows a high resilience to droughts and flood events and a low resilience to pumping and fires. This means that the system copes better with natural disturbances than anthropic disturbances. However, no measurements or estimations of the resilience change starting point could be made due to the lack of continuous flow recording data. Thus, good data series are key to having a strong conceptual understanding of the resilience of hydrogeological systems that in turn allow for a more adaptive style of management that better reflects that systems are not static but rather are constantly evolving. It is thus critical to understand this so that system resilience is in line with the protection of key hydrogeological system functions.

**Author Contributions:** Conceived, designed and original draft preparation: Á.d.l.H.-P. Formal analysis and review: J.L.-G. and P.Z.-M. English style and editing: E.L.-G., B.M., P.Z.-M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by European Commission H2020 program: NAIAD (Nature Insurance value: Assessment and Demonstration, grant number 730497).

**Acknowledgments:** This work has been funded by the H2020 Program of the European Commission Project "Nature Insurance value: Assessment and Demonstration" (NAIAD) No. 730497. The views expressed are those of the authors alone. The authors wish to thank to the anonymous reviewers for their useful and constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

#### *Article*
