Underway pCO2 Surveys Unravel CO2 Invasion of Lake Superior From Seasonal Variability

This study observed seasonal trends and inferred drivers of CO2 biogeochemistry at the air‐water interface of Lake Superior. Underway carbon dioxide partial pressure (pCO2) was measured in surface water during 69 transects spanning ice free seasons of 2019‐2022. These data comprise the first multiannual pCO2 time series in the Laurentian Great Lakes. Surface water pCO2 was closely tied to increasing atmospheric pCO2 over a 100‐day CO2 equilibration timescale, while seasonal variability was controlled equally by thermal and non‐thermal drivers during the ice‐free season. Comparison to previous modeling efforts indicated that Lake Superior surface pCO2 increased with two decades of rising atmospheric CO2. Spatial heterogeneity in CO2 dynamics was highlighted by a conductivity‐based delineation of “riverine” and “pelagic” regimes, each of which was associated with net CO2 influx over Julian days 100–300 on the order of 25 Gmol C. These findings refine previous estimates of Lake Superior C fluxes, support predictions of anthropogenic CO2 invasion during the ice‐free season, point to new observation strategies for large lakes, and highlight an urgent need for studies of changes to lacustrine C cycling.


Introduction
Measurements of carbon cycling in the Earth's hydrosphere are central to understanding global biogeochemical cycling and responses to perturbation (Le Quéré et al., 2013).Continuing anthropogenic emissions of carbon dioxide (CO 2 ) are increasing atmospheric concentrations at an unprecedented rate, which may force changes in carbonate equilibria in the oceans (Feely et al., 2001), in soils (Oh & Richter, 2004), in rivers (Raymond & Hamilton, 2018), and in lakes (Alin & Johnson, 2007).
Many studies of CO 2 dynamics in inland waters collect and analyze discrete water samples for inorganic C parameters including pH, dissolved inorganic carbon (DIC), total alkalinity (A T ), and partial pressure of carbon dioxide (pCO 2 ) (Cole et al., 1994).Direct measurements of CO 2 flux across the air-water interface are also made via floating chamber or eddy covariance methods (Podgrajsek et al., 2014).Calculation of one inorganic C parameter from two others remains fraught with uncertainty due to ongoing challenges associated with measurement and equilibrium calculations in low-ionic strength waters (Liu et al., 2020;Minor & Brinkley, 2022;Young et al., 2022).Furthermore, constructing time series of discrete water chemistry measurements is time-and labor-intensive and may not resolve inorganic C cycling in water bodies with high spatial or temporal variability (Schilder et al., 2013).To bridge these gaps in observational capabilities, instruments measuring inorganic C parameters continuously or autonomously have been developed and deployed in aquatic systems spanning the lacustrine-marine spectrum (Bushinsky et al., 2019;Lynch et al., 2010).Recent years have seen applications of pH and pCO 2 underway sensors that perform with uncertainties similar to those of discrete sample measurements (Ma et al., 2019;Takeshita et al., 2018).
Inorganic C chemistry remains less-studied in inland waters than in marine systems (Phillips et al., 2015), due in part to high heterogeneity within and among water bodies, each of which may require divergent approaches to inorganic C cycling analysis depending on hydrology, climate, or terrestrial influence.The approaches used in marine systems (where variations in inorganic C parameters are often smaller and require more precise measurements) might be applied to large lakes and provide new insights into aquatic systems.In this way, large lakes may serve as stepping-stones for application and further development of oceanographic chemical techniques in inland waters.Their great volume and relatively small terrestrial influences lend them a more constant chemistry and physics than their smaller peers.The largest of lakes share with oceans similar biogeochemical features and relative importance to local and global biogeochemical cycling (Sterner et al., 2017).On the other hand, large lakes respond more rapidly than the global ocean to perturbation.Their hydrologic residence times (c.190 years for Lake Superior) are shorter than that of the global ocean (millennia).Holomictic lakes experience full water column mixing at least annually, which represents a homogenizing driver not observed in oceans.For these reasons, large lakes can act as test systems for investigations of environmental variables, with responses occurring on more accessible spatial and temporal scales for research than the global ocean (Sterner, 2021).
The Laurentian Great Lakes lie on the border of the United States of America and Canada and within the historical and contemporary lands of Native American and First Nations.They constitute the largest contiguous aquatic ecosystem on Earth (Wetzel, 2001), yet C cycling in the Great Lakes is not well-understood (Minor & Oyler, 2021).It remains unclear to what extent the Great Lakes are net sources or sinks of CO 2 to the atmosphere (McDonald et al., 2013;Urban & Desai, 2009).Alin and Johnson (2007) concluded that they are annual net CO 2 sources, but Cotner et al. (2004) and Urban (2005) highlighted apparently unbalanced C budgets in Lake Superior, the best-studied of the five lakes in terms of carbon cycling.Bennington et al. (2012) and Atilla et al. (2011) noted that studies of CO 2 cycling in Lake Superior have been biased by sparse observations restricted to the ice-free period (generally April-November ±2 months), and could not close the cycle by modeling all C inputs and outputs.These pioneering studies were confounded by observations of inorganic C cycling that were sparse, irregular or unrepresentative of the lakes as a whole.This situation is similar to that of the Southern Ocean or South Pacific Ocean, in which limited observation hindered attempts to constrain biogeochemical budgets (Takahashi et al., 2009).Large lakes functioning as "sentinels, integrators, and regulators of climate change" (Williamson et al., 2009) exert significant influence on regional and global C budgets (Cole, 2013) and demand more detailed study.This research focuses on surface water pCO 2 measurements to illustrate the C cycle of Lake Superior in unprecedented spatial and temporal detail.pCO 2 in water responds to physical (temperature, pressure, salinity), chemical (pH, DIC, A T , CaCO 3 dissolution/precipitation), and biological (production, respiration) drivers (Zeebe & Wolf-Gladrow, 2001), such that a comprehensive understanding of pCO 2 variability sheds light on a suite of biogeochemical functions.
As a direct driver of CO 2 flux across the air/water interface, pCO 2 in surface waters acts as an important factor of atmospheric CO 2 accumulation.The surface ocean has maintained a lower mean pCO 2 than the atmosphere in recent decades, resulting in about a quarter of annual anthropogenic CO 2 emissions being absorbed by the ocean each year (Friedlingstein et al., 2022).In contrast, many inland waters are supersaturated with CO 2 recycled from the terrestrial C sink and emitted to the atmosphere or exported to the coastal ocean (Cole et al., 1994;Raymond et al., 2013).These marine and lacustrine C cycles are linked by C imports and exports along the terrestrialaquatic continuum which form important elements of comprehensive global C budgets, especially in light of their perturbation by anthropogenic global change (Tian et al., 2023).Accurate predictions of climate change and mitigation efforts require an improved understanding of surface waters as sources and sinks of CO 2 and other greenhouse gases (Cavallaro et al., 2018).Lake Superior is characterized by a small catchment-to-surface area ratio of 1.55 (Urban, 2005) and is underlain by a weathering-resistant igneous mineralogy leading to exceptionally dilute, soft, and carbonate-poor water chemistry (Weiler, 1978).Its primary production is dominated by a large proportion of picoplankton and Nano plankton, with identifiable phytoplankton species including an abundance of diatoms and small flagellated plankton (Fahnenstiel et al., 1986;Ivanikova et al., 2007;Kovalenko et al., 2019), and has been classified between oligo-and mesotrophic, with relatively high transparency and low phosphorus concentration (Sterner, 2010).Multiple aspects of the Lake Superior's limnology are changing: its waters are warming faster than the overlying atmosphere (Austin & Colman, 2008), the concentration of most of its major ions is increasing (Chapra et al., 2012), and its nutrient stoichiometry is shifting as the result of long-term nitrate accumulation (Sterner, 2011).Interannual trends in A T , pH, and pCO 2 have proven difficult to constrain due to covariation with lake level, influence from Dreissenid calcification in tributaries, large measurement uncertainty, and spatial heterogeneity (Minor & Brinkley, 2022).These poorly-understood trends, and especially the potential for acidification (Phillips et al., 2015), contribute to the need for a sustained campaign of spatially-and temporallycomprehensive measurements of the inorganic carbon system in Lake Superior.
In this work, underway pCO 2 measurements gathered by instrumentation aboard RV Blue Heron from cruises of opportunity during four consecutive field seasons (April-November 2019-2022) provided a survey of unprecedented spatial and temporal scope describing inorganic C cycling drivers and variability in Earth's largest lake by surface area.This information was used to illustrate ice-free season surface pCO 2 variability, infer trends in pCO 2 and CO 2 flux over space and time and establish the interplay of thermal and non-thermal drivers of pCO 2 .The results demonstrate a pathway toward comprehensive CO 2 budgets for the Laurentian Great Lakes via novel observation strategies and improved modeling efforts.

