At the laboratory, water samples were stored at 4°C and filtered within 24 h of collection. Sufficient water to visibly discolour a Whatman 1.2-µm GF/C filter was filtered was recorded, and filtrate was frozen (−10°C) in sealed darkened film canisters until analysis for chlorophyll pigments.

Greenhouse gas samples were analyzed at the Global Institute for Water Security, University of Saskatchewan. Headspace gas samples were analyzed for the dry molar fraction of CO₂, CH₄, and N₂O using a fully calibrated Scion 456 Gas Chromatograph (Bruker Ltd.) with Combipal autosampler (CTC Analytics–PAL System), using argon as the carrier gas. A flame ionization detector was used for methane samples (<100,000 ppmv). A thermal conductivity detector was used for CO₂ and high level CH₄ concentrations (>100,000 ppmv). N₂O was measured using a micro-electron capture detector and argon/methane (90/10) as a makeup gas (injector temperature 60°C, column temperature 60°C, detector temperature 350°C). All gases were calibrated using mixed gas standards (Praxair) with the addition of a single gas N₂O standard (0.1 ppmv; Scotty).

Water chemistry and water isotope samples were analyzed at the Institute for Environmental Change and Society at the University of Regina. Water chemistry samples were analyzed for soluble reactive phosphorus (SRP), total dissolved nitrogen (TDN), nitrate and nitrite (NO_x), ammonia (NH₃), dissolved organic carbon (DOC), and dissolved inorganic carbon (DIC). DOC and DIC were analyzed using standard analytic procedures on a Shimadzu model 5000 A total carbon analyzer (Finlay et al., 2009; Webb et al., 2019b). Dissolved nitrogen (TDN, NO_x, NH₃) and dissolved phosphorus (TDP, SRP) were measured using standard procedures on a Lachat QuikChem 8,500 (APHA-AWWA/WEF 1998; Bergbusch et al., 2021). Chl-a concentration was standard trichromatic spectrophotometric methods (Jeffrey and Humphrey 1975; Finlay et al., 2009). Water was analyzed for δ¹⁸O-H₂O using a Picarro L2120-I cavity ring-down spectrometer (CRDS), and δ¹⁸O inflow was calculated using coupled isotope tracer methods as an estimate of water source (winter vs. summer precipitation, groundwater, Yi et al., 2008; Webb et al., 2019a; Haig et al., 2020).

The squared Brunt-Väisälä buoyancy frequency (s⁻²) was calculated as a measure of maximum stratification strength within the water column. This approach uses the greatest density gradient from the vertical water temperature profiles taken at intervals of 0.5 m (0.25 m in shallower waterbodies) using the rLakeAnalyzer package (Read et al., 2011) in R (version 4.0.5; R Core Team, 2021) (Webb et al., 2019b).

The dry molar fractions of CO₂, CH₄, and N₂O were corrected for dilution and converted to concentrations according to the solubility coefficients and atmospheric partial pressure of each gas and the salinity of the water. Four replicate gas concentrations for each site were averaged for CO₂, CH₄, and N₂O to estimate mean concentrations for each sampling period. We removed an extreme CH₄ outlier for in one reservoir seen during April (93.5 µM) and assumed an ebullitive event that was captured in the associated water sample.

The concentration of each of the gases were used to estimate the diffusive flux for CO₂, CH₄, and N₂O. The diffusive flux (f _C) was calculated using the gas transfer velocity (k _c), gas concentration of the water (C_water), and the ambient air concentration (C_air) using the following equation (Eq. 1): f C = k c ( C w a t e r − C a i r ) (1)

The ambient air concentration of each gas was taken as the average over the sampling period at the Mauna Loa NOAA station from April to September 2018, inclusive. The average concentrations were 408.94 µatm (or ppm) for CO₂, 1.86 µatm (or ppm) for CH₄, and 0.33 µatm (or ppm) for N₂O, and corrected for the barometric pressure for each sampling effort. The diffusive flux was calculated using the gas transfer velocity as measured on similar agricultural reservoirs in 2017 (Webb et al. 2019a; Webb et al. 2019b). We did not estimate ebullitive CH₄ fluxes, although we recognize this will lead to an underestimation of the true total CH₄ efflux.

