Model-based assessment of interbasin groundwater flow in data scarce areas: the Gallocanta Lake endorheic watershed (Spain)

Abstract

Aquifer systems, because of the presence of frequently complex geological structures, may extend beyond watersheds limits. Interbasin groundwater flow is often identified among watersheds. Because geological systems are complex ones, modelling tools are needed for its estimation. In this paper, we quantify the outflows from the endhoreic Gallocanta watershed (Spain) by means of a MODFLOW numerical model in order to assess the robustness of the boundaries of the Gallocanta Groundwater Body. Our results show the Gallocanta watershed is hydrogeologically connected with the adjoining Piedra-Ortiz and Jiloca watersheds (discharging annually in these basins about 4 and 1 Mm³ respectively). Furthermore, we hypothesized the presence of geological features altering the groundwater flow. Additional simulations were run to analyse the changes in the water budget in the cases of: i) groundwater pumping no longer allowed by the authorities, and ii) a potential drought scenario. In the first case, the results forecast an increase in discharge to the Piedra-Ortiz and Jiloca watersheds, while in the second case a diminution of the outflows to the two neighboring basins is foreseen.

We then propose a larger and unique groundwater body, spanning from the Caminreal Springs on the east and the Piedra-Ortiz basin on the west, including a moving groundwater divide internal to the Gallocanta watershed. Monitoring the baseflow of the Piedra-Ortiz river and of the Caminreal Springs will allow to get information on the evolution of the groundwater resource availability in the Gallocanta watershed. Our results stress the importance of conjunctively using data and traditional geologic knowledge (i.e. surface geology maps) along with numerical modelling analyses. This holds especially true in areas, such those of hard-rock aquifers, where scarce hydrogeologic data are available, to test conceptual models, to derive and to infer information on water budgets and on the presence of relevant structural features driving the groundwater flow. This approach may lead to informed decision-making on groundwater body boundaries definition for the application of relevant groundwater management regulations.

1 Introduction

While topographic catchments are generally considered as self-contained, conservative, hydrological units, without exchange with adjacent ones [1, 2], aquifer systems, because of the presence of frequently complex geological structures may extend beyond watersheds limits. To this regard, interbasin groundwater flow (IGF) is often identified among watersheds. It occurs when groundwater flow crosses surface water drainage divides and discharges in an adjacent topographic basin [3].

The identification and assessment of IGF is of relevance in hard-rock aquifers in arid and semi-arid areas, especially for what regards endorheic basins, as it allows defining the whole hydrological cycle [4, 5] and addressing proper water resources management issues [6]. In some regions, IGF constitutes a significant part of the water budget, and being not measurable, it is estimated by means of modelling tools [7].

In hard-rock aquifers groundwater flows is difficult to analyse not only because of the heterogeneity of hydrostratigraphic units, but also due to often complex structural geological settings [8]. The presence of faults and thrusts may dislocate geological units and even facilitate or stop groundwater flow. Field hydrogeological studies complemented by physically-based, spatially-distributed groundwater flow numerical modelling allow to set a digital twin of the investigated domain and to gain knowledge of complex groundwater hydrodynamics including computation of IGF. The simulation based on a mathematical model is nowadays a tool commonly used to represent several processes of the hydrological cycle [9, 10].

However, few modelling studies aimed at quantitatively assessing IGF [1, 11] in hard-rock aquifers where data are commonly scarce. Model-based assessment is a valuable method, even in scarce data conditions, as it allows to quantify inflows and outflows through a system, to identify and to direct further investigations. Within this framework, in this paper, we assessed the groundwater budget of the Gallocanta Lake watershed (Spain) by means of numerical modelling considering the rate of IGF exchanges with adjoining basins, the influence of geological structures on the groundwater flow, and the relationships between surface- and ground-water bodies in order to support data-based groundwater resource management.

2 Material and methods

The Gallocanta watershed, approximately 35 km in length and 20 km wide, is an endhoreic basin located in the Iberian Range (NE Spain), an alpine chain surrounded by the Ebro (NE), the Duero, and the Tajo Basins (W) (Fig. 1). The Gallocanta Lake, giving name to the watershed and located at about 1000 m above mean sea level, constitutes the main hydrological feature of the area.

Fig. 1

Geographical setting of the investigated domain

2.1 Climate characteristics

