Spatial variability in surface-water pCO2 and gas exchange in the world's largest semi-enclosed estuarine system: St. Lawrence Estuary (Canada)

Abstract

The incomplete spatial coverage of CO${}_{\mathrm{2}}$ partial pressure ($p$CO${}_{\mathrm{2}})$ measurements across estuary types represents a significant knowledge gap in current regional- and global-scale estimates of estuarine CO${}_{\mathrm{2}}$ emissions. Given the limited research on CO${}_{\mathrm{2}}$ dynamics in large estuaries and bay systems, as well as the sources of error in the calculation of $p$CO${}_{\mathrm{2}}$ (carbonic acid dissociation constants, organic alkalinity), estimates of air–sea CO${}_{\mathrm{2}}$ fluxes in estuaries are subject to large uncertainties. The Estuary and Gulf of St. Lawrence (EGSL) at the lower limit of the subarctic region in eastern Canada is the world's largest estuarine system, and is characterized by an exceptional richness in environmental diversity. It is among the world's most intensively studied estuaries, yet there are no published data on its surface-water $p$CO${}_{\mathrm{2}}$ distribution. To fill this data gap, a comprehensive dataset was compiled from direct and indirect measurements of carbonate system parameters in the surface waters of the EGSL during the spring or summer of 2003–2016. The calculated surface-water$p$CO${}_{\mathrm{2}}$ ranged from 435 to 765 $\mathrm{\mu }$atm in the shallow partially mixed upper estuary, 139–578 $\mathrm{\mu }$atm in the deep stratified lower estuary, and 207–478 $\mathrm{\mu }$atm along the Laurentian Channel in the Gulf of St. Lawrence. Overall, at the time of sampling, the St. Lawrence Estuary served as a very weak source of CO${}_{\mathrm{2}}$ to the atmosphere, with an area-averaged CO${}_{\mathrm{2}}$ degassing flux of 0.98 to 2.02 mmol C m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ (0.36 to 0.74 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$. A preliminary analysis revealed that respiration (upper estuary), photosynthesis (lower estuary), and temperature (Gulf of St. Lawrence) controlled the spatial variability in surface-water $p$CO${}_{\mathrm{2}}$. Whereas we used the dissociation constants of Cai and Wang (1998) to calculate estuarine $p$CO${}_{\mathrm{2}}$, formulations recommended for best practices in open ocean environments may underestimate $p$CO${}_{\mathrm{2}}$ at low salinities, while those of Millero (2010) may result in overestimates.

Introduction

Although estuaries occupy a very small fraction (0.2 %) of the global ocean surface area, their CO${}_{\mathrm{2}}$ emissions are disproportionately large compared with CO${}_{\mathrm{2}}$ exchanges between the open ocean and the atmosphere (Bauer et al., 2013). With an estimated global efflux of 0.10–0.15 Pg C yr${}^{-\mathrm{1}}$ (Chen et al., 2013; Laruelle et al., 2013), estuarine CO${}_{\mathrm{2}}$ degassing is thought to counterbalance CO${}_{\mathrm{2}}$ uptake on the continental shelves (Chen and Borges, 2009; Laruelle et al., 2010; Cai, 2011). Almost every estuary on Earth, for which data are available, is generally supersaturated with CO${}_{\mathrm{2}}$ with respect to the atmosphere (Cai and Wang, 1998; Frankignoulle et al., 1998; Borges, 2005; Borges et al., 2005, 2006; Chen and Borges, 2009; Laruelle et al., 2010; Cai, 2011; Chen et al., 2012; Bauer et al., 2013; Chen et al., 2013; Regnier et al., 2013), with CO${}_{\mathrm{2}}$ partial pressures ($p$CO${}_{\mathrm{2}})$ ranging from 400 to 10 000 $\mathrm{\mu }$atm (in contrast, the atmospheric $p$CO${}_{\mathrm{2}}$ in coastal zones was approximately 360–385 $\mathrm{\mu }$atm in the year 2000) (Cai, 2011). Although estuaries are generally net sources of CO${}_{\mathrm{2}}$, there is considerable variability and uncertainty in estimates of their CO${}_{\mathrm{2}}$ emissions, reflecting the limited spatial and temporal coverage of $p$CO${}_{\mathrm{2}}$ measurements in estuaries as well as their heterogeneous nature (hydrological and geomorphological differences, differences in magnitude and stoichiometry of carbon and nutrient inputs) (Bauer et al., 2013; Regnier et al., 2013).

Estuaries are geochemical reaction vessels through which continentally weathered organic matter and inorganic nutrients must pass to enter the coastal ocean (Kaul and Froelich, 1984). Horizontal transport is controlled by a set of physical attributes (tides, wind, bathymetry, basin geography, and river flow) that determine the estuarine filter function (Cloern, 2001). The longer the freshwater flushing (or turnover) time of the estuary, the more opportunity there is for water-column biological activity, benthic exchanges, and particle-dissolved phase interactions to influence its biogeochemistry (Statham, 2012). Dissolved inorganic carbon (DIC) enrichments and $p$CO${}_{\mathrm{2}}$ supersaturations observed in estuaries can be mainly attributed to the in situ microbial degradation of internally and externally supplied organic carbon and the lateral transport of inorganic carbon from rivers, coastal wetlands, and groundwaters (Bauer et al., 2013).

In strongly tidal (macrotidal) systems, long water and particle residence times (on the order of weeks to months; Middelburg and Herman, 2007) allow for the extensive modification and degradation of particulate organic carbon during estuarine transport (Borges et al., 2006; Chen and Borges, 2009). In the absence of seasonal or permanent water stratification, the decoupling between production and degradation of organic matter at and below the surface, respectively, does not occur, resulting in less efficient export of DIC (Borges, 2005). Strongly tidal estuaries also tend to exhibit lower levels of photosynthetic activity (Monbet, 1992) and carry greater suspended particulate matter loads within their high-turbidity regions (Uncles et al., 2002; Middelburg and Herman, 2007) wherein suspended particles and organic-rich aggregates serve as hot spots of microbial recycling (Statham, 2012). Field measurements suggest that 10 % of the total CO${}_{\mathrm{2}}$ emissions from the inner estuary of macrotidal systems is sustained by the ventilation of riverine CO${}_{\mathrm{2}}$, whereas 90 % is due to local net heterotrophy (Borges et al., 2006) fueled by inputs of terrestrial- and riverine-algae-derived (planktonic) detritus and, in populated areas, sewage (Chen and Borges, 2009). In estuaries with long freshwater residence times, the riverine CO${}_{\mathrm{2}}$ will be fully ventilated to the atmosphere within the estuary, and the total CO${}_{\mathrm{2}}$ emissions can be attributed to net heterotrophy (Borges and Abril, 2011).

North American estuaries rank first in terms of global estuarine surface area (41 %) but account for the lowest numerically averaged CO${}_{\mathrm{2}}$ flux per unit area (12 %) among all continents (Chen et al., 2013). These estimates are subject to large uncertainties due to data paucity. A recent synthesis by Regnier et al. (2013) highlighted the meagre spatial coverage of estuarine $p$CO${}_{\mathrm{2}}$ measurements, particularly along the eastern Canadian seaboard. Ironically, the Estuary and Gulf of St. Lawrence (EGSL) in eastern Canada is the largest semi-enclosed estuarine system in the world and is among the world's most intensively studied estuaries (El-Sabh and Silverberg, 1990), but was left unmentioned in recent global (Cai, 2011; Chen et al., 2012, 2013) and regional (Laruelle et al., 2015) data compilations. Furthermore, previous estuarine CO${}_{\mathrm{2}}$ studies have focused on river-dominated estuaries, whereas there has been limited research on CO${}_{\mathrm{2}}$ dynamics in large estuaries and bays (Joesoef et al., 2015), i.e., marine-dominated systems, including the areas of mixing at sea (outer estuaries or river plumes) (Borges et al., 2005). A comparative study by Jiang et al. (2008) revealed large differences in CO${}_{\mathrm{2}}$ degassing between non-riverine and river-dominated estuaries and, more recently, Koné et al. (2009), Maher and Eyre (2012), and Cotovicz Jr. et al. (2015) reported low CO${}_{\mathrm{2}}$ uptake by strongly stratified and/or marine-dominated systems. On the US east coast, the ratio of non-riverine (flushed by tidal action and receiving minimum freshwater inputs) to river-dominated estuaries is nearly 1 : 1, demonstrating the geographic importance of coastal estuaries and bays on the eastern seaboard of North America (Cai, 2011).

Map of the greater St. Lawrence system, including the chain of Great Lakes, the St. Lawrence River, the upper St. Lawrence Estuary (USLE), the lower St. Lawrence Estuary (LSLE), and the Gulf of St. Lawrence (GSL) (Working Group on the State of the St. Lawrence Monitoring, 2014, with permission).

[Figure omitted. See PDF]

The large-scale (width often considerably greater than the internal Rossby radius; Cyr et al., 2015), macrotidal (mean tidal range greater than 2–4 m; Monbet, 1992) St. Lawrence Estuary is an excellent analogue of marine-dominated systems. Throughout its length, the full spectrum of oceanic variability can be found (Mertz and Gratton, 1990). Moreover, the basin characteristics and water transport timescales of the St. Lawrence Estuary (SLE) provide an almost ideal natural laboratory for geochemical studies. Its surface waters have a renewal time of several months while its bottom waters take several years to replenish, allowing for a comparison of spatial and temporal variations in physical and chemical properties (El-Sabh and Silverberg, 1990). Given its bimodal bathymetry, the St. Lawrence Estuary also permits the investigation of biogeochemical processes in two types of estuaries: (1) the shallow partially mixed upper estuary, where physical mixing and abiotic processes dominate and (2) the deep stratified lower estuary, where biological cycling and oceanic processes prevail (Yeats, 1990). As yet, no systematic study of the CO${}_{\mathrm{2}}$ dynamics in the SLE has been published and, hence, the present study provides the first comprehensive description of its mixed-layer carbonate chemistry, including (1) a multi-year compilation of springtime and summertime $p$CO${}_{\mathrm{2}}$ calculated from direct measurements of pH and alkalinity; (2) an area-averaged estimate of the air–sea CO${}_{\mathrm{2}}$ gas flux; and (3) an analysis of the relative importance of thermodynamic (temperature) and biological (photosynthesis and respiration) processes in controlling the spatial variability in surface-water $p$CO${}_{\mathrm{2}}$.

Materials and methods

Study area – St. Lawrence Estuary and Gulf

The greater St. Lawrence system (Fig. 1) connects the chain of Great Lakes, the second largest terrestrial freshwater reservoir in the world, to the Atlantic Ocean (Yang et al., 1996). With a drainage basin of approximately 1.32 million km${}^{\mathrm{2}}$, the St. Lawrence River channels the second largest freshwater discharge (11 900 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}})$ on the North American continent, second only to that of the Mississippi River (El-Sabh and Silverberg, 1990). The catchment area geology is dominated by silicate rocks of the Precambrian Shield and carbonates of the Paleozoic lowlands, whose components influence the downstream evolution of river chemistry (Yang et al., 1996). The erosion of the carbonate rocks of the drainage basin is practically constant, as the quantity of bicarbonate ions carried by the river varies little from season to season (Pelletier and Lebel, 1979). On a yearly basis, between 15 and 20 % of the outflow of dissolved inorganic carbon from the St. Lawrence River into its estuary originates from the system's tributaries (e.g., Ottawa River, Mascouche River) while 80–85 % is from the Great Lakes (Hélie et al., 2002). The SLE begins at the landward limit of the saltwater intrusion near Île d'Orléans ($\sim$ 5 km downstream of Québec City) and stretches 400 km seaward to Pointe-des-Monts where it widens into the Gulf of St. Lawrence (GSL) (El-Sabh and Silverberg, 1990), a semi-enclosed sea with an area of approximately 240 000 km${}^{\mathrm{2}}$ (Dufour and Ouellet, 2007) connected to the Atlantic Ocean through Cabot Strait and the Strait of Belle Isle (Coote and Yeats, 1979).

