Seasonality of pCO2 in a hard-water lake of the northern Great Plains: The legacy effects of climate and limnological conditions over 36 years

Biogeochemical processes are active year-round in ice-covered lakes, such that processes in one season can affect limnological conditions in subsequent seasons. However, the extent and nature of these legacy effects are poorly understood, particularly for the CO 2 content of lakes and when considering gas exchange with the atmosphere. Here, we used a unique 36-yr dataset of weekly limnological measurements of Buffalo Pound Lake in the northern Great Plains to assess seasonal changes in CO 2 concentration and ﬂ ux and determine how dependent lake pCO 2 is on limnological conditions of previous seasons. We found that the lake was a net source of CO 2 to the atmosphere (mean 18.5 (cid:1) 7.4 mol CO 2 m − 2 yr − 1 ), with spring potentially accounting for the majority (~ 64% (cid:1) 20%) of CO 2 ef ﬂ ux, assuming ice in spring was permeable to gas exchange (32.9% (cid:1) 19.8% if not). Analysis with generalized additive models (GAMs) demonstrated that current and antecedent seasonal conditions combined to explain 72.6% of deviance in spring pCO 2 , but that the strength of model predictions and the importance of antecedent conditions diminished in GAMs of summer (43.6%) and fall (23.3%) CO 2

It is now well established that inland waters contribute significantly to the global carbon budget (Cole et al. 2007;Prairie 2008;Tranvik et al. 2009), although many questions remain about the factors regulating variability in water-column pCO 2 at broad spatial and temporal scales. One such uncertainty relates to the legacy effects of antecedent water-column conditions on current ecosystem function. For example, biogeochemical cycling under ice can substantially alter the abundance and chemical form of macronutrients in spring (Kratz et al. 1987;Hampton et al. 2017) and, in the case of carbon (C), substantially increase CO 2 concentrations under ice (Kratz et al. 1987;Finlay et al. 2015). Additionally, although spring CO 2 flux has been shown to contribute significantly to total annual CO 2 flux in many lakes (Maberly 1996;Striegl et al. 2001;Ducharme-Riel et al. 2015), relatively few measurements of pCO 2 are available for shoulder seasons of summer, owing to logistical issues related to sampling during ice melt and formation. Given that lake pCO 2 is frequently elevated in spring and fall seasons relative to summer (Baehr and DeGrandpre 2002;Denfeld et al. 2015), it is important to better understand the magnitude and drivers of seasonal contributions to annual CO 2 fluxes to improve estimates of the role lakes in the global C cycle.
Seasonal variation in water-column pCO 2 in boreal lakes frequently follows predictable annual patterns of change in metabolic processes, particularly in ice-covered dimictic systems. In these lakes, CO 2 accumulates under ice in winter (Baehr and DeGrandpre 2002;Denfeld et al. 2015), causing a large efflux of CO 2 in spring when the ice melts and the water column circulates (overturn). pCO 2 levels are reduced in summer when the water column is stable and primary production increases, while pCO 2 often increases during fall when CO 2 from respired organic matter in the hypolimnion is mixed into the water column at fall overturn (Maberly 1996;Dillon and Molot 1997;Anderson et al. 1999;Baehr and DeGrandpre 2004;Ducharme-Riel et al. 2015). Deviations from this pattern can occur due to local variation in meteorological conditions (wind, atmospheric pressure, and storm runoff ), which affect lake stratification and gas solubility (Vachon and del Giorgio 2014) or which introduce labile allochthonous organic matter into the lake (López- Bellido et al. 2012).
Less is known about seasonal patterns of pCO 2 in hard water and saline lakes that account for nearly half of continental surface waters (Hammer 1986). In these hard-water systems, variation in pH, groundwater inputs, and calcite precipitation can uncouple lake pCO 2 from metabolically regulated processes (Striegl and Michmerhuizen 1998;Finlay et al. 2009;Stets et al. 2009). Moreover, the magnitude of atmospheric CO 2 exchange in spring and fall has not been widely quantified in the shallow polymictic lakes common in agricultural lowlands, but where summer CO 2 effluxes can be much less than that seen in dimictic lakes (Finlay et al. , 2015. In particular, hypolimnetic CO 2 accumulation should be relatively low in the absence of persistent thermal stratification, whereas frequent lake mixing should keep the vertical profiles of pCO 2 more uniform during the ice-free period (Anderson et al. 1999;Stets et al. 2009). Given the potential importance of such hard-water lakes in regulating atmospheric CO 2 exchange (Finlay et al. 2015) and the predominance of spring and fall CO 2 emissions in other lakes (Ducharme-Riel et al. 2015), further research is needed on the controls of seasonal and annual CO 2 content in polymictic hard-water lakes.
Under-ice processes can be influential for many biogeochemical cycles, including that of carbon (Kratz et al. 1987), and can thus result in legacy effects where antecedent conditions propagate into subsequent seasons (Meding and Jackson 2003;Hampton et al. 2017;Powers et al. 2017). For example, respiration rates can be high under ice (Denfeld et al. 2015), particularly near the sediments where the temperature tends to be warmer than in surface waters and where organic matter accumulation is high (Wetzel 2001). Photosynthesis can also be an important control of CO 2 immediately under ice when snow cover is limited DeGrandpre 2002, 2004;Pernica et al. 2017). As a result, the quantity of CO 2 accumulated under ice can be a function of the duration of ice cover (Finlay et al. 2015), the availability of nutrients and light for photosynthesis DeGrandpre 2002, 2004;Salmi and Salonen 2016;Pernica et al. 2017), and the quantity and quality of organic matter available for mineralization (Wetzel 2001;Hampton et al. 2017). Similarly, respiratory consumption of O 2 under ice is dependent on winter conditions, as well as previous seasons' primary production (Meding and Jackson 2003). However, while biogeochemical processes during winter conditions can affect limnological conditions in spring, it is less clear how these processes affect CO 2 flux at ice-off, or whether winter legacy effects continue through summer and fall.
Together, this evidence suggests that water-column pCO 2 at a given point in time is dependent on both present limnological conditions as well as those in preceding seasons. To evaluate this hypothesis, we used generalized additive models (GAMs) to quantify the magnitude and correlates of seasonal and annual CO 2 dynamics in a polymictic eutrophic hard-water lake that has been monitored year-round at weekly intervals for 36 yr. Our objectives were three-fold: (1) describe seasonal variation (spring, summer, and fall) in pCO 2 and potential atmospheric exchange in a polymictic lake; (2) quantify long-term (36 yr) trends in CO 2 dynamics and seasonality; and (3) evaluate the influence of antecedent environmental conditions (productivity and climate) on seasonal estimates of water-column pCO 2 . We predicted that spring pCO 2 would be influenced strongly by factors controlling the supply of labile organic matter and the duration of ice cover (Meding and Jackson 2003;Finlay et al. 2015) but that atmospheric CO 2 exchange in summer and fall would reflect the increasing influence of coeval meteorological and limnological conditions (Gerten and Adrian 2000;Winder and Schindler 2004). By integrating seasonal change with the importance of legacy effects, we hope to improve predictions of how future climate change may affect the contribution of lakes to the global carbon budget.