The climate of the area is sub-arid. The annual rainfall average for the period 1950–2020 was 440 mm yr⁻¹ (inferred following the method proposed by Serrano-Notivoli [12]; Fig. 2). For the same period of time, the yearly minimum precipitation was 173 mm (in 2019) and the maximum was 638 mm (in 1959). The major rain events are not evenly distributed in the year, with most rain falling between February and June, and the summer period characterised by low values of cumulative rainfall. The average temperature varies between 3 °C during the winter period and 21 °C during summer (January–July), with average annual temperatures of 12 °C. This results in cold winters and hot dry summers. Luzon et al. [13] calculated a mean annual potential evaporation of 926 mm.

Fig. 2

Monthly average rainfall and temperature in the Gallocanta basin (period 1950–2020; data from the Spanish Meterological Agency)

2.2 Geological setting

The Gallocanta watershed (Fig. 3) is seated in a vast area dominated by Mesozoic outcrops slightly folded along the SE sector of the Iberian Range, delimiting the largest fluvial basins in the Iberian Peninsula (i.e. Ebro Basin, Duero Basin and Tagus Basin). The Iberian Range is a fold-and-thrust range, whose structures are essentially NW–SE oriented, related to Mesozoic and Variscan orogeny [14]. Several authors studied the geology of the area from a lithostratigraphic and structural point of view [15, 16]. Paleozoic Sierras, made of metamorphic quasi-impervious rocks (quartzites and slates, alternated with meta-siltstone and meta-sandstone), whose depth is not known, limit the basin on the north-eastern and south-western boundaries. The central part of the basin (Figs. 3 and 4) is mostly flat, covered by Quaternary and Tertiary sediments (sand and clay up to 20 and 50 m thick in the center of the basin, respectively), above Triassic deposits of the Keuper Facies. The Keuper Facies includes clays and evaporites, intensely folded, whose thickness is variable due to its plastic behaviour. This level acts as a detachment layer for the overlying structures. Below the Quaternary sediments and the Keuper ones, between the Paleozoic rocks and the lake, a thin belt of Triassic rocks (Muschelkalk and Buntsandstein Facies) outcrops in very limited areas. Its geometry and depth are not well-known. The west and south parts of the basin are covered by Mesozoic dolomitic, clays and limestones (Jurassic and Cretaceous in age) reaching up to ~ 300 m maximum thickness.

Fig. 3

Geological setting of the Gallocanta Lake watershed

Fig. 4

Geological and hydrogeological cross-section (cross-section line in Fig. 3)

2.3 Main hydrological features

The Gallocanta Lake is a shallow (2 m maximum depth), ephemeral salt lake, and it is the largest of a regional lake system constituted by several small lakes, such as the Zaida Lake (north of the Gallocanta Lake). The Gallocanta Lake topographic basin is endorheic, so all the streams flow to the lake. The lake is the largest saline lake in Spain and all western Europe and is seated in an eminently rural area. Despite seasonal variations (the lake frequently dries up in the summer season) the lake extends over an area of about 14–15 km². Luzon et al. [13] reports salinity values ranging from 3600 to 136,000 µS cm⁻¹. Due to its relevant natural characteristics, the lake and its surroundings have been studied from different disciplines, such as geology, sedimentology [16], edaphology, or hydrology [13].

According to Ebro Hydrographic Confederation (CHE) [17], the water budget of the lake is constrained by high evapotranspiration and rainfall/runoff rates, while the contribution of groundwater inflow is of about 0.4 Mm³ yr⁻¹. Luzon et al. [16] stated the groundwater contribution to the lake water budget is about 10 time larger than CHE [17], with evaporation the only output term of the lake water budget.

One creek flows to the Zaida Lake (Cañada Real), and three to the Gallocanta Lake (Royo Creek, Acequia Madre, and Pozuelos Creek). The runoff of these creeks has been estimated in this study following the Curve Number method [18]. The estimated runoff (Table 1) was validated through the water level observed at two gauging points in the Royo Creek and Acequia Madre. Part of the flow infiltrates through the creek’s riverbed.

Table 1 Area, Precipitation, Potential Evapotranspiration (PET), Curve Number (CN) and runoff in the creeks under steady-state conditions

2.4 Hydrodynamics and hydrogeological characteristics

Thirty six boreholes, along with permeability tests, and hydrogeological information (from surface geology investigations and previous investigations) provided information on aquifer lithology, thickness and hydraulic conductivities for the groundwater body. The study area hosts a multi-layer aquifer system, and the main aquifers consist of carbonate materials from the Mesozoic and detrital sediments from Late Cenozoic-Quaternary [16]. From top to bottom the main aquifers may be identified (Figs. 3 and 4).

