Flow Numerical Modelling in Thermal Karst Systems: The Case of Alhama de Aragón and Jaraba Springs

Abstract

The underground flow of a karstic aquifer within one of Spain and Europe’s most important thermal systems (Alhama and Jaraba thermal springs, with a combined flow rate of 1200 L/s, 711 L/s at more than 30 °C) was simulated. In the simulation process, it was important to consider how temperature (a very sensitive parameter when calibrating the numerical model) and depth influence the variation in hydraulic conductivity in the aquifer. The location of previously unknown high recharge zones was also essential in the calibration. It was verified that some fault jumps break the hydraulic continuity of the aquifer, and the role of most of the existing faults in the regional flow is generally unimportant since they are incapable of explaining by themselves the large volume of water evacuated. It is relevant to highlight the importance of the orientation of the strata when calibrating the model, which become vertical in the area of the outcrops. In the end, the modelled regional flow as well as the simulated groundwater contour lines are consistent with the progressive increase in temperature, the age of the water, the mineralization, the piezometric values measured in the observation wells, and the springs’ flow through which the system discharges. The most significant finding is the validation of the conceptual hydrogeological model through regional flow simulations from numerical models, confirming the recharge area and supporting the inferred origins of the springs.

1. Introduction and Objectives

Modelling karst aquifers that sustain large springs is particularly compelling, yet it ranks among the most complex tasks relative to nearly any other groundwater study due to the multifaceted nature of these aquifers. The flow generated by a large spring that drains an extensive aquifer represents the cumulative result of all processes acting upon the aquifer, both natural and anthropogenic. Developing a physics-based quantitative model can be an invaluable tool for gaining deeper insights into the system’s behaviour and dynamics, supporting future research planning and ensuring effective management and protection of these springs [1].

One of the challenges in modelling the flow in highly karstified and heterogeneous aquifers is simulating two types of flow: one that follows Darcy’s law and another that occurs mainly along preferential flow pathways. Various research software programs are available for simulating dual-porosity aquifers (e.g., [2,3,4,5,6]), as cited in Zanini et al. (2021) [4], who evaluated the effectiveness of hybrid karst aquifer models that couple a single continuous medium with discrete features to accurately represent interactions between a vertical conduit and the surrounding unsaturated matrix.

Focussing on numerical models [7], Modflow-2005, Conduit Flow Process [8], emerged as an example developed to address this issue, providing the capability to simulate various flow types in karst aquifers using an explicitly physics-based approach.

However, our case involves a very extensive and minimally karstified aquifer, with slow flow velocities not exceeding 1 km/year, as verified through tritium isotopes, and where spring hydrographs, for example, appear nearly horizontal [9]. For this reason, the aquifer resembles an intergranular porosity system more than a karstic one. Consequently, we applied Modflow-2005 [7]. The successful simulation results further support this modelling choice.

The application of numerical simulation in thermal karst aquifers includes varied aspects and different scales of study. For example, it helps to understand how the thermal response of the springs is based on the depth of groundwater circulation [10]. For a deep regional flow system without additional heat input, Gong et al. (2019) present a coupled flow–heat–solute numerical model for the quantitative description of the factors that contribute to karstification [11]. They conclude that the conduits enlarge mainly in regions with lower temperature and higher water flow but remain unchanged in regions with higher temperature or lower water flow. Szijártó et al. (2021) carried out numerical simulations to examine the temperature field and flow pattern to highlight the role of different driving forces of groundwater flow in the Buda Thermal Karst (Hungary) [12]. Among other results, and on a basin-wide scale, they indicate that hydraulically conducting faults have only a minor influence on the temperature field and flow pattern.

