Evidence-based identification of integrated water quality systems

Identification of integrated models is still hindered by submodels’ uncertainty propagation. In this article, a novel identifiability and identification framework is applied to screen and establish reasonable hypotheses of an integrated instream (WASP) and catchment water quality (VENSIM) model. Using the framework, the models were linked, and critical parameters and processes identified. First, an ensemble of catchment nutrient loads was simulated with randomized parameter settings of the catchment processes (e.g. nutrient decay rates). A second Monte Carlo analysis was then staged with randomized loadings and parameter values mimicking insteam processes (e.g. algae growth). The most significant parameters and their processes were identified. This coupling of models for a two-step global sensitivity analysis is a novel approach to integrated catchment-scale water quality model identification. Catchment processes were, overall, more significant to the river’s water quality than the instream processes of this Prairie river system investigated (Qu’Appelle River).


Introduction
The requirements of recent policies for sustainable management of freshwater aquatic ecosystems in Canada and globally demand that large-scale models are used to provide answers for the emerging issues of aquatic ecosystem sustainability (EC 2000;Government of Canada 2007). Use-inspired modeling will, thus, effectively support informed decision-making by meeting modeling objectives.
The lack of evidence in the optimum performance of isolated managed systems, limited field data, the need to address connection issues between and within systems, and to meet the objectives of current policies such as sustainable resource allocation, have prompted integrated modeling solutions (Rauch et al. 2002;Willems and Berlamont 2002;Benedetti et al. 2013;Keupers and Willems 2017;Tranmer et al. 2018). Integrated modeling strategy has been applied to optimize field monitoring campaigns and to identify reasonable hypotheses of functional relationships governing the behavior of an existing system (Reichert 1994;Freni and Mannina 2012).
Mechanistic models with flexible structures, together with laboratory and field experiments, are typically used in establishing internal processes controlling the behavior of the environment system due to their complex and evolving nature (Reichert 1994). Two of the major issues currently confronting the advancement of large scale integrated water quality modeling are the use of very simple structures to represent complex catchment nutrient export and identification of practical parameters and structures of model constructs (Matott, Babendreier, and Purucker 2009;Frassl et al. 2019;Fu et al. 2019). Identification is used here as in the sense of the establishment of realistic hypotheses for the integrated system under consideration. The process is typically preceded by an assessment of whether system parameters can be uniquely determined from available datasets of the system, known as identifiability (Reichert 1994). Sparse monitored data at the catchment scale and the desire to extrapolate laboratory/field scale studies to catchment scale usually results in the use of either simple statistical relationships to represent complex catchment processes or the increment of model complexities that cannot be justified by available data. Either of these approaches, however, ends up with models with deep uncertainties regarding the system's structure and parameters. A useful integrated model will, thus, provide a cogent balance between a system's process representation and identifiability.
The issues of large-scale model structure and parameter identification are magnified and lagging in the cold regions. Model development in these regions is impeded by limited large-scale studies, over-parameterized models, and the ontology of a system's biogeochemistry (Costa et al. 2020). Land to water delivery of nutrients strongly influences snow mechanisms impacting overland flow, erosion, and diffuse loading in these regions (Eimers et al. 2009;Casson, Eimers, and Watmough 2012).
Identification of river or lake water quality models has received much attention in the literature. However, for integrated catchment-surface water quality models used for long-term management, less has been reported on parametric and process identification, and the effect of forcings on submodels (Yuan et al. 2015) at temporal, and spatial scales during model developments. Less progress has been reported in the development of such models during the last decade (Fu et al. 2019). Reliable predictions of large-scale water quality models require, among other things, the rigorous establishment of physical, biological, geological, and chemical bases of the system. This is especially the case for northern catchments with the unique characteristics of unexpected spring and winter blooms, and ice-covered river-lake systems. Iterative mathematical modeling and hypotheses refinement through field and laboratory investigation can assist in the successful development of such integrated water quality models. Careful parameterization can provide the means to fully understand the feedback of fundamental processes of a system at different scales (National Research Council 2005). Therefore, to advance ICWM, more attention is required on the identification of model constructs to justify the need for structural improvement or more data collection (Matott, Babendreier, and Purucker 2009;Marsili-Libelli et al. 2014;Guillaume et al. 2019). Due to the co-evolution of complex systems, continuously refining model constructs with a new system's trajectory will allow for precise predictions. Furthermore, delineation of the interaction and cross correlations between landscape characteristics and instream water quality has been limited (Lintern et al. 2018). There is currently an incomplete understanding of how the relationships between landscape characteristics and water quality responses can shift based on the other characteristics of the catchment especially in cold-climates.
System identification is iterative, and traditionally involves model parameter identifiability, initial state and parameter estimation, and model validation. During identifiability analysis, measured data are compared with a model to assess whether parameters can be uniquely and accurately determined (Reichert 1994).
The numerical identifiability criteria technique is proving to be robust in the practical identification of integrated urban water quality models due to its strong theoretical base of yielding naturally, the influence of parameters on model output (Weijers and Vanrolleghem 1997;Checchi and Marsili-Libelli 2005;De Pauw 2005;Freni, Mannina, and Viviani 2008;Marsili-Libelli and Giusti 2008). Conversely, the approach has not been prominent in non-urban (like the study site) integrated catchment-scale water quality modeling.
Consequently, a new stepwise practical/posteriori identifiability and identification framework (Lindenschmidt et al. 2020) is used to highlight critical parameters and processes underlying the integrated water quality model. Unlike the orthodox system identification, the new framework employs a multistage global sensitivity analysis in the identification process to limit error propagation between the coupled submodels and allows for the incorporation of new information about the system being studied.