The Tertiary-Quaternary shallow aquifer develops in the central part of the basin and the surroundings of the Gallocanta Lake (on average less than 40 m thick, and minimum thickness of about 4/5 m around the lake). It shows moderate to high permeability (with values of horizontal hydraulic conductivity ranging from 50 to 0.01 m d⁻¹). The Quaternary aquifer is in direct hydraulic connection with and channels the groundwater flows from other aquifers to the Gallocanta Lake. Groundwater flow is relatively radial towards the lagoon, but given the shape of the basin, the main flow direction is from west to east [20].

Mesozoic Units represent the main carbonate regional aquifers (Upper Cretaceous, Jurassic, and Triassic) of the area and also extend beyond the limits of the Gallocanta Lake watershed. The unconfined Upper Cretaceous aquifer extends along the western part of the study area. It is mainly composed by limestone and dolomites, so its permeability is moderate to high, similarly to that of the Tertiary-Quaternary shallow aquifer. Its thickness ranges between 200–300 m. The Jurassic aquifer, extending along the western part of the study area, can be considered as the main aquifer because of its hydraulic conductivity (the highest of the area), and its thickness (~ 250 m). The Lower Cretaceous unit (Utrillas facies), with moderate to low permeability compared to the Mesozoic aquifers, may be considered as an aquitard. It separates the high-permeability aquifers (Jurassic and Upper Cretaceous) along the western part of the study area.

The Triassic claystones from the Keuper facies, extend beneath the Quaternary and the Mesozoic aquifers, and show low permeability. A series of thrusts, in which the Keuper acts as the regional detachment level, halt the groundwater flow in some sectors (Figs. 3 and 4). Finally, the Triassic carbonates of the Muschelkalk Unit only outcrops as a thin band adjacent to the Paleozoic rock at the eastern part of the model domain, so scarce experimental hydrogeological information is available. This is a deep confined aquifer with locally piezometric level higher than that of the Gallocanta Lake [16].

The Paleozoic Unit constitutes the low permeability limit on the north-eastern side and at the bottom of the system. At the contact between the Paleozoic and the Triassic units along the foothills of the Santa Cruz Sierra some springs discharge enough water to supply water to small villages.

2.4.1 Boundaries of the Gallocanta Groundwater Body (GGB) and sink/source terms

The hydrogeological understanding of the GGB boundaries with the contiguous Jiloca River Basin and Piedra/Ortiz Rivers Basin undergone several changes through the years. The hydrogeological conceptual model of the adjoining High Jiloca Valley, prepared by SG DGOH [21] (Fig. 3), considered groundwater inflow from the adjacent GGB. Lately, geophysical works detected a change in the structural setting at depth of the large Mesozoic aquifers in the surroundings of the Caminreal thrust [22, 23] (Fig. 3). That is, an obstacle to the regional flow of the upstream part of the GGB. DGOH and CHE [24] estimated the Caminreal Springs (Ojos de Caminreal), a natural upwelling zone discharging from the GGB Mesozoic aquifer, average flux in about 9 Mm³ yr⁻¹ into the Jiloca Basin.

On the north-western part of the domain, groundwater outflows from the Jurassic and Cretaceous units of the GGB to the Piedra-Ortiz Rivers basin, and drains to the Cimballa Springs (Ojos de Cimballa; Fig. 3). Groundwater discharge into Ortiz and Piedra rivers was estimated to be 39.6 Mm³ yr⁻¹ using the hydrograph-separation analysis [25]. River flow hydrograph data are available for both rivers at gauging stations 9008 and 9129, operated the Ebro Hydrographic Confederation (www.chebro.es).

In 2003, the CHE defined the GGB, as it is now conceived, according to the EU Water Framework Directive (2006/118/EC), and set the first groundwater exploitation rules in the area, thanks to a first steady-state mathematical model implemented using the USGS MODFLOW code [17]. The CHE [17] model considered all the GGB model boundaries as no-flow, so no groundwater outflow to the Jiloca and the Piedra-Ortiz river basin was computed. Outflow from the GGB system was then only due to direct evaporation from the lake and evapotranspiration in the surroundings. The present GGB boundaries are included within those of the Gallocanta watershed (Fig. 3). Orellana-Macias et al. [20] reported the boundaries of the GGB roughly coincide at the eastern and southern side with the watershed boundaries. On the other hand, the western and northern boundaries are difficult to define due to the absence of physical boundaries [17].

However, it is known that a groundwater divide in the Mesozoic aquifers roughly stands along a line dividing the Gallocanta and the Zaida lakes (Fig. 3).

For what concerns sink/source terms in the GGB, rainfall directly recharges the Quaternary and Tertiary aquifers and the Cretaceous and Jurassic ones at the outcrops. Flows from the Cretaceous and Jurassic aquifers near to the lagoon also recharge the Quaternary and Tertiary ones. Loosing ephemeral streams also recharges the Cretaceous and the Triassic aquifers.

The Triassic aquifer discharges to springs and recharge the Quaternary aquifer through lateral flows. At the same time, the Jurassic aquifer laterally discharges to the Cretaceous and the Quaternary ones. Groundwater directly reaches the lagoon near the north-west shoreline. The Cretaceous aquifer also discharges by means of lateral flows to the Quaternary aquifer. Part of the groundwater is pumped locally.

Groundwater in the GGB mainly serves as: i) the main freshwater source for the villages in the surroundings of the Gallocanta Lake, and ii) irrigation source for several irrigated plots on the south and west shore of the lake. Pumping also supplies water to several pig farms across the study area.

2.5 Groundwater flow numerical model implementation

We implemented a steady-state groundwater flow numerical model for an active model domain extending over an area of 946 km² and including all of the GGB, most of the Gallocanta watershed, and parts of the adjacent Piedra-Ortiz and Jiloca watersheds (Fig. 5). We used the FREEWAT platform [26] and MODFLOW-2005 [27] as modelling code. MODFLOW solves the groundwater flow equations in three dimensions using a finite difference scheme:

$$\frac{\partial }{\partial x}({K}_{xx}\frac{\partial h}{\partial x})+\frac{\partial }{\partial y}({K}_{\mathit{yy}}\frac{\partial h}{\partial y})+\frac{\partial }{\partial z}({K}_{\mathit{zz}}\frac{\partial h}{\partial z})+W={S}_{s}\frac{\partial h}{\partial t}$$

where K_xx, K_yy, and K_zz, are values of hydraulic conductivity along the x, y and z axes, which are assumed as parallel to the major axes of hydraulic conductivity (L t⁻¹); h is the potentiometric head (L); W is a volumetric flux per unit volume representing sources and sinks of water, being W < 0 for outflows of the system and W > 0 for inflows of the system (t⁻¹); S_s is the specific storage of the porous material (L⁻¹); and t is time (t). In our case, being the model in steady-state, the right side of the equation is equal to zero, hence, values for the aquifer storage were not implemented.

Fig. 5

Model domain (in the figure model layer 1 is presented)

FREEWAT is a free and open-source composite plugin for QGIS and it allows performing pre-processing, post-processing and spatial analysis of the modelling results in a Geographic Information System. The FREEWAT modelling platform includes a whole suite of USGS modelling codes, such as MODFLOW, and it has been applied to a number of cases and hydrogeological conditions [28, 29].

The model domain is horizontally discretised using 250 × 250 m cells for a whole number of 171 rows and 88 columns (Fig. 5). The model is vertically discretised by means of eight vertical layers. Table 2 includes the initial implemented values of the hydrogeological parameters of the units, adapted from previous works [17]. Figures 5 and 6 show the distribution of hydraulic conductivity for each model layer. Because hard-rock aquifers are simulated in the model, we used the so-called equivalent porous medium approach considering their mainly fissured and fractured porosity as same as the intergranular one assumed by the MODFLOW code [30].

Table 2 Hydrodynamics parameters for the implemented hydrogeological units. Initial data adapted from CHE (2003) [17]

Fig. 6

Conceptual model and MODFLOW packages/modules (in yellow) used to simulate the groundwater flow terms at the GGB. Red and blue arrows represent outflows and inflows to the GGB, respectively

The use of the equivalent porous media approach in hard-rock aquifers (with specific reference to karst aquifers) has the following limitations [31]: i) it fails to capture groundwater hydrodynamics on a local scale, and ii) it fails to simulate transient turbulent flow through conduit network and its interaction with the matrix. Both two conditions are not relevant to our analyses.

Figure 5 shows the model boundaries and sink/source terms implemented in the model and the packages/modules used to represent them. In this study, we hypothesize two main outflows from the Gallocanta watershed. They are the above mentioned: i) outflow to the Piedra-Ortiz Rivers Basin (implemented by means of a General Head Boundary Package, GHB, set in layers 1 to 7, north-western side of the model), and ii) outflow towards the Jiloca River Basin at the Caminreal Springs (represented by means of drain cells in layer 2, Drain Package, east of the model; Fig. 5). All the other model boundaries are set as no-flow boundaries. We also set a no-flow boundary at the bottom model layer, assuming no-flow occurs between the bottom model layer and the lower hydrogeological structures.

Several MODFLOW packages have been used to represent the sink and source terms in the model area (Fig. 6). The MODFLOW Recharge Package mimicked rainfall direct recharge to the model domain. The recharge rates varied spatially and were applied depending on the infiltration characteristics of the various outcropping units. The recharge rate values come from CHE [17] and range from 0.0274 mm d⁻¹ in the Tertiary and Lower Cretaceous outcrops, to 0.0411 mm d⁻¹ in the Quaternary and 0.093 mm d⁻¹ in the Muschelkalk Unit. Recharge was considered null at the Keuper and Paleozoic outcrops. The infiltration rate at the Jurassic limestones (high hydraulic conductivity) and Upper Cretaceous (moderate hydraulic conductivity) outcrops were set at 0.496 and 0.225 mm d⁻¹, respectively. Those parameters were calculated by means of a hydrological balance validated on head levels from the monitoring network of the aquifers.

Three creeks exchanging with the Gallocanta Lake have been simulated using the MODFLOW Streamflow-Routing (SFR2) Package [32]: El Royo, Acequia Madre and Rambla de Pozuelos, and the Cañada Real flowing to the Zaida Lake. The SFR2 package allows performing a water balance, distinguishing input or output sections of the creeks.

This model also considers the connection between the aquifers and the lakes through the Lake package [33] (LAK). This package is used to estimate the water balance of the Zaida Lake and the Gallocanta Lake, and also simulates potential groundwater exchanges between the aquifers and both lakes. Table 3 shows the main parameters used in applying the LAK package. Direct precipitation on the lakes is assumed to be 1 mm d⁻¹ (AEMET 1970–2020), whereas the outflow through evaporation is considered as 5 mm d⁻¹, as per CHE [17]. Evaporation is the calculated by means of the LAK package by multiplying the evaporation rate and the area of the cells represented as lake.

Table 3 Main parameters used in applying the LAK package

Pumping wells were simulated by means of the Well Package. We considered 29 wells out of 99, representing 80% of the licensed yearly abstraction volumes, according to the official water volume licensed by water administrations.

For the calibration process, we used 31 head points from the official groundwater monitoring network and the MODFLOW Head Observation Package. This observation data is available at the CHE website (www.cheebro.es). The average head level for 1979–2020 was used and compared to the heads simulated by the model. Calibration was run in trial-and-error mode in order to focus on the hydrological functioning of the system. As shown in Table 2, we calibrated K_xx and K_yy parameters. Based on data coming from experimental fieldwork and literature data, we modified the hydraulic conductivity values with special regards to highly permeable layers (such us the Quaternary sediments and the Jurassic rocks). Those parameters were overestimated in previous models, so they had to be calibrated to get head levels best fit.

Once the model was validated, in order to verify and to quantify IGFs, the code Zone Budget [34] was implemented and run allowing calculation of zonal groundwater budgets. The model domain was divided into zones (Fig. 7) representing the Piedra-Ortiz Watershed (blue), the Gallocanta Watershed (orange), the Jiloca Watershed (yellow).

Fig. 7

Distribution of calibrated hydraulic conductivity values (K_xx) in model layers 2 to 7 (data for model layer 1 are presented in Fig. 5)

Additional simulations under steady-state conditions were then performed to analyse the changes in the water budget in the following cases: (1) a scenario in which groundwater pumping would be no longer allowed, and (2) a potential drought scenario, simulated by a 20% reduction in rainfall, in turn generating an analogue reduction in recharge and runoff. This reduction is a proxy to the potential consequence of climate change impacts on groundwater resources in the Mediterranean area [35].

3 Results and discussion

The analyses performed provided the whole model and the zonal groundwater budget for the steady state simulation. In the following section, we discuss such budgets along with model calibration.

3.1 Model validation

The analysis of the residuals showed a Mean Absolute Error of 4.2 m and a Root Mean Square Error of 5.74. The value of the Nash–Sutcliffe Efficiency Index, calculated using observed and simulated head levels at the 31 points, was 0.8, a reliable value according to Anderson et al. [29] for the characteristics of our model and our dataset. Figure 8 shows the map of heads residuals. The Zaida Lake area, represented by means of the LAK package, is the less calibrated in the investigated domain for all the Mesozoic and the Triassic aquifers. Head levels of the Quaternary aquifer show very limited changes.

Fig. 8

Water budget zones and results within the model domain. The balance includes recharge, pumping, surface- and groundwater exchanges, and groundwater outflow between the Gallocanta Basin, the Piedra-Ortiz Basin and the Jiloca Basin in Mm³ yr^-1 . The water balance was simulated in three scenarios: a) current scenario; b) no pumping; c) 20 % decrease in recharge, runoff and precipitation. Blue arrows represent inflows to the Gallocanta Lake watershed, whereas orange arrows represent outflows

Regarding head levels in the Upper Cretaceous, Jurassic and Muschelkalk units, simulation results fit well, with values along the 1:1 regression line (Fig. 9). The relation between all of the observed heads and simulated is statistically significant (p < 0.05), whereas according to R², 81% of the variance is justified by the model.

Fig. 9

Groundwater head residual map

3.2 Hydraulic head maps

We present here the groundwater head maps for the Jurassic and the Cretaceous aquifers. The equipotential lines of the Jurassic aquifer under steady-state conditions show a large hydraulic gradient in the GGB domain and map the presence of a groundwater divide west of the Gallocanta Lake (Fig. 10a). Northwards, groundwater flows to the Piedra-Ortiz system, and in the south, the hydraulic gradient increases towards the Caminreal Springs.

Fig. 10

Simulated and observed heads. Dashed grey line represent 1:1 line

The presence of some relevant discrepancies in the simulated heads tapping the Jurassic aquifers in the surroundings of Zaida Lake may not be justified by changes in the hydraulic conductivity of the medium. Because of the geological context, we inferred and tested in the model implementation the presence of a deep thrust that determines the groundwater flow of the area by isolating the northern part of the Jurassic aquifer. This deep thrust in the surroundings of the Zaida Lake should block the hydraulic connection and the flow direction. The thrust was simulated by reducing the hydraulic conductivity values along a line representing the thrust in the Jurassic unit model layer (model layer 5; Fig. 11).

Fig. 11

Equipotential lines in (a) the Jurassic Unit and (b) the Upper Cretaceous Unit

The flow pattern in the Upper Cretaceous aquifer (Fig. 10b) is similar to the Jurassic one in the northern part of the model, but with higher gradient. The groundwater divide is also located to the west of the Gallocanta Lake in a similar position to that of the Jurassic aquifer. Flow towards the Gallocanta Lake occurs. This flow channels through the Quaternary aquifer in the surroundings of the lake. In this case, the effect of the deep inferred thrust in the surrounding of Zaida Lake is significantly lower.

In both the Jurassic and the Cretaceous Units, the groundwater divide is located about 2.5 km to the south of the present GGB northern boundary. In this area, the simulated potentiometric lines tapping the Jurassic Unit and, especially those tapping the Cretaceous Unit, draw a 5.5 km length sector with a NW–SE orientation and very low piezometric gradient (< 0.02%). This area forms a dome with a wide apex, and within it the groundwater divide may move depending on recharge/discharge conditions. A relevant and not negligible outflow may be highlighted towards the Piedra-Ortiz Rivers Basin.

3.3 Water balance

We simulated zonal budget by dividing the model domain in three zones (Fig. 7) representing the Piedra-Ortiz Watershed (blue), the Gallocanta Watershed (orange), the Jiloca Watershed (yellow). In the Gallocanta watershed area, rainfall recharge is 5.74 Mm³ yr⁻¹, whereas 0.89 Mm³ yr⁻¹ (15.5% of the recharge) are pumped for irrigation and urban supply. A 70% of this recharge (4.05 Mm³ yr⁻¹) is transferred to the Piedra-Ortiz System, and 20% (1.19 Mm³ yr⁻¹) to the Jiloca Basin (Fig. 7 and Table 4). We compared the simulated discharge value to the Piedra-Ortiz Basin considering 12% of the outcrops from the Piedra, Ortiz and Gallocanta Basins are within our model domain. As such, we assume a proportional discharge to the Piedra-Ortiz basin, which is 4.7 Mm³ yr⁻¹. This value is similar to the simulated one (4.05 Mm³ yr⁻¹).

Table 4 Simulated and estimated groundwater outflow from the GGB to the Piedra-Ortiz Basin and the Jiloca Basin