In our case, a large-scale flow modelling of a large sedimentary basin with thermal water has been undertaken. After defining the conceptual model of hydrogeological functioning [9], the main objective of this study is to numerically model the regional flow of this geothermal system. The aim is to clear up uncertainties about the origin and hydrogeological functioning of the Alhama de Aragón and Jaraba thermal springs and other semi-thermal springs associated with the same geological area. Thus, the geothermal system has been reproduced for the first time by the numerical modelling of groundwater flow in 3D with hydrogeological modelling software. The main objectives of the modelling are the following: 1. To create a mathematical model of the area that validates the conceptual model explained. From this numerical model, the calibrated hydrogeological and geometric parameters are obtained so that they are consistent with the conceptual model. 2. Due to the great extension and complicated nature of the aquifer in question, it has been possible to orient the conceptual model from the calibration process of the numerical model, thus providing answers about the hydrogeological functioning which would not have been possible without the mathematical model.

It must be considered that this thermal aquifer, due to its largely deep location and confinement by the Tertiary cover, could be considered a naturally well-protected thermal water resource, with renewable and extensive resources. However, the nitrate contamination that lurks in some recharge areas could compromise, in the long term, the quality of the water of the spas that supply thermal and mineral water. In addition, the irrigation of large areas by extracting groundwater from this aquifer could compromise its resources in terms of quantity. Therefore, this article also aims to draw the attention of water managers to the challenges and problems associated with the sustainable use of thermal water resources from this aquifer. The assessment, evaluation, and modelling that has been carried out aims to provide a useful basis for the management of thermal water from the main thermal aquifer in Spain.

2. Hydrogeology of the Site

According to [9], the thermal aquifer is made up of about 350 m of Cretaceous limestone and dolomite. This calcareous aquifer extends in depth and in hydraulic continuity under the Tertiary up to 4500 m in the form of a large synclinorium, which outcrops along the edge of the Aragonese Branch of the Iberian Mountains, with the most impermeable Tertiary sediments of the Almazán Basin (Figure 1). The thermal and semi-thermal springs are aligned with this edge, with the most abundant and hottest springs emerging in the lowest areas, all of them belonging to the Ebro basin. The main recharge takes place in the calcareous outcrops of the aquifer of the Duero River basin, so it is a natural underground transfer. The Tertiary groundwater contour lines, however, adapt to the topography. On the other hand, an increase in the temperature, age, and flow of the springs is observed as we move away from the recharge zone.

3. Methodology

3.1. Conceptual Model

• Hydrothermal system domain boundaries: The definition of hydrothermal system domain boundaries was carried out in order to locate the hydrogeological divide that drains the springs located in the Ebro Basin of Alhama and Jaraba, where the boundary of the Cretaceous hydrothermal aquifer extends into the Duero Basin. There are several works and doctoral theses on the hydrogeology of the Cretaceous carbonate rocks located in the Almazán Basin belonging to the Duero Basin, informing us of a behaviour similar to the system studied here, although they are not such deep aquifers; rather, they are large springs with an extensive recharge area [14] and sometimes quite far from the springs and where the flow goes below the Paleogene and Neogene that confines it [15,16].
• Recharge calculation: Although the modelling has made it possible to estimate recharge in detail, the previous evaluations of [17] based on the hydrometeorological records of the study area were used as indicative starting values. These evaluations have been modified, mainly taking into account the variation in rainfall and evapotranspiration with altitude and the different topographical characteristics of the Jalón and Duero Basin.

3.2. Modelling of Hydrothermal System Flow

The hydrogeological conceptual model is known, and the subsurface flow situation was simulated with Modflow. At a regional or basin scale, it is customary to assume that the rock medium behaves as an equivalent porous medium with blocks whose permeabilities adequately reflect the hydraulic behaviour of the flow in the system. This modelling allowed us to test and refine the previous hydrogeological conceptual model. The model was subsequently calibrated in both steady and transient regimes and then subjected to a sensitivity analysis of the main parameters conditioning the behaviour of the system.

The mathematical code used is Modflow-2005, a numerical code that uses the three-dimensional finite difference method to simulate flow in the saturated zone. The flow regime can be stationary or transient. The Modflow-2005 version of Model Muse was used as the graphical analysis interface.

The study area can be represented by a three-dimensional finite difference grid, on which the system of differential equations of groundwater flow can be solved.

According to [9], the domain of the hydrogeological model of the Cretaceous hydrothermal aquifer covers the entire Tertiary of the Almazán Basin from the hydrogeological divide of the Ebro to the peripheral Cretaceous edges of the Aragonese branch of the Iberian Range and the Jurassic–Cretaceous of the Castilian branch. In addition, the domain of the model area is extended to the west, entering the Duero Basin. The hydrogeological divide of the hydrothermal aquifer does not coincide with the Duero–Ebro hydrographic divide, as justified above. The extension and geometry of the model are shown in Figure 2 and Figure 3.

The area shown in the model covers an area of approximately 2295.5 km².

The model considered two differentiated horizontal layers. Layer 1 is the upper layer of sedimentary deposits corresponding to the Tertiary that constitutes the Almazán Basin. It has a low permeability layer with variable thickness, ranging from a few meters in the outcrop area to 4000 m in the central area of the basin. Layer 2 is the lower layer, constituting the Cretaceous limestone aquifer, which is confined by the Tertiary. It has an average thickness of about 300 m, except in the outcrops and surrounding areas, where greater apparent thickness can be observed, reaching up to 800 m in some cases.

The overall dimensions of each cell in both layers were 500 × 500 m. In the areas of greatest interest, such as the edge of outcrops where springs are located or areas with a high density of contour lines, a spatial discretization (mesh refinement) was carried out to achieve greater accuracy in the simulation. The model was made up of a mesh of 49,621 cells, as shown in Figure 4.

The roof and wall elevations of each layer were modelled, matching the geometrical model as much as possible with the morphology of the geological contacts.

Both layers were discretised into three horizontal layers in order to improve the calculation of the software, as seen in Figure 5, Figure 6 and Figure 7.

3.3. Hydrogeological Parameters

The ranges of hydraulic permeability values established in the model were based on the initial horizontal hydraulic permeability values taken from the pumping tests [9] (

Table 1. Calculated vs. observed groundwater heads at calibration points.

Alhama de Aragón Piezometer		Embid de Ariza Piezometer		Deza Piezometer
Calculated Values	Observed Values	Calculated Values	Observed Values	Calculated Values	Observed Values
670	665–667.5	791	777–773	920	919–925

). These initial values were adjusted to calibrate the numerical model. It should be noted that according to the properties of water, a higher temperature produces a reduction in viscosity and therefore an increase in hydraulic conductivity. This translates in the model into a zonal variation in permeability depending on the depth, as shown in Figure 8. In addition, the system would be practically equivalent to an anisotropic homogeneous aquifer with a vertical permeability smaller than the horizontal permeability, according to [18]. As indicated in the introduction and objectives, the aquifer is extensive and minimally karstified, which makes it more like an intergranular porosity system than a karstic one; therefore, in general, a vertical permeability 10 times smaller than the horizontal permeability (K_z = K_xy/10) was defined in the model. This was considered where the strata are horizontal, but as the Cretaceous limestones reach the outcrops of the Aragonese branch, the strata become vertical, and therefore, the permeabilities K_x, K_y, and K_z had to be redefined so that the axes of the strata coincide, as shown in Figure 8.

Therefore, considering the initial values described in the conceptual model and taking into account the considerations described here, the following hydraulic conductivities were calibrated in the model. The hydraulic conductivities fit very well with those defined in the conceptual model.

3.4. Outline Conditions

According to the conceptual model defined in [9], four types of boundary conditions were defined in the model.

Concerning the first boundary condition, there are two types of water recharge by effective infiltration. The first corresponds to direct or autogenous recharge from precipitation recharge. Secondly, there is an indirect or diffuse recharge, coming from the infiltration of the Mesa River as it passes through the Cretaceous limestone between the village of Calmarza and Jaraba, as mentioned above.

The autogenous recharge of precipitation in the field was simulated in the model using the Modflow recharge package (Recharge).

As the study area is relatively large, around 2295.5 km², a spatial distribution of recharge was considered. As initial and indicative values, the previous assessments of [17], based on the hydrometeorological records of the study area, were used as starting values.

As shown in Figure 9, the annual recharge was separated into different recharge areas according to the initial values and the following considerations. This recharge corresponds to each complete hydrological year modelled in the stationary period.

• The higher the altitude, the higher the precipitation, as explained in [19]. Therefore, the northernmost areas of the model, being at a higher altitude, contribute to a higher recharge to the model.
• In addition, there are two poljes of special recharge. The first is Villaseca de Arciel and the head of the Rituerto (considered to be part of the recharge of the Sierra de la Pica and Sierra de Tajahuerce).
• Recharge in the Tertiary of the Duero Basin was considered in the model to be higher than in the Ebro basin. In the Duero Basin, there are flat endorheic areas with lower slopes, which favours recharge and reduces surface runoff. There are also two poljes with periodic flooding.

The direct or allogenic recharge through drains in the bed of the Mesa River as it passes through Calmarza was simulated using the “WELL” package, which allowed us to consider a flow of water leaving or entering a given cell. This water inflow was taken according to the values of [20,21], where it is indicated that this recharge is visualised mainly during low water periods and is on the order of 0.100 m³/s to 0.4 m³/s.

Therefore, a total allogenic recharge rate of 0.2 m³/s was considered, distributed homogeneously along the course of the Mesa River between the villages of Calmarza and Jaraba. This recharge rate was also calibrated by the numerical model.

Once the model was calibrated, the final estimated and calibrated recharge values were as shown in Figure 9.

As previously mentioned, the outflows associated with this model correspond to the discharge areas of the water flowing from the Cretaceous limestone aquifer, which are produced through the high-flow Alhama de Aragón and Jaraba springs, as well as through other springs associated with the same geological area, which have medium to low flows.

Concerning the second boundary condition, these springs were simulated in the numerical model using the Modflow drainage package (DRAIN).

The drainage package was also simulated along the outcrop edges of the Cretaceous limestone outcrops in the areas where the springs occur in order to have a greater degree of detail and not to limit it to only the cells where the springs are located. In this way, a better ability to calibrate the piezometric level was achieved. The drain boundary conditions are shown in detail in Figure 10.

The third boundary condition corresponds to the edge conditions (no flow). The definition of this behaviour was carried out through the boundary conditions imposed on the model. The following boundaries were considered closed in the model.

• Jaraba fault: As explained in [9], the Jaraba fault produces a disconnection in the Cretaceous limestones in a southeasterly direction. This disconnection justifies why there are no thermal water springs in the outcrops of Cretaceous limestone in the town of Ibdes, as these outcrops are at a lower altitude than the Jaraba springs. This edge has been modelled as a no-flow zone.
• Impermeable boundary to the west of the Cretaceous limestone outcrops, corresponding to the end of the Cretaceous limestone aquifer: This boundary was modelled with a no-flow zone.
• Hydrogeological divide boundary to the east of the model: non-zero flow boundary.
• Impermeable boundary to the south, corresponding to the contact with the Sierra del Solorio.

