**2. Methods**

#### *2.1. Study Site*

BPL is a shallow (mean depth 3.8 m), long (30 km), and narrow (<900 m) impounded lake. BPL is a cold polymictic, wind-driven lake that readily mixes when ice-free. Air temperatures range from an average daily minimum of −17.7 ◦C in January to an average daily maximum of 26.2 ◦C in July [33]. On average the lake is ice-covered from November to April, with approximately 30% of annual mean precipitation (365.3 mm) falling as snowfall.

The lake was impounded in 1939 and the dam was most recently upgraded in 2000 [34]. Since the 1950s water to BPL has been supplemented from the South Saskatchewan River Basin. At first, water was pumped but, after construction of Lake Diefenbaker (LDief) in 1967, water has been diverted from LDief via the Qu'Appelle River Dam through the Upper Qu'Appelle River channel to BPL (Figure 1). The main inflows into BPL are through these controlled inter-basin transfers of water. Mean daily discharges at flow gauge 05JG006 (Elbow Diversion), situated 3.3 kms below LDief (Figure 1), ranged from 1.4 m3/s to 5.9 m3/s over 20 years (min 0.013 m3/s/max 14.2 m3/s; 1999–2019, Water Security Agency hydrometric data). The Upper Qu'Appelle River channel consists of two sections, a channelized section (35 km) and the meandering natural river channel (62 km). Over recent years, the provincial Saskatchewan Water Security Agency (WSA) has worked to improve channel conveyance including deepening sections and undertaking various erosion control projects. Not considering water needs due to the expanded irrigation proposals, the current channel capacity is considered sufficient to meet present and future anticipated water demands. *Water* **2022**, *14*, x FOR PEER REVIEW 4 of 20 Lake residence time is highly variable, ranging from approximately 6 to 36 months [37]. Part of this variability is due to the occasional occurrence of backflows into the lake and variability of watershed runoff among years. When flows in Moose Jaw Creek, a tributary to the Qu'Appelle River downstream of BPL, exceed 50 to 60 m3/s, water levels increase such that water downstream of BPL flows backwards through the dam into BPL ('backflow'). This occurrence is infrequent (<20% of years) and to minimize the effect on water levels and water quality in BPL the dam's gates are closed; however, water can still backflow into BPL through the dam's fish passage channel and over its radial gates. Typically, water volumes that backflow are low, although in a few years they represent a significant proportion of the lake's volume.

**Figure 1.** Left: The Upper Qu'Appelle watershed. Right: Buffalo Pound Lake (BPL), Saskatchewan, Canada. The lake has a mean depth of 3.8 m, maximum depth of 6.0 m, average width of 890 m, and **Figure 1.** Left: The Upper Qu'Appelle watershed. Right: Buffalo Pound Lake (BPL), Saskatchewan, Canada. The lake has a mean depth of 3.8 m, maximum depth of 6.0 m, average width of 890 m, and total length of 30 kms. Red circles are the water quality sampling stations and green circles are the flow gauges used in this study. Hwy #2 and #19 are official highway numbers.

total length of 30 kms. Red circles are the water quality sampling stations and green circles are the

suitable for use in a range of waterbodies. The model is laterally averaged and includes a dynamic ice model, so is capable of multi-year simulations of a long, narrow, seasonally frozen waterbody such as BPL. The longitudinal segments and vertical layers are user defined and can be variable throughout the grid. Several types of inflow can be specified with individual temperature and constituent files. Numerous hydraulic structures and withdrawals can be placed throughout the grid. The W2 model has been in continuous development since 1986 and has been applied to numerous waterbodies worldwide. A complete description of the model development, hydrodynamic and ecological equations, and model assumptions are provided in the user manual [38]. This study used CE-QUAL-W2 version 4.2.2, the most recent public version at the time of undertaking the model

A 30 m resolution digital elevation model (DEM) representing the lake's bathymetry from Terry et al. [39] was used in the current project. The bathymetry grid was extended to

*2.2. Model Set-Up* 

simulations.

flow gauges used in this study. Hwy #2 and #19 are official highway numbers.

BPL's water quality is affected by its source waters and in very wet years water quality can be influenced by a surge in nutrients, organic carbon, and dissolved ions from overland run-off [30]. The region has a semi-arid to dry-subhumid climate with long-term ratios of precipitation to potential evapotranspiration (P/PET) of between 0.5 and 0.65 [35,36]. BPL's gross watershed area is 3340 km<sup>2</sup> ; however, the actual area that contributes flow in any given year is much lower. The area that contributes runoff to BPL in median runoff years, known as the effective drainage area, is 1300 km<sup>2</sup> . The reason for the variability of area contributing runoff to BPL includes low P/PET, the generally flat landscape with numerous pothole depressions that are often hydrologically isolated or poorly connected, and differences in landscape water content among years, notably fullness of wetlands and soil moisture content. This means that in dry years BPL receives little runoff from its watershed whereas in wet years it can receive substantial amounts.

The lake is split into separate waterbodies by an old highway (causeway) and the current highway (Figure 1). Inflows to BPL occur in the upper basin. Approximately 4 km from the inflow, the lake is constricted to a 53 m gap at the old highway before entering a triangular lake section between the two highways. The lake then constricts under a 45 m long bridge before entering the main lake body.

Lake residence time is highly variable, ranging from approximately 6 to 36 months [37]. Part of this variability is due to the occasional occurrence of backflows into the lake and variability of watershed runoff among years. When flows in Moose Jaw Creek, a tributary to the Qu'Appelle River downstream of BPL, exceed 50 to 60 m3/s, water levels increase such that water downstream of BPL flows backwards through the dam into BPL ('backflow'). This occurrence is infrequent (<20% of years) and to minimize the effect on water levels and water quality in BPL the dam's gates are closed; however, water can still backflow into BPL through the dam's fish passage channel and over its radial gates. Typically, water volumes that backflow are low, although in a few years they represent a significant proportion of the lake's volume.

### *2.2. Model Set-Up*

W2 is a complex two-dimensional coupled hydrodynamic and water quality model suitable for use in a range of waterbodies. The model is laterally averaged and includes a dynamic ice model, so is capable of multi-year simulations of a long, narrow, seasonally frozen waterbody such as BPL. The longitudinal segments and vertical layers are user defined and can be variable throughout the grid. Several types of inflow can be specified with individual temperature and constituent files. Numerous hydraulic structures and withdrawals can be placed throughout the grid. The W2 model has been in continuous development since 1986 and has been applied to numerous waterbodies worldwide. A complete description of the model development, hydrodynamic and ecological equations, and model assumptions are provided in the user manual [38]. This study used CE-QUAL-W2 version 4.2.2, the most recent public version at the time of undertaking the model simulations.

A 30 m resolution digital elevation model (DEM) representing the lake's bathymetry from Terry et al. [39] was used in the current project. The bathymetry grid was extended to include the additional 4 km long upper basin (upstream of Highway 2) (Figures 1 and 2). Previous studies did not include this section, making inflows less reliable, as the flow gauge was located above the omitted lake section. This led to uncertainty about the effect that this upper basin had on flows and nutrient transport. The most upstream tip of the lake was too shallow for the boat to navigate when collecting bathymetry and is excluded from the DEM. The shoreline around this section would also complicate the W2 model grid and is above the point of inflows (Figure 1). The DEM's lake extent was created when water levels were close to the operating full supply level, although water levels can occasionally surpass this during very wet years. The W2 model grid extends 1 m vertically from the surface of the DEM to account for these high-water events. The final model grid consists of 100 longitudinal segments around 300 m and up to 28 vertical layers of 0.25 m depth.

*Water* **2022**, *14*, x FOR PEER REVIEW 5 of 20

W2 has the ability to reflect constrictions in the model's grid, which were included at the highway locations (Figure 2). W2 has the ability to reflect constrictions in the model's grid, which were included at the highway locations (Figure 2).

include the additional 4 km long upper basin (upstream of Highway 2) (Figures 1 and 2). Previous studies did not include this section, making inflows less reliable, as the flow gauge was located above the omitted lake section. This led to uncertainty about the effect that this upper basin had on flows and nutrient transport. The most upstream tip of the lake was too shallow for the boat to navigate when collecting bathymetry and is excluded from the DEM. The shoreline around this section would also complicate the W2 model grid and is above the point of inflows (Figure 1). The DEM's lake extent was created when water levels were close to the operating full supply level, although water levels can occasionally surpass this during very wet years. The W2 model grid extends 1 m vertically from the surface of the DEM to account for these high-water events. The final model grid consists of 100 longitudinal segments around 300 m and up to 28 vertical layers of 0.25 m depth.

**Figure 2.** Left: Model grid with longitudinal segmentation. The upstream area in the red circle is the newly extended section of the model. Right: New digital elevation model of BPL including the upper shallow section of lake. Depths are shown relative to operating full supply level. **Figure 2.** Left: Model grid with longitudinal segmentation. The upstream area in the red circle is the newly extended section of the model. Right: New digital elevation model of BPL including the upper shallow section of lake. Depths are shown relative to operating full supply level.

#### *2.3. Data & Calibration 2.3. Data & Calibration*

The W2 simulation period covered 6.5 consecutive years between April 2013 and December 2019. The first years were wet with flood events occurring in 2013–2015. The latter The W2 simulation period covered 6.5 consecutive years between April 2013 and December 2019. The first years were wet with flood events occurring in 2013–2015. The latter half of the study years were relatively dry.

half of the study years were relatively dry. The hydrological data used in the model were sourced from several datasets. Previous W2 models for BPL had been calibrated using data ending in 1993, after which Environment and Climate Change Canada (ECCC) discontinued the flow gauges above and below BPL. ECCC reinstated the upstream flow gauge (05JG004) in June 2015. Herein, gauge 05JG004 is referred to as the Keeler Bridge gauge (Figure 1). The ECCC downstream The hydrological data used in the model were sourced from several datasets. Previous W2 models for BPL had been calibrated using data ending in 1993, after which Environment and Climate Change Canada (ECCC) discontinued the flow gauges above and below BPL. ECCC reinstated the upstream flow gauge (05JG004) in June 2015. Herein, gauge 05JG004 is referred to as the Keeler Bridge gauge (Figure 1). The ECCC downstream gauge has not been reinstated to date; however, the WSA began recording operating logs for the BPL dam structure outflows in January 2014. The WSA conducted a nutrient mass-balance study of the length of the Qu'Appelle River between April 2013 and March 2016 and derived inflows and outflows to BPL [40].

Giving preference to the gauged data, inflows were taken from the WSA nutrient mass-balance study (1 April 2013–31 May 2015), and ECCC flow gauge 05JG004 (1 June 2015–31 December 2019). Outflows from the BPL dam were provided by the WSA nutrient mass-balance study (1 April 2013–31 December 2013), and WSA dam operations log (1 January 2014–31 December 2019). Withdrawal data were provided by the Buffalo Pound Water Treatment Plant (BPWTP) for the main intake pipes to its treatment facility. Withdrawal amounts for a smaller intake pipe supplying surrounding industries were extrapolated from 1995–2013 data, as up-to-date data could not be obtained. Backflows into BPL during a downstream flood event occurred in spring 2013 and spring 2015 and were

derived by the WSA. These were included in the model as a tributary entering the most downstream segment. ECCC in-lake water level observations—5-min intervals averaged over 24-h—were used to close the water balance using the W2 water balance tool. This tool is recommended by the developers of W2 and prevents instabilities in the model during simulation due to abrupt changes in water levels. Flows may be positive or negative and are generally added to the model as a distributed tributary entering all segments equally. Adjustment flows were greatest during spring freshet (April 2013, 2014, 2015, 2018) and large precipitation events (June 2014, July 2015) and were assumed to be principally related to ungauged runoff.

Hourly meteorological data (air temperature, dewpoint temperature, wind speed and direction) were downloaded for ECCC climate station Moose Jaw A, located approximately 40 kms south of the lake. The landscape is flat and open between the two locations, and it is expected that the weather conditions will not have differed significantly for modelling purposes. Cloud cover data were not always available and any gaps in cloud cover data at Moose Jaw A were filled in from the nearby ECCC climate station at Regina International Airport. Daily total precipitation was from the ECCC climate station 'Buffalo Pound Lake' located at the site of the BPWTP. Precipitation temperature was set at dewpoint temperature as per W2 recommendations, or zero if the dewpoint was negative. Evaporation was calculated internally by the model based on provided meteorological data and geographic location.