Concerning the water balance of the lakes (Table 5), apart from precipitation and runoff, the simulation quantifies a small groundwater flux to the Gallocanta Lake (0.02 Mm³ yr⁻¹), which only represents 1.4% of the inflows to the lake. On the other hand, no groundwater flow reaches Zaida Lake, but outflows from the lake to the aquifers represent 23% of the total outflux. Additionally, according to the model, using the LAK package, we estimated the volume of the Gallocanta and the Zaida Lake in the whole simulated period to be 0.58 and 0.07 Mm³, respectively. The piezometric level of the Quaternary aquifer is directly related to the surface water body of Gallocanta Lake. Regarding the connection between the creeks and the aquifers, the model simulated the same flux (0.02 Mm³ yr⁻¹) of inflows and outflows respectively.

Table 5 Water balance of the Gallocanta and Zaida Lakes for a) current scenario, b) no pumping scenario, and c) 20% decrease in recharge, runoff and precipitation

3.4 Simulated scenarios

Aside from the base case scenario, we simulated two alternative scenarios to compare the current water management and hydrogeological conditions with potential future changes in water management and new environmental contexts associated with incoming climate change. In scenario b (no pumping), outflows to the Piedra-Ortiz and Jiloca Basins increase (4.55 and 1.26 Mm³ yr⁻¹, respectively), whereas the water balance of the creeks indicate higher runoff and inflows to the Gallocanta lake from the creeks (1.27 Mm³ yr⁻¹). Additionally, the absence of pumping also implies higher groundwater fluxes to the Lake (0.07 Mm³ yr⁻¹). These higher inflows lead to higher water volume in the lake (0.72 Mm³ yr⁻¹), thus a larger extension of the lake and hence higher direct rainfall to the lake surface (0.34 Mm³ yr⁻¹) and, in turn, higher evaporation (1.72 Mm³ yr⁻¹).

On the contrary, scenario c, simulating a 20% decrease in recharge, runoff and precipitation, would lead to lower groundwater flows to the Piedra-Ortiz and Jiloca Basins (13 and 16% lower, respectively). In addition, pumping would also decrease (18%) due to lower heads, represented by more cells gone dry in the model. However, infiltration from the creeks would remain constant, although they would no longer be fed by the aquifers, also due to lower heads. In this scenario, the lakes are especially impacted since inflows from runoff are practically null (Table 1). Finally, simulation’s uncertainty in scenario c is high, mainly due to the uncertainties related to the estimation of evaporation. This is due to the relation between the water volume and the lake’s surface. This relation will need to be specifically addressed in future simulations under transient conditions.

4 Discussion

Data in Table 4 highlight the persistence and scope of the questions related to the assessment of the IGF occurring from the Gallocanta Lake watershed boundaries. The model presented was successfully validated thanks to the outflows to the Piedra-Ortiz Basin, which comply with the estimations of the rivers’ baseflow, and to the Jiloca Basin, also in agreement with the estimated flow at the Caminreal Springs. Depending on the conceptual model and the methodology used, the rate of IGF ranges from 0 to 8 Mm³ yr⁻¹ to the Piedra-Ortiz Basin and from 0 to 15 Mm³ yr⁻¹ to the Jiloca Basin, the latter varying for more than one order of magnitude. The large IGF values presented in IGME [36] may be due to the use of data from a very wet period, a larger investigated domain, and a methodology using surface geology and groundwater head maps for evaluating the water budget. DGOH and CHE [24] did not considered transfer to the Piedra-Ortiz Basin. Luzon et al. [13] did not specify the methodology used (probably an analytical calculation), and CHE [17] choose no-flow boundaries for the numerical modelling, hence not allowing any outflow through the geological media.

Regarding the groundwater discharges to the Jiloca River Basin, the Caminreal thrust has to be taken into account, since it performs a key role in the regional hydrogeology. According to the simulation, groundwater from the Gallocanta Basin is about 17% of the Caminreal Spring average flux (7 Mm³ yr⁻¹ according to DGOH and CHE [24]).

The outflows to the Piedra-Ortiz Basin are higher than those to the Jiloca Basin. This is likely due to the structural control of the Jurassic aquifers to the northern part of the Gallocanta Lake, as well as the presence of quasi-impermeable rocks from the Keuper unit, and the topographical control on groundwater flow. Although the hydraulic connection between the surroundings of Zaida Lake and Gallocanta Lake is not completely interrupted, the flux towards Gallocanta Lake is severely hindered. The potential lines of the Jurassic aquifer to the north of Gallocanta Lake follow a different pattern and have lower heads. This pattern was already addressed by CHE [17] in the previous steady-state model.