The hydrogeological disconnections are the last boundary condition defined in the numerical model. The Cretaceous, along the edge of the Cretaceous limestones of the Aragonese Branch as it passes through Embid de Ariza and extending southwards, is disconnected by a series of thrusts through which the impermeable materials of the Mesozoic and Palaeozoic substrata have been uplifted. In some cases, these fault jumps cause barriers and a hydraulic disconnection of part of the aquifer. The existence of these thrusts was fundamental in calibrating the model. Firstly, a quasi-impermeable barrier was created, causing isolation of the Embid de Ariza outcrops. These thrusts act as a disconnection of the preferential flow of water from the outcrops of La Pica and Cardejón towards Alhama, thus leading to a diversion of the flow from the north towards the centre of the basin, which would justify the heating of the water. Even so, due to the high hydraulic pressure, a certain amount of water seeps through the thrusts, which justifies the existence of the 20 L/s of the Embid de Ariza springs. Therefore, this boundary condition was modelled as a low-permeability zone, which was calibrated.

The boundary conditions are shown in Figure 10 below.

A schematic of the workflow diagram concerning the calibration method is shown in Figure 11. It summarises the main parameters used for the validation of the conceptual model from the numerical model.

4. Results

4.1. Simulation and Calibration of Groundwater Flow in the Thermal Aquifer

Based on the model information of Section 3.2, a simulation and calibration process was performed in order to confirm the conceptual model. Calibration is the process of modifying the input parameters in the model until the solutions fit and match the observed data.

The objective is to solve the so-called inverse problem, i.e., trying to estimate and adjust the system parameters (such as conductivity, specific storage, and recharge) in such a way that the model is able to reflect the real behaviour of the groundwater recorded at the site from direct measurements of both spring flows and observation points (piezometric levels).

Firstly, it is worth highlighting the difficulty of calibrating this mathematical model. The vast distance of the model (2295.5 km²), the geological complexity, the great variety of values of hydraulic conductivities, the high number of boundary conditions, etc., together with the scarce hydrogeological knowledge and scarcity of piezometers, meant significant difficulty when adjusting the mathematical model. In any case, the numerical model was successfully calibrated in accordance with the conceptual model.

Mathematical model calibration was based on the adjustment of the piezometric levels and the flow rates of the successive thermal upwellings of the springs.

The simulation of the numerical model resulted in the groundwater contour shown in Figure 12, where the state of the piezometric level of the Cretaceous limestone aquifer is observed.

The calibration elements used to adjust the mathematical model are shown in Figure 13 and correspond to those previously indicated: firstly, the piezometric levels corresponding to the observation points of the “Confederación Hidrográfica del Ebro [22]” and secondly, the flows of the springs.

4.2. Piezometric Calibration

Calibration consisted of carrying out a series of simulations in which the model parameters (recharge, permeability, etc.) indicated in the conceptual model were varied in order to adjust the simulation results to the values established at the observation points, reducing the residual error between the simulated and observed values as much as possible.

Observation points were obtained from the piezometers of [22]. These piezometers, as shown in Figure 13, are located in Alhama de Aragón, Embid de Ariza, and Deza. According to the data of the Confederación Hidrográfica del Ebro [22], the levels measured at the observation points correspond to those in Table 1, all measured in units of meters above sea level (m.a.s.l.).

As can be seen, the calibration of the observation points is considered acceptable, as the measured and simulated values are very close (Table 1 and Figure 14).

It should be pointed out that there are very few observation points available, so other indicative elements were used to help with calibration. For example, in the village of Almaluez there are three surge wells, and although no pressure data are available, they served as a semi-quantitative calibration point. In fact, care was taken to ensure that the piezometric levels in the calibrated model are always above the topographic surface. As a counterexample, there is a Cretaceous limestone quarry near the village of Almazúl in the Tertiary, which is a window to the aquifer at 1050 m and through which water never emerges. Thus, this point also served to ensure that the groundwater contour lines of the model could never be above the bottom of the quarry.

4.3. Calibration of Spring Flow Rates

As indicated above, the information about observation points is scarce, and to guarantee a better and complete fit of the model, the flow rates of the different thermal springs were used as a calibration tool. It should be pointed out that the values of the flow rates measured for the springs shown in

Table 2. Comparison of gauged and simulated springs in the calibration process.

