Multidecadal Trends in Organic Carbon Flux Through a Grassland River Network Shaped by Human Controls and Climatic Cycles

Grasslands are among the most disturbed systems on Earth due to extensive agricultural development and river regulation. To define how these factors impact the movement of carbon (C) through river networks, we quantified multi‐decadal, seasonal patterns of dissolved organic C (DOC) fluxes through a grassland river network. Despite urbanization and effluent release in recent decades, we did not observe trends in DOC fluxes over 20 to 40 year monthly sample records. High interannual variability in DOC flux was linked to climate teleconnections, most strongly the Pacific Decadal Oscillation. Human impacts were more apparent on the seasonality of flux, since spring export dominated DOC flux in most sub‐watersheds, but seasonal changes were absent at the river mouth due to downstream flow control. In grassland regions, human impacts on DOC flux are complex, and subtle long‐term changes may be difficult to identify in contemporary monitoring records due to increasing hydro‐climatic variability.

in riverine DOC flux are unclear. The Canadian Prairies undergo extreme seasonal climatic cycles; summers are warm and dry while winters are cold and punctuated by snowfall events. This region also accounts for over 80% of agriculture in Canada (Corkal et al., 2011). Climatic changes are projected to alter meteorological patterns associated with atmospheric changes, causing more prolonged summer heat and more frequent Polar jet stream divergences (Archer & Caldeira, 2008). These temperature excursions are associated with variability in precipitation, influencing drought, flood, and snow events, which, along with anthropogenic impacts such as water consumption, may become more common in coming decades (Cook et al., 2015;Schindler, 2001;Schindler & Donahue, 2006). Empirical patterns and global circulation models suggest that future climate scenarios for the western Canadian Prairies will cause water shortages due to the occurrence of drought conditions, and that increased rationing and storage capacity will be needed to enhance socio-economic resiliency (Corkal et al., 2011;Gan, 2000;Sauchyn et al., 2016;St. Jacques et al., 2014). It is unclear how human land use activities (i.e., land cultivation for agriculture and water management for irrigation) interact with long-term climatic patterns to shape biogeochemistry and DOC processing in heavily modified grassland river networks.
Cyclical climate oscillations importantly shape the hydrologic cycle of rivers (Enfield et al., 2001), and can both complicate the interpretation of temporal trends in shorter riverine monitoring data sets, and obscure human effects on nutrient and organic matter cycling (Noacco et al., 2019). Climate oscillations have been linked to long-term discharge data sets across a range of environments and these may intensify with climate change (Foley et al., 2002;Marcé et al., 2010;Regier et al., 2016;Toohey et al., 2016). In the Canadian Prairies, Pacific Decadal Oscillation (PDO) has a typical length of 10 or more years and has been shown to enhance streamflow during negative phases associated with colder, wetter winters. The El Nino Southern (ENSO) and Arctic Oscillations (AO) are less influential on streamflow patterns, but exert some influence on patterns of summer drought and cold air intrusion, respectively and cycle on shorter timespans compared to PDO (Gobena & Gan, 2009;Newton et al., 2014aNewton et al., , 2014b. ENSO had a lag effect of 2-3 years on regional meteorological conditions, with extended dry spells occurring in the second summer following an ENSO event (Bonsal & Lawford, 1999;Shabbar et al., 1997). The effects of the PDO, and to a lesser extent ENSO and AO, have been shown to influence the position of the jet stream and thus tracking of moist Pacific systems over southern Alberta that influence winter snowpack and subsequent spring flooding events (Foley et al., 2002;Newton et al., 2014aNewton et al., , 2014bRood et al., 2005Rood et al., , 2008. It is presently unclear whether these climatic controls extend to shape the long-term patterns of riverine DOC flux. Information linking these climate oscillation cycles to river biogeochemistry patterns over multi-decadal timescales is limited (Noacco et al., 2019), but select evidence suggests the connection may be important in predicting river DOC fluxes. For example, in the US Everglades, long-term monitoring showed that the wet phase of the Atlantic Multidecadal Oscillation (AMO) positively correlated with DOC fluxes, playing an important, secondary role in regulating DOC cycling (Regier et al., 2016). These findings highlight the utility of climate oscillations in explaining hydrology and therefore DOC flux variability through time, and similar relationships may extend to grassland rivers of the Canadian Prairies.
Here we use a unique data set (spanning 2-4 decades) of DOC and discharge data for the South Saskatchewan River Basin (SSRB), a Canadian grassland river network that has been extensively developed for dryland and irrigation agriculture, and includes multiple expanding urban areas (Kerr, 2017;Pentney & Ohrn, 2008). Our goal was to produce a generalizable framework of the drivers of DOC flux in grassland aquatic networks, by addressing two objectives: (a) Quantifying the flux of water and DOC moving through the SSRB across multiple spatial scales (sub-basin to entire watershed); and (b) Assessing how land cover, hydrologic regulation, and interannual climatic patterns shape trajectories of DOC flux.

The South Saskatchewan River Basin
The SSRB originates in the North American Rocky Mountains and covers approximately 146,000 km 2 of southern Alberta, Montana, and Saskatchewan ( Figure 1). The SSRB transitions from the mountains through the foothills to the prairies, accompanied by substantial shifts in land cover and development, and water usage and allocation. Southern Alberta is heavily agriculturalized and has high proportions of water withdrawal from its rivers (Schindler & Donahue, 2006). Three major tributaries, the Red Deer (RDR), Bow (BR), and Oldman (OR) Rivers, are sub-basins within the SSRB. Watersheds above each sample location (n = 12) were delineated using publicly available 3 arc s digital elevation models (HydroSHEDS; Lehner et al., 2008) in QGIS 3.12. Land cover for each watershed was extracted using 2015 land cover rasters at 30 m resolution (North American Land Change Monitoring System) and was reported as a percent of the total watershed area. Environment and Climate Change Canada and the Province of Alberta have measured discharge (Q; daily measurements) and water chemistry (monthly measurements) over multiple decades at upstream, midstream, and downstream sites in each of these sub-basins ( Figure 1). Discharge (HYDAT, https://wateroffice.ec.gc.ca/search/historical_e.html) and DOC data from the Governments of Alberta and Saskatchewan for each site extend back to 1998-1978.

Fluxes and Yields
We calculated fluxes and yields of DOC and water at the up-, mid-, and downstream sites for the three subbasins of the SSRB, the SSR at Medicine Hat at the confluence of the BR and OR, the confluence of the SSR and RDR (SSR at Leader site; SSR-L) in Saskatchewan, and the SSR in Saskatoon (SSR-Sask; Figure 1). We used an online LOADEST platform (Park et al., 2015;Runkel et al., 2004) to calculate fluxes and yields with 175-462 DOC observations (mean = 330 observations). LOADEST computes fluxes using centered Q and DOC data to eliminate collinearity and applying the adjusted maximum likelihood estimation method (Runkel et al., 2004). A significant gap (6 years) in the DOC data record for the RDR downstream site was filled by interpolation (we confirmed no difference in flux estimates with and without interpolation). Water and DOC fluxes were normalized to watershed area to compute watershed yields.

Statistical Analyses
Statistical analyses were performed using R and Sigmaplot 13.0 (SyStat Software). Outliers were calculated as values lying outside the 10th and 90th percentiles. Linear regression models predicting fluxes and climate oscillation indices were calculated using annual and monthly oscillation values with and without a 1-year lag effect, regression models with an r 2 above 0.05 were reported. Analyses of variance (ANOVAs) are reported as significant when p < 0.05. For temporal trend analyses, annual time series data were analyzed using Mann Kendal and Sen slope (SS) tests in R (R Core Team, 2019) with packages "Kendall" and "zyp."

River Discharge
Our data sets capture the extreme range of river conditions for this region including severe droughts (i.e., [2001][2002] and flooding events (1995 and 2005; Figure S1 in Suppporting Information S1), and thus reflect the full range of flux and discharge. Average annual discharge from these sites ranged from 1.2 to 6.8 km 3 y −1 (Figure 2a), similar to other arid river systems (Spencer et al., 2013), but less than other large Canadian river systems (Hudon  Tank et al., 2016). Annual discharge generally increased at the downstream sites corresponding to watershed area (p < 0.01), however, water yield decreased moving downstream (Figure 2b; p < 0.01), with elevation of the sampling site strongly predicting mean water yield (r 2 = 0.75, p < 0.001).

DOC Fluxes and Yields
Fluxes of DOC ranged from an average of 2.9 ± 0.7 Gg C yr −1 at the BR Upstream site to 38.6 ± 17.2 Gg C yr −1 at SSR-L ( Figure 2c). Consistent with annual discharge, interannual variability was high across these sites, with fluxes varying over an order of magnitude in all sites except upstream reaches. The ice-free season from April to late September accounted for most of the DOC flux ( Figure 3), with June alone accounting for 20%-33% of annual DOC flux. This dominant flux corresponds with peak DOC export in the spring freshet ( Figure S1 in Suppporting Information S1; Pederson et al., 2013). From May to July across all sites DOC flux remained high, on average accounting for >50% of annual DOC flux.
Average DOC yields were highest at upstream sites (Figure 2; p < 0.01) and ranged from 0.20 to 0.81 g m −2 yr −1 (Table S1 in Suppporting Information S1). Across all sites we observed interannual variability in DOC fluxes and yields with high (i.e., year 1995) and low (i.e., year 2000) outlier values (Figures 2c and 2d). Water and DOC yields decrease across sites transitioning from mountains with greater forest cover to grasslands with lower precipitation and higher agricultural land usage (Figure 4; Table S2 in Suppporting Information S1).

Interannual Trends in DOC Flux
Few trends in annual discharge were observed across these sites (Table S1 in Suppporting Information S1). The BR Up and RDR Down sites both had trends in DOC flux through time, however they were directionally different (Table S1 in Suppporting Information S1; BR Up: p = 0.04, tau = −0.27, SS = −0.002; RDR Down: p = 0.05, tau = 0.23, SS = 0.018). The long record of both DOC and discharge data allows us to assess the drivers of multi-year climate oscillations on discharge and DOC flux ( Figure 5). Across most sites the December (winter) PDO index from the previous year was significantly correlated with DOC flux (r 2 = 0.07 to 0.27, p < 0.05; Table  S3 in Suppporting Information S1), while other comparisons of non-lagged and multi-year lagged climate data showed no clear correlation with DOC fluxes (Table S3 in Suppporting Information S1). Correlations with other climate oscillations (monthly and annually) varied, with somewhat weaker correlations between Arctic Oscillation from July of the previous year and river discharge and DOC flux (r 2 = <0.05 to 0.20 and < 0.05 to 0.25, respectively; Table S3 in Suppporting Information S1). Previous year May ENSO was weakly correlated at some sites as well with river discharge exhibiting stronger correlations than DOC flux (r 2 = <0.05-0.24 and < 0.05-0.16, respectively; Table S3 in Suppporting Information S1).

Discussion
In stark contrast to other regions of the world (e.g., De Wit et al., 2016;Tank et al., 2016), we found that the flux of water and DOC through the SSRB is highly variable on both seasonal and annual timescales and shows no clear directional shift in recent decades (Figures 2 and 3). The high interannual hydro-climatic variability observed here may mask long-term directional changes in DOC export (Figures 2 and 4). Multi-year climatic oscillations were correlated to fluxes of water and DOC, especially wintertime PDO ( Figure 4) and summer AO and ENSO. By establishing the link between DOC flux and climate oscillations, relative to other potential long-term drivers, we provide insight into the drivers of aquatic DOC cycling in grassland rivers. This is especially important because most studies of DOC cycling are limited to short (1-2 years) timescales of observation, making it difficult to discern long-term directional changes in DOC dynamics from cyclical patterns (Noacco et al., 2019). Given the large variation between sites and through time, our results provide a well-defined estimate of DOC flux that will help to constrain continental C budgets (Butman et al., 2018), and trajectories of change in grassland river C cycling, or lack thereof for the SSRB for the given measurement period.

High Interannual Variability in DOC Flux Partly Explained by Hydro-Climatic Oscillations
Few long-term studies exist for grassland riverine DOC flux, thus our new analyses provide important biogeochemical context from an understudied but socio-economically important region of the world. Across all sites we observed large variations in DOC flux (Figure 2) that was most pronounced at the downstream sites where flux was high. For example, DOC flux at SSR-L ranged 67.4 Gg C yr −1 across the discharge record (Figure 2). These large variations in DOC flux indicate that interannual variability can be extreme ( Figure S1 in Suppporting Information S1) and that short periods of record could lead to under or over estimations of DOC fluxes and limits on how trends are interpreted. For instance, Ren et al. (2016) found that the Mississippi River displays high interannual variability related to years of record droughts and floods, and that environmental change throughout this large watershed area was responsible for increasing DOC fluxes related to land management practices and land use change. It may be that similar, interactive natural and anthropogenic controls are shaping the complex patterns and interannual variability observed for the SSRB. Continental scale climate oscillations that drive regional temperature and precipitation patterns appear to explain a portion of the interannual variability observed for DOC flux ( Figure 5). Our observation extends those from previous work in the SSRB that demonstrated a link between multi-decadal streamflow patterns and climate oscillations, especially PDO (Philipsen et al., 2018;Rood et al., 2016). The DOC flux was negatively correlated to previous year December PDO at all 12 study sites, with more negative PDO values corresponding to periods of higher DOC flux. This indicates that the well-defined increase in winter snowpack during more negative phases of the PDO (Gurrapu et al., 2016;Whitfield et al., 2010) likely enhances river flows in the subsequent spring flood pulse, and enhances DOC export. This mechanism may explain some of the interannual variability in DOC flux across these sites ( Figure 5). These results highlight the potential importance of large scale climate variability on DOC flux consistent with a handful of previous studies of DOC export and climate oscillations in other regions (Jones et al., 1996;Regier et al., 2016). These relationships have not, to our knowledge, been shown for a grassland river system over multi-decadal timescales. Establishing this connection both provides mechanistic insight into a key driver of interannual variability in DOC flux in river networks and provides new information to help predict future patterns in grassland DOC flux. Moving forward, however, if PDO becomes less predictable with climate warming (Li et al., 2020), the trajectories of change to DOC fluxes could become more difficult to predict.

Trends in DOC Flux Differ From Other Regions
Grasslands appear to be distinct from many boreal and Arctic catchments where increasing DOC fluxes reflect greater terrestrial inputs through time. In more northern regions, increasing DOC concentrations, flux, and water color have been observed in many rivers (de Wit et al., 2016;Erlandsson et al., 2008;Tank et al., 2016). We did not observe consistent and significant directional changes in DOC fluxes through time (Table S1 in Suppporting Information S1). The fraction wetland cover (Figure 4) is low compared to other regions, which could explain the low responsiveness of DOC to precipitation-driven changes in DOC export, since wetland-derived DOC is a major driver of catchment DOC loading in other regions (Spencer et al., 2013). This may imply that management concerns related to increasing DOC (i.e., aquatic light regime, dystrophy in drinking water reservoirs, DOC-fueled biological oxygen demand) may be less relevant here compared to other regions where DOC fluxes are clearly increasing. These observations are particularly important in a management context, for the maintenance of downstream ecosystems such as Diefenbaker Reservoir where there is concern about the frequency of anoxia (Hudson & Vandergucht, 2015), which is partially related to organic matter delivery and processing.

Reservoir Creation and Operation a Major Impact on Seasonal DOC Fluxes
It is clear that the installation of dams are modifying the global organic matter cycle worldwide through the enhancement of water retention and aquatic OM processing (Maavara et al., 2017). This is reflected in our study by both the divergent seasonal DOC flux patterns below the downstream Diefenbaker Reservoir relative to upstream sites, and the apparent decline in DOC flux below the reservoir. At all upstream sites, DOC fluxes were highest during the spring freshet period (Figure 3) since these watersheds receive a majority of their precipitation in the form of snow during the winter months (Toth et al., 2009). This is consistent with observations from similar snowmelt driven systems in mountain and high latitude environments (Finlay et al., 2006;Holmes et al., 2012;Raymond & Saiers, 2010). Conversely, we found an absence of clear seasonality at the most downstream site (SSR below Saskatoon), where DOC ( Figure 3) and water ( Figure S2 in Suppporting Information S1) flux consistently accounted for 20%-30% of the annual flux, lacking a prominent freshet peak. This may be attributed unique reservoir management that prioritizes augmented winter flow with lowered summer flow (Gober & Wheater, 2014). Interestingly, in the BR Basin similar flow regulation is not observed despite the proximity of the BR Midstream site to the City of Calgary, and the presence of flow control structures above Calgary (Sinnatamby et al., 2020). The decline in DOC flux that we observed ( Figure 2) could be attributed to hydrologic regulation that lengthens water residence time and drives the removal of ∼9.6 Gg DOC yr −1 by Diefenbaker Reservoir during the ice-free season (Finlay et al., 2010). Taken together, we conclude that Diefenbaker Reservoir acts to modify the natural DOC cycle of the SSR river network in multiple ways, making it an important control point (Bernhardt et al., 2017) for DOC cycling in the broader SSR network.

Land Cover Is a Strong Spatial Predictor of Water and DOC Yields
Across the sites and similar to other biomes, land cover was strongly correlated with yields of water and DOC (Figure 4). Percent forest and grassland cover were the strongest predictors of mean DOC yield across the watersheds. Forest cover was correlated (p < 0.01) with higher DOC yields characteristic of headwater, mountainous watersheds while grassland cover was negatively correlated with DOC yields (p < 0.01; Figure 4). Similarly, agricultural land cover was weakly, negatively correlated with mean DOC yield (p < 0.05), suggesting that forest removal during conversion of land for agricultural uses has the potential to lower the long-term yield of DOC in some watersheds, particularly in headwater river reaches, however, these regions are unlikely to experience this conversion (Wijesekara et al., 2012). The natural DOC yields from grassland dominated systems is likely to remain low despite conversion to agricultural land cover (Figure 4). This suggests that the headwaters are supplied with water (Pomeroy et al., 2009) and DOC that fuels the bulk of downstream microbial and biogeochemical DOC cycling.

Conclusions
Our modeling approach using multi-decadal time series demonstrates the interplay between land cover, hydrology, and climate teleconnections that collectively generate complex patterns of DOC flux through grassland watersheds. The flux of DOC and water through the SSRB may become more variable as climate oscillations change in response to increased global temperatures. These changes are unlikely to be confined solely to the SSRB, and comparable grassland ecosystems worldwide may see similar changes to riverine DOC export patterns through time. By taking a network-scale approach to quantifying DOC and water flux, we demonstrated that large grassland reservoirs can represent a critical control point on the landscape for DOC flux regulation, which may impact downstream C delivery and possibly drainage to coastal margins, including Hudson Bay downstream of the SSRB. As seen here, despite widely varying hydro-climatic conditions, mountain headwaters remain critical sources of water and carbon to sustain downstream grassland ecosystems.