Assessing 1D Hydrodynamic Modeling of Júcar River Behavior in Mancha Oriental Aquifer Domain (SE Spain)

Abstract

In times of population growth, climate change, and increasing water scarcity around the world, it is important to take an objective look at water, a fundamental resource for life. Hydrodynamic modeling makes possible the research of different aspects of the water cycle and the evaluation of different hydrological and hydrogeological forecasting scenarios in the short and medium terms. The present research offers a more detailed scope at the hydrodynamic processes and their space-time distributions on a UE pilot in the Júcar River Basin, providing a calibrated and validated hydrodynamic model of 121 km river reach for 45 years period (1974–2019) on a daily scale. The obtained information is about discharge and water depths along the Júcar River reach within the hydrogeological boundaries of the Mancha Oriental Aquifer (MOA). The river–aquifer interactions have been represented as dynamic boundary conditions expressed as a difference between observed discharges measured in 3 gauging stations. The obtained calibration error performance evaluations of observed and simulated values cover two periods, according to observed data availability from gauging station 08036 with resulting R² for both discharges and water depths over 0.96. The model validation results were obtained for a different gauge 08132 and the determination coefficients R² also perform very well with value of 0.90. The model developed might be useful for decision making in water resources management and can be used to generate simulated time series of water depths, levels, discharges, and velocities in reaches where gauging measurements are not available with a desired space-time resolution (from meter/second to kilometer/month). Estimation of critical discharge value (1.973 m³s⁻¹) for system equilibrium, based on the balance between losing and gaining sub-reaches of the river, is also made with a statistical significance at 95% for hydrologic years 2007–2010, period influenced by restrictions in groundwater withdrawals. The results of the present research are important for the proper and objective management of the scarce water resources on a watershed scale in Júcar River Basin, a complex case study representing semiarid climate, growing anthropogenic pressures, and complex river–aquifer interactions. The used approach of dynamic representation of the river–aquifer interactions as distributed source boundary condition in the one-dimensional hydrodynamic model might be applied in another study case on similar scale.

1. Introduction

Continental water is a natural resource necessary to satisfy vital human needs. In the coming decades, humanity will face significant challenges, not only in meeting these needs, but also in encountering a balance between different uses and their impacts on the associated natural ecosystems [1,2,3,4]. In view of the estimated population growth that will increase the world’s population from the current 8 billion to 9.2 billion by 2050, according to projections by the United Nations [5], there are reasons for serious concerns in several regions of the world about the availability of water and food resources to ensure the development of future generations. According to [6], water degradation is now a widespread phenomenon in many river basins around the world, largely due to increased human activities. In the case of river basins with arid and semi-arid climates, the problems are even more acute. In these regions, the principal water use is for agricultural irrigation mainly supplied by groundwater [7], which can generate serious problems in surface water bodies [8].

It is well known that as part of the hydrological cycle, surface water and groundwater are hydraulically connected as a unitary resource. However, in many areas of the world they have been managed as two separate entities (for example, in California [9,10]). The joint use of surface water and groundwater is not a new concept; since the last quarter of the 20th century, it has been analyzed by different authors [11,12,13] and is successfully applied in different parts of the world [13,14,15]. In this regard, it should be noted that organizations such as the World Bank or the works of [16,17] propose joint use as a form of water management to fight poverty through agricultural development. Effective (real) joint use of surface and groundwater resources requires understanding and comprehension of the complex relationships between the two phases of the resource.

In fact, water management faces significant challenges in cases where available water resources are at the limit of safe yield, and where surface and groundwater are highly interconnected. Moreover, complex local dynamical processes exist such as pore bottom clogging, which causes changes in the hydraulic connection between the surface and the groundwater bodies [18,19,20,21,22]. In these cases, short- and medium-term hydrological and hydrogeological forecasting scenarios are necessary for proper planning and management of water resources.

A paradigm of the problems described above can be observed in the semi-arid basin of the Júcar River (SE Spain) (Figure 1) as it crosses the Mancha Oriental aquifer domain (one of the most extensive carbonate aquifers (7520 km²) in southwestern Europe). In this area, agriculture is currently (data from 2018) the largest water demands user (94.5%), followed by residential (3.8%) and industrial (1.6%) water supply. This demand is almost entirely supplied by groundwater which, due to the negative balance (for the river), has caused a decrease in groundwater levels that has caused Júcar River to experience a change from a gaining to a losing river in certain reaches [23,24]. The decrease in downstream discharges of Júcar River and the growing demand in the entire system are generating political disputes between stakeholders in the two main regions where Júcar River flows: Castilla-La Mancha (upper and middle sub-basins) and Valencia (low sub-basin) (Figure 1).

The system is managed by the Júcar River Basin Authority, in Spanish Confederación Hidrográfica del Júcar (acronym CHJ). As part of the management actions and based on the recommendations of the Water Framework Directive [25], a series of measures have been implemented, including policies to control water use, annual exploitation plans, different reservoir-release scenarios and public water rights acquisition offers (in Spanish Oferta Pública de Adquisición de Derechos, acronym OPAD) during droughts. Proposals such as increase of wastewater reuse, artificial aquifer recharge, reducing water use in agriculture by improving irrigation efficiency, new water allocation mechanisms, etc., are also considered. In addition, a groundwater flow model of the Mancha Oriental Aquifer has been developed and is available for more than 10 years [7,23,26], which is periodically updated and serves as a management tool for decision making regarding the evolution of groundwater levels and the determination of water balances in the context of pumping replacement scenarios, climate change, etc.

Nevertheless, groundwater flow model is a regional model with a cell size of 1 km × 1 km and a 1-month time step, while Júcar’s riverbed as it crosses through the MOA domain does not exceed few tens of meters in width and few meters in depth, where water may be released from the Alarcón reservoir on a daily scale. For this reason, and although a good approximation is obtained of the river–aquifer interactions with this model on a monthly time scale, it is impossible to reproduce the discharge peaks generated by the water releases’ policies from the upstream reservoir. Furthermore, the clogging of the bottom pores could lead to more uncertainties in the estimation of the high water level peaks, because of the accumulation of suspended solids, as well as the detritus on the river bottom [19].

Because of the above mentioned reasons, a hydrological forecasting tool is needed to quickly and unambiguously answer questions such as: what is the spatio-temporal distribution of discharge and water level along the river, what is the hydrological balance per river reach and period, what are the most likely reservoir water releases that will affect the river–aquifer relations on a daily time scale, or what is the minimum water discharge to be released from the reservoirs upstream? In addition, such a tool can provide completely simulated time series with a high degree of statistical reliability, both of discharge (Q) and water depths (h) for the studied period that can serve as input to models that need spatio-temporal information adjusted to the needs of more specific studies—hydrological, hydrogeological, ecological, economic, and/or environmental management.

In this context, the present study aims to create a numerical model and apply a hydrological simulation to investigate the space-time dynamics of the Júcar River stream flow, applying the one-dimensional hydrodynamic module (HD module) of the MIKE 11 software (hereinafter 1D-HD model). The model, based on daily data from the gauging stations, was calibrated and validated for a period of 45 years, covering important events from hydrological and water management perspective, caused by the incremented use of surface and groundwater. In that sense, this work represents a significant advance in the effort to create a model that is more plausible to the hydrological processes occurring in the Júcar River basin, in order to improve planning and the successful implementation in practice of future joint use management scenarios. The results of this study will, firstly, help decision makers to understand the current situation and, secondly, facilitate the design of future management options to improve the sustainability of water resources.

2. Study Area

The Mancha Oriental Aquifer system (already introduced as MOA) has a temperate Mediterranean climate that is highly influenced by a continental component because it belongs to the high plateau of the Iberian Peninsula with an average elevation above the sea level of approximately 700 m. It is usually characterized by very dry summers, with rainfall in spring in the western sector and in autumn in the eastern sector, being scarce in winter. The continental climate is evidenced with hot summers (average monthly temperatures of 22 °C) and cold winters (average monthly temperatures of 6 °C) with important daily oscillations. The average precipitation is 350 mm year⁻¹, varying from 280 mm year⁻¹ in the southern zone and 550 mm year⁻¹ in the northern zone. During dry years, the average precipitation is approximately 150 mm year⁻¹ (very common), while in the wet years 750 mm year⁻¹ can be reached. Potential evapotranspiration values exceed 1200 mm year⁻¹, so the area is considered semi-arid to arid. Geomorphologically, the area is characterized by large depressions of intramiocene age filled with later materials (Tertiary–Plio–Quaternary) that conserve their horizontal disposition, causing the practically flat relief of the area. The high plateau is interrupted only by the valley excavated by Júcar River, the main fluvial course that crosses the MOA domain. The hydrological network is scarcely developed; therefore, the surface runoff (except for the Júcar River) is practically zero or insignificant, ranging between 0 and 5 mm year⁻¹. The Alarcón reservoir, a basic element in the regulation of the Júcar River, is located upstream on its entrance into MOA domain (Figure 1).

The MOA is a groundwater body formed by the superposition of three carbonate aquifer hydrogeological units of Jurassic, Upper Cretaceous, and Miocene age, separated by impermeable and semipermeable layers that constitute the Lower Cretaceous and Tertiary detrital materials [7,27]. The primary groundwater discharge zone of MOA is Júcar River. The aquifer system in natural regime yielded its accumulated water resources to the river at a rate of 320 hm³. However, over the last 45 years, approximately 1000 km² of traditionally rain-fed land has been transformed into irrigated surface through the widespread use of groundwater for irrigation, which has led to significant socioeconomic development in the region. In recent decades, MOA has been in a quasi-equilibrium state with a stabilized annual balance. The annual volume extracted for irrigation was 320 hm³, while the estimated renewable resources were 389 hm³ [28]. This stabilization is currently threatened due to the proliferation of wells with irrigation rights of up to 7000 m³year⁻¹ (treated differently in the concession process, without being subject to the same controls as the wells with bigger extraction volumes), used to irrigate crops in expansion such as almond and pistachio trees. This fact together with the volume of environmental restrictions set at 114 hm³ by the last Júcar Basin water management plan [29] puts MOA in a state of quantitative overexploitation with the consequent negative effects on the interactions with Júcar River.

3. Materials and Methods

3.1. Data Acquisition

The hydrological data used as input for the 1D-HD model of the middle reach of Júcar River within MOA domain, developed with MIKE 11, include information about the topography, as well as the measurements obtained at the different gauging stations existing in the study area.

The definition of the fluvial network was obtained in vector format at scale 1:25,000 [29], the course of the river being defined through interconnected linear sections of variable length denominated “chainage”. Along the studied river reach, 24 cross topographic sections were surveyed at specific locations, 19 were used in the present study. For this purpose, based on spatial information on the location of each section and the geodetic vertices of the National Geodetic Network [30], the geometry of the Júcar River cross sections was obtained with differential GPS total station model Trimble S6 3” and the corresponding accessories with a distance accuracy of 3 mm ± 2 ppm. Based on the data obtained in the field with the use of Trimble Geomatic Office software, the information was further processed for its adequacy and subsequent use in other computer programs. The geometry of the gauging stations was obtained from the reports available in the yearbooks of the Ministerio para la Transición Ecológica y el Reto Demográfico [31]. The specific geographic location of each of the cross sections used in the study is shown in Figure 1.

