Novel statistical analysis illustrates the importance of flow source for extreme variation in dissolved organic carbon in a eutrophic reservoir in the Great Plains

Dissolved organic carbon (DOC) trends, predominantly showing long-term increases in concentration, have been observed across many regions of the Northern Hemisphere. Elevated DOC concentrations are a major concern for drinking water treatment plants, owing to the effects of disinfection byproduct formation, the risk of bacterial regrowth in water distribution systems, and treatment cost increases. Using a unique 30-year data set encompassing both extreme wet and dry conditions in a eutrophic drinking water reservoir in the Great Plains of North America, we investigate the effects of changing source-water and in-lake water chemistry on DOC. We employ novel wavelet coherence analyses to explore the coherence of changes in DOC with other environmental variables and apply a generalized additive model to understand predictor–DOC responses. We found that the DOC concentration was significantly coherent with (and lagging behind) flow from a large upstream mesotrophic reservoir at long (> 18-month) timescales. DOC was also coherent with (lagging behind) sulfate and in phase with total phosphorus, ammonium, and chlorophyll a concentrations at short (≤ 18-month) timescales across the 30-year record. These variables accounted for 56 % of the deviance in DOC from 1990 to 2019, suggesting that water-source and in-lake nutrient and solute chemistry are effective predictors of the DOC concentration. Clearly, climate and changes in water and catchment management will influence source-water quality in this already water-scarce region. Our results highlight the importance of flow management to shallow eutrophic reservoirs; wet periods can exacerbate water quality issues, and these effects can be compounded by reducing inflows from systems with lower DOC. These flow management decisions address water level and flood risk concerns but also have important impacts on drinking water treatability.

1 Introduction

Organic matter (OM) within lakes is derived from both internal processes (autochthonous) and the transport of land-derived organic matter (allochthonous) (Wetzel, 2001). The dissolved component of organic matter (DOM), often measured as dissolved organic carbon (DOC), is an important water quality characteristic that regulates the structure and function of lake and reservoir ecosystems (Prairie, 2008). In lakes and reservoirs used as drinking water sources, elevated DOM concentrations (i.e., DOM quantity) and changes in DOM colour and composition (i.e., DOM quality) can be problematic and are a major concern for drinking water treatment plants. At high DOM levels, water treatment is more challenging, necessitating costly DOM removal during pre-treatment (Cooke and Kennedy, 2001) and other actions to limit harmful disinfection byproduct (DBP) formation (Williams et al., 2019) and bacterial regrowth in water distribution systems (LeChevallier et al., 1996). The importance of DOC in regulating the structure and function of aquatic systems and the attention that it receives in drinking water treatment processes (e.g., Morling et al., 2017; Parr et al., 2019) have prompted many studies into the drivers of change in DOC.

Increasing DOC concentrations have been reported for freshwater lakes and reservoirs, primarily in forested regions of the Northern Hemisphere; this includes regions of eastern Canada (Couture et al., 2011; Hudson et al., 2003), the northeastern United States (Rodríguez-Cardona et al., 2015), northern Europe (Futter et al., 2014; Pärn and Mander, 2012), central Europe (Hruška et al., 2009), and the United Kingdom (Evans et al., 2005). Others have reported no significant trends in the DOC concentration (Dillon and Molot, 2005; Jane et al., 2017), that trends in flux do not accompany trends in concentration (Eimers et al., 2008), or decreasing trends over time (Rodríguez-Murillo et al., 2015). The search for common drivers has pointed to changes in acidic deposition, including both anthropogenic sulfur (Monteith et al., 2007) and nitrogen (N) deposition (Evans et al., 2008; Musolff et al., 2017), and natural cycles in sea salt (i.e., chloride) deposition (Evans et al., 2006) as being important. Climatic and hydrological factors, such as changes in temperature, precipitation, or runoff, can also play a role (Hongve et al., 2004; Weyhenmeyer and Karlsson, 2009; Björnerås et al., 2017). Moreover, trends in DOC have been linked to land management and land use (Yallop and Clutterbuck, 2009). Importantly, these drivers can operate at varying temporal and spatial scales (Clark et al., 2010), leading to individual or multiple drivers emerging as dominant in different studies.

DOC influx from catchments to lakes is regulated by local hydrology and climate as well as by landscape features, which influence DOC production and transport (Pace and Cole, 2002; Porcal et al., 2009). At the landscape scale, allochthonous DOC leached from terrestrial soils is transported to streams and lakes, and DOC flux from soils is regulated by soil moisture and flow paths from the soil to streams or lakes (Clark et al., 2010; Sobek et al., 2007). Landscape features that influence the DOC concentration in lakes include the wetland area (Dillon and Molot, 1997), the drainage ratio (catchment area : lake area) and proportion of the catchment that is open water (Kortelainen, 1993), the hydrological connectivity of water sources (Laudon et al., 2011; Schiff et al., 1997), the land use (Williams et al., 2010; Wilson and Xenopoulos, 2009), and the catchment slope and topography (Sobek et al., 2007; Xenopoulos et al., 2003).

Ultimately, it is evident that landscape attributes influence DOC exports and that landscape diversity complicates our understanding of the drivers of changing DOC. Importantly, there has been relatively little work assessing the change in aquatic DOC in regions with limited effects of acidic deposition. The Canadian Prairies Ecozone, a vast area of 1.8×10⁶ km², features flat topography, poorly developed stream networks, abundant depressional storage in pothole wetlands, and climatic variability. Periods of deluge and drought are common (Pham et al., 2009; Pomeroy et al., 2007; Vogt et al., 2018), and soils of the region often have high concentrations of DOC and SO${}_{\mathrm{4}}^{\mathrm{2}-}$, which is reflected in pothole pond chemistry (Arts et al., 2000; Miranda and Whitfield, 2024; Labaugh et al., 1987; Waiser, 2006). Importantly, approximately half of the region is hydrologically noneffective (does not contribute runoff to streams in a typical year) (Godwin and Martin, 1975; PFRA, 1983; Martin, 2001), although this varies strongly at small spatial scales, even for catchments of the same class (Wolfe et al., 2019). Low noneffective areas mean that surface flows from the local catchment are irregular and typically limited to a few weeks in spring. As such, they may play an important, albeit irregular, role in delivering nutrients and DOC to larger water bodies. Although endorheic saline prairie lakes are known for their often extremely high DOC concentrations, owing to prolonged evapoconcentration (Arts et al., 2000; Osburn et al., 2011; Waiser and Robarts, 2000), DOC concentrations in freshwater lakes of the region have received less attention, and few long time series of DOC concentration exist. Long-term change in DOC is of particular interest to drinking water treatment providers relying on surface water in the region (Imtiazy et al., 2023). Understanding these DOC patterns can support safe drinking water treatment and limit the production of disinfection byproducts like trihalomethanes and haloacetic acids. Given that the region is relatively dry and that large reservoirs are used to store water for irrigation and municipal supplies, support hydropower generation, and facilitate recreation, flow management to meet multiple objectives is important (e.g., Suen and Eheart, 2006; Yang et al., 2012).

This study investigates long-term variation in DOC concentration in Buffalo Pound Lake (Saskatchewan, Canada); this lake has a high DOC concentration for a drinking water source and is also prone to cyanobacterial blooms (Painter et al., 2022). Our objectives were to understand relationships between DOC, flows from the local catchment, managed flows diverted from upstream Lake Diefenbaker, and a suite of in-lake water chemistry parameters over a 30-year period (1990–2019) encompassing both wet conditions and severe drought, with varied water flows and management. Our analysis centres on the novel use of wavelet coherence and phase analyses to identify temporally dependent patterns between time series that are typically masked by traditional correlation methods. These relationships were further investigated with generalized additive modelling (GAM) to understand the predictor–response relationship between DOC concentration, flows, and water chemistry. We hypothesized that flows into Buffalo Pound Lake would be important determinants of patterns in the DOC concentration at longer timescales. We also assessed chemical changes at shorter timescales, focusing on total phosphorus and DOC relationships, which can be related due to allochthonous transport or the stimulation of autochthonous production via increased primary production in a lake that can be phosphorus-limited.

2 Methods

2.1 Study site