Study site
Buffalo Pound Lake is a natural lake that was impounded in 1939 and 1952 by the damming of the outflow into the Qu'Appelle River in southern Saskatchewan, Canada (Hall et al. 1999). The lake is long and narrow (1 km by 29 km), with an average depth of 3 m. The shallow depth of the lake, combined with long fetch along the prevailing storm track, results in a polymictic system that only rarely establishes weak thermal stratification (Dröscher et al. 2009).
Buffalo Pound provides the drinking water supply for the cities of Moose Jaw (population 45,000) and Regina (population 216,000), has water levels managed by the Saskatchewan Water Security Agency, and is maintained in part by hydrologic transfer from the upstream Lake Diefenbaker reservoir (Hall et al. 1999). The lake receives runoff from a 3310 km 2 agricultural catchment area in which nutrient-rich soils favor high nutrient influx and eutrophic conditions. Although other high pH lakes in this region typically ingas CO 2 (Finlay et al. 2015), Buffalo Pound is typically oversaturated with CO 2 and exhibits net outgassing of CO 2 to the atmosphere over the last 20 yr (Finlay et al. 2015).
As the main urban drinking water supply, Buffalo Pound has been monitored on a weekly basis since 1979 for 65 water quality parameters. Raw water taken from an inflow pipe 1 m above the bottom of the lake, at 3 m depth, is pumped into the water treatment plant for analyses, treatment for human use, and distribution. Measured parameters include physical (temperature), chemical (pH, nutrients, and major ions), and biological (chlorophyll a [Chl a], algae, and bacteria) properties. In this study, we used measured conductivity (μS cm −1 ), bicarbonate and carbonate (mg L −1 ), pH, and temperature ( C) to calculate pCO 2 and potential CO 2 flux. Furthermore, we explored proxies of planktonic metabolism (Chl a and dissolved organic carbon [DOC]) and physicochemical processes (water temperature and ice-cover duration) as predictors of variation in pCO 2 and CO 2 fluxes across the 36-yr time series. pCO 2 and CO 2 flux calculations CO 2 concentration (μM), pCO 2 (μatm), and CO 2 flux (mmol CO 2 m −2 d −1 ) were estimated for each sampling date from conductivity, water temperature, and pH measurements taken from the inflow water as described in Finlay et al. (2015). Dissolved inorganic carbon (DIC) concentrations were estimated using a previously derived relationship between measured DIC and conductivity for Buffalo Pound Lake (r 2 = 0.98, p < 0.001; Finlay et al. 2009). Given the elevated pH of the system (average pH during the open water period = 8.3), chemically enhanced C flux was calculated on each sampling date.
CO 2 flux was calculated as where [CO 2 ] lake is the concentration of CO 2 in the water, [CO 2 ] sat is the concentration of CO 2 at equilibrium with the atmosphere, α is the chemical enhancement of CO 2 flux at high pH (Hoover and Berkshire 1969), and k is piston velocity (cm h −1 ) as determined from Model B in Vachon and Prairie (2013), relating k to wind speed and lake surface area. Hourly wind speed as measured each day at 10-m height was collected from publicly available Environment Canada records for the city of Moose Jaw (Sta. 2967; http://climate.weather.gc.ca/). Flux was interpolated between time points by multiplying daily flux rates by 7 d to get a total potential flux for each of 52 weeks. The concentration of CO 2 at saturation with the atmosphere was taken as the global mean annual CO 2 concentration measured at Mauna Loa observatory. Raw water was collected between 07:00 and 07:30 h, once each week, and thus diel variations in pH were not considered in this analysis. Although pH can vary considerably during the day (Maberly 1996), an evaluation of continuous sonde data (15 min resolution) from Buffalo Pound during summer 2014 suggested no systematic bias in pCO 2 estimates due to the use of morning pH measurements. Specifically, pCO 2 at 07:00 h was correlated positively with mean daily pCO 2 (r = 0.44, p < 0.001) and was not consistently elevated or depressed relative to the daily values. Similarly, the pCO 2 values calculated from water treatment plant samples should be elevated relative to those of surface waters (Finlay et al. 2015), as water was extracted 1 m above the bottom (2 m from surface) and is more affected by sedimentary respiration and less by photosynthesis in this turbid system (mean summer Secchi depth = 1.1 m). However, given that the lake is polymictic year-round, we assumed that these surface-deep differences did not greatly influence our analysis of temporal variability and causal relationships, at least during the open water period (Finlay et al. 2015).
Ice-on and ice-off dates were provided by the Buffalo Pound water treatment plant and were determined as the day of year when the lake was 90% covered with ice (ice-on) or 90% ice free (for ice-off ).

Definition of seasons
As the goal of this study was to quantify the magnitude of variation in seasonal pCO 2 , identify potential controls thereof, and evaluate the importance of antecedent seasonal conditions on observed CO 2 , we needed to establish functional definitions to delineate seasons, which accounted for interannual variation in winter severity and ice-cover duration. Complete ice melt can take weeks and may vary with spring meteorology (Finlay et al. 2015), so we defined seasons based on a combination of potential gas exchange with the atmosphere (winter and spring) and consistent static calendar dates (summer and fall) that define when CO 2 concentrations were stable, as recommended by Anderson et al. (1999).
Herein, winter was defined as the period when the lake was completely covered with ice, and atmospheric gas exchange was negligible, beginning with the date of ice formation in the fall and continuing until the date of maximum modeled CO 2 concentration in Buffalo Pound (see below). The start of spring was defined as the date when CO 2 concentrations begin to decline, assuming that this pattern arises from loss of CO 2 to the atmosphere (see Discussion section), even if this occurred before the recorded ice-off date. Spring continued until the minimum CO 2 concentration was recorded (within 100 d of CO 2 maxima). Spring was further divided into two phases for CO 2 flux analyses: "potential spring," which was defined as the period of time between maximum CO 2 concentration and ice-off date, and "openwater spring," the period from documented ice-off until the minimum CO 2 concentration. Hereafter, "spring" refers to the period that includes both "potential" and "open-water" spring periods, and annual flux rates include "potential" spring, unless otherwise indicated. Summer was then calculated as the date after the spring CO 2 minimum continuing until 31 August of that year, whereas fall was defined as the period from 01 September until ice-on, as we have done previously (Finlay et al. 2015).

Statistical analyses
The time series of in situ pCO 2 was modeled using a GAM (Hastie & Tibshirani 1990;Wood et al 2016;Wood 2017;Simpson 2018) comprising terms to account for both withinand between-year variation in the time series. We chose to model the data using a GAM, because this approach better Finlay et al.
Regulation of seasonal pCO 2 accounts for nonlinearity of trends relative to other protocols (e.g., Mann Kendall test) and GAMs uniquely allow us to estimate the magnitudes of within-and between-year trends in the data, derive secondary estimates from the model (e.g., the magnitude of efflux at ice-out), and quantify uncertainties. For example, the commonly used (seasonal) Mann Kendall test does not estimate the magnitude(s) of trends, tests only for monotonic trends (which were not indicated in preliminary data screening), and does not allow derivation of secondary estimates as above.
Similarly, estimation of trends using parametric linear or generalized linear models would require us to a priori state the functional form of the within-and between-year trends in time series or perform model selection from among a set of complex polynomial models. Using GAMs, we avoid this subjective element of model specification by allowing the functional form of the trends to be determined from the data, whereas the use of splines avoids well-known bias issues at the ends of series that plague polynomial models. Details of candidate model selection and estimates are included in the Supporting Information.
To determine the start and end dates of spring, the bestfitting GAM for the pCO 2 time series was used to estimate annual mean differences between the minimum and maximum pCO 2 between days 50 (19 February) and 160 (~09 June). In this procedure, an estimate for the expected difference for each year can be derived by predicting daily pCO 2 for each day in the specified interval, finding the peak pCO 2 during the period, and calculating the difference between the two pCO 2 extremes. Of the 36 yr of data collection, 4 yr did not have a pronounced pCO 2 peak and those years were not included in subsequent statistical analyses. Uncertainty in the estimated pCO 2 trend was evaluated using 10,000 simulations of the trend from the posterior distribution of the fitted GAM (details in Supporting Information).
Variables known from the literature to affect pCO 2 content of prairie hard-water lakes were selected a priori to develop individual GAMs for each season to predict pCO 2 in Buffalo Pound (Meding and Jackson 2003;Finlay et al. 2009Finlay et al. , 2010. Specifically, spring CO 2 flux was expected to be dependent on icecover duration (longer ice cover resulting in greater accumulation of respired CO 2 ), and the productivity of the previous summer (providing the material for respiration over winter) approximated as Chl a (Meding and Jackson 2003;Finlay et al. 2015). In contrast, we expected that summer and fall pCO 2 would be more heavily dependent on pCO 2 in the previous season and coeval limnological conditions (Finlay et al. , 2010. To test these hypotheses, we developed GAMs to evaluate the effects of coeval mean water temperature, Chl a, DOC, and ice-cover duration in the models and examined legacy effects on seasonal lake pCO 2 by including mean values from preceding seasons. Given that not all variables were measured in all years, direct comparisons of Akaike information criteria (AIC) were not appropriate for determining the bestfitting model for each season. Instead, we selected models that maximized deviance explained, adjusted R 2 (R 2 adj ), and sample size (n), in addition to a qualitative exploration of the model fits. More details of model selection are included in the Supporting Information.
GAMs were estimated using the mgcv package (version 1.8-22; Wood 2017), and graphics were plotted with package ggplot2 (Wickham 2009) for R (version 3.5.1; R Core Team 2018). R scripts and data to reproduce the analysis can be found at https://doi.org/10.5281/zenodo.2431144. "potential spring," the period between maximum pCO 2 and the ice-off date (black bars), and "open-water spring," the period between ice-off and minimum pCO 2 (dark gray bars). Summer and fall are represented by medium and light gray bars, respectively. Flux is calculated using summed calculated CO 2 flux based on weekly data, extrapolated to 7 d. Winter flux under ice was considered to be zero.

Results
Estimation of pCO 2 from water chemistry suggested that Buffalo Pound Lake should outgas CO 2 during the open-water season of most years (Fig. 1). Estimates of total annual CO 2 flux (including potential spring flux) ranged from a minimum of 4.36 mol m −2 yr −1 in 1988 to a maximum of 41.97 mol m −2 yr −1 in 1992, with a mean (AE SD) annual flux rate of 18.53 AE 7.38 mol m −2 yr −1 . Instantaneous CO 2 fluxes ranged dramatically from an efflux of 886.8 mmol m −2 d −1 to an influx of 49.1 mmol m −2 d −1 . Buffalo Pound Lake also exhibited ingassing of CO 2 in summer or fall seasons of four years (1979, 1987, 1991, and 2012), but none of these events resulted in the basin experiencing a net influx of CO 2 when calculated at an annual scale. Over the entire 36-yr dataset, total spring CO 2 efflux averaged 63.8% of total flux (AE 19.8%), but this value declines to 32.9% (AE 19.8%) when only open-water spring flux is considered. In contrast, CO 2 efflux was lowest in summer (14.0% of annual total) and increased slightly in fall (22.2%).
There were no pronounced decadal-scale trends in the estimated CO 2 content or effluxes from Buffalo Pound Lake (Fig. 2). Seasonal averages of pCO 2 (μatm) varied by year, with spring pCO 2 being the highest, averaging 1797 μatm (range 480.1-3334 μatm), summer pCO 2 is the lowest and with less variability (average 683, range 62.3-1387 μatm), and fall intermediate between the spring and summer (average 818.9, range 327-1981 μatm). Winter pCO 2 averaged 1730 μatm with a minimum of 499 and a maximum of 3687. We did not see any significant (p > 0.1) decadal-scale trends with pCO 2 in each season vs. year (regression of data in Fig. 2) or with annual pCO 2 averages vs. year.
The best-fitting GAM to model the pCO 2 time series was a tensor-product smooth of sampling date and day of year as a seasonal trend, which varied smoothly with the between-year trend. This model was also the most complex in terms of the effective degrees of freedom (EDF = 506.8), but provided better fit to observed data (AIC = 25,769) than did the next best model (AIC = 25,926). The best-fit model explained 97% of deviance in the pCO 2 data, with an adjusted R 2 of 0.90. Qualitatively, the best-fit model also better explained large annual peaks in pCO 2 , as well as yearto-year variation in the magnitude of that peak. Overall, the mean pCO 2 estimated by the best-fit model oscillated slowly over the 36 yr of study and did not exhibit sudden changes between years.
In most years, modeled pCO 2 increased under ice, declined substantially at spring ice melt, stayed low during summer, and, with a few exceptions, remained low until ice formation in the fall (Fig. 3). Ice-cover duration varied > 7 weeks across the 36-yr period, from a minimum of 133 d in 2000 to a maximum of 183 d in 1979 (mean 156.7 AE 12.5 d). The length of spring CO 2 decline varied from 3 weeks (in 1986 and 1988) to 15 weeks (1991) with a mean (AE SD) of 10.5 (AE 2.9) weeks. On average, the spring pCO 2 decline started 4.9 AE 2.3 weeks before the observed date of ice melt. In contrast, there were few indications of a sudden change in estimated pCO 2 in fall such as would be expected if CO 2 was released suddenly from hypolimnetic waters at fall mixis (Vachon and del Giorgio 2014;Ducharme-Riel et al. 2015). Instead, pCO 2 peaked during winter in most years (Fig. 3b), with only a few years showing limited CO 2 build-up under ice (1984, 1987, 1995, and 2012). There were no statistically significant trends in the relationship between fall pCO 2 and the timing of ice formation.
Comparison of GAMs developed independently for spring, summer, and fall seasons revealed that the influence of antecedent seasonal conditions on pCO 2 declined from spring to fall ( Fig. 4; Supporting Information). For example, the best-fit GAM for spring pCO 2 used the scaled t distribution for heavily tailed data and explained 72.6% of the deviance (R 2 adj = 0.64, n = 31). This model demonstrated that mean values increased with the previous summer's average Chl a concentrations and the duration of ice cover but declined with spring water temperature (Fig. 4A-C) and was comparable when only open-water spring pCO 2 was considered (77.5% deviance explained, data not shown). In contrast, GAMs were significant but less predictive for both mean summer and fall pCO 2 values. In these models, mean summer pCO 2 declined with increases in mean summer Chl a concentrations but tended to increase with pCO 2 recorded in the previous spring, particularly at high values (Fig. 4D,

Discussion
Weekly estimates of CO 2 content of Buffalo Pound Lake over 36 yr demonstrated that a eutrophic hard-water lake could remain a net source of CO 2 to the atmosphere (Fig. 1) despite elevated algal production and pH (Finlay et al. , 2015. However, we found clear evidence of legacy effects in all seasons Fall Temperature (°C) Fall pCO 2 (µatm) Fig. 4. GAM results for seasonally averaged pCO 2 in Buffalo Pound. Spring average pCO 2 was best explained with a combination of (a) Chl a concentration in the previous summer, (b) ice-cover duration, and (c) current water temperature (GAM deviance explained = 72.6%), while summer average pCO 2 was explained using (d) current summer Chl a concentration and (e) average pCO 2 in the preceding spring (deviance explained = 43.6%), and fall was best explained using (f ) current fall Chl a concentration, (g) preceding summer pCO 2 , and (h) current fall water temperature (deviance explained = 23.3%). Plots are partial plots of the smooth terms in the model, and the y axis is the intercept plus the partial effect of the individual smooths. (Fig. 4). Specifically, GAMs suggested that spring CO 2 content and efflux rates were derived mainly from winter metabolism of organic matter that was produced the previous summer, with longer ice-cover duration serving to increase both factors (Kratz et al. 1987;Baehr and DeGrandpre 2004). While CO 2 effluxes were lower during summer and fall than in spring (Fig. 3), pCO 2 values in summer and fall were also regulated by interactions between current lake production and the legacy of the previous season's CO 2 concentration. Overall, we found no evidence of a major release of CO 2 during fall, as often occurs in thermally stratified lakes (Hesslein et al. 1991;Cole et al. 1994;Ducharme-Riel et al. 2015) and conclude that annual CO 2 budgets are strongly influenced by spring efflux, and therefore antecedent limnological conditions (Fig. 4). Given the highly synchronous patterns in CO 2 and other limnological parameters of lakes in this region (Vogt et al. 2011;Finlay et al. 2015), we believe that these results will be representative of other lakes in the Northern Great Plains (e.g., Meding and Jackson 2003;Donald et al. 2015;Maheaux et al. 2016).

Temporal variation in CO 2
Total CO 2 flux in spring strongly influenced the magnitude of annual CO 2 flux in Buffalo Pound, owing to the elevated respiratory-derived CO 2 under ice. Total spring efflux of CO 2 averaged 9.69 mol CO 2 m −2 , which is comparable to that observed in DOC-rich Wisconsin and Finnish lakes (1.1--13.7 mmol CO 2 m −2 spring −1 ; Striegl et al. 2001), whereas vernal open-water release alone was 32.9% AE 19.8% of annual values. Summer pCO 2 remained relatively low and, with the exception of a few anomalous years, there was little indication of hypolimnetic CO 2 release in fall (Vachon and del Giorgio 2014;Ducharme-Riel et al. 2015). Mean fall values for pCO 2 and CO 2 flux were only slightly higher than in summer, consistent with the polymictic status of Buffalo Pound and the irregular occurrence of thermal stratification during summer (Dröscher et al. 2009). In general, Buffalo Pound Lake was a net annual source of CO 2 to the atmosphere, with mean total annual flux of CO 2 (18.53 mol m −2 yr −1 ) and mean daily flux rates (103.95 mmol m −2 d −1 ), comparable to some other hard-water systems (Striegl and Michmerhuizen 1998) but higher than boreal lakes (Rantakari and Kortelainen 2005;Abnizova et al. 2012;Ducharme-Riel et al. 2015).
Atmospheric exchange of CO 2 during spring accounted for 63.8% of annual CO 2 efflux from Buffalo Pound Lake, a value much higher than that found elsewhere (Anderson et al. 1999;Ducharme-Riel et al. 2015) and reported earlier for this site (Finlay et al. 2015). However, this result includes both the "potential" and "open-water" spring fluxes and assumes both that ice is highly permeable to gas exchange prior to its complete disappearance and that CO 2 efflux to the atmosphere is the main mechanism reducing CO 2 content during spring. If instead, gas exchange is limited through even fully fractured ice due to limited hydrologic and atmospheric exchange  and degassing occurs following formation of marginal (lateral) open water immediately prior to full ice melt (Loose and Schlosser 2011), then CO 2 efflux declines to~33% of annual flux, a value more in line with boreal systems. Furthermore, we note that CO 2 could have declined under ice in Buffalo Pound due to alternate mechanisms, including elevated primary production by attached and motile algae (Salmi and Salonen 2016;Hampton et al. 2017), redistribution of CO 2 -rich deep waters by convective water-column currents (Kelley 1997;Mironov et al. 2002;Pernica et al. 2017) or chemical dissolution of sedimentary CaCO 3 (Finlay et al. 2015). Taken together, these observations suggest that our spring estimates represent the maximum possible CO 2 efflux and illustrate that intensive studies of under-ice processes in the weeks prior to ice melt are needed to fully characterize the magnitude and importance of vernal CO 2 release.
Here, our GAM analyses reaffirmed that warmer winters with reduced ice-cover duration results in lower winter CO 2 accumulation and thus emissions in spring, as has been observed previously (Finlay et al. , 2015. We did not, however, see similar consistent decadal trends for annual or seasonal CO 2 flux. In part, these differences arise because the duration of ice cover declined during the most recent 20-yr interval studied by Finlay et al. (2015), but not during the entire 36-yr period included here (r 2 = −0.01, p = 0.5). Given the elevated contribution of spring CO 2 fluxes to the annual budgets, as seen here and in other lake districts (Kratz et al. 1987;Cole et al. 1994;Striegl and Michmerhuizen 1998), climate change is likely to profoundly alter future lake CO 2 fluxes. Specifically, future prairie climates will be warmer and drier (Sauchyn and Kulshreshtha 2008;Lapp et al. 2009;Newton et al. 2014), with less ice cover (Shuter et al. 2013), patterns that should reduce the magnitude of CO 2 emissions in spring. This reduced spring CO 2 flux will translate into a reduction of long-term annual flux provided there is an alternative loss pathway for this C, such as carbonate precipitation or organic matter sedimentation (Tranvik et al. 2009), both of which are common in productive, hard-water lakes. Although reduced ice cover can potentially affect annual primary production, we found no significant relationship between ice duration and summer Chl a (r 2 = 0.061, p = 0.15) during the past 36 yr when ice cover varied from 133 to 183 d. Instead, given that regional phytoplankton biomass is a complex function of temperature and nutrient influx (Vogt et al. 2018) and given the legacy effects seen herein, we infer that further changes in spring and annual CO 2 emissions will also depend heavily on the effectiveness of nutrient management strategies (Leavitt et al. 2006;Bunting et al. 2016).

Controls of seasonal pCO 2
Controls of CO 2 content and potential atmospheric exchange were strongly influenced by lake production but differed in form and function among seasons. During spring, pCO 2 was strongly and positively influenced by mean Chl a content of the previous summer, consistent with microbial respiration of autochthonous organic matter consuming oxygen (Meding and Jackson 2003;Powers et al. 2017) and producing CO 2 (Kratz et al. 1987;Finlay et al. 2015) under ice. While pCO 2 levels in spring were also enhanced by the duration of ice cover (longer time for CO 2 accumulation) and cool spring temperatures (enhanced gas solubility), the paramount effect of Chl a may reflect the highly eutrophic conditions in Buffalo Pound (Hall et al. 1999;McGowan et al. 2005). Although present at elevated concentrations, allochthonous DOC in Buffalo Pound is recalcitrant relative to other sources (Williamson et al. 1999;Guillemette et al. 2017), largely unrelated to rates of microbial production (Finlay et al. 2010), and was not included in our final models of CO 2 content. Instead, it appears that factors regulating mid-summer production may be unanticipated but important controls of subsequent spring CO 2 efflux.
GAMs also suggested that coeval Chl a was correlated negatively with pCO 2 during in both summer and fall, consistent with a strong role of phytoplankton uptake during photosynthesis as seen in other autotrophic lakes (del Giorgio and Peters 1994). Additional influences of microbial processes are indicated also by the presence of a positive relationship between temperature and pCO 2 in fall, a pattern consistent with the role of bacterial respiration of OM (del Giorgio and Peters 1994; Cole et al. 1994), rather than changes in gas solubility as the lake cools (Pinho et al. 2016). Again, we were unable to detect an effect of DOC on pCO 2 from Buffalo Pound in either summer or fall models, possibly because groundwater inputs, carbonate buffering, calcification, and anaerobic metabolism also decouple the relationship between allochthonous DOC influx and microbial metabolism in regional lakes (Bogard and del Giorgio 2016;Stets et al. 2017). Full carbon budgets in each season, including catchment loading of organic and inorganic carbon, would be required in order to fully evaluate these alternatives.

Legacy effects of climate and limnological conditions
Analysis of a 36-yr continuous time series of water chemistry demonstrated that instantaneous estimates of CO 2 content in lake waters are regulated by current limnological conditions as well as the persistent influence of lake conditions in earlier seasons. Such legacy effects are well known from studies of terrestrial biogeochemistry (Cuddington 2011) and land-water linkages (Martin et al. 2011) but are less well understood for in situ biogeochemical cycles (Meding and Jackson 2003;Hampton 2015). Recent studies suggest that spring water chemistry is strongly influenced by under-ice processes (Powers et al. 2017), and, consistent with that view, our GAM analysis showed that vernal pCO 2 could be predicted best (72.6% of deviance explained) from a combination of mean Chl a concentration during the previous summer, ice-cover duration during the antecedent winter, and coeval spring water temperature. Overall, historical parameters had a paramount effect on spring model performance (Fig. 4), while coeval biological parameters were nonsignificant (e.g., spring Chl a), suggesting that vernal CO 2 fluxes were controlled mainly by limnological and climatic conditions in earlier seasons. Given the potential importance of spring CO 2 efflux to the annual CO 2 budgets (see above), these findings suggest that atmospheric warming (Finlay et al. 2015) and surface water eutrophication (Leavitt et al. 2006;Bunting et al. 2016) will interact in complex manners to regulate the importance of hard-water lakes in global carbon budgets.
Comparison of GAMs developed for individual seasons demonstrates that the strength of legacy effects declines continuously through the ice-free period, both in terms of explanatory power and the influence of historical parameters on gas fluxes (Fig. 4). Specifically, the predictive power of GAMs declined bỹ 50% in successive seasons, from spring (72.6%) to summer (43.6%) and fall (23.3%), while only pCO 2 levels during the preceding season were retained as a secondary predictor in summer and fall models. Although speculative, we infer that these declines may reflect the progressive accumulation of effects of intervening meteorological events (e.g., wind and low pressure cells), which are known to influence CO 2 fluxes at the scale of days to weeks (Morales-Pineda et al. 2014). In addition, as coeval Chl a concentrations were retained in latter models, yet are strongly influenced by summer water temperatures and nutrient content in Buffalo Pound and other regional lakes (Vogt et al. 2018), we infer that immediate controls of these limnological parameters may override the importance of legacy effects on CO 2 content and flux.

Conclusion
These analyses used a 36-yr time series with weekly resolution from a hard-water lake to demonstrate that instantaneous pCO 2 is regulated by a combination of current limnological conditions and legacy effects from earlier seasons. This legacy effect was most pronounced in spring and declined throughout the ice-free season. The form and identity of physicochemical controls also changed through time, with climate (ice-cover duration) being the strongest predictor in spring model, and coeval estimates of lake metabolism (Chl a) mainly regulating pCO 2 in summer and fall. The strength of these relationships reflects our ability to predict CO 2 in the future. Specifically, the strong spring relationship suggests that future climate warming and reduction of ice cover will diminish the importance of antecedent conditions, and may reduce annual CO 2 emissions to the atmosphere, particularly if efforts to reduce nutrient loading in this region are successful. This relationship explained 72% of the variability in the data, which allows for predictions of how future climate change will affect lake carbon processing in this system.