Daily discharge and water depths data were obtained from the gauging stations available in the study area [31] and their exact geographic locations are represented on Figure 1. With this information, time series files compatible with the MIKE11 hydrodynamic modeling program were created. The study period was conditioned by data availability and covers the interval from 1 January 1974 to 30 September 2019. General description and descriptive statistics of the daily observed time series of discharge mean values (Q, (m³s⁻¹)) and water depths (h, (m)) obtained from the four gauging stations used for modeling purposes can be found in Table 1. Furthermore, as complementary analysis to the descriptive statistics for the observed time series of discharge and depths, Hurst coefficient estimates [32,33,34,35,36,37] were used, since they present non-stationary component due to the increasing underground water extractions starting in mid-1970s, causing changes in the river-aquifer interactions. The Hurst coefficient H is a measure of self-similarity or a measure of the duration of the long-term dependence of a given process. Time series characterized by long-term dependence exhibit a slow decrease of the autocorrelation function. For a stationary process with long-term dependence H $\in \left(0.5, 1\right]$, while H > 1 indicates a non-stationary unbounded process [33]. More comprehensive non-stationary analysis of the Jucar River–MOA interactions, their temporal evolution, and relation to climatic variables and drivers can be found in [24].

Table 1. General description and descriptive statistics of the observed daily time series with mean values for river discharge (Q, (m³s⁻¹)) and depths (h, (m)) obtained from four gauging stations with available information for the period 01.01.1974–30.09.2019 and used for the purposes of the modeling. The location of the gauging stations can be seen in Figure 1 and their cross sections in Figure 2. The Hurst coefficients H were estimated using regression on periodogram (RP), averaged wavelet coefficients (AWC) and the detrended fluctuation analysis (DFA) methods.

General Description
CHJ Code		08129		08132		08036		08144
Location		El Picazo		Puente Carrasco		Los Frailes		Alcalá del Júcar
UTM ETRS89 H30 X		578628		584728		608082		635980
UTM ETRS89 H30 Y		4368600		4341164		43327912		4339792
Elevation, (m.a.s.l.)		694		647		605		514
Studied variables		QPicazo	hPicazo	QCarrasco	hCarrasco	QFrailes	hFrailes	QAlcalá	hAlcalá
		(m3s−1)	(m)	(m3s−1)	(m)	(m3s−1)	(m)	(m3s−1)	(m)
Time periods with data availability *		01.01.1974– 30.09.2019	01.01.1974– 30.09.2019	01.01.1974– 30.09.1986	01.01.1974– 30.09.1986	01.01.1974– 30.09.1988   25.10.1991– 30.09.2019	01.01.1974– 30.09.1988   26.10.1991– 30.09.2019	01.10.1974– 30.09.1980   12.06.1984– 30.09.2019	13.06.1984– 18.04.2011   30.09.2015– 30.09.2019
Data use in the model		Boundary cond. in	Calibration	Validation	Validation	Calibration	Calibration	Boundary cond. out	Boundary cond. out
Descriptive statistics
Mean (m3s−1)		8.055	0.568	10.120	0.907	9.416	0.480	9.915	0.400
Minimum (m3s−1)		0.044	0.080	0.000	0.000	0.344	0.010	0.001	0.130
Maximum (m3s−1)		42.440	3.440	46.600	1.990	57.102	1.220	74.000	1.170
Standard deviation (m3s−1)		8.260	0.263	8.267	0.335	9.165	0.199	8.283	0.215
Coefficient of variation		1.025	0.463	0.770	0.369	0.973	0.4715	0.835	0.508
Variance		68.230	0.069	68.342	0.112	84.005	0.0395	68.606	0.046
H—RP		1.328	1.310	1.216	1.230	1.345	1.360	1.297	1.335
H—AWC		1.180	1.068	1.199	1.204	1.277	1.262	1.218	1.220
H—DFA		1.162	1.152	1.205	1.203	1.196	1.192	1.174	1.187

* Missing data < 25 days are not considered in the table.

Table 1. General description and descriptive statistics of the observed daily time series with mean values for river discharge (Q, (m³s⁻¹)) and depths (h, (m)) obtained from four gauging stations with available information for the period 01.01.1974–30.09.2019 and used for the purposes of the modeling. The location of the gauging stations can be seen in Figure 1 and their cross sections in Figure 2. The Hurst coefficients H were estimated using regression on periodogram (RP), averaged wavelet coefficients (AWC) and the detrended fluctuation analysis (DFA) methods.

General Description
CHJ Code		08129		08132		08036		08144
Location		El Picazo		Puente Carrasco		Los Frailes		Alcalá del Júcar
UTM ETRS89 H30 X		578628		584728		608082		635980
UTM ETRS89 H30 Y		4368600		4341164		43327912		4339792
Elevation, (m.a.s.l.)		694		647		605		514
Studied variables		QPicazo	hPicazo	QCarrasco	hCarrasco	QFrailes	hFrailes	QAlcalá	hAlcalá
		(m3s−1)	(m)	(m3s−1)	(m)	(m3s−1)	(m)	(m3s−1)	(m)
Time periods with data availability *		01.01.1974– 30.09.2019	01.01.1974– 30.09.2019	01.01.1974– 30.09.1986	01.01.1974– 30.09.1986	01.01.1974– 30.09.1988   25.10.1991– 30.09.2019	01.01.1974– 30.09.1988   26.10.1991– 30.09.2019	01.10.1974– 30.09.1980   12.06.1984– 30.09.2019	13.06.1984– 18.04.2011   30.09.2015– 30.09.2019
Data use in the model		Boundary cond. in	Calibration	Validation	Validation	Calibration	Calibration	Boundary cond. out	Boundary cond. out
Descriptive statistics
Mean (m3s−1)		8.055	0.568	10.120	0.907	9.416	0.480	9.915	0.400
Minimum (m3s−1)		0.044	0.080	0.000	0.000	0.344	0.010	0.001	0.130
Maximum (m3s−1)		42.440	3.440	46.600	1.990	57.102	1.220	74.000	1.170
Standard deviation (m3s−1)		8.260	0.263	8.267	0.335	9.165	0.199	8.283	0.215
Coefficient of variation		1.025	0.463	0.770	0.369	0.973	0.4715	0.835	0.508
Variance		68.230	0.069	68.342	0.112	84.005	0.0395	68.606	0.046
H—RP		1.328	1.310	1.216	1.230	1.345	1.360	1.297	1.335
H—AWC		1.180	1.068	1.199	1.204	1.277	1.262	1.218	1.220
H—DFA		1.162	1.152	1.205	1.203	1.196	1.192	1.174	1.187