Integrating seasonal GHG flux across all gases is challenging, owing to temporal variability and the assumptions required to extrapolate between sampling dates. Here, we calculated total CO₂-eq flux using several methods to bound true values. First, the “July only” procedure used only the daily CO₂-eq value observed in July multiplied over the total number of ice-free days to extrapolate over the entire open water period (e.g., metanalysis in DelSontro et al., 2018). Second, the “Seasonal” protocol used the CO₂-eq values measured in April, May, July, and September for the respective intervals summed over the entire open water period (e.g., Striegl and Michmerhuizen 1998; Finlay et al., 2019). Third, we calculated a “Predicted” flux using the generalized additive models (GAMs) of the smoothed seasonal splines to interpolate fluxes on each date between sampling periods. For this last method, we took both the mean smooth value and the upper 95% credible interval (“Upper CI”) to estimate mean and maximum seasonal CO₂-eq flux. Given the lack of consistent diel trend (see Section 3.4 below), we used measured flux rates from daytime determinations only. In all cases, the CO₂-eq flux was calculated by first converting the moles of diffusive flux value to grams of diffusive flux, and then multiplying by the sustained global warming potential (45 and 270 for CH₄ and N₂O, respectively) if the diffusive flux value was positive, and by the sustained global cooling potential (203 and 349 for CH₄ and N₂O, respectively) if the diffusive flux value was negative (Neubauer and Megonigal 2015).

Generalized additive models (Wood 2006; Wood et al., 2016) were used to assess the seasonal trends in CO₂, CH₄, and N₂O concentration in the reservoirs, as they are particularly useful for estimating linear, nonlinear, and nonmonotonic relationships between response variables and predictors (Webb et al., 2019a; b; Swarbrick et al., 2019). All modelling was done using the mgcv package (Wood 2011; Wood 2016) for R (version 4.0.5; R Core Team, 2021). GAMs were also used to link CO₂, CH₄, and N₂O concentrations to various biotic and abiotic predictive parameters to identify which factors might regulate changes in GHG concentrations in agricultural reservoirs. Covariates included those thought to control seasonal variation in GHG levels during the open-water period and were based on spatial analysis of these agricultural reservoirs (Webb et al., 2019a; Webb et al., 2019b). Variables included nutrient content (DOC, DIN, SRP), production (Chl-a, DO), and stratification strength (buoyancy frequency). The CO₂ model additionally included water sources (as δ¹⁸O of inflow), specifically groundwater, rainfall, or snowmelt (Yi et al., 2008; Haig et al., 2020). Conductivity was additionally included in the CH₄ model as a proxy for ionic content and, in this region, sulfate concentrations (r² _{cond- SO4} = 0.64, p < 0.01). Similarly, pH was included in the N₂O model because of its known effects on N₂O production and consumption processes (Quick et al., 2019). Site identity (ID) was included as a random effect in the model, to account for the differences between individual sites.

A Gamma distribution was used with a log-link function, for positive, continuous responses of CO₂ and CH₄. We used a gammals distribution for N₂O concentrations to model the log mean and log scale parameter (or standard deviation/variance) for the probability density function. The gammals distribution was used because two sites exhibited high, outlying N₂O concentrations compared to the rest, so this allowed for grouping (high and low) of sites separately from the rest of the sites to inform a non-normal distribution. Sites were put into the high group if 25% or more of their N₂O concentrations were outliers. Basis size, dispersion of residuals, homogeneity of variance, and the relationship between the observed and predicted response were assessed in each model to ensure assumptions were not violated. Residual marginal likelihood (REML) was used for selection of smoothness parameters (Wood, 2011). To help with model selection, the double penalty approach of Marra and Wood (2011) was used. An additional penalty is applied to the perfectly smooth parts of the basis (the functions in the penalty null space) of each smooth function in the model, which allows entire smooths to be effectively penalised out of the model, while accounting for the selection procedure in the statistical tests applied to model terms. Parameters predicting CO₂, CH₄, and N₂O concentrations were considered significant at 95% confidence level (α = 0.05) for each waterbody.

Physical and environmental parameters of the 20 reservoirs exhibited a wide range of values among sites and over the season (Table 1). The mean monthly max depth (mean = 1.88 m, standard deviation (SD) = 0.92), Secchi depth (0.61 ± 0.43 m), DO saturation (89.5 ± 31.2%), and Chl-a (44.8 ± 61.1 μg L⁻¹) of the reservoirs did not change significantly between any month (ANOVA, p-value > 0.05). Mean surface water temperature was low in April (mean = 7.98°C, SD = 4.38), reached a maximum in July (mean = 22.7°C, SD = 1.47), and declined again in September (mean = 6.99°C, SD = 0.87, ANOVA post-hoc Tukey test, p-value < 0.001). The pH of the surface water was lowest in April (mean = 8.29, SD = 0.84), reached a maximum in July (mean = 8.94, SD = 0.61), and decreased again in September (mean = 8.62, SD = 0.44, ANOVA post-hoc Tukey test, p-value = 0.014). Conductivity increased from April (mean = 737 μS cm⁻¹, SD = 591) to September (mean = 1824 μS cm⁻¹, SD = 1722, ANOVA post-hoc Tukey test, p-value = 0.029). Salinity also increased from April (mean = 0.53 ppt, SD = 0.42) to September (mean = 1.25 ppt, SD = 1.11, ANOVA post-hoc Tukey test, p-value = 0.033).

Water chemistry parameters were highly variable both between months and among reservoirs (Table 2). Soluble reactive phosphorus (seasonal mean = 320 μg P L⁻¹, SD = 515), TDN (2,793 ± 2,141 μg N L⁻¹) and NH₃ (seasonal mean = 316 μg N L⁻¹, SD = 654) concentrations did not have a significant change between monthly means (ANOVA post-hoc Tukey test, p-value = 0.001). DOC concentrations also increased from April (mean = 20.1 mg C L⁻¹, SD = 6.54) to September (mean = 31.3 mg C L⁻¹, SD = 12.0, ANOVA post-hoc Tukey test, p-value = 0.004). The concentration of NO_x was highest initially following ice-off in April (mean = 644 μg N L⁻¹, SD = 1,054, ANOVA post-hoc Tukey test, p-value = 0.009).

Temporal variation in concentrations of CO₂, CH₄, and N₂O (Figure 1A) was identified as significant seasonal trends when analyzed by GAMs for all gases (Figure 1B). The concentrations of each of the gases was highly variable between reservoirs, yet there were common trends for each gas. For example, concentrations of CO₂ were greatest immediately following ice-off in April (mean = 45.8 µM, SD = 39.3), declined during spring and summer (May–mean = 23.5 µM, SD 23.6, July–mean = 16.5 µM, SD 17.0), then increased again in fall (mean = 32.0 µM, SD 31.2, GAM, p-value = 0.013, deviance explained = 49%). In contrast, CH₄ concentration were initially low in the spring (April–mean = 0.849 µM, SD = 1.54; May–mean = 0.73 µM, SD = 1.28), but increased into a mid-summer peak (mean = 1.32 µM, SD = 1.70) before decreasing in fall (mean = 0.40 µM, SD = 0.47, GAM, p-value = 0.0005, deviance explained = 91.4%). Finally, concentration of N₂O exhibited a strong pulse during ice-off in April (mean = 35.2 nM, SD = 50.3), followed by a sharp decrease to low values in spring and summer (May–mean = 10.1 nM, SD = 1.5; July–mean = 15.4 nM, SD = 36.1), before rising slightly again prior to ice formation (mean = 13.3 nM, SD = 4.31, GAM, p-value < 0.0001, deviance explained = 89.5%).

Estimates of diffusive fluxes of CO₂, CH₄, N₂O were all highly variable among sites (Figure 2). Despite among-site variation, diffusive flux of CO₂ varied little between months (ANOVA, p-value = 0.161), with the overall mean (= SD) of 19.7 ± 56.6 mmol m⁻² day⁻¹. Similar patterns were recorded for the diffusive flux of CH₄, with no seasonal difference between monthly means (ANOVA, p-value = 0.654) and an overall flux of 2.9 ± 10.9 mmol m⁻² day⁻¹. Diffusive flux of N₂O was variable among sites but not months (ANOVA, p-value = 0.361), with an overall mean diffusive N₂O flux of 9.7 ± 52.9 μmol m⁻² day⁻¹.

Most water quality parameters did not change substantially during the diel cycle (Figure 3). In general, temperature was greater during the day and lower at night, although there were up to 4°C differences among sites (Figure 3D). Similarly, while DO exhibited some high frequency variation, concentrations were generally lower at night than in the day (Figure 3E). In contrast, few trends were observed for conductivity, and pH exhibited diel cycle (elevated in day) at only one site with elevated Chl-a values and high diel changes in DO. Nutrients (TN, TP, SRP, DOC, DIC) and Chl-a did not exhibit any observable changes during the diel cycle in any sites, while NO_x concentration and ammonia concentration only exhibited notable variation at one site with high DIN levels (Supplementary Table S1).

Concentrations of CO₂, CH₄, and N₂O showed no consistent trends during a 24-h period in the four reservoirs sampled (Figure 3). At three sites, CO₂ concentrations were consistently undersaturated over the diel cycle, while values were supersaturated at the fourth site yet no more variable through time. CH₄ concentrations at two sites remained low throughout the diel cycle, while the other two basins showed inconsistent but elevated variation during the 24 h study period. Overall, N₂O concentrations were consistently low and undersaturated during the diel cycle.

Seasonal changes in CO₂ concentration in the reservoirs were explained best by a GAM using the pH of surface water as the sole predictor (Supplementary Figure S2). This model explained 89.1% of deviance in CO₂ and was characterized by a strong negative relationship with CO₂ above pH of 8. When pH was not included in the model, variation in CO₂ was explained best by a combination of DO saturation (Figure 4A), DIN concentration (Figure 4B), DOC levels (Figure 4C), buoyancy frequency (Figure 4F), and water source (as δ¹⁸O; Figure 4G). Dissolved oxygen, buoyancy frequency, and δ¹⁸O exhibited significant negative relationships with variation in CO₂ concentration, while DIN and DOC had positive relationships. Overall, this model explained 63.9% of deviance in the CO₂ levels.

Seasonal changes in CH₄ concentrations were explained best by a combination of DO saturation (Figure 5A), concentrations of DIN (Figure 5B), DOC (Figure 5C), and SRP (Figure 5D), as well as surface water temperature (Figure 5F), and conductivity (Figure 5G). Dissolved oxygen, SRP, and conductivity exhibited significant negative relationships with the seasonal CH₄ concentration, while DIN, DOC, and surface water temperature exhibited significant positive relationships with methane levels. Model deviance explained was 74.5% of variation in CH₄ levels.

Regulation of seasonal variation in N₂O concentrations was also complex, with significant effects of changes in DO saturation (Figure 6A), DIN concentration (Figure 6B), SRP concentration (Figure 6D), Chl-a concentration (Figure 6E), buoyancy frequency (Figure 6F), and pH of the surface water (Figure 6G). Variation in N₂O levels were related positively to changes in DO, DIN, and SRP, and negatively to variation in Chl-a, buoyancy frequency, and pH, while the model explained 87.4% of the deviance in N₂O levels.

Estimates of total CO₂-eq flux over the open-water period varied dependent on the method used to estimate in situ GHG concentrations (Figure 7), but all were within an order of magnitude of each other. Mean seasonal CO₂-eq flux was comparable between the “July-only” (9,811 ± 10,269 mmol m⁻², Figure 7A) and “seasonal” methods (10,446 ± 17,893 mmol m⁻², Figure 7B), and were roughly twice the fluxes predicted using GAMs (5,027 ± 4,806 mmol m⁻², Figure 7C). All methods suggest that methane was the dominant GHG responsible for CO₂-eq flux, and contributed more than 50% to total GHG efflux in 16 of 20 sites (Figure 7, yellow bars). While the mean seasonal values were comparable among methods, there were some discrepancies within individual sites. For example, Site 23 A demonstrated extremely high CO₂-eq flux, attributable to one very high CH₄ measurement in April. This likely reflected an ebullition event at this time that was not captured when the “July-only” method was used, and which was constrained using the GAM-predicted method.

CO₂ concentration in agricultural reservoirs was high in the spring, low throughout the summer, and began to increase again in fall (Figure 1B). This seasonal pattern was consistent with that observed in boreal and other hard-water systems, where CO₂ accumulates under the ice in the winter, causing a pulse of CO₂ released in the spring following ice off (Baehr and DeGrandpre 2002; Finlay et al., 2015; Denfeld et al., 2016). Depressed CO₂ in summer is attributed to isolation of surface waters leading to elevated photosynthetic rates that increase pH and deplete dissolved pCO₂ to values below atmospheric values (Huotari et al., 2009; Holgerson 2015). Conversely, elevated fall CO₂ concentrations appear to reflect loss of water-column stratification, and release of hypolimnetic CO₂ which had been accumulated during the summer period of stratification (Huotari et al., 2009; Holgerson 2015; Webb et al., 2019b; Finlay et al., 2019). Agricultural reservoirs in this study have a high relative depth (5.6 ± 3.1%), and thus exhibit a high resistance to vertical mixing (Wetzel and Likens 1991; Holgerson et al., 2022) during summer. As elevated primary production sinks into the hypolimnion and decomposes during summer, CO₂ is released leading to supersaturation of this gas in deep water (Huotari et al., 2009; Finlay et al., 2019).

Seasonally, pronounced negative correlations between CO₂ concentration and DO saturation (Figure 4A) and positive correlations between CO₂ and both DIN and DOC (Figure 4C) are indicative of metabolic controls of CO₂ concentrations (Holgerson 2015; Balmer and Downing 2017; Vachon et al., 2019). The metabolic stoichiometry between CO₂ and DO is well documented (Vachon et al., 2019) and represent autotrophic and heterotrophic processes occurring both in the water column or emanating from the sediments (Roulet et al., 1997; Jonsson et al., 2003; Rantakari and Kortelainen 2005; Kortelainen et al., 2006; Holgerson 2015). Here, the significant positive relationship between CO₂ concentration and DIN levels (Figure 4B), as well as DOC levels (Figure 4C), and lack of significant relationship with Chl-a (Figure 4E) suggests that elevated nutrient content stimulated microbial respiration rather than primary production, such as seen in other artificial waterbodies and prairie lakes (Ollivier et al., 2018; Wiik et al., 2018; Webb et al., 2019b; Peacock et al., 2019), and natural temporary wetland ponds (Holgerson 2015). As shown in mesocosm fertilization experiments, moderate eutrophication of shallow prairie waters first elevates primary production and DO saturation, then leads to water-column anoxia due to build-up of dissolved and particulate organic matter (Donald et al., 2011; Bogard et al., 2020). Consistent with this mechanism, monitored reservoirs exhibited low ratios of dissolved N:P characteristic of N-limited regional lakes (Donald et al., 2011) (Table 2), while on-going whole-reservoir fertilization experiments show that N fertilization both increases surface water Chl-a levels and causes deep water supersaturation of CO₂ (C.A.C. Gushulak and P.R. Leavitt, unpublished data). At present, it is unclear whether the relationships observed in this or other studies reflect water column processes or the disproportionate effect of warm sediments in these shallow water bodies.

The negative correlation between CO₂ concentration and buoyancy frequency, a measure of stratification strength (Figure 4F), is consistent with the effects of seasonal thermal stratification in isolating deeper waters and constraining CO₂ loss to the atmosphere (Huotari et al., 2009; Holgerson 2015). Here, CO₂ concentrations are elevated both in spring after ice melting allows release of the winters’ accumulation of respired CO₂ (Baehr and DeGrandpre 2002; Finlay et al., 2015; Denfeld et al., 2016), as well as in fall after when cooling air temperatures erode thermal stratification and release hypolimnetic CO₂ (Huotari et al., 2009; Holgerson 2015). Between these intervals, agricultural reservoirs can exhibit thermal stratification that isolates respired CO₂ derived from sediments and prevents accumulation in surface waters. In general, this pattern contrasts findings from a spatial survey of prairie farm reservoirs in 2017 that found a significant positive relationship between CO₂ concentration and buoyancy frequency (Webb et al., 2019b). We speculate that this difference arose because sampling for Webb et al. (2019b) took place during the summer and would have mainly captured variation in the degree of thermal stratification among sites, rather than the sharp seasonal patterns of water column mixing exhibited in this study.

CO₂ concentration and δ¹⁸O exhibited a pronounced negative relationship (Figure 4G) suggesting that elevated inflow of snowmelt or groundwater resulted in higher CO₂ concentrations in these farm reservoirs, consistent with observations in the 2017 spatial study (Webb et al., 2019b). In general, regional farmers site the reservoirs in existing lowland or wetland areas, often with direct contact to shallow aquifers, or with abundant overland runoff. In this region, snowmelt in spring is both the main source of overland flow (Pomeroy et al., 2007), and the mechanism recharging shallow aquifers (Haig et al., 2020) and, in both cases, exhibited highly depleted δ¹⁸O values relative to rainwater sources (Haig et al., 2020; Haig et al., 2021). Further, this groundwater has high concentrations of both DOC and DIC, and is frequently supersaturated with CO₂ (Macpherson 2009; Webb et al., 2019b). Seasonal monitoring of piezometers for water isotopes, CO₂ levels and dissolved C concentrations will be required to further resolve this mechanism.

We did not observe consistent diel fluctuations of CO₂ concentrations in the four reservoirs sampled (Figure 3). Previous studies have shown that CO₂ concentrations are typically lower during the daytime when primary production exceeds respiration (P > R), but increase overnight due to continued respiration (Liu et al., 2016; Raymond et al., 2013; Wiik et al., 2019). The magnitude in the shifts of CO₂ concentration is often linked to the water-column productivity, with the largest diel changes observed in systems with the highest productivity (Hanson et al., 2003; Morales-Pineda et al., 2014; Shao et al., 2015; Wiik et al., 2018). Here Chl-a concentrations ranged from 2.9 to 296 μg L⁻¹, therefore, a strong diel change in CO₂ concentration would be expected in the more productive sites. Lack of pronounced variation through the day may reflect the high pH of these systems which would favour rapid conversion of CO₂ to bicarbonate (HCO₃ ⁻) and carbonate (CO₃ ^2-). Furthermore, the high total alkalinity (241.3 ± 159.4 mg L⁻¹) of the systems likely resisted diel variation in pH in all but one reservoir (Figure 3D), as diel changes in CO₂ were much less than the total reservoir content of DIC (Dodds 2002; Stets et al., 2017).

CH₄ concentrations showed a marked peak during July-August (Figure 1B) similar to those seen high latitude wetlands (Pickett-Heaps et al., 2011), aquaculture ponds (Yang et al., 2015), shallow urban basins (van Bergen et al., 2019), and prairie wetlands (Bansal et al., 2016). In general, this interval corresponds to the period of greatest water temperatures and is consistent with peak CH₄ levels occurring when water temperatures are greatest in this study (Figure 5F) and the observation that metabolism of methanogenic bacteria is greater in warmer waters (Segers 1998). Such a positive relationship of CH₄ and temperature has been observed in other shallow systems including prairie wetland ponds (Bansal et al., 2016), urban ponds (van Bergen et al., 2019), natural temporary wetlands (Holgerson 2015), and sites with abundant Phragmites growth (Kim et al., 1998). Methane ebullition was not included in our study; however, we acknowledge that this process can contribute significantly to GHG fluxes from inland water and may also vary seasonally.

Dissolved CH₄ is expected to be more abundant under anoxic conditions (Glass and Orphan 2012), therefore, the significant negative relationship with DO saturation was expected (Figure 5A). Such a trend is common in diverse ecosystems, including boreal lakes (Kankaala et al., 2013), natural temporary wetland ponds (Holgerson 2015), and our earlier spatial survey of prairie farm reservoirs (Webb et al., 2019b).

The strong negative relationship between CH₄ concentration and conductivity (Figure 5G) was also observed in the 2017 spatial study, which is likely attributed to sulfate limiting methanogenesis (Webb et al., 2019b). Sulfate-reducing bacteria have been shown to outcompete methanogenic bacteria at sulfate concentrations as low as 60 µM in freshwater sediments in oligotrophic lakes (Lovley and Klug 1983), and as low as 200 µM in eutrophic lakes (Winfrey and Zeikus 1977; Lovley et al., 1982). A linear regression showed that conductivity and sulfate are highly related in these agricultural reservoirs (r ² = 0.64, p-value < 0.01). Based on this calculation, the sulfate concentration ranges from 0.33 to 35.9 mM in these agricultural reservoirs. Surface water conductivity significantly increased over the open water season (Table 1) and is likely one of the mechanisms driving down the CH₄ concentrations in fall.

Nutrient (DOC, DIN, and SRP) levels were also important predictors of CH₄ concentrations in agricultural reservoirs (Figures 5B,C,E), although solute concentrations did not appear to be correlated to primary production as Chl-a (Supplementary Figure S2.1). As noted above, effects of allochthonous nutrients may be felt mainly through changes in the redox potential and dissolved oxygen content of these shallow ecosystems, rather than via augmentation of primary production. In particular, elevated production (Figure 5E) and sedimentation of labile organic matter may have favoured decomposition and anoxia in sediments, thereby improving habitat and substrate for methanogenesis (Smith et al., 1999; Anderson et al., 2002; Ollivier et al., 2019). The pronounced negative relationship between CH₄ and SRP may suggest that there is competition for substrate between methanogenesis, and denitrification at high P levels, consistent with the increase in N₂O concentration is seen at high SRP concentrations (Figure 6D). Such a mechanism has also been observed in anaerobic wetland sediment slurry incubations (Kim et al., 2015).

CH₄ concentrations did not show a consistent diel trend in the four reservoirs sampled, despite obvious variation in both temperature and DO in the surface water (Figure 3). In contrast, distinct diel patterns of CH₄ production have been seen in other systems where gas levels increase following sunrise, peak in the early afternoon, then decline quickly a stable low level at night (Neue et al., 1997; Kim et al., 1998; Xing et al., 2004). In general, mid-day maxima were seen here only in ponds with greatly elevated CH₄ levels, possibly suggesting that the absence of strong diurnal maxima in some sites reflects a generally low abundance of methanogenic bacteria in some reservoirs. While speculative, we suggest that differences in the influx of groundwater- or sub-surface sulfate (SO₄ ⁼) may explain the variation in absolute and relative CH₄ levels during a 24 h cycle (see below), where the presence of sulfur-reducing bacteria can inhibit methanogenesis (Lovley and Klug 1983; Pennock et al., 2010).

Although freshwaters are often significant emitters of N₂O (Beaulieu et al., 2008; Beaulieu et al., 2011; DelSontro et al., 2018), patterns and controls of the seasonality of N₂O concentrations are poorly understood. The seasonal trend observed here shows that following an initial pulse in N₂O production immediately following ice melt in the spring, these reservoirs were often undersaturated in N₂O for much of the ice-free period (Figure 1B). These low N₂O values may be attributed to complete denitrification of dissolved NO₃ ⁻ to di-nitrogen gas (N₂), similar to patterns seem in some lakes (Piña-Ochoa and Álvarez-Cobelas 2006), including Lake Kasumiguara, a shallow eutrophic lake that does not have an anaerobic zone (Hashimoto et al., 1993). Widespread undersaturation of N₂O was also observed in regional farm reservoirs during 2017 (Webb et al., 2019a); however, this study provides the first evidence that N₂O undersaturation last much of the summer, a finding which has not previously been recorded for nitrogen-rich agricultural landscapes. Better understanding of the mechanisms controlling this process may be important to regulating net CO₂-eq fluxes from small water bodies in the continental interiors that represent much of the global agricultural regions.

The significant relationships seen between DO saturation, DIN (NO_x + NH₃), Chl-a, and buoyancy frequency (Figure 6) were similar to those seen spatial surveys of agricultural reservoirs (Webb et al., 2019a). This finding suggests that mechanisms relating to productivity and mixing are controlling N₂O concentrations in reservoirs both seasonally and spatially (Webb et al., 2019a). The pulse of N₂O in the spring may be explained by the build-up of N₂O under ice over the winter season, which is a common observation in prairie and boreal lakes (Soued et al., 2016; Cavaliere and Baulch, 2018). The lack of atmospheric exchange (during ice cover) combined with ongoing N₂O producing processes such as nitrification and denitrification (Cavaliere and Baulch, 2018), without the competing factor of primary productivity for DIN (NO_x was elevated in April, Table 2), may explain the higher N₂O levels. The switch to N₂O undersaturation coincides with a decline in NO_x concentration owing to elevated algal uptake, along with the onset of temporary water column stratification and a decline in the DO saturation in the deep waters (Table 1) where complete denitrification can proceed. Complete denitrification leads to the consumption of N₂O in the process of producing N₂ (Quick et al., 2019), which suggests that the conditions in agricultural reservoirs favour complete denitrification during most of the open-water period. Additionally, the significant positive relationship between N₂O concentrations and SRP may be driven by direct, or indirect effects, via impacts of SRP on algal biomass and carbon production, or via direct effects on N₂O producing microbes.

N₂O concentration in the four reservoirs sampled on a diel timescale remained under-saturated throughout the 24 h cycle. This pattern is similar to that seen in an arid zone alpine pond which remained undersaturated in N₂O for most of the diel cycle (Molina et al., 2021). Previous studies have diel variation in N₂O concentrations, primarily in rivers and streams, although comparisons among these site exhibit little common diel pattern among habitats. For example, some studies show clear diel patterns with greater concentrations observed during nighttime relative to those during the daytime (Rosamond et al., 2011; Wu et al., 2018a; Wu et al., 2018b), while other locations show a peak during the day (Molina et al., 2021) as well as a peak during the nighttime (Rosamond et al., 2011). One study found a sharp decrease in N₂O concentration during the nighttime (Harrison et al., 2005) and others did not find a consistent diel trend (Baulch et al., 2012).

We calculated areal CO₂-eq flux for the farm ponds to range from an average 32.0 mmol m⁻² day⁻¹ using the GAM-predicted flux, to 66.5 mmol m⁻² day⁻¹ using the seasonal method, over the 157-days open-water season. These values are somewhat lower than the mean (129 mmol m⁻² day⁻¹) but well within the range (−10 to 1,462 mmol m⁻² day⁻¹) of the daily flux estimates in the 2017 spatial survey (Webb et al., 2019a; Webb et al., 2019b). Similar to findings from previous studies (DelSontro et al., 2018; Ollivier et al., 2019), we observed that CH₄ was the major driver of CO₂-eq flux (average 91.5% across all methods), while N₂O contributed the least (0.73%). The contribution of CH₄ to total CO₂-eq flux is likely to be even higher when ebullition is also considered (DelSontro et al., 2016; Rosentreter et al., 2021).

Findings from this study show that all three GHGs exhibited significant seasonal trends, thereby suggesting that extrapolation from a single date to the entire open-water period may produce inaccurate estimates of the total CO₂ equivalent (CO₂-eq) flux from these systems. In our study, the seasonal trends of each GHG (Figure 1), combined with the fact that CH₄ is typically the dominant gas contributing to CO₂-eq flux (Figure 7) would suggest that extrapolating measurements in mid-summer are likely to overestimate seasonal flux. The GAM-predicted seasonal fluxes were, on average, half of that of the other two methods, likely because the modeled fluxes constrained occasional extreme effluxes GHG and avoided excessive extrapolation through time series. For example, the exceptionally high estimate of seasonal flux at site 23 A in Figure 7B reflects a pulse of CH₄ in April (observed CH₄ = 90.7 µM) which was not included in the July-only flux calculations, and was constrained in the GAM models. Depending on whether this extreme point was interpolated to the next sampling date or not, estimated CO₂-eq flux differed during the spring period by up to 72,800 mmol m⁻² depending on whether elevated values were (∼75,600 mmol m⁻²) or were not (2,800 mmol m⁻²) used to estimate vernal fluxes. Although ebullition can occur due to changes in production, water depth, or meteorological conditions (DelSontro et al., 2016), observed elevated values are unlikely to persist over the 3-week spring interval. Conversely, GAM predicted values missed the likely-real pulsed efflux, and thereby likely underestimated seasonal fluxes. Overall, the “seasonal” calculation also likely overestimates CO₂-eq flux by inflating short-term emissions, while the July-only calculation will miss seasonal trends, and the use of the GAM smooths may dampen any short-term pulses of GHG efflux and lose the contribution of these to the seasonal flux rates. Reconciling these discrepancies to provide the best estimate of total open-water GHG fluxes will require continuous monitoring of GHG with sondes, now possible for CO₂, but still limited for CH₄ and N₂O due to more limited sensor sensitivity.

This study shows that agricultural reservoirs exhibit significant seasonal trends in the concentrations of all major GHGs, and that the concentrations of these gases are regulated by a combination of physical, chemical, and biological mechanisms. Additionally, it was found that inclusion of estimates of this seasonal variability may be essential to accurately total CO₂-eq fluxes to the atmosphere. However, reconciling accurate net GHG over the annual scale will only truly be overcome through increased GHG monitoring via in situ sensors to capture extreme changes in water conditions (such as at ice off). The Northern Great Plains, as the largest agricultural area of Canada, is an area where agricultural water management is ongoing. Reservoirs in the landscape mediate substantive fluxes, although more work is required to reduce error in estimation of annual efflux from these often-eutrophic reservoirs. Perhaps more importantly, the Northern Great Plains are a region of extremely high variability in factors known to influence GHG efflux, and upscaling efforts aimed at understanding the role of inland waters in GHG efflux need to account for spatial variability in solutes, including sulfate.