BPL inflow water quality constituent concentrations and water temperatures were from measurements collected at Marquis station (Figure 1). For the two backflow events, backflowing constituent concentrations and water temperatures were taken from measured data at the Buffalo Pound Dam. Precipitation water quality constituent concentrations were not modelled.

The WSA provided water quality measurement data for June 2015–December 2019 for 13 in-lake stations, plus one station at the dam outflow. The 13 in-lake stations were located within eight different model segments. Data were averaged for segments with more than one station if recorded on the same day. The BPWTP also provided weekly observations for April 2013–December 2019 from a long-term dataset. The weekly samples are collected in the early morning at the BPWTP pumping station, located on the downstream end of the lake. Water is pumped from the lake through two pipes approximately one metre off the lakebed. The water travels three kilometers to the pumping station, rising approximately 82 m in elevation. The BPWTP intake pipes are in segment 86 of the model grid. A WSA sampling station is also located in this segment.

The historical BPL WQ model was calibrated in-depth through sensitivity analyses and semi-automated calibration [32,39,41]. For this study, parameter calibration values were initially set based on the older model. With the new data added, the eight model segments corresponding with water quality stations were plotted and results visually compared against the June 2015–December 2019 WSA data. The BPWTP data were plotted in segment 86 to validate the period April 2013–May 2015. Algal taxonomy data did not exist during previous modelling, and in this study the algal rates were updated to better reflect group composition per BPWTP data. Modelled groups were diatoms, green-algae, and cyanobacteria. The organic matter decay rates, and sediment oxygen demand were also recalibrated using the new WSA longitudinal profile water quality data. Wet years and dry years were treated the same in the model. The principal coefficients for the final calibrated model are listed in Table 1. Plots showing model outputs are shown in the results section. Plots focus on model segments 5, located at the top section of the upper basin, segment 32, located near the upstream end of the main lake body downstream of the highway, and segment 86, where the BPWTP intake is located. A description of the ice module and relevant coefficients can be found in Terry et al. [41].


**Table 1.** Main coefficients used in the water quality model (previous model rates in parentheses if changed) from Terry and Lindenschmidt [32]).

\* Denotes nitrogen fixation–Dolichospermum spp. (formerly Anabaena) are one of the dominant cyanobacteria taxa in BPL.

#### *2.4. Scenarios*

Three scenarios were developed for modelling LDief releases, with each scenario being tested at flow rates of 4 m3/s, 8 m3/s, and 12 m3/s for a total of nine simulations. These rates are within the hydrological constraints of the current Upper Qu'Appelle system. Scenario one is the constant flow scenario. It simulates a constant release rate from LDief for the entire model period. This would be a situation where the augmented water transfers from LDief occurred on a permanent year-round basis. There may be practical constraints to transferring water in wet years; however, the simulation is of theoretical interest.

Scenario two is the immediate release scenario. This was designed to investigate BPL's response if LDief releases started immediately after the first high flow event occurred in 2014. Vandergucht et al.'s [40] nutrient mass-balance report found that inflows at Keeler Bridge and nutrient concentrations at Marquis increased rapidly on 7 April 2014. Flows peaked at 65.6 m3/s on 9 April 2014. Model simulated releases commenced as soon as flows returned to rates <10 m3/s on 16 April 2014. This scenario was designed to mimic a management response of transferring additional water from LDief in reaction to increased nutrient loading into BPL. Modelled transfers started when flow in the Upper Qu'Appelle subsided sufficiently in consideration of channel capacity and continued until the end of the simulation. As with the previous simulation there are practical limitations with this approach during subsequent high-flow events, but this scenario specifically focusses on the response of BPL on initiating LDief transfers.

Scenario three commences after the first 2014 high flow event but with a delayed initiation and for a defined duration, referred to as the timed release scenario. Water transfers from LDief started on 9 May 2014. This date coincides with the time when 2014 flood levels in BPL lowered to near operating full supply level. This scenario assumed that additional water from LDief could not be added to BPL until water levels were sufficiently low. Water transfers terminated on 12 March 2015, when the following spring freshet caused notable increases in flow.

The scenarios were included in the model structure as a tributary entering the same upstream segment as the main inflows, which reflects the nature of the system's hydrology. For each simulation, outflows of equal rates were added to the BPL dam outflow file to aid the water balance. All scenarios used the same tributary water temperature and constituent files.

#### Scenario Input Files

Flow, water temperature, and constituent concentrations were required for input files for simulating the water releases from LDief. Water temperatures were assumed to be the same as the temperatures observed at Marquis station. Although temperature data existed for Highway 19 (Hwy 19), located 1.9 kms of channel length below LDief's outflow (Figure 1), consistent temperature changes occurred along the length of the Upper Qu'Appelle channel. Data at the Marquis station were more representative of the temperature of water flowing into BPL regardless of the source.

As with water temperature, constituent concentrations change along the Upper Qu'Appelle channel. Nutrients in rivers and streams cycle between the water column and sediments during transit [42,43]. Sources and sinks include recycling by aquatic biota, settling and resuspension of the sediment bed, adsorption to suspended solids, and biochemical transformations such as denitrification. Reoccurring patterns in constituent changes from upstream to downstream were examined so that they could be accounted for in the simulated LDief releases. Constituent concentrations at both Hwy 19 and Marquis water quality stations were compared during periods when flows were similar at both stations, a time when overland flow and tributary contributions from the watershed would have been minimal. Transit time, as a function of flow at Elbow Diversion (Figure 1), was calculated using a logarithmic function Equation (1) produced by a one-dimensional hydrologic river model (HEC-RAS) applied previously to the Upper Qu'Appelle River for winter flow testing [44]:

$$\text{Transit Time (days)} = -1.178 \times \ln(\text{flow}) + 4.7491; \text{(R}^2 = 0.9742) \tag{1}$$

Using the transit time equation, the flow at Elbow Diversion and water quality at Hwy 19 were matched with their downstream counterparts at Keeler Bridge and Marquis. Where the absolute difference in flow was greater than 0.5 m3/s the corresponding pair of water quality observations were removed from the sample set. For the remaining pairs of water quality observations, changes in concentration were calculated. Ordinary least squares regression analysis was performed between these differences and flow at Elbow Diversion to test if the regression model could be used to estimate net change in LDief concentrations for a given flow rate. The analysis was restricted to observations after 1 June 2015, which was the date the Keeler Bridge flow gauge was reinstated. In April 2016 the WSA changed their water quality reporting laboratory. For orthophosphate (PO<sup>4</sup> <sup>3</sup>−-P), ammonium (NH4-N), and nitrate-nitrite (NO3-N + NO2-N, or NOx), many observations after this date were at or under new reporting limits of detection of PO<sup>4</sup> <sup>3</sup>−-P < 0.02 mg/L, NH4-N < 0.02 mg/L, and NOx < 0.04 mg/L. The change in concentrations between the two water quality stations could not be calculated, and these observations were removed from the paired samples for the regression analysis. A greater number of samples may have returned different relationships between the three constituents and flow.

The difference between upstream and downstream concentrations for total suspended solids (TSS), total phosphorus (TP), PO<sup>4</sup> <sup>3</sup>−-P, dissolved organic carbon (DOC), and total dissolved solids (TDS) were found to be significantly related to flow (Table 2); however, for DOC and TDS, the amount of variance explained was low. The strongest relationships were for TSS and phosphorus (P).

**Table 2.** Regression models of the net change in concentration values for constituents travelling along the Upper Qu'Appelle River, as a function of flow at the upstream Elbow Diversion gauging station. Constituents are total suspended solids (TSS), total phosphorus (TP), orthophosphate (PO<sup>4</sup> <sup>3</sup>−-P), dissolved organic carbon (DOC), total dissolved solids (TDS), ammonium (NH<sup>4</sup> -N), total nitrogen (TN) and nitrate-nitrite (NOx). Flow is average daily flow in m3/s and constituents are in mg/L.


The change in concentrations for TSS, TP, and PO<sup>4</sup> <sup>3</sup>−-P are plotted in Figure 3. The regressions show a positive significant relationship with flow. The intercepts of all three relationships were rejected at 5% significance. Soluble P is strongly sorbed to suspended solids and bed sediments [45,46]. As a result, both soluble and particulate P concentrations can be influenced by the same hydrological processes that resuspend and transport sediments at different flow rates. Changes in TP and PO<sup>4</sup> <sup>3</sup>−-P concentrations show similar outcomes to TSS in the regression models. At higher flow rates, linear relationships between suspended solids and flow can break down as mobilized sediments are depleted and increasing water volume dilutes concentrations [47]. This is evident in the pattern for TSS in Figure 3. *Water* **2022**, *14*, x FOR PEER REVIEW 10 of 20

**Figure 3.** Regression plots for the three constituents that displayed a stronger relationship (adjusted R square > 0.5) between change in concentrations during transport along the Upper Qu'Appelle River with flow rate. Flow is average daily flow in m3/s at the WSA gauging station Elbow Diversion. **Figure 3.** Regression plots for the three constituents that displayed a stronger relationship (adjusted R square > 0.5) between change in concentrations during transport along the Upper Qu'Appelle River with flow rate. Flow is average daily flow in m3/s at the WSA gauging station Elbow Diversion.

The change in concentrations for variables DOC and TDS had a significant relationship with flow, as seen by the *p*-values. However, for a physical process, the ability of the regression model to explain the variability was low based on the adjusted R square value, and not suitable for predicting the nutrient transformations for the LDief releases. TDS formation and composition is influenced by site specific biological and chemical processes The change in concentrations for variables DOC and TDS had a significant relationship with flow, as seen by the *p*-values. However, for a physical process, the ability of the regression model to explain the variability was low based on the adjusted R square value, and not suitable for predicting the nutrient transformations for the LDief releases. TDS formation and composition is influenced by site specific biological and chemical processes

such as decay of organic materials, pH, temperature, and mineral types [48]. DOC chemical composition can be altered through prolonged reservoir storage before downstream

NH4-N, total nitrogen (TN), and NOx did not have a significant flow relationship. Poor or unexpected relationships between concentrations of particulate and dissolved

The positive relationships for TP and TSS could not be incorporated into the scenario input files. The inflow constituent files in W2 require a split between organic and inorganic substances for the boundary data and the model requires inorganic suspended solids (ISS) rather than TSS. ISS were not included in the modelled constituents as observations just for the inorganic forms were not available. Likewise, variables such as Chlorophyll-a (CHLA), TP, TN, and DOC are not accepted by W2 as inputs but are derived by the model using inflow concentrations for PO43−-P, NH4-N, NOx, algae concentrations, and organic matter (OM). The OM concentrations need to be split into four compartments in the inflow files (labile and refractory components of dissolved and particulate organic matter: LDOM, RDOM, LPOM & RPOM). The organic N, P and carbon content of the OM is calculated internally by W2 based on user-specified fixed ratios. If the organic N and P fractions are known then they can be directly included in the timeseries data as OM subcompartments (e.g., LDOM-P, LDOM-N, RDOM-P, etc.). W2 does not double count the OM and OM-P/OM-N loading when added in this manner. Available laboratory results for OM did not meet the required W2 inputs. Available data included DOC, and these were used as the basis for estimating OM compartments based on knowledge of the source

The positive transformation relationship for PO43−-P could be accounted for in the scenario input files, and W2 assumes all phosphorus to be completely available for uptake in this form. The tributary inflow constituent file included PO43−-P concentrations after the regression model in Table 2 had been applied. For the remaining inflow constituents in the LDief files, TDS, NH4-N, NOx, LDOM, RDOM, LPOM & RPOM remained at Hwy 19 concentrations, while CHLA and DO used input concentrations from Marquis. No CHLA data were available for Hwy 19, and for DO we use the same assumptions as for water

forms of N with flow have been found previously [51].

water.

temperature.

such as decay of organic materials, pH, temperature, and mineral types [48]. DOC chemical composition can be altered through prolonged reservoir storage before downstream release [49], and can be highly variable during storm events [50]. Thus, the transformations and recycling occurring along a river stretch will be inconsistent.

NH4-N, total nitrogen (TN), and NOx did not have a significant flow relationship. Poor or unexpected relationships between concentrations of particulate and dissolved forms of N with flow have been found previously [51].

The positive relationships for TP and TSS could not be incorporated into the scenario input files. The inflow constituent files in W2 require a split between organic and inorganic substances for the boundary data and the model requires inorganic suspended solids (ISS) rather than TSS. ISS were not included in the modelled constituents as observations just for the inorganic forms were not available. Likewise, variables such as Chlorophyll-a (CHLA), TP, TN, and DOC are not accepted by W2 as inputs but are derived by the model using inflow concentrations for PO<sup>4</sup> <sup>3</sup>−-P, NH4-N, NOx, algae concentrations, and organic matter (OM). The OM concentrations need to be split into four compartments in the inflow files (labile and refractory components of dissolved and particulate organic matter: LDOM, RDOM, LPOM & RPOM). The organic N, P and carbon content of the OM is calculated internally by W2 based on user-specified fixed ratios. If the organic N and P fractions are known then they can be directly included in the timeseries data as OM sub-compartments (e.g., LDOM-P, LDOM-N, RDOM-P, etc.). W2 does not double count the OM and OM-P/OM-N loading when added in this manner. Available laboratory results for OM did not meet the required W2 inputs. Available data included DOC, and these were used as the basis for estimating OM compartments based on knowledge of the source water.

The positive transformation relationship for PO<sup>4</sup> <sup>3</sup>−-P could be accounted for in the scenario input files, and W2 assumes all phosphorus to be completely available for uptake in this form. The tributary inflow constituent file included PO<sup>4</sup> <sup>3</sup>−-P concentrations after the regression model in Table 2 had been applied. For the remaining inflow constituents in the LDief files, TDS, NH4-N, NOx, LDOM, RDOM, LPOM & RPOM remained at Hwy 19 concentrations, while CHLA and DO used input concentrations from Marquis. No CHLA data were available for Hwy 19, and for DO we use the same assumptions as for water temperature.

#### **3. Results and Discussion**

### *3.1. Model Calibration*

Results from the updated model generally correspond well with measured data for the parameters tested (Figure 4). Modelled TDS matched observed data well. This was expected because of the relatively conservative nature of TDS. The variability of inflow TDS within and among years is reflected in the segment 5 results but further down the lake this is muted due to travel length and mixing, as expected. Total DOC followed a similar pattern to TDS. Modelled DO concentrations were largely driven by seasonal changes in temperature. Measured and modelled nutrient concentrations did not match as well as some other parameters, notably in segment 32 for TN, and during the wet years in downstream segment 86. CHLA in segment 5 was modelled well, and the measured CHLA peaks in 2018 and 2019 were reflected by the model. CHLA measurements can be highly variable spatially and temporally, and water quality models should capture the seasonal range of values rather than the variability of individual observations. By segment 86, modelled CHLA had a consistent seasonal pattern, although the amplitude did not capture the range in measured values in all years. In particular, the model did not capture winter algal blooms in the downstream segments. Interestingly, the model results showed an early spring algal bloom occurring in segment 86 in April 2013. This was the timing of the first backflow event into BPL that brought elevated concentrations of nutrients as seen by the corresponding spike in TP and DOC at this time. TDS concentrations were lower in the backflowing water, and the dilution of in-lake concentrations can also be seen

*Water* **2022**, *14*, x FOR PEER REVIEW 12 of 20

in the plots. The second backflow event occurred in March 2015, which can again be seen on the plots of segment 86. of the older model did not require much change, the comparison of output in several segments at once allowed the model assumptions to be evaluated in greater detail. The final model was deemed suitable for running the management scenarios.

addition of the new longitudinal profile data means much of the uncertainty was removed from the previous modelling work of Terry and Lindenschmidt [32]. While the parameters

**Figure 4.** Model results for BPL for the calibration period April 2013–December 2019. Segment 5 is located upstream in the top section of lake. Segments 32 and 86 are both located in the main body of the lake. Water Security Agency observed in-lake sample data (red \*) are included from June 2015. Observed weekly data provided by the Buffalo Pound Water Treatment Plant (green +) are shown in segment 86. Inflow constituents DO = dissolved oxygen, and CHLA = Chlorophyll-a. **Figure 4.** Model results for BPL for the calibration period April 2013–December 2019. Segment 5 is located upstream in the top section of lake. Segments 32 and 86 are both located in the main body of the lake. Water Security Agency observed in-lake sample data (red \*) are included from June 2015. Observed weekly data provided by the Buffalo Pound Water Treatment Plant (green +) are shown in segment 86. Inflow constituents DO = dissolved oxygen, and CHLA = Chlorophyll-a.

*3.2. Model Simulations*  For all scenarios there is a general pattern change from the wet period, shortly after the start of the scenarios in the spring of 2014, to the dry period after around 2016. The wet period saw larger magnitude concentration changes and larger differences among the scenarios, whereas scenario results were more similar to each other during the dry period (Figure 5). For the constant flow scenario (Figure 5), of note are the peak concentrations during high flow events. At the upstream end, segment 5 was closely aligned with inflows and inflow concentrations, and the additional constant discharge had a dilution effect but did not affect the timing of the peaks. In segments 32 and 86, the scenario constituents regularly peaked at similar concentrations to the base model, yet the peaks occurred much earlier in the season and decreased more rapidly. The peak events also subsided after each freshet or high precipitation event and multiple peaks of shorter duration were evident. The maximum residence time according to W2 is 338.4 days for the baseline model. This reduces to only 82.2 days when a constant discharge of 12 m3/s is simulated. The lake is Spring freshet plus intense rain-events led to several high-flow and loading events in the spring and summer of 2014. High concentrations of TP, TN, and DOC can be seen in the plots for segment 5. TP, TN, and DOC are derived internally in W2, and useful for checking the model's overall treatment of the different sources and sinks of each key constituent. It appears that the nutrient loading in this period was overstated in the model's boundary data. The sampling frequency for the Marquis station water quality data ranged from 3 days to 1 month between April and September 2014, with high concentrations of PO<sup>4</sup> <sup>3</sup>−-P, for example, recorded in all but two of 14 samples. Stored nutrients in the watershed, riverbed, and riverbanks can be quickly mobilized once flows pass a critical threshold. This local mobilization begins on the rising limb of the hydrograph and often peaks in advance of the flow peak [52,53]. More distant nutrient sources in overland runoff and subsurface flow can arrive both before or after the flow peak dependent on the hydrological event and site [52,53]. Travel time in the Upper Qu'Appelle channel can be under two days in high flows, and nutrient spikes dissipate in BPL rapidly. The recorded nutrient concentrations around this time were unlikely to be sustained at high concentrations over the entire five months but, rather, capture several peak events over these wet months. WSA in-lake data were not available for verification of inflowing nutrient concentrations, and concentrations were therefore added as per the available Marquis observations. The

