QUAL-NET, a high temporal-resolution eutrophication model for large hydrographic networks

To allow climate change impact assessment of water quality in river systems, the scientific community lacks efficient deterministic models able to simulate hydrological and biogeochemical processes in drainage networks at the regional scale, with high temporal resolution and water temperature explicitly determined. The model QUALityNETwork (QUAL-NET) was developed and tested on the Middle Loire River Corridor, a sub-catchment of the Loire River in France, prone to eutrophication. Hourly variations computed efficiently by the model helped disentangle the complex interactions existing between hydrological and biological processes across different timescales. Phosphorus (P) availability was the most constraining factor for phytoplankton development in the Loire River, but simulating bacterial dynamics in QUAL-NET surprisingly evidenced large amounts of organic matter recycled within the water column through the microbial loop, which delivered significant fluxes of available P and enhanced phytoplankton growth. This explained why severe blooms still occur in the Loire River despite large P input reductions since 1990. QUALNET could be used to study past evolutions or predict future trajectories under climate change and land use scenarios.


Introduction
River eutrophication has become a rising problem over the past decades, especially in India, Asia and South America, constituting a major risk for ecosystems and human health (e.g., Braga et al., 2000;Némery and Garnier, 2016;Yin et al., 2016;Dixit et al., 2017).Significant efforts to reduce non-point and point sources of nitrogen (N) and phosphorus (P) were made in Europe and North America, leading to a eutrophication decline in several large rivers (Hartmann et al., 2007;Friedrich and Pohlmann, 2009;Howden et al., 2010;Hardenbicker et al., 2014;Minaudo et al., 2015Minaudo et al., , 2016;;Powers et al., 2016;Poisvert et al., 2017).However, eutrophication crises are still occurring in many freshwater areas.
Previous studies often tried to assess which controlling factor of eutrophication prevails over the others and often contrasted nutrient availability with supposedly favorable physical conditions.Conflicting results in the literature did not solve this issue.In some rivers chlorophyll a concentration could directly be assessed confidently from P concentration (e.g., Basu and Pick, 1996;Dodds, 2006), whereas river flow conditions in other systems clearly constrained and determined the algal biomass (Biggs and Smith, 2002;Istvánovics et al., 2009).A few studies identified a combination of variables co-controlling phytoplankton blooms like the association of river flow conditions, water temperature and sunshine duration over the preceding days (Bowes et al., 2016), flow and light intensity (Hardenbicker et al., 2014), and flow, temperature and nutrient availability (Van Vliet and Zwolsman, 2008).If reducing P inputs has proved to be efficient to limit phytoplankton blooms in rivers, many recent studies show that both N and P availability must be considered as key elements to determine the trophic state of streams and rivers (Dodds and Smith, 2016;Paerl et al., 2016).Apart from nutrient availability, numerous other factors control phytoplankton composition and abundance in rivers, such as water residence time (directly linked to the river morphology, with the potential presence of flow veloc-ity dead zones), penetration of solar radiation into the water column (depth and turbidity), water temperature variations (hydrological and climate forcing), invertebrate grazing from endemic and invasive species, and self-shading effects by the phytoplankton colony itself (Reynolds and Descy, 1996;Reynolds, 2006;Abonyi et al., 2018).
Disentangling the relative influence of so many chemical, biological and physical factors on the river biogeochemistry can hardly be captured confidently through simple water quality monitoring and often requires the help of numerical modeling.Many deterministic water quality models at the catchment scale were developed and used initially to estimate nutrient source inputs into receiving waterbodies and support watershed stakeholders and decision makers to tackle eutrophication issues (Wellen et al., 2015).Only a limited number of models propose a mechanistic module simulating phytoplankton community dynamics and its impact on eutrophication.One can cite Riverstrahler (Billen et al., 1994;Garnier et al., 2002), ProSe (Even, 1995;Even et al., 1998;Flipo et al., 2004;Vilmin et al., 2015), Pegase (Deliège et al., 2009), QSIM (Kirchesch and Schöl, 1999;Schöl et al., 1999), WaterRAT (McIntyre and Wheater, 2004), QUAL2KW (Pelletier et al., 2006), WASP7 (Ambrose and Wool, 2009), QUASAR (Whitehead et al., 1997) or RWQM1 (Reichert et al., 2001).However, many of these models are only able to simulate river stretches and not the entire river network.The main reason is that very few models work at the catchment scale with subdaily timesteps (Wellen et al., 2015), mostly because program developers have to face long calculation times and usually compromise between a large spatial scale and high temporal and/or spatial resolution.The use of a high temporal resolution is, however, required to account for hydrological and biogeochemical processes occurring over short periods of time (e.g., storm events or subdaily phytoplankton growth variations).Water temperature is also a key factor for phytoplankton abundance and assemblage (Reynolds, 2006) which needs to be simulated at a high temporal frequency to assess the impact of potentially drier streams and warmer summers under climate change (Quiel et al., 2010).Developing methods appropriate to the regional scale is also required to account for instream processes in large rivers that control N, P and carbon (C) variations and constrain water quality in estuarine and coastal zones.Finally, models need to be appropriate for regional studies, i.e., the scale at which actions are undertaken by water body stakeholders and catchment managers.
The objectives of our study were twofold: firstly, develop a model able to simulate hydrological and biogeochemical processes in drainage networks at the regional scale (over 10 4 km 2 ), with hourly resolution and water temperature explicitly determined to allow potential climate change impact assessment; secondly, disentangle the different processes involved in eutrophication in a large river and identify their main drivers.To achieve this, the model QUALity-NETwork (QUAL-NET) was developed based on the integration of a biogeochemical model, Rive (Garnier et al., 2002), in a thermal model, T-NET (Beaufort et al., 2016).This new model was tested on a selected portion of the Loire River basin, the Middle Loire River Corridor, draining 43 × 10 3 km 2 , where the river main stem (270 km long) is prone to eutrophication in summer (Lair and Reyes-Marchant, 1997;Descy et al., 2011;Minaudo, 2015;Minaudo et al., 2015).

Study site
The Loire River (110 × 10 3 km 2 ) is the largest river flowing in France.The selected Middle Loire River Corridor (MLRC) is an intermediate sub-catchment located in the lowland section of the river main stem (Fig. 1).It separates the Upper Loire (a mountainous area where anthropogenic pressures greatly impact the river water quality but where eutrophication is only visible in lakes and reservoirs; Jugnia et al., 2004) from the Lower Loire River where the river main stem meets its major tributaries (Cher, Indre, Vienne and Maine rivers).The MLRC starts 450 km from the headwaters and runs over 270 km.From its entrance to its outlet (stations S1 to S2 in Fig. 1), the cumulated catchment area in the MLRC increases by only 26 %.This section of the river has a high eutrophication potential, combining most of the conditions favoring phytoplankton growth: high N and P concentrations (Minaudo et al., 2015), a low water level in summer (∼ 1 m), and a morphology with multiple channels and numerous islands slowing down flow velocity which increases the water travel time (Latapie et al., 2014).Chlorophyll a concentration was often over 250 µg L −1 in the 1980s, and many efforts have been conducted since 1990 to limit phosphorus point and non-point sources and counteract eutrophication: since 1990, phosphorus concentrations have been divided by 2.5 and phytoplankton blooms have declined 3-fold (Floury et al., 2012;Minaudo et al., 2015;Oudin et al., 2009).Even if phytoplankton in the Loire system is now clearly Plimited, algal blooms still occur (Abonyi et al., 2012), questioning phosphorus as a source and suggesting potential recycling processes.

