Skip to main content

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


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 Mm3 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
figure 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−1 (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
figure 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
figure 3

Geological setting of the Gallocanta Lake watershed

Fig. 4
figure 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 km2. Luzon et al. [13] reports salinity values ranging from 3600 to 136,000 µS cm−1. 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 Mm3 yr−1. 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−1). 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 Mm3 yr−1 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 Mm3 yr−1 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 (

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 km2 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 Kxx, Kyy, and Kzz, 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−1); 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−1); Ss is the specific storage of the porous material (L−1); 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
figure 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
figure 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−1 in the Tertiary and Lower Cretaceous outcrops, to 0.0411 mm d−1 in the Quaternary and 0.093 mm d−1 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−1, 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−1 (AEMET 1970–2020), whereas the outflow through evaporation is considered as 5 mm d−1, 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 ( 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 Kxx and Kyy 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
figure 7

Distribution of calibrated hydraulic conductivity values (Kxx) 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
figure 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 Mm3 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 R2, 81% of the variance is justified by the model.

Fig. 9
figure 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
figure 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
figure 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 Mm3 yr−1, whereas 0.89 Mm3 yr−1 (15.5% of the recharge) are pumped for irrigation and urban supply. A 70% of this recharge (4.05 Mm3 yr−1) is transferred to the Piedra-Ortiz System, and 20% (1.19 Mm3 yr−1) 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 Mm3 yr−1. This value is similar to the simulated one (4.05 Mm3 yr−1).

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 Mm3 yr−1), 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 Mm3, 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 Mm3 yr−1) 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 Mm3 yr−1, respectively), whereas the water balance of the creeks indicate higher runoff and inflows to the Gallocanta lake from the creeks (1.27 Mm3 yr−1). Additionally, the absence of pumping also implies higher groundwater fluxes to the Lake (0.07 Mm3 yr−1). These higher inflows lead to higher water volume in the lake (0.72 Mm3 yr−1), thus a larger extension of the lake and hence higher direct rainfall to the lake surface (0.34 Mm3 yr−1) and, in turn, higher evaporation (1.72 Mm3 yr−1).

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 Mm3 yr−1 to the Piedra-Ortiz Basin and from 0 to 15 Mm3 yr−1 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 Mm3 yr−1 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. 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. 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. 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. 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.

Availability of data and materials

