Next Article in Journal
Sensitivity Analysis of the Catalytic Ozonation under Different Kinetic Modeling Approaches in the Diclofenac Degradation
Next Article in Special Issue
Impact of Sediment Layer on Longitudinal Dispersion in Sewer Systems
Previous Article in Journal
Review of Global Interest and Developments in the Research on Aquifer Recharge and Climate Change: A Bibliometric Approach
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Department of Environmental Engineering, Technical University of Denmark, Bygningstorvet 115, 2800 Kgs. Lyngby, Denmark
2
Swedish Meteorological and Hydrological Institute, Folkborgsvägen 17, SE-601 76 Norrköping, Sweden
*
Author to whom correspondence should be addressed.
Water 2021, 13(21), 3002; https://doi.org/10.3390/w13213002
Submission received: 14 September 2021 / Revised: 3 October 2021 / Accepted: 13 October 2021 / Published: 26 October 2021
(This article belongs to the Special Issue Water Quality Modeling and Monitoring)

Abstract

:
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.

1. 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 pathways and related response times, with faster drainage from impervious urban areas compared to more natural ones, results in high spatial and temporal variability at different scales, as reported in [13,14,15]. Moreover, measurements taken in these types of catchments are still generally insufficient to capture—spatially or temporally—the driving processes leading to these variations, which are needed to ensure effective mitigation strategies and sustainable decision-making [16]. Management decisions regarding stream flow and quality within the peri-urban landscape are thus extremely challenging. Integrated water quality models, i.e., considering the mutual interaction of flow and quality with potential for integration with the broader socio-economical-ecological system, are therefore essential tools to support water resources management and facilitate decision-making, as advocated in a recent review of water quality models ([17]).
Many hydrological and water quality models have been developed over the years. These models can be classified by increasing degree of complexity from statistical (e.g., regression-based) to more mechanistic models (physical description of processes), and spatially from lumped to physically distributed, chosen depending on the application purpose. Widely applied water quality models (based on a number of published studies over the last 20 years [18]) include the Soil and Water Assessment Tool (SWAT) [19], WASP [20], QUAL2E/K [21], AQUATOX [22], and the MIKE series [23]. However, their use to simultaneously address the hydrology and water quality of streams within peri-urban contexts is challenging. Reasons for this may lie in the reliance on a hydrodynamic model not yet fully adapted to the peri-urban context [24], the large data requirement, and the parametrization level necessary when using spatially distributed models. This can also cause long simulation times [25] and/or difficulty in accounting for the inherent uncertainty and high variability of input data [26].
Consequently, specific research-based models have been developed for investigating peri-urban streams. However, these models tend to be fragmented in sub-disciplines, e.g., focusing on the hydrology [27] or targeting a specific water quality aspect, e.g., identifying source contribution and pathways for pesticide pollution [28]. While they certainly have the capability to perform integrated simulations or to be further developed by being coupled together, their interfacing (data format, calculation time step) may not be straightforward. Machine learning techniques have become extremely popular for water quality (for instance, for the prediction of pH, suspended solids, and ammonium concentration using different algorithms in [29]), possibly coupled to other models (e.g., application of artificial neural networks coupled with SWAT in [30]) and could therefore be applied in peri-urban settings. Nevertheless, these models also rely heavily on data and may be limited in terms of direct (cause-effect) links to process understanding and thus decision-support capabilities, as experienced by [31]. Finally, and to this day, very few models offer the possibility and flexibility for stakeholder engagement, the inclusion of local knowledge, and integration capabilities at a higher level for socio-economical-ecological consideration, which is deemed necessary for better and more acceptable stream water management [26].
A promising solution that can fill the gap in terms of watershed modeling approaches, capable of offering users a holistic and reliable understanding of the systems involved, including uncertainties related to both linear and nonlinear dependencies within and across sub-systems, is system dynamics (SD) simulation. The approach is based on the implementation of an interconnected system of flows and stocks, whose structure gives rise to the dynamic behavior [32]. SD introduces flexibility lacking in many other methods, including faster model development and shorter simulation time, ability to simulate interactions and thus interdependencies between model sub-systems, and improved transparency for conceptual understanding of the system and thus communication of model results [33]. Moreover, the benefits of SD for water resources modeling are numerous, e.g., multidisciplinary aspects, causality analysis, stakeholder participation, as documented in [34,35], and it has been broadly applied to numerous environmental (e.g., climate change [36], waste management issues [37]) and water resource management issues at different spatial scales (see Section 2.1 below).
Although SD simulation continues to gain traction as a useful tool, there is still a gap with respect to studies utilizing this approach to provide holistic, integrated model insights within the field of water resources modeling and planning [34], and more specifically, within water quality modeling (see also Section 2.1). This study aims to fill this gap, presenting an integrated SD model for (1) water quantity and (2) water quality (i.e., physico-chemical conditions), which is used to (3) enhance our understanding of the highly dynamic processes governing water quality in a peri-urban catchment context.
Other SD models exploring water quality and quantity exist (e.g., [38]); however, these are limited in their exploration of the aforementioned dynamic processes. Modules have been developed, which can be implemented in a plug-and-play fashion, to handle the complexity related to both urban and agricultural/natural land-use features that may be present within a given stream’s reach, in terms of sources and pathways which may impact water quality. The model has been applied in a data-rich catchment to demonstrate its applicability. In this catchment, degraded water and ecological quality have been routinely observed. More specifically, important temporal variations, in terms of dissolved oxygen, temperature, flow, and nutrient discharge from both urban and agricultural areas are possible drivers of ecological degradation, although other stressors (e.g., habitat degradation) co-exist [39]. Comparison of simulation results with both available sensor and grab sampled data provide valuable insights into the processes affecting peri-urban stream quality and the importance of integrating water quantity and quality modeling to improve system understanding.

2. Methodology

2.1. System Dynamics Approach

The SD model has been developed using the visual object-oriented software Stella Architect [40]. SD models have been developed to address both water quantity and quality challenges, though typically separately. This includes models that simulate local hydrological processes [41] or capture the general water balance at catchment or regional scales (see, for instance, [42,43,44]). Specific applications for water quality include nutrient release from agriculture and urban areas with stakeholder participation [43,45,46], pollutant point discharge from raw sewage or bypass effluents [47], pathogens [48], salinity issues [49], contaminant spills [50], and contaminated groundwater transport including with discharge to a receiving stream [51]. To the authors’ knowledge, the model presented here is the first SD model to integrate stream water quantity and quality with the objective of simulating the major features present within a mixed land-use catchment context.

2.2. 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 physico-chemical parameters. Currently, dissolved oxygen (DO), stream temperature (temp), nitrate (NO3), ammonium/ammonia (NH4), soluble reactive phosphorus (PO4), 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.
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.

2.2.1. Hydrological Model

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 (fimp) (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:
Q out = Q in , routed + RR Streamflow j outflows
where Qout is the flow at the outlet of a reach; Qin,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).

2.2.2. 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 ]):
d ( V · C ) dt = Q in C in Q out C + S ( C , p )
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]:
ρ C p ( V · T w ) t = AH + ρ C p ( T i Q i   T w Q out )
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 [m2], 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, Equation (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.

2.3. 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 spatio-temporal 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.

3. Model Application

3.1. The Usserød Peri-Urban Catchment

The model was applied to a small low-land catchment (ca. 120 km2, 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].
Figure 2. Overview of the Usserød catchment with stream and tributary (a) general land use [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.
Figure 2. Overview of the Usserød catchment with stream and tributary (a) general land use [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.
Water 13 03002 g002

3.2. 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.

3.3. Model Set-Up and Calibration

3.3.1. 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 up-gradient 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.

3.3.2. 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 least-squares 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, 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.

3.3.3. 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].

4. 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).

4.1. 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 m3/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 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).
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 sub-catchment, 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.

4.2. Physico-Chemical and Water Quality Parameters

4.2.1. Stream Temperature

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.

4.2.2. Dissolved Oxygen Concentration