Buffalo Pound Lake is a shallow, eutrophic reservoir located near the headwaters of the Qu'Appelle River drainage system in southern Saskatchewan, Canada (50.65127° N, 105.508225° W; Fig. 1). Buffalo Pound Lake is the drinking water source for ∼25 % of the province's residents, servicing the cities of Regina (population of 249 000) and Moose Jaw (population of 34 000). The lake is the site of a provincial park and a popular recreation spot for activities that include fishing, swimming, and camping; however, agriculture is the dominant land use in the region. The climate in this region is semi-arid with long, cold winters and short, warm summers. The mean annual precipitation is 320 mm (McGowan et al., 2005), although it varies between dry and wet climate phases, and the average annual temperature is ∼3 °C (Haig et al., 2021). Buffalo Pound Lake was created by damming a natural lake's outflow to the Moose Jaw River to raise water levels (beginning in 1939). The lake is long and narrow (∼29 km × 1 km) with an average depth of 3 m (max depth of 5.6 m) and storage capacity of 9×10⁷ m³ (BPWTP, 2020), with a short residence time (8 months to 2.5 years; BPWTP, 2016; Vogt et al., 2018). Its shallow depth and exposure to regular wind-induced mixing result in a polymictic system that only periodically establishes thermal stratification (Finlay et al., 2019).

Figure 1 Map of Buffalo Pound Lake, Saskatchewan, Canada, highlighting its gross (grey) and effective (green) drainage areas, including upstream Lake Diefenbaker and the Qu'Appelle River. Red diamonds denote Water Survey of Canada (WSC) streamflow gauging stations in the Buffalo Pound Lake catchment with complete or reconstructed records over the study period (see Sect. 2.3 and Table S1). The yellow circle shows the location of two submerged intake pipes that deliver raw lake water 3 km to the Buffalo Pound Water Treatment Plant (BPWTP). The short (Ridge Creek) and segmented (Iskwao Creek) blue lines are an artefact of their transience. Eyebrow Lake, which is not a true lake, is ephemeral and periodically contributes ungauged flow to Buffalo Pound Lake.

Buffalo Pound Lake has two major water sources: the local catchment and upstream diversions from Lake Diefenbaker (Fig. 1), a reservoir with a volume of 9×10⁹ m³, approximately ∼100 times greater than that of Buffalo Pound Lake. Importantly, the local catchment water quality and inflows from Lake Diefenbaker differ with respect to multiple attributes, including the DOC and sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$) concentrations. For example, measured local catchment DOC and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are on the order of 20 mg L⁻¹ and > 1000 mg L⁻¹, respectively, whereas the respective concentrations in Lake Diefenbaker are typically < 5 mg L⁻¹ (Saskatchewan Water Security Agency, unpublished data) and < 100 mg L⁻¹ (John-Mark Davies, unpublished data).

Lake Diefenbaker lies outside of the natural catchment boundary for Buffalo Pound (BPWTP, 2020) and receives a continuous flow sourced from the South Saskatchewan River, originating in the Rocky Mountains. Water releases from Lake Diefenbaker maintain Buffalo Pound lake levels within an established operating range (30 cm) which has some seasonal variability. During the summer months, due to the importance of the water supply and recreational interests, levels are typically managed within the upper 10 cm of the operating range, but they are reduced by about 15–20 cm in winter. The second major water source to Buffalo Pound is the local catchment (3310 km², of which only 38% is effective), which is dominated by agriculture with nutrient-rich runoff that contributes to eutrophic conditions (Finlay et al., 2019; Hammer, 1971). These inflows are generally intermittent, occurring primarily during spring snowmelt and (infrequent) extreme summer precipitation. During wet years, which lead to higher catchment flows to Buffalo Pound, managed releases from Lake Diefenbaker are much lower than in dry years. In contrast, dry years, with little or no catchment flow, require larger releases from this upstream reservoir to maintain lake levels within the operating range. Groundwater flows are not thought to be important in terms of volume or DOC contributions to the lake.

2.2 Data

2.2.1 Water chemistry

Water samples were collected and analyzed on weekly, semi-monthly, or monthly intervals at the Buffalo Pound Water Treatment Plant (BPWTP) between 1990 and 2019. Samples were collected from raw lake water that was pumped 3 km from near the Buffalo Pound Lake outlet to the BPWTP via two submerged intake pipes (Fig. 1). The BPWTP laboratory is accredited to the ISO/IEC 17025 standard and adapts its analysis methods using the Standard Methods for the Examination of Water and Wastewater (APHA, 2012). The accredited BPWTP laboratory did change instruments over the long history of this data set; however, new instruments were cross-validated and subjected to the testing requirements associated with accreditation (Blair Kardash, personal communication, 2020). Here, we briefly describe the current analytical methods (Table S2 in the Supplement). Prior to 2015, DOC samples were analyzed via a nondispersive infrared detector on a Tekmar Dohrmann Phoenix 8000 TOC (total organic carbon) Analyzer. Since 2015, DOC concentrations have been analyzed on a General Electric (GE) Sievers M5310 C TOC Analyzer with a GE-patented Selective Membrane Conductometric design. Filtration (e.g., for DOC) was to a pore size of 0.7 µm. Sulfate and nitrate (NO${}_{\mathrm{3}}^{-}$) concentrations were determined via ion chromatography using a Dionex ICS-1100 ion chromatograph. Ammonium (NH${}_{\mathrm{4}}^{+}$) concentrations were determined colourimetrically by Nessler's reagent addition after samples were buffered in boric acid solution and distilled. Chlorophyll a (Chl a) pigment concentrations were extracted with acetone and measured via spectrophotometry, but they were not corrected for the presence of pheophytin. Total phosphorus (TP) analysis was carried out via the molybdenum blue method at 690 nm on a spectrophotometer after digestion. Soluble reactive phosphorus (SRP) concentrations were analyzed using the ammonium molybdate method and then measured at 690 nm on a spectrophotometer.

The frequency of water chemistry analyses permitted the use of monthly averages (n= 360) across the 30-year time series with few missing data points (0.28 %–3.6 %) for the seven parameters described above. Statistical analyses (described below) required a complete record (i.e., 360 monthly averages of weekly to monthly samples collected over 30 years), necessitating (very limited) imputation. Because many water chemistry parameters were highly variable over decadal scales, replacing missing values with the (seasonal) mean or median consistently produced values outside of the reasonable range for a given month and year based on visual inspection. To overcome this, the k-nearest neighbour regression approach (Altman, 1992; Fix and Hodges, 1951) was used, replacing missing values with the mean of k= 1 nearest neighbours.

2.2.2 Streamflow measurements

The streamflow monitoring network in the Buffalo Pound Lake catchment is operated by the Saskatchewan Water Security Agency (WSA) and comprises four streamflow gauges with sometimes incomplete records for 1972–2019 (Fig. 1, Table S1). Lake Diefenbaker outflow (Q_LD) and Buffalo Pound Lake inflow (Q_BP) are gauges at the upstream and downstream reaches of the Qu'Appelle River that represent inflow from Lake Diefenbaker and upstream inflow to Buffalo Pound Lake, respectively. Flows from Ridge Creek (Q_RC) and Iskwao Creek (Q_IC) are tributaries sourced from the local catchment that flow into the Qu'Appelle River between Lake Diefenbaker and Buffalo Pound Lake (see Sect. 2.3). Ephemeral flows from the effective areas running the length of either side of the catchment (Fig. 1) represent another (ungauged) source of water to the lake.

2.3 Streamflow reconstruction scheme

Conventional hydrological modelling is not an easy task in the Buffalo Pound Lake catchment due to the poorly defined drainage system distributed over a relatively flat terrain (e.g., Mekonnen et al., 2015; Pomeroy et al., 2005). This forms complicated hydrological conditions, under which the hydrologic responses heavily depend on landscape feature detail (Fang et al., 2007; Van der Kamp et al., 2003; Spence, 2010) and are subject to nonlinear hysteresis (Shook and Pomeroy, 2011). With few exceptions (Clark and Shook, 2022), current hydrological models do not represent dominant hydrological processes in prairie landscapes and fail to reproduce the observed runoff. Here, we attempted to maximize the use of observed gauge information (Table S1) by implementing a data-driven reconstruction scheme. This approach allows us to consider three water sources in our analysis: diversions from Lake Diefenbaker (Q_LD); upstream inflows to Buffalo Pound Lake, representing a mix of local catchment flows and Q_LD (Q_BP); and flows originating from the local catchment (Q_LC).

Considering the mass conservation and ignoring other potential losses, the upstream inflow to Buffalo Pound Lake (Q_BP) can be described as follows:

$$\begin{array}{}\text{(1)}& Q_{\mathrm{BP}}=Q_{\mathrm{LD}}+Q_{\mathrm{RC}}+Q_{\mathrm{IC}}+Q_{\mathrm{U}}-{\displaystyle \frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}}.\end{array}$$

We used available flow information at Ridge Creek (Q_RC) to estimate missing flows at Iskwao Creek ($Q_{\mathrm{IC}}^{\ast })$. Likewise, we used Ridge Creek to estimate flows for the ungauged portion of the catchment upstream of the Buffalo Pound Lake inflow ($(Q_{\mathrm{U}}-\frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}{)}^{\ast }$), an area that includes Eyebrow Lake (Fig. 1), which periodically attenuates flow. While records for Q_LD and Q_RC are complete, we used a set of linear and nonlinear transfer functions, including polynomial, exponential, and inverse exponential functions, along with power and tangential functions and their linear combinations (resulting in a total of 420 functional mappings) to link $Q_{\mathrm{IC}}^{\ast }$ and $(Q_{\mathrm{U}}-\frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}{)}^{\ast }$ to Q_RC. We note that the nature of landscape hydrologic response in the lake catchment is significantly different between frozen (December–March) and ice-free (April–November) months; therefore, we considered separate transfer functions for warm and cold seasons. The 1972–1992 data were used for model development, and the models were compared and falsified based on the Bayesian information criterion, which penalizes the number of model parameters. For Iskwao Creek, the following transfer functions performed well during the training period:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& Q_{\mathrm{IC}}^{\ast }\\ & =\left\{\begin{array}{ll}\mathrm{0.3112}{Q}_{\mathrm{RC}}^{\mathrm{0.4537}}& \text{if cold season,}\\ \mathrm{0.4957}{Q}_{\mathrm{RC}}+\mathrm{0.1185}& \text{if warm season and }\\ & Q_{\mathrm{RC}}>\mathrm{0},\\ \mathrm{0}& \text{if warm season and }\\ & Q_{\mathrm{RC}}=\mathrm{0}.\end{array}\right.\end{array}\end{array}$$

Here, $Q_{\mathrm{IC}}^{\ast }$ and Q_RC are in cubic metres per second. Similarly, the total upstream ungauged flux, in cubic metres per second, was estimated based on Q_RC as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\left(Q_{\mathrm{U}}-{\displaystyle \frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}}\right)}^{\ast }\\ & =\left\{\begin{array}{ll}\mathrm{0.01292}{Q}_{\mathrm{RC}}^{\mathrm{3}}-\mathrm{0.303}{Q}_{\mathrm{RC}}^{\mathrm{2}}& \text{if cold season, }\\ +\mathrm{1.249}{Q}_{\mathrm{RC}}-\mathrm{0.0956}& \\ -\mathrm{0.05958}{Q}_{\mathrm{RC}}^{\mathrm{2}}& \text{if warm season.}\\ +\mathrm{2.77}{Q}_{\mathrm{RC}}-\mathrm{0.1463}& \end{array}\right.\end{array}\end{array}$$

Here, the storage term $\left(\frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}\right)$ accounts for changes in landscape storage, including at Eyebrow Lake. The efficiency of estimated $Q_{\mathrm{IC}}^{\ast }$ and ${\left(Q_{\mathrm{U}}-\frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}\right)}^{\ast }$ was investigated by calculating the coefficient of determination (R²) for observed flows at Iskwao Creek and upstream of Buffalo Pound Lake during separate training (1972–1992) and testing periods. The R² measures the percentage of described variance and was 0.71 during the training period and 0.68 during the testing period of 2008–2011 for $Q_{\mathrm{IC}}^{\ast }$. For the total inflow to the lake (Q_BP), R² was similarly high: 0.87 during training and 0.78 during testing (1993–1995).

Finally, because Q_BP is dominated by Q_LD, we used a third metric as a means of assessing the role of flows from the local catchment (Q_LC). In this instance, we approximate Q_LC, including runoff from the catchment area downstream of the gauge at the Buffalo Pound Lake inflow, according to the following:

$$\begin{array}{}\text{(4)}& Q_{\mathrm{LC}}\cong \left\{\begin{array}{ll}{Q}_{\mathrm{DC}}& \mathrm{if}{Q}_{\mathrm{RC}}+Q_{\mathrm{IC}}\\ & +{\left(Q_{\mathrm{U}}-\frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}\right)}^{\ast }<\mathrm{0},\\ Q_{\mathrm{RC}}+Q_{\mathrm{IC}}& \mathrm{if}{Q}_{\mathrm{RC}}+Q_{\mathrm{IC}}\\ +{\left(Q_{\mathrm{U}}-\frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}\right)}^{\ast }& +{\left(Q_{\mathrm{U}}-\frac{\mathrm{\Delta }S}{\mathrm{\Delta }t}\right)}^{\ast }\ge \mathrm{0}.\\ +Q_{\mathrm{DC}}& \end{array}\right.\end{array}$$

Here, we scale combined flows from the Ridge and Iskwao creeks to estimate flows from the effective catchment area that lies downstream of the upstream inflow to Buffalo Pound Lake (Q_DC), including areas along either side of the lake (Fig. 1). While this approach can only serve as an estimate, as any attenuation of flows in the ungauged portion of the Qu'Appelle River can be a mix of source waters from Lake Diefenbaker and the local catchment, it nonetheless provides a proxy for the flow contribution from the local catchment.

2.4 Statistical analyses

There were two major goals of our statistical analysis. First, we aimed to understand frequency-dependent relationships between the DOC concentration and a suite of environmental predictors that include upstream flows (Q_LD, Q_BP, and Q_LC) and in-lake water chemistry (SO${}_{\mathrm{4}}^{\mathrm{2}-}$, TP, SRP, Chl a, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$) using wavelet coherence analysis (Grinsted et al., 2004; Sheppard et al., 2016, 2017; Walter et al., 2021). Wavelet-based methods can be used to measure synchrony and coherence between environmental and biological variables and to investigate relationships that are not readily detected by conventional correlation methods (Reuman et al., 2021). Specifically, wavelet coherence analysis can reveal the degree of association between two variables, while the phase relationships provide information on the patterns of behaviour between these two variables. Variables that are in phase or antiphase are analogous to a positive or negative correlation, respectively, while the phase relationships also provide information on where there is a time lag between the variables (e.g., Walter et al., 2021). Importantly, this approach allows for investigation across different timescales, as it is possible that these relationships manifest in different ways at different timescales. We also note that these phase relationships do not identify dependence among variables. Variables identified as significantly coherent with DOC (α= 0.05) by wavelet coherence analysis were subsequently used to achieve our second objective, wherein generalized additive modelling (GAM) was used to investigate predictor–response relationships between the DOC concentration, flows, and in-lake water chemistry. Predictor–response relationships were assessed using GAM because we anticipated complex nonlinear trends between DOC and predictor variables that would be masked by traditional (e.g., Mann–Kendall) methods or would require some degree of a priori trend selection or functional form (e.g., linear and generalized additive models). This approach allows for investigation of the relationships between each predictor and the response variable, but we note that we use this technique as a tool to analyze complex relationships, not as a forecasting tool or to infer that these predictors are drivers of DOC. All statistical analyses were completed using monthly averages (n= 360; see Sect. 2.2.1).

All statistical analyses were carried out in R (version 4.1.3; R Core Team, 2022), including “Wavelet approaches to synchrony” (wsyn) for wavelet analysis (version 1.0.4; Reuman et al., 2021), base (R Core Team, 2022), and spatstat.core (version 2.4.2; Baddeley et al., 2015) for wavelet plotting. Generalized additive models were fit using the R package mgcv (version 1.8.40; Wood, 2017), and graphics were plotted with ggplot2 (version 3.3.6, Wickham, 2016), gratia (version 0.7.3; Simpson, 2021), and patchwork (version 1.1.1; Pedersen, 2020). The graphics in Fig. 1 were plotted using the R package mapcan (McCormack and Erlich, 2019). The gross and effective drainage area shapefiles used in Fig. 1 were obtained from the Prairie Farm Rehabilitation Administration, and lake and river geospatial data were retrieved through the Government of Canada's Open Government License (Statistics Canada, 2006, 2011). Further details on wavelet analyses and GAMs are provided below.