All data generated or analyzed during this study will be made available upon request.


  1. Pellicer-Martinez F, Martinez-Paz JM. Assessment of interbasin groundwater flows between catchments using a semi-distributed water balance model. J Hydrol. 2014;519:1848–58.

    Article  Google Scholar 

  2. Le Mesnil M, Charlier JB, Moussa R, Caballero Y, Dorfliger N. Interbasin groundwater flow: Characterization, role of karst areas, impact on annual water balance and flood processes. J Hydrol. 2020;585:124583.

    Article  Google Scholar 

  3. Frisbee MD, Meyers ZP, Stewart-Maddox NS, Caffee MW, Bogeholz P, Hughes MN. What is the source of baseflow in agriculturally fragmented catchments? Complex groundwater/surface-water interactions in three tributary catchments of the Wabash River, Indiana, USA. Hydrol Process. 2017;31:4019–38.

    Article  Google Scholar 

  4. Belcher WR, Bedinger MS, Back JT, Sweetkind DS. Interbasin flow in the Great Basin with special reference to the southern Funeral Mountains and the source of Furnace Creek springs, Death Valley, California, U.S. J Hydrol. 2009;369:30–43.

  5. Alvarez-Campos O, Olson EJ, Welp LR, Frisbee MD, Zuniga Medina SA, Diaz Rodriguez J, et al. Evidence for high-elevation salar recharge and interbasin groundwater flow in the Western Cordillera of the Peruvian Andes. Hydrol Earth Syst Sci. 2022;26:483–503.

    Article  Google Scholar 

  6. Boutt DF, Corenthal LG, Moran BJ, Munk L, Hynek SA. Imbalance in the modern hydrologic budget of topographic catchments along the western slope of the Andes (21–25°S): implications for groundwater recharge assessment. Hydrogeol J. 2021;29:985–1007.

    Article  Google Scholar 

  7. Danapour M, Hojberg AL, Jensen KH, Stisen S. Assessment of regional inter-basin groundwater flow using both simple and highly parameterized optimization schemes. Hydrogeol J. 2019;27:1929–47.

    Article  Google Scholar 

  8. De Vargas T, Boff FE, Belladona R, Faccioni LF, Reginato PAR, Carlos FS. Influence of geological discontinuities on the groundwater flow of the Serra Geral Fractured Aquifer System. Groundw Sustain Dev. 2022;18:100780.

    Article  Google Scholar 

  9. De Schepper G, Therrien R, Refsgaard JC, Hansen AL. Simulating coupled surface and subsurface water flow in a tile-drained agricultural catchment. J Hydrol. 2015;521:374–88.

    Article  Google Scholar 

  10. Joodavi A, Izady A, Karbasi Maroof MT, Majidi M, Rossetto R. Deriving optimal operational policies for off-stream man-made reservoir considering conjunctive use of surface- and groundwater at the Bar dam reservoir (Iran). J Hydrol Reg Stud. 2020;31:100725.

    Article  Google Scholar 

  11. Le Moine N, Andreassian V, Perrin C, Michel C. How can rainfall-runoff models handle intercatchment groundwater flows? Theoretical study based on 1040 French catchments. Water Resour Res. 2007;43:W06428.

    Article  Google Scholar 

  12. Serrano-Notivoli R, de Luis M, Beguería S. An R package for daily precipitation climate series reconstruction. Environ Modell Softw. 2017;89:190–5.

    Article  Google Scholar 

  13. Luzon A, Perez A, Sanchez JA, Soria AR, Mayayo MJ. Evolution from a freshwater to saline lake: a climatic or hydrogeological change? The case of Gallocanta Lake (northeast Spain). Hydrol Process. 2007;21:461–9.

    Article  Google Scholar 

  14. Guimera J. Structure of an intraplate fold-and-thrust belt: The Iberian chain. A synthesis. Geol Acta. 2018;16:427–38.

    Google Scholar 

  15. Schutt B. Reconstruction of Holocene paleoenvironments in the endorheic basin of Laguna de Gallocanta, Central Spain by investigation of mineralogical and geochemical characters from lacustrine sediments. J Paleolimnol. 1998;20:217–34.

    Article  Google Scholar 

  16. Perez A, Luzon A, Roc AC, Soria AR, Mayayo MJ, Sanchez JA. Sedimentary facies distribution and genesis of a recent carbonate-rich saline lake: Gallocanta Lake, Iberian Chain, NE Spain. Sediment Geol. 2002;148:185–202.

    Article  Google Scholar 

  17. CHE. Establishing the Exploitation Rules of the Gallocanta Hydrogeologic Unit and Delineation of the Protection Perimeter of the Gallocanta Lake. Zaragoza: Ebro Hydrographic Confederation; 2003 [in Spanish].

  18. SCS USDA. National Engineering Handbook, Section 4. Hydrology. Washington, DC: Soil Conservation Service, United States Department of Agriculture; 1972.

  19. Hargreaves GH, Samani ZA. Reference crop evapotranspiration from temperature. Appl Eng Agric. 1985;1:96–9.

    Article  Google Scholar 

  20. Orellana-Macias JM, Merchan D, Causape J. Evolution and assessment of a nitrate vulnerable zone over 20 years: Gallocanta groundwater body (Spain). Hydrogeol J. 2020;28:2207–21.

    Article  Google Scholar 

  21. SG DGOH. Delimitation of the Hydrogeological Units in the Iberian Peninsula and Balearic Islands. Summary of Their Characteristics. Madrid: Dirección General de Obras Hidráulicas; 1988 [in Spanish].

  22. SG DGOH. Geophysical Study of the Gallocanta Area (Zaragoza). Madrid: Dirección General de Obras Hidráulicas; 1988 [in Spanish].

    Google Scholar 

  23. SGOP. Test Cores and Piezometers in the Upper Basin of the Jiloca River. Zaragoza: Servicio Geológico de Obras Públicas; 1995 [in Spanish].

  24. DGOH and CHE. Analysis of the hydraulic resources of the aquifers within the Zaragoza province. Madrid: Dirección General de Obras Hidráulicas, Zaragoza: Ebro Hydrographic Confederation; 1990 [in Spanish].

  25. Barlow PM, Cunningham WL, Zhai T, Gray M. U.S. Geological Survey Groundwater Toolbox, a Graphical and Mapping Interface for Analysis of Hydrologic Data (Version 1.0)—User Guide for Estimation of Base Flow, Runoff, and Groundwater Recharge From Streamflow Data. In: Techniques and Methods 3-B10. Reston: United States Geological Survey; 2015.

  26. Foglia L, Borsi I, Mehl S, De Filippis G, Cannata M, Vasquez-Sune E, et al. FREEWAT, a free and open source, GIS-integrated, hydrological modeling platform. Ground Water. 2018;56:521–3.

    Article  Google Scholar 

  27. Harbaugh AW. MODFLOW-2005: the U.S. Geological Survey Modular Ground-Water Model--the Ground-Water Flow Process. In: Techniques and Methods 6-A16. Reston: United States Geological Survey; 2005.

  28. De Filippis G, Pouliaris C, Kahuda D, Vasile TA, Manea VA, Zaun F, et al. Spatial data management and numerical modelling: demonstrating the application of the QGIS-integrated FREEWAT platform at 13 case studies for tackling groundwater resource management. Water. 2020;12:41.

    Article  Google Scholar 

  29. Cannata M, Neumann J, Rossetto R. Open source GIS platform for water resource modelling: FREEWAT approach in the Lugano Lake. Spat Inf Res. 2018;26:241–51.

    Article  Google Scholar 

  30. Anderson MP, Woessner WW, Hunt RJ. Chapter 1: Introduction. In: Anderson MP, Woessner WW, Hunt RJ, editors. Applied Groundwater Modeling: Simulation of Flow and Advective Transport. 2nd ed. Cambridge: Academic Press; 2015.

  31. Ghasemizadeh R, Yu X, Butscher C, Hellweger F, Padilla I, Alshawabkeh A. Equivalent Porous Media (EPM) simulation of groundwater hydraulics and contaminant transport in karst aquifers. Plos One. 2015;10:e0138954.

    Article  Google Scholar 

  32. Niswonger RG, Prudic DE. Documentation of the Streamflow-Routing (SFR2) Package to Include Unsaturated Flow Beneath Streams - A Modification to SFR1. In: Techniques and Methods 6-A13. Version 1.2, revised Aug 2009 ed. Reston: United States Geological Survey; 2005.

  33. Merritt ML, Konikow LF. Documentation of a Computer Program to Simulate Lake-Aquifer Interaction Using the MODFLOW Ground-Water Flow Model and the MOC3D Solute-Transport Model. In: Water-Resources Investigations Report 00–4167. Reston: United States Geological Survey; 2000.

  34. Harbaugh AW. A Computer Program for Calculating Subregional Water Budgets Using Results from the U.S. Geological Survey Modular Three-dimensional Finite-difference Ground-water Flow Model. Reston: United States Geological Survey; 1990.

  35. MedECC. Risks Associated to Climate and Environmental Changes in the Mediterranean Region. A Preliminary Assessment by the MedECC Network Science-Policy Interface – 2019. Marseille: Mediterranean Experts on Climate and Environmental Change; 2019.

  36. IGME. Hydrogeological Study of the n. 57 Aquifer System: Monreal-Gallocanta Mesozoic. Madrid: Instituto Geológico y Minero de España. 1981 [in Spanish].

  37. Cardenas MB. Surface water-groundwater interface geomorphology leads to scaling of residence times. Geophys Res Lett. 2008;35:L08402.

    Article  Google Scholar 

  38. Genereux DP, Jordan M. Interbasin groundwater flow and groundwater interaction with surface water in a lowland rainforest, Costa Rica: A review. J Hydrol. 2006;320:385–99.

    Article  Google Scholar 

  39. Markovich KH, Manning AH, Condon LE, McIntosh JC. Mountain-block recharge: a review of current understanding. Water Resour Res. 2019;55:8278–304.

    Article  Google Scholar 

  40. Baudron P, Alonso-Sarria F, Garcia-Arostegui JL, Canovas-Garcia F, Martinez-Vicente D, Moreno-Brotons J. Identifying the origin of groundwater samples in a multi-layer aquifer system with Random Forest classification. J Hydrol. 2013;499:303–15.

    Article  Google Scholar 

  41. Sanchez D, Carrasco F, Andreo B. Proposed methodology to delineate bodies of groundwater according to the European water framework directive. Application in a pilot Mediterranean river basin (Málaga, Spain). J Environ Manage. 2009;90:1523–33.

  42. EU. Common Implementation Strategy for the Water Framework Directive (2000/60/EC). Guidance Document No 2 Identification of Water Bodies. Brussels: European Communities; 2003.

  43. Mutua S, Ghysels G, Anibas C, Obando J, Verbeiren B, Van Griensven A, et al. Understanding and conceptualization of the hydrogeology and groundwater flow dynamics of the Nyando River Basin in Western Kenya. J Hydrol Reg Stud. 2020;32:100766.

    Article  Google Scholar 

  44. Inbar N, Rosenthal E, Magri F, Alraggad M, Möller P, Flexer A, et al. Faulting patterns in the Lower Yarmouk Gorge potentially influence groundwater flow paths. Hydrol Earth Syst Sci. 2019;23:763–71.

    Article  Google Scholar 

  45. Arauzo M, Valladolid M, Martínez-Bastida JJ. Spatio-temporal dynamics of nitrogen in river-alluvial aquifer systems affected by diffuse pollution from agricultural sources: Implications for the implementation of the Nitrates Directive. J Hydrol. 2011;411:155–68.

    Article  Google Scholar 

  46. Winston RB. ModelMuse. version 5.0. Reston: United States Geological Survey; Software Release: 18 Mar 2022. (Accessed 10 Aug 2023).

Download references


This research was funded by the Government of Aragon, grant number C137/2016, by the Spanish Ministry of Science and Innovation—FEDER funds [EU] via the Research Project Agro-SOS, grant number PID2019-108057RB-I00.

Author information

Authors and Affiliations



Arce, M.: research design, conceptual modelling, numerical modelling; Orellana-Macías, J.M.: research design, data gathering, conceptual modelling, numerical modelling, manuscript writing and revising; Causapé, J.: research design, conceptual modelling, supervision; Ramajo, J.: data gathering, conceptual modelling; Gale, C.: data gathering, conceptual modelling; Rossetto, R.: research design, conceptual modelling, numerical modelling, manuscript writing, editing and revising. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Rudy Rossetto.

Ethics declarations

Competing interests

The authors declare they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Arce, M., Orellana-Macías, J.M., Causapé, J. et al. Model-based assessment of interbasin groundwater flow in data scarce areas: the Gallocanta Lake endorheic watershed (Spain). Sustain Environ Res 33, 32 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: