Quantifying Spatial Changes in the Structure of Water Quality Constituents in a Large Prairie River within Two Frameworks of a Water Quality Model

Abstract: A global sensitivity analysis was carried out on a water quality model to quantify the spatial changes in parameter sensitivity of a model of a large prairie river, the South Saskatchewan River (SSR). The method is used to assess the relative impacts of major nutrient loading sources and a reservoir on the river’s water quality. The river completely freezes over during winter; hence, the sensitivity analysis was carried out seasonally, for winter and summer, to account for the influence of ice-covered conditions on nutrient transformations. Furthermore, the integrity of the river’s aquatic ecosystem was examined through the inter-relationship between variables and comparing hierarchy index values and water quality indices at four locations along the river. Sensitivities of model parameters varied slightly at different locations along the river, with the phytoplankton growth rate being the most influential parameter. Nitrogen and phosphorus transformation processes were more sensitive in winter, while chlorophyll-a and dissolved oxygen parameters showed higher sensitivity in summer. A more complicated correlation between variables was observed downstream of the junction of the Red Deer River. Our results reveal that the lower correlation between variables may suggest a more balanced and healthier system, although further analysis is needed to support this statement.


Introduction
Increased nutrient loads in rivers and reservoirs lead to eutrophication and, in turn, degrade the surface waters' ecosystem health [1,2].High nutrient enrichment can lead to major health issues in aquatic ecosystems, with the most common concern being excessive algal blooms [2].Monitoring nutrients and quantifying water quality conditions of surface waters are both vital when attempting to manage eutrophication.Utilizing measured data alone is often not informative enough for making water management decisions, due to the resolution coarseness of spatial and temporal data [3].Thus, water quality modeling and theoretical analysis are essential for water quality management [2,3].A wide variety of surface water quality modeling studies have been undertaken to estimate nutrient concentrations in waterbodies.However, only a few modeling studies have been carried out on surface waters under ice conditions.Ice covers affect the hydrology and the biological and chemical conditions of rivers, particularly evident in changes in dissolved oxygen concentrations [4,5].
The water quality analysis was augmented with sensitivity and uncertainty analyses to evaluate the model's behavior and identify dependencies within the model structure [6][7][8][9].A variety of techniques exist for sensitivity analysis, including local sensitivity andglobal sensitivity [6][7][8][9][10].Reusser and Zehe [11] proposed a temporal sensitivity analysis method to detect the periods of poor model performance and the dominant parameters.This temporal error analysis helped to better understand the model structure and its deficits [11][12][13].Lindenschmidt & Chun [14] carried out a global sensitivity analysis (GSA) to assess the impact of fluvial geomorphology on river ice cover formation.They used GSA to highlight important processes, identify sensitive regions within the parameter space, and investigate the strength of the interactions between parameters and variables.
There are several water quality models to quantify the substance transformations in river and streams, such as QUAL2E [15], QUAL2K [16], WASP7 [17], CE-QUAL-RIV1 [18], and SWAT [19,20].There are also several studies conducting sensitivity analysis with different water quality models to investigate the sensitivity of water quality variables to different parameters.Examples include the work of Drolc and Koncan [21], who studied parameter sensitivity to dissolved oxygen using the QUAL2E model and the work of Griensven et al. [22] that used SWAT to examine parameter sensitivity to water quality variables.In this study, we used WASP (Water Quality Analysis Simulation Program) [17], which has been widely used to assess the water quality of various waterbodies such as lakes, rivers, and estuaries (e.g., Franceschini & Tsai [23], Kaufman [24], Lindenschmidt [25]).The latest version, WASP7, accounts for ice-covered conditions and was selected for this study to predict algal-nutrient dynamics along the South Saskatchewan River (SSR), both in summer and under ice-covered conditions in winter.
The main objective of this paper is to determine the spatial changes of the sensitivity of different water quality variables to various parameters.The determination of such spatial changes is of significance to understand the impact of nutrient loading from tributaries and effluents from lagoons and wastewater treatment plants on different water quality variables along a river.It also helps to understand the impact of lakes and reservoirs on water quality along the river, as reportedly they increase retention of nutrients through increased deposition [26,27].Therefore, GSA was used to examine the potential impact of nutrient loading point sources (the RDR and effluents from the Saskatoon wastewater treatment plant) and the retention capacity of Lake Diefenbaker.Seasonality of parameter sensitivity is also taken into account to examine the effect of ice-covered conditions.Furthermore, this study examines the correlation between parameters and state variables and their interactions with water quality at four water quality stations along the river.

Site Description
The SSR originates at the confluence of the Bow and Oldman rivers in Alberta, Canada, flows into south-central Saskatchewan and joins with the North Saskatchewan River to flow onward as the Saskatchewan River.Lake Diefenbaker, a reservoir formed by the construction of the Gardiner and Qu'Appelle River dams, modulates the flow of the SSR to secure water for domestic, agricultural and municipal water supplies, provides flood protection and hydropower supply, and maintains aquatic and riparian ecosystems.
The SSR was divided into two major modeling domains: ‚ upper SSR (uSSR), originating at the confluence of Bow and Oldman rivers to the inlet of Lake Diefenbaker.
‚ lower SSR (lSSR), originating at the Gardiner Dam and extending to the confluence of the North and South Saskatchewan rivers at the Saskatchewan River Forks east of Prince Albert.
The RDR is a major tributary of the upper South Saskatchewan River and provides forty percent of the flow to the uSSR.The RDR is a nutrient-rich river, with a large amount of phosphorus originating from intensive agricultural activities [28] Rapid agricultural and population growth has elevated nutrient loading into the river, making monitoring and management of surface water quality more pertinent [28].
Water 2016, 8, 158 3 of 15 Lake Diefenbaker is the only reservoir along the SSR that receives a high amount of phosphorus from the uSSR inflow (3146 t TP/yr).The trophic level is classified as meso to eutrophic at the inlet [29,30].Due to sediment deposition and phosphorus retention in the lake, the water quality improves downstream of the Gardiner and Qu'Appelle dam arms, where the water quality is classified as oligo to mesotrophic [30,31].The SSR receives large lagoon effluents and treated effluent from the Saskatoon wastewater treatment plant.
To study the impact of Lake Diefenbaker and the nutrient loads from the RDR and Saskatoon wastewater treatment plant, four stations along the SSR were selected: Medicine Hat, Leader, Outlook, and Clarkboro Ferry.The locations of these sites are shown in Figure 1.
Figure 2 compares the measured concentrations of total nitrogen (TN), total phosphorus (TP), dissolved oxygen (DO), and chlorophyll-a (CHLA) at the selected stations along the SSR.Higher concentrations of TP, DO, and CHLA, but lower TN concentrations, were observed in the uSSR downstream of the RDR inlet.As reported earlier, Lake Diefenbaker is an efficient phosphorus trap and the measured concentrations show a reduction in the TP and CHLA concentrations between the inlet and outlet of the lake (difference between "b" and "c" in Figure 2).However, surprisingly, several algal bloom events were recorded in the Qu'Appelle Dam arm of Lake Diefenbaker [26] and at Outlook (difference between "c" and "d" in Figure 2).Lower CHLA and higher TN were observed in the lSSR after receiving effluents from Saskatoon.
A seasonal trend was observed for the water quality of the SSR.As Figure 2 shows, DO and nitrate (NO3) concentrations are higher in winter than in summer.CHLA and TP concentrations are typically higher in summer compared to winter.
Water 2016, 8, 158 3 of 16 [29,30].Due to sediment deposition and phosphorus retention in the lake, the water quality improves downstream of the Gardiner and Qu'Appelle dam arms, where the water quality is classified as oligo to mesotrophic [30,31].The SSR receives large lagoon effluents and treated effluent from the Saskatoon wastewater treatment plant.
To study the impact of Lake Diefenbaker and the nutrient loads from the RDR and Saskatoon wastewater treatment plant, four stations along the SSR were selected: Medicine Hat, Leader, Outlook, and Clarkboro Ferry.The locations of these sites are shown in Figure 1.
Figure 2 compares the measured concentrations of total nitrogen (TN), total phosphorus (TP), dissolved oxygen (DO), and chlorophyll-a (CHLA) at the selected stations along the SSR.Higher concentrations of TP, DO, and CHLA, but lower TN concentrations, were observed in the uSSR downstream of the RDR inlet.As reported earlier, Lake Diefenbaker is an efficient phosphorus trap and the measured concentrations show a reduction in the TP and CHLA concentrations between the inlet and outlet of the lake (difference between "b" and "c" in Figure 2).However, surprisingly, several algal bloom events were recorded in the Qu'Appelle Dam arm of Lake Diefenbaker [26] and at Outlook (difference between "c" and "d" in Figure 2).Lower CHLA and higher TN were observed in the lSSR after receiving effluents from Saskatoon.
A seasonal trend was observed for the water quality of the SSR.As Figure 2 shows, DO and nitrate (NO3) concentrations are higher in winter than in summer.CHLA and TP concentrations are typically higher in summer compared to winter.Water 2016, 8, 158 3 of 16 [29,30].Due to sediment deposition and phosphorus retention in the lake, the water quality improves downstream of the Gardiner and Qu'Appelle dam arms, where the water quality is classified as oligo to mesotrophic [30,31].The SSR receives large lagoon effluents and treated effluent from the Saskatoon wastewater treatment plant.
To study the impact of Lake Diefenbaker and the nutrient loads from the RDR and Saskatoon wastewater treatment plant, four stations along the SSR were selected: Medicine Hat, Leader, Outlook, and Clarkboro Ferry.The locations of these sites are shown in Figure 1.
Figure 2 compares the measured concentrations of total nitrogen (TN), total phosphorus (TP), dissolved oxygen (DO), and chlorophyll-a (CHLA) at the selected stations along the SSR.Higher concentrations of TP, DO, and CHLA, but lower TN concentrations, were observed in the uSSR downstream of the RDR inlet.As reported earlier, Lake Diefenbaker is an efficient phosphorus trap and the measured concentrations show a reduction in the TP and CHLA concentrations between the inlet and outlet of the lake (difference between "b" and "c" in Figure 2).However, surprisingly, several algal bloom events were recorded in the Qu'Appelle Dam arm of Lake Diefenbaker [26] and at Outlook (difference between "c" and "d" in Figure 2).Lower CHLA and higher TN were observed in the lSSR after receiving effluents from Saskatoon.
A seasonal trend was observed for the water quality of the SSR.As Figure 2 shows, DO and nitrate (NO3) concentrations are higher in winter than in summer.CHLA and TP concentrations are typically higher in summer compared to winter.The measured forcing data during the study period included daily flow, water temperature, and daily solar radiation and are presented as supplementary materials in Section A.

Model Description and Set-Up
The eutrophication module in the Water Quality Analysis Simulation Program (WASP7) was used to model water quality of the SSR.WASP7 is an enhanced MS Windows version of the original program that was developed by the U.S. Environmental Protection Agency [17,32,33].The model is continually extended and updated and has a user-friendly interface that allows ease of applicability and facilitates learning in an educational environment.Additionally, the model can be run in batch mode for sensitivity and uncertainty analysis within a Monte-Carlo framework.Therefore, for this study, the WASP model was selected in lieu of other compatible water quality models such as CE-QUAL-RIV1.The eutrophication module uses a series of mass balance equations to simulate the spatial and temporal transport and transformation of nutrients, phytoplankton, and dissolved oxygen in the aquatic environment.The variables studied in this paper are phytoplankton carbon represented by chlorophyll-a (CHLA), ammonium (NH4), nitrate (NO3), dissolved oxygen (DO), organic nitrogen (ON), organic phosphorus (OP) and orthophosphate (OPO4).
In this study, we used two WASP7-based water quality models already developed and calibrated by Akomeah et al. [34] and Hosseini et al. [35] for the upper and lower parts of the SSR, respectively.They calibrated these models using a trial-and-error approach to achieve an acceptable model performance with respect to measured values.They also validated their models with independent datasets (yielded r 2 values of 0.69 for the lSSR, and 0.68 and 0.94 for the uSSR at Leader and Medicine-Hat), demonstrating the predictability of the developed models to describe our river system.In the present study, we utilized these established water quality models to carry out a global sensitivity analysis to assess the strength of the interactions between water quality parameters and state variables.The sensitivity results were compared at four different locations to investigate the effect of the RDR, Lake Diefenbaker, and effluents of the Saskatoon wastewater treatment plant on the water quality of the river.
The setup for each model, for the uSSR and the lSSR, is briefly explained in this section, but more details of the model descriptions can be found in Akomeah [34] and Hosseini [35].The river was discretized horizontally along the longitudinal profile of the river.The uSSR is about 476 km long and was discretized using 693 segments, each 700 m long.Similarly, the 330-km long lSSR was divided into 673 segments, each 500 m in length.Morphological characteristics of the segments were estimated based on surveyed cross-sectional profiles.The boundary conditions were estimated from the nearest water quality sampling stations and are provided in the supplementary materials, Section B. The initial concentration conditions in segments were calculated by averaging and interpolating water quality concentrations between two adjacent sampling stations.The measured forcing data during the study period included daily flow, water temperature, and daily solar radiation and are presented as supplementary materials in Section A.

Model Description and Set-up
The eutrophication module in the Water Quality Analysis Simulation Program (WASP7) was used to model water quality of the SSR.WASP7 is an enhanced MS Windows version of the original program that was developed by the U.S. Environmental Protection Agency [17,32,33].The model is continually extended and updated and has a user-friendly interface that allows ease of applicability and facilitates learning in an educational environment.Additionally, the model can be run in batch mode for sensitivity and uncertainty analysis within a Monte-Carlo framework.Therefore, for this study, the WASP model was selected in lieu of other compatible water quality models such as CE-QUAL-RIV1.The eutrophication module uses a series of mass balance equations to simulate the spatial and temporal transport and transformation of nutrients, phytoplankton, and dissolved oxygen in the aquatic environment.The variables studied in this paper are phytoplankton carbon represented by chlorophyll-a (CHLA), ammonium (NH4), nitrate (NO3), dissolved oxygen (DO), organic nitrogen (ON), organic phosphorus (OP) and orthophosphate (OPO4).
In this study, we used two WASP7-based water quality models already developed and calibrated by Akomeah et al. [34] and Hosseini et al. [35] for the upper and lower parts of the SSR, respectively.They calibrated these models using a trial-and-error approach to achieve an acceptable model performance with respect to measured values.They also validated their models with independent datasets (yielded r 2 values of 0.69 for the lSSR, and 0.68 and 0.94 for the uSSR at Leader and Medicine-Hat), demonstrating the predictability of the developed models to describe our river system.In the present study, we utilized these established water quality models to carry out a global sensitivity analysis to assess the strength of the interactions between water quality parameters and state variables.The sensitivity results were compared at four different locations to investigate the effect of the RDR, Lake Diefenbaker, and effluents of the Saskatoon wastewater treatment plant on the water quality of the river.
The setup for each model, for the uSSR and the lSSR, is briefly explained in this section, but more details of the model descriptions can be found in Akomeah [34] and Hosseini [35].The river was discretized horizontally along the longitudinal profile of the river.The uSSR is about 476 km long and was discretized using 693 segments, each 700 m long.Similarly, the 330-km long lSSR was divided into 673 segments, each 500 m in length.Morphological characteristics of the segments were estimated based on surveyed cross-sectional profiles.The boundary conditions were estimated from the nearest water quality sampling stations and are provided in the supplementary materials, Section B. The initial concentration conditions in segments were calculated by averaging and interpolating water quality concentrations between two adjacent sampling stations.

Parameter Sensitivity
A GSA was conducted to study the effect of parameters on model output within their entire parameter space and to identify the interactions between parameters.Various approaches to sensitivity analysis have been used in the literature, including simple one-factor-at-a-time strategies [36], Monte-Carlo filtering [37], variance-based approaches [38], and a variogram-based framework [39,40].In this study, a Monte-Carlo Analysis (MOCA) framework was used with many model simulations, each using a different set of parameter values sampled randomly from a certain range of model parameters.For each simulation, the parameters were randomly generated from uniform distributions within the parameter ranges, which for the most part, were derived from the literature and the WASP7 model manual.The initial parameter settings are those calibrated from Akomeah [34] and Hosseini [35] and are listed in Table 1.[34] the reaeration rate was defined as 1.4 (1/day) for the entire simulation period.
To run the model automatically, the WASP7 model was linked to OSTRICH (Optimization Software Tool for Research In Computational Heuristics), an open-source and model-independent code that can automate the processes of model calibration and evaluation [41].
In this study, we conducted two analyses to assess (1) the sensitivity of "individual" model variables to model parameters, and (2) the "overall" (collective) sensitivity of the model response (based on the entire model variables) to model parameters and their interactions.For the first analysis, the simulation performance for each variable was based on the sum of squared errors between simulated and observed values.To evaluate the sensitivity of overall model response to the parameters, the sum of normalized square errors (SNSE) was used, as shown in Equation (1): where n is the number of variables and k is the number of observations considered for sensitivity analysis.C sim (i.j), C obs (i,j) , and C obs pjq are simulated, observed and mean observed values for each variable, respectively.In both of the analyses described above, those parameter sets yielding good or poor performing simulation results were classified respectively as behavioral and non-behavioral [42].The behavioral sets were deemed to be those simulations yielding the best 10% of the simulated output.The cumulative distributions of the behavioral parameters were compared to the cumulative distributions of the Water 2016, 8, 158 6 of 15 entire parameter sets (i.e., a priori parameter sets).The parameter with greater deviation between the two curves is deemed to have the greater (most sensitive) impact on the model output.The Kolmogorov-Smirnov test was used to quantify the distance between the cumulative distributions of the behavioral parameters of the a priori parameter sets.The parameter sensitivity was classified using the significance level (α) of differences between the behavioral and a priori distributions.The significance level is the level at which the null hypothesis, that the two distributions are the same, is accepted.The same top 10% of the parameter sets were used to examine the interactions between parameters and to generate correlation maps.

Nutrient and Light Limitations
Light and nutrients are the most important resources for phytoplankton growth.Phytoplankton growth will be limited under the conditions of light deficiency or inadequate amount of nutrients relative to each other.Determining the limiting factor is important for watershed management and eutrophication control approaches [3,43,44].Light limitation is a function of solar radiation, fraction of daylight, depth of the water column, and total light extinction.In WASP7, the limitation factors are calculated using the following equations: Nutrient limitation is a function of dissolved inorganic phosphorus or nitrogen and varies between 0 (complete limitation) and 1 (no limitation).

Variable Interactions
First, sensitive parameters and their ranges were identified based on the Monte Carlo simulations.To investigate the interactions between variables, Latin Hypercube sampling was used to create parameter sets, which could represent the parameter space in a small sample better than just a simple Monte Carlo [45].
In this study, we aimed to examine the potential relationship between state variable correlations and the integrity of an ecosystem.The method used in this study for measuring the structure of correlation maps was the hierarchy index (h) [46]: Water 2016, 8, 158 7 of 15 where N is the total number of variables and od (vi) or outdegree is the sum of absolute values of correlations between a variable and the adjacent variables.
In order to test this hypothesis, we compared the hierarchy index and water quality index (WQI) values.The WQI values were calculated from measured data using CCME (Canadian Council of Ministers of the Environment) water quality index formulas.

Global Sensitivity Analysis (GSA)
The GSA results for stations along the SSR are presented in Tables 2 and 3 for winter and summer.The parameter sensitivity was classified using the significance level (α) of differences between the behavioral and a priori distributions.High significance was attributed to parameters with α < 0.001, medium significance was attributed to parameters with 0.001 < α < 0.01, and no significance was attributed to α > 0.01.
Although there were some differences in parametric sensitivities at different locations along the river, similar trends in the dominant controls were found at these locations.For example, phytoplankton growth (K1C) and death (K1D) rates, TON mineralization (K71C), and phosphoruscarbon (CCHL) ratio were found to be sensitive at all the locations.Denitrification rate (K20C), reaeration rate (K2), 1/2 saturation for N-limitation on phytoplankton uptake (KMNG1), and the fraction of phytoplankton death recycled to ON (fON) did not show a large impact on variables.A smaller number of parameters were sensitive to variables simulated at the Medicine Hat station (prior to the RDR confluence) in summer and at the Outlook station (downstream of Lake Diefenbaker) in winter.
NH4, ON, and OPO4 were more sensitive under ice-covered conditions than in open water.In contrast, CHLA and DO were more sensitive in summer.The transformation of ammonium to nitrate was pronounced in summer with nitrification rate (K12C) being the most influential parameter for the NH4 and NO3 variables.The temperature dependency of nitrification was consistent with the results of the study by Shammas [47].NH4 and ON were sensitive to TON mineralization (K71C), both in summer and winter, which indicates that the bacterial decomposition of organic nitrogen to ammonium was not seasonally affected.Similarly, TOP mineralization rate (K83C) was the dominant control on OP in summer and winter.Other than phytoplankton related parameters, CHLA was found to be sensitive to the phosphorus-related parameters but not to the nitrogen-related parameters.CHLA and OPO4 showed sensitivities to various similar parameters, suggesting that a strong linage between these two variables may exist.For the DO variable, SOD shows a high sensitivity, suggesting a strong interaction between the water column and the sediments.
The results presented above were based on the Monte-Carlo analysis framework with 1000 model simulations.Additional simulations were not necessary since the model outcome was stabilized by this number of simulations runs.However, our preliminary analysis showed that 1000 simulations were sufficient to yield relatively robust results.
The uncertainty in water quality variable is described here using the coefficient of variation defined as the ratio of the standard deviation to the mean [48].In Figure 3, a comparison is made between the water quality variables distributions at two locations (1) in uSSR at Leader and (2) lSSR at Clarkboro Ferry in winter and summer.The uncertainty in uSSR is larger than the uncertainty in lSSR.NH4 in summer and OPO4 and CHLA in both winter and summer have larger uncertainties compared to the other water quality variables.

K2s
•very significant (α < 0.001); ○medium significance (0.001 <α < 0.      We studied the nutrient and light limitation results to determine the degree of light or nutrient deficiency for phytoplankton growth.The results indicate that phosphorus was the limiting factor for phytoplankton growth, except in winter when light was most limiting (see Figure 4).There is no evidence of nitrogen deficiency, which is consistent with the findings of Donald, et al. [29], who showed that phytoplankton growth in Lake Diefenbaker is phosphorus limited, specifically when light was sufficient.As expected, Outlook was the most phosphorus limited site (80% of the time) due to the phosphorus retentive effect of Lake Diefenbaker.Clarkboro Ferry was the least phosphorus limited site (71% of the time).The interactions between parameters were examined using correlation maps for the best performing parameter sets (top 10%) from the Monte Carlo parameters (as presented in Figure 5).There were relatively weak correlations between parameters at all locations along the river.Phytoplankton growth rate (K1C) was directly or indirectly related to phosphorus-related parameters (i.e., KMPG1, fOP, and PCRB) at three stations along the SSR, Medicine Hat, Leader, and Clarkboro Ferry, highlighting the sensitivity of phytoplankton to the carbon:phosphorus ratio.At Outlook, there were some interactions between nitrification (K12C) and phosphorus-carbon ratio (PCRB), TON mineralization rate (K71C), and carbon-chlorophyll ratio (CCHL), suggesting that nitrogen may be slightly more limiting at that site.The weak relationships between parameters may be a result of the large number of parameters used for our Monte-Carlo analysis (high degrees of freedom).Therefore, the parameter effects are easily compensated by other parameter settings leading to weaker interactions between parameters.

Interactions between State Variables
The most sensitive parameters were selected using the Kolmogorov-Smirnov test.The parameters were ranked in order of importance for each model (uSSR and lSSR models), as shown in Table 4. Based on the results, the nitrification rate (K12C), phytoplankton growth rate (K1C), phytoplankton death rate (K1D), phosphorus-carbon ratio (PCRB), fraction of phytoplankton death recycled to OP (fOP), half-saturation for P-limitation on phytoplankton (KMPG1) and sediment oxygen demand (SOD) were the most sensitive parameters on the overall model output.The dot plots of the parameter values versus model performance (SNSE) were used to define the optimal parameter ranges.As an example, the red ovals in the dot plots for phytoplankton death rate (K1D) (as presented in Figure 6) point to the phytoplankton death rate (K1D) range with the lowest error between simulations and measured data.
The most sensitive parameters and their optimal ranges were used for Latin Hypercube sampling.To assess the potential interaction between the variables, 2187 (3ˆ7, with 7 = number of parameters) parameter sets were selected.The simulation with the least sum of squared errors from the Monte Carlo simulations was used as the base simulation.
Figure 7 shows the correlation maps of water quality variables using a graph theory approach.Thicker lines between edges show stronger correlations between the variables.The gray and red lines are positive and negative correlations, respectively.
The correlation graphs at the lSSR at Outlook and Clarkboro Ferry illustrate two clusters; the first one includes DON, OP, OPO4, and CHLA, whereas the second includes inorganic nitrogen (NO3 and NH4).DO was weakly linked to inorganic nitrogen at Clarkboro Ferry.The results indicate that the biological production for lSSR was generally more phosphorus controlled, which is consistent with the limiting factor results in Figure 4.
oxygen demand (SOD) were the most sensitive parameters on the overall model output.The dot plots of the parameter values versus model performance (SNSE) were used to define the optimal parameter ranges.As an example, the red ovals in the dot plots for phytoplankton death rate (K1D) (as presented in Figure 6) point to the phytoplankton death rate (K1D) range with the lowest error between simulations and measured data.
The most sensitive parameters and their optimal ranges were used for Latin Hypercube sampling.To assess the potential interaction between the variables, 2187 (3^7, with 7 = number of parameters) parameter sets were selected.The simulation with the least sum of squared errors from the Monte Carlo simulations was used as the base simulation.At Medicine Hat, DON was disconnected from the other variables and two clusters (first cluster included OP, OPO4, CHLA and the second one included NO3+NH4) were both connected to DO.The relationships between variables were more complex at Leader.There were negative correlations between OPO4 and all other variables.The complexity of relationships may indicate the vulnerability of the uSSR water quality system downstream of the RDR inlet.A slight change in one of the water quality variables may affect all other variables.
The greater complexity in the correlation graph and poorer water quality (based on the observed data) at the Leader station motivated an examination of the relationship between the correlation between variables and water quality, by comparing the hierarchy index values (degree of connectivity) to the water quality indices at the four studied stations along the SSR. Figure 8 displays an inverse relationship between h values and WQIs.
The network representation of the correlation linkages between model state variable has been adapted and extended from ecological network theory, reviewed in Montoya et al. [49] and Inges et al. [50].Those reviews highlight the description of network complexity in terms of the number of linkages between variables and the strength of those linkages.In our case, the magnitude of the correlations between the variable interactions in the model simulations represents the number and strength of these linkages.Montoya et al. [49] mention that more complexity may make a network more "fragile" in terms of resilience and persistence of ecosystem variables.This may be the case for the more complex network of the uSSR shown in Figure 8, which coincides with the more polluted state of the uSSR downstream of the Red Deer River confluence.Although there are more interactions (linkages) between the variables in the uSSR network compared to the lSSR network, the interactions are weaker (generally lower correlations) imposing less of a chance that the system exists at a stable equilibrium [49].Moreover, weaker interactions tend to lead to more local equilibria that are less stable [50], which corroborates nicely with the greater equifinality of the model simulations.
Thicker lines between edges show stronger correlations between the variables.The gray and red lines are positive and negative correlations, respectively.
The correlation graphs at the lSSR at Outlook and Clarkboro Ferry illustrate two clusters; the first one includes DON, OP, OPO4, and CHLA, whereas the second includes inorganic nitrogen (NO3 and NH4).DO was weakly linked to inorganic nitrogen at Clarkboro Ferry.The results indicate that the biological production for lSSR was generally more phosphorus controlled, which is consistent with the limiting factor results in Figure 4.
At Medicine Hat, DON was disconnected from the other variables and two clusters (first cluster included OP, OPO4, CHLA and the second one included NO3+NH4) were both connected to DO.The relationships between variables were more complex at Leader.There were negative correlations between OPO4 and all other variables.The complexity of relationships may indicate the vulnerability of the uSSR water quality system downstream of the RDR inlet.A slight change in one of the water quality variables may affect all other variables.The greater complexity in the correlation graph and poorer water quality (based on the observed data) at the Leader station motivated an examination of the relationship between the correlation between variables and water quality, by comparing the hierarchy index values (degree of connectivity) to the water quality indices at the four studied stations along the SSR. Figure 8 displays an inverse relationship between h values and WQIs.
The network representation of the correlation linkages between model state variable has been adapted and extended from ecological network theory, reviewed in Montoya et al. [49] and Inges et al. [50].Those reviews highlight the description of network complexity in terms of the number of  linkages between variables and the strength of those linkages.In our case, the magnitude of the correlations between the variable interactions in the model simulations represents the number and strength of these linkages.Montoya et al. [49] mention that more complexity may make a network more "fragile" in terms of resilience and persistence of ecosystem variables.This may be the case for the more complex network of the uSSR shown in Figure 8, which coincides with the more polluted state of the uSSR downstream of the Red Deer River confluence.Although there are more interactions (linkages) between the variables in the uSSR network compared to the lSSR network, the interactions are weaker (generally lower correlations) imposing less of a chance that the system exists at a stable equilibrium [49].Moreover, weaker interactions tend to lead to more local equilibria that are less stable [50], which corroborates nicely with the greater equifinality of the model simulations.

Conclusion
This research studied parametric sensitivities to state variables of a water quality model for a large Prairie river, the SSR, in western Canada.A Monte-Carlo sensitivity analysis was carried out to study the global sensitivities of parameters on state variables at different locations along the SSR.

φ
L = light limitation factor I a = the average incident light intensity during daylight hours below the surface, (ly/day) I s = the saturating light intensity of phytoplankton, (ly/day) K e = the light extinction coefficient, (m ´1) D = depth of the water column or model segment, (m) f = fraction day that is daylight, (unitless)
d a b c d a b c D a b c d a b c d a b c d a b c d 01); blank-low significance (α > 0.01); (a) uSSR at Medicine Hat; (b) uSSR at Leader; (c) lSSR at Outlook, and (d) lSSR at Clarkboro Ferry.
d a b c d a b c D a b c d a b c d a b c d a b c d

c d a b c d a b c d a b c d a b c d a b c d a b c
(α < 0.001); -medium significance (0.001 <α < 0.01); blank-low significance (α > 0.01); (a) uSSR at Medicine Hat; (b) uSSR at Leader; (c) lSSR at Outlook, and (d) lSSR at Clarkboro Ferry.

c d a b c d a b c d a b c d a b c d a b c d a b c
significant (α < 0.001); -medium significance (0.001<α < 0.01); blank-low significance (α > 0.01); (a) uSSR at Medicine Hat; (b) uSSR at Leader; (c) lSSR at Outlook, and (d) lSSR at Clarkboro Ferry

Figure 4 .
Figure 4. Nutrient and light limitations for the best performing simulation.(a) uSSR at Medicine Hat; (b) uSSR at Leader; (c) lSSR at Outlook, and (d) lSSR at Clarkboro Ferry.

Figure 4 .
Figure 4. Nutrient and light limitations for the best performing simulation.(a) uSSR at Medicine Hat; (b) uSSR at Leader; (c) lSSR at Outlook, and (d) lSSR at Clarkboro Ferry.

Figure 4 .
Figure 4. Nutrient and light limitations for the best performing simulation.(a) uSSR at Medicine Hat; (b) uSSR at Leader; (c) lSSR at Outlook, and (d) lSSR at Clarkboro Ferry.

Figure 5 .
Figure 5. Parameter correlation: grey and red edges are positive and negative correlations, respectively.(a) uSSR at Medicine Hat, (b) uSSR at Leader, (c) lSSR at Outlook, and (d) lSSR at Clarkboro Ferry.

Figure 5 .
Figure 5. Parameter correlation: grey and red edges are positive and negative correlations, respectively.(a) uSSR at Medicine Hat, (b) uSSR at Leader, (c) lSSR at Outlook, and (d) lSSR at Clarkboro Ferry.

Figure 8 .
Figure 8.Comparison of hierarchy indices and water quality indices.

Figure 8 .
Figure 8.Comparison of hierarchy indices and water quality indices.

Table 1 .
Selected parameters for sensitivity analysis.

Table 2 .
Sensitivity analysis in winter.

Table 3 .
Sensitivity analysis in summer.

Table 2 .
Sensitivity analysis in winter.

Table 3 .
Sensitivity analysis in summer.

Table 4 .
Parameter ranking in order of significance and the optimal ranges.

Table 4 .
Parameter ranking in order of significance and the optimal ranges.