2.4.1 Wavelet coherence analyses

We applied the continuous Morlet wavelet transform W_σ(t) (Addison, 2017) to each time series following the approaches outlined by Sheppard et al. (2016, 2017, 2019). The wavelet transform W_σ(t) uses small localized wavelike functions to transform time series signals into a more useful form (Addison, 2017). For this study, scaling was applied such that one wavelet oscillation was equal to 2 months (σ= 2), as a 2-month period is the highest-frequency fluctuation that can be identified in monthly time series. Wavelets were generated across a range of timescales, from 2 to 120 months (10 years). Wavelet transforms require scalloping to remove poorly estimated values at the tails of time series (Addison, 2017). For a 30-year time series of monthly observations, scalloping limits wavelet periods to a maximum of ∼120 months, or approximately one-third of the length of the full time series (e.g., Fig. S1 in the Supplement).

Coherence between two variables is a measure of the strength of association between those variables in a timescale-specific way that is not confounded by lagged or phase-shifted associations (Reuman et al., 2021). Accordingly, wavelet coherence quantifies the degree to which two time series have correlated magnitudes of oscillation and consistent phase differences through time, as a function of timescale (Walter et al., 2021), with a magnitude ranging from 0 (no relationship) to 1 (perfect coherence). The “wavelet mean field” normalization method (Sheppard et al., 2016) was used to measure wavelet coherence because its coherence magnitude increases both with increasing synchrony between time series and when oscillations in time series at time t and timescale σ have a similar phase (direction and time lag). This method also permits one to test for significant coherence across multiple timescales. Two timescale bands were tested for significance: short timescales (≤ 18 months) were selected based on Buffalo Pound Lake intra-annual/seasonal dynamics and potential lagged relationships between DOC and environmental predictors; long timescales (> 18 months – i.e., up to 120 months) were selected based on multiyear to decadal patterns observed in DOC concentration.

Phase (direction and time lag) was also investigated for environmental predictors that were significantly coherent with the DOC concentration. Coherent variables may be in phase (positively correlated) or antiphase (negatively correlated) and are typically time-lagged (positively or negatively). To understand phase difference relationships between significantly coherent DOC–predictor pairs, we computed the average phase φ across the corresponding timescale band (Walter et al., 2021). Because phase φ is an angular measurement, cosine- or sine-transforming this value provides information about how close the relationship is to being in phase [cos(φ)] and whether the time-lagged relationship between time series tends to be positive or negative [sin(φ)] (Walter et al., 2021). Cosine transformation assigns in-phase relationships (φ=0) to 1, antiphase relationships ($\mathit{\phi }=\pm \mathit{\pi }$) to −1, and quarter-phase (i.e., time-lagged) relationships ($\mathit{\phi }=\pm \mathit{\pi }/\mathrm{2}$) to 0. Sine transformation assigns both in-phase (φ=0) and antiphase ($\mathit{\phi }=\pm \mathit{\pi }$) relationships to a value of 0 because they exhibit no time lag. When a change in DOC occurs in advance of a predictor variable ($\mathit{\phi }=-\mathit{\pi }/\mathrm{2}$), the relationship is lagged negative, whereas when a change in DOC lags behind a change in a predictor ($\mathit{\phi }=\mathit{\pi }/\mathrm{2}$), the relationship is lagged positive.

To test for significant coherence, all time series were transformed using standard optimal Box–Cox normalization prior to wavelet transformation (Sakia, 1992). Box–Cox transformation improves normality and ensures that variability in individual time series is not dominated by extreme values (Sheppard et al., 2019). Significance testing for wavelet methods relies on Fourier surrogate techniques where time series are normally distributed (Schreiber and Schmitz, 2000), so fair comparisons and statements of significance can only be made if underlying data are normally distributed (Sheppard et al., 2019). Box–Cox transformation removes the linear trend for each time series and rescales the variance to 1, producing transformed times series with a mean of 0 and approximately normal distributions. Fourier-transform-based methods for generating surrogate coherence data sets can be used to test for statistically significant coherence relationships between wavelet-transformed variables (Sheppard et al., 2017). These Fourier-transformed data retain the original characteristics (e.g., temporal autocorrelation) of the time series and test whether coherence values are likely to occur under the null hypothesis that no actual coherence or phase relationships exist between variables. For each DOC–predictor time series pair, 10 000 surrogate randomizations were run to facilitate more accurate significance testing results and reduce variability in repeat runs (Reuman et al., 2021; Sheppard et al., 2017).

2.4.2 Generalized additive models

Predictors identified as significant by coherence analysis were used to model the DOC time series using a GAM. This approach was chosen because it can account for nonlinearity in trends, whereas other methods are limited to identifying increasing or decreasing linear (monotonic) trends (e.g., (seasonal) Mann–Kendall test) or require a priori selection of the functional form of trends in time series or selection from less flexible polynomial models (e.g., parametric linear or generalized linear models) (Finlay et al., 2019; Simpson, 2018). The use of splines by GAMs is also advantageous in that it can reduce bias issues and over-fitting at the tails of data sets, which is a common problem with polynomial models (Finlay et al., 2019). In practice, a GAM takes the following general form:

$$\begin{array}{}\text{(5)}& y_i={\mathit{\beta }}_{\mathrm{0}}+f_{\mathrm{1}}\left(x_{i\mathrm{1}}\right)+f_{\mathrm{2}}\left(x_{i\mathrm{2}}\right)+\mathrm{\dots }+f_n\left(x_{in}\right)+{\mathit{ϵ}}_i,\end{array}$$

where y_i is the response variable; β₀ is the model intercept; x_i1, x_i2, and x_in are covariates; f₁, f₂, and f_n are nonparametric smoothing functions; and ϵ_i represents independent N(0,σ²) random errors (Wood, 2017). All predictors (i.e., f(x_i)) were estimated using thin plate regression splines and penalized using restricted maximum-likelihood–based smoothness selection procedures (Wood, 2011). The initial basis dimension (k) of each smooth function was checked following the procedure described in Pya and Wood (2016): if initial k was deemed too low (i.e., if the k index < 1 and estimated degrees of freedom was close to k^′), a larger basis size was used and the model was refitted. The Tweedie distribution (Tweedie, 1981) was used as the model conditional distribution, as histograms of most of the environmental predictors showed deviation from normality and resembled some form of the gamma distribution, which is also included in the Tweedie family of distributions and identified when the Tweedie power parameter, p, is 2. In our model, p was 1.99, indicating that the gamma distribution was a good fit; however, selecting the Tweedie distribution a priori avoided the potential of erroneously choosing gamma as the conditional distribution. Model fit was assessed through qualitative inspection of quantile–quantile residuals, residuals vs. linear predictor, histogram of residuals, and observed vs. fitted values plots. Uncertainty in the estimated DOC trend was simulated using 10 000 iterations from the posterior distribution of the fitted values. The simulated trends are consistent with the estimated trend but include the uncertainty in the estimates of the spline coefficients (Finlay et al., 2019). The posterior simulation involves drawing 10 000 samples from the multivariate normal distribution and then deriving the difference between the peak and minimum DOC for each sample (trend). The upper and lower 2.5 % probability quantiles of the distribution of 10 000 differences in trend for each year form a 95 % confidence interval on the difference estimated from the fitted trend (Finlay et al., 2019).

3 Results

3.1 Flow and water chemistry parameters

Flow sources varied over the 30-year period with no distinct long-term or multiyear patterns (Fig. 2a–c). Buffalo Pound Lake inflows (Q_BP) and Q_LD each showed seasonal patterns, peaking around late spring to early summer and decreasing during winter months when water demands were lower and the lake was covered with ice (Fig. S2). Control structures from Lake Diefenbaker to Buffalo Pound Lake permit flows to continue during winter months when the lakes are ice-covered. Local catchment flows (Q_LC) typically peaked in late spring to early summer and did not contribute water to Buffalo Pound Lake in winter months or dry years. Peaks in Q_LD, Q_BP, and Q_LC are typically associated with spring freshet, but large rainfall events can also contribute considerable flow within the region. Streamflow averaged 3.2±2.6 (range: 0–12.4) m³ s⁻¹ at Q_LD and 3.6±2.8 (range: 0.1–18.9) m³ s⁻¹ at Q_BP (Table S3), highlighting much lower catchment flows which averaged 0.7±2.3 (range: 0–25.4) m³ s⁻¹. In years where Q_LC was low or absent, particularly in the late 1990s and early 2000s, Q_BP generally followed patterns of water release from Lake Diefenbaker, and fluctuations in Q_BP were comparable to those of Q_LD. Several years showed evidence of wetter-than-normal conditions (e.g., 1997 and 2014) with short-term episodes of considerable discharge from the local catchment. These episodes correspond to peaks in Q_BP and were generally associated with low Q_LD flows (Fig. 2).

Figure 2 Monthly time series of flows of (a) Q_LD, (b) Q_BP, (c) Q_LC, (d) SO${}_{\mathrm{4}}^{\mathrm{2}-}$, (e) TP, (f) SRP, (g) Chl a, (h) NO${}_{\mathrm{3}}^{-}$, and (i) NH${}_{\mathrm{4}}^{+}$ concentrations from 1990 to 2019. Note that the y-axis scales for the Q_LD, Q_BP, and Q_LC plots are different and that all nitrogen concentrations are given as milligrams of N per litre (mg N L⁻¹).

The DOC concentrations fluctuated considerably over the 30-year observation period, ranging from 3.3 to 12.4 mg L⁻¹ (Fig. 3, Table S3). The mean (± standard deviation) DOC concentration was 6.8±1.8 mg L⁻¹. Several notable fluctuations occurred over multiyear periods. For example, after sustained high concentrations from 1997 to 2000 (mean 9.5 mg L⁻¹), DOC fell to 4.4 mg L⁻¹ by 2004 and remained near or below average until the mid-2010s. Between 1990 and 2000, there was a gradual increase in the DOC concentration where, alongside intra-annual increases and decreases, DOC concentrations increased from 3.9 mg L⁻¹ in 1990 to 9.3 mg L⁻¹ in 2000. Other notable changes occurred at shorter timescales: in 1991, DOC concentrations were ≤5.0 mg L⁻¹ from January to July, spiked to 10.5 mg L⁻¹ in August, and fell to 5.5 mg L⁻¹ by December. The mid-2010s also saw a sharp increase in DOC concentration followed by below-average levels within a 4-year span.

Figure 3 Monthly Buffalo Pound Lake DOC concentrations (1990–2019) by (a) year, showing the pattern in the DOC concentration over the full time series, and (b) day of year, showing the intra-annual variation in the DOC concentration.

Water chemistry was also highly variable over the 30-year period (Fig. 2d–i, Table S3), with evidence of dramatic changes in the chemistry of Buffalo Pound Lake over time. For example, SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations ranged from 56.8 to 340 mg L⁻¹ and were relatively stable from 1990 to 2009 before rapidly rising up to 250 mg L⁻¹ between 2011 and 2017 (Fig. 2d). By the end of 2019, SO${}_{\mathrm{4}}^{\mathrm{2}-}$ returned to levels observed during 1990–2010. Nitrate also showed two distinct patterns. After averaging 0.20±0.28 mg L⁻¹ from 1990 to 1999 and peaking at 1.5 mg L⁻¹ in 1997, NO${}_{\mathrm{3}}^{-}$ concentrations were nearly an order of magnitude lower from 2000 to 2019, averaging 0.05±0.08 mg L⁻¹ (Fig. 2h). Total phosphorus, SRP, Chl a, and NH${}_{\mathrm{4}}^{+}$ concentrations rose rapidly in 1991 (Fig. 2e–g, i), concomitant with the rise in DOC in that year, but this was transient. In the few years after 1991, chemical concentrations were average or below average, followed by increases in concentrations into the early 2000s. The mid-2010s also showed elevated levels of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, TP, Chl a, and NH${}_{\mathrm{4}}^{+}$, concurrent with elevated DOC concentrations and Q_LC during this time.

3.2 Wavelet coherence and phase relationships

Coherences between DOC and nine environmental predictors at short (≤ 18-month) and long (> 18-month) timescales were highly variable (Table 1, Fig. S3). At short timescales, coherence magnitudes ranged from 0.08 to 0.50 (median of 0.15), whereas at long timescales, they ranged from 0.05 to 0.74 (median of 0.54); this is consistent with potential bias in wavelet coherence, where greater coherence magnitudes tend to be returned at long timescales (Walter et al., 2021). Four DOC–predictor pairs were significant at short timescales (α=0.05 significance level): SO${}_{\mathrm{4}}^{\mathrm{2}-}$, TP, Chl a, and NH${}_{\mathrm{4}}^{+}$. At long timescales, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, TP, and Q_LD were significantly coherent with DOC (α= 0.05 significance level). Significant coherence relationships at these timescales indicate that the magnitude of changes that occur in DOC concentrations and the environmental predictors are closely related. Phase relationships, expounded upon next, describe the time lags in these changes. Pairs of variables that are in phase will increase (and decrease) together, and a time lag can tell us if changes in one variable precede the other. Antiphase behaviour indicates that one variable is increasing while another decreases.

Table 1 Coherence and phase relationships between DOC and environmental predictors at short (≤18-month) and long (> 18-month) timescales. Mean coherence p values are calculated over the reported timescale band for each significant relationship. Significant relationships shown using bold text and are denoted by ${}^{***}$ (p < 0.001), ${}^{**}$ (p < 0.01), and ^* (p < 0.05). Relationships that are marginally nonsignificant (p < 0.1) are shown using italic font and are denoted by ${}^{****}$ Cos(φ) and sin(φ) are transformations of the average phase φ over the reported timescale band; cos(φ) describes how close the relationship between DOC and each parameter is to being in phase, whereas sin(φ) focuses on whether the time-lagged relationship between DOC and each predictor tends to be positive or negative. Because phases are angular measurements, φ ranges between −π and π, and cos(φ) and sin(φ) range between −1 and 1. Relationships are identified as in phase when cos(φ) ≈1, antiphase when cos(φ) $\approx -\mathrm{1}$, and quarter-phase when cos(φ) ≈ 0. When sin(φ) $\approx -\mathrm{1}$, the relationship is time-lagged positive (i.e., a change in DOC precedes a change in the predictor), whereas when sin(φ) ≈1, the relationship is time-lagged negative (DOC lags behind the predictor). When relationships are perfectly in phase or antiphase (i.e., not time-lagged), sin(φ) = 0.

Phase relationships spanned the range of possible values (0 to ±π) (Fig. S4). For DOC–predictor pairs with significant coherence, most relationships were approximately ($-\mathit{\pi }/\mathrm{4}\mathit{<}\mathit{\varphi }\mathit{<}\mathit{\pi }/\mathrm{4}$) at both short and long timescales (Table 1). At short timescales, DOC was strongly in phase with TP, Chl a, and NH${}_{\mathrm{4}}^{+}$ ($\mathrm{0}\mathit{<}\mathit{\varphi }\mathit{<}\mathit{\pi }/\mathrm{4}$), whereas it tended to lag behind fluctuations in SO${}_{\mathrm{4}}^{\mathrm{2}-}$ ($-\mathrm{3}\mathit{\pi }/\mathrm{4}\mathit{<}\mathit{\varphi }\mathit{<}-\mathit{\pi }/\mathrm{2}$), indicating that DOC increases (or decreases) when these predictors increase (or decrease). Interestingly, the phase relationship between DOC and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ shifted to strongly in phase ($-\mathit{\pi }/\mathrm{4}\mathit{<}\mathit{\varphi }\mathit{<}\mathrm{0}$) at long timescales, while TP remained strongly in phase with DOC ($\mathrm{0}\mathit{<}\mathit{\varphi }\mathit{<}\mathit{\pi }/\mathrm{4}$). Unlike SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and TP, fluctuations in DOC were antiphase with and tended to lag behind fluctuations in Q_LD ($-\mathrm{3}\mathit{\pi }/\mathrm{4}\mathit{<}\mathit{\varphi }\mathit{<}-\mathit{\pi }/\mathrm{2}$), indicating that DOC decreases follow increases in Q_LD. Significant coherence and phase relationships described here provide insight into the timescale and time-lag relationships between DOC and the environmental predictors. Next, analogies to traditional positive and negative correlation are discussed, using GAM.

3.3 Generalized additive modelling

Five environmental predictors (SO${}_{\mathrm{4}}^{\mathrm{2}-}$, TP, Chl a, NH${}_{\mathrm{4}}^{+}$, and Q_LD) found to be significantly coherent with the DOC concentration, as identified by wavelet coherence analyses at either short (≤18-month) or long (> 18-month) timescales, were included in a simple GAM. A thin plate regression spline for each covariate explained 56 % of the deviance in the DOC concentration over the 30-year period, with an adjusted R² of 0.50 (effective degrees of freedom = 16.9). Sulfate, TP, NH${}_{\mathrm{4}}^{+}$, and Q_LD all explained the significant variation in the DOC concentration (p < 0.001 for NH${}_{\mathrm{4}}^{+}$ and Q_LD; p < 0.0001 for SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and TP). The Chl a smooth term was the only nonsignificant predictor (p= 0.15; Fig. 4a). Total phosphorus and Q_LD showed linear positive and negative relationships with the DOC concentration, respectively (Fig. 4b, e). At NH${}_{\mathrm{4}}^{+}$ concentrations <∼0.1 mg L⁻¹ (61 % of the observed NH${}_{\mathrm{4}}^{+}$ data), NH${}_{\mathrm{4}}^{+}$ had an approximately linear positive effect on the DOC concentration but did not have much predictive power at concentrations > 0.1 mg L⁻¹, as shown by the increasingly large confidence interval (CI; shaded region) at higher NH${}_{\mathrm{4}}^{+}$ levels (Fig. 4c). The relationship between DOC and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (Fig. 4d) was generally positive but was more variable than the other environmental predictors. For example, at an SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentration of ∼150–200 mg L⁻¹, increases in SO${}_{\mathrm{4}}^{\mathrm{2}-}$ were concurrent with decreases in DOC; however, at SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations < 150 mg L⁻¹ (71 % of the observed SO${}_{\mathrm{4}}^{\mathrm{2}-}$ data), DOC tended to increase with increases in SO${}_{\mathrm{4}}^{\mathrm{2}-}$. At higher concentrations of SO${}_{\mathrm{4}}^{\mathrm{2}-}$, TP, and NH${}_{\mathrm{4}}^{+}$ and high Q_LD, the GAM lost predictive power (large CI for the right tails of each smooth), but it performed well in the lower to moderate ranges (Fig. 4). Visually, the GAM provided a good fit between observed and fitted DOC at a concentration <∼ 7 mg L⁻¹, but considerable variation was apparent in the model fit at concentrations above a threshold of 8 mg L⁻¹ (Fig. S5). Fitting the GAM to the DOC time series (Fig. 5) revealed that using these predictors to estimate the DOC concentration was useful for identifying multiyear trends (e.g., in the 1990s and mid-2010s), but the model suffered from over-fitting intra-annual fluctuations, such that DOC minima and maxima are not always captured (Fig. 5).

Figure 4 Partial effects of (a) SO${}_{\mathrm{4}}^{\mathrm{2}-}$, (b) TP, (c) Chl a, (d) NH${}_{\mathrm{4}}^{+}$, and (e) Q_LD for the GAM fitted to the DOC time series with significantly coherent environmental predictor covariates. The x axes show observed values for each environmental predictor. The y axes show the partial effects of thin plate regression spline smooths (black lines) for each environmental predictor. Grey shaded regions are the 95 % confidence interval. The rug (the inset on the x axes) displays the distribution of observed values. Chlorophyll a was the only nonsignificant covariate (p=0.15); all other environmental predictors were significant at p < 0.001.

Figure 5 The DOC concentration time series (black points) overlain by the GAM fit (GAM model using SO${}_{\mathrm{4}}^{\mathrm{2}-}$, TP, Chl a, NH${}_{\mathrm{4}}^{+}$, and Q_LD; blue line). The shaded blue area is the 95% confidence interval of the GAM model fit using 10 000 samples drawn from the posterior multivariate normal distribution. This model explained 56 % of the deviance in the DOC concentration with an $R_{\mathrm{adj}}^{\mathrm{2}}$ of 0.50.

4 Discussion

4.1 DOC concentration and flow conditions

Semi-arid regions are subject to high inter-annual variation in moisture conditions. Buffalo Pound Lake is a shallow reservoir in a DOC-rich landscape (Labaugh et al., 1987; Waiser, 2006) with high climatic variability (Pham et al., 2009; Pomeroy et al., 2007; Vogt et al., 2018) where prolonged periods of drought are often followed by multiple high-precipitation years. This leads to very high variability in DOC, with a range of ∼10 mg L⁻¹ over the 30-year observation period (Fig. 3). Rapid changes in DOC were sometimes observed, likely related to the lake's short residence time (8 months to 2.5 years; BPWTP, 2016; Vogt et al., 2018) that makes the system highly responsive to changing inputs. In regard to our first hypothesis, the flow source appears to have an important influence on DOC. Wavelet analyses demonstrated that flow from the upstream reservoir (Q_LD) was coherent with DOC concentrations in Buffalo Pound Lake at long timescales, with changes in the DOC concentration lagging behind changes in flows. This suggests that “managed” flows from Lake Diefenbaker have a flushing (diluting) effect, which is consistent with the generally lower DOC concentrations of this lake compared with local catchment runoff (Saskatchewan Water Security Agency, unpublished data). Thus, we fail to reject the null hypothesis (that Q_LD is not related to DOC).

We expected that Q_LC would be strongly related to lake DOC for the following reasons: prairie catchment runoff tends to be high in solutes, including DOC, salinity, and nutrients (e.g., Casson et al., 2019; Fasching et al., 2019), and prairie wetlands become important flow sources during periods of high hydrologic connectivity (Ali and English, 2019; Nachshon et al., 2014) and can have extremely high DOC concentrations (sometimes in excess of 100 mg L⁻¹) (Arts et al., 2000; Waiser, 2006). Instead, coherence of Q_LC with DOC was low (0.10) at short timescales and moderate (0.51) at long timescales, with neither being significant (Table 1). Statistically, the smaller and intermittent contribution of the local catchment to flow makes signals challenging to ascertain without much longer time series. Specifically, Q_LC, although sometimes large, is infrequent (Fig. 2) due to the ephemeral nature of these sources and their tendency to have no or minimal flow in drier years. There is also some potential for the attenuation of upstream catchment flows and for unknown lags in these flows reaching the lake, which may impair statistical detection. In contrast, Q_LD is large, even in years with high Q_LC. In addition, Q_LC is also estimated from smaller areas of gauged flow and, therefore, has greater uncertainty than Q_LD. In addition, total flows into Buffalo Pound Lake are not coherent with DOC (at short or long timescales; Table 1). Because Q_BP combines flows from Lake Diefenbaker and parts of the local catchment, this indicator may obscure relationships by integrating different mechanisms by which the two flow sources impact chemistry (e.g., greater allochthonous inputs from the local catchment and more dilute waters from Lake Diefenbaker).

Ultimately, the lack of statistical significance between Q_LC and lake DOC does not negate the likely importance of catchment DOC inputs to Buffalo Pound Lake. Due to the fact that there are only two flow sources to the lake and our analysis revealed a flushing effect of Q_LD, we would not expect this flushing to occur without the input of water of a different chemical character from the catchment. The wettest multiyear period in the long-term record (2010–2015) coincided with a prolonged period of DOC increase in Buffalo Pound Lake (Fig. 3), suggesting that a pulse of inflows from the local catchment influenced DOC and that the local catchment is likely an important source of allochthonous DOC when the Q_LC magnitude is large. Existing streamflow data from the local catchment indicate much higher DOC concentrations than in Buffalo Pound Lake. The Ridge Creek and Iskwao Creek DOC concentrations are each on the order of 20 mg L⁻¹ (Saskatchewan Water Security Agency, unpublished data), while water diverted from Lake Diefenbaker is similar to the minimum observed concentrations in Buffalo Pound Lake (< 5 mg L⁻¹; Saskatchewan Water Security Agency, unpublished data). Further, observed increases in lake DOC due to autochthony and other in-lake processes are modest (Baron, 2023): DOC concentrations measured at eight locations along the length of Buffalo Pound Lake from 2016 to 2019 showed a 1.2 mg L⁻¹ increase from near the lake's inflow to near its outflow. In contrast, DOC concentrations from all lake locations over this period ranged from 3.7 to 8.1 mg L⁻¹ (Baron, 2023), As such, increases in DOC from autochthonous sources are unlikely to be the driving factor in the observed across-year increases in DOC. It is important to emphasize that Q_LD is reduced when Q_LC is high. This is a management-driven relationship that stems from the need to maintain stable water levels and limit downstream water transfers during wet periods (due to flooding effects). However, reducing the inputs from Lake Diefenbaker (Q_LD) when the inputs from the local catchment are high reduces the likelihood that higher-DOC waters will be flushed.