Dissolved oxygen concentrations in the stream are relatively well captured by the model with a coefficient of determination R2 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, Figure 7 and Figure 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 Sjæl 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.
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 (Figure 6a and Figure 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.

4.2.3. 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 NO3-N concentrations rule out a possible underestimation of the nitrification process, suggesting instead an assimilation preference for NH4-N (rather than NO3-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 NO3 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 (Figure 6a and Figure 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].

5. 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 (NH4-N, NO3-N, PO4-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].

5.1. 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 peri-urban 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.
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.

5.2. 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.

6. 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 RMSE-observations 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: Parameter probability distributions used for setting up the MC simulation regarding physico-chemical conditions and shown in Figure 7.

Author Contributions

Conceptualization, G.G.L. and U.S.M.; methodology, G.G.L., S.A.C., S.G., V.M., U.S.M.; software, G.G.L., S.A.C., S.G. and V.M.; validation, G.G.L.; formal analysis, G.G.L.; investigation, G.G.L., U.S.M., S.A.C., S.G. and V.M.; resources, U.S.M. and P.L.B.; data curation, G.G.L.; writing—original draft preparation, G.G.L.; writing—review and editing, All authors; visualization, G.G.L.; supervision, U.S.M. and P.L.B.; project administration, U.S.M.; funding acquisition, U.S.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded partially by a PhD grant from the Department of Environmental Engineering, Technical University of Denmark. The RIVERSCAPES project (funded by the Innovation Fund Denmark under file number 7048-00001B) financially supported the measurements (grab samples) used for verification.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The model equations are presented in full in Appendix A and the data presented in this study are available on request from the corresponding author.

Acknowledgments

The authors thank the water utility company NOVAFOS for information and all data regarding the WWTPs, drainage system, and water abstraction in the catchment.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Appendix A.1. List of Variables Describing the System Considered in Figure 2

Table A1. Catchment and reach.
Table A1. Catchment and reach.
ParametersUnitDescriptionData TypeValue
A m 2 Subcatchment areaInput
L m Reach lengthInputAll inputs are
b m Average stream reach widthInputreach dependent
d m Average stream reach depthCalculatedSee Table S4
s -Average stream reach slopeInput
V m 3 Water volume in reachCalculated
U m / s Stream flow velocityCalculated
Table A2. Hydrology—natural flow component.
Table A2. Hydrology—natural flow component.
ParametersUnitDescriptionData TypeValue
P mm · day 1 PrecipitationInput[73]
Ta °CAir temperatureInput[74]
Re J / ( m 2 · d ) Extraterrestrial radiationInput[102]
n S/m1/3Manning’s coefficientCalculated
SM mm Soil moisture reservoirCalculated
ROFF mm Overland flow reservoirCalculated
GW mm Groundwater reservoirCalculated
STREAM mm Stream reservoirCalculated
I mm · day 1 InfiltrationCalculated
AET mm · day 1 Actual EvapotranspirationCalculated
PET mm · day 1 Potential EvapotranspirationCalculated
Perc mm · day 1 Percolation Calculated
Q runoff mm · day 1 Overland flowCalculated
Q ui mm · day 1 Upper interflowCalculated
Q li mm · day 1 Lower interflowCalculated
Baseflow mm · day 1 BaseflowCalculated
Streamflow m 3 · day 1 Stream flowCalculated
FC mm Field capacity equivalentParameter
LPET mm Soil moisture threshold for Potential Evapotranspiration Parameter
β-Evapotranspiration reduction factorParameter
W -Wetness indexCalculated
SMT mm Soil moisture threshold—upper interflowParameter
RTC day Runoff time constantParameter
UITC day Upper interflow time constantParameter
LITC day Lower interflow time constantParameter
PRTC day Percolation time constantParameter
BTC day Baseflow time constantParameter
SDR day Stream discharge rateParameter
GWAR m 3 · day 1 Groundwater abstraction rateInput[103]
contrib GWAR -Scaling factor for GWARParameterReach dep. Table S4
GSI -Factor for loosing stream sectionParameterReach dep. Table S4
Table A3. Hydrology—urban component.
Table A3. Hydrology—urban component.
ParametersUnitDescriptionData TypeValue
Q wwtp m 3 · day 1 WWTP effluent flowInput[103]
f imp -Fraction of impervious areaInput[103], Figure S4
f cs -Combined system fraction impervious areasInput[103], Figure S4
SS mm Separated system reservoirCalculated
Q urban _ cs + Q dry m 3 · day 1 WWTP effluent flow in combined systemInput[103]
Q urban _ ss m 3 · day 1 Separated system flowInput[103]
SSTC day Separated system time constantParameter0.05 (calibrated)
CSTC day Drainage outflow time constant Parameter0.5 (calibrated)
CSOTC day Overflow combined system time constantParameter0.95 (calibrated)
Q dry ¯ m 3 · day 1 Average dry discharge WWTP (black grey water)Parameter[103]
Q MAX m 3 · day 1 Max flow capacity towards WWTP (combined system)Parameter50,000 (Calibrated)
V threshold m 3 Max volume capacity—retention basin (combined system)Parameter10,000 (Calibrated)
Table A4. Temperature.
Table A4. Temperature.
ParametersUnitDescriptionData TypeValue
Tw °CStream water temperatureCalculated
Te °CEquilibrium temperatureCalculated
ρ kg / m 3 Mass densityParameter1000
Cp J / ( ° C · kg ) Specific heat capacity waterParameter4182
K J / ( m 2 · ° C · kg · day ) Heat transfert coefficientParameter1.9 × 106 [104]
Table A5. Oxygen and nutrients.
Table A5. Oxygen and nutrients.
ParametersUnit *DescriptionData TypeValue
Alt m Average altitudeInput8 [71]
Sal ppt   ( g · L 1 ) SalinityInput0
DO mg · L 1 DO concentrationCalculated
DOsat mg · L 1 DO saturation concentrationCalculated
pH - Input
ka day 1 Reaeration rate Calculated
θa-Temperature correction factor reaerationParameter1.024 [60]
BOD mg · L 1 Carbonaceous Biochemical Oxygen DemandCalculated
k d s day 1 Degradation rate for organic matterCalculated
k d 20 day 1 Reference degradation rate for organic matter (20 °C)Parameter0.5 [60]
θBOD-Temperature correction factor BODParameter1.047 [62]
ks day 1 Settling rate, organic matterCalculated
V s , BOD m · day 1 Settling velocity, organic matterParameter0.1 [60]
NBOD mg · L 1 Nitrogenous Biochemical Oxygen Demand (nitrification)Calculated
r on g · gN 1 Stoichiometric ratio DO consumed during nitrificationParameter4.57 [60]
k nit day 1 Nitrification rateCalculated
k nit 20 day 1 Reference nitrification rate (20 °C)Parameter0.3 [62]
θ nit -Temperature correction factor nitrificationParameter1.085 [62]
pH nit _ corr -Correction factor. pH effect on nitrification Calculated
f DO _ nit _ corr -Correction factor. DO concentration on nitrification Calculated
NH 4 mgN · L 1 Ammonium concentrationCalculated
NO 3 mgN · L 1 Nitrate concentrationCalculated
PO 4 mgP · L 1 Orthophosphate concentrationCalculated
F NH 4 , pref -Preference ratio for nitrogen assimilationCalculated
k NH 4 , pref mgN · L 1 Half saturation constant for ammonium preferenceParameter0.05 [60]
k denit day 1 Denitrification rateCalculated
k denit 20 day 1 Reference denitrification rate (20 °C)Parameter0.1 [62]
θ denit -Temperature correction factor denitrificationParameter1.045 [62]
k s , denit _ lim mg · L 1 Half saturation constant oxygen limitationParameter0.1 [62]
SOD mg · L 1 · day 1 Sediment Oxygen DemandCalculated
SOD BG g · m 2 · day 1 Background value Sediment Oxygen DemandCalculated
SOD enhanced , mac mg · L 1 · day 1 Enhanced heterotrophic respiration, macrophyte effectCalculated
SOD o 20 g · m 2 · day 1 Reference Sediment Oxygen Demand (20 °C)Parameter1 [62]
θ SOD -Temperature correction factor Sediment Oxygen DemandParameter1.065 [62]
k s , o _ lim mg · L 1 Half saturation constant for SOD oxygen limitationParameter1.4 [62]
* All units of concentration, e.g.,   mg · L 1 corresponds to mg O 2 by default.
Table A6. Freshwater plant and algae biomass.
Table A6. Freshwater plant and algae biomass.
ParametersUnit *DescriptionData TypeValue
Chla sus μ gChla · L 1 Suspended algae concentrationCalculated
Chla mac μ gChla · L 1 Plant biomass concentration (macrophyte—epiphyton)Calculated
k growth , sus / mac day 1 Algae/plant biomass growth rateCalculated
k resp , sus / mac day 1 Mortality/respiration/excretion rateCalculated
k set , sus day 1 Settling rate suspended algaeCalculated
V set , sus m · day 1 Settling velocity suspended algaeParameter0.1 [62]
G max , sus day 1 Optimum growth rate under no limitation, suspended algae Parameter2 [62]
G max , mac day 1 Optimum growth rate (macrophyte) no limitation Parameter0.1 [105]
f θ -Temperature effect function on biomass growth rateCalculated
f l -Light effect function on biomass growth rateCalculated
f n , p -Nutrient effect function on biomass growth rateCalculated
k s , P _ lim mgP · L 1 Half saturation constant phosphorus limitation of biomass growthParameter0.025 [62]
k s , N _ lim mgN · L 1 Half saturation constant nitrogen limitation of biomass growthParameter0.0005 [62]
L photo -Photoperiod of the dayInput[102]
γ l , background m 1 Light attenuation coefficient—non biomass materials Parameter0.5 [106]
γ l , bio _ sus m 1 · L · mgchla 1 Light attenuation coefficient factor—suspended algae Parameter0.035 [60]
I A ly · day 1 Average light radiation over daylight hoursCalculated
I opt ly · day 1 Optimal light radiation for plant growthParameter200 [62]
R e ¯ ly · day 1 Mean extraterrestrial radiation over the daylight hoursInput
atm -Atmospheric absorption (cloud free)Parameter0.38 [107]
PAR -Ratio PAR/global horizontal radiationParameter0.47 [108]
shadow -Radiation attenuation—shadow effect at ground levelInput0.1 (estimate)
reflec -Radiation attenuation due to reflection on the stream surfaceParameter0.1 [106]
cloud -Radiation attenuation due to could coverCalculated
N okta Cloudiness dataInput[74]
roa mg · μ gChla 1 Oxygen production rate per plant biomassParameter1.5 [60,63]
P ¯ mg · L 1 · day 1 Mean daily photosynthesis rateCalculated
AR mg · L 1 · day 1 Autotrophic respiration rateCalculated
AR ratio -Autotrophic respiration ratioParameter0.45 [64]
P MAX mg · L 1 · day 1 Max photosynthesis rateCalculated
EM ratio -Ecosystem metabolism ratioParameter0.26 (Calibrated)
* All units of concentration, e.g.,   mg · L 1 corresponds to mg O 2 by default.
Additional details regarding the input data (catchment and reaches, as well as physico-chemical conditions of the different flow components) can be found in Table S4.

Appendix A.2. Hydrological Model

Appendix A.2.1. Pervious Area Flow Component

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 ]:
dSM dt = I AET +   Perc +   Q ui +   Q li + Q runoff
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:
I   =   P ( 1 W ) ( 1 f imp )
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:
W = ( SM FC ) β
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:
dROFF dt = P · W ( 1 f imp ) Q runoff
Q runoff = 1 RTC · ROFF
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:
AET   = SM LPET PET if   SM LPET < 1
AET   =   PET if   SM LPET > 1
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]:
PET   = Re γ ρ · T a + 5 100
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:
Q li = SM · 1 LITC
Q ui = MIN ( 0 , ( SM SMT ) · 1 UITC )
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:
dGW dt = Perc Baseflow Abstraction
Perc = SM · 1 PRTC
Baseflow = GW · 1 BTC
Abstraction = contrib GWAR · GWAR
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 / day ].
dSTREAM dt = streamflow
streamflow = STREAM · 1 SDR
where SDR is a reservoir time constant [ days ] .