Traditionally, the SLE is divided into two segments based on its bathymetry and hydrographical features (Ingram and El-Sabh, 1990). The upper St. Lawrence Estuary (USLE), from Île d'Orléans, where the estuarine circulation begins, to Tadoussac, near the mouth of the Saguenay Fjord, covers an area of 3470 km${}^{\mathrm{2}}$. It is relatively narrow (2 to 24 km wide) and mostly shallow (depths less than 30 m; d'Anglejan, 1990) and features an uneven, fairly complex bottom topography characterized by several disconnected channels and troughs separated by ridges and islands (El-Sabh and Murty, 1990). Topographically modified flows give rise to very large tidal ranges and currents (up to 10 m and 3 m s${}^{-\mathrm{1}}$, respectively; Mertz and Gratton, 1990). In this tidally energetic region, wind mixing is 1 to 2 orders of magnitude smaller than tidal mixing (Painchaud et al., 1995). Owing to the resuspension of bottom sediments (tide, wind, or wave generated) and the net non-tidal estuarine circulation (d'Anglejan and Smith, 1973), a well-developed turbidity maximum stretches between Île d'Orléans and Île-aux-Coudres (Painchaud and Therriault, 1989), where suspended particulate matter concentrations vary from 10 to more than 200 mg L${}^{-\mathrm{1}}$ (Silverberg and Sundby, 1979). The sources of particulate organic matter (POM) in the estuary are still debated (Gearing and Pocklington, 1990). Carbon isotope studies indicate that less than half of the POM is derived from terrestrial sources (Pocklington and Leonard, 1979) and is quite refractory to biodegradation (Lucotte et al., 1991), whereas the major contributor to POM is believed to be fresh organic matter, i.e., living or recently living material, of river-borne origin (Tan and Strain, 1983; Hélie and Hillaire-Marcel, 2006). During the spring freshet in April–May, when freshwater discharge delivers 40 % of the annual solid inputs to the estuary, the input of terrigenous POM is equivalent to the average POM kept in suspension in the turbidity maximum (Lucotte, 1989).

The lower St. Lawrence Estuary (LSLE) is fairly unique in that its character is more oceanic than most estuaries due to its grand size in all three dimensions and unimpeded connection with continental shelf and slope waters of the northwest Atlantic Ocean (El-Sabh and Silverberg, 1990). Relative to the USLE, the LSLE is much larger (9350 km${}^{\mathrm{2}}$; d'Anglejan, 1990), wider (30 to 50 km), and deeper ($\sim$ 300 m), and it displays a smoother, less variable bottom topography. Tidal currents are weaker (on the order of 30 cm s${}^{-\mathrm{1}}$ or less; Mertz and Gratton, 1990) and, under these less turbulent conditions, the lower estuary is the major sink of continental inputs to the St. Lawrence system. Most ($\sim$ 75 %) of the terrigenous POM carried by the St. Lawrence River is deposited on the LSLE floor (Lucotte et al., 1991). The dominant bathymetric feature of the LSLE is the Laurentian Channel (or Trough), a deep, central U-shaped glaciated valley that extends 1240 km from the eastern Canadian continental shelf break through the GSL and into the LSLE (d'Anglejan, 1990). The termination (head) of the Laurentian Channel at an abrupt and shallow sill near Tadoussac marks the region of transition between the upper and lower estuaries and is an area of complex tidal phenomena (Gratton et al., 1988). Due to rapid shoaling, tidal movements (e.g., internal tides and strong flows over the steep sill) locally generate significant mixing of surface freshwater with cold, nutrient-rich waters from the intermediate and deep layers of the GSL, resulting in a fertile surface layer that flows continuously seaward (Coote and Yeats, 1979; Saucier and Chassé, 2000) and sustains important feeding habitats for several large marine mammals (Dufour and Ouellet, 2007). The lower estuary's seaward outflow, together with the Gaspé Current, a rapidly moving coastal jet, are a major input of nutrients and zooplankton to the near-surface waters of the GSL (Coote and Yeats, 1979; Plourde and Runge, 1993). Mesoscale features such as coastal jets, internal Kelvin waves, baroclinic eddies, and unstable waves are all possible due to strong Coriolis effects (Ingram and El-Sabh, 1990).

The lower St. Lawrence Estuary is one of the most laterally stratified estuaries in the world (Larouche et al., 1987) and is also strongly vertically stratified. During summertime, the SLE can be described as a three-layer system on the basis of its thermal stratification (Gratton et al., 1988). Each spring, a new surface layer flow is initiated by the freshwater runoff from the St. Lawrence River, Saguenay Fjord, and rivers on the north shore of the estuary (Dufour and Ouellet, 2007; see Fig. 1). Discharge from the St. Lawrence River (mean annual discharge of 10 000 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$, peaking at 15 000 m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ during the spring freshet; Painchaud and Therriault, 1989) provides about 80 % of the total freshwater input to the estuary (Ingram and El-Sabh, 1990), whereas the combined runoff from the Saguenay and Manicouagan rivers accounts for most of the remainder (Tee, 1990). The warm and relatively fresh surface layer (0 to 30 m) overlies the cold intermediate layer or CIL (30–150 m deep; $S_P=\mathrm{32.0}$ to 32.6) that is formed by advection of the GSL's wintertime surface mixed layer (Galbraith, 2006). Below the CIL, a warmer (2 to 6 ${}^{\circ }$C) and saltier ($S_P=\mathrm{33}$ to 35) bottom layer ($>$ 150 m deep), originating from the mixing of western-central Atlantic and Labrador shelf waters that intrude at depth primarily through Cabot Strait, flows sluggishly landward ($\sim$ 0.5 cm s${}^{-\mathrm{1}}$; Bugden, 1988) toward the head region of the Laurentian Channel (Saucier et al., 2003; Gilbert et al., 2005).

Water-column sampling and analytical procedures

Water samples were collected aboard the RV Coriolis II during ten research cruises within the EGSL in the ice-free spring or summer season between 2003 and 2016. Water sampling was conducted mainly along the central axis of the SLE and the Laurentian Channel. The sampling locations are shown in Fig. 2. Samples were taken from discrete depths throughout the water column, typically at 3, 20, 50, 70, 100 m, and 50 m intervals to the bottom (or within 10 m of the bottom). A comprehensive dataset was compiled from field or laboratory measurements of the following physical–chemical properties: temperature ($T$), practical salinity ($S_P)$, pH${}_{\mathrm{NBS}}$ and/or pH${}_{\mathrm{T}}$, total alkalinity (TAlk), dissolved inorganic carbon (DIC), soluble reactive phosphate (SRP), and dissolved silicate (DSi), nitrate (NO${}_{\mathrm{3}})$, and oxygen (DO).

Map showing the four principal subdivisions of the study area and the sampling locations (markers). Water samples were collected during 10 spring and summer research cruises: July 2003, June 2006, May 2007, July 2007, June 2009, July 2009, July 2010, May 2011, June 2013, and May 2016. The estuary, from the landward limit of the saltwater intrusion near Île d'Orléans ($\sim$ 5 km downstream of Québec City) to Pointe-des-Monts, where the coastline diverges, extends $\sim$ 400 km and covers a total surface area of $\sim$ 12 820 km${}^{\mathrm{2}}$. The solid gold line follows the 200 m isobath of the Laurentian Channel. The vertical dotted lines delineate the five segments of the estuary used for the calculation of the area-averaged CO${}_{\mathrm{2}}$ flux.

[Figure omitted. See PDF]

$T$ and $S_P$ were determined in situ using the conductivity, temperature, and depth (CTD) probe (Sea-Bird SBE 911) mounted on the sampling rosette. The temperature probe was calibrated by the manufacturer, whereas the conductivity sensor was calibrated by the manufacturer and recalibrated using discrete salinity samples collected throughout the water column and analyzed on a Guildline Autosal 8400 salinometer calibrated with IAPSO standard seawater. Water samples destined for pH and TAlk measurements were transferred directly from the 12 L Niskin bottles mounted on the CTD-rosette system to, respectively, 125 mL plastic bottles without headspace and 250 mL glass bottles as soon as the rosette was secured onboard. In the latter case, a few crystals of HgCl${}_{\mathrm{2}}$ were added before the bottle was sealed with a ground-glass stopper and Apiezon^® Type-M high-vacuum grease.

pH was determined onboard at 25 ${}^{\circ }$C, potentiometrically on the NBS/NIST scale (infinite dilution convention, pH${}_{\mathrm{NBS}})$ for low-salinity waters ($S_P$ $<$ 5) and potentiometrically and/or colorimetrically on the total hydrogen ion concentration scale (constant ionic medium convention, pH${}_{\text{T}})$ for higher-salinity waters. Potentiometric pH measurements were carried out using a Radiometer Analytical^® GK2401C combination glass electrode connected to a Radiometer Analytical^® PHM84 pH/millivoltmeter. Prior to each measurement, the electrode was calibrated against three NIST-traceable buffer solutions: pH-4.00, pH-7.00, and pH-10.00 at 25 ${}^{\circ }$C. The electrode response to these buffers was then least-squares fitted to obtain the Nernstian slope. For $S_P$ $>$ 5, pH measurements were converted to the pH${}_{\mathrm{T}}$ scale using TRIS buffer solutions prepared at $S_P=\mathrm{5}$, 15, 25, or 35, for which the pH${}_{\mathrm{T}}$ was assigned at 25 ${}^{\circ }$C (Millero, 1986). Colorimetric pH measurements were carried out using a Hewlett-Packard UV-visible diode array spectrophotometer (HP-8453A) and a 5 cm quartz cell after thermal equilibration of the plastic sampling bottles in a constant temperature bath at 25.0 $\pm$ 0.1 ${}^{\circ }$C. Phenol red (Robert-Baldo et al., 1985) and $m$-cresol purple (Clayton and Byrne, 1993) were used as color indicators. The pH${}_{\mathrm{T}}$ of the water samples and buffer solutions were calculated according to the equation of Byrne (1987). The reproducibility of the pH measurements was typically better than $\pm$0.003.

TAlk was measured at McGill University using an automated radiometer (TitraLab865^®) potentiometric titrator and a Red Rod^® combination pH electrode (pHC2001). The dilute HCl titrant was calibrated prior, during, and after each titration session using certified reference materials provided by Andrew Dickson (Scripps Institution of Oceanography). Raw titration data were processed with a proprietary algorithm specifically designed for shallow end-point detection. The reproducibility of the method was better than 0.5 %.

Direct DIC measurements were carried out during an additional cruise in 2014 using a SciTech Apollo DIC analyzer. After being thermostated at 25 ${}^{\circ }$C, 1–1.5 mL of the sample was injected into the instrument's reactor, where it was acidified with 10 % H${}_{\mathrm{3}}$PO${}_{\mathrm{4}}$ and the evolved CO${}_{\mathrm{2}}$ was carried by a stream of pure nitrogen to a LI-COR infrared analyzer. A calibration curve was constructed using gravimetrically prepared Na${}_{\mathrm{2}}$CO${}_{\mathrm{3}}$ solutions, and the accuracy of the measurements was verified using certified reference material solutions provided by Andrew Dickson (Scripps Institute of Oceanography). The reproducibility of the measurements was typically on the order of 0.2 %. Results of the direct DIC measurements were used to assess the contribution of organic alkalinity to the total alkalinity in the upper estuary, as discussed in Sect. 2.3.2.

DO concentrations were determined by Winkler titration (Grasshoff et al., 1999) on distinct water samples recovered directly from the Niskin bottles. The relative standard deviation, based on replicate analyses of samples recovered from the same Niskin bottle, was better than 1 %. These measurements further served to calibrate the SBE-43 oxygen probe mounted on the rosette. For the determination of nutrient concentrations, aliquots of the water samples taken from the Niskin bottles were syringe filtered through a 0.45 $\mathrm{\mu }$m Millipore polycarbonate (MA) filter. DSi was measured onboard on the same day of sampling using the method described in Grasshoff et al. (1999). Water samples destined for NO${}_{\mathrm{3}}$ and SRP measurements were transferred into acid-washed 15 mL polyethylene and borosilicate tubes, respectively, quickly frozen, and stored at $-$20 ${}^{\circ }$C. Their concentrations were determined using standard colorimetric methods adapted from Grasshoff et al. (1999) with a SEAL AutoAnalyzer III at the Institut des Sciences de la Mer de Rimouski. The analytical detection limit was 0.04 $\mathrm{\mu }$M for NO${}_{\mathrm{3}}$, 0.05 $\mathrm{\mu }$M for SRP, and 0.1 $\mathrm{\mu }$M for DSi. Based on replicate analyses of the standards, the reproducibility of these measurements was typically 1 %.

The in situ pressure and density of the samples were calculated from the Thermodynamic Equation of Seawater – 2010 (TEOS-10) using the Gibbs SeaWater (GSW) Oceanographic Toolbox (MATLAB version 3.05; McDougall and Barker, 2011). All field measurements reported in $\mathrm{\mu }$mol L${}^{-\mathrm{1}}$ were converted to $\mathrm{\mu }$mol kg${}^{-\mathrm{1}}$ using the in situ density data.

Calculation of aqueous $p$CO${}_{\mathrm{2}}$

$p$CO${}_{\mathrm{2}}$ in mixed-layer waters

Aqueous $p$CO${}_{\mathrm{2}}$ ($p$CO${}_{\mathrm{2}}$(water)) is defined as the partial pressure of carbon dioxide in wet (100 % water-saturated) air that is in equilibrium with the water sample. Because direct $p$CO${}_{\mathrm{2}}$ measurements were not available from the RV Coriolis II cruises, $p$CO${}_{\mathrm{2}}$(water) ($\mathrm{\mu }$atm) and DIC ($\mathrm{\mu }$mol kg${}^{-\mathrm{1}})$ were calculated from the measured pH (total or NBS scale) and TAlk ($\mathrm{\mu }$mol kg${}^{-\mathrm{1}})$, at in situ temperature (${}^{\circ }$C), salinity ($S_P)$, and pressure (dbar), using the program CO2SYS (MATLAB version 1.1; van Heuven et al., 2011) and the carbonic acid dissociation constants ($K_{\mathrm{1}}$, $K_{\mathrm{2}})$ of Cai and Wang (1998) for estuarine waters. Wherever data were available, the contributions to TAlk from phosphate and silicate were included in the calculations. Although the $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$ formulations from Lueker et al. (2000) are recommended for best practices by Dickson et al. (2007), they are not suitable for the low-salinity conditions found in estuaries ($S_P$ $<$ 19) (Orr et al., 2015). The revised equations for $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$ from Cai and Wang (1998) are applicable over a larger range of salinities (0 to 40) and, thus, were used to examine the carbonate system in the estuarine waters of our study area.

This study focuses on the CO${}_{\mathrm{2}}$ dynamics in near-surface waters. To obtain individual data points of surface-water $p$CO${}_{\mathrm{2}}$ at each sampling location, the $p$CO${}_{\mathrm{2}}$ data in the surface mixed layer (SML) were averaged. The SML is the site of active air–sea interaction where heat and gases are exchanged directly with the atmosphere and within which physical (temperature, salinity, and density) and chemical (dissolved gases) properties are vertically homogeneous due to turbulent mixing (Sprintall and Tomczak, 1992). The lower limit to air–sea interaction, i.e., the mixed-layer depth (MLD), is demarcated by a pycnocline, a sharp density gradient that generally coincides with both a temperature gradient (thermocline) and a salinity gradient (halocline). Here, we determine the thickness of the SML using a density-based criterion which defines the MLD as the depth at which a threshold difference of 0.03 kg m${}^{-\mathrm{3}}$ from the sea surface occurs (de Boyer Montégut et al., 2004). In the following sections, surface-water $p$CO${}_{\mathrm{2}}$ will be taken to mean the SML-averaged $p$CO${}_{\mathrm{2}}$.

Sources of error in $p$CO${}_{\mathrm{2}}$ calculation

In this study, pH $/$ TAlk was used as the input combination to study the consistency between calculations of $p$CO${}_{\mathrm{2}}$ using different sets of carbonic acid dissociation constants. The importance of using appropriate formulations of $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$ in estuarine waters is shown by the discrepancies in the calculated $p$CO${}_{\mathrm{2}}$ values ($p$CO${}_{\mathrm{2}}$ at 15 ${}^{\circ }$C) at low salinities (Fig. 3). The percent difference between values calculated using the dissociation constants of Cai and Wang (1998) and those calculated using the best-practice constants of Lueker et al. (2000) was on average 3.07 %. At $S_P$ $<$ 19, differences between the calculated $p$CO${}_{\mathrm{2}}$ values were as large as 18.0 % (average difference of 6.88 %), whereas, at $S_P$ $>$ 19, the calculated values were in better agreement (only $\sim$ 1.48 % difference). The $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$ formulations of Millero (2010), the most recent set of constants proposed for estuarine waters ($S_P=\mathrm{1}$ to 50), yielded $p$CO${}_{\mathrm{2}}$ values that differed substantially from those derived using the constants of Cai and Wang (1998) at $S_P$ $<$ 19, with the largest divergence reaching 34.4 % (average difference of 13.6 %). At $S_P=\mathrm{0}$, the $p$CO${}_{\mathrm{2}}$ values calculated using the constants of Cai and Wang (1998) compared very well with those given by the Millero (1979) constants for freshwater (difference of only $\sim$ 0.08 %), whereas the $p$CO${}_{\mathrm{2}}$ values calculated using the constants of Millero (2010) showed very poor agreement with the freshwater results (differing by $\sim$ 34.3 %). These discrepancies highlight the need for new or revised measurements of the carbonic acid dissociation constants under estuarine conditions (in brackish waters) especially at $S_P$ $<$ 5. Studies that use the best-practice formulations of $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$ to calculate estuarine $p$CO${}_{\mathrm{2}}$ may underestimate CO${}_{\mathrm{2}}$ emissions at low salinities, whereas those that implement the Millero (2010) formulations may produce overestimates.

Differences in $p$CO${}_{\mathrm{2}}$ ($\mathrm{\Delta }$$p$CO${}_{\mathrm{2}})$ calculated using different published formulations of $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$, including Cai and Wang (1998) (C&W), Lueker et al. (2000) (Luek), Roy et al. (1993) (Roy), Millero (2010) (Mill), and Millero (1979) for pure water only ($S_P=\mathrm{0}$) (Fresh). All calculations were carried out at 15 ${}^{\circ }$C ($p$CO${}_{\mathrm{2}}$ at 15 ${}^{\circ }$C) with measured pH and TAlk. The constants of Lueker et al. (2000) are recommended for best practices by Dickson et al. (2007) and are applicable for waters at $S_P$ $>$ 19 (delineated by the vertical dashed line), whereas those of Roy et al. (1993) ($S_P=\mathrm{5}$ to 45) are recommended by Dickson and Goyet (1994). Both the constants of Cai and Wang (1998) and Millero (2010) have been proposed as more appropriate for the study of the carbonate system in estuarine waters.

[Figure omitted. See PDF]

Another potential source of error in the calculation of $p$CO${}_{\mathrm{2}}$ (pH $/$ TAlk) in low-salinity estuarine waters is the contribution of dissolved organic compounds to the total alkalinity. The contributions of borate, phosphate, and silicate species are taken into consideration in CO2SYS, whereas the magnitude of organic alkalinity (Org-Alk, or excess alkalinity) is usually assumed to be small or negligible, and it is simply ignored when using TAlk to calculate $p$CO${}_{\mathrm{2}}$ in open ocean waters. In riverine and coastal waters, however, the contribution of organic species to the TAlk can be significant (Yang et al., 2015). Rivers draining organic-rich soils and non-carbonate rocks have low DIC concentrations (a few hundred $\mathrm{\mu }$mol L${}^{-\mathrm{1}})$ that are often exceeded by dissolved organic carbon (DOC) concentrations (Abril et al., 2015). As discussed by Hunt et al. (2011), a significant contribution of Org-Alk (the organic acid anions in DOC) leads to an overestimation of calculated $p$CO${}_{\mathrm{2}}$ using any algorithm that accounts only for the contributions of inorganic species to TAlk. A comparison of the calculated TAlk (DIC $/$ pH) and the measured TAlk from the 2014 cruise reveals that Org-Alk is on the order of $-$20 $\mathrm{\mu }$mol kg${}^{-\mathrm{1}}$ for the St. Lawrence River end-member, whereas it is as high as $-$120 $\mathrm{\mu }$mol kg${}^{-\mathrm{1}}$ for the Saguenay River end-member (A. Mucci, personal communication, 2017). Given that TAlk exceeds $\sim$ 1000 $\mathrm{\mu }$mol kg${}^{-\mathrm{1}}$ throughout our study area and the Saguenay River contribution to the surface waters of the SLE is limited (at most $\sim$ 6 % at the head of the lower estuary; Mucci et al., 2017), consideration of the Org-Alk in the calculation of $p$CO${}_{\mathrm{2}}$ (pH $/$ TAlk) yielded values that were at most 1.9 % different from those uncorrected for Org-Alk. Bearing in mind the uncertainties in the $K_{\mathrm{1}}$ and $K_{\mathrm{2}}$ formulations as well as the analytical uncertainties, the influence of Org-Alk on the calculated $p$CO${}_{\mathrm{2}}$ (pH $/$ TAlk) did not represent a significant source of error.

Temperature normalization of $p$CO${}_{\mathrm{2}}$

The effect of temperature on aqueous $p$CO${}_{\mathrm{2}}$ is primarily the manifestation of changes in the solubility of CO${}_{\mathrm{2}}$ gas in water (Takahashi et al., 1993). The temperature dependence of $p$CO${}_{\mathrm{2}}$ in seawater, i.e., $\partial \ln \left(p{\mathrm{CO}}_{\mathrm{2}}\right)/\partial T=\mathrm{0.0423}$ ${}^{\circ }$C${}^{-\mathrm{1}}$, was determined experimentally by Takahashi et al. (1993) on a single North Atlantic surface-water sample with $S_P=\mathrm{35.380}$ under isochemical conditions. As this oft-used approximation for thermally induced changes in $p$CO${}_{\mathrm{2}}$ was derived from direct measurements in open ocean waters, we use a different approach to remove the temperature effect on the estuarine $p$CO${}_{\mathrm{2}}$ in our study area. The in situ $p$CO${}_{\mathrm{2}}$ values were normalized to the average surface-water temperature ($p$CO${}_{\mathrm{2}}$($\stackrel{\mathrm{‾}}{\mathrm{SST}}$), $\stackrel{\mathrm{‾}}{\mathrm{SST}}=\mathrm{7.82}$ ${}^{\circ }$C), using the temperature normalization method of Jiang et al. (2008) in which $p$CO${}_{\mathrm{2}}$ values are recalculated from the TAlk and DIC data at a common temperature. The results yielded a temperature coefficient of $\partial \ln \left(p{\mathrm{CO}}_{\mathrm{2}}\right)/\partial T=\mathrm{0.0402}$ ${}^{\circ }$C${}^{-\mathrm{1}}$ ($R^{\mathrm{2}}=\mathrm{0.99}$), in excellent agreement with that of Takahashi et al. (1993). The $p$CO${}_{\mathrm{2}}$ changes due to temperature deviations from $\stackrel{\mathrm{‾}}{\mathrm{SST}}$ ($\mathrm{\Delta }$$p$CO${}_{\mathrm{2}}$(temp)) were calculated as $$\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{temp}\right)=p{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{obs}\right)-p{\mathrm{CO}}_{\mathrm{2}}\left(\stackrel{\mathrm{‾}}{\mathrm{SST}}\right),$$ where $p$CO${}_{\mathrm{2}}$(obs) is the in situ $p$CO${}_{\mathrm{2}}$ and $p$CO${}_{\mathrm{2}}$($\stackrel{\mathrm{‾}}{\mathrm{SST}}$) is the temperature-normalized $p$CO${}_{\mathrm{2}}$. Since changes in $p$CO${}_{\mathrm{2}}$ at a common temperature primarily reflect changes in DIC, the spatial variations in $p$CO${}_{\mathrm{2}}$($\stackrel{\mathrm{‾}}{\mathrm{SST}}$) can be attributed to the combined influences of non-thermal processes that affect DIC (water mass mixing and biological activity). Theoretically, $p$CO${}_{\mathrm{2}}$($\stackrel{\mathrm{‾}}{\mathrm{SST}}$) can be further partitioned into the $p$CO${}_{\mathrm{2}}$ change due to biology and the $p$CO${}_{\mathrm{2}}$ change due to mixing through an analysis of the water mass structure, e.g., an optimum multiparameter (OMP) water mass analysis. Results from its application will be presented in a subsequent study.

Air–sea CO${}_{\mathrm{2}}$ flux estimation

Air–sea CO${}_{\mathrm{2}}$ gas exchange ($F$, mmol C m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ at each sampling location was estimated as follows: $$F=k\cdot K_{\mathrm{0}}\cdot \left(p,{\mathrm{CO}}_{\mathrm{2}},\left(\mathrm{water}\right),-,p,{\mathrm{CO}}_{\mathrm{2}},\left(\mathrm{air}\right)\right),$$ where $k$ (cm h${}^{-\mathrm{1}})$ is the gas transfer velocity of CO${}_{\mathrm{2}}$, $K_{\mathrm{0}}$ (mol kg${}^{-\mathrm{1}}$ atm${}^{-\mathrm{1}})$ is the solubility coefficient of CO${}_{\mathrm{2}}$ at in situ surface-water temperature and salinity (Weiss, 1974), and $p$CO${}_{\mathrm{2}}$(water) and $p$CO${}_{\mathrm{2}}$(air) ($\mathrm{\mu }$atm) are the partial pressures of CO${}_{\mathrm{2}}$ in the water and the air, respectively. The difference between $p$CO${}_{\mathrm{2}}$(water) and $p$CO${}_{\mathrm{2}}$(air) ($\mathrm{\Delta }p$CO${}_{\mathrm{2}})$ determines the direction of gas exchange across the air–sea interface. Positive values of $F$ indicate CO${}_{\mathrm{2}}$ release by the surface water, whereas negative values indicate CO${}_{\mathrm{2}}$ uptake. Conversion factors were applied to express the final $F$ with the aforementioned units.

Atmospheric $p$CO${}_{\mathrm{2}}$ ($p$CO${}_{\mathrm{2}}$(air)) was calculated using the monthly averages of the measured mole fraction of CO${}_{\mathrm{2}}$ in dry air ($x$CO${}_{\mathrm{2}}$, at the greenhouse gas observational station in Fraserdale, Ontario) obtained from the Climate Research Division at Environment and Climate Change Canada. The mean $p$CO${}_{\mathrm{2}}$(air) in the sampling month was computed using the following relationship (Takahashi et al., 2002): $$p{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{air}\right)=x{\mathrm{CO}}_{\mathrm{2}}\cdot \left(P_{\mathrm{b}},-,P_{\mathrm{w}}\right),$$ where $x$CO${}_{\mathrm{2}}$ is in ppm, $P_{\mathrm{b}}$ (atm) is the atmospheric (or barometric) pressure at the sea surface, and $P_{\mathrm{w}}$ (atm) is the equilibrium (or saturation) water vapor pressure at in situ surface-water temperature and salinity (Weiss and Price, 1980). Barometric pressures, 1-month averaged per sampling period, were calculated using the hourly station pressure data from Environment Canada at the following weather observing stations: Québec/Jean Lesage International Airport (upper estuary), Mont-Joli Airport (lower estuary), and Gaspé Airport (GSL). The $P_{\mathrm{b}}$ at station elevation was converted to mean sea level pressure using the formula of Tim Brice and Todd Hall (NOAA's National Weather Service, http://www.weather.gov/epz/wxcalc_stationpressure).

The formulation of the gas transfer velocity, $k$, is the largest source of error in the computation of air–sea CO${}_{\mathrm{2}}$ fluxes (Borges et al., 2004a, b). Properly constraining values of $k$ in estuaries is problematic (Raymond and Cole, 2001) due to their hydrodynamic and geomorphologic complexity compared to the open ocean (Abril et al., 2000). Gas transfer is thought to be regulated by turbulence at the air–water interface (Wanninkhof, 1992). Wind stress plays a key role in the generation of turbulence at the ocean surface through the transfer of momentum to waves and currents (Ho et al., 2011), whereas, in estuarine environments and especially macrotidal estuaries, surface turbulence can be created by interactions of wind forcing, tidal currents, and boundary friction (Zappa et al., 2003, 2007; Borges et al., 2004a, b). In turbid estuaries, surface turbulence can be attenuated by suspended material (Abril et al., 2009). The turbulence generated from bottom stress varies with water depth and tidal velocity (Raymond et al., 2000), and it is important only in shallower estuaries with high current velocities (Cerco, 1989). Raymond and Cole (2001) have shown that wind stress controls turbulence at the air–water interface for all systems with depths greater than 10 m (at depths $<$ 10 m, either wind or bottom stress may dominate).

Several different predictive relationships between wind speed and gas transfer velocity of CO${}_{\mathrm{2}}$ have been proposed based on laboratory and field studies. Here, we estimate the latter from short-term (or steady) wind speed measurements using the equations of Wanninkhof (1992) revised by Wanninkhof (2014) and Raymond and Cole (2001):

$k$ from Wanninkhof (2014), denoted as $k_{\text{W-14}}$: $$k_{\text{W-14}}=\mathrm{0.251}{u}^{\mathrm{2}}{\left(\mathit{\text{Sc}},/,\mathrm{660}\right)}^{-\mathrm{0.5}},$$ and $k$ from Raymond and Cole (2001), denoted as $k_{\text{R}\mathit{\&}\text{C-01}}$: $$k_{\text{R}\mathit{\&}\text{C-01}}=\mathrm{1.91}{e}^{\mathrm{0.35}u}{\left(\mathit{\text{Sc}},/,\mathrm{660}\right)}^{-\mathrm{0.5}},$$ where $u$ is the wind speed (m s${}^{-\mathrm{1}})$ and Sc is the Schmidt number (Sc $=$ $\mathit{\mu }/D$, where $\mathit{\mu }$ is the kinematic viscosity of the water and $D$ is the diffusion coefficient) for CO${}_{\mathrm{2}}$ gas in solution. The Schmidt number for CO${}_{\mathrm{2}}$ in seawater at 20 ${}^{\circ }$C is 660 and was adjusted to Sc $=$ 600 for freshwater. Hourly wind speed data were obtained from Environment Canada at the aforementioned weather observing stations and averaged over the sampling month to obtain short-term wind speeds. The correction to a common Schmidt number was performed using the equations of Wanninkhof (1992) for the temperature dependence of Sc for CO${}_{\mathrm{2}}$ gas in seawater ($S_P=\mathrm{35}$) and freshwater and assuming that $k$ is proportional to Sc${}^{-\mathrm{0.5}}$.

Because of increased turbulence, one would expect $k$ values calculated from estuarine parameterizations to be higher than those predicted from oceanic parameterizations at equivalent wind speeds (Abril et al., 2000). Within the confines of the SLE, estimates of $k$ using the Wanninkhof (2014) relationship ranged from 1.6 to 4.5 cm h${}^{-\mathrm{1}}$, whereas those calculated from Raymond and Cole (2001) were between 3.8 and 8.1 cm h${}^{-\mathrm{1}}$. Hence, we take the air–sea CO${}_{\mathrm{2}}$ flux values calculated with $k_{\text{W-14}}$ to be the theoretical lower limit of gas exchange ($F_{\text{W-14}})$, whereas those computed from $k_{\text{R}\mathit{\&}\text{C-01}}$ represent the upper limit of gas exchange ($F_{\text{R}\mathit{\&}\text{C-01}})$.

In order to estimate the area-averaged CO${}_{\mathrm{2}}$ flux in the SLE, the estuary proper was divided into five segments, with each section containing at least one sampling location. Given that the lower estuary occupies $\sim$ 75 % of the total estuarine surface area and encompasses a fairly wide range of $p$CO${}_{\mathrm{2}}$ values (standard deviation of 119 $\mathrm{\mu }$atm), the SLE was divided into longitudinal sections (Fig. 2) rather than being segmented by salinity. The fluxes in each segment were normalized to the sectional surface area and then summed to obtain a spatially integrated air–sea CO${}_{\mathrm{2}}$ flux ($F_{\text{area-avg}})$ for the whole estuary, as follows (Jiang et al., 2008): $$F_{\text{area-avg}}={\displaystyle \frac{\sum F_i\cdot S_i}{\sum S_i}},$$ where $F_i$ is the average of all the fluxes within segment $i$, and $S_i$ is the surface area of segment $i$. Sectional surface areas were tabulated in MATLAB using the land mask of eastern Canada obtained from Fisheries and Oceans Canada. An area-averaged CO${}_{\mathrm{2}}$ flux was obtained for both the upper and lower limits of gas exchange in the SLE. These two final estimates are assumed to bracket the real areal CO${}_{\mathrm{2}}$ flux.

Conceptual framework for the analysis of variations in biogenic gas concentrations

A comparison of the distribution of biologically reactive dissolved gases, i.e., CO${}_{\mathrm{2}}$ and O${}_{\mathrm{2}}$, can provide useful information about the physical (thermal) and biological processes controlling their concentrations (Richey et al., 1988). Temperature-related gas solubility effects occur in the same direction for CO${}_{\mathrm{2}}$ and O${}_{\mathrm{2}}$, whereas biological production and respiration affect CO${}_{\mathrm{2}}$ and O${}_{\mathrm{2}}$ in opposite directions. Following the approach of Carrillo et al. (2004), the saturation states (or percent saturation) of $p$CO${}_{\mathrm{2}}$ and DO, with respect to the atmosphere, were compared in order to determine the relative importance of temperature effects (heating or cooling) and biological activity (photosynthesis or respiration) in the surface waters at each sampling location. The $p$CO${}_{\mathrm{2}}$ percent saturation ($p$CO${}_{\mathrm{2}}$(% sat)) was calculated as follows: $$p{\mathrm{CO}}_{\mathrm{2}}\left(\mathit{\%},,\text{sat}\right)=\left(p,{\mathrm{CO}}_{\mathrm{2}},\left(\text{water}\right),/,p,{\mathrm{CO}}_{\mathrm{2}},\left(\text{air}\right)\right)\cdot \mathrm{100}.$$ The DO percent saturation (DO(% sat)) was calculated as $$\text{DO}\left(\mathit{\%},,\text{sat}\right)=\left(\text{DO},/,{\text{DO}}^{\ast }\right)\cdot \mathrm{100},$$ where DO${}^{\ast }$ is the equilibrium DO concentration ($\mathrm{\mu }$mol kg${}^{-\mathrm{1}})$ at in situ surface-water temperature and salinity (Benson and Krause, 1984). The relationship between DO(% sat) and $p$CO${}_{\mathrm{2}}$(% sat) is roughly analogous to that of the apparent oxygen utilization (AOU) and excess DIC (eDIC) (Abril et al., 2000). The former is defined as the difference between DO${}^{\ast }$ and DO, whereas the latter is defined as the difference between the observed DIC and a theoretical DIC at atmospheric equilibrium.

Mean, standard deviation, and range of the surface-water temperature ($T$), practical salinity ($S_P)$, dissolved inorganic carbon (DIC), total alkalinity (TAlk), and in situ partial pressure of CO${}_{\mathrm{2}}$ ($p$CO${}_{\mathrm{2}})$ in the St. Lawrence River (near Québec City), USLE (Île d'Orléans to Tadoussac), LSLE (Tadoussac to Pointe-des-Monts), and GSL (Pointe-des-Monts to Cabot Strait) during all sampling months. Numbers in parentheses indicate the ranges observed.

	$T$ (${}^{\circ }$C)	$S_P$	DIC ($\mathrm{\mu }$mol kg${}^{-\mathrm{1}})$	TAlk ($\mathrm{\mu }$mol kg${}^{-\mathrm{1}})$	$p$CO${}_{\mathrm{2}}$ ($\mathrm{\mu }$atm)
River	14.2 $\pm$ 3.9	0.03 $\pm$ 0.05	1242 $\pm$ 132	1204 $\pm$ 99	604 $\pm$ 76
($N=\mathrm{3}$)	(9.8–17.2)	(0–0.09)	(1148–1335)	(1124–1314)	(550–658)
Upper	9.6 $\pm$ 3.6	10.9 $\pm$ 8.0	1514 $\pm$ 242	1492 $\pm$ 272	571 $\pm$ 72
($N=\mathrm{46}$)	(4.2–17.4)	(0–24.5)	(1081–2005)	(969–2030)	(435–765)
Lower	6.2 $\pm$ 2.2	26.2 $\pm$ 2.1	1837 $\pm$ 82	1957 $\pm$ 82	394 $\pm$ 119
($N=\mathrm{60}$)	(2.7–12.6)	(21.2–30.4)	(1634–2005)	(1752–2088)	(139–578)
GSL	8.8 $\pm$ 3.1	30.1 $\pm$ 1.5	1936 $\pm$ 64	2096 $\pm$ 61	352 $\pm$ 80
($N=\mathrm{30}$)	(3.9–13.7)	(25.5–31.5)	(1761–2032)	(1921–2175)	(207–478)

According to the method of Carrillo et al. (2004), data points fall into one of four quadrants on a graph of DO(% sat) versus $p$CO${}_{\mathrm{2}}$(% sat), with the origin at 100 % saturation for both gases. Quadrant I (upper left; supersaturated DO, undersaturated $p$CO${}_{\mathrm{2}})$ suggests net photosynthesis, Quadrant II (upper right; supersaturated DO and $p$CO${}_{\mathrm{2}})$ indicates the effects of heating, Quadrant III (lower right; undersaturated DO, supersaturated $p$CO${}_{\mathrm{2}})$ implies net respiration, and Quadrant IV (lower left; undersaturated DO and $p$CO${}_{\mathrm{2}})$ represents the effects of cooling. Although general patterns become apparent, we urge caution in the interpretation of these results as significant limitations apply. Surface-water CO${}_{\mathrm{2}}$ and O${}_{\mathrm{2}}$ may be acted upon by other forcings such as air–sea gas exchange. The net transfer of CO${}_{\mathrm{2}}$ and O${}_{\mathrm{2}}$ gases occurs across the air–sea interface whenever their partial pressures in the SML differ from those in the atmosphere. Because of their differential gas exchange rates (i.e., O${}_{\mathrm{2}}$ exchanges $\sim$ 19 times faster than CO${}_{\mathrm{2}}$; Peng et al., 1987), CO${}_{\mathrm{2}}$ and O${}_{\mathrm{2}}$ dynamics may be decoupled in surface waters, causing an asymmetry in the observed CO${}_{\mathrm{2}}$ : O${}_{\mathrm{2}}$ relationship (Carrillo et al., 2004).

Results and discussion

Spatial variability in surface-water $p$CO${}_{\mathrm{2}}$

Data were compiled from all 10 cruises to describe the inorganic carbon chemistry in the mixed-layer waters of the St. Lawrence River, USLE, LSLE, and GSL (Table 1). Large spatial variations in surface-water $p$CO${}_{\mathrm{2}}$ were observed within the EGSL system, with values ranging from 139 to 765 $\mathrm{\mu }$atm (452 $\pm$ 134 $\mathrm{\mu }$atm) during the spring and summer sampling periods. Overall, the $p$CO${}_{\mathrm{2}}$ were higher in the USLE (571 $\pm$ 72 $\mathrm{\mu }$atm) than in the LSLE (394 $\pm$ 119 $\mathrm{\mu }$atm) and GSL (352 $\pm$ 80 $\mathrm{\mu }$atm), whereas the atmospheric $p$CO${}_{\mathrm{2}}$ showed less variability, ranging from 372 to 405 $\mathrm{\mu }$atm during the sampling years. As shown in Fig. 4, the USLE was always a CO${}_{\mathrm{2}}$ source (i.e., surface-water $p$CO${}_{\mathrm{2}}$ values were above atmospheric level), while the LSLE and GSL were generally either a CO${}_{\mathrm{2}}$ sink or nearly neutral (i.e., surface-water $p$CO${}_{\mathrm{2}}$ values were below or close to atmospheric level).

Spatial distributions of surface-water $p$CO${}_{\mathrm{2}}$ (circles, diamonds) and practical salinity ($\times$ symbols) in the St. Lawrence River, SLE, and GSL during spring and summer cruises. Horizontal lines show the mean $p$CO${}_{\mathrm{2}}$(air) in the sampling months. The $p$CO${}_{\mathrm{2}}$ data points above atmospheric level are sources of CO${}_{\mathrm{2}}$ to the atmosphere, whereas those below atmospheric level are sinks of atmospheric CO${}_{\mathrm{2}}$.

[Figure omitted. See PDF]

Within the confines of the SLE, the surface-water $p$CO${}_{\mathrm{2}}$ values generally decreased with increasing distance from the head of the estuary (Île d'Orléans) and along the salinity gradient (Figs. 4 and 5). The highest values of $p$CO${}_{\mathrm{2}}$ were observed near the landward limit of the saltwater intrusion in the SLE's upper reaches, in the vicinity of the Cap Tourmente intertidal flats and marshes. This area ($\mathrm{3}\times {\mathrm{10}}^{\mathrm{6}}$ m${}^{\mathrm{2}})$ is located along the core of the estuary's maximum turbidity zone (MTZ) (Lucotte and d'Anglejan, 1986). The lowest surface-water $p$CO${}_{\mathrm{2}}$ values were found downstream of the MTZ in the lower reaches of the SLE near Pointe-des-Monts, where the channel widens into the GSL. Due to favorable environmental conditions (nutrients, light, and stratification), phytoplankton blooms typically occur in late spring or early summer in the LSLE (Zakardjian et al., 2000), with maximal biological production occurring in its downstream portion due to the mixing of cold nutrient-rich waters, which upwell at the head of the Laurentian Channel, with warmer freshwaters flowing in from the north-shore rivers (Savenkoff et al., 1994). Seaward of the estuary–gulf boundary, the $p$CO${}_{\mathrm{2}}$ gradually increased from 207 to 478 $\mathrm{\mu }$atm, coinciding with rising surface-water temperatures ($T=\mathrm{3.9}$ to 13.7 ${}^{\circ }$C).

Spatial distribution of surface-water $p$CO${}_{\mathrm{2}}$ ($\mathrm{\mu }$atm) in the EGSL during all spring and summer cruises. Linear interpolation was used between the sampling locations in the estuary ($+$ symbols). Circles show calculated $p$CO${}_{\mathrm{2}}$ (pH $/$ TAlk), whereas diamonds show $p$CO${}_{\mathrm{2}}$ measured by the underway system (General Oceanics model 8050) aboard the CCGS Amundsen in June 2016 (T. Papakyriakou, University of Manitoba, personal communication, 2017). For neighboring locations sampled in May 2016 aboard the RV Coriolis II, measured and calculated $p$CO${}_{\mathrm{2}}$ differed by, on average, $\sim$ 4.2 %. The mean atmospheric $p$CO${}_{\mathrm{2}}$ during the sampling months ranged from 372 to 405 $\mathrm{\mu }$atm.

[Figure omitted. See PDF]

The spatial variability in surface-water $p$CO${}_{\mathrm{2}}$ due to temperature variations was removed by normalizing the $p$CO${}_{\mathrm{2}}$ data to a common temperature ($T=\mathrm{7.82}$ ${}^{\circ }$C). From a comparison of the in situ and corresponding temperature-normalized $p$CO${}_{\mathrm{2}}$, spatial variations in surface-water temperature lowered or raised the $p$CO${}_{\mathrm{2}}$ by $-$122 to 181 $\mathrm{\mu }$atm within the EGSL system. The maximum (minimum) values of $\mathrm{\Delta }p$CO${}_{\mathrm{2}}$(temp), expressed as a percent change, were 38 % ($-$14 %) in the USLE, 24 % ($-$20 %) in the LSLE, and 29 % ($-$17 %) in the GSL. Temperature normalization, however, removed only a small part of the overall spatial variability in surface-water $p$CO${}_{\mathrm{2}}$ (Fig. 6). Given that the spread of the $p$CO${}_{\mathrm{2}}\left(\stackrel{\mathrm{‾}}{\mathrm{SST}}\right)$ data remained large (153–668 $\mathrm{\mu }$atm; 447 $\pm$ 133 $\mathrm{\mu }$atm), most of the spatial variability in surface-water $p$CO${}_{\mathrm{2}}$ can be explained by nonthermal physical and biological processes that affect DIC concentrations in the surface mixed layer.

Surface-water in situ $p$CO${}_{\mathrm{2}}$, temperature-normalized $p$CO${}_{\mathrm{2}}$, and $\mathrm{\Delta }p$CO${}_{\mathrm{2}}$(temp) averaged over salinity bins of 3. The open circles show the average temperature for each salinity bin. To correct for the increasing or decreasing effect of temperature on surface-water $p$CO${}_{\mathrm{2}}$, the in situ $p$CO${}_{\mathrm{2}}$ values were normalized to the average surface-water temperature of the study area ($T=\mathrm{7.82}$ ${}^{\circ }$C). The $\mathrm{\Delta }$$p$CO${}_{\mathrm{2}}$(temp) values are the thermally induced $p$CO${}_{\mathrm{2}}$ changes due to temperature deviations from $T=\mathrm{7.82}$ ${}^{\circ }$C, whereas variations in temperature-normalized $p$CO${}_{\mathrm{2}}$ are due to water mass mixing and/or biological activity.

[Figure omitted. See PDF]

Mean, standard deviation, and range of $\mathrm{\Delta }$$p$CO${}_{\mathrm{2}}$, $k_{\text{W-14}}$, $k_{\text{R}\mathit{\&}\text{C-01}}$, $F_{\text{W-14}}$, and $F_{\text{R}\mathit{\&}\text{C-01}}$ in the surface waters of the St. Lawrence River, USLE, LSLE, and GSL during all sampling months. $k_{\text{W-14}}$ is the gas transfer velocity given by the Wanninkhof (2014) relationship, whereas $k_{\text{R}\mathit{\&}\text{C-01}}$ is that given by Raymond and Cole (2001). Values of $F_{\text{W-14}}$ are taken to be the theoretical lower limit of air–sea gas exchange, whereas values of $F_{\text{R}\mathit{\&}\text{C-01}}$ are the upper limit. Numbers in parentheses indicate the ranges observed; the extreme $F$ data points are shown in bold.

	$\mathrm{\Delta }p$CO${}_{\mathrm{2}}$	$k_{\text{W-14}}$	$k_{\text{R}\mathit{\&}\text{C-01}}$	$F_{\text{W-14}}$	$F_{\text{R}\mathit{\&}\text{C-01}}$
	($\mathrm{\mu }$atm)	(cm h${}^{-\mathrm{1}})$	(cm h${}^{-\mathrm{1}})$	(mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$	(mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$
River	217 $\pm$ 99	3.0 $\pm$ 1.4	6.1 $\pm$ 2.0	5.8 $\pm$ 3.2	12.7 $\pm$ 6.4
($N=\mathrm{3}$)	(147/287)	(1.9/4.5)	(4.3/8.2)	(3.5/8.0)	(8.2/17.3)
Upper	184 $\pm$ 72	2.8 $\pm$ 0.8	5.6 $\pm$ 1.1	6.1 $\pm$ 3.0	12.3 $\pm$ 5.4
($N=\mathrm{46}$)	(43/386)	(1.6/4.5)	(3.8/8.1)	(2.0/14.7)	(3.6/28.4)
Lower	9.2 $\pm$ 116	3.2 $\pm$ 0.4	5.9 $\pm$ 0.6	0.6 $\pm$ 4.9	1.0 $\pm$ 9.0
($N=\mathrm{60}$)	($-$266/188)	(2.0/3.8)	(4.3/6.9)	($-$12.1/8.3)	($-$21.9/15.1)
GSL	$-$31.0 $\pm$ 90	1.2 $\pm$ 0.3	3.4 $\pm$ 0.3	$-$0.8 $\pm$ 1.5	$-$1.7 $\pm$ 3.9
($N=\mathrm{30}$)	($-$178/107)	(0.8/1.7)	(2.8/4.1)	($-$3.6/1.1)	($-$8.4/3.6)

Air–sea CO${}_{\mathrm{2}}$ flux and spatial integration

Large spatial variations in the air–sea CO${}_{\mathrm{2}}$ flux were observed within the EGSL system during the spring and summer sampling months, with fluxes ranging from $-$21.9 to 28.4 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ (Fig. 7). Values of $F$ were always positive in the USLE (2.0 to 28.4 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ and either negative or positive in the LSLE ($-$21.9 to 15.1 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ and GSL ($-$8.4 to 3.6 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$. As expected, $F_{\text{R}\mathit{\&}\text{C-01}}$ (estuarine parameterization of $k$) were larger than $F_{\text{W-14}}$ (oceanic parameterization of $k$) due to the inherently greater surface turbulence in estuarine systems. The average difference between CO${}_{\mathrm{2}}$ fluxes calculated using the two formulations of the gas transfer velocity (Eqs. 4 and 5) was 71.7 %. Details of the $k$ and $F$ values given by each parameterization of $k$ are shown in Table 2. Irrespective of the parameterization, the calculated CO${}_{\mathrm{2}}$ fluxes were more positive in the USLE (9.2 $\pm$ 5.3 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ than in the LSLE (0.8 $\pm$ 7.2 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ and GSL ($-$1.2 $\pm$ 3.0 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$.

Sectional and area-averaged air–sea CO${}_{\mathrm{2}}$ fluxes (mmol C m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ in the SLE during all sampling months. To obtain the area-averaged CO${}_{\mathrm{2}}$ flux, the SLE was divided into five segments at equal intervals. The first row of the table shows the surface area (km${}^{\mathrm{2}})$ of each segment. The flux data in each segment were numerically averaged to obtain sectional fluxes, which were then area weighted and summed to obtain the spatially integrated whole-estuary flux (in bold). The $F_{\text{W-14}}$ and $F_{\text{R}\mathit{\&}\text{C-01}}$ data provide the lower and upper estimates, respectively.

	Seg 1	Seg 2	Seg 3	Seg 4	Seg 5	Whole
	($N=\mathrm{17}$)	($N=\mathrm{23}$)	($N=\mathrm{21}$)	($N=\mathrm{17}$)	($N=\mathrm{8}$)	estuary
Surface area (km${}^{\mathrm{2}})$	1098	2250	2726	3404	3303	12 781
$F_{\text{W-14}}$ (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$	7.2	5.4	4.3	$-$1.8	$-$4.0	0.98
$F_{\text{R}\mathit{\&}\text{C-01}}$ (mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$	14.5	11.0	7.9	$-$3.5	$-$7.4	2.02

Box plot showing the variability in air–sea CO${}_{\mathrm{2}}$ fluxes in the four principal subdivisions of the study area (St. Lawrence River, USLE, LSLE, and GSL). The box spans the interquartile range (25–75 percentiles), the dashed line is the median, and the whiskers extend to the extreme data points. The $F_{\text{W-14}}$ data were combined with the $F_{\text{R}\mathit{\&}\text{C-01}}$ data from all spring and summer sampling months to depict the upper and lower limits of gas exchange.

[Figure omitted. See PDF]

The SLE was divided into five segments to obtain an area-averaged CO${}_{\mathrm{2}}$ flux for the whole estuary. The data used to calculate the $F_{\text{area-avg}}$ are listed in Table 3. Overall, the SLE served as a weak source of CO${}_{\mathrm{2}}$ to the atmosphere at the time of sampling, with an area-averaged degassing flux of 0.98 to 2.02 mmol C m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ (0.36 to 0.74 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$ during the late spring and early summer. This efflux compares favorably with that of the Delaware Estuary (2.4 $\pm$ 4.8 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$; Joesoef et al., 2015), another large estuarine system with a long water residence time, but is significantly lower than estimates in the marine-dominated Sapelo Sound and Doboy Sound estuaries (10.5 to 10.7 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$; Jiang et al., 2008). From a compilation of 165 estuaries worldwide, almost all systems, with the exception of those in the Arctic ($-$1.1 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$, serve as sources of CO${}_{\mathrm{2}}$ to the atmosphere (Chen et al., 2013). Chen et al. (2013) concluded that the world's upper estuaries ($S_P$ $<$ 2) are strong sources (39.0 $\pm$ 55.7 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$, mid-estuaries (2 $<$ $S_P$ $<$ 25) are moderate sources (17.5 $\pm$ 34.2 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$, and lower estuaries ($S_P$ $>$ 25) are weak sources (8.4 $\pm$ 14.3 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$. Predictably, with its maritime region occupying almost three-fourths of the total surface area, the SLE behaves like an outer estuary with only small CO${}_{\mathrm{2}}$ evasion. The lack of temporal coverage of surface-water $p$CO${}_{\mathrm{2}}$ data, however, prevents us from reliably synthesizing an annual air–sea CO${}_{\mathrm{2}}$ flux.

Major drivers of estuarine $p$CO${}_{\mathrm{2}}$ variability

The $p$CO${}_{\mathrm{2}}$ in the surface mixed layer is a function of its $T$, $S_P$, DIC, and TAlk, as described by the following relationship (Takahashi et al., 1993): $$\begin{array}{cc}& {\displaystyle }\text{d}p{\mathrm{CO}}_{\mathrm{2}}& {\displaystyle }=\left(\partial ,p,{\mathrm{CO}}_{\mathrm{2}},/,\partial ,T\right)\mathrm{d}T+\left(\partial ,p,{\mathrm{CO}}_{\mathrm{2}},/,\partial ,S_P\right)\text{d}{S}_P{\displaystyle }& {\displaystyle }+\left(\partial ,p,{\mathrm{CO}}_{\mathrm{2}},/,\partial ,\text{DIC}\right)\text{dDIC}+\left(\partial ,p,{\mathrm{CO}}_{\mathrm{2}},/,\partial ,\text{TAlk}\right)\text{dTAlk}.\end{array}$$ Through changes in $T$, $S_P$, DIC, and TAlk, variations in surface-water $p$CO${}_{\mathrm{2}}$ are mainly controlled by dynamic processes (water mass mixing), thermodynamic processes (temperature and salinity changes), air–sea gas exchange, and biological processes (photosynthesis and respiration) (Poisson et al., 1993). Among these, the effects of temperature and DIC, i.e., the addition or removal of DIC through biological activity and mixing processes, are generally the most important drivers of estuarine $p$CO${}_{\mathrm{2}}$ variability. In the absence of a significant source or sink of TAlk (e.g., calcium carbonate formation and/or dissolution, anaerobic organic matter decomposition) and after normalizing the data to a common temperature, changes in DIC determine the buffer capacity (DIC : TAlk ratio) and the $p$CO${}_{\mathrm{2}}$ saturation of the water. Whereas the physically and biologically induced changes of DIC $/$ $p$CO${}_{\mathrm{2}}$ will be quantified in a future study, using a modified OMP water mass analysis, here we evaluate the relative importance of thermal and biological processes in controlling the spatial distribution of $p$CO${}_{\mathrm{2}}$ in the EGSL.

(a) Comparison of the saturation states (or percent saturation) of $p$CO${}_{\mathrm{2}}$ and DO in the mixed-layer waters of the St. Lawrence River, SLE, and GSL. Dashed lines delineate the 100 % saturation levels for both gases. Surface-water samples (markers) fall into one of four quadrants representing the dominant controls on CO${}_{\mathrm{2}}$ $/$ O${}_{\mathrm{2}}$ dynamics. Quadrants I and III indicate the effects of photosynthesis and respiration, whereas Quadrants II and IV indicate heating and cooling. Panel (b) presents a comparison of the apparent oxygen utilization (AOU) and excess DIC (eDIC). Respiration and remineralization processes are reflected in positive values of AOU and eDIC, whereas the effects of photosynthesis are reflected in negative values.

[Figure omitted. See PDF]

Spatial distributions of surface-water $p$CO${}_{\mathrm{2}}$ (circles) and temperature (triangles) in the St. Lawrence River, SLE, and GSL during spring and summer cruises. Temperatures ranged from 4.2 to 17.4 ${}^{\circ }$C (generally decreasing) in the USLE, 2.7 to 12.6 ${}^{\circ }$C (generally increasing) in the LSLE, and 3.9 to 13.7 ${}^{\circ }$C (generally increasing) in the GSL. The horizontal line shows the mean atmospheric $p$CO${}_{\mathrm{2}}$, $p$CO${}_{\mathrm{2}}$(air), during all sampling months. The dashed line is the smoothed $p$CO${}_{\mathrm{2}}$(water) data using a moving average filter with a span of 50 % of the total number of data points, whereas the red line is the smoothed temperature data.

[Figure omitted. See PDF]

To disentangle the biological and temperature effects on the spatial variability in $p$CO${}_{\mathrm{2}}$, the DO(% sat) values were plotted against the $p$CO${}_{\mathrm{2}}$(% sat) values, with the origin at 100 % saturation for both gases. This simple approach uses the four possible combinations of $p$CO${}_{\mathrm{2}}$(% sat) $/$ DO(% sat) as integrated measures of thermally and biologically induced changes. As shown in Fig. 8a, microbial respiration was the major driver of $p$CO${}_{\mathrm{2}}$ variability in the USLE, whereas photosynthesis and temperature were the dominant controls in the LSLE and GSL. Similar results were found from a comparison of eDIC and AOU (Fig. 8b). In the strongly stratified lower estuary, as well as near the estuary–gulf boundary, the biological drawdown of CO${}_{\mathrm{2}}$ counteracted the decrease in CO${}_{\mathrm{2}}$ gas solubility due to increasing temperature (Fig. 9). Its waters were mostly undersaturated with CO${}_{\mathrm{2}}$ with respect to the atmosphere (values of $p$CO${}_{\mathrm{2}}$ were below atmospheric level) despite a general trend of surface-water warming ($T=\mathrm{2.7}$ to 12.6 ${}^{\circ }$C). This pattern is consistent with the finding that, in spring and summer, the increasing effect of warming on $p$CO${}_{\mathrm{2}}$ is counteracted by the photosynthetic utilization of CO${}_{\mathrm{2}}$, particularly in a strongly stratified shallow mixed layer (Takahashi et al., 1993). Whereas direct measurements of chlorophyll-$a$ concentrations were not carried out during the research cruises, a fluorescence sensor was mounted on the CTD probe. As shown in Fig. 10, maximum fluorescence values, as well as high values of transmission (percentage of light transmission approaching 100 %), were observed in the eastern lower estuary and the western GSL, where the system appears to shift from net heterotrophy to net autotrophy. Farther into the GSL (near Anticosti Island), the temperature dependence of $p$CO${}_{\mathrm{2}}$ exerted a stronger influence, causing values of surface-water $p$CO${}_{\mathrm{2}}$ to increase concomitantly with temperature (Fig. 9).

Spatial distributions of maximum fluorescence values (circles) and mean transmission values ($\times$ symbols) in the euphotic zone of the St. Lawrence River, SLE, and GSL during the May 2016 cruise. Fluorescence is a primary production proxy, whereas transmission is an excellent proxy for turbidity (low transmission values are due to light absorption by suspended particulate matter and/or colored dissolved organic matter).

[Figure omitted. See PDF]

Conclusions

Because of its large physical dimensions and unimpeded connection to the Atlantic Ocean, the St. Lawrence Estuary encompasses both a river-dominated inner estuary, where physical mixing and abiotic processes dominate, and a marine-dominated outer estuary, where biological cycling and oceanic processes prevail. The physical and biogeochemical processes of these contrasting environments are reflected in the spatial distribution of surface-water $p$CO${}_{\mathrm{2}}$ (139–765 $\mathrm{\mu }$atm). The shallow partially mixed upper estuary, with a turbidity maximum controlled by tide- and wind-induced turbulence, was, during our sampling period, a net source of CO${}_{\mathrm{2}}$ to the atmosphere due to microbial respiration (low biological productivity), whereas the deep stratified lower estuary, with its stable, summertime three-layer vertical structure, was generally a net sink of atmospheric CO${}_{\mathrm{2}}$ due to the enhanced biological drawdown of $p$CO${}_{\mathrm{2}}$ (light availability, nutrient supply, and strong stratification).

Overall, the large subarctic St. Lawrence Estuary was a weak source of CO${}_{\mathrm{2}}$ to the atmosphere, with an area-averaged CO${}_{\mathrm{2}}$ degassing flux of 0.98 to 2.02 mmol C m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}}$ (0.36 to 0.74 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$. This efflux is somewhat smaller than the numerically averaged CO${}_{\mathrm{2}}$ flux per unit area (2.19 mol C m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}})$ reported for North American estuaries by Chen et al. (2013), highlighting their relatively small contribution ($\sim$ 12 %) to global estuarine CO${}_{\mathrm{2}}$ emissions. The pronounced shift in source and sink dynamics in the St. Lawrence Estuary, between its river-dominated (9.2 $\pm$ 5.3 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ and marine-dominated (0.8 $\pm$ 7.2 mmol m${}^{-\mathrm{2}}$ d${}^{-\mathrm{1}})$ regions, is consistent with the conclusions of the comparative study carried out by Jiang et al. (2008) that revealed large differences in CO${}_{\mathrm{2}}$ degassing between riverine (inner) and maritime (outer) estuaries. Given the limited research on CO${}_{\mathrm{2}}$ dynamics in large estuaries and bay systems, which cover approximately one-half of the estuarine surface area on the US east coast, as well as the large uncertainties in the indirect measurement of $p$CO${}_{\mathrm{2}}$ (carbonic acid dissociation constants and organic alkalinity contribution), current global-scale estimates of estuarine CO${}_{\mathrm{2}}$ degassing may be overestimated. To better constrain the role of large estuaries and bays in the coastal ocean carbon cycle, more extensive spatial and temporal coverage of direct $p$CO${}_{\mathrm{2}}$ measurements across estuary types is needed.

Data presented in this paper (Figs. 4 and 8) are available upon request from one of the authors ([email protected]).

AD and AM conceived the project. AM acquired and processed the data prior to 2016. AD conducted the data analysis and wrote the first draft of the paper and AM provided editorial and scientific recommendations.

The authors declare that they have no conflict of interest.

Acknowledgements

We wish to thank the captains and crews of the RV Coriolis II for their unwavering help over the years. We also wish to acknowledge Gilles Desmeules and Michel Rousseau for their dedicated electronic and field sampling support as well as Constance Guignard for her help in cruise preparation and field data acquisition. Most of the data acquisition was carried out opportunistically on research cruises funded by grants to Alfonso Mucci or Canadian colleagues by the Natural Sciences and Engineering Research Council of Canada (NSERC), whereas the work was funded by a Regroupement Stratégique grant from the Fonds Québécois de Recherche Nature et Technologies (FQRNT) to GEOTOP as well as NSERC Discovery and MEOPAR grants to Alfonso Mucci. Ashley Dinauer wishes to thank the Department of Earth and Planetary Sciences at McGill University for financial support in the form of scholarships and assistantships. Edited by: Gwenaël Abril Reviewed by: Chen-Tung Arthur Chen, Helmuth Thomas, and one anonymous referee