Study area
The Qu'Appelle River (QR) drains a 52,000 km 2 catchment with its headwaters below the Qu'Appelle Dam and finally joins the Assiniboine River in Manitoba ( Figure 1). The main land-use in the valley is agriculture, with cattle, swine, canola, barley, and wheat farming covering over 95% of the catchment area. Productivity in the reaches and lakes of the Qu'Appelle increases downstream from Buffalo Pound Lake in the upper reaches (eutrophic) to Katepwa Lake (middle) and Crooked Lake in the lower Qu'Appelle (hyper-eutrophic). In addition to diffuse loadings from the catchment, Pasqua Lake, the first of four lakes in series on the middle Qu'Appelle, receives point source loads from Regina and Moose Jaw. These lakes are polymictic, except Katepwa and subsaline (Hall et al. 1999). Cyanobacteria and cyclopoid copepods dominate the plankton community in the summer, while bigmouth buffalo, cisco, yellow perch, northern pike, and walleye represent fish communities in the region (Patoine, Graham, and Leavitt 2006;Vogt et al. 2011). Some of the major tributaries contributing to the middle reaches include Wascana Creek, Last Mountain River, Jumping Deer Creek, and Loon Creek.
The river and lakes are under ice for almost 6 months with an average temperature of À20.1 C, in January and 25.8 C in July (Environment Canada 2017).

