Assessment of Water Quality Regulation Functions in Southwestern Europe Watersheds

: Estimation at large scale of the water quality regulation services is still lacking. It is essential to develop methodological approaches to quantify nutrient-related functions’ distribution. The present study aims to quantify nitrate-related ecological functions through nitrate net balance (NNB), nitrate removal (NR), and nitrate production (NP). This study explores the spatiotemporal dynamics of these indicators in South-Western Europe (SUDOE, 216 subsystems over 81 basins) at a monthly scale from 2000 to 2010. We use the Soil and Water Assessment Tool model to simulate nutrient transfer at the subsystem scale (~3000 km 2 ) and calculate ecological functions. The modeled NNB is validated at the subsystem scale by comparing with NNB predicted at the water body scale (~60 km 2 ) over the Garonne watershed (France). Hot spots of NR are located in the south of SUDOE, characterized by a warmer and dryer climate, whereas NP hot spots are located in the most anthropized streams. The mean NNB (the balance between NP and NR) at the subsystem scale for the SUDOE territory reaches − 3.5 and − 2.8 gN m − 2 day − 1 during the most active seasons. The results highlight drivers inﬂuencing NR such as streamﬂow, river slope, and hydrological alteration. Getting an overview of where and when these nitrate regulation functions (NP and NR) occur is essential for socio-ecosystem sustainability and this study highlights its sensitivity to anthropogenic stressors.


Introduction
The increase in the world's population together with living standard improvements, consumption pattern changes, and the expansion of irrigated agriculture is driving rising global water demand [1,2]. The associated increase of water withdrawals from rivers and aquifers causes river streamflow decreases [3] which lead to the loss of cultural and amenity services linked to the river [4][5][6], the loss of river ecosystems, and modification of their trends [7]. Freshwater scarcity is known to have impacts on nutrient cycles [8,9], often leading to water quality degradation due to a lower dilution of pollutants [10,11]. To prevent water scarcity and pollution increase, the Water Framework Directive (EU-WFD) has been set up in Europe [12]. Each European country applies a set of laws at various scales of the watersheds in order to comply with the EU-WFD.
However, little is known about the impacts of water scarcity on river ecological functions related to water quality regulation. This is especially true at high spatial and temporal resolution over regional scales, despite the essential role of these functions for sustaining the ability of the environment to supply ecosystem services [13]. If different definitions of ecological function exist, we consider here that ecological function gathers all the physical, biological, and biogeochemical processes that contribute to the ecosystem The hydrology is highly impacted by human activities such as urban, agricultural, and industrial activities or river management structures. More than 50% of the 60 million inhabitants are mostly settled in major cities such as Madrid, Lisbon, and Toulouse ( Figure  1). During droughts, cities such as Murcia and Valencia are frequently facing drinking water supply issues [46,47]. The SUDOE is the major agricultural region in Europe with 50% of farming lands (CLC 2012). Agriculture is dominated by cereal, maize, and forage crops. Vegetables, olives, and wine production occupy a large surface area, especially in Spain (29%) and Portugal (31%). These crops are potentially major water consumers. Irrigated cotton crops are also widespread in the Southern part of SUDOE, especially in the Andalusian region. The other half of the SUDOE territory is 46% forests and semi-natural areas, 3% artificial areas, and 1% wetlands and water bodies. Water bodies are highly anthropized by reservoirs and canals for irrigation, city water supplies, and hydroelectricity. The hydrology is highly impacted by human activities such as urban, agricultural, and industrial activities or river management structures. More than 50% of the 60 million inhabitants are mostly settled in major cities such as Madrid, Lisbon, and Toulouse ( Figure 1). During droughts, cities such as Murcia and Valencia are frequently facing drinking water supply issues [46,47]. The SUDOE is the major agricultural region in Europe with 50% of farming lands (CLC 2012). Agriculture is dominated by cereal, maize, and forage crops. Vegetables, olives, and wine production occupy a large surface area, especially in Spain (29%) and Portugal (31%). These crops are potentially major water consumers. Irrigated cotton crops are also widespread in the Southern part of SUDOE, especially in the Andalusian region. The other half of the SUDOE territory is 46% forests and semi-natural areas, 3% artificial areas, and 1% wetlands and water bodies. Water bodies are highly anthropized by reservoirs and canals for irrigation, city water supplies, and hydroelectricity. For example, the SUDOE territory gathers at least 1974 dams which are commonly used for hydroelectricity and irrigation. Indeed, 13.5% of agricultural areas are permanently irrigated (https://ec.europa.eu/eurostat accessed on 10 October 2021). SUDOE soils are Water 2021, 13, 2980 4 of 21 mainly composed of Calcisol, Cambisol, and Regosol [48]. The present study focuses on the seven main watersheds.

Subsystem Scale Model
The model used in this study is the Soil and Water Assessment Tool, SWAT [49], an hydro-agronomic model simulating the nitrogen cycle and human stressors on watersheds. The SWAT model has been widely used at regional [50,51] and continental [36] scales. Many studies used this model to simulate surface water quantity and quality and to predict the environmental impact of land use, human stressors, and climate change on water quantity and quality [41] and also on ecosystem services [42]. The modeling approach has been set up at the scale of subsystems (from 203 to 24,963 km 2 ; Figure 1), which is the scale used for basin management planification and at the political decision-making level, i.e., the level of law application as described in [10]. Each subsystem is divided into Hydrological Response Units (HRU), characterized by a unique combination of soil, land cover, and slope information. Water balance, erosion, and nutrient cycles are calculated in each HRU, and summed at subsystem level. The SWAT model also integrates urban sewage, agriculture management, and dams. The implementation of the model at this subsystem scale, including the conceptualization, the calibration, and the validation steps of streamflow and nitrate loads, has been described in [10]. The present paper uses this model to quantify and validate nitrate-related ecological functions (NNB, NR, and NP) at the subsystem scale. This model will be referred as the "subsystem model" in the rest of the paper.
The subsystem model has been calibrated from 2000 to 2010 and validated from 1990 and 1999 over the territory with an average NSE, R 2 , and PBIAS, respectively, equal to 0.52, 0.78, and ±12.52 for streamflow, and 0.51, 0.67, and ±19.59 for the nitrate load. The Nash-Sutcliffe efficiency (NSE) is calculated as one minus the ratio of the error variance of the simulation divided by the variance of the observation. The coefficient of determination (R 2 ) is the variation between observation and simulation. The percent bias (PBIAS) is calculated to measure the interannual monthly average tendency of the simulated values to be larger or smaller than observations. According to the criteria provided by [52], the streamflow and nitrate load simulation of the subsystem model can be considered as satisfactory over the entire territory.
The river nitrogen cycle is modeled by the QUAL2E model [53]. Dissolved oxygen concentration drives the biological processes as in-stream algae production, carboneous oxygen uptake, and nitrogen transformation. The QUAL2E module also integrates benthic algae processes. Nitrogen concentrations are controlled by biological-mineralization, ammonification, nitrification, and denitrification-and physical processes-settling and dilution. The physical processes consist in atmospheric aeration, storage of stream water and nutrients during a time step. All the processes are represented by differential equations described in the QUAL2E manual [53].
Four water quality regulation functions are investigated: the nitrate net balance (NNB), the nitrate removal (NR), the nitrate production (NP), and the nitrate net balance rate (NNBR).
The NNB indicator is the nitrate net balance that gathers the NR and NP processes with the following equation: where NO 3 load is gN day −1 and "Reach wetted area" is the length of the reach (m) multiplied by the width (m). The NNB (gN m −2 day −1 ) is the amount of nitrate that stays in or leaves the water bodies of the reach during each time step. If NNB is negative, nitrate is removed from the stream by denitrification, retention in sediment, or plant and algae assimilation, indicating a NR higher than the NP component. The areas with negative NNB define sink zones where the self-purification capacity of the river is high. When NNB is positive, the overall nitrate load is increasing in the free water indicating that anthropic and/or natural inputs are dominant compared to retention processes thus resulting in nitrate production (NP). Subsystems with positive NNB enable us to underline areas and periods when high nitrate production occurs and to identify nitrate sources (natural or not). A comparison between modeled and measured NNB enabled to validate the SWAT modeling of NNB at the water body resolution [54]. This study aims at validating the definition and quantification of NNB by changing the spatial resolution of the model from the water body resolution [54] to the subsystems resolution, before applying this methodology to all the SUDOE territory.
Another indicator is the nitrate net balance rate (NNBR): is the ratio between NNB (gN m −2 day −1 ) and the nitrate load that enters into the reach (gN day −1 ). The streamflow effects are erased in this indicator, which can have an important effect on reach efficiency and seasonal analysis.
Two runs of the model are performed: a first run with activated QUAL2E module, and a second run with inactive QUAL2E module (suppressed instream biological processes). The difference between the predicted values of NNB, NR, NP, and NNBR obtained with the two runs allows us to quantify the contribution of biological processes to ecological functions (NNB, NR, NP, and NNBR). More details about the activation and deactivation of the module are presented in [54]. The Garonne River basin is the principal watershed in South-Western France with an area of 60,000 km 2 and a river length of 600 km ( Figure 1). The watershed has a high temporal and/or spatial heterogeneity in its topographic, geologic, and climatic conditions. The climate is Mediterranean in the East of the watershed, and Oceanic in the West. The annual rainfall average is 900 mm, with a maximum of 2000 mm in the upper part of the watershed. The Strahler order is 8 at the Garonne River outlet [55]. The Garonne flow regime is a pluvionival regime characterized by two high flood periods in February and May and low flow from August to September. Tonneins is considered to be the outlet of the watershed (monitoring point "Outlet" in Figure 1) because it is the last gauging station with no tidal influence. The streamflow at Tonneins varies from 53 to 4740 m 3 s −1 with an annual mean of 496 m 3 s −1 .
The Garonne watershed is composed of 60% agricultural lands, 32% forest, and 2.5% artificial areas. Agriculture is dominated by maize and wheat crops and is located in the lowland plain. The area is dominated by Cambisols, favoring agriculture practices, and the floodplain is characterized by Fluvisols and Luvisols [56].

Methodology: Comparing Subsystem and Water Body Scale Models
To validate the quantification of nitrate-related ecological functions at the SouthWestern Europe scale, model outputs (river flow, nitrate loads, and NNB) at subsystem scale (this study;~3000 km 2 ) were compared to those obtained at water body scale (~60 km 2 ) on the Garonne watershed in the South of France [54]. Even if data on nitrate-related ecological functions are sparse and rare, results from [10] were validated by in situ measurements of nitrate removal intensities in different streams.
The subsystem and water body models used the same input data. The main difference between the two models is the spatial resolution and how the subbasin division is achieved based on river network and topography. The Garonne basin model, explained in [54], is discretized at the administrative level of water bodies, the level used by stakeholders to monitor stream water quality and quantity, with a surface area varying from 1 to 208 km 2 ( Figure 1). This model will be referred to as "water body model" in the following work. The streamflow and nitrate load simulated with the water body model are calibrated from 2000 to 2010 and validated from 1990 to 1999 using a similar procedure as the subsystem model. According to the criteria provided by [52], the streamflow and nitrate load simulation of the water body model can be considered as satisfactory over the entire Garonne basin with NSE, R 2 , and PBIAS, respectively, equal to 0.51, 0.62, and ±6.0 for the streamflow, and 0.41, 0.68, and ±9.69% for the nitrate load.
Simulations of streamflow, nitrate load, and NNB were compared on 25 subbasin outlets over 10 years of simulations (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). Model outputs were also compared with observations at four gauging stations: Tonneins (Outlet), Villefranche (K), Portet (R), and Roquefort (W). When the dataset distribution was normal (following a Gaussian law), a Student t-test was used to investigate significative differences between groups. When non-normal, the nonparametric Mann-Whitney U test was used. The comparison of the simulated NNB from the subsystem and water body models were compared at the node of the subsystem scale ( Figure 2). The NNB comparison was made for annual and monthly datasets. A good estimation procedure should result in PBIAS values close to zero [52]. The similarity of the results obtained using the two resolution models was used to validate the application of the SWAT model at subsystem scale to quantify nitrate-related ecological functions over the SUDOE territory.
The subsystem and water body models used the same input data. The main difference between the two models is the spatial resolution and how the subbasin division is achieved based on river network and topography. The Garonne basin model, explained in [54], is discretized at the administrative level of water bodies, the level used by stakeholders to monitor stream water quality and quantity, with a surface area varying from 1 to 208 km² ( Figure 1). This model will be referred to as "water body model" in the following work.
The streamflow and nitrate load simulated with the water body model are calibrated from 2000 to 2010 and validated from 1990 to 1999 using a similar procedure as the subsystem model. According to the criteria provided by [52], the streamflow and nitrate load simulation of the water body model can be considered as satisfactory over the entire Garonne basin with NSE, R 2 , and PBIAS, respectively, equal to 0.51, 0.62, and ±6.0 for the streamflow, and 0.41, 0.68, and ±9.69% for the nitrate load.
Simulations of streamflow, nitrate load, and NNB were compared on 25 subbasin outlets over 10 years of simulations (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010). Model outputs were also compared with observations at four gauging stations: Tonneins (Outlet), Villefranche (K), Portet (R), and Roquefort (W). When the dataset distribution was normal (following a Gaussian law), a Student t-test was used to investigate significative differences between groups. When non-normal, the nonparametric Mann-Whitney U test was used. The comparison of the simulated NNB from the subsystem and water body models were compared at the node of the subsystem scale ( Figure 2). The NNB comparison was made for annual and monthly datasets. A good estimation procedure should result in PBIAS values close to zero [52]. The similarity of the results obtained using the two resolution models was used to validate the application of the SWAT model at subsystem scale to quantify nitrate-related ecological functions over the SUDOE territory. Figure 2. Example illustrating the difference between the subsystem and water body (WB) scales of the models applied to the Garonne River watershed. One subsystem integrates several water body units. Comparisons of nitrate net balance (NNB) are performed on common nodes between subsystem and WB scale models to validate the calculation of NNB at the subsystem scale over SouthWestern Europe.

NP and NR Driver Selection
A recent study highlighted that environmental factors such as river slope and streamflow might impact NP and NR functions along the Garonne River [54]. According to previous studies and local knowledge [3,24,52], a panel of qualitative and quantitative environmental variables were selected to proceed to the investigation of potential drivers of NNB (Table 1). These influencing factors were chosen to illustrate both environmental variability and human stressors. A principal component analysis of mixed data (PCAMix), Figure 2. Example illustrating the difference between the subsystem and water body (WB) scales of the models applied to the Garonne River watershed. One subsystem integrates several water body units. Comparisons of nitrate net balance (NNB) are performed on common nodes between subsystem and WB scale models to validate the calculation of NNB at the subsystem scale over SouthWestern Europe.

NP and NR Driver Selection
A recent study highlighted that environmental factors such as river slope and streamflow might impact NP and NR functions along the Garonne River [54]. According to previous studies and local knowledge [3,24,52], a panel of qualitative and quantitative environmental variables were selected to proceed to the investigation of potential drivers of NNB (Table 1). These influencing factors were chosen to illustrate both environmental variability and human stressors. A principal component analysis of mixed data (PCAMix), which is a mix between an ordinary principal component analysis and a multiple correspondence analysis, was used on the drivers (variables list in Table 1) and NNB over the SUDOE territory at the subsystem scale. The PCAMix analysis aims to find correlations between input data, named "observations", quantitative variables, named "numerical variables", and qualitative variables, named "level". R version 3.5.1 [57] with the PCA mixdata package [58] was used to complete this statistical analysis. Part of SUDOE territory ( Figure 1) -

Riparian
Percentage of riparian area in the subsystem % Slope Average slope of rivers in the subsystem % Urb Percentage of the artificial areas in the subsystem % Wetl Percentage of the wetland areas in the subsystem % Width River width m Wtr Percentage of the water body areas in the subsystem % The drivers may change with a larger panel of pedoclimatic conditions and land occupation. For these reasons, clustering using the K-means algorithm [59] on the NNB values obtained in each subsystem allowed us to split streams into three groups. Nonparametric statistical analysis was used due to non-normal datasets. The Mann-Whitney U statistical test was used to identify significant differences in drivers (river slope, streamflow, stream width, stream length, and altitude) variations between groups.
As the level of anthropization of a subsystem might influence the NNB function, this study also looked for the link between NNB and an indicator of human stressors, the widely used Index of Hydrological Alteration (IHA) [60][61][62][63]. A lot of literature has explored the impact of humans on stream hydro alterations, including through this indicator. The IHA characterizes the impact of river regulation on the flow regime. It was developed by the Nature Conservancy organization [60][61][62] and modified by [63]. The IHA is the sum of 12 indicators that reflect the modification of the current human-impacted flow regime from the natural one, with different statistics (magnitude index, variability index, flood and drought intensity, and frequency indexes) [60][61][62]. The IHA was calculated at a monthly scale following the methodology of [10]. The IHA varied between 0 and 100 with the maximum value for unaltered river flow. The correlation between IHA and NNB has been investigated on 54 subsystems over SUDOE territory.   The spatial distribution of NNB in the Garonne watershed is compared in Figure S1 (Supplementary material). The PBIAS values between water body model (high-resolu- The simulated monthly interannual average of NNB by subsystems ( Figure 4) are not significantly different between the two models throughout the year (p-value = 0.47, W = 16,368 and PBIAS = 23%). The similarity of the NNB monthly interannual average between both models demonstrate that the modeling resolution change does not affect the nitrate-related ecological function representation in the SWAT model. The annual average of N removal is 8.5 gN m −2 y −1 with the subsystem model and 7.4 gN m −2 y −1 with the water body model. The highest quantity of nitrate is removed during the summer (from June to August). The only slight difference between the two models is observed between July and September and November-December.

Comparison between High and Low Resolution SWAT Models for the Garonne Watershed
The spatial distribution of NNB in the Garonne watershed is compared in Figure S1   Comparison of (A) monthly streamflow (m 3 s −1 ) and (B) monthly nitrate load (kg) outputs between two model resolutions (the water body model (olive, high resolution) and subsystem model (grey, low resolution)) and observations (white) over Garonne watershed at the 25 monitoring points. Letters correspond to the different stations in Figure 1. The spatial distribution of NNB in the Garonne watershed is compared in Figure S1 (Supplementary material). The PBIAS values between water body model (high-resolution) and subsystem model (low-resolution) vary between −62.0% and +49.3% with an average of +27.4%. The subsystem model underestimates this indicator in small streams of the upper watershed part in the Pyrenees and Massif Central chains, but has a good estimation in the rest of the basin.
The differences of NNB simulations over the Garonne watershed at the monthly time step between water body and subsystem resolution models follow a normal distribution and are not significant (Student's t-test, t = 0.71, p-value = 0.48). Considering that streamflow and nitrate load are validated with the high resolution water body model on the Garonne watershed [54] and with the low resolution subsystem model [10], the absence The differences of NNB simulations over the Garonne watershed at the monthly time step between water body and subsystem resolution models follow a normal distribution and are not significant (Student's t-test, t = 0.71, p-value = 0.48). Considering that streamflow and nitrate load are validated with the high resolution water body model on the Garonne watershed [54] and with the low resolution subsystem model [10], the absence of differences between the NNB calculated with the two model versions for almost all months enables us to analyze the spatial and temporal variation of NNB over the entire SUDOE territory. Table 2 shows the seasonal variability of NNB and NNBR in the SUDOE territory, and the physical and biological components for each variable (see Section 2.4). NNB values are lower in summer and spring with −2.8 and −3.5 gN m −2 day −1 , respectively, indicating a stronger retention during these seasons. Physical processes for NNB only cause NR, with higher values also predicted during spring and summer, and lower values in winter. The variability of biological processes for NNB is not significant between seasons (t [−1.01, 0.136], p-value > 0.05). The highest values of NNBR, NP, and NR are predicted during summer with −631.5 km −2 day −1 , 30.0 km −2 day −1 , and −186,185 km −2 day −1 , respectively. For NNBR, physical processes occur mainly during summer (−626.5 km −2 day −1 ) and spring (−528.7 km −2 day −1 ), and biological processes during autumn (with an average of −8264 km −2 day −1 ). The highest production rate and the lowest removal rates occur during winter.

Spatial and Seasonal Variations of Nitrate-Related Ecological Functions over SUDOE Territory
The interannual spatial variation of NNB in each subbasin calculated between 2000 and 2010 is shown in Figure 5A. Some rivers remove nitrate with more or less efficiency whereas a minority of rivers in SUDOE exhibit a production of nitrate. Rivers producing nitrate are mainly located in the North of the territory in Occitanie and Cataluña regions. The highest NP occurs in the main large streams whereas the highest NR happens in the intermediate reaches of the streams. Physical processes retain more nitrate at the outlet of major rivers such as Tajo, Garonne, and Guadiana ( Figure 5C), whereas biological processes produce nitrates in these same locations ( Figure 5B). Streams located in the Southern part of the territory retain more nitrate. Intermediate river reaches in the North of the territory are hot spots of NR for both physical and biological processes, which is consistent with the results obtained at the waterbody resolution for the Garonne River [54].

NNB and Drivers in the SUDOE Territory
Several principal component analyses of mixed data, called PCAMix, were performed to investigate correlation between NNB, NR, and influencing factors. First, a PCAMix was set up considering all the variables of Table 1. The first five dimensions gathered only 33% of the information. NNB does not influence the first dimension (0.03) and 22% of the NNB is explained on dimension 3. Dimension 3 also explains streamflow (81%) and nitrate load (67%). As expected, the streamflow and the nitrate load are directly linked to NNB. However, NNB is not only explained by these two variables.
A second analysis is run without streamflow, nitrate load, and variables contributing less than 10% on dimension 3 of the first PCAMix. The variables selected for the second PCAMix are in bold in Table 1. This second analysis gathered 31% of the information. The first and the third dimensions explain most of the variables. Figure 6 illustrates the trends on the first and third dimensions of the observation ( Figure 6A), the levels ( Figure 6B), the numerical variables ( Figure 6D), and all the variables ( Figure 6C). The three geographic groups (illustrated in Figure 1) contributed differently on the PCAMix axis and are correlated to specific land covers ( Figure 6B). The East part (E) is characterized by almonds, olive, orchard, and orange crops, the West part (W) is composed of cotton, rice, and intensive cereal cultures, and the North part (N) is associated with other land uses including forests, seminatural areas, urban areas, and wetlands. These results seem consistent with the SUDOE crop repartition with intensive crops in the South of the territory. NR seems to be positively correlated with elevation, river slope, and forest land use percentage ( Figure 6D) whereas NNB is linked to agriculture and urban land use proportions ( Figure 6D). Finally, it appears that the denitrification rate in soil simulated by the SWAT model is linked to the fertilization intensity ( Figure 6D). Figure 6C shows the PCAMix considering all qualitative and quantitative variables. Denitrification rate seems related to land use whereas NR and NNB have the same trends and are influenced by the percentages of land use (water bodies, agricultural, and forest, urban), river slope, and elevation. However, the dimensions 1 and 3 gather only 16% of the dataset information. Figure 7a shows the distribution of NNB values for the three geographical regions of the SUDOE territory: East, North, and West. The Eastern part of the SUDOE territory removes more nitrate than the Western part, while the Northern part has a high intersubsystem heterogeneity. groups (illustrated in Figure 1) contributed differently on the PCAMix axis and are correlated to specific land covers ( Figure 6B). The East part (E) is characterized by almonds, olive, orchard, and orange crops, the West part (W) is composed of cotton, rice, and intensive cereal cultures, and the North part (N) is associated with other land uses including forests, seminatural areas, urban areas, and wetlands. These results seem consistent with the SUDOE crop repartition with intensive crops in the South of the territory.   Figure 7b. Group 1 is characterized by NP (Figure 7b). Group 2 gathers rivers presenting low processes with NP and NR values close to 0. NNB is significantly different depending on groups (p-values < 0.005, X = 119 960). Group 3 is characterized by high NR with an average NR value of −4.65 gN m −2 day −1 for physical processes and −0.25 gN m −2 day −1 for biological processes. Group 1 is characterized by a production of nitrate by biological processes with an average NP of 11.24 gN m −2 day −1 , significantly different from groups 2 and 3. Group 2 removes nitrate in lower quantities (average NR value lower than −1.28 gN m −2 day −1 ) than group 3 (average NR value = −24.25 gN m −2 day −1 ). These differences between groups are significant and are shown in Figure 7b. Figure 7c-g illustrates the analysis of key factors (streamflow, stream length and width, river slope, and altitude) by groups in order to find a common tendency between groups. Streamflow seems to be the key driver for NNB with a significant difference between groups (Figure 7c), which confirms the first PCAMix analysis. Group 1 (NP) gathers streams with high streamflow whereas group 3 (higher NR) is characterized by streams with low discharge. This study analyzes river width ( Figure 7d) and length (Figure 7e). It appears that NR is higher for narrower and smaller streams (group 3) and this difference is significant with the other groups. Flat reaches are gathered in group 1 and sloping reaches in group 3 (Figure 7f). The river slope difference is not significant between groups 2 and 3, both characterized by NR function (Figure 7f). For elevation (Figure 7g), group 1 gathers streams with the lowest altitude and river slope, group 3 contains the highest and sloping reaches, while group 2 has the highest altitude (Figure 7g).

Relation between IHA and Nitrate Net Balance (NNB)
This study explores the variations of the IHA according to the logarithm of the averaged NNB on each subsystem. The 54 studied subsystems are dominated by nitrate removal (NR) with only negative NNB. Figure 8 shows a negative correlation between IHA and the absolute NNB logarithm in 54 subsystems over the SUDOE territory (R 2 = 0.35 and p-value < 0.005). The direct representation of NNB by the IHA values is not possible but there is a clear correlation between the two variables. Nitrate removal increases with decreasing IHA. The spatial distribution of subsystems over the SUDOE territory shows that the Eastern part ( Figure 1) has subsystems with the highest NR but the lowest IHA values.

Relation between IHA and Nitrate Net Balance (NNB)
This study explores the variations of the IHA according to the logarithm of the averaged NNB on each subsystem. The 54 studied subsystems are dominated by nitrate removal (NR) with only negative NNB. Figure 8 shows a negative correlation between IHA and the absolute NNB logarithm in 54 subsystems over the SUDOE territory (R 2 = 0.35 and p-value < 0.005). The direct representation of NNB by the IHA values is not possible but there is a clear correlation between the two variables. Nitrate removal increases with decreasing IHA. The spatial distribution of subsystems over the SUDOE territory shows that the Eastern part ( Figure 1) has subsystems with the highest NR but the lowest IHA values.

Large Scale Application: Model Performances
The similar performance of the high and low resolution models in simulating streamflow, nitrate load, and NNB on the well-studied Garonne River watershed in the South of France constitutes a validation step to quantify NNB at the subsystem scale for pluri-regional applications. Only few studies have analyzed hydrology, water quality, and related ecosystem service indicators at such large scales [3,32,36]. For the first time, this study focuses on the quantification of NP and NR in several watersheds at a pluri-regional scale and includes a validation step by comparing two model resolutions on a well-studied watershed. The Garonne river watershed has been chosen as an ideal site to validate the upscaling approach because (i) NR has been measured in situ in different streams which constitute a rare dataset, (ii) a model has already been applied at a high-resolution scale and is available for comparison, and (iii) the surface area of the watershed is consistent (60,000 km 2 ). Even if this study shows a good correlation between NNB values at subsystem and water body scales, differences can still be notified during low flow (July-September) and high flow periods (November-December). Two reasons could explain these discrepancies. The first reason is that the spatial and temporal resolution of the process calculation over one subsystem is different if the model used is the water body (sum of NNB calculated for each water body included in the subsystem) or the subsystem model (directly based on the difference between output and input to the subsystem). This can lead to more spatial heterogeneity, but also to a mathematical underestimation of nitrate-based processes at the water body scale compared to the subsystem scale. The second reason is

Large Scale Application: Model Performances
The similar performance of the high and low resolution models in simulating streamflow, nitrate load, and NNB on the well-studied Garonne River watershed in the South of France constitutes a validation step to quantify NNB at the subsystem scale for pluriregional applications. Only few studies have analyzed hydrology, water quality, and related ecosystem service indicators at such large scales [3,32,36]. For the first time, this study focuses on the quantification of NP and NR in several watersheds at a pluri-regional scale and includes a validation step by comparing two model resolutions on a well-studied watershed. The Garonne river watershed has been chosen as an ideal site to validate the upscaling approach because (i) NR has been measured in situ in different streams which constitute a rare dataset, (ii) a model has already been applied at a high-resolution scale and is available for comparison, and (iii) the surface area of the watershed is consistent (60,000 km 2 ). Even if this study shows a good correlation between NNB values at subsystem and water body scales, differences can still be notified during low flow (July-September) and high flow periods (November-December). Two reasons could explain these discrepancies. The first reason is that the spatial and temporal resolution of the process calculation over one subsystem is different if the model used is the water body (sum of NNB calculated for each water body included in the subsystem) or the subsystem model (directly based on the difference between output and input to the subsystem). This can lead to more spatial heterogeneity, but also to a mathematical underestimation of nitrate-based processes at the water body scale compared to the subsystem scale. The second reason is related to the different resolution and origin of the forcings between water body and subsystem scales, e.g., the different scale of point sources datasets. Indeed, releases from seven waste water treatment plants (WWTP) have been implemented based on local datasets for the high resolution model, while only three WWTP were implemented based on a European dataset for the low resolution model. These differences can lead to discrepancies in direct nitrogen inputs to the river, which can be associated to the stimulation or reduction of nitrogen retention. Even if variations can be expected due to the different resolution of the forcings in the low and high resolution models, the relative agreement between NNB calculated at the two resolutions on the Garonne watershed shows that the upscaling of the SWAT model from water body scale to subsystem scale is possible. Nevertheless, further works are encouraged to pursue the validation of models at local and basin catchment scales to improve upscaling application of pluri-regional models.

NNB Variation in Time and Space
The analysis of NNB variations in time and space highlights hot spots and hot moments of both components of NNB: NP and NR. When only the Garonne basin is considered, this study shows that NR occurs mainly in intermediate reaches of the watershed during summer and spring when the biological processes are the most active, whereas NP occurs mostly during winter when in-stream processes are the least active, and in major rivers, downstream of important cities. These trends of NR and NP at the subsystem scale are consistent with those found at water body scale in the Garonne River watershed [54].
When all the SUDOE territory is considered, the NNB is higher during the spring and summer seasons when the ecological functions are active [64] and where processes are more active [38]. The hot spots of NR are indeed located in the South and the East of the SUDOE territory, where the climate is warmer and dryer. A warmer climate is known to activate elimination processes such as denitrification [65]. The hot spots of NP are located in the main streams where human pressures are higher [45] especially on biological processes. The Tajo, the Garonne, and the Guadiana basins are known to be highly impacted by human stressors such as intensive farming or arable lands, the massive use of irrigation, and the presence of multiple large dams on rivers [66]. The streamflow in Southern SUDOE rivers is lower compared to the North due to the specific pedoclimatic conditions with lower precipitations and larger water withdrawal for human consumption and uses. This induces an increase in residence time and thus in exchanges between the water column and the sediment [67]. NNB processes in these streams might be explained by the nutrient spiraling phenomena in the fluvial continuum [6].

Drivers for NR and NP with Different Resolutions
This study offers an initial key to understand the main drivers influencing NNB and NR. The hydrogeomorphological factors determine the residence time, depth, and exchanges with the hyporheic zone [29], which has a strong impact on the oxygenation of the environment and the kinetics of nitrate removal processes. The present study explores the relationship between these factors and NR obtained with low resolution. NR seems to be linked to streamflow, nitrate load, land cover, river slope, and elevation in agreement with results obtained at high-resolution water body scale [54]. Our results also show that nitrate retention does not occur in large streams with low elevation, which are net nitrate producers.
Some of the environmental drivers of NR will also influence NP function. NP functions are highly dependent on human stressors such as agriculture and urban effluents. Large cities and intensive agriculture are located downstream of the watershed in flatted sedimentladen zones [68]. Our study confirms this hypothesis with NP production occurring in large and long rivers located in areas with low elevation, low river slopes, and high discharge, which confirms the results of [54]. There is a clear distinction between streams characterized by NR and NP functions. The functional and structural gradient along the river network has an influence on NNB functioning. Mountain areas characterized by high elevation and high river slopes do not support NP and NR functions whereas NR are favored in intermediate streams and NP in large rivers. These results are similar to the ones found in [54], which show that influencing factors do not change by changing the resolution of the study.
Land cover seems to be an influencing factor of NR in SUDOE. The SUDOE territory, especially in the South of Spain and Portugal and particularly in the East of Spain, characterized by the highest NR (Figure 7a), is known to be more impacted by human stressors [10,69] and subject to water scarcity due to different pedoclimatic conditions such as Oceanic, Mediterranean, and Mountain climates [7,43]. The human-intensive use of water and land resources through agriculture, urbanization, dam regulation, water withdrawals, and channelization, induces water stress in freshwater ecosystems which might have implications on the water quality, biological structure, and functioning [4,8]. This approach at the SUDOE territory scale with more diversified and intense human stressors, explored for the first time, highlights the response of NR and NP to drought climates and related withdrawals and land uses.
Some drivers such as river hydrogeomorphology and the availability of reaction products (nitrate, carbon, and hydrogen) influence NR function. Several multidisciplinary works have shown that physical features like hydrogeomorphology are just regulation drivers of biological transformations, such as denitrification, which play the main role in nitrogen balance at the watershed scale [24,70]. Further studies are thus encouraged to continue to identify the main influencing factors of NNB.

The Human Impact on NR Function Determination
Some environmental factors seem to partly explain the NR simulated at the SUDOE scale. However, the results also showed that a factor integrating human stressors is necessary to evaluate NR functions at large scale since NR is influenced by land cover. This study found a correlation between NR and the IHA, an indicator integrating the level of anthropization. Human activities have a direct impact on nitrate-related ecological functions, particularly during low water periods [8]. Indeed, agriculture, cities, industries, and the management of water resources shape the landscapes and influence natural processes including nitrate production and removal [27,71,72]. NR variability caused by anthropogenic influence depends also on the soil-climate variability and resource management policies. These variabilities and policies greatly influence ecological functions [73][74][75].
NR functions are impacted by the availability of nitrate and organic carbon which are essential to activate denitrification processes even if an excess of nitrate can limit denitrification [76]. It could be interesting in further studies to focus on the nutrient range activating NR and on the potential inhibition limit of this concentration for NR [77]. The aquatic plant diversity and abundance will have a role in NR by plant and algae uptake [78]. Plant diversity and abundance will also be dependent on stakeholder policies [20,79], especially on the restoration of aquatic ecosystems [80]. Finally, the NR might be explained by environmental factors coupled with anthropization indicators, ecology diversity indicators, and nutrient concentrations. This study encourages further studies to focus on the model development to forecast NR functioning.

Setup of EF Indicators to Quantify Water Quality Regulation Ecosystem Services
Quantifying ecological functions make it possible to estimate indicators for ecosystem services assessment, such as the service of water quality regulation. The modeling approach allows an integrative view of these indicators with spatial and temporal trends inside the watershed. Analyzing simultaneously NP, NR, and NNB permits us to understand where and when the maximal potentials of nitrate production, total nitrate removal, and the final balance of these functions (NNB) occur in the watershed. The quantification of ecosystem services is relevant for public and private institutions to understand the importance of water treatment capacities provided by natural processes in parallel to the water treatment plans. This knowledge is essential for the regional planning of economic development with a sustainable approach by taking into account these ecosystem services and preserving them [81]. Moreover, the information about spatial and temporal indicators at the reach scale allows a global vision of the stream functioning, i.e., to identify the reaches that help purifying and the ones that cannot accept more pollution inputs. This overview of the different reach abilities is supplied as a functional support to orientate and help stakeholders in their management decisions such as the location of new implementation for water treatment structures. Water resource management, with this information in hand, could help a sustainable development that combines natural and human-made treatment possibilities, with respect to the EU-WFD policies. The spatial and temporal quantification of NNB capacity could also be included into the design of a new methodology to assess the ecological status of European rivers by combining regular and functional bioindicators for a better completion of the EU-WFD.
Furthermore, the NP, NR, and NNB indicators can be used to estimate the cost avoided for water treatment thanks to the free processes provided by the environment, i.e., the total amount that the taxpayer could not avoid paying in order to preserve the society standard of living and health in case of ecosystem service loss. A methodology based on SWAT outputs has recently been proposed to evaluate the ecosystem services of water quality regulation at the scale of the Ebro watershed in the North of Spain [82]. This methodology should now be extended to understand the spatial and temporal dynamics of ecosystem services at pluri-regional scales. For example, using the results of the present study would permit the quantification of water quality regulation services at the scale of the South-West of Europe. Based on this modeling completion, it could be possible to estimate the biophysical and economic values of this service at the regional scale, as economic arguments for the conservation of natural river services.
Finally, further studies should focus on human stressors and environmental drivers that will affect these NP and NR functions. In parallel, the practices and conservation of the river compartments that favor these service deliveries should be encouraged. Since it is well know that the nitrogen cycle is related to instream algae biomass, to hyporheic and river wetlands, and to the riparian forest health conditions, these compartments should be conserved or restored according to their real values. In this way, increasing wetlands and riparian forests will create a natural barrier effect on leaching and consequently will lower NP and increase NR. The restoration of this barrier is also favored by the conservation of the logjams that create islands inside the stream where biodiversity finds habitats and sediment flow is slowed down.

Conclusions
This paper assesses the water regulation functions related to nitrate at the pluriregional scale of the South-West of Europe, based on SWAT modeling. The quantification of nitrate-related ecological functions at the subsystem scale is validated based on the similarity between the modeled streamflow, nitrate load, and NNB outputs and those obtained with a higher resolution model (at the water body scale). Hot spots and hot moments of nitrate production and removal are identified in South-Western Europe. The hot spots of nitrate removal are located in the South of the territory, characterized by a warmer and dryer climate, whereas the hot spots of nitrate production are located in the main streams more impacted by human activities. Summer and spring are the seasons where the biological processes are the most active for nitrate production and removal. Finally, this study also highlights the existence of environmental and anthropic drivers of nitrate retention, such as streamflow and river slope, or the Index of Hydrologic Alteration. This study highlights the importance of considering other indicators such as anthropic indicators and biodiversity abundance indicators to explain the nitrate removal functions.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/w13212980/s1, Figure S1: PBIAS performance of nitrate net balance comparison between water body and subsystem models over Garonne watershed. The simulation has been performed from 2000 to 2010; Figure S2: Spatial repartition of three groups defined by NNB value simulated by the subsystem model over SUDOE territory between 2000 and 2010.