* Missing data < 25 days are not considered in the table.

Figure 2. The geometries of Júcar River gauging stations are shown as introduced in the Cross Section file of MIKE11 for (a) El Picazo, (b) Puente Carrasco, (c) Los Frailes, and (d) Alcala del Júcar. Minimum and maximum water levels are indicated with green and red dashed lines, respectively, while intermediate water levels are drawn in blue. The rest of the cross sections with irregular geometries were also introduced but are not shown here. Their locations can be found on Figure 1.

Figure 2. The geometries of Júcar River gauging stations are shown as introduced in the Cross Section file of MIKE11 for (a) El Picazo, (b) Puente Carrasco, (c) Los Frailes, and (d) Alcala del Júcar. Minimum and maximum water levels are indicated with green and red dashed lines, respectively, while intermediate water levels are drawn in blue. The rest of the cross sections with irregular geometries were also introduced but are not shown here. Their locations can be found on Figure 1.

3.2. Model Setup

3.2.1. Hydrodynamic Modeling Setting

In order to model the flow of Júcar River, the MIKE 11 hydrodynamic module was used [38] to obtain the discharges (Q) and water depths (h) in the middle catchment of the JRB, in the reach between El Puente Picazo (chainage 0.00 m) and Alcalá del Júcar (chainage 120,990.29 m) gauging station (Figure 1). The software was developed for detailed modeling of river networks including flood plains, road overflows, culverts, and dam control, including dam breach. The numerical model is based on kinematic, diffusive, or fully dynamic, vertically integrated mass conservation and fluid momentum conservation (Saint-Venant equations, see Equations (1) and (2)). The numerical solution of the equations is an implicit finite difference scheme, structured in such a way as to guarantee independence of the wave type (i.e., kinematic, diffusive or dynamic). The solution of the equations is based on the following assumptions: (i) water is incompressible and homogeneous (i.e., differences in density are negligible); (ii) the slope of the bottom is small, so that the cosine of the angle towards the horizontal can be assumed to be 1; (iii) the wavelength is large compared to the water depth, i.e., it can be assumed that the flow is parallel to the bottom (vertical accelerations are neglected and hydrostatic pressure variation in the vertical direction can be assumed); (iv) the flow is subcritical (supercritical flow can be modeled with MIKE HD, but more restrictive conditions apply).

$$\frac{\partial \mathrm{Q}}{\partial \mathrm{x}}+\frac{\partial \mathrm{A}}{\partial \mathrm{t}}=\mathrm{q}$$

(1)

$$\frac{\partial \mathrm{Q}}{\partial \mathrm{t}}+\frac{\partial \left(\frac{{\mathsf{\alpha }\mathrm{Q}}^2}{\mathrm{A}}\right)}{\partial \mathrm{x}}+\mathrm{gA}\frac{\partial \mathrm{h}}{\partial \mathrm{x}}+\frac{{\mathrm{n}}^2\mathrm{gQ}\left|\mathrm{Q}\right|}{{\mathrm{AR}}^{4/3}}=0,$$

(2)

where Q is the discharge (m³s⁻¹); A is the cross section flow area (m²); q is the lateral inflow (m²s⁻¹); h is the water level above a reference datum (m); x is the downstream path (m); t is the time (s); n is the Manning coefficient (s m^−1/3); R is the hydraulic or resistance radius (m); g is the gravity acceleration (m²s⁻¹); and α is the momentum distribution coefficient [38]. The MIKE 11 HD module solves numerically the river networks and floodplains as a system of interconnected branches. Water levels and discharges are calculated at alternating points along the river branches as a function of time.

3.2.2. Conceptual Model, Discretization, and Boundary Conditions

The first phase in building a 1D-HD model is to define the conceptual model, i.e., how it is initially considered that the system works in reality. From a hydrological point of view, the study area can be divided according to the location of the upstream reservoir and the gauging stations along the river stream, and four reaches can be distinguished (Figure 1). Section 1: from the Alarcón reservoir to the El Picazo gauging station, Section 2: from the El Picazo gauging station to the Los Frailes gauging station, Section 3: from the Los Frailes gauging station to the Alcalá del Júcar gauging station, and Section 4: from the Alcalá del Júcar gauging station to the MOA domain exit. The hydrological and hydrogeological functioning of the system is closely related to the geomorphology of the Júcar River as it flows through MOA domain, which is conditioned by the lithology, the structural arrangement of the geological materials and the topographic slopes that affect the river–aquifer relations [39].