BPWTP weekly sample set indicated that the model was overstating TP, TN, and DOC as far downstream as segment 86 in that year.

The most significant outcome from the development of the BPL model is that the model appeared to handle both extended wet and dry conditions without predicted results deviating too far from actual observations. A major challenge in multi-year simulations is having only one set of coefficients to represent a range of environmental conditions. The model has been calibrated to relatively dry periods, including the original 1986–1993 period. The same coefficients were then used for model validation of two extremely wet years (April 2013–June 2015) with reasonable success, per the segment 86 comparison with the BPWTP data. Benefits of updating the model to the most recent data available were that the water quality scenarios became more relevant to current lake conditions. In addition, the extension of the model grid to include the top section of the lake, and the addition of the new longitudinal profile data means much of the uncertainty was removed from the previous modelling work of Terry and Lindenschmidt [32]. While the parameters of the older model did not require much change, the comparison of output in several segments at once allowed the model assumptions to be evaluated in greater detail. The final model was deemed suitable for running the management scenarios.

#### *3.2. Model Simulations*

For all scenarios there is a general pattern change from the wet period, shortly after the start of the scenarios in the spring of 2014, to the dry period after around 2016. The wet period saw larger magnitude concentration changes and larger differences among the scenarios, whereas scenario results were more similar to each other during the dry period (Figure 5).

For the constant flow scenario (Figure 5), of note are the peak concentrations during high flow events. At the upstream end, segment 5 was closely aligned with inflows and inflow concentrations, and the additional constant discharge had a dilution effect but did not affect the timing of the peaks. In segments 32 and 86, the scenario constituents regularly peaked at similar concentrations to the base model, yet the peaks occurred much earlier in the season and decreased more rapidly. The peak events also subsided after each freshet or high precipitation event and multiple peaks of shorter duration were evident. The maximum residence time according to W2 is 338.4 days for the baseline model. This reduces to only 82.2 days when a constant discharge of 12 m3/s is simulated. The lake is shallow, and materials easily suspend in the water column through wind stress and water flow. The additional discharge velocity from the LDief releases transported the constituents more rapidly and further through the lake model.

Outside of the peak events, TDS, DOC and TN concentrations decreased along with additional LDief flows. TN showed substantial reduction in concentrations in the later drier years at all three flow rates. TDS and DOC also showed reduction in concentrations in these same years, although results were closer to baseline.

TP demonstrated the opposite behaviour, and outside of the large 2014 and 2015 loading events, TP concentrations increased from baseline with each rise in LDief discharge. CHLA largely followed the trend of the nutrients with reduced concentrations in the wet years, when nutrients were smaller than the base model, and higher concentrations in the dry years when TP was increased by the scenarios.

The difference in behaviour of TP can be explained by the PO<sup>4</sup> <sup>3</sup>−-P boundary data. The percentage of PO<sup>4</sup> <sup>3</sup>−-P in TP according to the WSA data between April 2013–March 2016 was 4–100% at Hwy 19 (avg. 32%), and 5–86% (avg. 39%) at Marquis. The W2 model assumes all available P for uptake is of this form with PO<sup>4</sup> <sup>3</sup>−-P a large proportion of derived TP. PO<sup>4</sup> <sup>3</sup>−-P is the only constituent to which a regression model transformation was applied to estimate nutrient transformations in the Upper Qu'Appelle River. The relationship with flow was significant; however, the regression equation was not strong and based on only eight observation pairs, though the regression relationship with TP and flow was also positive (Figure 3), implying a transformation equation of some kind was

required. Conversely, if linear relationships break down at higher discharge, as per TSS, then the equation may not have applied equally at all three flow rates across time. Clearly, further work is needed to establish a method for estimating nutrient transformations along the Upper Qu'Appelle. This could be achieved through targeted sampling at different times of year, and/or model chaining or extension of the W2 lake model through the Upper Qu'Appelle channel. *Water* **2022**, *14*, x FOR PEER REVIEW 13 of 20 shallow, and materials easily suspend in the water column through wind stress and water flow. The additional discharge velocity from the LDief releases transported the constituents more rapidly and further through the lake model.

**Figure 5.** Output for five model-derived water quality constituents of interest at segments 5, 32 and 86 for the constant flow scenario. Scenario flow rates for the water transfers are 4, 8 and 12 m3/s. **Figure 5.** Output for five model-derived water quality constituents of interest at segments 5, 32 and 86 for the constant flow scenario. Scenario flow rates for the water transfers are 4, 8 and 12 m3/s.

Outside of the peak events, TDS, DOC and TN concentrations decreased along with additional LDief flows. TN showed substantial reduction in concentrations in the later drier years at all three flow rates. TDS and DOC also showed reduction in concentrations in these same years, although results were closer to baseline. TP demonstrated the opposite behaviour, and outside of the large 2014 and 2015 loading events, TP concentrations increased from baseline with each rise in LDief discharge. CHLA largely followed the trend of the nutrients with reduced concentrations in the wet years, when nutrients were smaller than the base model, and higher concentrations in the dry years when TP was increased by the scenarios. For the immediate release scenario (Figure 6), the model responded quickly to the initiation of the LDief transfers. Concentrations in segment 5 immediately returned to the concentrations modelled in the constant flow scenario, as would be expected so close to the inflow point. In segment 32 concentrations were similar, although TN peaked at slightly higher concentrations in this scenario when the transfers started—most obviously at the transfer rate of 12 m3/s. By segment 86 there were no noticeable differences in concentration between the scenarios on initiation of transfers. As per the constant flow scenario, constituent peaks decreased in segment 5 with the LDief transfers during the period of high overland runoff. This led to lower concentrations overall in the downstream segments as the model traced the water quality signal from the releases.

The difference in behaviour of TP can be explained by the PO43−-P boundary data. The percentage of PO43−-P in TP according to the WSA data between April 2013–March 2016 was 4–100% at Hwy 19 (avg. 32%), and 5–86% (avg. 39%) at Marquis. The W2 model assumes all available P for uptake is of this form with PO43−-P a large proportion of derived TP. PO43−-P is the only constituent to which a regression model transformation was applied to estimate nutrient transformations in the Upper Qu'Appelle River. The relation-For the timed release scenario (Figure 7), transfers were terminated in March 2015 prior to spring freshet. The model responded similarly to the other scenarios on commencement of the water transfers. Of most interest in this scenario was the impact that terminating the transfers ahead of the freshet had on 2015 concentrations. By pre-lowering the concentrations through the LDief transfers, and then reducing overall loading during the freshet itself by stopping the transfers, TDS, DOC, TP, and TN all peaked at far lower

ship with flow was significant; however, the regression equation was not strong and based on only eight observation pairs, though the regression relationship with TP and flow was also positive (Figure 3), implying a transformation equation of some kind was required.

concentrations that year in segment 86. Constituent concentrations did not return to higher concentration baseline levels until mid-July 2015, four-months after LDief water transfers were terminated. Segment 86 is of interest as model results are more removed from input values, and it is located at BPWTP intake, so is of most management relevance. It should be noted that concentrations returned to baseline levels much more rapidly in the upstream segments and constituents peaked at similar concentrations to the base model for the 2015 freshet. The timed release is most likely to have returned better results than the constant flow and immediate release scenarios in segment 86 that year as the reduced flow velocity, from the termination of transfers, increased residence time and allowed more materials to settle in the model. For the immediate release scenario (Figure 6), the model responded quickly to the initiation of the LDief transfers. Concentrations in segment 5 immediately returned to the concentrations modelled in the constant flow scenario, as would be expected so close to the inflow point. In segment 32 concentrations were similar, although TN peaked at slightly higher concentrations in this scenario when the transfers started—most obviously at the transfer rate of 12 m3/s. By segment 86 there were no noticeable differences in concentration between the scenarios on initiation of transfers. As per the constant flow scenario, constituent peaks decreased in segment 5 with the LDief transfers during the period of high overland runoff. This led to lower concentrations overall in the downstream segments as the model traced the water quality signal from the releases.

Conversely, if linear relationships break down at higher discharge, as per TSS, then the equation may not have applied equally at all three flow rates across time. Clearly, further work is needed to establish a method for estimating nutrient transformations along the Upper Qu'Appelle. This could be achieved through targeted sampling at different times of year, and/or model chaining or extension of the W2 lake model through the Upper

*Water* **2022**, *14*, x FOR PEER REVIEW 14 of 20

Qu'Appelle channel.

**Figure 6.** Output for five model-derived water quality constituents of interest at segments 5, 32 and 86 for the immediate release scenario. **Figure 6.** Output for five model-derived water quality constituents of interest at segments 5, 32 and 86 for the immediate release scenario.

For the timed release scenario (Figure 7), transfers were terminated in March 2015 prior to spring freshet. The model responded similarly to the other scenarios on commencement of the water transfers. Of most interest in this scenario was the impact that terminating the transfers ahead of the freshet had on 2015 concentrations. By pre-lowering the concentrations through the LDief transfers, and then reducing overall loading during The total spread of predicted concentrations at segment 86 illustrates the overall impact on water quality concentrations per scenario (Figure 8). For DOC, TP, and TN results for all scenarios and the base model appear positively skewed. This data distribution infers that most predicted concentrations were small (the box plus left tail includes 75% of the predicted concentrations). The whiskers in Figure 8 are plotted at interquartile range (IQR) ∗ 1.5, with the remaining data, beyond the whiskers, representing the infrequent high concentration loading events (e.g., freshet). TDS demonstrated increasing skewness with the duration of LDief releases and increased transfer rates.

lowed more materials to settle in the model.

the freshet itself by stopping the transfers, TDS, DOC, TP, and TN all peaked at far lower concentrations that year in segment 86. Constituent concentrations did not return to higher concentration baseline levels until mid-July 2015, four-months after LDief water transfers were terminated. Segment 86 is of interest as model results are more removed from input values, and it is located at BPWTP intake, so is of most management relevance. It should be noted that concentrations returned to baseline levels much more rapidly in the upstream segments and constituents peaked at similar concentrations to the base model for the 2015 freshet. The timed release is most likely to have returned better results than the constant flow and immediate release scenarios in segment 86 that year as the reduced flow velocity, from the termination of transfers, increased residence time and al-

**Figure 7.** Output for five model-derived water quality constituents of interest at segments 5, 32 and 86 for the timed release scenario. **Figure 7.** Output for five model-derived water quality constituents of interest at segments 5, 32 and 86 for the timed release scenario.

The total spread of predicted concentrations at segment 86 illustrates the overall impact on water quality concentrations per scenario (Figure 8). For DOC, TP, and TN results for all scenarios and the base model appear positively skewed. This data distribution infers that most predicted concentrations were small (the box plus left tail includes 75% of the predicted concentrations). The whiskers in Figure 8 are plotted at interquartile range (IQR) ∗ 1.5, with the remaining data, beyond the whiskers, representing the infrequent high concentration loading events (e.g., freshet). TDS demonstrated increasing skewness The constant flow and immediate release scenarios show a decrease in the median value and the IQR range for both TDS, DOC, and TN over the base model. This indicates that predicted concentrations over the model period were lowered by the LDief water transfers. TP displayed opposite behaviour to the other constituents, which may in part be a result of the transformation equation applied to PO<sup>4</sup> <sup>3</sup>−-P in the LDief constituent file. Overall, CHLA did not change greatly from the base model due to scenario predictions showing both increases and decreases in concentrations over the simulation period, as per the time series plots.

with the duration of LDief releases and increased transfer rates. The constant flow and immediate release scenarios show a decrease in the median value and the IQR range for both TDS, DOC, and TN over the base model. This indicates In the timed release scenario, the LDief water was released for only 306 days (12.4% of days) with the rest of the inflows the same as the base model. This shorter duration of water transfers did not alter the overall median and IQR as markedly as the other two scenarios. Segment 86 results from the timed release scenario (Figure 7) suggest that, if the water transfers had been repeated once the 2015 freshet had subsided and continued between the subsequent freshets, then the medians and IQR range may have been much lower. Such a scenario would have been similar to the immediate release scenario but with transfers stopped temporarily each year during flood periods.

Overall, the scenario results implied that adding additional transfers of water from LDief resulted in a greater number of days with lower levels of nutrients. As we anticipated, the influx of better-quality source water from LDief improved the general water quality of BPL. The constant flow scenario appeared to be the most effective at lowering total concentrations in the plotted segments at all three transfer rates.

**Figure 8.** Boxplots comparing model output at segment 86 for the base model and the nine simulation runs for April 2013–December 2019. The box represents 50% of the ranked data (interquartile range, or IQR) with the median marked as a horizontal line within the box. Whiskers represent upper and lower bounds at (IQR ∗ factor of ±1.5). **Figure 8.** Boxplots comparing model output at segment 86 for the base model and the nine simulation runs for April 2013–December 2019. The box represents 50% of the ranked data (interquartile range, or IQR) with the median marked as a horizontal line within the box. Whiskers represent upper and lower bounds at (IQR ∗ factor of ±1.5).

that predicted concentrations over the model period were lowered by the LDief water transfers. TP displayed opposite behaviour to the other constituents, which may in part be a result of the transformation equation applied to PO43−-P in the LDief constituent file. Overall, CHLA did not change greatly from the base model due to scenario predictions showing both increases and decreases in concentrations over the simulation period, as per

In the timed release scenario, the LDief water was released for only 306 days (12.4% of days) with the rest of the inflows the same as the base model. This shorter duration of water transfers did not alter the overall median and IQR as markedly as the other two scenarios. Segment 86 results from the timed release scenario (Figure 7) suggest that, if the At the same time, the increased discharge transported constituents more rapidly and further through the lake during peak events. This created shorter, more frequent peaks at the location of the BPWTP intake pipe than occurred with the base model—with the freshet nutrient peaks arriving earlier in the season. These results are consistent with the decrease in residence time and increased flow velocity from adding the additional discharge volume in the scenarios. The timed release scenario results, at segment 86 for the 2015 freshet, suggest that pausing the water transfers during periods of high flow may reduce the intensity of the nutrient peak arriving at the BPWTP intake location. It may be that trying to push as much water as one can at any time is not always effective at attaining water quality goals. In addition, pausing the water transfers during high flow events will need to occur to ensure flows remain within the Upper Qu'Appelle's channel capacity.

#### **4. Concluding Remarks**

Stakeholder engagement at the model design stage helps give confidence in the modelling process and end results. Workshops have shown that end users have a better understanding of what is being modelled, modelling possibilities and limitations, e.g., [54,55]. Stakeholders are more satisfied with output because model scenarios have been tailored to meet their own specific objectives. In this study, the provincial water management agency, the WSA, was provided with the opportunity to refine the scope of the model development and scenarios in order to address pertinent questions about water diversion strategy. For example, the WSA were keen to incorporate the extension to the model bathymetry to remove uncertainties in results brought to light in earlier studies, e.g., [32]. In turn, the enhanced in-lake water quality sampling by the WSA allowed this bathymetry to be implemented and tested in the model. The result was a more robust model that satisfies both stakeholder and scientific objectives. With this study we are bridging the needs of the end user and water quality modelling.

The modelled results show that, dependent on the timing and quantity of water transferred, the increased volume of water is predicted to decrease some water quality parameters in BPL. A flow rate of 12 m3/s was most effective at lowering total concentrations for all scenarios. BPL's residence time is highly variable and scenario results indicate that the lake will respond rapidly to increased discharge velocity. The expansion of the model grid, and inclusion of the newer profile data, removed the uncertainties existing in the original BPL model with regards to hydrology and transport. There is strong support that the timing of LDief water transfers is important. Results for the timed release scenario suggest the optimum timing would be to wait until overland runoff is at a minimum before commencing water transfers, followed by terminating transfers before runoff resurges. Constituent concentrations near the treatment plant intake are lowest during the 2015 freshet for this scenario (Figure 7). Another benefit to this approach is that the risk of flooding the Upper Qu'Appelle system during natural high runoff events is reduced.

The difference in TP behaviour with and without the transformation equation underscores the need to understand how constituents change along the Upper Qu'Appelle before releasing any transfers. Such a study would require higher-frequency sampling data along the Upper Qu'Appelle than is available at present. Options include extending the W2 model upstream to LDief, or the coupling of two or more water quality models in a river–lake chain. The simulations in this study were based on current environmental conditions and should be treated with caution for long-term planning. Climate is changing in the Prairies and the climate change signal may eventually outweigh the benefits of releasing additional water from LDief. Testing the flow management options under climate change simulations will be the next step for the BPL model.

**Author Contributions:** Conceptualization, J.T., J.-M.D. and K.-E.L.; Formal analysis, J.T.; Funding acquisition, K.-E.L.; Investigation, J.T. and J.-M.D.; Methodology, J.T. and J.-M.D.; Supervision, K.-E.L.; Visualization, J.T.; Writing—original draft, J.T.; Writing—review & editing, J.-M.D. and K.-E.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** Funding for this project was provided by the University of Saskatchewan's Global Water Futures Program.

**Acknowledgments:** The authors thank the Buffalo Pound Water Treatment Plant and the Saskatchewan Water Security Agency for providing water quality data used to build and verify the model. The University of Saskatchewan's School of Environment and Sustainability and Global Institute for Water Security also provided data and technical support. Analyses presented in this paper do not necessarily reflect the opinions or official position of the Saskatchewan Water Security Agency or the Government of Saskatchewan.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