4.2 DOC and chemical predictor relationships

Dissolved organic carbon was not the only chemical variable that underwent dramatic change within Buffalo Pound over the study period. Indeed, 6-fold variation in the concentration of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ was observed (Fig. 2d). At short timescales (≤18 months), DOC lagged behind SO${}_{\mathrm{4}}^{\mathrm{2}-}$ on average, but this relationship was strongly in phase at long timescales (> 18 months), suggesting that within-year or seasonal DOC and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ dynamics may differ; however, the patterns in these chemical variables become more synchronous at longer timescales. The rapid rise and fall in DOC and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations in the 2010s may be linked to both internal and external processes. The majority of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ fluxes to lakes come from soils and wetlands in the catchment, often in mineral form or bound in allochthonous organic matter containing oxygen-bound S (e.g., as ester-SO₄ or other organo-SO₄ species) (Couture et al., 2016), helping to explain coherence between SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and DOC. Much like prairie catchments have high DOC, they can also have high SO${}_{\mathrm{4}}^{\mathrm{2}-}$, again, due to soil and wetland pools. Pothole ponds have a wide range of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations (median levels of 163 (mean: 519) mg L⁻¹), with some ponds reaching concentrations of 5500 mg L⁻¹ (Miranda and Whitfield, 2024), typically much higher than in many other regions, owing to sulfur-rich geologic sources and the potential for evapoconcentration in surface waters. Prairie soils can be similarly high in SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (Fennell and Bentley, 1998). In contrast, annually measured SO${}_{\mathrm{4}}^{\mathrm{2}-}$ samples at Lake Diefenbaker from 1990 to 2019 showed median concentrations of only 58 mg L⁻¹ and maxima of < 90 mg L⁻¹ (John-Mark Davies, unpublished data), which are values near the low range of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ seen in Buffalo Pound Lake. The much higher median Ridge Creek SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations (estimated at ∼1094 mg L⁻¹; unpublished data) combined with below-average Q_LD (and higher Q_LC) in the early 2010s under more “natural” flow conditions mean that catchment SO${}_{\mathrm{4}}^{\mathrm{2}-}$ sources must have contributed to the higher Buffalo Pound Lake SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations during the latter part of the observation record. Extreme flooding across the prairies in 2011 could have created an exceptional fill-and-spill (Shaw et al., 2012) scenario (see Nachshon et al., 2014; Shook and Pomeroy, 2012), during which water containing elevated levels of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and DOC was transported via runoff from spilling pothole ponds to Buffalo Pound Lake. Similar mechanisms were likely important in 2014 and 2015, 2 years with periods of high precipitation in the region (WSA, 2018) and above-average catchment inflows, when SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and DOC concentrations further increased. There is also the potential that groundwater may be infusing minerals, such as SO${}_{\mathrm{4}}^{\mathrm{2}-}$, into the lake (BPWTP, 2022), and analogous patterns have been observed during multiyear wet periods that have led to the activation of subsurface flow pathways that transport solutes from the soil to surface water (Draude, 2025). Prairie soils and lake sediments rich in gypsum (CaSO₄) and pyrite (FeS₂) deposited during glaciation 10 000 years ago may also play a role. The reoxidation of lake sediments rich in H₂S, CaSO₄, FeS₂, and other sulfides can also increase surface water SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations (Holmer and Storkholm, 2001); moveover, in freshwater sediments, high rates of oxidation of reduced S compounds can shift the sediments from net sinks to sources of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (Holmer and Storkholm, 2001), particularly if groundwater SO${}_{\mathrm{4}}^{\mathrm{2}-}$ sources to the lake are not high.

Lake DOC concentrations in cold regions tend to vary both seasonally and inter-annually. Our results highlight in-phase relationships of DOC and TP, Chl a, and NH${}_{\mathrm{4}}^{+}$ (Table 1). Seasonal cycles alter the rates of biological DOC production and hydrological export (Clark et al., 2010), and DOC influxes often peak during spring snowmelt, whereas biological DOC production peaks in summer (Buffam et al., 2007; Clark et al., 2010; Laudon et al., 2004). In cold regions where lakes are ice-covered for a portion of the year, DOC exclusion from ice can be important for driving winter metabolism and DOC concentrations (Guo et al., 2012; Kurek et al., 2022), whereas photo- and biodegradation processes (Hansen et al., 2016) and autochthonous DOC production (Baron, 2023) may dominate during the open-water season. In many lakes, the seasonal variation in the DOC concentration can be 1–2 orders of magnitude greater than the long-term rate of change over decadal timescales (Clark et al., 2010). Within Buffalo Pound Lake, we observed both high inter- and intra-annual variability in the DOC concentration, with fluctuation as high as 7 mg L⁻¹ within a single year.

Total phosphorus, Chl a, and NH${}_{\mathrm{4}}^{+}$ concentrations were strongly in phase with DOC at short timescales, indicating that increases in these three variables are accompanied by increases in DOC. Total phosphorus remained strongly in phase at long timescales (Table 1). Possible mechanisms for these relationships include (1) internal production and transformations and (2) synchronous inputs of nutrients and allochthonous DOC from the local catchment, as expected from SO${}_{\mathrm{4}}^{\mathrm{2}-}$ behaviour. We anticipate that the latter plays a more prominent role, given the modest DOC increases as water moves through the lake (above). Using SO${}_{\mathrm{4}}^{\mathrm{2}-}$, TP, NH${}_{\mathrm{4}}^{+}$, Chl a, and Q_LD as predictors in a GAM describing DOC explained 56 % of the deviance in the DOC concentrations observed over this 30-year data set, where all predictors were significant except for Chl a. This result presents evidence that the DOC concentration in Buffalo Pound Lake is strongly linked to in-lake nutrient and solute chemistry in addition to its major flow source, Lake Diefenbaker. We attempted to incorporate time-lagged relationships of DOC with SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and Q_LD into the GAM approach; however, a feature of wavelet coherence analysis is that phase relationships can only be determined on average over the timescale of interest (J. Walter personal communication, 2022), and thus specific time lags (e.g., 6 months or 24 months) cannot be accurately measured with this approach. The true time-lag relationship between Q_LD and DOC also likely varies depending on flow magnitude. To assess concerns over potential collinearity between DOC and TP, we ran a GAM with SRP in place of TP but found no discernable differences in model results.