Appendix A.2.2. Impervious Areas Flow Component (Urban)

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 sub-catchment 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:
A · P · f imp = Q urban _ cs + Q dry + Q urban _ ss + [ CSOs ]
where A is the sub-catchment area [m2], P [ mm / day ] is the precipitation, f imp [-] is the aggregated imperviousness of the sub-catchment, Q urban _ cs + Q dry [ m 3 / day ] correspond to a WWTP effluent (wet and dry discharge), if relevant for the investigated reach,   Q urban _ ss is the flow component generated by precipitation and routed to the stream through the separate system. Finally, [ CSOs ] are possible overflow discharges from a combined sewer overflow system. Q urban _ cs + Q dry is currently input as a time series of WWTP effluent in this version of the model.
Separate systems flow Q urban _ ss   is estimated by a linear reservoir approach:
dSS dt = A · P · f imp ( 1 f cs ) Q urban _ ss .
Q urban _ ss = SS · 1 SSTC
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:
dV dt =   Q in , basin + Q out , basin + CSOs
And
Q in , basin = A · P · f imp · f cs + Q dry ¯
Q out , basin = MIN ( Q max , V · 1 CSTC )
CSO = MAX ( 0 , ( V V threshold ) · 1 CSOTC )
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 Nash-cascade 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:
loss = GSI · GWAR

Appendix A.2.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 cross-section 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:
d = ( n · Q out s 1 / 2 · b ) 3 / 5
U = 1 n · s 1 / 2 ·   d 2 / 3  
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]:
n = α n Q β n
where α n , β n are regression parameters.

Appendix A.3. Water Quality Model

Appendix A.3.1. Stream Temperature

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]:
ρ C p ( V · T w ) t = AH + ρ C p ( T i Q i   T w Q out )
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]:
H = K ( T e T w )
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]:
T e = 1 If   T a < 6
T e = 0.83 ·   T a 5 6 < T a < 30
T e = 25 if   T a > 30
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:
( V · DO ) t = i Q i DO i Q out DO + rea + carboneous   Biological   Oxygen   Demand nitrification Sediment   Demand + photosynthesis
The reaeration rea process is driven by the hydrological conditions, and the oxygen deficit with respect to oxygen saturation in the stream:
rea = k a · ( DO DO sat )
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):
k a = 5.3 U 0.67 d 1.85   θ a T w 20
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]:
DO sat = Alt eff · exp ( 139.3 + 1.57 × 10 5 T w + 6.64 × 10 7 T w 2 + 1 × 2410 8 T w 3 + 1 × 57.10 5 T w 4 Sal eff )
And
Alt eff = 100 ( 0.0035 × 3.28 × Alt ) 100
Sal eff = Sal × ( 0.01767 10.75 T w + 2.14 × 10 3 T w 2 )

Appendix A.3.3. Carbonaceous Biochemical (Biological) Oxygen Demand

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):
carboneous   Biological   Oxygen   Demand =   k d s · BOD
k d s = k d 20 DO DO sat · θ BOD Tw 20
where k d 20   [ 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:
( V · BOD ) t = i Q i BOD i Q out BOD ( k d s + k s ) BOD
where i refers to the different flow component entering the reach (Equations (A16), (A19) and (A23) + WWTP effluent), k d s 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]:
k s = V s , BOD d

Appendix A.3.4. Nitrification

The nitrification process corresponds to the double step oxidation of ammonium resulting in the combined formation of nitrate and oxygen consumption:
NBOD = k nit · NH 4 · r on
where k nit 20 is the nitrification rate [ day 1 ] , NH 4   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]:
k nit = k nit 20 · θ nit Tw 20 · pH nit _ corr · f DO _ nit _ corr
And
f DO _ nit _ corr = 1 e 0.6 DO
where   k nit 20 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.
Figure A1. Correction factor pH nit _ corr   for the effect of pH on the nitrification rate, from [22].
Figure A1. Correction factor pH nit _ corr   for the effect of pH on the nitrification rate, from [22].
Water 13 03002 g0a1

Appendix A.3.5. Macronutrients

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:
( V · NH 4 ) t = i Q i · NH 4 i Q out · NH 4 k nit · NH 4 assimilation
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 NH 4 , pref for ammonium compared to nitrate for the assimilation, estimated by the following formula [60]:
F NH 4 , pref = NH 4 k NH 4 , pref + NH 4
where k NH 4 , pref is a half-saturation constant for the ammonium preference [ mgN / L ] .
Furthermore, NH 4 can also be found in the un-ionized form of ammonia NH 3   in stream water, which cannot be assimilated directly by plant biomass and is therefore considered not available for assimilation. Ammonia concentration NH 3   is estimated by the acid-based equation system:
R = 1 1 + 10 pKa pH
And
NH 3 = R 1 R ·   NH 4
where pKa is the acid-dissociation constant for ammonia (and is temperature-dependent).