Methods
Underway instrument data sets from 69 transects of RV Blue Heron were compiled.These efforts included singleday endeavors near the vessel's home port of Duluth Minnesota, as well as multi-week transects across the Laurentian Great Lakes (Figure 1).Water was directed from the ship's water intake line at 2 m depth through a suite of sensors measuring parameters including dry molar fraction of carbon dioxide (xCO 2 ), sea surface temperature (SST), and sea surface conductivity.These were combined with wind velocity, barometric pressure, and air temperature collected from an onboard meteorological station.The multi-year span considered in this study permitted evaluation of interannual variability in inorganic C biogeochemistry despite limited cruises in 2020 and 2021 due to challenges associated with the Coronavirus pandemic.
xCO 2 was measured in water from the underway system at 2 s intervals by a Sunburst Sensors Super CO 2 instrument equipped with a dual showerhead equilibrator which recirculates air at a rate of 1 L min 1 through a twostage system as described by (Pierrot et al., 2009), which has been found to have a response time around 2 min (Evans et al., 2019).The sample water flow rate c. 4 L min 1 through two c. 2 L showerhead equilibrators ensured complete equilibration with air-sea xCO 2 gradients never exceeding the manufacturer's guideline of 1,000 ppm.This instrument has proven robust to wide ranges of temperature and sample pH in intercomparison experiments in which it was used as the reference instrument (Shangguan et al., 2022) and its equilibrator design aligns with established best practices (Pierrot et al., 2009).
Measurements from four standard gases with CO 2 concentrations between 0 and 1,018 ppm were performed every 2 hr (Figure S1 in Supporting Information S1) and the 60 s before and after calibration removed from the time series to prevent memory effects.The slope and intercept values from a type-I linear regression of measured versus standard xCO 2 were used to correct surface water xCO 2 before conversion to pCO 2 (Equation 1) A nearly-identical instrument demonstrated pCO 2 measurement uncertainty of ±5 μatm (DeGrandpre et al., 2020).SST and conductivity were obtained from a SBE21 thermosalinograph every 2 s.Conductivity was converted to practical salinity using the equations of Hill et al. (1986).Previous work found such an extrapolation from seawater salinity to Lake Superior's ionic strength fit overconstrained carbonate system equilibrium calculations better than assuming freshwater carbonate dissociation constants (Minor & Brinkley, 2022).Salinity as used in this work cannot be quantitatively compared to ocean salinity due to differenes in major ion ratios which invalidate the seawater equation of state.Salinity is used here only for equilibrium calculations.
Wind velocity was measured with a Young 05106 wind monitor on a mast 10 m from the sea surface.Air temperature was obtained from a Young 41372VC thermometer.It was assumed that mast-measured windspeed (corrected for vessel velocity) approximated neutral wind speed at 10 m (U 10n ) sufficiently well for the determination of instantaneous CO 2 flux.Measured pCO 2 and calculated CO 2 flux were averaged for each day of each transect in 0.01°× 0.01°boxes (approximately 1.1 × 0.8 km at the latitude of Lake Superior) to normalize distributions of pCO 2 and CO 2 flux on an areal basis, prevent biases due to vessel idling, and provide observations on a scale greater than the equilibration time of the instrument.pCO 2 was calculated as a product of ambient atmospheric pressure (p atm , in atm) and xCO 2 both measured by the SuperCO 2 instrument and corrected for water vapor partial pressure (p H 2 O , in atm) calculated as a function of temperature assuming saturation (Dickson et al., 2007): pCO 2 was also corrected for the temperature difference between SST and the equilibrator as in Pierrot et al. (2009).CO 2 flux (in mol m 2 d 1 ) was calculated as the difference between aqueous and atmospheric pCO 2 , multiplied by the gas transfer velocity (k, in m d 1 ), a function of the unitless Schmidt number Sc (Ho et al., 2006), mean squared neutral wind speed at 10 m above the sea surface (<U 2 10 n >, in m 2 s 2 ), and K o , the solubility of CO 2 in water (in mol m 3 atm 1 ) (Weiss, 1974).Positive values of CO 2 flux indicate sea-to-air efflux.
CO 2 Flux = kK o (pCO 2 aq pCO 2 atm ) (2) We compared two sources of atmospheric CO 2 concentrations for calculation of CO 2 flux: underway-measured atmospheric pCO 2 measured every 2 hr by the SuperCO 2 instrument and atmospheric pCO 2 as measured at the WLEF/Park Falls Wisconsin tower (Desai, 2022).The WLEF/Park Falls time series was chosen for CO 2 air-sea flux calculations, as detailed in the Results.
There is considerable disagreement among parameterizations and measurements of k in inland waters.Varying empirical formulations of the wind speed dependence of k have been derived and applied in lakes (e.g., Atilla et al., 2011;Cole & Caraco, 1998;Oxburgh et al., 1991;Wanninkhof, 1992), while comparison with eddy covariance-derived values of k indicates that heat and water fluxes also play a role in determining gas flux rates in lakes (Heiskanen et al., 2014).The k parameterization in this study (Ho et al., 2006) was chosen on the grounds that Lake Superior can be understood similarly to marine environments, with a high range of wind speeds and large fetch which merit the quadratic wind dependence discussed by Wanninkhof (1992) (D. Ho, personal communication).The extent to which this parameterization may bias calculated CO 2 flux in Lake Superior is unclear, and may benefit from direct measurement of k in future studies.

Results
More than 6 × 10 6 measurements of xCO 2 in Lake Superior surface waters were assembled into a pCO 2 and CO 2 flux timeseries.These data spanned the lake's most significant hydrological regions, including shallow coastal zones, deep (maximum 406 m) waters, riverine outlets, and regions bordering significant human development (Figure 1).The most heavily-observed regions included the far western arm of Lake Superior and a cross-lake transect from Duluth to Sault Ste.Marie.Binning of pCO 2 and CO 2 flux data by grouping observations by date and 0.01°boxes yielded 1.3 × 10 4 observations.

Underway Timeseries Overview
Mean observed SST was 11.4°C with a median of 12.7°C.SST varied widely among and within cruises, ranging from a maximum of 23.5°C in July 2019 near the center of the Far Western Arm to a minimum of 0.45°C in April 2022 in the plume of the St. Louis River Estuary.Specific conductance (SC) ranged from a near-constant 99.44 µS cm 1 in unstratified offshore waters to values exceeding 170 µS cm 1 in the plume of the St. Louis River Estuary, displaying a mean of 99.6 µS cm 1 , a median of 99.3 µS cm 1 , and a standard deviation of 2.2 µS cm 1 .The timing of thermal stratification in Lake Superior varied widely among locations and years (Austin et al., 2022), so observations within 0.5°C of the temperature of maximum density of freshwater (3.98°C) were designated as unstratified.Stratification onset occurred between late June and August, depending on year and location (Figure 2a).Interannual meteorological variability, including the previous winter's thermal state, exerts considerable influence on Lake Superior's summer stratification (Austin & Allen, 2011), as indicated by the historically late stratification of Lake Superior in August 2022 (J.Austin, personal communication), as well as maximum winter ice coverage which varied from 94.9% in early 2019 to 22.6%, 50.6%, and 79.5% the following years (T.-Y.Yang et al., 2020).
Surface-water DIC and pH (free scale) were calculated from measured pCO 2 , SST, and an assumed A T of 840 μmol kg 1 (Figures 2d and 2e) with PyCO 2 SYS, using the carbonate constants of Waters et al. (2014).A T is given as μmoles of hydrogen ion equivalent kg 1 in this work, and is largely invariant in Lake Superior (Minor & Brinkley, 2022;Sandborn et al., 2023) except in regions with significant terrestrial influence; no A T -conductivity relationship for Lake Superior has been published, so A T was not estimated from underway data.
Calculated pH free exhibited a mean of 8.075 and standard deviation of 0.093, while calculated DIC exhibited a mean of 855.0 μmol kg 1 and standard deviation of 8.8 μmol kg 1 .Seasonal variation in DIC was evident as a summertime decrease on the order of 20 μmol kg 1 , followed by an autumn increase of c. 10 μmol kg 1 .The pH free distribution fell within the range of values given in Minor and Brinkley (2022), while the mean calculated DIC was 10-40 μmol kg 1 higher than observations given in Zigah et al. (2011) and Sandborn et al. (2023).The discrepancy may be due to interannual DIC increases, sampling bias in the latter two studies favoring regions or periods of lower DIC, contribution of organic alkalinity, uncertainty associated with equilibrium calculation, or error propagated from uncertainty in A T .Minor et al. (2019) Reported mean annual A T values for the 1996-2014 EPA GLNPO timeseries ranging from approximately 810 to 870 μmol kg 1 .Propagation of this range in A T through our pCO 2 observations yielded a calculated DIC range of 825-885 μmol kg 1 and a calculated pH range of 8.059-8.090.This sensitivity analysis produced DIC and pH ranges overlapping with previous measurements, such that some of the discrepancy between calculated and previously-reported DIC could be due to error in assumed A T .

Atmospheric CO 2
Daily mean shipboard atmospheric xCO 2 varied seasonally in concert with the CO 2 timeseries observed at the Park Falls/WLEF tower (Desai, 2022)  observed between the underway and Park Falls/WLEF time series within years, yet the underway atmospheric signal displayed a much larger variability.Several anomalies emerged in the underway atmospheric data.Atmospheric xCO 2 measurements in several cruises were consistently higher than expected despite nominal measurements of standard gases and sea surface xCO 2 .These cruises included extended periods of idling, and presumably detection of exhaust CO 2 by the underway system.In another two cruises in September 2022, atmospheric (but not sea surface) xCO 2 was depressed over a period of weeks for reasons related to a filter on the air inlet.Due to these discrepancies, we chose to use daily means of nearby Park Falls/WLEF tower hourly measurements of atmospheric xCO 2 with the expectation of a well-mixed atmosphere over these scales.The occurrence of most atmospheric underway xCO 2 measurements within a close approximation of the Park Falls/ WLEF timeseries validated this expectation.If atmospheric CO 2 were not well-mixed on these scales, an uncertain amount of error in calculated CO 2 flux could result, which further work should explore.

Wind Speed
Wind speed observed on Lake Superior (corrected for direction of travel) exhibited a skewed unimodal distribution with a peak at 4.5 m s 1 (Figure S2 in Supporting Information S1).Some bias may have been incurred by intentional planning of transects around inclement weather and targeting the ice-free season, so it was unclear how well these transects represented the true distribution of wind velocity above Lake Superior.The underwayobserved wind speed distribution in 2020 stood out from other years with a lower and irregular distribution; these transects were limited in time and space (Figure S1 in Supporting Information S1) and are less likely to represent the true distribution of wind speed over Lake Superior.Comparison of the underway wind speed distributions with those measured offshore at the Stannard Rock Lighthouse over the same periods (Figure S3 in Supporting Information S1) indicates that the underway-observed wind speed distribution closely approximated that of the whole season.
The wind speed distribution peaks observed from either source were lower than the global U 10n distribution peak of approximately 7 m s 1 in M. Yang et al. (2022), which may imply an underestimation of CO 2 flux as parameterized by dual-tracer models as in this research.The present scarcity of research on gas flux parameterization validity in large lake systems for which size, morphometry, and variable winds greatly influence gas flux magnitude and timing (Perolo et al., 2021;Schilder et al., 2013) does not yet allow exploration of similar biases in this research.
Gas transfer velocities (k) calculated from the underway wind distribution displayed a mean of 1.6 m d 1 , about half the mean ocean value of 3.3 m d 1 given by Broecker and Peng (1982) and supported by revised gas transfer velocity parameterizations (e.g., Ho et al., 2006;Wanninkhof, 2014).Given k, along with the 147 m mean depth of Lake Superior (Fuller & Shear, 1995), DIC and aqueous CO 2 concentration [CO * 2 ] (both in molal units), and Revelle Factor (RF, the unitless ratio of relative change in [CO * 2 ] per relative change in DIC), all produced via PyCO 2 SYS equilibrium calculation, the characteristic timescale, or e-folding time, of CO 2 equilibration in Lake Superior (τ CO 2 ) can be estimated (Zeebe & Wolf-Gladrow, 2001): τ CO 2 functions as an local gas exchange relaxation timescale describing a system seeking equilibrium; the extent to which the period in which a system undergoes gas exchange approaches τ CO 2 dictates the efficiency of the exchange (Jones et al., 2014).In the case that surface waters remain in contact with the atmosphere for periods meeting and exceeding τ CO 2 , equilibration approaches completion barring other processes (primary production, respiration, etc.) maintaining air-sea disequilibrium.For systems like Lake Superior in which these other processes are minimal during mixing (see Section 3.5), τ CO 2 gains further relevancy as a descriptor of the actual process of air-sea CO 2 equilibration as observed in this work, and may be used to compare the physicochemical drivers of mixing among systems.
During unstratified periods, mean RF was 26.9 ± 0.6, mean DIC was 867.0 ± 0.9 μmol kg 1 , and mean [CO 2 *] was 29.6 ± 0.8 μmol kg 1 (all ± s.d.).The resulting τ CO 2 during the unstratified period was 100.± 4 days.This period is much smaller than that of most of the surface ocean mixed layer, indicating relatively fast CO 2 equilibrium despite Superior's deeper mixed layer.This period is similar in magnitude to the duration of the twiceannual unstratified periods in December-January and May-July (though stratification phenology varies among years; Austin and Colman (2008); Woolway et al. ( 2021)), so it is reasonable to expect that on multiannual timescales, Lake Superior maintains near-atmospheric CO 2 equilibrium.This inference depends on lake stratification and wind velocity, both of which may shift with the changing climate (Xue et al., 2022).Climate change effects on lake thermal state and atmospheric circulation are likely to have complex effects on lake biogeochemistry which extend to CO 2 flux behavior changes (Desai et al., 2009).This brief calculation also implies that stratified waters and high wind events may exhibit a much smaller τ CO 2 ; assuming a mixing depth of 10 m and k of 6.4 m d 1 , which might be observed in a strong gale during the winter stratified season, τ CO 2 would decrease by a factor of c. 40.This thought experiment informed by our data emphasizes the importance of thermal structure and ephemeral wind events for determining large-scale C cycling in Lake Superior.

pCO 2 Variability
A continuous multiannual cycle of observed pCO 2 could not be constructed due to large gaps in the time series, so the 4 years of observations were collated by Julian day of year in order to illustrate seasonal CO 2 cycling variability, facilitate comparison to previous observations and models, and contrast Lake Superior with other systems.
LOESS regression of observations grouped by 0.01°boxes and date of observation highlighted seasonal cycling of sea surface pCO 2 which resembled the sinusoidal behavior modeled by Bennington et al. (2012) (Figure 3a).This regression used a data window spanning the 5% of observations on either side of regression points and 10 iterative fits.A local maximum was observed in the regression at 427 μatm in early June, while a local minimum was observed at 310 μatm in early September.Individual observations ranged from 199 to 1,056 μatm, with a mean value of 379 ± 60 μatm (±s.d.).
Heterogeneity was visible in the range of pCO 2 values observed on a given date, with super-and under-saturated conditions observed throughout the year.This high degree of spatial and/or temporal heterogeneity obscured the seasonal pCO 2 cycle.Additionally, the high concentration of transects in the riverine-influenced Western Arm of Lake Superior may not have represented open-water conditions prevailing in the remainder of the lake.Diel variability was examined as a potential source of bias, but no significant difference between daytime and nighttime pCO 2 was found (Text S1 in Supporting Information S1) in contrast with observations of large diel CO 2 cycling variability in many smaller lakes (Golub et al., 2023).
Confounded spatial and seasonal variability was partly separated by SC into "riverine" and "pelagic" regimes in order to isolate open-water seasonal variability.A cutoff SC value was defined by statistically significant departure from the surface SC distribution observed in unstratified periods.In every year of observation, springtime unstratified surface SC observations formed a narrow distribution with a mean of 99.44 µS cm 1 and a standard deviation of 0.64 µS cm 1 .This value was taken to represent the mean SC of the well-mixed freshwater lake.
Observations with SC 3 standard deviations greater than the unstratified period mean were considered riverinfluenced, as in this system, river inputs have higher ionic strength or SC due to interactions with the rocks and soil in their watersheds and due to anthropogenic inputs such as road salts (Weiler, 1978).This scheme decreased the noise around the seasonal trend of surface water pCO 2 in pelagic observations (Figure 3b) and highlighted heterogeneity associated with riverine-influenced observations (Figure 3c).Potential interferences with this classification included evaporation and precipitation, which would be expected to increase and decrease surface SC, respectively.For this reason, we elected not to construct a quantitative mixing relationship based on underway-measured surface SC and merely used it as a rough proxy for riverine influence.
In pelagic waters of Lake Superior during April-November the mean observed pCO 2 was 380 ± 53 μatm, while in river-influenced waters, the mean observed pCO 2 was 343 ± 38 μatm (±s.d.).The depression of riverine regime mean pCO 2 pCO 2 may have been due to promotion of primary production and CO 2 drawdown in nutrient-rich riverine-influenced waters (Minor et al., 2014;Sterner et al., 2020), yet contemporaneous observations of highly elevated pCO 2 point to a gradient of heterotrophy and autotrophy associated with riverine plumes (Delvaux, 2017).The pelagic pCO 2 cycle regression displayed a greater seasonal variability than the modeled time series of The two models presented by Bennington et al. differed in their treatment of primary production limitation, which resulted in the greatest differences after spring mixing, when this study's observations also displayed high variability.These differences in the shape of the post-mixing pCO 2 cycle may indicate a model shortcoming or a change in inorganic carbon system drivers over the preceeding two decades.
The observed increase in spring mixing period pCO 2 relative to that modeled by Bennington et al. (2012) was consistent with the magnitude of atmospheric CO 2 concentration increase (c. 2 ppm yr 1 , Keeling & Keeling, 2017) over the approximately 22 years separating the modeled period of Bennington et al. and these observations, as well as the direction of increase in Lake Superior surface water pCO 2 calculated from pH and A T over the period 1992-2019 by Minor and Brinkley (2022).The precise rate of increase of Lake Superior surface water pCO 2 over decadal timescales remains difficult to constrain, but its repeated acheivement of nearatmospheric equilibrium, along with radiocarbon measurements indicating rapid (<3 years) recycling of the DIC pool (Zigah et al., 2011), and a relatively brief characteristic time of CO 2 equilibration (Section 3.3), indicate SANDBORN AND MINOR that it mirrors atmospheric pCO 2 over interannual periods and may continue to do so.This atmospheric CO 2 invasion of Lake Superior occurs regardless of its net CO 2 air-sea flux due to its dynamic chemical equilibration with atmospheric CO 2 .
The magnitude of seasonal variability in Lake Superior surface pCO 2 was comparable to that of subtropical ocean regions (Bates, 2001), but shifted in the year.In terms of pCO 2 phenology, Lake Superior resembled the Arctic Ocean most closely, despite exhibiting a larger amplitude (Orr et al., 2022).Neighboring Lake Michigan displayed a minimum surface pCO 2 2-3 months earlier than Superior in a process model study (Pilcher et al., 2015).Scarcity of data from November-April prevented extrapolation to those periods, but previous, model-based studies suggest Lake Superior pCO 2 likely remains supersaturated or near-atmospheric equilibrium throughout that period (Bennington et al., 2012).Interannually-variable wintertime ice cover (White et al., 2012) may modify the expected CO 2 efflux.

Competing Drivers of pCO 2
Deconvoluting the pelagic pCO 2 cycle (Figure 3b) into inferred drivers sheds light on biogeochemical cycling in Lake Superior by isolating direct temperature control of pelagic surface pCO 2 from control by other drivers.
Changes in water temperature directly control the speciation of the inorganic carbon system via changes to chemical equilibria and solubility.By removing the direct influence of water temperature variability, the remaining variation can be ascribed to biological and physical processes including production, respiration, sea-air CO 2 gas flux, and allochthonous C inputs.The indirect influence of temperature on these processes was not normalized by this analysis to preserve the seasonality of the biological and physical drivers.
The method of Takahashi et al. (1993) was used to separate measured pCO 2 into thermal (pCO 2 T ) and nonthermal (pCO 2 non-T ) signals using observed SST (T ) and the temperature-driven variation of pCO 2 : These signals have been termed "temperature" and "biological"/"biophysical" effects in other studies for example, Minor et al. (2019) and Fassbender et al. (2018).This work uses the T/non-T terminology of Atilla et al. (2011).CaCO 3 dissolution and precipitation were neglected in this analysis because Lake Superior is vastly undersaturated with respect to CaCO 3 so precipitation and dissolution processes are not relevant.Overbars indicated arithmetic mean values in the literature source, but this study analyzed an incomplete annual time series of pCO 2 , so mean temperature T) and mean pCO 2 pCO 2 ) were adjusted to 1°C and 400 μatm to ensure convergence of the driver signals at the beginning of the observed period.The temperature partial derivative of ln (pCO 2 ) was calculated via PyCO2SYS, yielding an average value of 0.03606°C 1 for Lake Superior over the temperature range 0-20°C (Text S2 in Supporting Information S1).This temperature dependence is in good agreement with values used in previous studies (0.038°C 1 Atilla et al., 2011; 0.0384°C 1 Lynch et al., 2010).
Plotting the measured, thermal, and non-thermal pCO 2 signals illustrated the interplay of these competing drivers of pCO 2 in Lake Superior (Figure 5).Seasonal temperature effects were visible as the springtime increase and autumn decrease in pCO 2 T , opposed by the summertime dip in pCO 2 non-T .Measured pCO 2 lay suspended between the curves.The degree to which thermal versus non-thermal drivers control pCO 2 can be conceptualized as the vertical distance between the measured curve and its two drivers; in spring, measured pCO 2 was closely tied to pCO 2 T , indicating that most of the spring trend in pCO 2 was driven by seasonal warming.pCO 2 moved equidistant between drivers before dipping with the non-thermal curve through the summer, resulting in CO 2 influx presumably maintained by net biological C drawdown, which persisted until the end of the ice-free season.
Quantitatively, the ratio of thermal to non-thermal control of pCO 2 can be calculated (Fassbender et al., 2018;Takahashi et al., 2002) as which yielded a value of 1.1 using the regressions in Figure 5, indicating roughly equal thermal and non-thermal driver magnitudes over the ice-free period.Interestingly, this value aligns with that of the Atlantic Ocean at the approximate latitude of Lake Superior (Fassbender et al., 2018), which raises questions about latitudinal gradients in R T non-T 1 in inland waters compared to marine systems.Minor et al. (2019) found majority non-thermal control of calculated pCO 2 from discrete samples of Lake Superior surface water in 2014-2016 in agreement with the findings of Atilla et al. (2011).This deconvolution analysis highlights the importance of non-thermal drivers of C drawdown in the latter stages of the ice-free season acting to depress pelagic surface pCO 2 despite opposition from temperature effects.

CO 2 Flux Variability
CO 2 air-sea flux rates calculated from observed pCO 2 displayed seasonal variation similar to that of pCO 2 , but with a greater degree of variability within individual cruises.The mean calculated CO 2 flux rate was 2.7 ± 6.4 mmol C m 2 d 1 .In pelagic waters, the mean calculated flux rate was 1.9 ± 5.7 mmol C m 2 d 1 , while in riverine-influenced waters a lower value of 4.5 ± 7.5 mmol C m 2 d 1 was calculated (all ± s.d.).These values fall the range of instantaneous CO 2 fluxes observed or modeled in previous work (Bennington et al., 2012;Delvaux, 2017).The most extreme air-sea fluxes were calculated in mid-summer, when high wind speeds coupled with CO 2 -undersaturated surface waters to create high rates of CO 2 drawdown exceeding 70 mmol C m 2 d 1 .Differences between pelagic and riverine-associated observations were clearest in spring, when pelagic observations tended to favor CO 2 air-sea influx while riverive observations favored efflux.This correlated with observations of springtime allochthonous C loading associated with nearshore CO 2 evasion in Lake Superior (Minor et al., 2019).
Net CO 2 air-sea flux over the observed seasons was estimated via trapezoidal integration of median daily observations of instantaneous CO 2 flux collated by Julian day across the observed time domain: Julian day 100 (April 9 or 10) through 300 (November 26 or 27).The resulting values (Table 1) were both multiplied by the total area of Lake Superior (8.21 × 10 10 m 2 ) to yield total fluxes as it was unclear what fraction of the lake was considered pelagic versus riverine.Uncertainty in integrated CO 2 fluxes was determined by bootstrap random resampling of observations with replacement for 100 repetitions of the integration and given as the standard deviation of the ensemble net fluxes.These values serve as rough bounds for the mean net CO 2 flux of Lake Superior throughout the ice-free seasons of 2019-2022.The contrast between spring riverine regime instantaneous CO 2 efflux and integrated ice-free season CO 2 influx speaks to temporal variability in riverine effects on inorganic C cycling in Lake Superior.
The resulting ice-free period CO 2 influx on the order of 25 Gmol C (300 Gg C, 0.30 mol C m 2 ) was similar in magnitude but opposite in sign to the only fully-annual estimate of Lake Superior air-sea CO 2 flux: a mean net annual efflux of 16 Gmol C yr 1 (190 Gg C yr 1 , 0.19 mol C m 2 yr 1 ) over the period 1997-2001 (Bennington et al., 2012).The discrepancy is accounted for by winter supersaturation of surface pCO 2 .Assuming the veracity and comparability of the above values (perhaps an invalid assumption due to the shift in model curves suggested by Figure 4), an efflux of 41 Gmol C (492 Gg C, 0.50 mol C m 2 ) during Julian days 301-399 would be implied.The rough approximations of carbon budgets allowed by available observational CO 2 flux data continues to prohibit integration of Lake Superior into regional and global C budgets.There remains the possibility that the modeled annual CO 2 flux and this study's calculated sub-annual flux are not comparable due to two intervening decades of ecological and climate change, an under-constrained modeled pCO 2 cycle, and ongoing uncertainty about comparisons of measured versus calculated pCO 2 in Lake Superior.An updated observation-based and/or process model constrained by spatially-and temporally-comprehensive direct observations of pCO 2 and CO 2 flux is required for substantive comparisons of observed and modeled C cycling.

Discussion
Four years of surface pCO 2 measurements gathered on transects across Lake Superior were used to describe inorganic carbon system variability across temporal and spatial scales.Ice-free season (April-November) observations yielded a detailed account of the seasonal pCO 2 cycle, driven by thermal and non-thermal drivers acting in opposition to perturb surface pCO 2 from its interannually-repeated atmospheric equilibrium, resulting in sustained periods of CO 2 influx and efflux.Increases in surface pCO 2 over the last two decades illustrate that Lake Superior is undergoing CO 2 invasion in agreement with Phillips et al. (2015).Separating the lake into pelagic and riverine regimes on the basis of SC highlighted riverine influences on the inorganic C system which tended to increase spring pCO 2 observations and instantaneous air-sea CO 2 efflux, likely through a combination of DIC loading and promotion of respiration.
Comparison of this study's calculated areal CO 2 flux rates demonstrates the importance of Lake Superior relative to regional C budget terms.Previous eddy covariance work indicated Lake Superior's ice-free season CO 2 sink is likely less than that of the surrounding terrestrial ecosystem (Vasys et al., 2011).The mean CO 2 instantaneous airsea influx of 2.7 mmol C m 2 d 1 found in this study supports this finding when compared to the terrestrial net ecosystem production (NEP) of 19.6 mmol C m 2 d 1 (86.1 g C m 2 yr 1 ) for the Great Lakes hydrologic unit (Butman et al., 2016).This terrestiral NEP is still nearly three times larger than the mean annual Lake Superior CO 2 efflux of 6.3 mmol C m 2 d 1 modeled by Bennington et al. (2012).The absolute values of these three C flux rates remain less than an order of magnitude apart, demonstrating the importance of Lake Superior as an major player in the regional C cycle.This point is further emphasized by the sheer size of Lake Superior, which has a surface area 24% the size of the Great Lakes hydrologic unit, and a net annual Lake Superior CO 2 efflux (Bennington et al., 2012) 17% that of the Great Lakes hydrologic unit net lake CO 2 efflux of 9.2 Gmol C (Butman et al., 2016).Integration of measured CO 2 air-sea fluxes rates over the observed period indicated net CO 2 influxes of 24.3 ± 2.9 Gmol C (pelagic) and 28.0 ± 1.4 Gmol C (riverine), which are considered bounding values for the mean whole-lake ice-free period mean CO 2 flux during 2019-2022.The results of this study are in agreement with previous work indicating sustained periods of CO 2 influx during the summer stratified period contrasting with efflux during winter, (Bennington et al., 2012), yet the contribution of the Great Lakes to regional and global C budgets remains uncertain, as recent regional and global estimates of CO 2 evasion from lakes and reservoirs have presented a broad range of values with large margins of error.Cavallaro et al. (2018) produced a North Future work should seek to address some of the limitations of this work.In particular, a paucity of early Spring and late Fall data hindered analysis of periods at the extremes of the observation period, which could shed light on the effects of ice-off as a driver of ice-free season CO 2 flux (cf., Ahmed et al., 2019).As previously noted, there may be some bias in wind-parameterized gas transfer velocities associated with dual-tracer experiments (M.Yang et al., 2022), such that the gas transfer velocities calculated here may be underestimates by as much as 20%.
Small-scale variations in turbulence and eddy scales are also important and underexplored mechanisms affecting air-water gas transfer uncertainty in low-wind environments (Wang et al., 2015).Wind speed gas flux parameterization applications in large lakes may be an important area for future study.

Consequences of Increasing pCO 2
Among the most impactful findings of this research is the observation that Lake Superior surface pCO 2 maintains near-equilibrium with overlying atmosphere over multi-year periods.Temperature variability and biogeochemical processes drive seasonal departures of pCO 2 from atmospheric equilibrium (creating a likely net annual CO 2 efflux), yet surface water pCO 2 returns to a baseline state near atmospheric equilibrium on timescales shorter than a year.This fact has several significant consequences in a world of increasing atmospheric CO 2 concentration: First, the solubility pump of Lake Superior acts as a partial CO 2 sink which can be approximated by an equilibrium calculation: Assuming A T = 840 μmol kg 1 , SST = 3.98°C (temperature of maximum density during destratification), an initial pCO 2 = 400 μatm, and atmospheric ΔpCO 2 Δt 1 = 2.50 μatm yr 1 , then equilibrium calculation indicates ΔDIC Δt 1 = 0.184 μmol kg 1 yr 1 , which is multiplied by the approximate mass of Lake Superior (1.21 × 10 17 kg) to give a CO 2 storage of 22.3 Gmol C yr 1 (267 Gg C yr 1 ) due to increasing atmospheric CO 2 alone.This storage is characteristic of any body of water maintaining CO 2 equilibrium with a nonsteady-state atmosphere.It acts alongside C sources (e.g., DIC loading) and sinks (e.g., C burial) to compose the net annual C budget of Lake Superior.Development of an annual net CO 2 flux using expanded observational and modeling capabilities may yield insights on all of these contributors.If atmospheric pCO 2 were stable, then Superior's annual net CO 2 efflux could be larger than it is today, mirroring the case of the pre-industrial global ocean, which likely acted as a CO 2 source instead of a sink (Cartapanis et al., 2018).
Second, Lake Superior's water chemistry will undergo changes as a result of consistently-higher pCO 2 .Its weak CO 2 buffer (RF 25-30 in this work, compared to marine values 8-16; Sarmiento and Gruber (2006)) and absence of sediment carbonate buffer (unlike neighboring Lake Michigan) result in relatively high sensitivity to atmospheric CO 2 acidification.The outcomes of hypothesized lake acidification mirror those in the ocean: decreasing pH and CaCO 3 saturation states, impacts on primary producer communities and carbonate-shelled organisms such as Dreisseina mussels, changes to metal ion activities, and other phenomena with potentially detrimental ecosystem effects (Doney et al., 2009).Trends in A T and temperature may modify the speciation (e.g., [CO 2 3 ] , pH) of the inorganic carbon system as well as the seasonal and spatial expression of the surface water pCO 2 cycle, but not the surface pCO 2 of a system at equilibrium with the atmosphere.
Third, efforts to observe Lake Superior's inorganic C system must capture a greater fraction of the annual cycle and spatial variability to constrain these changes.The twice-annual time series of chemical parameters (including glass electrode pH and Gran titration alkalinity) collected by US EPA Great Lakes National Program Office includes samples over a broad spatial scale, during periods of mean CO 2 efflux (April-May) and influx (August-September) but fails to observe intervening periods which provide necessary context for interannual variability of the annual pCO 2 cycle.Undersampling a complex signal like inorganic C chemistry delays detection of climate change effects (Carter et al., 2019).A more complete picture of biogeochemical parameters is sorely needed during the present period of climate change and ecological disruption.This gap in observational capabilities can be addressed by a sustained campaign of higher-quality, higher-frequency measurements of inorganic C parameters in the Laurentian Great Lakes.

Observational Challenges and Opportunities
Environmental and instrumental challenges limit deployment of underway pCO 2 systems as tools for biogeochemical observation on large lakes like Superior.These instruments describe only a small fraction of a water body at any given time, which complicates efforts to generalize results to the system as a whole.A network of similar sensors equipped on moorings, vessels of opportunity, and other vehicles (e.g., drifters, saildrones, wavegliders) may be suited for more comprehensive observation.Seasonal ice cover limits winter deployment of autonomous sensors, and has long acted as a blinder focusing scientific attention on more accessible seasons (Block et al., 2018).Novel observation platforms designed to observe under-ice pCO 2 (DeGrandpre et al., 2019;Lee et al., 2022) demonstrate the potential to expand the horizons of inorganic C observation in seasonally icecovered lakes.Direct measurements of gas flux may also be obtained by eddy covariance towers in the vicinity of the Great Lakes (Shao et al., 2015).This research grappled with problems of bias in transect data due to overrepresentation of certain regions in space (the far western lake) and time (summer).Although these problems were partially addressed by regression analysis and separation of pelagic and riverine regimes, future work should consider other drivers of spatial and temporal heterogeneity, for example,: dissolved organic matter and chlorophyll measured by in-situ instruments or remote sensing (e.g., Lohrenz et al., 2018;Sims et al., 2023).Expanded monitoring of pCO 2 and related chemical properties in the Laurentian Great Lakes provides a fruitful avenue for observation and modeling of CO 2 budgets in the world's largest surface freshwater resource.

Conclusions
This study provided the most comprehensive observations to date of surface pCO 2 variability in Earth's largest freshwater lake by area and demonstrated techniques for inferring C cycling drivers in an understudied system.As the present perturbation of Earth's C cycle continues, the need for such knowledge to inform water and climate policy will grow apace, requiring continuing innovation of observational and modeling capabilities.This is as true for the Laurentian Great Lakes as for the African Rift Lakes and other understudied surface waters of the world.
A spatially-comprehensive, fully annual CO 2 flux budget is not achievable with the data presented here because of spatial and temporal gaps in the time series presented.Future work must perform more observation of neglected regions in space and time, extrapolation to unobserved domains, and generalization of observed fluxes and drivers by modeling efforts.To this end, we recommend further development of observational strategies such as underway data collection, moored and autonomous instrumentation, remote sensing, and winter limnology techniques to better constrain CO 2 flux in Superior and other large lake systems.Efforts to resolve the modeled C budgets of the Great Lakes will benefit from a greater number of CO 2 measurements to constrain and correct models (Gloege et al., 2022).Insights into the balance of productivity and respiration may result from pairing a large pCO 2 survey with measurements of other biogeochemical tracers such as dissolved oxygen (Evans et al., 2022) or primary productivity (Sterner, 2010).As ice cover of temperate lakes declines with climate change, the period amenable to transects of seasonally ice-covered lakes will grow.This disappearance of the ice cover regime is among driving forces of the sub-discipline of winter limnology, which studies a vanishing environment (Ozersky et al., 2021).It is unclear how changes in ice cover will affect annual pCO 2 fluxes in these changing lakes systems.Spatially-and temporally-comprehensive observations of element cycling in these large lakes hint at the depth and complexity of biogeochemical functions responding and feeding back to a changing planet.

Figure 1 .
Figure 1.Underway measurement density transects 2019-2022, visualized as the number of occupations of approximately 5 km squares.The number of days of observation ranged from 0 to nearly 600.The cities of Duluth and Sault Ste.Marie, between which multi-lake transects traverse, are indicated by red triangles.The Park Falls/WLEF tower is denoted by a black square.

Figure 2 .
Figure 2. Sea surface temperature, pCO 2 , calculated CO 2 flux, calculated dissolved inorganic carbon (DIC), and calculated pH free observed in 0.01°boxes on transects of Lake Superior, 2019-2022.Median values for each day of observation are connected by a gray line.(a) The 3.98°C temperature of maximum density is indicated by a dotted line, along which lie unstratified conditions, highlighted in red.Depressed springtime surface temperatures of 2022 are visible as a delayed warming trend.(b) The Park Falls/WLEF time series is displayed as a dotted line separating observations of CO 2 supersaturation and undersaturation.(c) The division of CO 2 efflux versus influx is indicated by a dotted line.(d) DIC as calculated from pCO 2 and assumed total alkalinity (A T ) = 840 μmol kg 1 (e) pH (free scale) as calculated from pCO 2 and assumed A T = 840 μmol kg 1 .

Figure 3 .
Figure 3. Seasonal variation collated from pCO 2 observations grouped by 0.01°squares and date during transects of Lake Superior over 2019-2022.Black dashed lines represent LOESS regressions of each time series.

Figure 4 .
Figure 4. Median daily observations of pelagic surface water pCO 2 observed during 2019-2022 compared with Model 1 from Bennington et al. (2012), which described mean lake surface pCO 2 1997-2001.A 44 μatm adjustment of Model 1 to account for 22 years' atmospheric CO 2 increase (assuming 2 μatm yr 1 ) aligned spring and mixing season modeled results with contemporary observations.

Figure 5 .
Figure 5. Deconvolution of median daily measured sea surface pCO 2 (blue dashed line) into non-thermal (green dash-dot line) and thermal (red dotted line) drivers.Seventh-order power function regressions are shown as visual aids, and their equations are given in Supporting Information S1.
Bennington et al. (2012)(Figure4).Annual pCO 2 summer minima and spring maxima were approximately 330 and 400 μatm in Model 1 of that work, compared to 313 and 428 μatm in this study's regression of collated pelagic observations from 2019 to 2022.Bennington et al. modeled surface water equilibrium with an atmospheric pCO 2 of 360 μatm at the end of a mixing period spanning late April-late June 1997-2001.At the end of destratification in this (2019-2022) study, a mean surface water pCO 2 of 430 ± 30 μatm (±s.d.) was observed, which was indistinguishable from contemporaneous atmospheric pCO 2 .

Table 1
Lauerwald et al. (2023)of CO 2 Over the Air-Water Interface of Lake Superior Ascribed to Pelagic and Riverine Chemical Regimes for Julian Days 100-300 Note.Total flux values assume the whole lake surface area is either similar to pelagic or to riverine values and thus give the range within which total CO 2 lake flux is likely to fall.Uncertainties are given as standard deviations propagated via bootstrap resampling with replacement for 100 repetitions.Negative signs indicate influx.lakeandreservoir CO 2 efflux of 122 Tg C (10.1 Tmol C) yr 1 ± 100% which neglected the Laurentian Great Lakes, whileLauerwald et al. (2023)reported a North American natural lake CO 2 evasion of 337 Tg C (28.1 Tmol C) yr 1 with a range of 202-459 Tg C yr 1 .These recent studies represent significant advancements in the science of inland water carbon cycling, but their uncertainties indicate the potential for improvements in CO 2 cycling observation as demonstrated in this work to contribute to improved regional and global C budgets. American