The integrated catchment-instream model
Two submodels were linked to hypothesize the transport, transformation, and delivery of diffuse loadings from contributing catchments and coupled cycling of dissolved oxygen balance, nitrogen, phosphorous, and phytoplankton in the main stem of the QRB (Figure 2). The first submodel (Figure 2), the Qu'Appelle nutrient loading (QR_NL) model, is an extension of the SD-Qu'Appelle (SQ), initially developed to estimate mean annual loadings from the QRB (see Hassanzadeh et al. 2019). In the SQ model, mobilization of nutrients, total nitrogen and total phosphorus (TN and TP) from sources is characterized by their export coefficients. Delivery and attenuating factors are used to predict the integrated effect of landscape features on the fate of nutrients per unit length of streams, wetland, and / or reservoirs. Geomorphologic features such as soil properties, average overland flow characteristics, wetlands, and management practices influence the transport and transformation of catchment nutrients. These landscape characteristics represent the net effect of prevailing processes such as nutrient uptake, storage, mineralization, settling, and denitrification. Estimated fluxes at the outlet of each subcatchment represent all attenuated nutrients from sources upstream of the outlet location (Alexander et al. 2002;Kim et al. 2014).
The steady-state assumption inherent in the SQ model limits its ability to accurately define seasonal variation in diffuse loading. However, seasonality has been found to be a strong predictor of nutrient export in the Prairies (Gomi, Sidle, and Richardson 2002). In particular, meltwater plays a key role in the land to water delivery of nutrients, contributing up to 80% of their export (Corriveau et al. 2011;Uusi-K€ ampp€ a et al. 2012;Cade-Menun et al. 2013). Compared to summer rainfall-runoff, snowmelt runoff is marked by frozen ground, unsaturated flow, and relatively less nutrient uptake in the Canadian Prairies. In effect, spring runoff is relatively faster with accentuated loading due to restricted infiltration and lower temperatures (Gray, Landine, and Granger 1985;Cade-Menun et al. 2013). Diffuse loadings in regions with a continental climate exhibit peak rates in spring and minimum values in the winter.
In this study, the base SQ model was extended to include seasonality in the estimated long-term nutrient loads. Equation (1) represents the modified model structure.
Y ¼ Estimated daily load at subcatchment i (kg/day) a j ,b j ,c j ,and e j ¼ Non-point source N export coefficient, land to water delivery coefficient, point source N export coefficient, wetland/lake/reservoir decay coefficient, respectively, The second submodel which represents eutrophication in the main stem of QR was modeled using the advanced EUTRO module of Water quality Analysis Simulation Program (WASP) (Akomeah, Lindenschmidt, and Chapra 2019). The transport, loading, and transformation cycles of nitrogen, phosphorus, phytoplankton, and dissolved oxygen balance were integrated in WASP. A one-dimensional mass balance equation was solved for the conversion of each of the eight state variables according to Equation (2).
Direct and diffuse loading rate, g/m 3 -day S B ¼ Boundary loading rate (including upstream, downstream, benthic, and atmospheric), g/m 3 -day S K ¼ Total kinetic transformation rate; positive is source, negative is sink, g/m 3 -day The cycling of phytoplankton in the WASP submodel is linked to nutrient cycles, and DO, and has a critical effect on pH. Forms of nitrogen simulated include organic N (DON, phytoplankton N), and inorganic N (NH 4 þ -N, NO 3 --N). NH 4 þ -N and NO 3 --N are taken up and incorporated into the biomass of phytoplankton during photosynthesis. Nutrient uptake is based on Monod growth kinetics, with growth following saturation rate given temperature, light, and sufficient substrate concentration. As the concentration of ammonium accentuates, however, uptake is solely based on ammonium. Respiration and death of phytoplankton recycle and transform cellular N to NH 4 þ -N (1-fon) and detrital N (fon). The model formulation also includes the N 2 fixation and sediment diagenesis describing processes such as sedimentation, mineralization, and surface exchanges at the sediment-water interface (Wool et al. 2006). Catchment-scale nutrient export and instream nutrient, phytoplankton, and DO transformation processes were loosely coupled using a python wrapper. The wrapper calls the VENSIM model executable and runs it. The results of each run are then decomposed into daily loads and partitioned to the various species of TN (ON, NH 4 þ -N, and NO 3 À -N) and TP (OP and PO 4 3À -P). Partitioning ratios were determined from measured data. The resulting VENSIM output file is further formatted according to WASP's non-point source (NPS) file format and passed to WASP. WASP then uses the NPS as an input together with other forcings to simulate eutrophication in the QRS.

Evidence-based identifiability and identification modeling framework
A newly stepwise multivariate identification framework was used to highlight critical parameters and processes that regulate eutrophication in the integrated model ( Figure  3). The framework is based on a Bayesian scheme using Monte Carlo filtering and model update using new measured data. Details of the framework can be found in (Lindenschmidt et al. 2020).
To implement the framework, first, uncertainties associated with model parameters, structures, and forcings are constrained within practical space in the catchment nutrient model (Osidele and Beck 2001). At this stage, site measured parameter ranges are preferred as much as possible. For this study, to adequately cover parameter space and minimize model computation time, 10,000 parameter sets were generated for the catchment model using uniform random sampling (Sarrazin, Pianosi, and Wagener 2016).
See Table 1 for the integrated model parameters and ranges. Next, the randomized parameterization for the catchment nutrient model is then staged in a Monte Carlo analysis to yield an ensemble of feasible nutrient loadings. This step within the framework is designed to limit error propagation that ensues from using only a calibrated/fixed parameterized catchment model. Third, the ensemble of nutrient loadings is then passed to the instream eutrophication model. In this study, the ensembles were passed to the non-point source (NPS) file of the instream WASP model to participate in a second global sensitivity analysis using randomized instream parameters, constrained within practical limits (Osidele and Beck 2001). The posteriori distribution of this coupled Monte Carlo simulation was then compared with instream observations. Estimated output performance metrics and input factors of the integrated model were then ranked and partitioned into equal sample sizes (Wagener et al. 2001). A cumulative Distribution Function (CDF) was then generated for each parameter. A sensitivity index, which measures the sensitivity of the model performance metric to the integrated input factor was then estimated by taking the absolute vertical difference between the CDFs for each group (Noacco et al. 2019).
In this study, integrated parameter population, posteriori performance metric, and the relative root mean square error (RRE) (Akomeah, Lindenschmidt, and Chapra 2019), were ranked from best to worst and grouped into ten of equal size (that is 1000 parameters in each group). A cumulative distribution (CDF) of the likelihood of RRE associated with each of the parameter sets versus the parameters was then plotted (Wagener and Kollat 2007). A sensitivity index for each parameter was computed as  the maximum vertical distance between two CDFs using two-sample Kolmogorov-Smirnov (KS) statistic, Equation (3)  . The higher the index, the higher the sensitivity.
where F i () is the CDF of x in the ith group and F j () is the CDF in the jth group. Estimated maximum vertical distance between CDFs was then compared with critical difference values of the KS test. Because of noise, error and the filtering role played by RRE (Guillaume et al. 2019), critical parameters and processes were identified at a significance level of a ¼ 0.001. This level was selected to discard the null hypothesis that two CDFs are from similar distributions. Values above this critical value were therefore considered very significant, and hence, identified as parameters and processes controlling eutrophication in the studied system. Moreover, moderate and less significant parameters were identified at significance levels 0.01< a < 0.001 and a < 0.01, respectively.
After parameters and processes are identified, there is a need to field verify them. Based on the field studies, the coupled systems can then be used for forecasting purposes (probabilistic forecasting) or the model updated with new found information.

Results and discussion
3.1. Global sensitivity analysis (GSA) and identification of eutrophication controls for an annual average condition From the CDF plots, sensitivity indices, and identified parameters, net DO concentration in the QR is sensitive to critical parameters: N fertilizer export coefficient(CNfert), P fertilizer export coefficient ( Meanwhile CNfert, K12C, K71C, Organic Phosphorus Mineralization Rate (K83), CCHL, K1C, and K1T exert influence on DO in Pasqua Lake. These identified parameters were estimated by combining the sensitivity indices with an error margin of 0.1% (that is, a significance level a ¼ 0.001; Figure A.1 (bottom) in Appendix (online supplementary data), and Figure 6).
The parametric sensitivities in the model are associated with landscape-receiving water processes that induce variation in DO concentration in the QR. The underscored DO parameter subset posits that DO flux in the river results from photosynthetic DO production by phytoplankton which is enhanced by uptake of nutrients (NH 4 þ -N, NO 3 À -N, and PO 4 3À -P) in the presence of adequate light and temperature, and to a lesser extent, by reaeration. Nutrient availability for primary productivity in the river appears to be linked largely to mineralization and nitrification of available organic matter (OM) from the contributing catchment, and then, substrate. N, P, dead leaves, insects, and plants exported from agricultural land, corrals (upstream of Pasqua), and forest catchments to the QR supply inorganic nutrients (NO 3 À -N, NH 4 þ -N, and PO 4 3À -P), and ON and OP which are transformed to NH 4 þ -N and NO 3 À -N for plant growth in the water column. In addition to utilizing oxygen to decompose ON and oxidize NH 4 þ -N in the water column, benthic metabolic activities on settled organic matter also draw on epilimnetic DO, with the river metabolism weighing heavier on DO than catchment processes. From Appendix Figure A.1 and Figure 6, it can be deduced that the degree of influence of instream processes on DO availability does exceed catchment processes by an average of 42%. In Pasqua Lake, DO availability is further  mediated by mineralization of organic phosphorus, in lake temperature, and the lake's polymixis ( Figure 6).
Critical parameters regulating net NH 4 þ -N concentration in the river at HWY6 and upsofPasqua stations include CNfert, Coefficient of Permeability (CPER), CPFert, Temperature Coefficient (CT), Coefficient of Decay for wetlands (CW), and K71C. Furthermore, while Precipitation Coefficient (CP) is a critical parameter at HWY 6, it has a moderate influence on NH 4 þ -N at upsofPasqua (Figures 6 and 7 and Appendix Figure A.2).
The nutrient pathway delineated by these identified parameters includes nutrient mobilization from areal extents in the subcatchments with excess N and P fertilizers applied in the growing season and intervened by heavy summer rainfall and spring snowmelt. Nutrient flux delivered to the QR is impacted by landscape characteristics including soil type, temperature, streams, and wetlands. The mass of nutrients being transported in the shallow streams of the QRB is significantly influenced by their contact with the stream bed through retention of particulate portions (PO 4 3À -P and NH 4 þ -N) and denitrification. The morphometry of reservoirs, lakes, and wetlands in the surrounding catchments also aid in sinking portions of nutrient being transported. Particulate portions of the nutrient tend to settle and retain in these surface water bodies. Other mass balance modeling studies have shown how these landscape features retain nutrients (Chapra 1975;Kirchner and Dillon 1975;Molot and Dillon 1993;Westphal et al. 2019). The loss rate in the water column of these systems is a net of average subcatchment loadings, N fixation and benthic processes including Figure 6. Identified parameters for each variable in the river and lake. Critical parameters are highlighted red, moderately sensitive parameters yellow and less sensitive parameters green. Blue bar represents sensitivity index. Colour online. sedimentation and denitrification. Once delivered to the QR, inorganic nutrients become readily available for uptake, while the organic portion of the exported OM undergoes mineralization and sedimentation. To a lesser extent, decomposition of substrate also contributes to the nutrient pool in the pelagic zone with increasing influence from upstream to downstream (Figure 6).
Similarly, net NO 3 À -N concentration at these locations on the river traverses the identified path with an additional process, nitrification following mineralization of allochthonous and autochthonous organic matter, especially in the lake. From Figure 6, it is apparent that landscape processes exert dominant control on nutrient enrichment in the river than instream processes. In addition to the nutrient mechanisms in the riverine stretch, mobilization of nutrients from corrals also moderately governs nutrient availability in the lake. Inductions from this study's results are consistent with studies by Rattan et al. (2017), who also found a strong correlation between landscape activities and stream nutrients. Higher nutrient concentrations were measured during the snowmelt period.
Proliferation of the phytoplankton community in the lake was mediated by mobilization of nutrients from agricultural lands and wetlands in the catchment, phytoplankton stoichiometric make-up, irradiance, and seasonality. Variation in C:Chla ratio in response to nutrient and light availability drives the primary productivity and cycling of nutrients in Pasqua Lake (Appendix Figure A.4, and Figure 6).

Parameter and process modulation by seasonality
3.2.1. Spring-summer eutrophication mechanism Terrestrial N mobilization and delivery processes appear to exert greater control on net ammonium and nitrate concentrations in the river and lake. Prevalence of these mechanisms is during the freshet and rainfall in the early part of the growing season. The degree of instream microscopic process influences on nutrients, however, accentuates þ -N at HWY 6. The first twenty-two subplots counting from the left (CcorrN to RsCP) (see subplots enclosed in the green box) represent catchment submodel parameters while the next eighteen subplots (K12C to KMNG1) represent instream model parameters.
downstream. Toward late fall, when soil temperature is about 10 C, applied N fertilizer (e.g. granular urea, anhydrous ammonia, and urea ammonium nitrate) tend to accumulate (Davies et al. 1987;Pomeroy et al. 2006) due to reduced/no rainfall events (less dilution), decreasing temperature, and terrestrial microbial activities including denitrification. In the cold winter season, the frozen ground prevents leaching and retards biogeochemical reaction rates. This antecedent nutrient condition together with winter cattle feeding sets the stage for higher nutrient delivery in the spring and earlier summer. Studies by Corriveau, Chambers, and Culp (2013) show how elevated levels of TN were recorded during the winter in a Prairie catchment. The rising temperature in the spring causes previously frozen ground to thaw and snow to melt. With the large volume of snowpack accumulated over the winter period, the melt results in an expanse of runoff over less restricted terrain. The nutrients are transported to the river through preferential flow paths with load magnitude as a function of flow rate and nutrient concentration (Costa and Pomeroy 2019). The nutrients retained in the snowpack are, thus, released and delivered to the river. Spring nutrient export accounted for more than 80% of annual N export (Figures 8-10) (Corriveau, Chambers, and Culp 2013).
During summer rainfall, the erosive potential of the runoff increases with a resulting mobilization of dissolved and particulate nutrient to the river. Instream nutrient enrichment depends on the constituent delivery ratio and connectivity (Lintern et al. 2018). These influence how much nutrient is ultimately delivered to the river. Processes such as sorption of ammonium onto stream sediments, sedimentation of particulate organic N in wetlands, denitrification, and uptake by plants reduce nutrient connectivity. NH 4 þ -N and NO 3 À -N received from the catchment during spring snowmelt and summer rainfall render nutrient readily available to phytoplankton in the pelagic zone of the river compared to the longer instream N recycling through mineralization of detrital and substrate. Meanwhile, organic N and NH 4 þ -N sorbed onto particles that is delivered to the river undergo mineralization, nitrification, and Figure 8. Spring-summer cumulative distribution functions (CDF) of the coupled model parameters for the variable NH 4 þ -N at HWY 6. The first twenty-two subplots counting from the left (CcorrN to RsCP) (see subplots enclosed in the green box) represent catchment submodel parameters while the next eighteen subplots (K12C to KMNG1) represent instream model parameters. Figure 9. Spring-summer cumulative distribution functions (CDF) of the coupled model parameters for the variable NO 3 À -N at HWY 6. The first twenty-two subplots counting from the left (CcorrN to RsCP) (see subplots enclosed in the green box) represent catchment submodel parameters while the next eighteen subplots (K12C to KMNG1) represent instream model parameters. Figure 10. Sensitivity of NH 4 þ -N, NO 3 À -N, DO, and Chl a to variation in parameters for the seasons: starting from the leftspring-summer, fall, and winter (December). then uptake. These internal processes accelerate phytoplankton growth rate ( Figure 10).
Nutrient uptake is also dependent on the phytoplankton carbon to chlorophyll ratio. The instream processes are profound at the more sluggish part of the river, just upstream of the lake where residence time for organic N is longer (Figure 10).
Availability of DO in the river during the spring/summer is dependent primarily on instream primary productivity, nitrification, N export, and mineralization processes. The rate of atmospheric DO dissolution declines with rising temperature, especially in the summer months. Photosynthetic oxygen production in the river is made possible by the rising temperature during spring and summer with available irradiance and nutrients. The available oxygen raises the redox potential of the river, allowing for the oxidation of organic N and NH 4 þ -N ( Figure 10). In addition to landscape nutrient export mechanisms, mineralization, nitrification, and uptake that interact to govern N cycling in the river, N fertilization in the lake is further influenced by water temperature, oxidizing/reducing conditions, sediment oxygen demand, and phytoplankton stoichiometry. With the influx of ON and nutrients from the catchment and overturn, the rising of spring and summer temperatures spurs microbial activities such as mineralization, nitrification, and uptake in the water column and sediments.
Consequently, growth and oxygen production are speeded up as a result of the influx of NH 4 þ -N and NO 3 À -N from the catchment and sediments with the typical consequential spring maximum phytoplankton biomass (dominated by diatoms). In the sediments, settled organic matter is mineralized and at the sediment-water interface, oxygen is drawn to aid nitrification of NH 4 þ -N into NO 3 À -N which diffuses into the water column during circulation. During anoxic conditions in the sediment, NH 4 þ -N is released into the water column ( Figure 10).
From the foregoing, terrain processes are a more dominant influence in nutrient enrichment of the river and lake during spring/summer than instream nutrient cycling.

Fall eutrophication mechanism
Instream processes including mineralization, nitrification, phytoplankton growth rate, water temperature, and reaeration exert a strong influence on ammonium and nitrate availability in the river during fall. Increased solubility, longer residence time due to reduced flow, available organic matter accumulated over the growing season cause mineralization and subsequent NH 4 þ -N oxidation to maintain N transformation in the river. Reaeration rate exceeds photosynthetic oxygen production rate due to the declining temperature. Draw on DO for water column and sediment mineralization and nitrification is immense with limited photosynthesis due to the falling temperature, limited diffuse loading. Phytoplankton growth rate slows down compared to summer rates ( Figure 10).

Winter eutrophication mechanism
In the winter (March), snowmelt N mobilization and delivery processes control nutrient availability (NH 4 þ -N and NO 3 À -N) in the river. In December, however, at HWY 6, organic N mineralization, nutrient uptake, and temperature regulate NH 4 þ -N ( Figures  10 and 11).
Although winter sampling frequency was low, above zero temperature conditions accounted for the landscape N export (see Appendix Figure A.5 and Figure 10). When temperatures are above zero in the winter, meltwater carries and delivers concentrated levels of nutrients accumulated from livestock farming areas close to the river and fertilizer applications in the fall. Positive temperatures during the winter create a grade in the snow crystals which results in snow fractionation and densification. These transformations lead to ion separation from the crystals and eventual elution. Rivers receive episodic loadings from such landscape mechanisms (Appendix Figure A.5 and Figure  10) (Pomeroy et al. 2006;Lilbaek and Pomeroy 2008;Bormann et al. 2013).
A warm winter correspond to NO 3 À -N export in the winter (Casson, Eimers, and Watmough 2012). Elsewhere, external N flux was attributed to the riverine accumulation of NO 3 À -N in the winter (Knowles and Lean 1987;Pettersson et al. 2003;George et al. 2010). In another study, NO 3 À -N export in the winter contributed up to 40% of annual exports (Eimers, Buttle, and Watmough 2007). This was mainly due to events such as rain-on-snow and flow over frozen ground. In the winter, snow cover and frozen ground limit NO 3 À -N contact with the soil and retard denitrification rates. Microbial activities and temperature in the sediment play prominent roles in DO availability in winter (see also Akomeah and Lindenschmidt 2017;Hosseini et al. 2018;Terry et al. 2018). External mechanisms such as reaeration become limited due to the presence of ice cover (Figure 10).

Conclusions
This study has demonstrated the utility of applying a novel stepwise practical identification of an integrated catchment-river-lake water quality model in the light of parsimony. Two models, VENSIM and WASP were coupled using a python wrapper to hypothesize catchment nutrient export and instream eutrophication. A stepped Figure 11. Winter (December) Cumulative Distribution Functions (CDF) of the coupled model parameters for the variable NH 4 þ -N at HWY 6. The first twenty-two subplots counting from the left (CcorrN to RsCP) (see subplots enclosed in the green box) represent catchment submodel parameters while the next eighteen subplots (K12C to KMNG1) represent instream model parameters.
posteriori identification analysis was used to glean key parameters and processes modulating river eutrophication: catchment versus instream processes.
Findings from the study show the impact of diffuse loading on incremental instream eutrophication. Although the river is eutrophic, diffuse loading from the agricultural lands (crop and livestock farming) is exacerbating nutrient enrichment in the river. This is especially the case in the spring/summer months. Surprisingly in this study, nutrient export in the winter accounted for elevated instream NO 3 À -N concentration. Winter warming due to climate change will, thus, worsen the trophic state of the system. In the fall and subzero winter conditions, instream N cycling processes sustain the eutrophication.
Our results also show that there is a threshold of landscape processes until they start to dominate instream nutrient enrichment (Akomeah, Lindenschmidt, and Chapra 2019). In the Qu'Appelle river basin, this threshold is above 20% of instream nutrient levels. The stepwise practical identifiability analysis is, thus, another tool to further the identification of integrated catchment scale water quality models, especially in largescale data poor catchments. To capture the evolution of the components of the study site, identified parameters and processes would have to be validated in the field before forecasting the future behavior of the river. The developed framework will be key in the sustainable management of ecological systems.