Nitrate

The mass balance for nitrate in stream water is:
( V · NO 3 ) t = i Q i · NO 3 i Q out · NO 3 + k nit · NH 4 assimilation denitrification
where denitrification is modeled as a first-order removal dependent on temperature and DO concentration [62]:
  denitrification = k denit 20 · θ denit Tw 20 k s , denit k s , denit _ lim + DO ·   [ NO 3 ]
where k denit 20 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:
( V · PO 4 ) t = i Q i · PO 4 i Q out · PO 4 assimilation
The assimilation process for orthophosphate is similar to nitrogen, described in Appendix A.3.5.

Appendix A.3.6. Sediment Oxygen Demand

In this model version, SOD is built around two main contributions. First, a background contribution SOD BG with temperature and DO concentration as co-factors. And second, an extra dynamic term accounts for the possible enhanced heterotrophic respiration SOD enhanced , mac   associated with the indirect effects of aquatic plant biomass occurring during the spring and summer seasons and trapping fine sediments [56]:
SOD = SOD BG + SOD enhanced , mac
SOD BG = SOD o 20 d × DO DO + k s , o _ lim θ SOD T w 20
where SOD o 20   [ 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)).

Appendix A.3.7. Aquatic Macrophytes, Algal Biomass and Autotrophic Metabolism

Aquatic Macrophytes and Algal Biomass

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 chlorophyll-a 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
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
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:
k set , sus = V set , sus d
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:
k resp = k r 20 ·   f θ
where k r 20 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
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]:
f θ = e ( 2.3 · ( T w T opt T x T opt ) 2 ) With   T x = T min   if   T w < T opt   and   T x = T max   if   T w > T opt
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 ( PO 4 k s , P _ lim + PO 4 , N k s , N _ lim + N )
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]:
f l = 2.718 × L photo γ l · d ( e α 1 e α 0 )
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:
α 0 = I A I opt
α 1 = I A I opt e γ l · d
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:
I A = R e ¯ · ( 1 atm ) · cloud · PAR · ( 1 shadow ) ·   ( 1 reflection )
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]:
cloud = ( 1 0.75 ( N / 8 ) 3.4 )
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
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]:
AR = AR ratio × P ¯
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 ) + 2 L photo 1 )
and
P MAX = P ¯ · π 2 L photo
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 T w 20
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)).

