Data-Driven System Dynamics Model for Simulating Water Quantity and Quality in Peri-Urban Streams

Holistic water quality models to support decision-making in lowland catchments with competing stakeholder perspectives are still limited. To address this gap, an integrated system dynamics model for water quantity and quality (including stream temperature, dissolved oxygen, and macronutrients) was developed. Adaptable plug-n-play modules handle the complexity (sources, pathways) related to both urban and agricultural/natural land-use features. The model was applied in a data-rich catchment to uncover key insights into the dynamics governing water quality in a peri-urban stream. Performance indicators demonstrate the model successfully captured key water quantity/quality variations and interactions (with, e.g., Nash-Sutcliff Efficiency ranging from very good to satisfactory). Model simulation and sensitivity results could then highlight the influence of stream temperature variations and enhanced heterotrophic respiration in summer, causing low dissolved oxygen levels and potentially affecting ecological quality. Probabilistic uncertainty results combined with a rich dataset show high potential for ammonium uptake in the macrophyte-dominated reach. The results further suggest phosphorus remobilization from streambed sediment could become an important diffuse nutrient source should other sources (e.g., urban effluents) be mitigated. These findings are especially important for the design of green transition solutions, where single-objective management strategies may negatively impact aquatic ecosystems.


Introduction
Surface water ecosystems worldwide are deteriorating at an alarming rate under ever-increasing human pressures and climate change [1][2][3], threatening biodiversity and ecosystem services that link to human water security and public health [4][5][6][7]. Notably, increases in agricultural productivity, urbanization, and their associated impacts have multiplied the number and severity of stressors to surface waters since the middle of the 20th century. Anthropogenic changes to the soil surface and sub-surface (e.g., drainage network) have resulted in important flow alterations that can have direct ecological consequences or indirect consequences via water quality impairment stemming from enhanced loads of nutrients or pollutants [8][9][10]. Bioeconomy-related pressures resulting from the drive towards green transition solutions in response to climate change threats may pose additional threats to stream water quality [11].
Streams draining peri-urban landscapes (as defined in [12]) are a prime example of such systems, with great potential to be impacted by the heterogeneous sprawl of urban, industrial, and agricultural activities coexisting with natural areas that may be found sporadically throughout a given catchment. The combination of various hydrological

Model Description
The SD model is comprised of a scalable, combined hydrologic and in-stream water quality model. Overall, the model aims to simulate the water flow and associated stream depth, incorporating key natural and urban-related hydrological processes as described further below. Fundamentally, it is based on a representation of fully mixed connected reaches, through which different estimated hydrological flow components (groundwater, urban-related flows, tributaries) are aggregated and routed along ( Figure 1). The geometry of the stream reach is simplified to an ideal rectangular cross-section, which is deemed sufficient considering the bathymetric data available. The natural and anthropogenic features of the hydrological cycle will combine to influence the various in-stream physicochemical parameters. Currently, dissolved oxygen (DO), stream temperature (temp), nitrate (NO 3 ), ammonium/ammonia (NH 4 ), soluble reactive phosphorus (PO 4 ), as well as chlorophyll-a (chl-a)-used as a proxy for suspended algae and benthic plant biomass (macrophytes and benthic algae), are simulated in terms of water quality parameters. ment model can thus be quickly assembled by first creating the required number of reaches and connecting all relevant input and output variables between modules of interest within a reach, and then between the reaches, providing a flexible environment for fast and transparent model development ( Supplementary Information, Figure S1). A condensed description of the model and simulated processes is provided in more detail below and in Appendix A. Conceptual model structure for a single reach encompassing the key hydrological (upper) and in-stream water quality (lower) processes. Key interconnecting processes for water quality given as numbers include: (1) carbonaceous biochemical oxygen demand, (2) reaeration, (3) nitrification, (4) sediment oxygen demand, (5) photosynthesis/autotrophic respiration, (6) nutrient assimilation, (7) denitrification, (8) settling, (9) ammonium/ammonia partitioning (pH-dependent), and (10) water-atmosphere heat exchange. All processes are temperature-dependent. Abbreviations are defined in the related Sections 2.2.1 and 2.2.2. Figure 1. Conceptual model structure for a single reach encompassing the key hydrological (upper) and in-stream water quality (lower) processes. Key interconnecting processes for water quality given as numbers include: (1) carbonaceous biochemical oxygen demand, (2) reaeration, (3) nitrification, (4) sediment oxygen demand, (5) photosynthesis/autotrophic respiration, (6) nutrient assimilation, (7) denitrification, (8) settling, (9) ammonium/ammonia partitioning (pH-dependent), and (10) water-atmosphere heat exchange. All processes are temperature-dependent. Abbreviations are defined in the related Sections 2.2.1 and 2.2.2.
All processes and associated sub-processes were developed as separate modules within the provided software interface, enabling the plug-n-play environment. A catchment model can thus be quickly assembled by first creating the required number of reaches and connecting all relevant input and output variables between modules of interest within a reach, and then between the reaches, providing a flexible environment for fast and transparent model development ( Supplementary Information, Figure S1). A condensed description of the model and simulated processes is provided in more detail below and in Appendix A. The hydrological cycle found in peri-urban streams is comprised of both a natural component, as well as regulated (e.g., wastewater effluent outlets) and unregulated (e.g., separate system outlets) components representing the anthropogenic modification of the water cycle. The hydrological model is driven by precipitation, air temperature and extraterrestrial radiation, equivalent to a lumped rainfall-runoff (RR) formulation. When applied as several modules in series as presented here, however, it can be seen as a semi-distributed model (or series of lumped models), capable of representing key features specific to smaller sub-catchments that may have different governing parameters. To account for the urban component, precipitation can then be partitioned between the natural/agricultural (or pervious) and urban (or impervious) components estimated within the catchment, which is allocated using an aggregated coefficient of imperviousness (f imp ) (Equations (A2) and (A17)) varying between 0 (natural sub-catchment) and 1 (fully impervious/urban catchment). The general water balance for a given reach is: (1) where Q out is the flow at the outlet of a reach; Q in,routed is the routed outflow from the previous reach and routed flow components (tributaries discharging within the reach, urban water flows); RR streamflow represents the flow component generated by the precipitation falling over pervious areas and handled by the RR model; and finally, outflows, j, represents the sum of any water withdrawals (e.g., direct water abstraction, losing reach section). More specifically, the RR model is inspired by existing lump formulations [52], accounting for soil moisture, groundwater storage, and an extra runoff stock ( Figure 1; Equations (A1), (A4) and (A13)). The contribution resulting from the impervious surfaces, corresponding to the urban-related outflows (i.e., separate systems and combined sewer overflows) are estimated using an equivalent drainage area estimated from spatially aggregated data of imperviousness, network extent in the sub-catchment, and dedicated reservoirs (Equations (A17)-(A19)). Combined sewer overflows are specifically modeled as a simple reservoir but with a nonlinear outflow that becomes active above a certain volume threshold; this is incorporated in the reach as a point inflow (other, Figure 1; Equation (A23)). Wastewater treatment plants (WWTP) effluents, constituting point discharges to the stream, are currently introduced as a time series. All the urban inflows are aggregated at the inlet of the reach and routed using a 3rd-order delay function (equivalent to a 3rd-order Nash-cascade model) before being merged with the RR model output. This enables hydrograph shape prediction when the inflows and channel characteristics are known. Finally, a leaking term is introduced for each reach to account for a possible losing stream section (outflow from the stream to the groundwater, Equation (A24)). Water depth is inferred from the calculated water flow using Manning's equation with a shallow depth approximation and a variable Manning's coefficient, n, estimate (Equation (A25)). The latter is dynamically updated on a daily basis in the model using a power regression law, with the estimated streamflow as the independent variable (Equation (A27)) based on the findings from [53]. This approach is deemed reasonable in small streams for which Manning's coefficient (and consequently depth) are strongly affected by the influential seasonal macrophyte coverage occurring at low flow in the summer but less at high discharge ( Figure S5).