Both Reach 1 and Reach 4 were excluded from the model due to the diversity of anthropic elements that affect the flow of the river. Section 1, with a slope of approximately 2.50‰ and 28 km in length, is formed by a narrow valley over Cretaceous materials (limestone and dolomites). Its functioning depends on both the releases from the Alarcón reservoir (including contributions from the Tajo–Segura water transfer channel, a hydraulic infrastructure outside the JRB) and the turbining of the El Pizaco hydroelectric power plant. Section 4, 36 km long and with a slope of 2.00‰, runs from the village of Alcalá del Júcar to the exit of the MOA, where Mesozoic materials (mostly Cretaceous) emerge. It is affected by the reservoirs-drainage of the El Molinar reservoir and the outflows derived to the Cofrentes hydroelectric power plant located outside of the hydrogeological system (Figure 1). The reaches, subjects of the present research, are 2 and 3. Reach 2 is approximately 74 km long and two sectors can be delimited: (a) between El Picazo and the village of Villalgordo del Júcar, where the river flows through Plio–Quaternary materials and the slope is reduces to 1.60‰ and, (b) from the latter to the vicinity of the gauging station of Los Frailes, which is characterized by the presence of Tertiary materials (limestones and marls) and the slope decreases to 1.30‰. Section 3 of approximately 47 km in length continues crossing tertiary materials and extends to the town of Alcalá del Júcar increasing the slope up to 2.00‰, forming a narrow valley, enclosed between deep gorges.

In relation to the interactions between Júcar River and the MOA, in the natural regime the Júcar River stream presented a clearly gaining behavior (it received contributions from the aquifer until 1985). At present, the Júcar River has a gaining–losing behavior: its stream is losing from its entry into the groundwater body domain to the vicinity of the Cuasiermas location, where it becomes a gaining river. However, in the entire Júcar River Basin hydrologic unit (its upper, middle, and low sub-basins are shown on Figure 1), there is a predominance of the diffuse connection conceptual model in the river-aquifer interactions, both direct and indirect, along the course of the river [40].

The construction of HD model is based on reference information about the topography of the river and the floodplain. Furthermore, the 1D dynamic wave equations for continuity and momentum conservation (Equations (1) and (2)), are solved numerically by vertical integration. A computational grid of alternating Q (discharge) and h (water level) points is automatically generated on the basis of the user requirements; respectively Q and h are computed at each time step. Q-points are always placed midway between two h-points, while the distance between h-points may differ. Based on the river network, the cross sections obtained with a topographic study have been allocated along the studied river reach (h points) and the corresponding chainages are introduced in the network file. The minimum distance between two cross sections is 50 m and the maximum is 11,210 m (see Figure 1).

The discretization of river networks for modeling purposes is an important step in order to ensure cost-effective simulations, especially at watershed scale. For 1D model described with Equation (2), the convergence is conditioned by Courant–Friedrichs–Lewy coefficient (Cr), given with the Expression (3):

$$\mathrm{Cr}=\mathrm{u}\frac{\mathsf{\Delta }\mathrm{t}}{\mathsf{\Delta }\mathrm{x}}$$

(3)

where the dimensionless coefficient Cr needs to be less than 1 for stability, u is a characteristic velocity, Δt is the time discretization and Δx is the space discretization. In the case of a shallow water flow where advection is ignored the characteristic velocity is $\sqrt{\mathrm{gh}}$, g and h being as defined in Equation (2). The maximum acceptable time step is determined following the simplification proposed by [41] and commonly used to estimate the step at the moment t + Δt for values of Cr usually assumed in the interval 0.2–0.7:

$$\mathsf{\Delta }\mathrm{t}=\mathrm{Cr}\frac{\mathsf{\Delta }\mathrm{x}}{\sqrt{\mathrm{gh}}}$$

(4)

In the case of the present study, considering the model aims and the information available Δx = 1000 (m) and Δt = 15 (s), sufficient to keep the Courant number less than 1 for the whole simulation. The initial condition was set at a water depth of 0.15 m, which is important for the simulation to run smoothly and avoid unrealistic results for the stream-flow velocities along the river.

The correct definition of boundary conditions is one of the most important steps in the construction of a numerical river model. In the case of the present model, these conditions must include the configuration for upstream, downstream, and distributed sources (named leakage). In MIKE 11, the boundary conditions that can be defined include water depth (h), discharge (Q), differential discharge between reaches (ΔQ), and rating curve (Q/h ratio). The water depth boundary condition can be applied upstream or downstream of the river. The water discharge boundary condition Q can be used at both upstream or downstream and as a lateral flow. The differential discharge between reaches boundary condition require an upstream and downstream chainage between which the inflow is equally distributed. A boundary condition set as a rating curve (Q/h ratio) can only be applied downstream of the system.

Regarding the boundary condition upstream of the river average daily (24 h) discharge Q (m³s⁻¹) data from the El Picazo gauging station were used for the studied time period, which were entered as a time-varying (dynamic) open boundary condition upstream of the river. Observed discharge at the input varied from 0.04 m³s⁻¹ to 42.44 m³s⁻¹ during the study period, and the average discharge was 8.04 m³s⁻¹ (Table 1). Downstream of the river; (gauging station Alcalá del Júcar), the rating curve of the cross section is entered in table form as a boundary condition.

About a distributed source boundary condition (leakage), these are used to represent the river–aquifer interactions between the surface water body (the Júcar River) and the groundwater body (Mancha Oriental Aquifer). That boundary condition was dynamically represented by obtaining the differences of the observed discharges in the three gauging stations, namely Leakage₁ and Leakage₂, expressed as follows:

Leakage₁(t_k) = Q_LosFrailes (t_k) − 1/2 ((Q_Picazo(t_k) + Q_Picazo(t_k+1))

(5)

Leakage₂(t_k) = Q_Alcalá (t_k) − 1/2 ((Q_LosFrailes (t_k) + Q_LosFrailes (t_k+1))

(6)

where Q is the discharge, Leakage is the difference between the observed discharges in two consecutive gauging stations, where t_k and (t_k + 1) are the kth and the ((k + 1)th) times of the observation. The obtained time series were introduced as distributed source boundary condition at the corresponding river reaches—Leakage1 is between El Picazo and Los Frailes gauging stations and Leakage2 is between Los Frailes and Alcalá del Júcar gauging stations.

3.2.3. Calibration and Validation Model Evaluation Criteria

The calibration and the validation processes consisted in adjusting the model boundary conditions and parameters (within reasonable ranges) so that both the Q and h were close to the gauging stations observed values. As an indicator of the accuracy of the calibration-validation process, the error performance is being evaluated followed statistics’ recommendations for hydrodynamic modeling given in [42] and references therein. These goals were achieved by a trial-and-error method, represented schematically on Figure 3.

Manning coefficient n is a highly sensitive parameter used in the Saint-Venant formulation that describes the flow below a pressure surface in a fluid (Equation (2)), therefore it is common to use it for model calibration [43], and in most cases, it is verified versus observed discharges and/or water depths.

The Manning coefficient plays a crucial role in natural river’s modeling and represents the roughness of the riverbed or the friction between the streamflow and different types of surfaces. In the case of one-dimensional hydrodynamic numerical modeling Manning coefficient n is considered a key parameter. The works of [44,45] are fundamental for establishing the parameter for different materials, flow conditions and geographic locations (mountains, plains, artificial channels, etc.). Its magnitude depends on particle size, vegetation cover, channel alignment, channel roughness, meanders, and other river characteristics.

The use of global and spatially varying values for n in the 1D hydrodynamic modeling could be particularly useful in long river reaches when the number of gauging stations is limited [46], which is the case of the present research. Global values are suitable when comparing similar cross sections of different rivers and with the assumption of low to mild water level changes.

In the case of water depths with high temporal variability, it is important to integrate a dynamic component of the value n, in order to improve the model performance. The dynamic component could be water depth dependent and could be implemented by using a different value for n along of the river cross section, based on the assumption that the bed material is heterogeneous, which could lead to minor changes in the velocity field and occurrence of secondary currents [47]. This approach is preferable when modeling long river systems and could provide a solid base for further calibration and verification of the numerical model.

Given the importance of correct model performance for the intended model use [48,49,50] and following the recommendations for evaluation criteria given in [42] and references therein, various statistical methods for error estimation were used. In this study, model performance is assessed by statistical measures and characterization of the studied variable’s error including the coefficient of correlation (R), coefficient of determination (R²), Nash–Sutcliffe Efficiency (NSE) [51] and the root mean square error (RMSE). Their equations, ranges, and optimal values are given in

Table 2. Model performance statistics, their equations, ranges, and optimal values, modified from [42]. O and S are observed and simulated values for a time series with length n and $\overline{\mathrm{O}}$ and $\overline{\mathrm{S}}$ are the mean values, respectively.

Statistic	Symbol	Equation	Range	Optimal Value	(Eq)
Coefficient of correlation	R	$\frac{{\sum }_{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{O}}_{\mathrm{k}}- \overline{\mathrm{O}})({\mathrm{S}}_{\mathrm{k}}- \overline{\mathrm{S}})}{\sqrt{{\sum }_{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{O}}_{\mathrm{k}}- \overline{\mathrm{O}})}\sqrt{{\sum }_{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{S}}_{\mathrm{k}}- \overline{\mathrm{S}})}}$	−1 to 1	−1(negative slope) 1(positive slope)	(7)
Coefficient of determination	R2	${\left(\frac{{\sum }_{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{O}}_{\mathrm{k}}- \overline{\mathrm{O}})({\mathrm{S}}_{\mathrm{k}}- \overline{\mathrm{S}})}{\sqrt{{\sum }_{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{O}}_{\mathrm{k}}- \overline{\mathrm{O}})}\sqrt{{\sum }_{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{S}}_{\mathrm{k}}- \overline{\mathrm{S}})}}\right)}^2$	0 to 1	1	(8)
Nash–Sutcliffe Efficiency	NSE	$1-\frac{{\sum }_{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{O}}_{\mathrm{k}}-{\mathrm{S}}_{\mathrm{k}})}{{\sum }_{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{O}}_{\mathrm{k}}- \overline{\mathrm{O}})}$	−∞ to1	1	(9)
Root mean square error	RMSE	$\sqrt{\frac{1}{\mathrm{n}}\sum _{\mathrm{k}=1}^{\mathrm{n}}({\mathrm{O}}_{\mathrm{k}}-{\mathrm{S}}_{\mathrm{k}})}$	0 to ∞	0	(10)

.

The correlation coefficient (R) and the coefficient of determination (R²) describe the degree of collinearity between simulated (S) and observed data (O). They are widely used in hydrological modeling; however, they are not sensitive to additive and proportional differences between model simulations and observed data [52]. The Nash–Sutcliffe Efficiency, NSE, accounts for the model’s ability to predict variables different from the mean and gives the magnitude of residual variation in comparison to observed data variance. The RMSE is suitable for continuous long-term simulations and is used in model performance evaluation as a measure of the difference between simulated and observed values.

Model performance criteria defined in [42] recommend the use of various statistical measures, among them: R², NSE, and RMSE. The performance criteria on daily, monthly, and annual temporal scales for watershed scale models can be considered “satisfactory” for flow simulations if R² > 0.60 and NSE > 0.50, “good” for 0.75 < R² ≤ 0.85 and 0.70 < NSE ≤ 0.80, and “very good” for R² > 0.85 and NSE > 0.80 [53]. In regard to the water table depth (on a daily scale), the reported statistical evaluation criteria in [54] rate the model “acceptable” for NSE >0.40, “very good” for 0.60 < NSE ≤ 0.75, and “excellent” for NSE > 0.75.

4. Results and Discussion

4.1. Calibration and Validation Results of HD Modeling

During the calibration and validation procedure, the model was initially simulated using the default value of Manning’s roughness coefficient (n = 0.033 s m^−1/3). Observed and simulated Q y h at the Los Frailes gauging station (see Figure 1 and Figure 2) had a fit with an R² considerably lower than 0.80, therefore as suggested by the performance evaluation criteria reported for hydrodynamic models [42], a more precise calibration was needed, in this case, of the value of n along the river. Consequently, n was adjusted to a variable value that would produce the best match between observed and simulated values. The averaged values of n along the river vary from 0.027 (s m^−1/3) to 0.042 (s m^−1/3) and are introduced for each of the sections defined in the model. The characteristics of the riverbed, such as its geology, geomorphology, and surrounding vegetation, were taken into account [44]. The resulting values of n along the studied sub-reaches of the Júcar River, are shown in Figure 4a.

Specific focus was put on the calibration of the n values for Los Frailes gauge cross section where the values of Q and h are used for model calibration. When operating with a constant roughness coefficient as a function of depth, in the sections where observed and simulated values are to be compared, the model results will only be adjusted for a constant Q/h. However, the variation of Q/h over time, mostly due to reservoir-release policies, makes it necessary to assume a depth varying Manning’s coefficient n. Based on the type of construction of this gauging station (geometry, materials, etc.), n was adjusted with respect to the depth following the criteria given in [44,47] (Figure 4b). By varying the roughness as a function of depth, the error is within 8–10% covering the entire range of observed water depths, or ±6 cm over the full range of depths (up to 3.44 m at the model input).

Regarding the results of the calibration at Los Frailes gauging station, simulated and observed discharges and water depths are shown in Figure 5. As it can be seen in Figure 5a, a good fit between the observed and simulated values is achieved with the 1D-HD model. However, with respect to water depths the simulated water depths are higher than those observed in the first period (Figure 5b). In that sense, most probably this is due to the modification of the existing geometry of the gauging station (Los Frailes) from year 1911 and its adoption to the new national automated hydrological information system (acronym SAIH), which was pioneered in the Júcar River basin; built and put into operation in the period 1985–1991. In fact, it is also reflected in the statistic NSE which measures the model’s ability to predict variables different from the mean, for this first period NSE(Q) = 0.99 for the discharge and NSE(h) = 0.45 for the water depths (

Table 3. Calibration and validation statistical performances at Los Frailes and Puente de Carrasco stations (simulated versus observed time series), represented by correlation coefficient (R), determination coefficient (R²), Nash–Sutcliffe efficiency (NSE), and root mean square error (RMSE) for discharge Q (m³s⁻¹) and depths h (m) time series, calculated within periods with availability for observed data.

Error Statistics	R	R2	NSE	RMSE
Calibration—08036 Los Frailes
Q (01.01.1974–30.09.1988)	0.997	0.994	0.994	0.798
h (01.01.1974–30.09.1988)	0.983	0.967	0.446	0.181
Q (25.10.1991–30.09.2019)	0.998	0.996	0.997	0.375
h (25.10.1991–30.09.2019)	0.998	0.995	0.812	0.064
Validation—08132 Puente Carrasco
Q (01.01.1974–30.09.1986)	0.946	0.896	0.723	4.345
h (01.01.1974–30.09.1986)	0.946	0.895	0.791	0.153

). It should be noted that the model maintains the same geometry for the Los Frailes station for the entire simulation, being the one that corresponds to the current one.

In contrast, for the second period, NSE(h) = 0.81 (see Table 3), meaning that the model performance has shifted from satisfactory to very good for the water depths, while for the discharge the values are close to the optimal value 1. For the second period, NSE for Q and h can be rated as very good. The rest of statistics with optimal value 1 (R and R²) for both periods are also performing very well. In practically more than 96% of the cases, the simulated discharges and water depths represent reality. The optimal value for RMSE is zero; therefore, the statistic indicates tolerable agreement between observed and simulated values for Q (0.30 m³s⁻¹) and h (0.06 m).

After the model calibration, the MIKE 11 HD River model was validated using the daily discharge and water depths from Puente de Carrasco (chainage 36,910 m) over the period with data availability (01.01.1974–30.09.1986). The simulated and observed discharges and water depths are shown on Figure 6. Apparently, water depths h fit better than discharges. On Figure 6a can be observed that simulated discharges are higher than those observed for the entire periods. In fact, a constant shift of approximately 4 m³s⁻¹ exists. The error performance statistics (see Table 3) also reflect this situation, the RMSE for the discharge departs from the optimal value zero, however the rest of statistics with optimal value 1 can be rated as very good (R and R²) and good (NSE). Considering that in Figure 6b the water depths adjustment is satisfactory, the reason for this shift could be the way MIKE HD module represents distributed sources (denominated leakage in our study) boundary condition which represents the interactions between the river and any type of inflow-outflow (e.g., aquifer–river interactions). Indeed, distributed sources can only be established between points where complete historical series are available. Therefore, although Leakeage1 is well defined in the El Picazo–Los Frailes reach, it should be different for the El Picazo–Puente de Carrasco and Puente de Carrasco–Los Frailes reaches. A highly probable explanation is that the model overestimates the leakage in the reach El Picazo–Puente de Carrasco and in the following reach Puente de Carrasco–Los Frailes it is underestimated; as a result, they compensate each other. This issue is interesting and serves to advertise to the water authority about the need for flow control at intermediate points.

In view of the results of the calibration and validation of the Júcar River 1D-HD model, there is overall agreement between observed and simulated daily discharges and water depths for the entire studied period. The error performance statistics between simulated and observed data indicate that the 1D-HD model is capable to represent complex and unsteady river flows at a watershed scale. Time series can be obtained for a 45-year period with a required space-time discretization (from seconds/meters to months/kilometers) for discharge, water level, depths, velocities, and shear stress. Consequently, it is possible to determine both in space and time the spatio-temporal distribution of different variable (discharge, depth, water level, velocities, etc.) from the El Picazo to the Alcalá del Júcar gauging station, approximately 121 km in length (Figure 1). As an example, in Figure 7 are shown mean monthly values of discharge Q and water depths h for hydrologic year 10.2007–09.2008. They are represented with a space distribution of 1000 m along the studied river reach. The represented information was used to obtain the values of the river gains/losses in low waters (<5 m³s⁻¹).

4.2. Discussion (Implications for the Water Resources Management)

Despite the ubiquitous availability of data from digital elevation models and the development of 2D/3D hydrodynamic models, 1D modeling remains useful because of its accuracy in describing the hydraulic behavior of river networks from scarce field data, with the guarantee of low computational times [55]. In view of the current research goals, aiming to calibrate and validate hydrodynamic model and obtain discharge and water depths at un-gauged sections of the river for 45 years period, rather than flood inundation analysis, the 1D-HD modeling approach has been chosen. The main reasons are (i) the relatively long river reach (121 km) with scarce spatial information obtained with topographic study and (ii) the high computational time needed for such a long period of time (45 years with a time discretization for the differential equations integration of 15 s).

In fact, there are other river basin systems worldwide, using 1D-HD river models (including coupled models) developed with MIKE DHI software for different purposes: data sensitivity analysis [56], water quality in large river delta [57], and superficial streams–lake interactions on a catchment scale [58]. In all these works, as well as in our own study have been demonstrated the reliability, simplicity and effectiveness of 1D-HD models for regional-scale water resource evaluation, planning and management.

The obtained completely simulated time series of Q and h for the chosen time period along the studied Júcar River reach with a very high degree of statistical reliability can serve as inputs for other models (hydrogeological for instance) and can also be used for the analysis of their fractal and spectral properties through the estimation of their Hurst coefficient. This kind of analysis serves as a measure of the so-called long-range dependence of the time series. In this sense, the complete series studied reflect changes in the Q/h behavior of the river and, therefore, in the river–aquifer relationship. The Hurst coefficient (H) values are shown on Figure 8. They are higher than 1.0 and approaching 1.5, which confirms the non-stationary component in the time series mainly due to the characteristics of the environment (described above, see conceptual model) and to variations in flow caused by water withdrawals from the river (mainly due to groundwater pumping in the vicinity of the river). Indeed, the values of H indicate a process far from stochastic behavior (H~0.5) and resemble the Red noise (also called Brownian noise) which gives additional information about the nature of the hydrological processes in the studied basin.

The Júcar River in a practically natural regime (prior to 1974) was a clearly gaining river mainly due to the significant water amounts stored in the MOA. The groundwater levels were quite shallow and the connection between the groundwater and the Júcar River was upstream of the Puente de Carrasco site, being the gaining river until its exit from the MOA domain. In the following decades, groundwater extraction for irrigation began to increase significantly [59]. The groundwater levels in the MOA decreased in the area surrounding the Júcar River and, therefore, the natural inflows to the Júcar River. As a result, the river flows decreased in this section, making it necessary to increase the flows released from the Alarcon reservoir to avoid a significant discharge decrease in certain sections, i.e., Reach 2 close to Cuasiermas, years 1995 and 2008 [60,61]. Different simulations of the 1D-HD model indicate that with discharges released from Picazo gauging station lower than 1.83 m³s⁻¹ the Júcar River may dry up between chainages 71,000 m and 80,000 m. This effect was observed during model calibration for the period influenced by OPAD measures. During this period, due to the environmental restrictions imposed by the river basin authorities (CHJ), no direct water abstractions were made from Júcar River and a perimeter of influence was established around it where groundwater pumping that could influence river flows was not allowed. In addition to this situation, and because of the scarcity of water in the headwaters of the system (Alarcón reservoir), the discharges released to Júcar river were relatively low, more specifically between 2 and 5 m³s⁻¹ at El Picazo gauging station (for a summary, see [23]). These are the cases when hydrodynamic modeling is considered essential to establish the minimum discharge’s stream flow so that the river would not dry up. As can be seen in the results of the model shown in Figure 9, for Reach 2 (net losing river), in 95% of the cases when the discharge at El Picazo gauging station was less than 1.97 m³s⁻¹, the river would dry out close to Cuasiermas. In contrast, reach 3 in 90% of the cases the river would gain a base flow of 1.68 m³s⁻¹.

5. Conclusions

The 1D-HD model developed in the present research covers a period of 45 years in 121 km Júcar River reach as it passed within the MOA domain, built with MIKE 1D-HD module, it contributes significantly to the understanding of the hydrological functioning of the middle sub-basin under different discharge policies, including in conditions of water extractions from the river, either by direct withdrawals or by the influence of groundwater pumping. In this sense, the model provides simulated time series for the studied period with daily data of discharge and water depths of Júcar River, additionally, variables such as water level, velocities, and shear stress can be obtained along the studied reach. Statistical evaluations indicate that the model calibration for daily Q data (based on [42]) can be rated as “Very Good” and for h data (based on [54]) can be rated as “Excellent”.

Therefore, the 1D-HD model presented here becomes an essential tool for the management of water resources in the study area, since it not only allows hydrological balances to be made for shorter reaches of river (between gauging stations) but can also be very useful for estimating or predicting with more accurate spatio-temporal precision the critical flow for the river to maintain adequate ecological flow, established in relation to the reach where the river and the aquifer are disconnecting.

Additional information as for example riverbed permeability, Manning coefficient, and other parameters are needed to better determine the interaction between the river and the aquifer. The study of Jucar River behavior within the MOA domain offers an opportunity to further improve the modeling process. From a broader perspective, gaining knowledge about the interactions between rivers and groundwater bodies and their long-term evolution can be critical and can offer useful information for the managers of water resources, where similar problems exist (semiarid climate effects, overexploitation, and stakeholders’ conflicts). Future environmental management decisions can be supported in their efforts to minimize extreme climate impacts and adverse anthropogenic pressures, aiming to maintain environmentally beneficial hydrogeological balance.