Strong in-phase coherence between DOC and TP at both short and long timescales and the strong linear effect of TP on DOC (Fig. 4e) suggest that there is synchrony in the processes influencing DOC and TP in Buffalo Pound Lake. As with SO${}_{\mathrm{4}}^{\mathrm{2}-}$, prairie ponds can also be important phosphorus (P) sources to downstream lakes during periods of elevated catchment flow, particularly because ponds high in SO${}_{\mathrm{4}}^{\mathrm{2}-}$ tend to have greater sediment P release rates (Jensen et al., 2009). Spring TP concentrations measured at ∼150 ponds across the prairie region in 2019 revealed mean TP concentrations in excess of 500 µg P L⁻¹ (McFarlan, 2021); when water is mobilized from the catchment through wetland fill-and-spill behaviour (Shaw et al., 2012), Q_LC is elevated and increased nutrient export from the catchment to the lake is expected. Periods of water flow from the catchment also prompt lower Q_LD, as less water is diverted from Lake Diefenbaker, during the natural flow events. This would reduce the dilution effects of Lake Diefenbaker inputs for both DOC and TP. Sediment processes may also support this coherence. Larger pulses of internal P released from sediments during transient stratification events could be linked to mid-summer algal blooms (D'Silva, 2017; Painter et al., 2023), potentially supporting pulses of greater autochthonous DOC production. In addition, anoxia and organic P are important to internal P loading in lakes further downstream (Wauchope-Thompson et al., 2024) which could similarly support in-phase DOC and TP relationships at shorter timescales.

The Chl a concentration, a measure of algal abundance, was not significant in the GAM but showed in-phase coherence with DOC concentrations at short timescales, suggesting a seasonal effect associated with cyanobacterial/algal biomass in the lake. Inefficient conversion of light energy to organic molecules during photosynthesis means that a portion of the organic materials produced by algal cells are released as DOM/DOC to ambient water over the lifetime of the algae (Fogg, 1966; Myklestad, 1995). Extracellular DOC release, via the overflow model (Fogg, 1966, 1983; Williams, 1990) or the passive diffusion model (Bjørnsen, 1988; Fogg, 1966), has been shown to be an important source of autochthonous DOC, where up to 5 %–35 % of fixed organic carbon may be release immediately as DOM, often with high proportions of DOC (Carlson and Hansell, 2015). Other mechanisms for the release of autochthonous carbon, such as release by grazers and following sedimentation, are also likely important (Keller and Hood, 2011). Additionally, DOC concentrations can increase by 1–2 mg L⁻¹ between the Buffalo Pound Lake inflow and outflow in years dominated by flow from Lake Diefenbaker (Baron, 2023), consistent with autochthonous production and the potential for extracellular release.

Observed in-phase coherence between DOC and NH${}_{\mathrm{4}}^{+}$ at short timescales and the significance of NH${}_{\mathrm{4}}^{+}$ as a predictor of DOC in the GAM may be related to several mechanisms. The lake is sometimes N-limited (Swarbrick et al., 2019), suggesting that increased NH${}_{\mathrm{4}}^{+}$ could support increased productivity as well as autochthonous production. This temperate lake is also highly seasonal, with major summer blooms (Painter et al., 2022, 2023), followed by long, cold winter periods when the lake is ice-covered and NH${}_{\mathrm{4}}^{+}$ dynamics undergo rapid change (Cavaliere and Baulch, 2020). The early-winter phase is marked by increasing NH${}_{\mathrm{4}}^{+}$ concentrations associated with decreased autotrophy, oxygen depletion, and organic matter mineralization. We observed several years during which DOC concentrations decreased in the early-winter period. The late-winter phase shows a decline in NH${}_{\mathrm{4}}^{+}$ likely associated with biotic uptake. Seasonal variation in the DOC : DON ratios are also apparent (Fig. S6). Higher C : N ratios occur under ice cover and then decrease beginning around spring ice-off (∼100th day of year), reaching a minima in late summer, concurrent with seasonal maximum Chl a concentrations.

4.3 Flow management and effective water treatment

Flood risk and drought underpin water management concerns in the Prairie region; cyclic wet–dry phases mean that these systems frequently alternate between years of having too much or too little water associated with elevated flood and drought risk, respectively. Examples include extensive flooding in 2011 and 2014, a major drought from 1999 to 2004, and extreme drought conditions in recent years (2021–2023). Natural and managed flow scenarios associated with wet–dry cycles each present challenges for Buffalo Pound Lake water quality and for the Buffalo Pound Water Treatment Plant (BPWTP). Our analysis highlights the importance of flows from Lake Diefenbaker as a means of “flushing” the lake and contributing to lower DOC concentrations. While our analysis did not yield significant relationships between Q_LC and DOC at long timescales, this may be partly attributed to the intermittent nature of this flow source. Given that there are only two flow sources in this system, we argue that the flushing mechanism identified for Q_LD would not likely exist if years with high local catchment flows (natural flow scenario) were not contributing to high DOC and poor water quality. These conditions can take years to dissipate. For example, during a prolonged drought from 1999 to 2005, specific ultraviolet absorbance at 254 nm normalized to DOC concentration (SUVA₂₅₄) values were very low (∼0.9–1.6 L mg C⁻¹ m⁻¹ over this period), but SUVA₂₅₄ exceeded 2.5 L mg C⁻¹ m⁻¹ after a very wet 2011 that resulted in extensive flooding and then remained above 1.5 L mg C⁻¹ m⁻¹ into 2019 (Baron, 2023; Fig. S7).

There are important trade-offs here that include balancing the provision of water from Buffalo Pound Lake to meet the needs of ∼300 000 people. While Q_LD has a flushing effect, diverting additional water during wetter periods when more source water originates from the local catchment, to manage DOC concentrations (and support drinking water treatment), can be expected to exacerbate downstream flooding and is currently avoided. It should be noted that dry conditions that prompt elevated flows from Lake Diefenbaker to Buffalo Pound Lake also represent a trade-off, as this diversion of water means that less flow is available for hydroelectric power generation via water exiting Lake Diefenbaker to the South Saskatchewan River at the Gardiner Dam. Diversions to maintain Buffalo Pound Lake levels during dry phases have also been associated with worsening algal blooms in Buffalo Pound Lake, in terms of magnitude (Painter et al., 2022), and coincide with autochthonous DOC production in the lake (e.g., Baron, 2023). Periods of elevated DOC under either flow scenario can create added challenges and costs for the BPWTP to maintain safe drinking water standards. Disinfection byproduct formation is a major concern for the BPWTP (Williams et al., 2019), particularly when the DOC concentration is high. Under these circumstances, the BPWTP has used a pre-chlorination step to reduce algal growth and prevent rising floc (organic materials that coalesce and resist coagulation) that can impact treatment capacity by accumulating on filtration media (Painter et al., 2022). Pre-chlorination has typically been used in summer months when algae are abundant in Buffalo Pound Lake (Painter et al., 2022); however, this can make it difficult for the BPWTP to meet regulatory limits for trihalomethanes in years in which DOC is also elevated (BPWTP, 2016). Untreated DOC can also contribute to poor taste and odour problems (Matilainen et al., 2011).

Environmental flows represent a framework that could be used to better manage the competing needs for drinking water and the hydrological functions of the reservoir, as well as its critical fish and wildlife habitat. The variable magnitude, timing, and water quality of flows from two water sources and the impacts on DOC underscore some of the complexity that would accompany the development of environmental flow rules for Buffalo Pound Lake. Maintaining stable water levels within the lake to limit valley slumping and erosion is important (SERM, 2001), and competing regional priorities for freshwater, including hydropower production, flood control, and irrigation, are also key factors in water management (WSA, 2013; Wheater and Gober, 2013). Water quality concerns are broad and include increasing nutrient inputs to agricultural lands; agricultural drainage of wetlands, which can have major impacts on catchment streamflow (Spence et al., 2022); effluent disposal practices on agricultural land; leaking septic fields; and pollution from tourism and recreational activities. Ultimately, concurrent efforts to consider environmental flow rules and advance source-water protection could help ensure that the maximum benefit is attained from current investments in upgraded water treatment.

5 Conclusions

Our novel analysis of a rare long-term observational record for a prairie reservoir demonstrates the relationship of DOC with the flow source and several water chemistry variables. Both the water quality indicators, including DOC, and flow source show high variability at both long and short timescales. Total phosphorus and NH${}_{\mathrm{4}}^{+}$ concentrations are synchronous with DOC and in phase, an observation which could reflect synchronous inputs (driven by high solute inputs from the local catchment), autochthony, and seasonal nutrient and organic matter dynamics (across ice-covered and open-water seasons). Sulfate, which was also highly variable over the 30-year period, was synchronous and in phase with DOC, revealing the close link between solutes and DOC sourced from the local catchment. Most critically for management, the flow source (strongly influenced by climatic and hydrologic variability) is an important control on the DOC concentration via the diluting effects of an upstream reservoir. This reservoir, with its waters primarily originating in mountain headwaters, plays an important role in flushing Buffalo Pound Lake, contributing to lower DOC concentrations at longer timescales (>18 months). In this system, where wet periods stimulating more flow from the local catchment lead to less water being diverted from Lake Diefenbaker to the Qu'Appelle River, this also suggests that catchment source waters can contribute to poorer-quality lake water. If the prairies experience higher precipitation that yields more catchment runoff under a changing climate, higher allochthonous inputs of DOC, as well as inputs of nutrients and other solutes, can be expected to further degrade water quality in the lake, making drinking water treatment increasingly difficult. Importantly, management to control the water level and flood risk can be in tension with managing flows for water quality, due to the need to reduce flushing inflows to manage the water level and downstream flooding. Advancing environmental flow rules may help balance competing goals in this water-insecure region that is already subject to major climatic and water stressors.