Methods
The model QUALity NETwork (QUAL-NET) was developed based on a deterministic approach.It is the coupling between a thermal model, T-NET (Beaufort et al., 2016), and a biogeochemical model, Rive (Garnier et al., 2002).
Model T-NET is a physically based model able to estimate the water temperature in each reach of a large hydrographical network (10 5 km 2 ) with an hourly resolution and low errors, especially in the lowland area (Beaufort et al., 2015(Beaufort et al., , 2016;;Loicq et al., 2018).It has previously been developed specifically for the Loire River Basin (110 × 10 3 km 2 and over 50 × 10 3 river reaches from headwaters to the estuary).The temperature in the river network is computed as follows: (i) resolution of the heat budget in a given reach and estimation of the equilibrium temperature (Bustillo et al., 2014); (ii) longitudinal propagation downstream of the thermal signal according to the estimated water velocity throughout the river reach; (iii) discharge-weighted mix of the thermal signal when two or more streams meet in one node.
The model Rive is a mechanistic model describing many of the biogeochemical interactions that occur in the river between the water column and the benthos.It simulates the dynamic of dissolved and particulate organic matter, nutrients (N, P, Si), dissolved oxygen, the phytoplankton biomass (three algae groups: green algae, diatoms and cyanobacteria), zooplankton and bacteria.Equations from model Aquaphy (Lancelot et al., 1991) were used to describe primary producer variations.The model Rive is the biogeochemical module of the Riverstrahler (Billen et al., 1994) and ProSe (Even et al., 1998) models.Riverstrahler was largely used in past studies to simulate with a 10-day time step biogeochemical variations in large lowland eutrophic rivers under varying climate conditions, e.g., the Seine Basin, the Danube River, the Red River in Vietnam, and over long periods of time (Garnier et al., 1995(Garnier et al., , 2002(Garnier et al., , 2005;;Billen and Garnier, 2000;Billen et al., 2001;Quynh et al., 2010).Equations and variables included in the model are extensively described in Billen et al. (1994) and Garnier et al. (2002).Both the water and the benthic components were considered, including chemical and physical exchanges in-between these two components, according to the Billen et al. (2014) formulation.Equations in this formulation provided estimates of nitrogen, phosphorus, silica and dissolved oxygen fluxes across the water-sediment interface.The sediment layer was split into two sub-layers.The one at the bottom was considered compact and not erodible, the other one could potentially be resuspended.Nutrient fluxes between these two sediment layers were also considered in our model.For each river reach and at each time step, the model estimated quantities of eroded particles or those settled on the riverbed.Particles were considered both inorganic and organic with three levels of lability.Resuspension could potentially fuel the water column with soluble reactive phosphorus via desorption processes from suspended matter.
The temporal resolution of QUAL-NET was hourly.QUAL-NET was coded in C++ language and allowed parallel computing, i.e., the simultaneous use of several processors in order to reduce computation time as much as possible.Simulating hourly biogeochemical evolutions of 3361 stream segments over a 3-year period took nearly 4 h on a two-processor platform (Intel(R) Xeon(R) CPU E5-2670 0 @ 2.60 GHz) with 16 cores (64 Go, DDR3 = 1600 MHz).

Data inputs and main spatialization choices
Hydrological, geomorphological and meteorological forcing variables were determined and used on the basis of T-NET model implementation on the Loire Basin (Fig. 2).Thus, a more detailed description is available in Beaufort et al. (2016), except for nutrient sources forcing variables.

Meteorological variables
Hourly meteorological variables were taken from the SAFRAN atmospheric reanalysis (Quintana-Segui et al., 2008) T-NET thermal model (Beaufort et al., 2016) Δx = segment ≈ 1 km Δt = 1 h RIVE biogeochemical model (Billen et al., 1994) Δx = segment ≈ 1 km Δt = 1 h Water temperature C, N, P, Si, O 2 , TSS, phytoplankton, zooplankton, bacteria rological variables were used to compute the hydrological model (see below) for both thermal and biogeochemical modules.Air temperature, specific humidity, wind velocity and atmospheric radiation were used to compute water temperature.Most biogeochemical variables were water temperature dependent, and phytoplankton photosynthesis processes were directly linked to atmospheric radiation variations.

Hydrology
Daily mean discharge and groundwater flows were simulated by the semi-distributed hydrological model EROS (Thiéry and Moutzopoulos, 1995) at the outlet of 17 subwatersheds.Within each of these sub-watersheds, flows were redistributed into the hydrographic network according to the corresponding drainage area of each river reach.This approach proved its efficiency and reliability at the regional scale in the Loire Basin (Beaufort et al., 2016).Discharge and groundwater flows were considered constant over 24 h even if the water quality model output was hourly.

Geomorphology
The hydrographical network was determined from the Carthage ® database (Carthage, 2012), after transforming multiple channels into single channels.In the MLRC, we counted 3361 reaches, every one of them being defined as the river section between two confluences.Slopes for each reach were assessed based on a 25 m resolution digital terrain model (ALTI, 2012).Transversal morphology in streams were assumed to be rectangular, while depth and width were assessed on a daily time step but differently for the Loire River main stem and other streams: (i) depths in the Loire River main stem reaches were assessed based on field measurements conducted during both low-and highflow periods (Latapie, 2011;Latapie et al., 2014) combined with Manning-Strickler formulation (Strickler coefficient was calibrated numerically); (ii) in all other rivers and streams, where no field measurements were available, depth and width were assessed on a daily basis on the ES-TIMKART application (Lamouroux et al., 2010) which uses stream slope, watershed area, and daily and interannual discharge to estimate streams morphology.

Non-point sources
Non-point sources of nutrients and exports of total suspended solid concentrations (TSSs) were defined based on land use (European Corine Land Cover dataset, 2006), climate characteristics, lithology (LITHO ® , 2008) and previous observations conducted on 108 streams located in the Loire headwaters, upstream of any potential point sources (Blan-chard, 2007).Overall, land use categories were grouped into seven large categories (urban, arable land, cultivated land, prairie, forest, wetland, other types) and associated with a corresponding non-point source concentration for the following variables: nitrate, ammonium, total inorganic phosphorus, biogenic silica, dissolved and particulate organic carbon for three different biodegradability classes, total suspended solids, and fecal matter.The MLRC Basin was divided into 479 small sub-catchments (the average was 27 km 2 ), and diffuse source concentrations were applied homogenously for all streams located in a given sub-catchment as a combination of concentrations originating from all the different land use types.Land use was considered constant over time, leading to constant nutrient concentrations for non-point sources.Thus, it was hypothesized that the hydrological variability alone could be responsible for temporal variations in nonpoint nutrient fluxes.

Point sources
Industrial and domestic point sources of nutrients and TSS fluxes originated from Loire Basin water authorities (AELB) surveys conducted in 2010.In the MLRC Basin, 641 wastewater treatment plants (WWTPs) were recorded.Datasets provided total organic carbon, total nitrogen and total phosphorus fluxes.It was estimated in 2010 that point sources in the Middle Loire sub-catchment (our study) represented 322 kg P day −1 and 1.9 t N day −1 .These fluxes were divided into the different chemical forms for C, N and P, according to Servais and Billen (2007), depending on the type of point sources and the characteristics of WWTP.Fluxes of point sources were considered constant over time.

Upstream boundary in the Loire River and validation dataset at catchment outlet
A daily survey was conducted at S1 (Saint-Satur-sur-Loire) and S2 (Cinq-Mars-la-Pile) in the Loire River during the period August 2011-July 2014 (Minaudo, 2015).Data collected at S1 were used as data input for the model, and data at S2 were used for both calibration and model performance assessment.Samples were collected every day from a bridge using the same procedure at each station.TSSs were measured every day.The following parameters were analyzed on a 3-day frequency basis: dissolved and particulate organic carbon (DOC and POC), total and soluble reactive phosphorus (TP and SRP), nitrate (NO − 3 ), dissolved silica (Si), and chlorophyll a concentrations.Filtrations were immediately made on-site using 0.45 µm cellulose acetate membrane filters for chemical parameters and a 0.70 µm glass filter (Whatman GFF) previously burned at 500 • C for 6 h for chlorophyll a and POC analysis.Total suspended solid concentrations were determined by filtration of a precise volume of each water sample through pre-weighed filters and by drying them at 105 • C.After filtration, water samples and filters were stored at −80 • C in polypropylene tubes after acidification of aliquots for NO − 3 , SRP and DOC analysis.Tubes and filters were unfrozen on the day of the analysis.DOC concentrations were measured with a carbon analyzer (Shimadzu TOC-V CSH/CSN).The NO − 3 concentration was determined by ionic chromatography.Phosphorus was measured by colorimetry after solid digestion in the case of TP analysis (potassium-persulfate digestion).Dissolved silica (Si) was measured by colorimetry.For POC analyses, filters were first treated with HCl 2N to remove carbonates, dried at 60 • C for 24 h and then measured with a C / S analyzer (LECO C-S 200).Chlorophyll a was measured by fluorimetry at a wavelength > 665 nm after an excitation step between 340 and 550 nm.Chlorophyll a concentrations were expressed in mg C L −1 considering a C : Chl a ratio of 31, according to Minaudo et al. (2016), and constituted the variable hereafter named "PHY".

Computation steps in the model
Computation in the model was based on a network topology: each reach in the hydrographic network corresponded to the stream segment between two confluences.Each reach had an upstream (or upper) and a downstream (or lower) node (Fig. 1).Except for Strahler first-order streams in the headwaters, upper nodes were always connected to two lower nodes.

Initialization at upper node and boundary conditions
All variables were initialized at the upper node of Strahler first-order streams.Water component variables were initialized according to non-point sources estimated for hillslope catchments located upstream of the upper nodes.Variables in the sediment component were initialized homogenously everywhere in the stream network, assuming that the model should quickly reach its equilibrium based on the interactions with variables from the water component.The upstream boundary in the Loire River (S1) was determined based on the daily survey conducted at S1 (see Sect. 3.1.6).

Propagation downstream
All variables computed at one reach in the water component were transferred downstream according to travel time estimated from discharge and stream morphology.Variables from the benthic component interacted with the water component but were not transferred downstream.For a given time step, a given reach was discretized depending on the estimated travel time.If travel time was shorter than 1 h, the reach was not segmented: thermal and biogeochemical equations were solved at the downstream node.If travel time exceeded 1 h, the reach was segmented into as many subsegments as needed to get 1 h travel time subsegments.with a 15 min sub time step to avoid potential numerical resolution drifts.When two streams met, biogeochemical signals were mixed with respect to stream discharge.This determined the values for the next downstream upper node.Because the exact location of potential WWTP input within a segment was not always known, it was assumed that the location of point sources occurred at the downstream node.

Calibration step
The thermal model was fully deterministic and no calibration step was needed.Even if Rive was built as a universal representation of the mechanisms occurring in rivers, some processes were based on empirical relationships.Nearly 150 coefficients were counted overall (Fig. 3), the majority of them were used to describe bacteria and phytoplankton dynamics depending on light intensity, water temperature, and nutrient availability.Most coefficients are currently accepted as universal constants, but several studies pointed out that hydrosedimentary and P sorption-desorption processes needed experimental or numerical calibration (Vilmin et al., 2015), especially because processes involved highly impacted performances on phytoplankton and water quality predictions (Aissa-Grouz, 2015).The phosphorus dynamic in the water compartment was based on the Langmuir equilibrium concept (Limousin et al., 2007), a description largely found in the literature for water quality models (e.g., Chao et al., 2010;Rossi et al., 2012;Vilmin et al., 2015).Very different values for P sorption-desorption coefficients were found experimentally or numerically in the literature, with up to 5 orders of magnitude differences from one study to another (Vilmin et al., 2015).No specific laboratory experiments were conducted in the Loire River, leading us to deploy numerical calibration methods to calibrate TSS and SRP dynamics.Because SRP computation relies on TSS dynamic, the first variable calibrated was TSS.Calibration was conducted manually by changing the values of the different coefficients to be calibrated over a range of values found in the literature.In total, five coefficients were calibrated.The best set of coefficients was selected when results minimized root mean square errors (RMSEs) of the calibrated variable.Among the recorded time series (1 August 2011 to 31 July 2014), the period selected for calibration was the first year, i.e., 1 August 2011 to 31 July 2012, and the remaining data served for validation.

Calibration of TSS dynamic
Total suspended solid concentration increments (dTSS) were computed based on a simple difference between eroded matter from the riverbed (eros TSS ) and settled particles (sedim TSS ), as described in Eqs.(1-4).Erosion was defined as a power law function of flow velocity (Eqs. 2 and 3).
where Vs TSS was the sedimentation velocity, depth(t) was the water depth at time t, Cap TSS was the erosion capacity depending on coefficients Veli 0 , Veli 1 and flow velocity V (t), SED was the height of the layer of sediments potentially erodible, and SED 0 was the layer of sediments set during initialization step.Thus, TSS concentration depended on the coefficients Veli 0 , Veli 1 and Vs TSS .These three coefficients were calibrated.

Calibration of P dynamic
SRP concentration was estimated based on sorptiondesorption equations originating from the Langmuir equilibrium displayed by Eqs. ( 5) and ( 6).This formulation requires the maximal sorption capacity of P onto suspended solids (Pac, in mg P g −1 ) and a half-saturation constant (Kpads, in mg P L −1 ) that needed to be defined.
where TIP corresponded to total inorganic phosphorus concentration at time step t and Kpads and Pac were the two parameters needing to be calibrated.

Model performance criteria for validation
To estimate model performances and define criteria for model validation, bias and standard deviation errors were used, following Eqs.( 7) and ( 8): where SD was the standard deviation and n the total number of observations.These metrics were calculated for each variable observed at S2 over the entire period of validation (1 August 2012 to 31 July 2014), and were also computed seasonally: "summer" corresponded to the bloom season from April to October; "winter" corresponded to the remaining part of the year.

Lagrangian representation and fluxes budgets
In addition to more common ways of presenting results longitudinally, we proposed two other graphical representations of transfers and biogeochemical transformations from S1 to S2 along the Loire River.One representation consisted in following the same water body transferred from S1 to S2, i.e., a Lagrangian representation.Lagrangian profiles were estimated from the matrix of travel time computed for each reach and at each time step from measured discharge and river morphology estimates (see Sect. 3.1.3).This representation was both spatial and temporal since it displayed longitudinal variations according to travel time going downstream.It was used for two typical situations: one in winter (starting on 9 February 2013 during a high-flow period) and another one during a phytoplankton bloom (starting on 10 July 2012).Additionally, average seasonal flux budgets of all the main processes (inputs and outputs) simulated between S1 and S2 over winter or summer periods were computed for a selection of variables (TSS, NO − 3 , total inorganic P, Si, PHY, DOC, POC and O 2 ).In those graphs, arrow widths were proportional to the corresponding calculated flux, allowing the comparison between the two different seasons.

Constant flow and constant water temperature simulations
The sensitivity of phytoplankton variations to constant flow conditions in the Loire River for both low-flow and high-flow conditions (200 and 1000 m 3 s −1 , respectively) was assessed.
A similar approach was tested with constant water temperature (13.7 • C, i.e., the average temperature).Lagrangian profiles during a phytoplankton bloom (starting on 10 July 2012) of these simulations can be found as Fig. S1 in the Supplement.

Calibration step
The best set of coefficients that minimized errors over the period are displayed in Table 1.RMSE on calibrated variables were 15 mg L −1 for TSS and 14 µg P L −1 for SRP.The selected values for TSS coefficients largely differed from other values found in the literature, justifying the need for this calibration step.Compared to the Seine River, it appeared necessary to increase the erosion capacity (Veli 1 ) but also to reduce considerably suspended solids sedimentation rates (Vs TSS ), which resulted on an increased sediment reactivity within the Loire system.Values calibrated for P sorption processes were close to the values found experimentally in the neighboring Seine Basin (Aissa-Grouz, 2015

Model performances at station S2
Over the study period, discharge variations at S2 presented highly seasonal variations (Fig. 4): Q ranged between 60 and 150 m 3 s −1 in summer low flows and peaked over 1200 m 3 s −1 in winter high flows.Water temperature simulated by T-NET was highly seasonal and fluctuated between 0 and 30 • C. In summer, the amplitude of water temperature diel cycles ranged between 0.2 and 1.5 • C. Phytoplankton concentrations presented three clearly delimited bloom events between March and September.The maximum observed each year was 60 to 70 µg chl a L −1 corresponding to 1.6 and 1.9 mg C L −1 .Observed TSS concentrations were correlated with discharge and ranged between nearly 0 in summer and 150 mg L −1 during high flows.Nitrate concentrations presented a clear seasonal signal, fluctuating between ∼ 1.5 mg N L −1 in summer to ∼ 3.5 mg N L −1 in winter.Dissolved silica concentrations ranged between nearly 0 and 8 mg Si L −1 .Concentrations always peaked in winter during high flows and dropped in spring, concomitantly with the start of phytoplankton activity.Soluble reactive P con-centrations presented a clear seasonal cycle, with very low concentrations reached during summer (< 10 µg P L −1 ) and relatively high concentrations in winter (∼ 60 µg P L −1 ).Particulate organic carbon concentrations ranged between 0.4 to 5 mg C L −1 , with strong correlations between POC and TSS in winter and between POC and phytoplankton biomass during algae blooms (Minaudo et al., 2016).Dissolved organic carbon concentrations ranged between 4 and 10 mg C L −1 .The highest concentrations were observed during high-flow periods, but no clear seasonal variations could be deciphered.QUAL-NET provided reasonable estimations for the main variables (see Table 2 for bias and standard deviation errors).Seasonal variations were correctly simulated for all variables.At the scale of the storm event, a few events were observed with the daily survey but were not represented by the model, especially for several storm events that occurred under low-flow conditions.A phytoplankton bloom event at the end of summer 2012 was simulated but did not correspond to our observations.Dissolved oxygen was not measured, but concentrations simulated by QUAL-NET presented a clear seasonal cycle, with ∼ 12 mg O 2 L −1 estimated in winter, and Figure 4. Results at station S2 after calibration for the main variables in the model: discharge (Q), phytoplankton (PHY), nitrate (NO − 3 ), dissolved silica (Si), soluble reactive phosphorus (SRP), water temperature, total suspended solids (TSSs), particulate organic carbon (POC), dissolved organic carbon (DOC), dissolved oxygen (O 2 ).Last row zooms in on SRP and O 2 concentrations in July 2013 to show simulated diel fluctuations.6 to 9 mg O 2 L −1 in summer.During phytoplankton blooms, the model provided interesting diel fluctuations for PHY, SRP and O 2 concentrations.For instance, SRP concentration fluctuated between 0 and 15 µg P L −1 and O 2 concentrations presented a minimum at midnight and a maximum at noon.Unfortunately, the reliability of these variations could not be verified with our measurements.Performances appeared similar between seasons (Table 2) with approximately the same range of errors in winter or summer, except for dissolved sil-ica whose simulated concentrations in winter were subject to higher imprecisions (2.1 against 1.3 mg Si L −1 in summer) and for PHY with lower absolute errors in winter (a period with very low PHY concentrations).

Lagrangian views of winter versus summer dynamics
The Lagrangian views of the evolution of the different biogeochemical species highlighted different hydro- 3 , SRP, Si, PHY and O 2 during two selected events.For the summer event, P input from mineralization processes, P uptake from phytoplankton, phytoplankton growth rate (availability of intracellular carbon and nutrients), sedimentation and mortality rates are also displayed on the right axis.biogeochemical functioning depending on the season, (Fig. 5).The selected winter event corresponded to a highflow period: Q at S1 was 940 m 3 s −1 and increased to 1110 m 3 s −1 by the time the water arrived at S2.It took almost 2 days for the water to travel between S1 and S2 (∼ 250 km).Most elements were simply transferred downstream, with no significant transformation or alteration between S1 and S2.The concentration of TSS presented a decreasing evolution from 33 mg L −1 at S1 to 25 mg L −1 at S2. Nitrate concentration slightly increased from 2.8 to 3.1 mg N L −1 (+11 %), and so did SRP (+40 %).Dissolved silica concentration decreased (−12 %).Phytoplank-ton activity remained very low and declined steadily (5 to 2 µg chl a L −1 ).Dissolved oxygen slightly increased (+8 %).
During the selected summer event, discharge was much lower: Q was 330 m 3 s −1 when the water left S1 on 10 July 2012 and increased to 340 m 3 s −1 when the water reached S2.
The model estimated that it took nearly 3 days for the water to cover the distance from S1 to S2, and all biogeochemical variables were largely modified when the water moved downstream.Two steps were identified: -The first 2.5 days, total phytoplankton concentration increased from 0.5 to 1.7 mg C L −1 .Simultaneously, SRP was dramatically depleted from 50 to nearly 0 µg P L −1 .Nitrate, silica and oxygen concentrations slightly de- creased (∼ −10 %).The amount of P released from organic matter mineralization remained limited but reached a first peak concomitantly with a large P uptake from the phytoplankton colony.Phytoplankton mortality rates kept increasing and peaked when growth rate reached its maximum (0.15 mg C L −1 h −1 ), i.e., when travel time from S1 was 2.3 days.
-Then, during the next 24 h (the time needed for the water to reach S2), phytoplankton concentration started to decrease (−15 %); SRP remained very low under 5 µg P L −1 and presented a diurnal fluctuation with a minimum reached during the afternoon and rising concentrations when the water arrived at S2 by night.During this phase, organic matter mineralization as a source of inorganic P increased substantially from 2 to 13 µg P L −1 h −1 , and phytoplankton growth rates first dropped from 0.15 to near 0 mg C L −1 h −1 and then rose again to 0.1 mg C L −1 h −1 when SRP input from mineralization counteracted phytoplankton uptake.

Storm event disturbance during a phytoplankton bloom
A storm event occurred in August 2013, during a phytoplankton bloom.Over 5 days (9 to 14 August), discharge at S2 increased from 200 to 406 m 3 s −1 and then declined to reach 230 m 3 s −1 on 19 August.This largely disturbed TSS, SRP and PHY dynamics (Fig. 6).This storm event caused a suspended solids peak which propagated over the entire river stretch.TSS concentration peak amplitude decreased from 120 to 50 mg L −1 when the water moved downstream from S1 to S2, and the peak width widened.At the beginning of the event, the SRP concentration profile was showing a complete P depletion starting approximately 80 km downstream S1.This P limitation threshold progressively moved further downstream when the storm event hit.SRP slightly

Si (t day -1 )
Figure 7. Average winter and summer budgets between S1 and S2 for TSS, nitrate, inorganic P and dissolved silica.All arrow widths are proportional to calculated fluxes, allowing the visual comparison between winter and summer periods.increased at S2, but concentrations remained very low.When the discharge peak hit S2 (14 August), SRP concentrations presented a steady longitudinal decline from 50 µg P L −1 down to nearly 0. Before the storm event, phytoplankton concentrations showed a limited longitudinal increase, from 0.5 to 1.2 mg C L −1 , but when the discharge peak event occurred, PHY concentrations decreased in the upper part of the MLRC but clearly increased in the lower part suggesting that phytoplankton was flushed away by the storm event.PHY concentrations during discharge recession presented an increasing longitudinal profile from 0.1 to 1.1 mg C L −1 and began to increase again everywhere along S1 to S2 when hydrological conditions stabilized.

Fluxes, transfers and transformations in the Middle Loire River Corridor
Proportions of the different contributions or biogeochemical transformations largely depended on the season (Figs. 7 and 8).In winter, most of the biogeochemical species entering the MLRC at S1 were transferred downstream, with nonsignificant interactions with the biological component.Suspended solids and particulate P showed an almost balanced budget between erosion and sedimentation processes.Lateral contribution between S1 and S2 remained small compared to the upstream flux at S1, except for nitrate because tributaries and lateral non-point source inputs contributed to 25 % of the total NO − 3 flux at S2. Reaeration of the water body represented a significant portion of the dissolved oxygen budget at S2 (14 %).
In summer low flows, the biological component largely modified the river biogeochemistry in the studied sector.Nitrate fluxes were 15% higher at S2 (38 t N day −1 ) than at S1 (28 t N day −1 ) despite N uptake by phytoplankton (3.2 t N day −1 ∼ 11 % S1 flux) and a moderate contribution from the lateral streams (12 t N day −1 ).Diffuse sources in the tributaries contributed to 94 % of total lateral inputs.Inorganic phosphorus loads were by 3 between S1 and S2 (from 1 to 0.3 t P day −1 ) due to phytoplankton and bacteria uptakes (respectively, 2.6 and 0.4 t P day −1 ).P recycling from organic matter mineralization (phytoplankton dead cells) supplied 1.3 t P day −1 , i.e., more available phosphorus than both upstream and lateral P inputs.Inorganic P inputs from WWTPs within the MLRC subbasin represented less than a third of P load in the Loire at S1 (0.3 t P day −1 compared to 1 t P day −1 ).Particulate inorganic P constituted a very small amount of total inorganic P, most of it was balanced between erosion and sedimentation processes.The riverbed acted like a source of inorganic P (299 kg P day −1 ).Dissolved silica fluxes were slightly affected by phytoplankton activity: 20 % of the flux at S1 was assimilated by diatoms.Lateral stream contributions represented 13 % of the flux quantified at S2. Phytoplankton increased 4-fold between S1 and S2 during summer blooms (Fig. 8), from 4.3 to 17.1 kg C day −1 , but this calculation only took into account the surviving cells when the water body reached S2.A larger proportion of phytoplankton actually grew but part of it decayed and was eventually recycled: the model estimated that 50 % of green algae and 25 % of diatoms that grew between S1 and S2 decayed.Additionally, approximately 25 % diatoms were deposited on the riverbed.The lateral contributions by the Loire river tributaries were not significant (only 0.1 kg C day −1 ), indicating that phytoplankton grew only within the river main stem.Approximately 100 t of organic C entered the system at S1 every day in low-flow periods (see also Minaudo et al., 2016).Approximately 80 % of it was dissolved, the rest was particulate.It was estimated that 16 t C day −1 in summer of DOC was bioavailable and consumed by heterotrophic bacteria.Part of this organic matter was eventually mineralized, depending on oxygen content.The dissolved oxygen budget was balanced between S1 and S2 (respectively, 192 and 208 t O 2 day −1 ), with oxygen inputs from primary producers (phytoplankton, 136 t O 2 day −1 ) similar to oxygen depletion by bacteria and zooplankton respiration processes (137 t O 2 day −1 ).

Discussion
Interannual, annual and seasonal variations in the main water quality variables simulated by QUAL-NET corresponded to the observations, proving the efficiency of the model at both transferring the different biogeochemical species and also modeling the main processes instream.term (but highly impacting) events such as a storm during a phytoplankton bloom.These performances were considered good enough to allow us investigate confidently the different processes and their drivers.This is highlighted in Sect.5.1 and 5.2.Additionally, QUAL-NET was subject to several weaknesses, and potential improvements will be made; this is detailed in Sect.5.3 and 5.4.
5.1 Drivers of eutrophication in the Middle Loire River Corridor

Biological versus hydrological control of the river biogeochemistry
The model showed that the Loire River biogeochemistry was the result of complex interactions between nutrient availability and hydrological variations.In winter, the MLRC was mainly controlled by hydrological processes, and nutrients were simply transferred downstream, with no noticeable control by biological processes.Under low-flow conditions and warmer water temperature, C, N, P and oxygen dynamics were dominated by biological processes.Stream algae were clearly P-limited and never reached N or Si limitations, supporting previous studies (Descy et al., 2011;Minaudo et al., 2015).In the MLRC, lateral inputs during summer were not significant compared to the magnitude of fluxes within the Loire River main stem.The highest phytoplankton concentration was not necessarily observed at the catchment outlet: during phytoplankton blooms, P was often depleted before the water could reach S2, and when this occurred, lower phytoplankton growth and higher mortality rates started to cause a decline in phytoplankton concentration.As soon as phytoplankton started to be P-limited, and bacterial activity caused the decrease in oxygen concentration (Li et al., 2014).When a storm event entered the Middle Loire system, the phytoplankton colony was flushed downstream, and, as long as physical conditions for phytoplankton growth remained degraded (shorter transit time, increased turbidity), available P was not totally assimilated.Thus, the SRP concentration increase in the lower section during such an event was rather the consequence of lower phytoplankton activity than increased inputs.Storm events in summer simply move the available P-exhaustion point further downstream.

P recycling within the MLRC
In summer, most of the inorganic P entering the MLRC was assimilated by phytoplankton and bacteria biomasses.However, mineralization of organic matter constituted a significant source of bioavailable P. The model estimated that P releases from mineralization was equivalent to all fluxes entering the MLRC (point and non-point sources), i.e., ∼ 1.3 t P day −1 .Besides, the phytoplankton concentration peak in-between S1 and S2 corresponded to inorganic P exhaustion.This caused a 15 % decrease in PHY con-centration when the water moved further downstream.Thus, SRP was most of the time fully assimilated by phytoplankton in summer, but phytoplankton was also subject to mortality and could partly be recycled to eventually constitute an autochthonous source of available P. Remineralization of autochthonous labile organic particulate P, known as part of the "microbial loop", is described in the literature of phytoplankton ecology (Reynolds, 2006;Li et al., 2014) and mostly identified in lakes, reservoirs or estuarine systems (Jossette et al., 1999;James and Larson, 2008;Song and Burgin, 2017) and rarely in rivers (Descy et al., 2002;Withers and Jarvie, 2008).On the one hand, bacteria compete with phytoplankton for SRP availability, and on the other hand, bacterial mineralization recycles P and supports phytoplankton growth.These observations confirm the necessity of considering bacterial activity as a major driver of carbon cycling in large eutrophic rivers.

High temporal resolution is needed in water quality models
The high temporal resolution in QUAL-NET enabled us to disentangle the interactions between hydrological and biogeochemical processes when a storm event occurred in summer low flow.Besides, the model identified diel fluctuations of O 2 or SRP (daily variations oscillated between 0 and 15 µg P L −1 during phytoplankton blooms).Diel fluctuations of O 2 were often observed and described in previous studies, directly linked to primary producers' activity (Moatar et al., 2001;Wade et al., 2012;Rode et al., 2016).Sub-daily fluctuations of inorganic phosphorus are rarely observed, but this is due to limited measurements of high-frequency variations in P concentration.Similar diel fluctuations were found in some other lowland eutrophic rivers, but these cycles were mostly explained as a balance between P contributions from direct sources and non-point sources (Wade et al., 2012).In the case of the Loire River, QUAL-NET simulates these diel fluctuations as likely being the result of complex interplay between biological uptake, P mineralization instream, lateral inputs and diffusion from the benthos.During the night, phytoplankton growth was zero, while lateral contributions from point and non-point sources still occurred.Additionally, P kept being diffused from the benthic compartment, resulting in an increased SRP concentration in the water column.After sunrise, as soon as the biological compartment started to assimilate more P than the amount of P originating from the different P sources, the SRP concentration decreased again.These subtle variations, revealed by the model, could not be seen based on the daily-scale survey and need to be confirmed with higher-frequency sampling measurements.

Sensitivity to phosphorus sorption-desorption representation
During the calibration step, QUAL-NET showed a high sensitivity to the formulation of phosphorus sorption-desorption processes.Compared to other studies using the same formulation, the optimized values found manually for our study appeared relatively close to those determined experimentally in the Seine River (Table 1).However, the large variability in the results when one of these two coefficients was modified challenged the use of the model: if modifications are conducted on the model (in terms of data inputs and/or processes), these coefficients should be recalibrated.This appears as an important weakness in the model until an experimental survey is deployed to assess spatial and temporal heterogeneities of P sorption-desorption characteristics in the Loire River.
5.4 Issues and potential improvements for model QUAL-NET

Conflicting time steps between forcing variables and output resolution
The use of a high temporal resolution in QUAL-NET proved its usefulness to model processes that occur over short temporal scales.However, the only forcing variables with such a fine resolution were the meteorological variables, allowing us to compute hourly water temperature and light availability in the water column.Flows, and therefore water depth and velocity, were measured on a daily basis, and spatial discretization for discharge was based on catchments that were on average 27 km 2 .Flows within each of these 17 catchments were redistributed into the hydrographic network according to the corresponding drainage area of each river reach, assuming a simultaneous temporal dynamic within each of the 17 catchments.This could provoke conflicting signal propagation through the hydro-system during storm events.A semi-distributed hydrological model could address some potential propagation issues during storm events, even if the output frequency remains daily because of the lack of discharge estimation on a sub-daily basis.Nutrients fluxes discharged from point sources were considered constant through time.Wastewater treatment plants efficiency in treating sewage can be seasonal (biological processes, variation in population in tourist areas) and sometimes highly impacted by storm events.Therefore, we urge local and national water basins authorities to provide at least monthly concentrations and fluxes for the different wastewater treatment plants, especially for plants treating sewage from the biggest cities.
Non-point source concentrations were constant through time.Therefore, it was assumed that only hydrological variability could modify non-point fluxes.This representation proved its reliability with a 10-day time step (Garnier et al., 2002) but misses many processes occurring at least at the seasonal scale, such as for instance nutrient retention by riparian vegetation during spring and summer (Peterjohn and Correll, 1984;Pinay et al., 1993), denitrification increases during warmer conditions, and peaks of nutrient concentrations during soil-rewetting events and when groundwater connects with streams (Dupas et al., 2015a, b).QUAL-NET proved to be efficient to model in-stream processes and would certainly benefit if coupled with land use models that predict more reliably nutrient non-point inputs such as SWAT (Douglas-Mankin et al., 2010) or HSPF (Fonseca et al., 2014).This would allow us to model the biogeochemical variations for the whole drainage system, instead of forcing the system with daily-scale measurements at S1.To upscale the model to the entire Loire Basin, the influence of lakes and reservoirs have to be considered since they largely modify nutrient transfers downstream.This raises another issue, because the connection between streams/rivers with lakes/reservoirs is hardly considered in water quality models at the catchment scale.

New eco-hydrological issues
In eutrophic rivers, several recent studies clearly showed the increasing concern with Asian Corbicula spp.clams that have invaded river networks in South and North America over the past decades and later in Europe (Cohen et al., 1984;Phelps, 1994;Cataldo and Boltovskoy, 1998;Pigneur et al., 2014).This clam plays a significant role in the dynamic of phytoplankton for several rivers in Europe.Pigneur et al. (2014) estimated for instance that Corbicula was responsible for a 70 % decrease in the phytoplankton biomass of the Meuse River.The main challenge with Corbicula is the lack of datasets, both spatially and temporally.In the Loire River, Descy et al. (2011) determined that a population density of 2.5 to 10 g C m −2 was needed to explain the phytoplankton variations, but clam density was then uniformly distributed depending on the river.This has not been tested in QUAL-NET yet, since very few surveys have been conducted, and spatial distributions of Corbicula spp.population remain unknown.In addition, aquatic fixed vegetation is able to extract nutrients from the sediment and might keep growing even if phytoplankton has reached its phosphorus limitation in the water column (Carignan and Kalff, 1980;Hood, 2012).Despite low P availability, macrophytes might keep growing, especially under a high P legacy in the riverbed sediments.We lack data about macrophytes in the Loire River, but a few unpublished observations in the MLRC presented very significant densities of Ranunculus fluitans, Myriophyllum spicatum and Elodea nuttallii (Michel Chantereau, personal communication, 2015).Their impact on the Loire River's biogeochemistry is likely significant; further developments in the model QUAL-NET should consider this biological compartment, and macrophyte biomass within the MLRC needs to be surveyed properly.The deterministic modeling approach we developed helped disentangle the interactions between hydrological and biological processes in the Loire River.Results from QUAL-NET fitted the available daily observations, and the main driving processes could be identified.The Middle Loire River Corridor functions as a biogeochemical reactor in summer during the low water period.The system clearly reaches a P limitation, and our model indicates that internal loadings of P due to bacterial mineralization enhance phytoplankton blooms.The use of a high temporal resolution enabled us to study the impact of a storm event during a phytoplankton bloom, and identified large diel fluctuations for C, P and O 2 , but these variations still need to be compared to highfrequency in situ measurements.QUAL-NET simulated realistic sub-daily variations from low-frequency forcing variables and could be applied at a larger scale (e.g., the entire Loire Basin, 110 × 10 3 km 2 ).It could be used to study past evolutions using a low-frequency dataset as data input or predict future evolutions under climate change and land use scenarios.
Data availability.The data used in this work are part of a wider database collected within the project "Eutrophication Trends" and cannot be publicly accessed for the moment.The global database will be properly deposited in a reliable public data repository as soon as it is fully analyzed and validated by all contributors.For more information, the reader is kindly asked to send requests to camille.minaudo@univ-tours.fr or florentina.moatar@univ-tours.fr.
Author contributions.CM and FC designed the model structure.
CM and YJ developed the model code.FM, NG and CM designed and conducted the daily sampling and the chemical analysis.CM prepared the manuscript with contributions from all co-authors.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Study area, the Middle Loire River Corridor sub-catchment defined between stations S1 and S2, and network topology concept used in the model.

Figure 2 .
Figure 2. Architecture of QUAL-NET and data sources.

Figure 3 .
Figure3.Main variables interdependency in the biogeochemical model Rive and associated coefficients.Gray plain rectangles identify variables describing the benthos component; red rectangles are generic functions often called within the code.Variables and function names and their extensive definitions can be found inGarnier et al. (2002).

Figure 5 .
Figure5.Lagrangian profiles from S1 to S2 of TSS, NO − 3 , SRP, Si, PHY and O 2 during two selected events.For the summer event, P input from mineralization processes, P uptake from phytoplankton, phytoplankton growth rate (availability of intracellular carbon and nutrients), sedimentation and mortality rates are also displayed on the right axis.

Figure 6 .
Figure 6.Longitudinal evolution of discharge Q, TSS, SRP and PHY concentrations when a storm event occurred between 8 and 19 August 2013.
-p o in t so u rc e s T ri b s.W W T P D if fu si o n fr o m b e n th o s P h y .u p ta k e B a c t. u p ta k e M in e ra l. u ta ri e s W W T P N o n -p o in t so u rc e s T ri b s.W W T P D if fu si o n fr o m b e n th o

Figure 8 .
Figure8.Average winter and summer budgets between S1 and S2 for phytoplankton biomass, dissolved and particulate organic carbon, and dissolved oxygen.All arrows widths are proportional to calculated fluxes, allowing the visual comparison between winter and summer periods.

Table 2 .
Model performances (bias ± SD errors) for different timescales: over the entire period of validation (1 August 2012 to 31 July 2014), in summer (April to October) and in winter (November to March).
High temporal resolutions provided reasonable daily variations and made it possible to estimate biogeochemical variations during short-