Springs	Gauged Flows (L/s)		Simulated Flow Rates (L/s)	Deviation
	Maximum	Minimum
Jaraba	500	600	609.81	2%
Alhama de Aragón	434	520	551.85	6%
Embid de Ariza	20	40	42.5	6%
Deza (South)	120	140	102.13	0%
Deza (North)			29.28
San Roquillo	5	10	10	0%
Almazúl	5	10	6.37	0%

have a certain degree of uncertainty since, although they are extraordinarily constant over time, it is not easy to measure them in the field. These are, however, the ones that appear officially and in public studies [23].

As with the observation points, the objective was to adjust the simulation results to the real flow values observed in the springs. In this way, the calibration of the model was completed by adjusting both the simulated flow and piezometric values, reducing as much as possible the residual error between the simulated and observed values. The springs with which the model was calibrated are shown in Table 2 and Figure 15.

As can be seen, the calibration of the flow rates is more than acceptable, presenting a really good fit (Table 2 and Figure 15). Therefore, according to the piezometric and spring flow calibration, the mathematical model can be considered adjusted.

4.4. Transit Time Estimation

Owing to mathematical modelling of the hydrothermal aquifer, it was possible to approximate the effective porosity of the hydrothermal system as well as the Darcy velocity.

In order to carry out this simulation, the Modpath package was used to simulate the path of a contaminant plume along the aquifer. This allowed us to estimate the time it takes for a pollutant to travel through the aquifer until it emerges at the discharge points.

To achieve this, a contaminant located in the northernmost area of the calcareous outcrops of the Aragonese branch was introduced, and another contaminant was introduced in the outcrops of this same formation corresponding to the catchment area of the Deza springs.

In Figure 16, we can see the progress of the tracer after 20 years and 80 years. It can be seen in the figure on the left how, after around 20–25 years, it comes out through the Deza springs. In the figure on the right, it comes out after 80 years. It leaves through the springs of Alhama de Aragón, then through Jaraba, having made a similar route concentrated in the deepest and most conductive area, although the peripheral flow that flows into Jaraba makes it slightly longer and older. Although these have a semi-quantitative and approximate value, it is observed that the ages obtained in the numerical model and those obtained by tritium are of the same order of magnitude. Indeed, it is observed that the simulation shows that the waters of the Deza springs are more recent than those of Alhama and Jaraba, and that the age obtained from the model coincides with the time estimated by tritium from the latest analyses, which is about 20–25 years. In the case of Alhama and Jaraba, we observe that the ages obtained by tritium from the most recent analyses give us an age of more than 61 years, so 80 years could fit.

According to the Modpath simulation, it was observed that with a porosity of 0.3%, the water takes approximately one century to reach the thermal springs of Alhama de Aragón and Jaraba. In the case of Deza, a duration of 20 years has been estimated, which coincides perfectly with the time estimated by tritium, as can be seen in [9].

4.5. Discussion

This thermal system is characterised by the hydraulic continuity of a calcareous formation located at great depth, as occurs in other examples of large karst aquifers; see [24,25,26].

The circulation of water in this thermal karst system is driven by gravity and caused by topographic gradients [27]. The depth of the aquifer influences the circulation depth and, therefore, the resulting water temperature, according to normal geothermal gradients, without additional heat input. Temperature-induced density gradients and reduced viscosities act simultaneously and facilitate and accelerate the flow further along the aquifer as well as the upwards flow of hot water towards the springs. Indeed, the distribution of the heat source and its depth have influenced the variation in hydraulic conductivity in the aquifer.

