*2.2. Water Quality Models*

PATRICAL [28,30] is a large-scale, conceptual model, with a monthly time step, that discretizes the territory with a resolution of 1 km × 1 km. The water quality component simulates nitrate transport through the hydrological cycle in the entire basin. This model includes the SW–GW interaction as it takes into account irrigation returns that recharge aquifers, lateral transfers among aquifers, and water movement through the river network. However, PATRICAL only reproduces part of the altered hydrological cycle, as it does not include the management of infrastructure or the modifications produced in the flow regime. Inputs to PATRICAL are monthly pluvial precipitation; air temperature; urban and industrial discharges to the GW bodies; nitrogen surplus in the soil; and GW withdrawals [39]. The data set employed in the PATRICAL model are shown in Table 1. Outputs of PATRICAL are streamflow-accumulated time series, GW levels, and total nitrate loads from diffuse pollution in rivers and aquifers. The schematic with the steps carried out by the model is shown in Figure 2a and detailed in the Appendix A. A more extensive description of PATRICAL model and the parameters used is provided by the Refs. [28,30]. **Table 1.** Data set employed in the PATRICAL and RREA models. NO3 −-SW: nitrate concentration in surface water (mg/L); NO3 <sup>−</sup>-GW: nitrate concentration in groundwater; Q: streamflow (m3/s); P: pluvial precipitation (mm); T: temperature (◦C); N-soil: nitrogen surplus in soil (KgN/ha); V discharge: point discharge volume (m3/year); PE: population equivalent.


**Figure 2.** PATRICAL (**a**) and RREA (**b**) models' structure and variables. Evaluation of simulation performance metrics (**c**). Rectangles with smoothed edges are the input variables to the model, the rectangles represent storages, and the document flowchart symbol represents the outputs of the process. DP: diffuse pollution; GW: groundwater; K: pollutant degradation constant N-water: nitrogen in the water resources; SW: surface water; t: time of residence of nitrate.

The information obtained from PATRICAL is the starting point for the second largescale surface-water-quality model, RREA. The two models complement each other, as RREA allows to include reservoir management and measurement regulation, agricultural and urban demands, and changes in the streamflow regime effects. A series of algorithms were developed using Python software [43] to automate decumulation loads and streamflow processes (PATRICAL output).

RREA estimates concentrations of pollutants in surface water bodies considering the load contributed to each SW rivers, the pollution coming from upstream and the possible degradation occurring in the water body itself. Input variables to RREA are: physical characteristics of the hydrographic network; water demands [39]; streamflow records of rivers and reservoirs; diffuse nitrate load (output PATRICAL); streamflow time series (output PATRICAL); point discharge; and degradation constant by pollutant. The data set employed in the RREA model is shown in Table 1. Point sources of nitrate were entered into the model by linking the authorized discharge location of wastewater treatment plants (WWTP) and the SW rivers into which they discharge. Output variables are the time series of streamflow and nitrate concentration circulating through the SW rivers under conditions altered by human activities. The general scheme of RREA is shown in Figure 2b and detailed in the Appendix A.

#### *2.3. Calibration*

The parameters were calibrated by an iterative process taking into account the following: (1). To assess the skill of the models to simulate the nitrate status of the water bodies; (2). to estimate the statistical error in order to obtain a greater number of SW rivers with satisfactory performance in the simulation of streamflow and nitrate concentration; and (3). to represent nitrate load generated by point and diffuse pollution.

In previous works, the hydrological component of PATRICAL was calibrated and validated by Pérez-Martín et al. [30], who reported satisfactory behaviour of the model for all evaluated water bodies. In addition, improvements have been made to the groundwater component and the SW–GW interactions, finding better fits between the simulated and observed flows with respect to the previous calibration [28].

The results of the model under an altered regime were compared with the observed streamflow and nitrate concentration (SIA Júcár, Available online: ps.chj.es/siajucar/, accessed on 26 March 2021) in the calibration process. The Python software [43] was used to calculate the main descriptive statistics (25%, 50%, and 75% quantiles, mean and standard deviation). The evaluations used the median of observed and simulated data in the SW rivers for greater robustness and to avoid outliers. To check the consistency of the data, automatic graphs were generated with the time series of the nitrate concentrations and streamflow in each SW river.

The statistical error was calculated using three indicators. First, the relative bias (PBIAS) shows the simulation deviation expressed as a percentage. In addition, it differs from other indicators because it has a specific classification for streamflow and water quality components. The second is the Nash-Sutcliffe efficiency (NSE), which determines the relationship between the error variance of the simulated data and the variance of the observed data [44]. The NSE ranges from −∞ to 1 and the optimal value is 1. Finally, the indicator Modified Kling-Gupta Efficiency (*KGEM*) (Equation (1)) decomposes the bias into three different terms, *r* represents the correlation coefficient between the simulated and observed time series, *β* is the ratio between the simulated and observed means (μ) (Equation (2)), and *γ* is the ratio of the coefficients of variation of both time series (Equation (3)). The optimal value for each of the three components of the *KGEM* is 1 [45,46].

$$KGEM = 1 - \sqrt{(r - 1)^2 + (\beta - 1)^2 + (\gamma - 1)^2} \tag{1}$$

$$
\beta = \frac{\mu\_{sim}}{\mu\_{obs}} \tag{2}
$$

$$\gamma = \frac{CV\_{\text{sim}}}{CV\_{\text{obs}}} \tag{3}$$

#### *2.4. Nitrate Status Classification Performance*

To assess the skill of the models to simulate the nitrate status, a 2 × 2 contingency table for dichotomous events was used [35]. This table allows assessing the performance of the models to evaluate the status of water bodies based on the nitrate concentration values. For this purpose, nitrate status was classified in the complete time series of simulated and observed data for each SW river, considering the threshold value of 25 mg NO3 −/L [47], and using the same length of data in both series. In this way, a matrix of discrete nonprobabilistic values was obtained as shown in Table 2.

**Table 2.** Contingency table to assess nitrate status classification performance of PATRICAL/RREA models.


Four different measures were used to assess the skill of the models to simulate nitrate status: the Accuracy (*ACC*) assesses the model performance to reproduce an event correctly and was calculated using Equation (4), *ACC* ranges from 0 to 1, and 1 is the best value; the bias measured is the ratio of the simulated mean and observed mean Equation (5). Bias ranges from 0 to infinite, and 1 is the best value; the Success Ratio (*SR*) provides information on the proportion of *TP* in the whole time series (Equation (6)) [35,48]; and in contrast, specificity (*SP*) which is the proportion of *TN* correctly classified in the simulation (Equation (7) [49]. For the indicators *SR* and *SP*, the best value is 1 and the worst is 0.

$$\text{ACC} = \frac{TP + TN}{TP + FN + FP + TN} \tag{4}$$

$$IAS(TC) = \frac{TP + FP}{TP + FN} \tag{5}$$

$$SR = \frac{TP}{TP + FP} \tag{6}$$

$$SP = \frac{TN}{TN + FP} \tag{7}$$

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

### *3.1. Calibration*

Streamflows and nitrate concentrations were jointly calibrated in the six main water resource systems of the Júcar RBD. The values obtained for the three statistical indicators are shown in Figure 3. According to the PBIAS indicator, the streamflow calibration provided a good fit between simulated and observed values in the Mijares, Turia, Júcar, and Vinalopó; and satisfactory fit in Palancia and Serpis. For nitrate concentration, a very good fit was obtained in Turia, Júcar, Serpis; a good fit in Palancia, and a satisfactory fit in Mijares and Vinalopó.

Based on the NSE values for the monthly streamflow, the fit was satisfactory in Mijares, Turia, and Júcar, whereas in Palancia, Serpis, and Vinalopó the performance was unsatisfactory. NSE values for the nitrate concentration in Mijares, Palancia, and Vinalopó were below zero; whereas in Júcar, Turia, and Serpis, positive values were obtained, which indicates better behaviour of the model in the simulation of nitrate concentration.

The *KGEM* indicator and the three components in the streamflow performance was close to the optimum in most of the systems evaluated, except in Vinalopó (Figure 3c), which also presented a ratio between coefficients of variation (*γ*) close to zero.

**Figure 3.** Evaluation parameters of the calibration process of the streamflow (altered regime) (**a**) and nitrate concentration (**b**). The *KGEM* components for streamflow (**c**) and nitrate concentration (**d**) in the main systems of the Júcar RBD. r = correlation coefficient; *β* = bias ratio; *γ* = ratio of the coefficients of variation; *KGEM* = Modified Kling-Gupta Efficiency.

*KGEM* values for nitrate concentration were between 0.3 and 0.7 in the Júcar, Mijares, Palancia, Turia, and Serpis (Figure 3d); whereas in Vinalopó, a value close to zero was obtained, with similar behaviour to that found in the streamflow. Analysing *KGEM* components (Figure 3d), the correlation coefficient (r) was 0.81 for Mijares and 0.28 for Vinalopó, meaning that simulated and observed data series are more correlated in Mijares than in Vinalopó. The bias ratio (*β*) was 1.59 in Mijares and 0.40 in Vinalopó, so nitrate concentrations are overestimated in Mijares, while it is underestimated in Vinalopó. Júcar, Palancia, Turia, and Serpis have a bias relation close to the optimum. The ratio between the coefficients of variation (*γ*) are close to optimal in Júcar, Palancia, and Serpis, and presented values between 0.6 and 0.52 in Mijares and Vinalopó, respectively. The NSE index for Mijares was not satisfactory but there was a high correlation between simulated and observed data as a satisfactory *KGEM* value was obtained.

The models performed well in the simulation of water resources in basins with large surface areas (such as Júcar and the Turia), but in small basins, with less surface area and less flow (such as the Vinalopó), the fit was less satisfactory. This is influenced by the greater number of gauging stations and measurements in the basins with a larger area.