Physico-Chemical Conditions and Water Quality
All dissolved substance concentrations, C, in the stream model are simulated using fully-mixed mass balance equations with relevant source and sink terms, S, that are often temperature-dependent (all terms expressed in [mass/time]): Water 2021, 13, 3002 6 of 40 where V represents the volume of water in a given reach, Q in C in are the substance loading terms, C is the fully-mixed concentration of the substance in the reach, and Q in and Q out represent, respectively, the inflow and outflow in the reach. The substance loadings are calculated from estimated/input flow components combined with measured or representative concentration levels at the stream interface, i.e., the source-pathway and possible attenuation is not explicitly modeled.
The stream water temperature model is based on an instantaneous mixing model using temperature loading in the investigated reach, the temperature loading estimated from the previous reach, and the equilibrium temperature concept for the heat exchange at the stream/air interface [54,55]: where ρ is the water density, C p is the water-specific heat capacity [MJ·kg −1 · • C −1 ], A is the area of river/atmosphere interface [m 2 ], H is the net atmospheric flux [MJ·m 2 ·day −1 ], T i , Q i are the water temperature and inflows to the reach respectively, and Q out is the outflow of the reach. Temperature loadings, T i Q i , are assessed from the different calculated or input flow components combined to temperature estimate provided as time series or assessed using a regression model based on air temperature ( Figure S6). DO in the stream is depleted by carbonaceous biological oxygen demand (cBOD, Equation (A37)), nitrification (NBOD, Equation (A43)), a constant "background" sediment oxygen demand (SOD, Equation (A54)), and an additional dynamic heterotrophic respiration term (Equation (A72)). This term accounts for enhanced heterotrophic activities, stemming from enhanced settling of fine organic matter or plant exudation in streams with a high density of water plants, as highlighted in [56]. DO is also affected by photosynthesis/autotrophic plant respiration (Equations (A68)-(A72)) and balanced by a reaeration process. To prevent negative DO concentrations from occurring, we implemented a simple feedback control mechanism by linking the DO to the oxidation rate for organic matter (modeled as cBOD) following [57] (Equation (A40)). The potential DO diurnal variation stemming from plant biomass photosynthesis, and autotropic respiration is simulated on an hourly basis and determined by the photoperiod, estimated max. photosynthesis rate on a daily basis and idealized diurnal cycles [58] (Equation (A70)). Finally, DO saturation is computed from temperature, salinity, and altitude [59] (Equations (A36)-(A38)), and reaeration is driven by the Owens-Gibbs formulation (Equation (A35)).
The model simulates dissolved orthophosphate and inorganic nitrogen (in nitrate and ammonium forms) as nutrients received in the stream water. The nitrification process is simplified to a one-step process in which the first oxidation to nitrite is omitted due to its fast reaction rate (Equation (A46)), while on the other hand, nitrate is potentially removed by denitrification [60] (Equation (A51)). Ammonium/ammonia partitioning is evaluated using an equilibrium reaction based on pH data [61] (Equation (A48)). pH is currently incorporated directly as input data to the model and is not calculated. The nutrient N, P pool is potentially depleted via assimilation by the plant biomass based on mass stoichiometric ratios or simply transported further downstream. Nutrient fluxes from the sediment are currently not accounted for in this model version.
The plant biomass in the stream water contributing to photosynthesis and respiration is split into two different stocks, corresponding to two separate autotroph groups according to their mobility: phytoplankton, which is transported in the water column, and fixed aquatic plant communities (e.g., macrophytes, benthic algae), which are pooled together. For both stocks, the plant biomass is reduced by a combined respiration/excretion and non-predatory mortality term (Equations (A55)-(A56)). Phytoplankton, additionally, can settle (Equation (A57)). The respiration/excretion and mortality terms are only temperature-dependent [62] (Equation (A58)), whereas gross growth is controlled and limited by environmental conditions including temperature (Equation (A59)), nutrient availability (minimum of N, P pool following a Michaelis-Menten formulation, Equa-tion (A61)), and light (Beer-Lambert light attenuation integrated over time and depth (Equation (A62)). The light attenuation coefficient is dynamically estimated and linearly dependent on both suspended particles (non-algal suspended solids as background attenuation) and the phytoplankton concentration, but also non-linearly related to macrophyte biomass, following the regression model built by [63] on various plant communities (Equation (A65)).
The daily photosynthetic rate is a linear function of the gross growth rate for the different aquatic plant stocks, identical for all stocks (Equation (A68)), and the corresponding autotrophic respiration is estimated as a fraction of this daily photosynthetic rate [64] (Equation (A69)).
It is worth noting that the decaying of plant materials after death is so far not included. Hence, the resulting nutrient flux from this decomposition, and the nutrient fluxes from the streambed in general, are not accounted for. However, the oxygen required for the decomposition of the plant matter is essentially accounted for in the SOD calculations in the current version of this model.

Model Uncertainty
Stream water quality is characterized by high spatio-temporal variability, and associated models are often applied in data-scarce environments [26]. These variations stem, for example, from variations in landscape characteristics and catchment topography in the spatial dimension [18], whereas changes in environmental conditions coupled to variations in the source, active pathways, and related processes will influence the temporal variability at hourly, daily, and seasonal scales [14]. Notably, considering these variations and related uncertainty in any model output may ultimately become crucial, especially if the model results are used in any kind of decision-making process [65].
Due to the limited availability of associated data, the uncertainty and potential spatiotemporal variations of the physico-chemical conditions is addressed using the Monte Carlo (MC) simulation capabilities offered in the STELLA software. Specifically, Latin Hypercube Sampling technique [66] was used to ensure a proper spread of possible parameter values, including the extremes, are sampled within the defined parameter space. For each of the tracked physico-chemical conditions, we carried out 200 realizations, i.e., >10 times the number of uncertain parameters (Tables S5 and S6) following the recommendations for this sampling technique [67]. The uncertainty arising from the model structure itself is not evaluated here, considering the structure presented in this study was developed based on well-accepted processes. However, possible omitted processes will be discussed in Section 5.

The Usserød Peri-Urban Catchment
The model was applied to a small low-land catchment (ca. 120 km 2 , mean elevation 30.8 m.a.s.l (DVR.90) located on Sjaelland, 25 km north of Copenhagen, Denmark ( Figure 2). This catchment is drained by Usserød Stream, flowing from its origin at Sjael Lake (via automatic sluice control) before merging with the Nivå Tributary and discharging in the Baltic Sea ca. 8000 m to the north. The stream is considered as a small-to-medium stream at the national level, with width ranging from ca. 2 to 13 m and depth from 0.2 to 1.0 m. The median flow is ca. 348 L/s −1 for the period 2017-2019 [68] (St.5005). The geology in the catchment consists of clay tills with embedded sand lenses/layers, which is typical for this part of Denmark. Groundwater is abstracted from the underlying Danien limestone aquifer ( Figure S2). A detailed description of the catchment can be found in [39].
This catchment is a typical peri-urban area, with a complex mixture of land-use activities spreading along the stream corridor. Urban (dwelling and industrial), agricultural, and natural-like areas (mostly secondary forest and surface water) cover 22, 57, and 21% of the catchment, respectively. The Usserød Stream receives additional flow from the Donse Tributary, draining mostly natural areas and agricultural lands, before running through more agricultural lands in the north direction further downstream. Urban areas are predominant in the upstream part of the catchment, where two sluices with overflow and one low head dam with a side derivation ensure a possible flow regulation. These urban areas are drained by both combined sewer overflow systems (CSO) and separate storm (rainwater) overflow systems (SSO), the latter discharging stormwater directly into the stream ( Figure S2). Approximately 11 CSO structures are found along the Usserød stream for a direct discharge of excess sewage water in case of overload of the combined system network ( Figure S3). Wastewater in the catchment is treated by three wastewater treatment plants releasing effluent into the stream.
Former monitoring activities in the catchment revealed that some physico-chemical and water quality parameters (e.g., dissolved oxygen, temperature) were below the local guideline values, with potential impairment of the ecological status of the stream, as documented in [69,70].
Water 2021, 13, x FOR PEER REVIEW 8 of 40 more agricultural lands in the north direction further downstream. Urban areas are predominant in the upstream part of the catchment, where two sluices with overflow and one low head dam with a side derivation ensure a possible flow regulation. These urban areas are drained by both combined sewer overflow systems (CSO) and separate storm (rainwater) overflow systems (SSO), the latter discharging stormwater directly into the stream ( Figure S2). Approximately 11 CSO structures are found along the Usserø d stream for a direct discharge of excess sewage water in case of overload of the combined system network ( Figure S3). Wastewater in the catchment is treated by three wastewater treatment plants releasing effluent into the stream. Former monitoring activities in the catchment revealed that some physico-chemical and water quality parameters (e.g., dissolved oxygen, temperature) were below the local guideline values, with potential impairment of the ecological status of the stream, as documented in [69,70].  [71] and abstraction well location along the stream, (b) key urban features (WWTP) and available monitoring network. The reach boundaries are also displayed and set by the monitoring stations for verification purposes [72]. Each reach is named after the corresponding downstream monitoring station in this study.

Available Input Data and Monitoring
A network of continuous monitoring stations deployed in the catchment provided hourly to daily time series for several of the streams' hydrological and physical-chemical properties used for calibration and verification (Table S1, Figure 2). Furthermore, the dual measurement flow and water level/depth by individual sensors enabled the estimation of an equivalent Manning's coefficient and the regression parameters employed to characterize its dynamic variations in our model ( Figure S5).
The different sub-catchments and associated reach characteristics (i.e., elevation, land cover, reach length) were retrieved from [71] for delineation, and zonal statistics operations were performed in QGIS (v2.18.28, TauDEM package). Meteorological data (air temperature and precipitation time series) were collected from a weather station located in  [71] and abstraction well location along the stream, (b) key urban features (WWTP) and available monitoring network. The reach boundaries are also displayed and set by the monitoring stations for verification purposes [72]. Each reach is named after the corresponding downstream monitoring station in this study.

Available Input Data and Monitoring
A network of continuous monitoring stations deployed in the catchment provided hourly to daily time series for several of the streams' hydrological and physical-chemical properties used for calibration and verification (Table S1, Figure 2). Furthermore, the dual measurement flow and water level/depth by individual sensors enabled the estimation of an equivalent Manning's coefficient and the regression parameters employed to characterize its dynamic variations in our model ( Figure S5).
The different sub-catchments and associated reach characteristics (i.e., elevation, land cover, reach length) were retrieved from [71] for delineation, and zonal statistics operations were performed in QGIS (v2.18.28, TauDEM package). Meteorological data (air temperature and precipitation time series) were collected from a weather station located in the catchment [73] (St.5622), while cloud cover was collected from the closest station with available data [74] (St.06188, ca. 12 km from the catchment).
Details and technical information about the drainage of urban areas were gathered from the local drinking water and wastewater company with GIS shape files displaying the average degree of imperviousness and sewer network type at a spatial resolution close to cadastral units ( Figures S2-S4). Groundwater extraction data were received as time series on a daily basis and aggregated per reach. The different WWTP influent and effluent flows, as well as some associated physico-chemical conditions (e.g., temperature, oxygen saturation, nutrient concentrations), were provided as daily time series for the largest WWTP and at a lower and variable frequency for the two smallest ones (Table S2). Almost all CSOs in this catchment are located in the reach Grønnegade-Ådalsvej (except one in the upstream reach Sjaelsø-Grønnegade), and their overflow volume is closely monitored by the local municipality [75]. Other data specifically focusing on water quality (nutrients, chl-a) were available from temporal monitoring (grab samples at monthly to bi-monthly intervals) at 11 different sampling stations during 2018 and 2019 [39]. The remaining data were extracted from literature sources: nutrient concentrations in some specific compartment, BOD, and DO saturation for urban drainage compartment (SSOs and CSOs), as well as most data related to stream aquatic plant biomass. All input data and model parameters are provided in Tables S2 and S3.

General Set-Up
The Usserød Stream was discretized into three reaches in this catchment so that the reach outlets coincided with the available monitoring stations and judiciously included the various different features and input discharges. The reach lengths range from 1500 to 3500 m and are designated (named) according to their outlet points through the rest of this paper (i.e., Grønnegade, Ådalsvej, Parallelvej, Nivemølle, Figure 2b). The length range theoretically ensures fully mixed conditions of any substances discharging at the upgradient boundary and are short enough to consider the reach geometry as uniform. It is worth noting that an additional reach was calculated, i.e., the sub-reach Parallelvej (Figure 2b), to facilitate verification of the DO model result with the only available DO concentration measurement (sensor) in the catchment. The Donse Tributary contribution to the streamflow at Brønsholmsdalsvej was estimated using the developed RR module. Reach characteristics can be found in Table S4.
We used three consecutive years of monitoring station data (2017-2019) for calibration and verification: the years 2017 and 2018 for the hydrological calibration (these years were respectively wet and dry, with annual precipitation of 777 and 512 mm, average cumulated = 712 mm) and year 2019 for the verification (cumulated precipitation = 791 mm). The results from the aforementioned field investigation provided a dataset to assess the model's ability to simulate the dynamic variations of the physico-chemical conditions [39]. The pH monitored at the most upstream part of the catchment was used as input and assumed uniform throughout all reaches. Such an assumption allows a simulation of extreme cases as we can suspect the highest variations of pH stemming from the eutrophic conditions at the source of the catchment and more subtle variation downstream due to the more controlled variation imposed by the WWTP operations. The model is run using a daily time-step.

Model Calibration Procedure and Application
The RR model and the CSO sub-module model were automatically calibrated using the dedicated STELLA module based on a differential evolution algorithm and leastsquares method for error minimization. The calibration was first run on the tributary and its associated sub-catchment (Fredtoften, Figure 2b), i.e., the most pervious part of the catchment and relatively free from any urban disturbance and for which the least data were available. In a second step, this simulation was run for the entire tributary sub-catchment (Brønsholmsdalsvej, Figure 2b) with the inclusion of the urban areas to verify its satisfactory performance in peri-urban settings. The calibrated parameter values were used as default parameter values in the RR model for the Usserød stream reaches, assuming relatively uniform geohydrological conditions over such a small catchment size. The groundwater storage time constant required an extra calibration operation (carried out on the most upstream reach and assumed equal in the other reaches) due to its strong influence on the hydrological response ( Figure S14). The corresponding initial stock value was here manually adjusted to fit the start of the simulation period. All other stock initial values were simply initialized to non-null values, as their influences were negligible beyond the first few days of simulation. Finally, a non-null leaking term was introduced and calibrated in the most downstream reach to account for a losing section (flow from the stream towards groundwater) in periods of low flow. The CSO module calibration was carried out using the aggregated time series of all CSOs events (arithmetic summation) in this reach, considering the lump character of our model. The single CSO structure present in the Grønnegade reach was neglected at this stage. Calibrated values can be found in Table S4.
The automatic calibration for all physico-chemical parameters of interest (i.e., DO, temp, NH 3 , NO 3 , PO4, chl-a) is hindered by the absence of a multi-objective calibration procedure, complex feedback mechanisms, and limited measurement data. Instead, we used a multiple steps calibration to get a deterministic solution for the DO concentration, based on the identified most sensitive parameters ( Figure S18): DO saturation of the different flow components, heterotrophic respiration, followed by BOD degradation rate and nitrification. The oxygen yield parameter for the plant biomass was set to realistic values using the value range retrieved from [63]. Uncertainty and time variation of parameters and inputs were then addressed by Monte Carlo Simulation. No specific calibration was used for the temperature simulation that is mostly relying on the input regression model built from the measurements in the catchment ( Figure S6). Finally, the ability of the model to handle nutrient and chl-a concentrations were carried out via a forcing input function at Grønnegade station, based on simulated flow, measured concentration data [39], and MC simulation to address the parameter uncertainty and the current capabilities of the model. An overview of the overall set-up and calibration procedure of the model can be found in Figure 3. were available. In a second step, this simulation was run for the entire tributary sub-catchment (Brønsholmsdalsvej, Figure 2b) with the inclusion of the urban areas to verify its satisfactory performance in peri-urban settings. The calibrated parameter values were used as default parameter values in the RR model for the Usserød stream reaches, assuming relatively uniform geohydrological conditions over such a small catchment size. The groundwater storage time constant required an extra calibration operation (carried out on the most upstream reach and assumed equal in the other reaches) due to its strong influence on the hydrological response ( Figure S14). The corresponding initial stock value was here manually adjusted to fit the start of the simulation period. All other stock initial values were simply initialized to non-null values, as their influences were negligible beyond the first few days of simulation. Finally, a non-null leaking term was introduced and calibrated in the most downstream reach to account for a losing section (flow from the stream towards groundwater) in periods of low flow. The CSO module calibration was carried out using the aggregated time series of all CSOs events (arithmetic summation) in this reach, considering the lump character of our model. The single CSO structure present in the Grønnegade reach was neglected at this stage. Calibrated values can be found in Table  S4.
The automatic calibration for all physico-chemical parameters of interest (i.e., DO, temp, NH3, NO3, PO4, chl-a) is hindered by the absence of a multi-objective calibration procedure, complex feedback mechanisms, and limited measurement data. Instead, we used a multiple steps calibration to get a deterministic solution for the DO concentration, based on the identified most sensitive parameters ( Figure S18): DO saturation of the different flow components, heterotrophic respiration, followed by BOD degradation rate and nitrification. The oxygen yield parameter for the plant biomass was set to realistic values using the value range retrieved from [63]. Uncertainty and time variation of parameters and inputs were then addressed by Monte Carlo Simulation. No specific calibration was used for the temperature simulation that is mostly relying on the input regression model built from the measurements in the catchment ( Figure S6). Finally, the ability of the model to handle nutrient and chl-a concentrations were carried out via a forcing input function at Grønnegade station, based on simulated flow, measured concentration data [39], and MC simulation to address the parameter uncertainty and the current capabilities of the model. An overview of the overall set-up and calibration procedure of the model can be found in Figure 3. . Flowchart depicting the data-driven application of the model. Calibration and verification for the water quantity module is carried out first, followed by the water quality module, where time-series data are available. Sensitivity analyses can then be conducted using the calibrated model to determine key water quantity (see Figure S17) and water quality ( Figure S18) variables. These can be used further in a probabilistic application of the model, where the Monte Carlo (MC) method is applied. All values for these simulations have been derived either from the literature or site-specific (grab sample) data (see Tables S5 and S6 for values and sources).

Statistic Indicators Used for Comparison between Model Results and Data
The model's ability to capture stream water quantity and quality trends was evaluated using the comparison of plotted time series combined with standard statistical indicators: RMSE values and coefficient of determination R 2 were computed for all quantities and specifically only when discrete measurement data were available. Nash-Sutcliffe Efficiency (NSE), Percent BIAS (PBIAS), and RMSE-observations standard deviation ratio (RSR) for hydrological quantities monitored continuously [76]. Finally, we used a flashiness index to characterize the impact of the urban compartment on the peri-urban stream flow, Q, and alteration of the flow regime [77].

Results
To demonstrate the model applicability, results for both the calibration and verification period are shown for various water quantity and quality indicators in the following sections. Taking advantage of the numerous stream gauging stations in the catchment, results could be shown for four locations along the Usserød stream (Grønnegade; Ådalsvej, Parallelvej, Nivemølle) and two locations along the Donse Tributary (Fredtoften, Brønsholmsdalsvej), listed here in order from up-to downstream flow locations along each water course ( Figure 2). Results have been selected to show the versatility and transferability within the catchment, as well as the model performance. In terms of water quality, Parallelvej (Usserød Stream) is shown for DO, the location of the sensor; else results are shown according to the outlets for the three reaches comprising the Usserød catchment (Grønnegade, Ådalsvej, Nivemølle).

Stream Flow and Water Depth
The model performs quite well in terms of flow simulation for Usserød Stream (Figure 4a,b and Figure S7), with the performance indicators categorized as good or very good for all reaches in terms of flow (NSE, RSR, PBIAS) for both the calibration and verification periods (Table 1). In terms of the stream water level (depth), the model was also found to simulate the dynamics well, both on a daily and seasonal basis, with most indicators in the satisfactory to good range (except at Grønnegade), and with RMSE values below 0.08 m for both calibration and verification periods (Table 1, Figure 4a and Figure S8), indicative for the capability of the model in capturing seasonal variations of stream roughness.
Notably, during the entire simulation period, 2017-2019, 49 CSOs events were recorded in the Ådalsgade reach (15 below 10 m 3 /day), of which 11 were dynamically captured by our model without prior specific knowledge of the network structure (actual number of outlets) nor actual active control measures (related to CSO release requirements). Deviation of estimated and measured overflows was, however, relatively high (59% relative deviation for the events simulated) ( Figure S9). Importantly, the simulation underlines the variable contribution and dynamics of the different flow components. In winter, an important part of the flow originates from the lake outflow (input) and the pervious areas. In summer, the tributary at Brønsholmdalsvej becomes negligible, while the WWTP effluents contribute significantly to the overall stream flow. Stormwater becomes an important contributor during rain events via the numerous separate system outlets to the stream. CSOs are sporadically active, but their respective short and intense discharges contribute relatively little to the flow when aggregated on a daily basis (see discussion in Section 5.1). The RR model application comprising the full sub-catchment area (all feature types) of the Donse Tributary (Brønsholmdalsvej) is shown in Figure 4c for both the calibration and verification period. Satisfactory agreement between simulated and measured data could be documented for the calibration period, based on the NSE and RSR criteria, while PBIAS reached an unsatisfactory value of 30%, caused by an underestimation of the flow in the second year of the calibration period (Table 1). These all improved to very good for the verification period, where the effect of the urban drainage system on the flow is particularly visible. Even on a daily time-step, sharp peaks and an increase in the flashiness index (RB) at Brønsholmdalsvej are representative of the fast drainage of impervious areas via separate systems in this part of the sub-catchment (RB = 0.13 and 0.15 at Fredtoften The RR model application comprising the full sub-catchment area (all feature types) of the Donse Tributary (Brønsholmdalsvej) is shown in Figure 4c for both the calibration and verification period. Satisfactory agreement between simulated and measured data could be documented for the calibration period, based on the NSE and RSR criteria, while PBIAS reached an unsatisfactory value of 30%, caused by an underestimation of the flow in the second year of the calibration period (Table 1). These all improved to very good for the verification period, where the effect of the urban drainage system on the flow is particularly visible. Even on a daily time-step, sharp peaks and an increase in the flashiness index (RB) at Brønsholmdalsvej are representative of the fast drainage of impervious areas via separate systems in this part of the sub-catchment (RB = 0.13 and 0.15 at Fredtoften and Brønsholmdalsvej, respectively). In comparison, for the upstream most pervious (agricultural/natural) part of this sub-catchment (Fredtoften), the flow dynamics could be classified as good or very good (RSR, NSE) for the calibration period, and good (NSE, RSR) for the verification period (Table 1, Figure S10). Table 1. Performance indicators for flow and depth simulations, including the Nash-Sutcliffe Efficiency (NSE), Percent BIAS (PBIAS), RMSE-observations standard deviation ratio (RSR), and Flashiness index (RB, flow only) for all reaches, for both the calibration and verification period. Performance rating values for all indicators, except RB are indicated as follows [78]: very good (in bold); good (in italic); satisfactory; unsatisfactory (marked with an x). Additional simulation results can be found in Figures S7-S9 The PBIAS indicator, however, highlights specific periods of flow underestimated and overestimated. For example, the flow at Fredtoften station is overestimated at the end of the summer-fall season in 2018 and 2019 of the verification period ( Figure S10). Such a discrepancy likely stems from the source of the tributary, consisting of an overflow from a water impoundment not accounted for in the model, which becomes insignificant in the late summer time. Conversely, PBIAS indicates a flow underestimation during the calibration period at Brønsholmdalsvej. This deviation could potentially be caused by an extension of the separate system's network, not aligned with the natural delineation of the subcatchment, and therefore collecting and adding stormwater from an adjacent sub-catchment. Nevertheless, we believe these results confirm the implemented approach is suitable for capturing the flow dynamics from urban and agricultural/natural catchment features.

Physico-Chemical and Water Quality Parameters
The simulation of daily stream water temperature compares very well with the measurements for both the calibration and verification period at Parallelvej ( Figure 5; R 2 = 0.90, RMSE = 1.6 • C). Moreover, the narrow percentile spread for the MC simulation results (results not shown) suggests that the uncertainty related to the relevant model parameters for simulating stream temperature is not critical (i.e., low sensitivity parameters) and that the daily stream temperature variations are mostly driven by the daily variation in air temperature. Similar observations were made in other studies ( [54,79]) where the small stream size and relatively shallow depth resulted in a low water thermal inertia and fast equilibrium of water temperature, leading in turn to important seasonal water temperature variations in the investigated streams. urements for both the calibration and verification period at Parallelvej ( Figure 5; R 2 = 0.90, RMSE = 1.6 °C). Moreover, the narrow percentile spread for the MC simulation results (results not shown) suggests that the uncertainty related to the relevant model parameters for simulating stream temperature is not critical (i.e., low sensitivity parameters) and that the daily stream temperature variations are mostly driven by the daily variation in air temperature. Similar observations were made in other studies ( [54,79]) where the small stream size and relatively shallow depth resulted in a low water thermal inertia and fast equilibrium of water temperature, leading in turn to important seasonal water temperature variations in the investigated streams.

Dissolved Oxygen Concentration
Dissolved oxygen concentrations in the stream are relatively well captured by the model with a coefficient of determination R 2 for the daily average concentration ranging from 0.58 to 0.63 for the calibration and verification period, respectively, compared to the sensor data at Parallelvej (Figure 6a,b). The computation at a sub-daily time-step also compares relatively well with observations in terms of the overall dynamic trend (Figure 6c), albeit with some short-term deviations. The observed discrepancies in amplitude at such a time discretization is likely related to the daily input data (not available at higher frequency), shadow effects not properly captured (either riparian or in-stream), and the uncertainty arising from the aggregation of aquatic plant biomass into two groups, as corroborated by the MC simulation results (most of the variations and the grab samples at Parallelvej and Nivemølle lie within the estimated percentiles, Figures 7 and 8b). The most important discrepancy between simulated and observed daily amplitude is found during the period March-May 2019 (with grab sample measurements lying outside of the 10-90percentile for this period, Figure 7). This deviation is possibly connected to a suspended algae or benthic/epiphytic community bloom stimulated by favorable spring light conditions that particular year and would require an additional dedicated biomass stock in the model. Similar ecological events have been witnessed in other streams with high nutrient levels at similar periods of the year, with, for instance, phytoplankton blooms during the spring period in [80,81] in England. Another algae bloom was also observed at the stream's source Sjael Lake in August 2019 (visual inspection). In terms of time variation, the diel oscillation and max. DO concentration resulting from photosynthetic activity occur relatively late compared to the simulation for point C (Figure 6c) and may possibly

Dissolved Oxygen Concentration
Dissolved oxygen concentrations in the stream are relatively well captured by the model with a coefficient of determination R 2 for the daily average concentration ranging from 0.58 to 0.63 for the calibration and verification period, respectively, compared to the sensor data at Parallelvej (Figure 6a,b). The computation at a sub-daily time-step also compares relatively well with observations in terms of the overall dynamic trend (Figure 6c), albeit with some short-term deviations. The observed discrepancies in amplitude at such a time discretization is likely related to the daily input data (not available at higher frequency), shadow effects not properly captured (either riparian or in-stream), and the uncertainty arising from the aggregation of aquatic plant biomass into two groups, as corroborated by the MC simulation results (most of the variations and the grab samples at Parallelvej and Nivemølle lie within the estimated percentiles, Figures 7 and 8b). The most important discrepancy between simulated and observed daily amplitude is found during the period March-May 2019 (with grab sample measurements lying outside of the 10-90-percentile for this period, Figure 7). This deviation is possibly connected to a suspended algae or benthic/epiphytic community bloom stimulated by favorable spring light conditions that particular year and would require an additional dedicated biomass stock in the model. Similar ecological events have been witnessed in other streams with high nutrient levels at similar periods of the year, with, for instance, phytoplankton blooms during the spring period in [80,81] in England. Another algae bloom was also observed at the stream's source Sjael Lake in August 2019 (visual inspection). In terms of time variation, the diel oscillation and max. DO concentration resulting from photosynthetic activity occur relatively late compared to the simulation for point C (Figure 6c) and may possibly stem from an unknown sensor timestamp shift, local variability in environmental conditions, or high macrophyte density in the most upstream part of the catchment combined with a much lower reaeration rate than estimated (results not shown). This point will be further discussed in the next section. stem from an unknown sensor timestamp shift, local variability in environmental conditions, or high macrophyte density in the most upstream part of the catchment combined with a much lower reaeration rate than estimated (results not shown). This point will be further discussed in the next section. Overall, the model satisfactorily captures the seasonal trend for dissolved oxygen and highlights a potentially critical state in summer (NSE = 0.52 and 0.59 for calibration and verification period, respectively). Low average DO concentrations corresponding to low saturation levels are indeed observed and simulated in the period July-September for both the calibration and verification periods and possibly driven by sustained heterotrophic respiration (Figures 6a and S11). Interestingly, the model, as applied in this study (i.e., in data-driven mode), was unable to capture the low daily concentrations between May and July 2019. Such a dynamic trend does not seem to be related to input data uncertainty, as the MC simulation did not capture these events either (Figures 6a and 7). The overestimation of the simulated DO concentration from May to July could be related to an important heterotrophic respiration event, stemming from (1) scouring, settling, and decomposition of the possible algae bloom in the previous months, or/and (2) a high dissolved organic matter (DOM) load during high flow prior to a sustained low flow period in spring (see the simulated/observed high flow in April 2019, followed by a low flow period, Figure 4). Such a phenomenon was recently reported in [82], where an important DOM load following a high flow and episode was the main explanatory factor for the low DO concentrations observed in the river Thames, UK.  (Figures 8c and S12). The ammonium concentrations seem to  Table S5. underestimation of the suspended chl-a in spring (March-May), coinciding with the aforementioned underestimation of DO concentrations by our model (Figures 6a and 7). These remarks (combined with the high density of emergent aquatic plants observed in this reach) support the assumption of an epiphyte or benthic algae community development at that time, possibly washed away by the stream water flow. Similar observations were made in shallow stream systems in the same country, with epiphyte and benthic microalgae bloom from March to April lost later in the season during high flow episodes [88].   Table S6. Overall, the model satisfactorily captures the seasonal trend for dissolved oxygen and highlights a potentially critical state in summer (NSE = 0.52 and 0.59 for calibration and verification period, respectively). Low average DO concentrations corresponding to low saturation levels are indeed observed and simulated in the period July-September for both the calibration and verification periods and possibly driven by sustained heterotrophic respiration (Figure 6a and Figure S11). Interestingly, the model, as applied in this study (i.e., in data-driven mode), was unable to capture the low daily concentrations between May and July 2019. Such a dynamic trend does not seem to be related to input data uncertainty, as the MC simulation did not capture these events either (Figures 6a and 7). The overestimation of the simulated DO concentration from May to July could be related to an important heterotrophic respiration event, stemming from (1) scouring, settling, and decomposition of the possible algae bloom in the previous months, or/and (2) a high dissolved organic matter (DOM) load during high flow prior to a sustained low flow period in spring (see the simulated/observed high flow in April 2019, followed by a low flow period, Figure 4). Such a phenomenon was recently reported in [82], where an important DOM load following a high flow and episode was the main explanatory factor for the low DO concentrations observed in the river Thames, UK.

Nutrient Species and Chlorophyll-a
In terms of inorganic nitrogen, the estimated nitrate concentrations are relatively well captured, with most of the grab sample concentrations falling within the 10-90 percentile of the MC simulation results (Figure 8c and Figure S12). The ammonium concentrations seem to be well estimated in winter but overestimated in summer (Figure 8d). This overestimation seems to affect only the last simulated reach of the stream ( Figure S13), where an important coverage of macrophytes and emergent plants is observed (data not shown). The relatively good simulation results for the NO 3 -N concentrations rule out a possible underestimation of the nitrification process, suggesting instead an assimilation preference for NH 4 -N (rather than NO 3 -N) by the aquatic plant biomass. Such an assimilation preference has also been reported in [83,84], for instance, using controlled releases of nutrient pulses and stable isotope analysis, respectively.
Orthophosphate concentrations exhibit the opposite trend with a significant underestimation of concentrations during an important part of the summer at the beginning of autumn, specifically in the last simulated reach (June-October, Figure 8e and Figure S14). This result highlights a possible omission of an important source-pathway (e.g., wash-off from agricultural lands, septic tank drainage). Alternatively, and considering the continuous overestimation in the summer period for all associated sampling periods, a diffuse release from phosphorous stored in the streambed sediments in the summertime cannot be excluded. Such processes have been observed and described in lentic environments [85], but also in streams where phosphorous becomes available by desorption and fueled by enhanced decomposition of organic matter and lower oxygen levels at higher temperatures [56,86]. P-release from the sediment in the urban stream environment has also been documented when coupled with low DO concentrations ( [87]), although often in connection with lower NO 3 concentrations compared to the present study. Finally, the model simulation for suspended chl-a compares relatively well with the few available observations, and despite the lack of data and characterization of the suspended algae species assemblage in the stream (Figure 8f). The simulations in the different reaches show that an important part of the suspended chl-a actually originates from the most upstream reach, likely as an output from the eutrophic lake ( Figure S15). At the most downstream station, Nivemølle, the percentile range for the simulations clearly shows an underestimation of the suspended chl-a in spring (March-May), coinciding with the aforementioned underestimation of DO concentrations by our model (Figures 6a and 7). These remarks (combined with the high density of emergent aquatic plants observed in this reach) support the assumption of an epiphyte or benthic algae community development at that time, possibly washed away by the stream water flow. Similar observations were made in shallow stream systems in the same country, with epiphyte and benthic microalgae bloom from March to April lost later in the season during high flow episodes [88].

Discussion
This SD model and lumped module formulation, implemented as a semi-distributed model structure to capture key spatial variability, enabled the computationally quick and effective simulation of hydrology and water quality variations in a small peri-urban stream. Notably, one of the great advantages of the lumped formulation is a fast simulation time, providing a good opportunity to account for data uncertainty and their inherent variability (i.e., through MC simulation).
In terms of the hydrological module, the model and the implemented structure captured well the very dynamic contribution of flows stemming from urban and more pervious areas that will influence water quality (with performance indicators ranging from satisfactory to very good). The good results in terms of water depth (RMSE < 0.1 m) highlight the importance of a dynamic and variable Manning's coefficient, especially for small stream systems like the one modeled here, where shallow water depths combined with underwater vegetation have a strong influence. On the water quality side, the representation of in-stream processes and simulation results for stream temperature, oxygen, and macronutrients (NH 4 -N, NO 3 -N, PO 4 -P) were deemed acceptable compared to the available measured data (based on RMSE and coefficient of determination).
Importantly, the observed discrepancies allowed us to identify key processes affecting the hydrology or water quality in a mixed land-use environment, due in part to the transparency of the model structure allowing the direct inspection of cause-effect system feedbacks (built using an SD approach), as well as potential model improvements and limitations (both discussed further below). This could facilitate a more in-depth analysis of the investigated stream system, which can sometimes be challenging with more complex model structures [89].

Process Understanding and Model Application
The combined simulation supported by a rich dataset for our application gives some insightful information on the dynamics processes at stake in peri-urban catchments with strong implications on the general stream water quality.
Measurement of nutrient fluxes on a seasonal basis in this catchment showed very dynamic contributions (e.g., nutrients discharge) in the different reaches, confounding the identification of detrimental impacts from sources and land use [39]. The simulations highlight this important dynamic contribution of the different flow components (Figure 9), and consequently, the challenges with respect to maintaining good water quality in periurban stream systems. There is not one single dominant component, and the effects from the different contributions and associated loadings are naturally very dependent on the local catchment attributes (e.g., degree of urbanization, land-use, type and extent of drainage system). These local differences need to be accounted for when designing restoration strategies that target stream water quality in peri-urban systems [90]. In our case, the stream reaches could be further divided to check that the current model setup appropriately captures the key attributes governing the conditions in the stream system.
In terms of CSO contribution, it is worth mentioning that the simulation on a daily basis shows that a single CSO event has a quite limited effect in terms of flow due to its relatively short and transient properties (Figure 9b,d). Nevertheless, the impairment related to particle load (and related degradation) and pollutants (dissolved and particle-bounded) is non-negligible for these structures and underlined as a possible reason for ecological degradation in this catchment, but it is out of the scope of the model at this time. In terms of CSO contribution, it is worth mentioning that the simulation on a daily basis shows that a single CSO event has a quite limited effect in terms of flow due to its relatively short and transient properties (Figure 9b,d). Nevertheless, the impairment related to particle load (and related degradation) and pollutants (dissolved and particlebounded) is non-negligible for these structures and underlined as a possible reason for ecological degradation in this catchment, but it is out of the scope of the model at this time.
The results from the stream water quality simulation highlight important seasonal variations of stream water temperature that may be amplified in peri-urban settings. Notably, here, temperatures above 21 degrees Celsius (local threshold value) are estimated during the summer season, with potential adverse effects in terms of DO concentration (due to decrease in DO saturation) and on the stream biota [91]. The model indicates a limited stream thermal inertia (driven by the shallow depth and fast equilibrium rate with the atmosphere, not shown), resulting in the stream water temperature being strongly correlated to the air temperature as the main explanation. Furthermore, the extent of urban and impervious areas and relative contributions to the peri-urban stream system can exacerbate this problem. Generally, air temperatures in the urban areas are higher, and runoff on heated impervious areas drained via a separated system can contribute as well [92].
The simulation of DO concentrations shows high heterotrophic respiration events during the summer and beginning of autumn (May-October) that appear as the cause of sustained periods at relatively low daily DO levels witnessed in this catchment. This observation is in agreement with freshwater ecosystem metabolism being generally heterotrophic [64]. However, this heterotrophic state seems highly dynamic in time and space and does not always result in significant DO depletion (see, e.g., the difference in measured DO saturation levels between two consecutive years, 2018 and 2019, in Figure S11). Several studies highlighted that urbanization leads to increased and highly dynamic loads of more labile dissolved organic matter fraction and enhanced heterotrophic respiration [93,94], which will constitute a challenge in peri-urban settings considering the myriad of The results from the stream water quality simulation highlight important seasonal variations of stream water temperature that may be amplified in peri-urban settings. Notably, here, temperatures above 21 degrees Celsius (local threshold value) are estimated during the summer season, with potential adverse effects in terms of DO concentration (due to decrease in DO saturation) and on the stream biota [91]. The model indicates a limited stream thermal inertia (driven by the shallow depth and fast equilibrium rate with the atmosphere, not shown), resulting in the stream water temperature being strongly correlated to the air temperature as the main explanation. Furthermore, the extent of urban and impervious areas and relative contributions to the peri-urban stream system can exacerbate this problem. Generally, air temperatures in the urban areas are higher, and runoff on heated impervious areas drained via a separated system can contribute as well [92].
The simulation of DO concentrations shows high heterotrophic respiration events during the summer and beginning of autumn (May-October) that appear as the cause of sustained periods at relatively low daily DO levels witnessed in this catchment. This observation is in agreement with freshwater ecosystem metabolism being generally heterotrophic [64]. However, this heterotrophic state seems highly dynamic in time and space and does not always result in significant DO depletion (see, e.g., the difference in measured DO saturation levels between two consecutive years, 2018 and 2019, in Figure S11). Several studies highlighted that urbanization leads to increased and highly dynamic loads of more labile dissolved organic matter fraction and enhanced heterotrophic respiration [93,94], which will constitute a challenge in peri-urban settings considering the myriad of delivery pathways. Another possible driving mechanism could be algae blooms phenomena followed by settling and decomposition (as discussed in Section 4.2.2) or aquatic plant biomass, with macrophytes strongly affecting the local flow conditions, resulting in flow velocity reduction, fine sediment accumulation, and ultimately, enhanced metabolism in lowland streams, as shown in [56] for another small shallow stream. To this day, the modeling of this type of interaction and understanding of heterotrophic respiration is still limited (but see [95] and a coupled model of stream metabolism and microbial biomass, calibrated on long-term monitoring series ) and constitute an important source of uncertainty for any DO simulation in water quality models [82].
Finally, the sustained period of low dissolved oxygen previously described may also create favorable conditions for the enhanced release of nutrients. The simulations revealed, for example, an underestimation of phosphorus concentrations compared to the observed data. While unknown sources such as agricultural drainage or septic tanks cannot be excluded in peri-urban catchments and should be investigated further, several studies pointed out an enhanced release of reactive and legacy phosphorus in summer periods under low DO conditions for both mixed and single land-use streams [56,85,87]. Such a process can become an important feedback mechanism. The released nutrients sustain aquatic plant biomass and favor increased organic matter settling and decomposition and, thus, oxygen depletion as previously mentioned, especially if P is a limiting factor or other sources of reactive phosphorus (e.g., urban effluents) are removed or controlled. Notably, in the investigated catchment, the push towards green transition solutions for more resource and energy-efficient water treatment systems may result in centralized systems, with the suppression of urban effluent outlets [96]. While the benefit of such a solution is obvious in terms of overall environmental impact (e.g., carbon reduction and energy savings), is evident, it may come with trade-offs for the receiving waters (impacting biodiversity) that should be considered holistically, e.g., enhanced residence time due to flow reduction, lower depth, reduced thermal inertia, and more particle settling triggering some of the impairment mechanisms previously described.

Future Model Development and Data Needs
The results from the hydrological model were generally simulated well; however, some of the deviations observed in the period of extreme flows could not be entirely captured. The verification period for the three reaches comprising Usserød Stream finished 10 October 2019 (vertical dashed line in Figure 5), when the simulation started to deviate from the observations significantly. Floodplain inundation (flooding likely associated with a reduced channel capacity and high flows) strongly altered the flow regime and could not be handled with the current model structure (1-D system). Secondly, the introduction of a leaking term in the most downstream reach (see Section 3.3), based on local water abstraction data, was deemed necessary to improve the results during the period of low flows in summer.
Although a local change in geology has not been considered (we assume the calibrated parameters were similar between all reaches), this correction highlights the dynamic interaction between both groundwater and surface water especially important in small stream systems [97,98]. This interaction is challenging to handle in hydrological modeling, though critical in lowland catchments, and should be addressed to better understand the flow dynamics in low flow periods (see, e.g., the lump formulation and the coupling between surface water and groundwater using a lump hydrological model formulation in [99]). Finally, the period of overestimation for the flow at the beginning of the autumn period in the Donse Tributary (in connection with a possible water impoundment), the deviation between simulated and measured CSOs or the use of some time series (e.g., WWTP effluent; lake sluice data at the most upstream point) underline the importance of anthropogenic features affecting the hydrology and the required knowledge of their active control and stakeholder engagement for improved modeling, as also pointed out in [26].
The simulations for water quality in the stream are currently limited to the in-stream process. Therefore, a recommended next step would be to address the land source dynamics (e.g., fertilizer and nutrient applications, dissolved organic matter source) and the multiple flow contributions. The simulation of DO concentrations at the sub-daily time step highlights the need for improved characterization of the influence of the aquatic plant biomass, both in terms of spatial coverage but also in terms of autotroph group dynamics, to fully capture the overall dynamics of DO in small peri-urban streams. This task, however, may need to rely on more spatially detailed and high-frequency data, very often not available [81,100].
Furthermore, the observed time shift for the diel oscillation cycle points towards potential uncertainty in the reaeration formulation. This uncertainty is especially high for shallow streams but is inherent to any water quality model. It has been observed in [101] that errors between measured (gas tracer studies) and estimated reaeration coefficient can be up to 30-50%, and even >100% depending on the used formulation and depth/water velocity. The aquatic plant biomass will influence the hydrological model and should ideally have a feedback effect on the hydraulic roughness relationship currently implemented, in addition to their potential to force the settling of suspended sediments. Finally, the deviation in terms of orthophosphate concentration in summertime highlights the potential need for a more detailed process description of the streambed compartment in terms of nutrient cycling and role as a pollutant-bound stock, but also as the scene for complex heterotrophic respiration processes.

Conclusions
In this study, we developed an integrated (water quantity and quality) model using a system dynamics approach to investigate the variations and interactions between the flow and physico-chemical parameters (stream temperature, dissolved oxygen, nutrients, and chlorophyll-a) in a peri-urban stream on a daily basis. To our knowledge, this is the first SD model investigating stream hydrology and water quality within a mixed land-use catchment context.

•
In terms of hydrology, the model performs as well as other integrated models (see, for instance, the flow simulations results in [82]) with performance indicators, e.g., NSE/RSR, ranging from very good to satisfactory for all reaches in the verification period for both flow and depth variation. It also gives satisfactory results in terms of physico-chemical conditions (stream temperature and dissolved oxygen, NSE/RMSE performance indicators). Notably, the model combined with a rich dataset highlighted the very dynamic contribution of flows stemming from urban and more pervious areas that will greatly influence water quality: inflow of labile dissolved organic matter from the numerous flow pathways triggered in periods of high flow, followed and/or combined with settling and decomposition of algae in periods of low flow may fuel an important heterotrophic activity leading to low dissolved oxygen levels. • The developed model allows the use of probabilistic Monte Carlo simulations to account for the high temporal variability and uncertainty inherent to water quality parameters. Notably, with respect to nutrients, it showed a potential preference for ammonium uptake by macrophytes compared to nitrate. The coupled investigation between simulation and measurements also indicates a potential for remobilization of phosphorus from the sediment in peri-urban streams, which may constitute a critical source for nutrients should other sources (e.g., WWTP effluent) be removed or reduced. Such a process should consequently be accounted for in water quality modeling.

•
The use of the model in combination with a rich dataset demonstrates, for the first time, that water quality impairments could be further exacerbated by the implementation of a green transition solution (here, rerouting the urban effluents out of the stream system). The sustainability of these solutions should therefore be holistically evaluated to fulfill multi-objective strategies and limit adverse environmental trade-offs.
Overall, the integrated model was valuable in uncovering key insights into the dynamics of a peri-urban stream and, in combination with measured data, revealed that these types of lowland stream catchments may have a high potential to be impacted by green transition solutions. We expect the model's flexibility, simple structure, and use of system dynamics will constitute a strong advantage for additional stakeholder engagement activities in this catchment and facilitate its ease of transferability and application in others.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10 .3390/w13213002/s1, Figure S1: Example for an SD module developed for simulating stream temperature, Figure S2: Sewer network extent and type for the Usserød catchment, Figure S3: Locations for documented water abstraction, separate systems outlets and combined sewer overflows for the Usserød catchment, Figure S4: Imperviousness distribution data for the Usserød catchment, Figure  S5: Scatterplot for the estimated Manning's number vs. flow rate measured at the outlet of the three modelled reaches in the catchment, Figure S6: Linear regression model for the water temperature estimation of the natural water and urban system components, Figure S7: Flow simulation results (blue line) compared to continuous monitoring station measurements (red line) for all three Usserød stream reach outlets, Figure S8: Water depth simulation results (blue line) compared to continuous monitoring station measurements (red line) for all three Usserød stream reach outlets, Figure S9: Combined Sewer Overflow events in Ådalsvej reach. Simulated events (green) are compared to measured events (red). Figure S10: Flow simulation results (blue line) compared to continuous monitoring station measurements (red line) at two points, Fredtoften and Brønsholmdalsvej, in the Donse Tributary, Figure S11: Simulation result for DO daily average concentration (blue line) compared to the online sensor data (red line) at Parrallelvej monitoring station, Figure S12: Simulation results (blue line) for NO3-N concentration in Ådalsvej and Nivemølle reach, compared to grab sample measurements (red circles), Figure S13: Simulation results (blue line) for NH4-N concentration in Ådalsvej and Nivemølle reach, compared to grab sample measurements (red circles), Figure S14: Simulation results (blue line) for PO4-P concentration in Ådalsvej and Nivemølle reach, compared to grab sample measurements (red circles), Figure S15: Simulation results (blue line) for suspended chl-a concentration in Ådalsvej and Nivemølle reach, compared to grab sample measurements (red circles), Figure S16: Forcing input concentrations at Grønnegade (red circle) and at the outlet of the tributary (purple diamond) for NO3-N; NH4-N; PO4-P and susp. Chl-a used for the simulations presented in Figures S12-S15 and Figure 7, Figure S17: Rainfall-Runoff model sensitivity testing (Nivemølle station), Figure S18: Sensitivity testing for the DO concentration simulation (Parallelvej Station), Table S1: Summary of available input and monitoring data for the simulated peri-urban catchment, Table S2: Overview of the model input data for the Usserød catchment, data type (time series or constant value), associated time range, and source for the baseline simulation, Table S3: Model performance rating for Nash-Sutcliffe Efficiency (NSE), Percent BIAS (PBIAS) and RMSEobservations standard deviation ratio (RSR), Table S4: Model parameters, all reaches, Table S5: Parameter probability distributions used for the MC simulation for dissolved oxygen concentration shown in Figure 6, Table S6 [103], Figure S4 f cs -Combined system fraction impervious areas Input [103], Figure S4 SS    Additional details regarding the input data (catchment and reaches, as well as physicochemical conditions of the different flow components) can be found in Table S4. The natural RR component of the streamflow is calculated using a lump model structure inspired by [63], with an additional fast runoff component included. It consists of 3 main stocks (or reservoirs): soil moisture, surface runoff (overland flow), and groundwater, which contribute to the hydrological response of the investigated catchment (see also Figure 2).
The first stock represents the soil moisture and is controlled by the following water balance equation with all parameters expressed in [mm/day]: where I is the infiltration, AET is the actual evapotranspiration, Perc is the flow towards the groundwater reservoir, Q ui , Q li represent an upper and lower interflow, and Q runoff is an overland flow ( Figure 2). The infiltration I is dependent on the saturation degree of the soil, defined by: where P [mm/day] is the precipitation, f imp [-] is the aggregated imperviousness of the sub-catchment, and W [-] an index representing the "wetness" based on the estimated soil moisture, SM, following the non-linear relationship: where FC [mm] represent a field capacity, and β a calibrated index.
The aggregated imperviousness f imp factor is estimated as a weighted mean value per sub-catchment based on spatial analysis of imperviousness per cadastral units (0: natural catchment, to 1: fully impervious catchment).
The part of precipitation not infiltrating corresponding to a saturated soil moisture stock, P·W 1 − f imp , is routed to the stream as an overland flow Q runoff via a runoff stock ROFF with a relatively short time constant: where RTC [days] is a time constant associated to the ROFF stock.
The evapotranspiration term AET is dependent on the degree of saturation of the soil and is controlled by a moisture ratio defined by SM LPET ; then, evapotranspiration is varying following a piecewise linear function: where LPET is a calibrated moisture threshold [mm], and PET is potential evapotranspiration [mm·day −1 ]. PET is estimated in our model using the following formulation [109]: where Re is the daily extraterrestrial radiation [MJ·m −2 ·day −1 ], T a is the air temperature [ • C], ρ is the mass density of water, and γ the latent heat of vaporization [MJ·kg −1 ]. Q ui , Q li are two flows out of the soil moisture stock corresponding to shallow groundwater flow, typically interflows. Q ui is only active when the soil moisture is above a certain threshold, according to: where UITC, LITC are the stock time constants used for calibration [days], and SMT [mm] a moisture threshold for activation of the upper interflow. Finally, percolation Perc is a downwards water movement from the soil moisture stock towards the groundwater reservoir GW, which can be reduced by potential deep groundwater abstractions GWAR, and contributes to the stream Baseflow: where PRTC, BTC are the reservoir time constants used for calibration [days]. The coefficient contrib GWAR accounts for the fact that abstraction may occur directly from this stock (contrib GWAR = 1), but could also be impacted only by an abstraction from a deeper aquifer, one that may not be bounded by the sub-catchment delineation (contrib GWAR < 1). All these flow components (Equations (A5), (A9), (A10), and (A13) are scaled by the sub-catchment area and summed within a stream water reservoir STREAM to obtain the streamflow [m 3  The part of precipitation P·f imp not naturally drained by the sub-catchment, as described in Section A.1.1, is routed to the stream through the urban compartment, either by separate or combined systems. We assume that the drainage network follows the subcatchment delineation, which is deemed reasonable if the network is drained by gravity. Thus, any precipitation collected by the urban compartment within the sub-catchment will be discharged into its associated reach. The general water balance equation for the urban compartment, expressed in [m 3 /day] is: where f cs is the fraction of combined sewer network (=1 if only combined sewer network is implemented, 0 if fully separated) and SSTC [days] is a reservoir time constant. CSOs are simulated using a simple reservoir equivalent to a basin with two non-linear outflows, collecting water prior to the transfer, treatment, and discharge from the WWTP, inspired from [31], and was used in this paper. These flows consist of an overflow CSO activated when the water volume V reaches a maximum volume threshold V threshold and varies linearly with the excess volume of water V − V threshold . Simultaneously, the basin is emptied by a flow Q out,basin proportional to the volume of water V in this basin, up to a maximum flow capacity Q max . Only the overflow component CSO is considered here, the other one Q out,basin having already been accounted for by the WWTP effluent time series. The water balance for the basin is: And Q in,basin = A·P·f imp ·f cs + Q dry (A21) where Q dry [m 3 /day] is the daily average dry flow of the WWTP in this combined system network, Q max [m 3 /day] is the maximum flow capacity of the combined system, CSTC is a time constant for the drainage outflow, V threshold m 3 is the maximum volume capacity for the basin, and CSOTC is a time constant for the overflow [day]. Flow components (Equations (A19) and (A23)) and WWTP effluent time series are routed through the reach using a 3rd-order delay function (equivalent to a 3rd-order Nashcascade model) and added to the pervious component (Equation (A16)). A final possible flow transfer from the stream to the groundwater (loosing stream section) is accounted for by use of a leaking term Loss assumed linearly driven by the water abstraction via a coefficient GSI and possibly reducing the overall streamflow:

.3. Equations for Stream Water Depth and Velocity
The water depth d and average stream water velocity U for a given reach is estimated using Manning's equation and the computed stream flow. We assume a rectangular crosssection and a relatively low depth compared to the width W of the channel (d b), resulting in the following depth and flow velocity U estimates: where n is the Manning's coefficient, Q out the stream flow [m 3 /s], s the average channel slope [-], and b the average channel width [m]. The Manning's coefficient n represents the resistance or roughness of the channel to the flow and is a critical parameter for the depth estimation with a significant impact also on water quality (especially dissolved oxygen). We use a dynamic flow-dependent Manning's coefficient in the model formulation, following observations in lowland streams with high macrophyte densities, where Manning's coefficient decreases at higher discharges when submerged plants are flexible and bend with the flow [53]: where α n , β n are regression parameters. The stream temperature T w is predominantly influenced by the temperature of the different flow components discharging to the reach and by the heat exchange at the interface stream/air [110]. We use a lump heat balance over the water volume V in a reach according to the following equation [54]: here ρ is the water density, C p is the water-specific heat capacity [MJ·kg −1 · • C −1 ], A is the area of river/atmosphere interface [m 2 ], H is the net atmospheric flux [MJ·m 2 ·day −1 ], T i , Q i are the water temperature and inflows to the reach, respectively, and Q out is the outflow of the reach. At a daily time step, the thermal exchange with the streambed is considered as having a minor influence, and the net atmospheric flux H can be estimated using the concept of equilibrium temperature T e . Notably, H is proportional to the temperature difference between equilibrium and water temperature [55]: where K is a coefficient of heat transfer [MJ·m −2 · • C −1 ·day −1 ]. We assume a linear correlation between equilibrium temperature and air temperature, deemed appropriate for temperate regions, leading to the piecewise function for the equilibrium temperature estimation [55]: where T a [ • C] is the air temperature.

Appendix A.3.2. Dissolved Oxygen
Dissolved oxygen in the stream is affected by carbonaceous biochemical oxygen demand (BOD), nitrification, background sediment oxygen demand processes, and daily oscillations stemming from photosynthesis/autotrophic respiration activities while being simultaneously replenished by reaeration driven by the DO deficit to saturation level: The reaeration rea process is driven by the hydrological conditions, and the oxygen deficit with respect to oxygen saturation in the stream: where k a day −1 is the reaeration rate constant and DO sat [mg/L] the oxygen saturation concentration. The reaeration rate constant k a is essentially driven by the hydrological conditions (output from the hydrological model) and the stream water temperature. We use the Owens-Gibbs formulation in our model, suitable for relatively low flow velocities (<0.5 m·s −1 ) and shallow streams (<0.75 m): where U is the stream flow velocity [m/s], d water depth, and θ a temperature correction factor θ a ∼ = 1.024 [60]. The oxygen saturation concentration DO sat is dependent on the stream water temperature, altitude Alt, and salinity Sal, as described in [73]: The carbonaceous biochemical oxygen demand (BOD) induces an oxygen depletion caused by the aerobic degradation of organic matter. BOD is described in our model by a Streeter-Phelps-Shishkin equation system, according to [57]. This latter component introduces feedback between the degradation rate and the concentration of dissolved oxygen (DO), i.e., the rate of degradation is limited by how much dissolved oxygen is available, preventing the occurrence of negative DO values (a common problem when using the Streeter-Phelps formulation, alone): where k 20 d day −1 is the BOD decomposition rate at 20 • C, BOD [mg/L] is the carbonaceous biochemical oxygen demand, and a temperature correction factor θ BOD ∼ = 1.047 [62].
Currently, BOD is represented as a single stock based on the general Equation (3), with a decrease in BOD from the aerobic degradation of organic matter, as well as a possible settling for large organic particles: where i refers to the different flow component entering the reach (Equations (A16), (A19) and (A23) + WWTP effluent), k s d is the BOD degradation rate as previously defined, and k s is the settling rate day −1 dependent on a particle settling velocity V s,BOD m·day −1 and water depth d [60]: The nitrification process corresponds to the double step oxidation of ammonium resulting in the combined formation of nitrate and oxygen consumption: where k 20 nit is the nitrification rate day −1 , NH4 is the ammonium concentration [mgN/L] and r on is the stoichiometric ratio of mass oxygen consumed per mass nitrogen (=4.57 g·gN −1 ).
The nitrification rate k nit is highly dependent on environmental conditions co-limited primarily by temperature, pH and the oxygen concentration [60]: where k 20 nit is the nitrification rate at 20 • C, pH nit_corr is a correction factor dependent on the pH ( Figure A1, following [22]), f DO_nit_corr accounts for the influence of the DO concentration, and θ is a temperature correction factor θ nit ∼ = 1.085. where k nit 20 is the nitrification rate at 20 °C , pH nit_corr is a correction factor depen the pH ( Figure A1, following [22]), f DO_nit_corr accounts for the influence of the centration, and θ is a temperature correction factor θ nit ≅ 1.085. The nutrient concentrations currently simulated in the model include inor trogen (in terms of nitrate and ammonium) and dissolved reactive phosphorou phosphate).

Ammonium
The mass balance for ammonium in the stream water is: The nutrient concentrations currently simulated in the model include inorganic nitrogen (in terms of nitrate and ammonium) and dissolved reactive phosphorous (orthophosphate).

Ammonium
The mass balance for ammonium in the stream water is: Overall, nutrient assimilation by the aquatic plant biomass is estimated by computing a net carbon assimilation from the estimated daily net photosynthesis rate [111] and from the mass stoichiometric ratios C:N and C:P using the Redfield ratio. Nitrogen will be assimilated from both the ammonium and nitrate present in the stream water. We assume a preference ratio F NH4,pref for ammonium compared to nitrate for the assimilation, estimated by the following formula [60]: where k NH4,pref is a half-saturation constant for the ammonium preference [mgN/L]. Furthermore, NH4 can also be found in the un-ionized form of ammonia NH3 in stream water, which cannot be assimilated directly by plant biomass and is therefore considered not available for assimilation. Ammonia concentration NH3 is estimated by the acid-based equation system: where pKa is the acid-dissociation constant for ammonia (and is temperature-dependent).

Nitrate
The mass balance for nitrate in stream water is: where denitrification is modeled as a first-order removal dependent on temperature and DO concentration [62]: where k 20 denit is the denitrification rate defined at 20 • C [day −1 ], θ denit is a temperature correction factor θ denit ∼ = 1.047, and k s,denit_lim [mg/L] is the half-saturation constant.

Ortho-Phosphate
The mass balance for soluble reactive phosphorus in stream water is: The assimilation process for orthophosphate is similar to nitrogen, described in Appendix A.3.5. where SOD 20 o [mg/L] is a constant average value for sediment oxygen demand at 20 • C, d [m] is the stream depth, DO [mg/L] is the dissolved oxygen concentration, k s,olim [mg/L] is the half-saturation value for the oxygen dependency, and θ SOD is a temperature correction factor θ SOD ∼ = 1.065 [62]. The SOD enhanced,mac term is detailed further in the following section (Equation (A72)). The stream plant biomass is split into two main stocks (reservoirs) in this version of the model. One represents the aquatic plants attached or fixed in the stream (e.g., macrophytes and periphyton), and the second represents algae suspended in the water column and thus transported (e.g., phytoplankton). These stocks aggregate the plant biomass without specific species distinction. The biomass for both groups is expressed in terms of chlorophylla as a proxy and follows a similar mass balance, except for the transport and potential settling of suspended algae. Any loss by predation, grazing, sloughing, or scouring (for benthic algae and macrophytes) is currently neglected for both stocks.
For aquatic plants attached to the streambed (macrophytes), it is assumed that no transport of plant matter occurs and that the decay and decomposition of organic matter will be accounted for in the SOD pool; then the mass balance for plant biomass can be calculated by: ∂(V·chla mac ) ∂t = k growth,mac − k resp,mac ·chla mac (A55) where chla mac is the fixed biomass chlorophyll-a concentration [mg/L] k growth,mac is the associated growth rate, and k resp,mac is the a mortality/respiration/excretion rate. The mass balance for phytoplankton in a given reach is: ∂(V·chla sus ) ∂t = ∑ i Q i ·chla sus,i − Q out ·chla sus,i + k growth,sus − k resp,sus − k set,sus ·chla sus,i (A56) where chla sus is the suspended chlorophyll-a concentration [mg/L], Q i , Q out are the inflow and outflow from the reach, respectively, chla sus,i is the suspended chlorophyll-a concentration flowing into the reach, k growth,sus [day −1 ] is the growth rate, k resp,sus [day −1 ] is a combined mortality/respiration/excretion rate, and k set,sus [day −1 ] is a settling rate. The settling term k set is strongly dependent on the shape, size, and hydrological conditions, but will be simply defined in this model by: where V set,sus [m·day −1 ] is a constant settling velocity and d [m] is the water depth. As previously mentioned, the term k resp,sus/benthic (biomass stock-dependent) accounts for all losses affecting plant growth, i.e., plant maintenance respiration, excretion, and decay. This process is temperature-dependent and defined by: where k 20 r is a respiration/excretion rate at the reference temperature of 20 degrees, and f θ is a temperature dependency function (see below).
The growth rate k growth, sus/mac (biomass stock-dependent) is based on an optimal growth rate of the plant biomass, and limited by some environmental conditions in terms of light, nutrients and temperature: k growth = G max ·f θ ·f n,p ·f l (A59) where G max is an optimal growth rate (stock dependent), f l [-] is a function to account for the light dependency, f n,p to account for the nutrient concentration dependency, and f θ [-] is a temperature correction function. The temperature dependency function, f θ , is defined as a skewed normal distribution to account for the optimal growth of plant biomass at a specific temperature, and suboptimal conditions at higher or lower temperatures [62]: With T x = T min if T w < T opt and T x = T max if T w > T opt (A60) where T opt is the optimum temperature corresponding to the optimal growth rate, and T min , T max are the minimum and maximum extreme temperatures at which growth ceases.
The nutrient concentration dependency f n,p follows a Michaelis-Menten formulation. It is assumed that only inorganic nitrogen (in the form of ammonium or nitrate) and orthophosphate are limiting factors: f n,p = min PO4 k s,P_lim + PO4 , N k s,N_lim + N (A61) where k s,P_lim and k s,N_lim are half-saturation constants [mgP·L −1 or mgN·L −1 ].
The light correction f l accounts for the photoinhibition of growth at high light levels, as well as the light attenuation over depth due to particles and resulting turbidity. This correction factor is integrated over time and depth to get a mean daily correction for a well-mixed stream reach [60]: where L photo [-] is the photoperiod, d [m] is the stream water depth, γ [m −1 ] is the light attenuation coefficient, and α 0 and α 1 are functions of the light radiation and defined by: where I opt [ly·day −1 ] is the optimal light radiation for plant growth (dependent on plant species but assumed here simply dependent on the autotroph group) and I A is the average light radiation over daylight hours. The light attenuation coefficient γ l is function of the turbidity of the stream water, caused by all non-biomass matter that may be present in the water γ l,background (referred here as background), as well as any suspended algae, and/or self-shadowing effects from any fixed macrophytes. This coefficient is defined in our model by the combining the formulation by [60,63]: γ l = γ l,background +γ l,bio_sus × Chla sus + 10 0.57 log 10 (Chla mac −0.95) where γ l,background [m −1 ] is the light attenuation coefficient due to non-biomass materials, and γ l,bio_sus [m −1 ·L·mgchla −1 ] is the attenuation coefficient factor for the suspended algae, Chla mac is the concentration of water plant (macrophyte). I A , in the absence of data, is estimated from the mean extraterrestrial radiation and successive attenuation terms: where R e is the mean extraterrestrial radiation over the daylight hours [converted in ly·day −1 ], ∅ atm is a mean atmospheric absorption under cloud-free conditions [-], ∅ cloud represents the cloud absorption, ∅ par is the ratio between global horizontal radiation and photosynthetically active radiation, [-] is a reduction factor to account for any shadow effects at ground level and reach dependent, e.g., riparian vegetation, and ∅ reflection is a reduction factor to account for the light reflection at the stream surface. ∅ cloud is assessed using mean cloudiness data on a daily basis [112]: where N is the cloudiness data [okta].

Photosynthesis and Autotrophic Respiration
The daily mean gross photosynthesis rate P [in mg·L −1 ·day −1 ] is estimated by the following formula [60]: P = r oa · ∑ sus+ mac G max,sus/mac ·f l ·f θ ·chla sus/mac (A68) where r oa [mg·µgChla −1 ] is the oxygen yield per unit biomass and correction terms defined in Equations (A60) and (A62). The associated autotrophic respiration rate AR [mg·L −1 ·day −1 ] is estimated as a fraction of the daily mean gross photosynthesis rate [64]: where AR ratio = 0.44.
We emulate the diel variation of dissolved oxygen resulting from photosynthesis P(t) from the daily mean gross photosynthesis rate using the photoperiod length and the idealized half sinus profile for available light during the day [58]: P(t) = P MAX (cos(−2πt) + 2L photo − 1) (A70) and P MAX = P· π 2L photo (A71) Many streams are heterotrophic ecosystems, i.e., the gross primary production, GPP, from the autotrophic biomass is less than the overall ecosystem respiration, ER, (sum of autotrophic respiration and heterotrophic respiration, HR, from organic matter decomposition). Notably, such conditions have been documented in streams with important macrophyte coverage, enhancing the settling of fine particles fueling heterotrophic respiration and resulting in oxygen consumption (see, for instance, [56,83]). We defined the SOD enhanced,mac corresponding to this enhanced heterotrophic respiration to account for this effect, with temperature and oxygen limitations, as follows: SOD enhanced,mac = P/EM ratio − AR DO DO + k s,o_lim θ SOD where EM ratio [-] is the ecosystem metabolism ratio, P is the daily mean gross photosynthesis rate (Equation (A68)), and AR is the autotrophic respiration term (Equation (A69)); the other terms representing the temperature and oxygen limitation functions have already been defined (Equation (A54)).