The simulation carried out with the calibrated porosity values of 0.3% provided a flow model with average velocities of about 3 m/day (1 km/year), which are entirely consistent with the age of the different thermal springs determined with tritium isotopes [8]. This study draws attention to the importance of topography and other different driving forces of groundwater flow in deep carbonate sequences located in adjacent river basins. The subsurface flow in the Tertiary that adjusts to the dividing lines between rivers does not exist in the deep thermal aquifer, where a huge transfer of groundwater between neighbouring river basins is observed. Thus, it is a thermal system in a natural regime associated with a calcareous aquifer that occupies the dividing line of two neighbouring basins (Duero–Ebro), allowing underground hydraulic communication between them, which can be explained by the flow of water directed towards the springs located in the lower-altitude basin (Ebro). It should be noted that we are located on the western edge of the Castilian Duero plateau, with altitudes of around 1000–1100 m, surrounded by altitudes of around 600 m belonging to the Ebro basin. The slopes are moderate and justify the fact that it is a low-temperature thermal system. In addition, the fact that it is the largest plateau in Europe explains the large volume of water transferred, despite having modest recharge and precipitation (500–600 mm).

To calibrate the model, it was necessary to assign recharge values both for the poorly permeable Tertiary (k = 2 × 10⁻⁶ m/s) and for the Cretaceous thermal aquifer. In the Tertiary, the resulting recharge is somewhat greater in the flat areas of the tributaries of the Duero Basin (16.5 mm/year) than in those of the Ebro, which are much steeper (2.75 mm/year). However, one of the most important aspects of the calibration process was the need to review the behaviour of the Cretaceous aquifer in the areas where it emerges in the Rituerto basin (Duero), where it was not known where the water was going and the hydraulic balances were unknown. It has been proven that the headwaters of the Rituerto are a large recharge polje where the flow is directed below the Tertiary towards the Alhama and Jaraba springs. Both here and in the La Pica mountain range, the recharge values are the maximum (228 mm/year) compared with other areas feeding the aquifer (Sierra de Cardejón and Miñana: 171 mm/year).

On the other hand, as is known, faults can form conduits or barriers, but their hydraulic function depends on several factors and is often difficult to predict [28,29]. Hydraulically conducting faults are crucial for the development of other examples of thermal systems, especially at small scales [30,31,32,33]. In this system, it has been shown that, at the scale of such a large basin, faults have a minor influence on the flow pattern and the temperature field. Although existing faults were considered in the modelling, they do not explain by themselves the large volume of thermal water evacuated. The influence of the Cretaceous limestone series has been recognised as greater, since near the borders with the Tertiary and near the discharge springs of the system, it constitutes a set of vertical and quite karstified limestone strata of about 350 m thick. In the opposite direction, faults have also been located whose jumps (faults) act as barriers or hydraulic disconnections. The careful review of the lack of hydraulic continuity in certain sectors of the aquifer has shown that it is due to fault jumps that break the aforementioned continuity.

The successful calibration of the mathematical model from all the existing observation piezometers, the flow rates of the simulated springs with respect to the real values, as well as the calibration of the velocities that coincide with those measured with tritium, have allowed us to confirm and justify the conceptual model proposed in [9].

4.6. Conclusions

Owing to the integration of diverse hydrogeological information in a previous study [9], it has been possible to know the origin of the most voluminous thermal springs in Spain. The modelling of the regional flow of the studied thermal system has proven to be an integrative and fundamental tool that has mainly served to confirm and quantify various aspects of the initial hydrogeological conceptual model.

This thermal system is associated with a large sedimentary basin and is clearly related on a large scale to a deep regional flow system characterised by the hydraulic continuity of a hydrostratigraphic calcareous formation of Cretaceous age located at great depth.

The regional flow modelling of the thermal system studied has served to confirm and quantify various aspects of the initial hydrogeological conceptual model and to determine the origin of the most voluminous thermal springs in Spain. Among other aspects, the calibration of the model has confirmed the low permeability of the Tertiary and has also confirmed the recharge zone, which had been somewhat disputed, by quantifying it. But it has also required a careful review of the lack of hydraulic continuity in certain sectors of the aquifer due to fault jumps and the low importance of the existing faults as transport routes to explain, by themselves, the large volume of water evacuated.