IGF is an important hydrological process to the understanding of water and chemical fluxes and budgets for water resources planning and management [3, 37, 38]. For instance, groundwater coming to an aquifer from an adjoining mountain-block recharge area [39] may constitute a relevant part of the storage to be exploited or a relevant baseflow component in rivers external to that area [1, 40], as it is in our case. This additional storage has to be taken into account in the design of Integrated Water Resources Management Programs when aquifers extend across watersheds.

The delineation of groundwater bodies boundaries is not an easy task [41]. EU [42] suggests delineating groundwater bodies in such a way that any groundwater flow from one groundwater body to another (a) is so minor that it can be ignored in water balance calculations; or (b) can be estimated with adequate precision. Moreover, in hard-rock aquifers not well-known geological structures may heavily alter groundwater flow dynamics and hence the delineation of groundwater bodies boundaries [43, 44]. In our case, we inferred the presence of geological structure in the form of a deep thrust, having a large impact on the groundwater flow at depth, but not at shallow levels. The thrust needs to be confirmed by new investigations (i.e. geophysical studies).

The performed simulations provide key elements in addressing the areal extent of groundwater bodies and the hydraulic relationships with adjoining bodies. In the investigated case, the presence of an internal (to the presently defined GGB) moving groundwater divide (due to meteo-climatic variability and climate change), and large groundwater outflows to the two adjoining watersheds of Piedra-Ortiz on the west and Jiloca River on the east, makes the areal extent of the GGB inconsistent to the EU definition of groundwater body [42]. This must be taken into account to revise the current delimitation and the laws associated with it. Not in vain, previous studies addressed the lack of efficiency of the environmental programs due to the failure when delimiting both groundwater bodies and protected areas [45].

The small groundwater inflows to the Gallocanta Lake (1.4%) shows the Lake high dependence on runoff and rainfall, and justifies the low nitrate concentration of the lake. Given the current quality status of the aquifers, if groundwater inflows were higher, pollution could be a severe issue in the Gallocanta Lake. This because the only outflow from the lake is evaporation, and water renewal is practically null, which simultaneously enhance salinization. In the regional hydrogeological context, groundwater fluxes towards the Gallocanta Lake are relatively low, although they are essential for the ecological system associated with the wetland. Those fluxes should be specifically analysed under transient modelling conditions.

5 Conclusions

In this paper we numerically quantified at the Gallocanta watershed (Spain) the importance of interbasin groundwater flow, the presence of a groundwater divide, and the inferred influence of geological structures on groundwater flow dynamics. We then highlight the following conclusions:

• 1. Interbasin groundwater flow is an important term to be considered in the water balance of the studied watersheds. In our case, except for a small flux from the Upper Cretaceous aquifer to the Gallocanta Lake, most of the groundwater flux from the Jurassic and the Upper Cretaceous aquifers flows to the Piedra-Ortiz Basin and the Jiloca Basin. The groundwater outflow may also imply a significant flow of pollutants since both aquifers have high nitrate concentrations. Not in vain, most of the GGB is declared as a Nitrate Vulnerable Zone.

• 2. The model indicates a similar water divide under steady-state conditions for both the Jurassic and the Upper Cretaceous aquifers, confirming the hydraulic connection of the aquifers. However, this water divide must be considered non-fixed depending on recharge/discharge conditions. We then propose, for groundwater resource management purpose, a unique groundwater body, spanning from the Caminreal Springs on the east and the Piedra-Ortiz basin on the west, hence including the groundwater divide internal to the endorheic Gallocanta watershed. Monitoring the baseflow of the Piedra-Ortiz river and of the Caminreal Springs will allow to get information on the evolution of the groundwater resource availability in the Gallocanta watershed.

• 3. Our model is a dynamic, growing tool in order to manage the groundwater resource in the Gallocanta Lake. We calibrated the model using head data from different aquifers levels. Although, our set of parameters is not unique, we identified areas where further investigations should be run in order to further test the conceptual model under varying hydrological conditions.

• 4. Finally, the results of our work show that the definition of groundwater body boundaries according to present EU regulations for protecting the groundwater resources should be based on thorough analyses of the geological structures and hydrogeological information. Numerical modelling analyses, especially in areas where scarce data are available, could provide valuable help in testing conceptual model hypothesis and inferring information on the presence of relevant structural features driving the groundwater flow.

Nowadays, the easy availability of free and open-source numerical codes, graphical user interfaces and educational/training material (FREEWAT [26]; ModelMuse [46]) gives no excuses for not applying such tools in groundwater resource analyses and management.