References

  1. Vörösmarty, C.J.; McIntyre, P.B.; Gessner, M.O.; Dudgeon, D.; Prusevich, A.; Green, P.; Glidden, S.; Bunn, S.E.; Sullivan, C.A.; Liermann, C.R.; et al. Global threats to human water security and river biodiversity. Nature 2010, 467, 555–561. [Google Scholar] [CrossRef]
  2. Dudgeon, D.; Arthington, A.H.; Gessner, M.O.; Kawabata, Z.-I.; Knowler, D.J.; Lévêque, C.; Naiman, R.J.; Prieur-Richard, A.-H.; Soto, D.; Stiassny, M.L.J.; et al. Freshwater biodiversity: Importance, threats, status and conservation challenges. Biol. Rev. Camb. Philos. Soc. 2006, 81, 163–182. [Google Scholar] [CrossRef]
  3. Reid, A.J.; Carlson, A.K.; Creed, I.F.; Eliason, E.J.; Gell, P.A.; Johnson, P.T.J.; Kidd, K.A.; MacCormack, T.J.; Olden, J.D.; Ormerod, S.J.; et al. Emerging threats and persistent conservation challenges for freshwater biodiversity. Biol. Rev. 2018, 94, 849–873. [Google Scholar] [CrossRef] [Green Version]
  4. Heal, K.V.; Bartosova, A.; Hipsey, M.R.; Chen, X.; Buytaert, W.; Li, H.-Y.; McGrane, S.J.; Gupta, A.B.; Cudennec, C. Water quality: The missing dimension of water in the water–energy–food nexus. Hydrol. Sci. J. 2020, 66, 745–758. [Google Scholar] [CrossRef]
  5. von Schiller, D.; Acuña, V.; Aristi, I.; Arroita, M.; Basaguren, A.; Bellin, A.; Boyero, L.; Butturini, A.; Ginebreda, A.; Kalogianni, E.; et al. River ecosystem processes: A synthesis of approaches, criteria of use and sensitivity to environmental stressors. Sci. Total Environ. 2017, 596–597, 465–480. [Google Scholar] [CrossRef]
  6. Green, P.A.; Vörösmarty, C.J.; Harrison, I.; Farrell, T.; Sáenz, L.; Fekete, B.M. Freshwater ecosystem services supporting humans: Pivoting from water crisis to water solutions. Glob. Environ. Chang. 2015, 34, 108–118. [Google Scholar] [CrossRef]
  7. Tickner, D.; Opperman, J.J.; Abell, R.; Acreman, M.; Arthington, A.H.; Bunn, S.E.; Cooke, S.J.; Dalton, J.; Darwall, W.; Edwards, G.; et al. Bending the Curve of Global Freshwater Biodiversity Loss: An Emergency Recovery Plan. Bioscience 2020, 70, 330–342. [Google Scholar] [CrossRef] [PubMed]
  8. Van Meter, K.; Thompson, S.E.; Basu, N.B. Human Impacts on Stream Hydrology and Water Quality; Elsevier Inc.: Amsterdam, The Netherlands, 2016; ISBN 9780124058903. [Google Scholar]
  9. Pinto, U.; Maheshwari, B.L. River health assessment in peri-urban landscapes: An application of multivariate analysis to identify the key variables. Water Res. 2011, 45, 3915–3924. [Google Scholar] [CrossRef]
  10. Meyer, J.L.; Paul, M.J.; Taulbee, W.K. Stream ecosystem function in urbanizing landscapes. J. N. Am. Benthol. Soc. 2005, 24, 602–612. [Google Scholar] [CrossRef]
  11. Marttila, H.; Lepistö, A.; Tolvanen, A.; Bechmann, M.; Kyllmar, K.; Juutinen, A.; Wenng, H.; Skarbøvik, E.; Futter, M.; Kortelainen, P.; et al. Potential impacts of a future Nordic bioeconomy on surface water quality. Ambio 2020, 49, 1722–1735. [Google Scholar] [CrossRef] [PubMed]
  12. Piorr, A.; Ravetz, J. Peri-Urbanisation in Europe. Towards European Policies to Sustain Urban-Rural Futures; Synthesis Report; University of Copenhagen: Copenhagen, Danmark, 2011. [Google Scholar]
  13. Aurélio, M.; Cruz, S.; De Azevedo, A.; Ricardo, G.; Roberto, J.; De Amorim, A.; Vinicius, P.; Vajapeyan, M.; Alexandre, C.; Garcia, B.; et al. Spatial and seasonal variability of the water quality characteristics of a river in Northeast Brazil. Environ. Earth Sci. 2019, 78, 1–11. [Google Scholar]
  14. Guo, D.; Lintern, A.; Angus Webb, J.; Ryu, D.; Bende-Michl, U.; Liu, S.; William Western, A. A data-based predictive model for spatiotemporal variability in stream water quality. Hydrol. Earth Syst. Sci. 2020, 24, 827–847. [Google Scholar] [CrossRef] [Green Version]
  15. Lintern, A.; Webb, J.A.; Ryu, D.; Liu, S.; Waters, D.; Leahy, P.; Bende-Michl, U.; Western, A.W. What Are the Key Catchment Characteristics Affecting Spatial Differences in Riverine Water Quality? Water Resour. Res. 2018, 54, 7252–7272. [Google Scholar] [CrossRef]
  16. Pinto, U.; Maheshwari, B. A framework for assessing river health in peri-urban landscapes. Ecohydrol. Hydrobiol. 2014, 14, 121–131. [Google Scholar] [CrossRef]
  17. Fu, B.; Merritt, W.S.; Croke, B.F.W.; Weber, T.R.; Jakeman, A.J. A review of catchment-scale water quality and erosion models and a synthesis of future prospects. Environ. Model. Softw. 2019, 114, 75–97. [Google Scholar] [CrossRef]
  18. Burigato Costa, C.M.; da Silva Marques, L.; Almeida, A.K.; Leite, I.R.; de Almeida, I.K. Applicability of water quality models around the world—a review. Environ. Sci. Pollut. Res. 2019, 26, 36141–36162. [Google Scholar] [CrossRef]
  19. Neitsch, S.L.; Arnold, J.G.; Kiniry, J.R.; Williams, J.R. SWAT: Soil & Water Assessment Tool Theoretical Documentation Version 2009. Texas Water Resour. Inst. 2011, 1–647. Available online: https://www.researchgate.net/publication/312462968_Soil_and_water_assessment_tool_theoretical_documentation (accessed on 10 February 2020).
  20. Di Toro, D.M.; Fitzpatrick, J.J.; Thomann, R.V. Documentation for Water Quality Analysis Simulation Program (WASP) and Model Ver-Ification Program (MVP); Environmental Research Laboratory, Office of Research and Development, US Environmental Protection Agency: Washington, DC, USA, 1983.
  21. Chapra, S.C.; Pelletier, G.; Tao, H. QUAL2K: A modeling framework for simulating river and stream water quality, Version 2.12. Doc. User Man. 2012, 97. Available online: www.ecs.umass.edu/cee/reckhow/courses/577/Qual2/Q2KDocv2_11b8%20v211.pdf (accessed on 10 February 2020).
  22. Park, R.A.; Clough, J.S.; Wellman, M.C. AQUATOX: Modeling environmental fate and ecological effects in aquatic ecosystems. Ecol. Modell. 2008, 213, 1–15. [Google Scholar] [CrossRef]
  23. DHI. ECOLAB User Guide; DHI: Hørsholm, Denmark, 2017. [Google Scholar]
  24. Braud, I.; Fletcher, T.D.; Andrieu, H. Hydrology of peri-urban catchments: Processes and modelling. J. Hydrol. 2013, 485, 1–4. [Google Scholar] [CrossRef] [Green Version]
  25. Thuy, T.; Keupers, I.; Willems, P. Environmental Modelling & Software Conceptual river water quality model with fl exible model structure. Environ. Model. Softw. 2018, 104, 102–117. [Google Scholar]
  26. Fu, B.; Horsburgh, J.S.; Jakeman, A.J.; Gualtieri, C.; Arnold, T.; Green, T.R.; Quinn, N.W.; Volk, M.; Hunt, R.J.; Vezzaro, L.; et al. Modeling water quality in watersheds: From here to the next generation. Water Resour. Res. 2020, 56, e2020WR027721. [Google Scholar] [CrossRef] [PubMed]
  27. Jankowfsky, S.; Branger, F.; Braud, I.; Rodriguez, F.; Debionne, S.; Viallet, P. Assessing anthropogenic influence on the hydrology of small peri-urban catchments: Development of the object-oriented PUMMA model by integrating urban and rural hydrological models. J. Hydrol. 2014, 517, 1056–1071. [Google Scholar] [CrossRef]
  28. Wittmer, I.K.; Bader, H.P.; Scheidegger, R.; Stamm, C. REXPO: A catchment model designed to understand and simulate the loss dynamics of plant protection products and biocides from agricultural and urban areas. J. Hydrol. 2016, 533, 486–514. [Google Scholar] [CrossRef]
  29. Najah Ahmed, A.; Binti Othman, F.; Abdulmohsin Afan, H.; Khaleel Ibrahim, R.; Ming Fai, C.; Shabbir Hossain, M.; Ehteram, M.; Elshafie, A. Machine learning methods for better water quality prediction. J. Hydrol. 2019, 578, 124084. [Google Scholar] [CrossRef]
  30. Noori, N.; Kalin, L.; Isik, S. Water quality prediction using SWAT-ANN coupled approach. J. Hydrol. 2020, 590, 125220. [Google Scholar] [CrossRef]
  31. Schmidt, L.; Heße, F.; Attinger, S.; Kumar, R. Challenges in Applying Machine Learning Models for Hydrological Inference: A Case Study for Flooding Events Across Germany. Water Resour. Res. 2020, 56, e2019WR025924. [Google Scholar] [CrossRef]
  32. Simonovic, S.P. Managing Water Resources: Methods and Tool for System Approach; Taylor & Francis: Abingdon, UK, 2012; Volume 9781849771. [Google Scholar]
  33. McKnight, U.S.; Finkel, M. A system dynamics model for the screening-level long-term assessment of human health risks at contaminated sites. Environ. Model. Softw. 2013, 40, 35–50. [Google Scholar] [CrossRef]
  34. Winz, I.; Brierley, G.; Trowsdale, S. The use of system dynamics simulation in water resources management. Water Resour. Manag. 2009, 23, 1301–1323. [Google Scholar] [CrossRef]
  35. Zomorodian, M.; Lai, S.H.; Homayounfar, M.; Ibrahim, S.; Fatemi, S.E.; El-Shafie, A. The state-of-the-art system dynamics application in integrated water resources modeling. J. Environ. Manag. 2018, 227, 294–304. [Google Scholar] [CrossRef] [PubMed]
  36. Randers, J.; Golüke, U.; Wenstøp, F.; Wenstøp, S. A user-friendly earth system model of low complexity: The ESCIMO system dynamics model of global warming towards 2100. Earth Syst. Dyn. 2016, 7, 831–850. [Google Scholar] [CrossRef] [Green Version]
  37. Dianati, K.; Schäfer, L.; Milner, J.; Gómez-Sanabria, A.; Gitau, H.; Hale, J.; Langmaack, H.; Kiesewetter, G.; Muindi, K.; Mberu, B.; et al. A system dynamics-based scenario analysis of residential solid waste management in Kisumu, Kenya. Sci. Total Environ. 2021, 777, 146200. [Google Scholar] [CrossRef]
  38. Liu, H.; Benoit, G.; Liu, T.; Liu, Y.; Guo, H. An integrated system dynamics model developed for managing lake water quality at the watershed scale. J. Environ. Manag. 2015, 155, 11–23. [Google Scholar] [CrossRef] [PubMed]
  39. Lemaire, G.G. Assessing the Spatio-Temporal Dynamics and Environmental Impacts in Peri-Urban Stream Systems. Ph.D. Thesis, Technical University of Denmark, Lyngby, Denmark, 2021. [Google Scholar]
  40. ISEE. Stella and Ithink v. 1.2.2. Technical Documentation. 2020. Available online: https://www.iseesystems.com/resources/help/v1-2/ (accessed on 2 August 2020).
  41. Khan, S.; Yufeng, L.; Ahmad, A. Analysing complex behaviour of hydrological systems through a system dynamics approach. Environ. Model. Softw. 2009, 24, 1363–1372. [Google Scholar] [CrossRef]
  42. Sehlke, G.; Jacobson, J. System dynamics modeling of transboundary systems: The river basin model. Ground Water 2005, 43, 722–730. [Google Scholar] [CrossRef]
  43. Ghashghaei, M.; Bagheri, A.; Morid, S. Rainfall-runoff Modeling in a Watershed Scale Using an Object Oriented Approach Based on the Concepts of System Dynamics. Water Resour. Manag. 2013, 27, 5119–5141. [Google Scholar] [CrossRef]
  44. Tian, Y.; Li, C.; Yi, Y.; Wang, X.; Shu, A. Dynamic model of a sustainable water resources utilization system with coupled water quality and quantity in Tianjin city. Sustainability 2020, 12, 4254. [Google Scholar] [CrossRef]
  45. Rivers, M.R.; Weaver, D.M.; Smettem, K.R.J.; Davies, P.M. Estimating farm to catchment nutrient fluxes using dynamic simulation modelling—Can agri-environmental BMPs really do the job? J. Environ. Manag. 2013, 130, 313–323. [Google Scholar] [CrossRef] [PubMed]
  46. Teegavarapu, R.S.V.; Tangirala, A.K.; Ormsbee, L. Modeling Water Quality Management Alternatives for a Nutrient Impaired Stream Using System Dynamics Simulation. J. Environ. Inform. 2005, 5, 73–81. [Google Scholar] [CrossRef]
  47. Brown Gaddis, E.J.; Vladich, H.; Voinov, A. Participatory modeling and the dilemma of diffuse nitrogen management in a residential watershed. Environ. Model. Softw. 2007, 22, 619–629. [Google Scholar] [CrossRef]
  48. Elshorbagy, A.; Ormsbee, L. Object-oriented modeling approach to surface water quality management. Environ. Model. Softw. 2006, 21, 689–698. [Google Scholar] [CrossRef]
  49. Venkatesan, A.K.; Ahmad, S.; Johnson, W.; Batista, J.R. Systems dynamic model to forecast salinity load to the Colorado River due to urbanization within the Las Vegas Valley. Sci. Total Environ. 2011, 409, 2616–2625. [Google Scholar] [CrossRef] [PubMed]
  50. Zhang, B.; Qin, Y.; Huang, M.; Sun, Q.; Li, S.; Wang, L.; Yu, C. SD-GIS-based temporal-spatial simulation of water quality in sudden water pollution accidents. Comput. Geosci. 2011, 37, 874–882. [Google Scholar] [CrossRef]
  51. McKnight, U.S.; Funder, S.G.; Rasmussen, J.J.; Finkel, M.; Binning, P.J.; Bjerg, P.L. An integrated model for assessing the risk of TCE groundwater contamination to human receptors and surface water ecosystems. Ecol. Eng. 2010, 36, 1126–1137. [Google Scholar] [CrossRef] [Green Version]
  52. Bergström, S. The HBV model—Its structure and applications. Swed. Meteorol. Hydrol. Inst. Norrköping 1992, 4, 1–33. [Google Scholar]
  53. de Doncker, L.; Troch, P.; Verhoeven, R.; Buis, K. Deriving the relationship among discharge, biomass and Manning’s coefficient through a calibration approach. Hydrol. Process. 2011, 25, 1979–1995. [Google Scholar] [CrossRef]
  54. Toffolon, M.; Piccolroaz, S. A hybrid model for river water temperature as a function of air temperature and discharge. Environ. Res. Lett. 2015, 10, 1–22. [Google Scholar] [CrossRef]
  55. Caissie, D.; Satish, M.G.; El-Jabi, N. Predicting river water temperatures using the equilibrium temperature concept with application on Miramichi River catchments (New Brunswick, Canada). Hydrol. Process. 2005, 19, 2137–2159. [Google Scholar] [CrossRef]
  56. Alnoee, A.B.; Levi, P.S.; Baattrup-Pedersen, A.; Riis, T. Macrophytes enhance reach-scale metabolism on a daily, seasonal and annual basis in agricultural lowland streams. Aquat. Sci. 2021, 83, 1–12. [Google Scholar] [CrossRef]
  57. Gotovtsev, A.V. Modification of the Streeter-Phelps system with the aim to account for the feedback between dissolved oxygen concentration and organic matter oxidation rate. Water Resour. 2010, 37, 245–251. [Google Scholar] [CrossRef]
  58. Simonsen, J.F.; Harremoës, P. Oxygen and pH fluctuations in rivers. Water Res. 1978, 12, 477–489. [Google Scholar] [CrossRef]
  59. APHA. Standard Methods for the examination of Water and Wastewaters, 18th ed.; APHA: Washington, DC, USA, 1992. [Google Scholar]
  60. Chapra, S.C. Surface Water-Quality Modeling; Waveland, P., Ed.; Waveland Press: Long Grove, IL, USA, 1997. [Google Scholar]
  61. Nazaroff, W.W.; Alvarez-Cohen, L. Environmental Engineering Science; Wiley: Hoboken, NJ, USA, 2001. [Google Scholar]
  62. Bowie, G.; Mills, W.; Porcella, D.; Campbell, C.; Pagenkopf, J.; Rupp, G.; Johnson, K.; Chan, P.; Gherini, S. Rates, Constants, and Kinetic Formulations in Surface Water Modeling; US Environmental Protection Agency: Athens, GA, USA, 1985.
  63. Krause-Jensen, D.; Sand-Jensen, K. Light attenuation and photosynthesis of aquatic plant communities. Limnol. Oceanogr. 1998, 43, 396–407. [Google Scholar] [CrossRef]
  64. Hall, R.O.; Beaulieu, J.J. Estimating autotrophic respiration in streams using daily metabolism data. Freshw. Sci. 2013, 32, 507–516. [Google Scholar] [CrossRef]
  65. Uusitalo, L.; Lehikoinen, A.; Helle, I.; Myrberg, K. An overview of methods to evaluate uncertainty of deterministic models in decision support. Environ. Model. Softw. 2015, 63, 24–31. [Google Scholar] [CrossRef] [Green Version]
  66. McKay, M.D.; Beckman, R.J.; Conover, W.J. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 2000, 42, 55–61. [Google Scholar] [CrossRef]
  67. Heebner, D.; Toran, L. Sensitivity Analysis of Three-Dimensional Steady-State and Transient Spray Irrigation Models. Ground Water 2000, 38, 20–28. [Google Scholar] [CrossRef]
  68. Danmark Miljøportal, Environmental Database. Available online: https://arealinformation.miljoeportal.dk/html5/index.html?viewer=distribution (accessed on 1 February 2020). (In Danish).
  69. Krüger. Measurements in Usserød Stream. 2010 (Målinger i Usserød Å—2010). 2011. (In Danish) [Google Scholar]
  70. Gørtz, P.; Schultz, J.R. Biological Conditions in Usserød Stream. (Den Biologiske Tilstand i Usserød Å). 2020. (In Danish) [Google Scholar]
  71. DSFE. Danish Ministry Data Supply and Efficiency. Map service (Kortforsyningen). 2020. Available online: https://kortforsyningen.dk/ (accessed on 10 February 2020).
  72. Rudersdal, Hørsholm, and Fredensborg. Hydroinform. 2020. Available online: http://hydroinform.dk/UsseroedIntern.html (accessed on 1 February 2020).
  73. DMI. SVK Bestilling—Rain Monitoring Stations. 2020. Available online: http://svk.dmi.dk/dmi/RainEvents/*.login. (accessed on 1 January 2018).
  74. DMI. DMI Open Data—Meteorological Observation. 2020. Available online: https://confluence.govcloud.dk/display/FDAPI/Meteorological+Observation (accessed on 15 October 2020).
  75. Kommune, H. Hørsholm Vand. Nedbør og Overløb Til Usserød å [In Danish]. Available online: http://hydroinform.dk/HorsholmVandAa.html (accessed on 1 January 2020).
  76. DMoriasi, N.; Arnold, J.G.; van Liew, M.W.; Bingner, R.L.; Harmel, R.D.; Veith, T.L. Model Evaluation Guidelines for Systematic Quantification of Accuracy in Watershed Simulations. Trans. ASABE 2007, 50, 885–900. [Google Scholar] [CrossRef]
  77. Baker, D.B.; Richards, R.P.; Loftus, T.T.; Kramer, J.W. A new flashiness index: Characteristics and applications to Midwestern rivers and streams. J. Am. Water Resour. Assoc. 2004, 40, 503–522. [Google Scholar] [CrossRef]
  78. Pérez-Sánchez, M.; Sánchez-Romero, F.J.; Ramos, H.M.; López-Jiménez, P.A. Calibrating a flow model in an irrigation network: Case study in Alicante, Spain. Spanish J. Agric. Res. 2017, 15, 1–3. [Google Scholar] [CrossRef]
  79. Hamid, A.; Bhat, S.U.; Jehangir, A. Local determinants influencing stream water quality. Appl. Water Sci. 2020, 10, 1–16. [Google Scholar] [CrossRef] [Green Version]
  80. Hutchins, M.G.; Williams, R.J.; Prudhomme, C.; Bowes, M.J.; Brown, H.E.; Waylett, A.J.; Loewenthal, M. Projections of future deterioration in UK river quality are hampered by climatic uncertainty under extreme conditions. Hydrol. Sci. J. 2016, 61, 2818–2833. [Google Scholar] [CrossRef] [Green Version]
  81. Bowes, M.J.; Loewenthal, M.; Read, D.S.; Hutchins, M.G.; Prudhomme, C.; Armstrong, L.K.; Harman, S.A.; Wickham, H.D.; Gozzard, E.; Carvalho, L. Identifying multiple stressor controls on phytoplankton dynamics in the River Thames (UK) using high-frequency water quality data. Sci. Total Environ. 2016, 569–570, 1489–1499. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  82. Hutchins, M.G.; Harding, G.; Jarvie, H.P.; Marsh, T.J.; Bowes, M.J.; Loewenthal, M. Intense summer floods may induce prolonged increases in benthic respiration rates of more than one year leading to low river dissolved oxygen. J. Hydrol. X 2020, 8, 100056. [Google Scholar] [CrossRef]
  83. Riis, T.; Tank, J.L.; Reisinger, A.J.; Aubenau, A.; Roche, K.R.; Levi, P.S.; Baattrup-Pedersen, A.; Alnoee, A.B.; Bolster, D. Riverine macrophytes control seasonal nutrient uptake via both physical and biological pathways. Freshw. Biol. 2020, 65, 178–192. [Google Scholar] [CrossRef]
  84. Cohen, R.A.; Fong, P. Nitrogen uptake and assimilation in Enteromorpha intestinalis (L.) Link (Chlorophyta): Using 15N to determine preference during simultaneous pulses of nitrate and ammonium. J. Exp. Mar. Bio. Ecol. 2004, 309, 67–77. [Google Scholar] [CrossRef]
  85. Jarvie, H.P.; Pallett, D.W.; Schäfer, S.M.; Macrae, M.L.; Bowes, M.J.; Farrand, P.; Warwick, A.C.; King, S.M.; Williams, R.J.; Armstrong, L.; et al. Biogeochemical and climate drivers of wetland phosphorus and nitrogen release: Implications for nutrient legacies and eutrophication risk. J. Environ. Qual. 2020, 49, 1703–1716. [Google Scholar] [CrossRef] [PubMed]
  86. Jaiswal, D.; Pandey, J. Human-driven changes in sediment-water interactions may increase the degradation of ecosystem functioning in the Ganga River. J. Hydrol. 2021, 598, 126261. [Google Scholar] [CrossRef]
  87. Blaszczak, J.R.; Delesantro, J.M.; Urban, D.L.; Doyle, M.W.; Bernhardt, E.S. Scoured or suffocated: Urban stream ecosystems oscillate between hydrologic and dissolved oxygen extremes. Limnol. Oceanogr. 2019, 64, 877–894. [Google Scholar] [CrossRef] [Green Version]
  88. Sand-Jensen, K.; Borg, D.; Jeppesen, E. Biomass and oxygen dynamics of the epiphyte community in a Danish lowland stream. Freshw. Biol. 1989, 22, 431–443. [Google Scholar] [CrossRef]
  89. Clark, M.P.; Bierkens, M.F.P.; Samaniego, L.; Woods, R.A.; Uijenhoet, R.; Bennet, K.E.; Pauwels, V.R.N.; Cai, X.; Wood, A.W.; Peters-Lidard, C.D. The evolution of process-based hydrologic models: Historical challenges and the collective quest for physical realism. Hydrol. Earth Syst. Sci. Discuss. 2017, 21, 3427–3440. [Google Scholar] [CrossRef] [Green Version]
  90. Booth, D.B.; Roy, A.H.; Smith, B.; Capps, K.A. Global perspectives on the urban stream syndrome. Freshw. Sci. 2016, 35, 412–420. [Google Scholar] [CrossRef] [Green Version]
  91. Fernando, A.M.E.; Súarez, Y.R. Resource use by omnivorous fish: Effects of biotic and abiotic factors on key ecological aspects of individuals. Ecol. Freshw. Fish 2021, 30, 222–233. [Google Scholar] [CrossRef]
  92. Wenger, S.J.; Roy, A.H.; Jackson, C.R.; Bernhardt, E.S.; Carter, T.L.; Filoso, S.; Gibson, C.A.; Hession, W.C.; Kaushal, S.S.; Martí, E.; et al. Twenty-six key research questions in urban stream ecology: An assessment of the state of the science. J. N. Am. Benthol. Soc. 2009, 28, 1080–1098. [Google Scholar] [CrossRef] [Green Version]
  93. Khamis, K.; Bradley, C.; Hannah, D.M. High frequency fluorescence monitoring reveals new insights into organic matter dynamics of an urban river, Birmingham, UK. Sci. Total Environ. 2020, 710, 135668. [Google Scholar] [CrossRef]
  94. Xenopoulos, M.A.; Barnes, R.T.; Boodoo, K.S.; Butman, D.; Catalán, N.; D’Amario, S.C.; Fasching, C.; Kothawala, D.N.; Pisani, O.; Solomon, C.T.; et al. How humans alter dissolved organic matter composition in freshwater: Relevance for the Earth’s biogeochemistry. Biogeochemistry 2021, 3, 1–26. [Google Scholar]
  95. Segatto, P.L.; Battin, T.J.; Bertuzzo, E. Modeling the coupled dynamics of stream metabolism and microbial biomass. Limnol. Oceanogr. 2020, 65, 1573–1593. [Google Scholar] [CrossRef]
  96. Forum, D.W. Unlocking the Potential of Wastewater; State of Green: Copenhagen, Denmark, 2016. [Google Scholar]
  97. Lemaire, G.G.; McKnight, U.S.; Schulz, H.; Roost, S.; Bjerg, P.L. Evidence of Spatio-Temporal Variations in Contaminants Discharging to a Peri-Urban Stream. Groundw. Monit. Remediat. 2020, 40, 40–51. [Google Scholar] [CrossRef]
  98. Sophocleous, M. Interactions between groundwater and surface water: The state of the science. Hydrogeol. J. 2002, 10, 52–67. [Google Scholar] [CrossRef]
  99. Brauer, C. Modelling Rainfall-Runoff Processes in Lowland Catchments. Ph.D. Thesis, Wageningen University, Wageningen, The Netherlands, 2014. [Google Scholar]
  100. Fones, G.R.; Bakir, A.; Gray, J.; Mattingley, L.; Measham, N.; Knight, P.; Bowes, M.J.; Greenwood, R.; Mills, G.A. Using high-frequency phosphorus monitoring for water quality management: A case study of the upper River Itchen, UK. Environ. Monit. Assess. 2020, 192, 3. [Google Scholar] [CrossRef] [Green Version]
  101. Bojanowski, J.S.; Brown, L.C. Assessing the Performance of Reaeration Prediction Equations. J. Environ. Eng. 2014, 140, 04013013. [Google Scholar]
  102. Bojanowski, J.S. Sirad: Functions for calculating daily solar radiation and evapotranspiration. R Package Vers. 2013, 140, 1–35. [Google Scholar]
  103. NOVAFOS. Survey of Public Perception and Activities in Usserød å [Technical Meeting and Data Transfer]; NOVAFOS: Hørsholm, Denmark, 2020. [Google Scholar]
  104. Bogan, T.; Mohseni, O.; Stefan, H.G. Stream temperature-equilibrium temperature relationship. Water Resour. Res. 2003, 39, 1–12. [Google Scholar] [CrossRef] [Green Version]
  105. Nielsen, S.L.; Sand-Jensen, K. Variation in growth rates of submerged rooted macrophytes. Aquat. Bot. 1991, 39, 109–120. [Google Scholar] [CrossRef]
  106. Julian, J.P.; Doyle, M.W.; Stanley, E.H. Empirical modeling of light availability in rivers. J. Geophys. Res. Biogeosci. 2008, 113. [Google Scholar] [CrossRef] [Green Version]
  107. Wild, M.; Hakuba, M.Z.; Folini, D.; Dörig-Ott, P.; Schär, C.; Kato, S.; Long, C.N. The cloud-free global energy balance and inferred cloud radiative effects: An assessment based on direct observations and climate models. Clim. Dyn. 2019, 52, 4787–4812. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  108. Carr, G.M.; Duthie, H.C.; Taylor, W.D. Models of aquatic plant productivity: A review of the factors that influence growth. Aquat. Bot. 1997, 59, 195–215. [Google Scholar] [CrossRef]
  109. Oudin, L.; Hervieu, F.; Michel, C.; Perrin, C.; Andréassian, V.; Anctil, F.; Loumagne, C. Which potential evapotranspiration input for a lumped rainfall-runoff model? Part 2—Towards a simple and efficient potential evapotranspiration model for rainfall-runoff modelling. J. Hydrol. 2005, 303, 290–306. [Google Scholar] [CrossRef]
  110. Caissie, D. The thermal regime of rivers: A review. Freshw. Biol. 2006, 51, 1389–1406. [Google Scholar] [CrossRef]
  111. Bott, T.L. Primary Productivity and Community Respiration, 2nd ed.; Elsevier Inc.: Amsterdam, The Netherlands, 2007; ISBN 9780123329080. [Google Scholar]
  112. Piedallu, C.; Gégout, J.-C. Multiscale computation of solar radiation for predictive vegetation modelling. Ann. For. Sci. 2007, 64, 219–228. [Google Scholar] [CrossRef] [Green Version]
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 Section 2.2.1 and Section 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 Section 2.2.1 and Section 2.2.2.
Water 13 03002 g001
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).
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).
Water 13 03002 g003
Figure 4. Selected results of flow and depth simulations for Usserød stream (a,b) and the Donse tributary (c) compared to available measured data. Other simulation results can be found in Figure S7 (Usserød flow), Figure S8 (Usserød water depth), Figure S9 (Donse tributary flow).
Figure 4. Selected results of flow and depth simulations for Usserød stream (a,b) and the Donse tributary (c) compared to available measured data. Other simulation results can be found in Figure S7 (Usserød flow), Figure S8 (Usserød water depth), Figure S9 (Donse tributary flow).
Water 13 03002 g004
Figure 5. Simulation results for the stream temperature compared to grab sample measurements at Parallelvej (red marked circle in map). The grey ribbon indicates the calibration period for the RR model, the rest of the time series being the verification.
Figure 5. Simulation results for the stream temperature compared to grab sample measurements at Parallelvej (red marked circle in map). The grey ribbon indicates the calibration period for the RR model, the rest of the time series being the verification.
Water 13 03002 g005
Figure 6. DO concentration simulation (blue line) compared to an online sensor measurement (red line) at Parallelvej, with (a) simulated and measured daily mean DO concentration for the calibration and verification period. The grey ribbon indicates the calibration period. The pale blue and red shadings correspond to the amplitude of simulated and measured daily extremes, respectively. (b) 1:1 scatter plot of simulated and measured daily average DO concentration (black cross-marker corresponds to the verification period, grey dot marker to the calibration period); and (c) simulated and measured hourly DO for a selected summer period showing the diurnal variation induced by plant photosynthesis/respiration.
Figure 6. DO concentration simulation (blue line) compared to an online sensor measurement (red line) at Parallelvej, with (a) simulated and measured daily mean DO concentration for the calibration and verification period. The grey ribbon indicates the calibration period. The pale blue and red shadings correspond to the amplitude of simulated and measured daily extremes, respectively. (b) 1:1 scatter plot of simulated and measured daily average DO concentration (black cross-marker corresponds to the verification period, grey dot marker to the calibration period); and (c) simulated and measured hourly DO for a selected summer period showing the diurnal variation induced by plant photosynthesis/respiration.
Water 13 03002 g006
Figure 7. MC simulation results at Parallelvej for the DO concentration (gray and black shading/lines) compared to the online sensor (red line; located under a bridge in the shade) and grab sample measurements (red circles, measured in the afternoon) j. The results are shown for the verification period only. The gray and black shadings show the 10–90 percentiles of estimated daily DO concentrations. The dashed grey lines show the 10–90 percentiles of maximum daily amplitude resulting from photosynthesis/respiration. Input parameter distributions can be found in Table S5.
Figure 7. MC simulation results at Parallelvej for the DO concentration (gray and black shading/lines) compared to the online sensor (red line; located under a bridge in the shade) and grab sample measurements (red circles, measured in the afternoon) j. The results are shown for the verification period only. The gray and black shadings show the 10–90 percentiles of estimated daily DO concentrations. The dashed grey lines show the 10–90 percentiles of maximum daily amplitude resulting from photosynthesis/respiration. Input parameter distributions can be found in Table S5.
Water 13 03002 g007
Figure 8. MC simulation results at Nivemølle. (a) temp; (b) DO; (c) NO3-N; (d) NH4-N; (e) PO4-P and (f) susp. chla compared to available grab sampled measurements (red circles). The dashed grey lines for the DO concentration (b) highlight the 10–90 percentile for the sub-daily DO variations. Results are displayed for the verification period only, with grab samples for nutrients and chla combined with simulated flow used as forcing input at the most upstream station Grønnegade and at Donse tributary (Figure S16). Parameter distributions for the MC simulation are given in Table S6.
Figure 8. MC simulation results at Nivemølle. (a) temp; (b) DO; (c) NO3-N; (d) NH4-N; (e) PO4-P and (f) susp. chla compared to available grab sampled measurements (red circles). The dashed grey lines for the DO concentration (b) highlight the 10–90 percentile for the sub-daily DO variations. Results are displayed for the verification period only, with grab samples for nutrients and chla combined with simulated flow used as forcing input at the most upstream station Grønnegade and at Donse tributary (Figure S16). Parameter distributions for the MC simulation are given in Table S6.
Water 13 03002 g008
Figure 9. Simulations of the contribution of the different flow components at Nivemølle. (a) Measured precipitation in the model area; (b) Simulated flow in Usserød Stream, (c) Simulated flow contribution for the key features captured by the model; and (d) zoom to highlight a CSO event. The contribution is based on the sum of all inflows (lake (blue); flow from pervious areas (dark green); tributaries (light green); various urban features (WWTP: dark grey; SS: light gray)) and does not account for possible loss along the reaches. (d) Selected time period showing the contribution of a CSO event (purple shading).
Figure 9. Simulations of the contribution of the different flow components at Nivemølle. (a) Measured precipitation in the model area; (b) Simulated flow in Usserød Stream, (c) Simulated flow contribution for the key features captured by the model; and (d) zoom to highlight a CSO event. The contribution is based on the sum of all inflows (lake (blue); flow from pervious areas (dark green); tributaries (light green); various urban features (WWTP: dark grey; SS: light gray)) and does not account for possible loss along the reaches. (d) Selected time period showing the contribution of a CSO event (purple shading).
Water 13 03002 g009
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.
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.
Output
Indicator
CalibrationVerification
[2017–31 October 2018][1 November 2018–10 October 2019 */31 December 2019]
GrønnegadeÅdalsvejNivemølleFredtoftenBrønshomsdalsvejGrønnegade *Ådalsvej *Nivemølle *FredtoftenBrønshomsdalsvej
Usserød StreamDonse TribUsserød StreamDonse Trib
Flow
NSE0.760.730.760.800.570.900.760.870.660.87
RSR0.480.510.490.450.650.320.480.350.580.35
PBIAS3.3−9.30.013.330 x10.916.48.3−42.7 x−1.3
RB0.190.170.200.140.150.180.180.200.130.15
Water depth
NSE0.45 x0.590.67 0.610.700.60
RSR0.74 x0.640.57 0.620.550.63
PBIAS7.5−6.1−7.8 4.4−4.8−8.3
RMSE0.070.050.07 0.080.050.08
*: indicator given for the period [1 November 2018–10 October 2019, see Section 5].
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lemaire, G.G.; Carnohan, S.A.; Grand, S.; Mazel, V.; Bjerg, P.L.; McKnight, U.S. Data-Driven System Dynamics Model for Simulating Water Quantity and Quality in Peri-Urban Streams. Water 2021, 13, 3002. https://doi.org/10.3390/w13213002

AMA Style

Lemaire GG, Carnohan SA, Grand S, Mazel V, Bjerg PL, McKnight US. Data-Driven System Dynamics Model for Simulating Water Quantity and Quality in Peri-Urban Streams. Water. 2021; 13(21):3002. https://doi.org/10.3390/w13213002

Chicago/Turabian Style

Lemaire, Gregory G., Shane A. Carnohan, Stanislav Grand, Victor Mazel, Poul L. Bjerg, and Ursula S. McKnight. 2021. "Data-Driven System Dynamics Model for Simulating Water Quantity and Quality in Peri-Urban Streams" Water 13, no. 21: 3002. https://doi.org/10.3390/w13213002

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop