Spatial variations in CO2 fluxes in the Saguenay Fjord (Quebec, Canada) and results of a water mixing model

The Saguenay Fjord is a major tributary of the St. Lawrence Estuary and is strongly stratified. A 6–8 m wedge of brackish water typically overlies up to 270 m of seawater. Relative to the St. Lawrence River, the surface waters of the Saguenay Fjord are less alkaline and host higher dissolved organic carbon (DOC) concentrations. In view of the latter, surface waters of the fjord are expected to be a net source of CO2 to the atmosphere, as they partly originate from the flushing of organic-rich soil porewaters. Nonetheless, the CO2 dynamics in the fjord are modulated with the rising tide by the intrusion, at the surface, of brackish water from the Upper St. Lawrence Estuary, as well as an overflow of mixed seawater over the shallow sill from the Lower St. Lawrence Estuary. Using geochemical and isotopic tracers, in combination with an optimization multiparameter algorithm (OMP), we determined the relative contribution of known source waters to the water column in the Saguenay Fjord, including waters that originate from the Lower St. Lawrence Estuary and replenish the fjord’s deep basins. These results, when included in a conservative mixing model and compared to field measurements, serve to identify the dominant factors, other than physical mixing, such as biological activity (photosynthesis, respiration) and gas exchange at the air–water interface, that impact the water properties (e.g., pH, pCO2) of the fjord. Results indicate that the fjord’s surface waters are a net source of CO2 to the atmosphere during periods of high freshwater discharge (e.g., spring freshet), whereas they serve as a net sink of atmospheric CO2 when their practical salinity exceeds ∼ 5–10.


Introduction
Anthropogenic emissions of carbon dioxide (CO 2 ) have recently propelled atmospheric CO 2 concentrations above the 410 ppm mark, the highest concentration recorded in the past 3 million years (Willeit et al., 2019). The oceans, the largest CO 2 reservoir on Earth, have taken up ca. 30 % of the anthropogenic CO 2 emitted to the atmosphere since the beginning of the industrial era (Feely et al., 2004;Brewer and Peltzer, 2009;Doney et al., 2009;Orr, 2011;Friedlingstein et al., 2019), mitigating the impact of this greenhouse gas on global warming . On the other hand, the uptake of CO 2 by the oceans has led to modifications of the seawater carbonate chemistry and a decline in the average surface ocean pH by ∼ 0.1 units since pre-industrial times, a phenomenon dubbed ocean acidification (Caldeira and Wickett, 2005). According to the Intergovernmental Panel on Climate Change (IPCC) "business as usual" emissions scenario IS92a and general circulation models, atmospheric CO 2 levels may reach 800 ppm by 2100, lowering the pH of the surface oceans by an additional 0.3-0.4 units, a rate that is unprecedented in the geological record (Caldeira and Wickett, 2005;Hönisch et al., 2012;Rhein et al., 2013). The growing concern about the impacts of anthropogenic CO 2 emissions on climate, as well as marine and terrestrial ecosystems, calls for a meticulous quantification of organic and inorganic carbon fluxes, especially in coastal environments, including fjords, a major but poorly quantified component of the global carbon cycle and budget (Bauer et al., 2013;Najjar et al., 2018). Meaningful predictions of the effects of climate change on future fluxes are intricate given the very large uncertainty associated with present-day air-sea CO 2 flux estimates in coastal waters, including rivers, estuaries, tidal wetlands and the continental shelf (Bauer et al., 2013;Najjar et al., 2018). The coastal ocean occupies only ∼ 7 % of the global ocean surface area but plays a major role in biogeochemical cycles because it (1) receives massive inputs of terrestrial organic matter and nutrients through continental runoff and groundwater discharge; (2) exchanges matter and energy with the open ocean; and (3) is one of the most geochemically and biologically active areas of the biosphere, accounting for significant fractions of marine primary production (∼ 14 % to 30 %), organic matter burial (∼ 80 %), sedimentary mineralization (∼ 90 %), and calcium carbonate deposition (∼ 50 %) (Gattuso et al., 1998).
Although the carbon cycle of the coastal ocean is acknowledged to be a major component of the global carbon cycle and budget, accurate quantification of organic and inorganic carbon cycling and fluxes in the coastal ocean -where land, ocean and atmosphere interact -remains challenging (Bauer et al., 2013;Najjar et al., 2018). Constraining the exchanges and fates of different forms of carbon along the land-ocean continuum is so far incomplete, owing to limited data coverage and large physical and biogeochemical variability within and between coastal subsystems (e.g., hydrological and geomorphological differences, differences in the magnitude and stoichiometry of organic matter inputs). Hence, owing to limited data coverage and suspicious upscaling due to the large physical and biogeochemical variability within and between coastal subsystems, there remains a debate as to whether coastal waters are net sources or sinks of atmospheric CO 2 . Recent compilations of worldwide CO 2 partial pressure (pCO 2 ) measurements indicate that most open shelves in temperate and high latitudes are sinks of atmospheric CO 2 , whereas low-latitude shelves and most estuaries are sources (Chen and Borges, 2009;Cai, 2011;Chen et al., 2013). As noted by Bauer et al. (2013), estuaries are transitional aquatic environments that can be riverine or marine dominated and, thus, they typically display strong gradients in biogeochemical properties and processes as they flow seaward. Chen et al. (2013) reported that the strength of estuarine sources typically decreases with increasing salinity. However, marsh-dominated estuaries, in which active microbial decomposition of organic matter occurs in the intertidal zone, are strong sources of CO 2 (Cai, 2011).
High-latitude waters such as the Arctic Ocean have recently been the foci of much research, while coastal, seasonally ice-covered aquatic environments, such as the Saguenay Fjord, that display comparable inter-annual and climatic seaice cover variabilities but are much more accessible, have been neglected (Bourgault et al., 2012). Characteristics of Arctic coastal ecosystems are found in the Saguenay Fjord, including the presence of many species of plankton, fish, birds and marine mammals, as well as important freshwater inputs and the presence of seasonal ice cover (Bourgault et al., 2012). Fjords stand amongst the most productive ecosystems on the planet, while they have a yet unexplored role in regional and global carbon cycles as part of the estuarine family (Juul-Pedersen et al., 2015). They are crucial hotspots for organic carbon (mostly terrestrial) burial and account for nearly 11 % of the annual organic carbon burial flux in marine sediments, while covering only 0.12 % of oceans' surface (Rysgaard et al., 2012;Smith et al., 2015). In other words, organic carbon burial rates in fjords are 100 times faster than the average rate in the global ocean. Rates of organic carbon burial provide insights into the mechanism that controls atmospheric O 2 and CO 2 concentrations over geological timescales (Smith et al., 2015).
This study presents (1) the relative contribution of known source waters to the water column in the Saguenay Fjord, estimated from the solution of an optimization multiparameter algorithm (OMP) using geochemical and isotopic tracers, and (2) results of a conservative mixing model, based on results of the OMP analysis, and from which theoretical surface water pCO 2 values are derived and then compared to field measurements. The latter comparison serves to identify the dominant factors, other than physical mixing (i.e., biological activity, gas exchange), that impact the CO 2 fluxes at the air-sea interface and modulate their direction and intensity throughout the fjord (i.e., whether it is a source or a sink of CO 2 to the atmosphere).

Study site characteristics
Located in the subarctic region of Quebec, eastern Canada, the Saguenay Fjord is up to 275 m deep, 110 km long and has an average width of 2 km, with a 1.1 km wide mouth where it connects to the head of the Lower St. Lawrence Estuary (Fig. 1a). The fjord's bathymetry includes three basins bound by three sills (Fig. 1b). The first one, at a depth of ∼ 20 m, is located at its mouth near Tadoussac and controls the overall dynamics of the fjord. The second is located 18 km further upstream and sits at a depth of 60 m, while the third one is found another 32 km further upstream and rises to a depth of 115 m. The fjord's drainage basin is 78 000 km 2 and is part of the greater St. Lawrence drainage basin (Smith and Walton, 1980), forming a hydrographic system, along with the Great Lakes, of more than 1.36 million km 2 .
Tributaries to the Saguenay Fjord include the Saguenay, Éternité and Sainte-Marguerite rivers (Fig. 1a). The Saguenay River is the main outlet from Saint-Jean Lake and flows into the north arm of the fjord near St. Fulgence (Fig. 1a) with a mean freshwater discharge of ∼ 1200 m 3 s −1 (Bélanger, 2003). Two other local, minor tributaries, the Rivière-à-Mars (95 km long, mean discharge ∼ 8 m 3 s −1 ) and the Rivière des Ha! Ha! (35 km long, mean discharge ∼ 15 m 3 s −1 ) discharge into the Baie des Ha! Ha!, a distinct feature of the Saguenay Fjord (Fig. 1a). Finally, the fjord receives denser marine waters from the St. Lawrence Estuary, filling the bottom of the three basins, as these waters episodically overflow the entrance sill (Therriault and Lacroix, 1975;Stacey and Gratton, 2001;Bélanger, 2003;Belzile et al., 2016). According to Seibert et al. (1979), the tidal amplitude at the mouth of the fjord near Tadoussac averages 4.0 m and increases slightly toward the head of the fjord (4.3 m near Port Alfred). Spring tides may reach an amplitude of 6 m.
The overflow and the intrusion of marine waters from the St. Lawrence Estuary generate a sharp halocline, leading to a simplified two-layer stratification in the fjord (Fig. 1b). The tidally modulated intrusion of marine waters from the St. Lawrence Estuary into the Saguenay Fjord, as well as the outflow of the fjord into the estuary, have a major influence on the water column stratification and circulation in the Saguenay Fjord and at its mouth (Belzile et al., 2016;Mucci et al., 2017). In other words, the properties of the uppermost 100 m of the water column in the adjacent estuary are critical in determining the water stratification in the Saguenay Fjord, since salinity and temperature control the density of waters that spill over the sill and fill the fjord's deep basins (Belzile et al., 2016).
During most of the ice-free season, the St. Lawrence Estuary is characterized by three distinct layers: (1) a relatively warm and salty bottom layer (LSLE, 4 • C < T < 6 • C, 34 < S P < 34.6, where T stands for temperature and S P refers to practical salinity) that originates from mixing on the continental shelf of Northwestern Atlantic Current and Labrador Current waters, (2) a cold intermediate layer (CIL, 30-150 m deep; −1 • C < T < 2 • C, 31.5 < S P < 33) that forms in the Gulf of St. Lawrence in the winter and flows landward, and (3) a warm brackish surface layer (0-30 m deep, −0.6 • C < T < 12 • C, 25 < S P < 32) that results from the mixture of freshwater from various tributaries (mostly the St. Lawrence and Saguenay rivers but also north-shore rivers such as the Betsiamites, Romaine and Manicouagan) and seawater and flows seaward to ultimately form the Gaspé Current (Dickie and Trites, 1983;El-Sabh and Silverberg, 1990;Gilbert and Pettigrew, 1997). Seasonal variations greatly affect the properties of the surface layer, which merges with the intermediate layer during winter, as temperature and salinity change with atmospheric, and buoyancy forcing and the contribution from tributaries decreases during winter months (Galbraith, 2006).
Likewise, the Saguenay Fjord is characterized by a strongly stratified water column that includes at least two water masses: (1) a warm, shallow layer, the Saguenay Shallow Water (SSW; 0 • C < T < 16.8 • C, 0.2 < S P < 26.9), that lies above (2) the Saguenay Deep Water (SDW; 0.9 • C < T < 4.0 • C, 27.3 < S P < 29.8). The SDW most likely forms from a mixture of surface fjord water, St. Lawrence River waters and the St. Lawrence Estuary cold intermediate layer (CIL), when the latter spills over the entrance sill at the mouth of the fjord (Bourgault et al., 2012;Belzile et al., 2016). Nonetheless, our study shows that, because the Saguenay Fjord is a relatively deep fjord with multiple sills, the vertical structure of the water column is far more complex than described above.

Water column sampling
The data presented in this paper were gathered on five cruises, between the years 2014 and 2018 aboard the R/V Coriolis II in late spring (May 2016 and May 2018) and early summer (June 2017), as well as early and late fall (September 2014 and November 2017). Sampling of the water column was carried out with a rosette system along the central axis of the Saguenay Fjord, between St. Fulgence and the mouth of the fjord, including the Baie des Ha! Ha!. Stations in the St. Lawrence Estuary, near the mouth of the fjord, were also sampled. The sampling locations are identified in Fig. 1a. The surface water of the Saguenay River was sampled, with a rope and bucket in 2013 and 2017, from the Dubuc Bridge that joins Chicoutimi and Chicoutimi-Nord, to determine the chemical characteristics of the freshwater Saguenay River end-member.
The rosette system (12 × 12 L Niskin bottles) was equipped with a Seabird 911Plus conductivity-temperaturedepth (CTD) probe, a Seabird ® SEB-43 oxygen probe, a WETLabs ® C-Star transmissometer and a Seapoint ® fluorometer. The Niskin bottles were closed at discrete depths as the rosette was raised from the bottom, typically at the surface (2-3 m), 25, 50, 75, 100 and at 50 m intervals to the bottom (or within 10 m of the bottom). Samples were taken directly from the bottles for dissolved oxygen (DO), pH NBS and/or pH T , total alkalinity (TA), dissolved inorganic carbon (DIC), dissolved silicate (DSi), practical salinity (S P ), and the stable oxygen isotopic composition of the water (δ 18 O water ). Water samples destined for pH measurements were transferred to 125 mL plastic bottles without headspace, whereas TA and TA/DIC samples were stored in, respectively, 250 and 500 mL glass bottles. TA and TA/DIC samples were poisoned with a few crystals of mercuric chloride (HgCl 2 ), and bottles were sealed using a ground-glass stopper and Apiezon ® Type-M high-vacuum grease. δ 18 O water and S P samples were stored in 13 mL plastic screw-cap test tubes.
Direct measurements of surface water (∼ 2 m) pCO 2 were carried out using a CO 2 -Pro CV (Pro-Oceanus, Bridgewater, NS) probe in May 2018. The CO 2 -Pro CV probe operates through rapid diffusion of gases through a supported semipermeable membrane to a thermostated cell in which the CO 2 mole fraction is quantified by a nondispersive infrared detector (NDIR) that was factory calibrated using standard trace gas mixtures. The instrument was operated in continuous mode, with measurements taken nearly every 7 s. Stable pCO 2 values were achieved after a 15 min equilibration period and averaged over the next 20 min. Relative standard deviations over this period were typically on the order of 0.2 to 6 % but were on the order of 0.1 % in a stable water mass at 220 m depth, implying that deviations recorded at the surface likely reflected natural variations over the period of sampling as the ship drifted with the current. The manufacturer claims a 1 % accuracy, but the performance of the instrument may be even better (Hunt et al., 2017).
Total freshwater discharge data of the Saguenay River were provided by Rio Tinto Alcan (a multinational aluminum smelter and producer that manages its own hydroelectric dam on the Saguenay River) from their bank stabilization program. Data for the relevant sampling days in September 2014, May 2016, June 2017, November 2017 and May 2018 were taken from the Shipshaw and Chute-à-Caron monitoring stations.

Analytical procedures
T and S P were determined in situ using the CTD probe. The conductivity probe was calibrated by the manufacturer over the winter prior to the cruises. In addition, the S P of surface waters was determined by potentiometric argentometric titration at McGill University and calibration of the AgNO 3 titrant with IAPSO standard seawater. The reproducibility of these measurements is typically better than ±0.5 %.
pH T was determined spectrophotometrically on board, on the total hydrogen ion concentration scale for saline waters (S P > 5), using phenol red and purified m-cresol purple as indicators and a Hewlett-Packard UV-visible diode array spectrophotometer (HP-8453A) with a 5 cm quartz cell, after thermal equilibration of the sample in a constant temperature bath at 25 • C ±0.1. The salinity-dependence of the dissociation constants and molar absorptivities of the indicators were taken from Robert-Baldo et al. (1985) for phenol red and from Clayton and Byrne (1993) for m-cresol purple. The salinity-dependence of the phenol red indicator dissociation constant and molar absorptivities was extended (from S P = 5 to 35; Bellis, 2002) to encompass the range of salinities encountered in this study, but computed pH T values from the revised fit were not significantly different from those obtained with the relationship provided by Robert-Baldo et al. (1985). Results computed from these parameters yielded values that were more similar to each other as well as to potentiometric glass electrode measurements than the revised equation for the purified m-cresol purple provided by Douglas and Byrne (2017). The pH of low-salinity waters (S P < 5) was determined potentiometrically on board at 25 • C, on the NIST (formerly NBS) scale (pH NBS ), using a Radiometer Analytical ® (GK2401C) combination glass electrode connected to a Radiometer Analytical ® pH/millivoltmeter (PHM84). A calibration of the electrode was completed prior to and after each measurement, using three NISTtraceable buffer solutions: pH 4.00, pH 7.00 and pH 10.00 at 25 • C. The Nernstian slope was then obtained from the least-squares fit of the electrode response to the NIST buffer values. For waters with S P of between 5 and 35, pH NBS was converted to pH T according to the electrode response to TRIS (tris(hydroxymethyl)aminomethane) buffer solutions prepared at S P = 5, 15, 25 and 35 and for which the pH T was assigned at 25 • C (Millero, 1986). Reproducibility of pH measurements based on replicate analyses of the same sample or at least two of the three methods used was typically better than ±0.005.
Dissolved oxygen (DO) concentrations were determined on board by Winkler titration on distinct water samples recovered directly from the Niskin bottles, following the method described by Grasshoff et al. (1999). The relative standard deviation, based on replicate analyses of samples recovered from the same Niskin bottle, was 0.5 %. These measurements served to calibrate the SBE-43 oxygen probe mounted on the rosette sampler.
The stable oxygen isotopic composition of the water samples (δ 18 O water ) was determined using the CO 2 equilibration method of Epstein and Mayeda (1953). Aliquots (200 µL) of the water samples and three laboratory internal reference waters were transferred into 3 mL vials stoppered with a septum cap. The vials were then placed in a heated rack maintained at 40 • C. Commercially available 99.998 % pure CO 2 gas (Research Grade) was introduced in all the vials using a Micromass AquaPrep and allowed to equilibrate for 7 h. The headspace CO 2 was then sampled by the Micromass AquaPrep, dried on a −80 • C water trap and analyzed on a Micromass Isoprime universal triple collector isotope ratio mass spectrometer in dual inlet mode at the GEOTOP-UQAM Stable Isotope Laboratory. Data were normalized against the three internal reference waters, themselves calibrated against Vienna Standard Mean Ocean Water (V-SMOW) and Vienna Standard Light Arctic Precipitation (V-SLAP). The results are reported on the δ scale in ‰ relative to V-SMOW: Based on replicate analyses of the samples, the average standard deviation of the measurements was better than 0.05 ‰. TA was measured using an automated Radiometer (TitraLab865 ® ) potentiometric titrator and a Red Rod ® combination pH electrode (pHC2001) at McGill University. The diluted HCl titrant was calibrated prior, during and after each titration session using certified reference materials (CRMs) provided by Andrew Dickson (Scripps Institution of Oceanography). Raw titration data were processed with a proprietary algorithm designed for shallow endpoint detection. Surface water samples from the Saguenay Fjord and the Upper St. Lawrence Estuary were also analyzed at Dalhousie University using a VINDTA 3C ® (Versatile Instrument for the Determination of Titration Alkalinity, by Marianda) following the method described in Dickson et al. (2007). A calibration of the instrument was performed against CRMs, and the reproducibility of the measurements was better than 0.1 %.
The DIC concentration of samples, recovered in 2016, 2017 and 2018 in the Saguenay Fjord and surface waters of the Upper and Lower St. Lawrence Estuary, were determined at Dalhousie University using the VINDTA 3C ® . In 2014, DIC was determined on board using a SciTech Apollo DIC analyzer. Once thermally equilibrated at 25 • C, 1-1.5 mL of the sample was acidified with 10 % H 3 PO 4 after being injected into the instrument's reactor. The evolved CO 2 was carried to a LI-COR infrared analyzer by a stream of pure nitrogen. A calibration curve was constructed using gravimetrically prepared Na 2 CO 3 solutions, and the accuracy of the measurements was verified using a CRM. Reproducibility was typically on the order of 0.2 %.

Water mass distribution analysis
A combination of transport processes associated with ocean circulation and biogeochemical cycles generally controls the distribution of tracers in the ocean (Chester, 1990). Resolving the effects of mixing and biogeochemical cycling is imperative if one is to evaluate the movement of nutrients and tracers in a water body. An optimum multiparameter (OMP) analysis allows for the determination of the relative contributions of pre-defined source water types (SWTs), representing the parameter values of the unmixed water masses in one specific geographic location, by optimizing the hydrographic data gathered in a given system (Tomczak, 1981). The original OMP algorithm is a linear inverse model that assumes all hydrographic tracers are conservative. The algorithm has since been modified to handle nonconservative properties such as DIC and nutrients by taking into consideration the stoichiometry of microbial respiration and photosynthesis (Dinauer and Mucci, 2018;Karstensen and Tomczak, 1998).
OMP calculates the SWT fractions, x i , for each data point by finding the best linear mixing combination defined by parameters such as T , S P , δ 18 O water , DO, TA and DIC. The contributions from all SWT must add up to 100 % and cannot be negative. Assuming that four SWT (a, b, c and d) are sufficient to characterize the water column structure, and six parameters (T , S, δ 18 O water , DO, TA and DIC) characterize each of these, the following set of linear equations is solved in the classical OMP analysis (MATLAB -version 1.2.0.0; Karstensen, 2013): where T obs , S obs , δ 18 O obs , DO obs , TA obs and DIC obs are the observed values in any given parcel of water and R is their respective associated fitting residual.
To account for potential environmental variability and measurement inaccuracies and allow for the comparison of parameters with incommensurable units, a weighting proce-dure based on covariances between tracers is typically applied. In this study, weights were assigned arbitrarily based on their conservative behaviors and variability (Lansard et al., 2012). Conservative tracers (i.e., S P , TA, δ 18 O water ) were assigned heavy weights, while nonconservative tracers (i.e., T , DO, DIC) were given low weights according to their seasonal variability. For instance, temperatures in the surface waters of the Saguenay River range from 3.1 • C in the winter to 21 • C in the summer. Dissolved oxygen was also considered a nonconservative tracer as it is heavily reliant on temperature and salinity, as well as biological activity. DIC was given an intermediate weight given that it is relatively conservative except in the surface waters, where photosynthesis and air-sea gas exchange take place. Several OMP analyses were carried out using different weights for each parameter, while weighing their conservative behavior appropriately (i.e., highly conservative vs. lightly conservative). Results were not affected significantly.

Source water type definitions
A water mass is, by definition, a body of water having its origin in a particular source region (Tomczak, 1999). An OMP analysis requires the user to define the major water masses contributing to the structure of the water column in the study area. In the context of biogeochemical cycles, a SWT should be defined where the water mass enters the basin, before it enters the mixing region (Karstensen, 2013). Parameter values are preferably extrapolated from hydrographic observations in the water mass formation region or can be found in the literature.
In this study, source water type definitions were derived from property-property diagrams (see Appendix, Fig. A1) of an observational dataset relevant to the Saguenay Fjord: the Saguenay River (SRW), the St. Lawrence Estuary summertime cold intermediate layer (CIL), the Lower St. Lawrence Estuary bottom waters (LSLE) and the St. Lawrence River (SLRW). Each definition was captured relative to the fjord, i.e., each source water type is only appropriate for the fjord and for the period of study. Definitions and weights are reported in Table 1. A seasonality analysis was carried out to ensure SWT definitions were appropriate for the period of study. Insignificant variations were observed in tracers such as δ 18 O, DIC, TA, DO and S P . The only highly variable tracer was T , which was given the lowest possible weight in the OMP analysis.

CO 2 partial pressures
The CO 2 partial pressure in seawater (pCO 2(SW) ) is defined as the pCO 2 in water-saturated air (pCO 2(air) ) in equilibrium with the water sample or the ratio of the CO 2 concentration in solution to the equilibrium concentration at T , P , and S P , multiplied by the actual pCO 2(air) . As direct measurements of the surface mixed layer pCO 2 were not available in September 2014, May 2016, June 2017 and November 2017, it was calculated (pCO 2(SW-calc) ) using CO2SYS (Excel v2.1; Pierrot et al., 2006) and the measured pH (total or NBS/NIST scale; see Appendix B, Tables B1 and B2), DIC (µ mol kg −1 ), in situ temperature ( • C), practical salinity (S P ) and pressure (dbar) as input parameters. When available, soluble reactive phosphate (SRP) and dissolved silicate (DSi) concentrations were also included in the calculations, but their inclusion did not affect the results significantly because their concentrations are relatively low in surface waters (0.49 and 37.0 µM, respectively) and introduce an insignificant error. DIC rather than TA was used as an input parameter to CO2SYS, since the fjord surface waters are enriched in colored dissolved organic carbon (> 4 mg L −1 ) delivered by the Saguenay River and are characterized by a negative organic alkalinity (positive organic acidity) (see below). The carbonic acid dissociation constants (K * 1 and K * 2 ) of  were used for the calculations, as the latter were found to be more suitable for the low-salinity waters encountered in estuarine environments, such as the Saguenay Fjord (S P < 20) (Dinauer and Mucci, 2017). pCO 2(SW-calc) values were computed for the surface mixed layer located above the sharp pycnocline (∼ 10 m), where most physical and chemical properties are directly impacted by biological activity (photosynthesis and respiration), as well as heat and gas exchange across the air-sea interface (Table 2). Direct measurements of pCO 2 (pCO 2(SW-meas) ) were acquired in May 2018, and pCO 2(SW-calc) was also calculated from pH and DIC for this sampling month for comparison purposes, following the aforementioned procedure.

CO 2 flux across the air-sea interface
The difference between the air and sea surface pCO 2 values ( pCO 2 = pCO 2(SW) -pCO 2(air) ) determines the direction of gas exchange and whether the surface mixed layer of a body of water is a source or a sink of CO 2 for the atmosphere. The air-sea CO 2 gas exchange, or CO 2 flux, can be estimated 554 L. Delaigue et al.: Spatial variations in CO 2 fluxes in the Saguenay Fjord at each station using the following relationship: where F is the flux of CO 2 across the air-sea interface in mmol m −2 d −1 , k is the gas transfer velocity of CO 2 in cm h −1 (Wanninkhof, 1992), K 0 is the solubility of CO 2 in mol kg −1 atm −1 at the in situ temperature and salinity of the surface waters (Weiss, 1974), and pCO 2 is the difference between the air and sea surface pCO 2 values in µatm. Whereas, formally, Fick's first law of diffusion should be written as F = −D δC/δx (where F is the diffusion flux in mole s −1 m −2 , D is the diffusion coefficient in m 2 s −1 , C is the concentration of CO 2 in mole m −3 and x is the distance in m), as commonly expressed by Eq. (3a), positive values of F indicate the release of CO 2 to the atmosphere by surface waters, whereas negative values imply that surface waters serve as a sink of atmospheric CO 2 . The flux of CO 2 was computed for each sampling month, using the pCO 2(air) for each sampling date (395 µatm for September 2014, 407 µatm for May 2016, 408 µatm for June and November 2017, and 411 µatm for May 2018; see below for details).
The gas transfer velocity of CO 2 was calculated using the revised relationship of Wanninkhof (2014): where u is the wind speed (m s −1 ) and Sc is the Schmidt number (Wanninkhof, 2014). Wind speed was estimated using the hourly station wind speed data from Environment Canada at the La Baie weather station (Fig. 1a) for each sampling month. The Schmidt number is defined as the kinematic viscosity of water divided by the diffusion coefficient of CO 2 . Sc was corrected for the temperature dependence of CO 2 in freshwater (S P = 0), assuming that k is proportional to Sc −1/2 (Wanninkhof, 1992). In the case of CO 2 , the increase in Sc −1/2 (and k) with increasing temperature is compensated for by a decrease in solubility; therefore, k was considered nearly temperature independent (Wanninkhof, 1992). Sc was computed using the following equation: where t is the temperature ( • C) and A, B, C, D, and E are fitting coefficients for seawater (S P = 35) and freshwater (S P = 0), for water temperatures ranging from −2 to 40 • C (Wanninkhof, 2014). The uncertainty in Sc ranges from 3 % to 10 % and is mainly due to the imprecision of diffusion coefficients (Wanninkhof, 2014). Estimates of k, calculated at each sampling point using the equation of Wanninkhof (2014), ranged from 0.36 to 3.38 cm h −1 for the fjord, compared to 1.6 to 4.5 cm h −1 in the St. Lawrence Estuary (Dinauer and Mucci, 2017). Atmospheric pCO 2 values (pCO 2(air) ) were computed using the daily averages of measured mole fractions of CO 2 in dry air, obtained at the La Baie weather station and retrieved from the Climate Research Division at Environment and Climate Change Canada. The mean pCO 2(air) was then calculated for each year using the following equation: where xCO 2 is the measured mole fraction of CO 2 in dry air in ppm, P b is the barometric pressure at the sea surface in atm, and P w is the saturation water vapor pressure at in situ temperature and salinity in atm. P b was obtained using the conversion formula of Tim Brice and Todd Hall (from NOAA's National Weather Service -https://www.weather. gov/epz/wxcalc_wxcalc2go, last access: 3 November 2019), using the La Baie weather station's elevation (152 m). P w was calculated using the Rivière-à-Mars properties (i.e., closest body of water to the weather station), and the P w calculated from its relationship to T and S P provided by Weiss and Price (1980). The area-averaged CO 2 flux (F area-avg ) was computed for the whole fjord, following the procedure described by Jiang et al. (2008): where F i is the average of all the fluxes within segment i and S i is the surface area of segment i. The fjord was divided into two segments, one including the inner basin and the other encompassing the two outer basins, as each segment often displays distinct behaviors. Segments are identified in Fig. 1b. The fjord's surface area (∼ 290 km 2 ) was computed using a land mask in MATLAB.

Water mixing model
A two end-member mixing model was constructed based on the chemical properties of the freshwater delivered to the fjord (Saguenay River) and marine bottom waters entering the fjord from the St. Lawrence Estuary (Fig. 2a). As shown in the results of the OMP analysis (Sect. 3.1), the LSLE and SLRW have a negligible influence on the fjord's water structure and thus were not included in the model. Given that the carbonate chemistries of the CIL and LSLE waters are similar, the bottom waters were assumed to be well mixed and constitute a single end-member. This is illustrated in Fig. 2, as the high S P end-member alkalinity extends linearly beyond that of the CIL end-member (Table 1). The measured surface TAs were strongly correlated to S P (R 2 = 0.999) in the fjord waters. Therefore, end-member properties were obtained by extrapolating the surface water (above the pycnocline) data to S P = 0 and bottom water data to the highest measured salinity (Fig. 2a). The extrapolated TA (meas) (Fig. 2b; 154 µmol kg −1 ) is in good agreement with the average TA (meas) of samples taken directly from the Saguenay River in 2013 and 2017 (157 µmol kg −1 ). The organic alkalinity of the fjord waters was estimated from the difference between the measured and calculated TA (TA (calc) ; Fig. 2b).  (calc) and Org Alk definitions for the Saguenay River (SRW), using surface water data from all sampling months, with standard error. The Org Alk (positive) contribution to the TA of the CIL is not considered, as it accounts for less than 0.1 % of its TA.
The latter was calculated using CO2SYS (Excel v2.1;Pierrot et al., 2006) and pH and DIC as input parameters. The end-member source waters were then mixed, assuming that TA (calc) and DIC behave conservatively. Hence, the salinity, total alkalinity (TA (mix) ) and dissolved inorganic carbon (DIC (mix) ) of the mixed solutions were calculated using the following equations: where m i is the mass contribution of each end-member to the mixture. pCO 2(SW-mix) was then computed from TA (mix) and DIC (mix) for practical salinities ranging from 0 to 33 at four different temperatures (0, 5, 10 and 15 • C) using CO2SYS. Results of the model (Fig. 8) show that, at the lower and higher salinities, the pCO 2(SW-mix) is elevated and the fjord serves as a net source of CO 2 to the atmosphere, but at intermediate salinities (5 < S P < 15) or mixing ratios, the fjord may serve as a net sink of atmospheric CO 2 when surface water temperatures are close to freezing. The data from the various cruises are superimposed on the model results, after correction for the organic alkalinity.

Salinity normalization of DIC in surface waters
To quantitatively evaluate the impact of biological activity on the DIC budget in the surface waters of the fjord, DIC and TA (calc) were normalized to the average surface salinity of each sampling month (S P = 12.4 for September 2014; S P = 2.58 for May 2016; S P = 7.61 for June 2017; S P = 10.9 for November 2017; and S P = 5.9 for May 2018) following the procedure of Friis et al. (2003): where DIC meas is the measured DIC, DIC S=0 is the DIC extrapolated to S P = 0, S meas is the measured practical salinity and S ref is the average measured practical salinity per sampling month (Friis et al., 2003). The change in NDIC (i.e., NDIC) along the fjord, relative to the waters at the head of the fjord, was then computed for each sampling month. These values reveal how DIC evolves along the fjord beyond what is expected based on conservative mixing.

Oxygen saturation and apparent oxygen utilization in the surface waters
To further account for the biological activity in the surface waters, the oxygen saturation index was calculated for each sampling month in the surface waters of the fjord using the following equation: where [O 2 ] meas is the dissolved oxygen concentration measured in the fjord waters and [O 2 ] equil is the equilibrium dissolved oxygen concentration (or solubility) at in situ conditions (i.e., temperature and salinity) for each sample.
The oxygen saturation index indicates if the system is autotrophic (i.e., production of oxygen, dominated by photosynthesis) or heterotrophic (consumption of oxygen, dominated by microbial respiration). The oxygen saturation remains a qualitative proxy, as O 2 exchange at the air-sea interface is about 9 times faster than CO 2 exchange (Zeebe and Wolf-Gladrow, 2001 (Tomczak and Large, 1989;Tomczak, 1981;Mackas and Harrison, 1997). A variational analysis (DIVA) interpolation was applied between field data points in Ocean Data View.

Water mass analysis
Relative contributions (mixing ratios, f ) of the Saguenay River (SRW), the St. Lawrence Estuary summertime cold intermediate layer (CIL) and Lower St. Lawrence Estuary (LSLE) bottom waters throughout the Saguenay Fjord's water column for the sampling month of June 2017 are shown in Fig. 3. As expected, the SRW and CIL are dominant contributors, with the SRW forming a brackish surface layer (f = 1 in surface waters) and the CIL replenishing the bottom waters of the fjord (0.7 < f < 1). According to the OMP analysis, the LSLE bottom waters have a small contribution to the fjord's bottom waters (f = 0.2), adding to the complexity of the water structure. Although somewhat unexpected, this can readily be explained by tidal upwelling, internal waves and intense turbulent mixing of the water column resulting from the rapid shoaling at the head of the Laurentian Channel (Gratton et al., 1988;Saucier and Chassé, 2000). The relative contribution of the LSLE bottom waters in the deep waters of the fjord is small and could only be detected because of the suite of geochemical and isotopic tracers used in the OMP analysis, especially the difference in the δ 18 O water signature of the CIL and LSLE waters. The contribution from the St. Lawrence River Water (SLRW) is negligible, as it intrudes slightly at the surface at the mouth of the fjord and is thus not shown here. Although the water column structure is similar throughout the year, seasonal variations do occur and will be addressed in a forthcoming paper.

Aqueous pCO 2 and CO 2 flux
Variations in the inorganic carbon chemistry in the Saguenay Fjord water column are described using field data ac- . It was inversely proportional to the salinity of the surface waters of the fjord and became positive, yet a negligible fraction (< 0.1 %) of TA (corr) when S P > 25, like in the St. Lawrence Estuary. The negative organic alkalinity of the Saguenay River water most likely originates from soil humic acids that are flushed by percolation with groundwaters that drain the metamorphic and igneous rocks of the Canadian Shield. Surface water pCO 2 (pCO 2(SW-calc) ) values were higher at the head of the fjord (i.e., near the Saguenay River mouth) and lower at the mouth of the fjord, although large variations (315 to 740 µatm -average 503 µatm) were observed on a seasonal and yearly basis ( Table 2). Values of pCO 2(SW) were higher in May 2018 (623 µatm), June 2017 (506 µatm), and May 2016 (563 µatm) than in November 2017 (418 µatm) and September 2014 (406 µatm). This can be explained by the larger freshwater discharge from the Saguenay River in the spring (i.e., spring freshet, average of 1856 ± 21 m 3 s −1 for spring periods of 1998-2018) compared to the fall (1470 ± 10 m 3 s −1 for fall periods of 1998-2018). As atmospheric pCO 2(air) varied marginally between September 2014 (395 µatm) and May 2018 (411 µatm), the fjord was generally a source of CO 2 to the atmosphere near its head (i.e., surface pCO 2 values above atmospheric level), while the zone near its mouth was most often a sink (i.e., surface pCO 2 values below atmospheric level) (Fig. 4). An anomaly was observed in November 2017, with a high pCO 2(SW-calc) value (> 550 µatm) near the mouth of the fjord. Given the statistics of the box plot presented in Fig. 7, this value appears to be erroneous.
Air-sea CO 2 fluxes within the fjord ranged from −2.4 to 10.0 mmol m −2 d −1 (Fig. 6). Near the head of the fjord, fluxes were mostly positive, while values decreased when  approaching its mouth. Overall, the total area-averaged degassing flux of the fjord adds up to 2.14±0.43 mmol m −2 d −1 or 0.78 ± 0.16 mol m −2 yr −1 . In comparison, the degassing flux in the adjacent St. Lawrence Estuary was estimated at between 0.36 and 0.74 mol m −2 yr −1 during the late spring and early summer (Dinauer and Mucci, 2017). This discrepancy can be explained by the low carbonate alkalin-ity (and buffer capacity) of the Saguenay River waters that flow through the Grenvillian metamorphic and igneous rocks of the Canadian Shield (Piper et al., 1990), as with most rivers on the north shore of the St. Lawrence Estuary (e.g., Betsiamites, Manicouagan, Romaine; Paul del Giorgio, personal communication, 2016), and the low productivity of the fjord surface waters because of very limited light penetration due to their high chromotrophic dissolved organic matter (CDOM) content (Tremblay and Gagné, 2009;Xie et al., 2012). In contrast, waters of the St. Lawrence River have an elevated carbonate alkalinity (∼ 1200 µM), inherited from the Ottawa River that drains through limestone deposits (Telmer and Veizer, 1999). Furthermore, the estuary is host to multiple seasonal phytoplankton blooms (Levasseur and Therriault, 1987;Zakardjian et al., 2000;Annane et al., 2015) that strongly modulate its trophic status (Dinauer and Mucci, 2018).
The correlation between pCO 2(SW-meas) and pCO 2(SW-calc) is presented in Fig. 5. The average difference between pCO 2(SW-meas) and pCO 2(SW-calc) is 48 µatm, implying that calculations underestimate pCO 2(SW) values by approximately 7 % and thus contribute to the uncertainty associated with CO 2 fluxes. This discrepancy most likely originates from uncertainties associated with the carbonic acid dissociation constants (K * 1 and K * 2 ) in low-salinity estuarine environments, particularly those affected by strong organic alkalinities or acidities such as in the Saguenay Fjord Ko et al., 2016). This concurs with the results of Lueker et al. (2000), who showed that, depending on the choice of K * 1 and K * 2 , computed pCO 2(SW) values from other carbonate system parameters (TA, DIC, pH) can be up to 10 % lower than those of direct measurements. Consequently, although the constants of  are the most suitable for this study, direct measurements of the pCO 2(SW) should preferentially be carried out whenever possible.

Water mixing model approach
As results of the OMP analysis reveal, LSLE and SLRW have a negligible influence on the water properties in the fjord, except for the latter near the mouth. Additionally, given the relatively small contribution of the LSLE deep waters and their similarity to the carbonate chemistry of the CIL, their influence is considered inconsequential on the properties of the mixture. Hence, a conservative mixing model was constructed based on the chemical properties of the two main source water masses in the fjord (i.e., SRW and the CIL mixture for bottom waters) and the relationship between practical salinity and TA (corr) /DIC, respectively (Fig. 8). pCO 2(SW-calc) were normalized at each station to the average surface water temperature per sampling month (i.e., T = 10.4 • C for September 2014, T = 5.04 • C for May 2016, T = 11.9 • C for June 2017, T = 7.13 • C for November 2017 and T = 5.08 • C for May 2018) to account for the effects  of temperature on the CO 2 solubility in water, following the procedure described in Jiang et al. (2008). The temperaturenormalized pCO 2(SW-calc) values, pCO 2(SW-SST) , from the various cruises were superimposed on the model results in Fig. 8.
Field measurements follow the trend displayed by the mixing model. The fjord appears to be a net source of CO 2 to the atmosphere during periods of high freshwater discharge (i.e., spring freshet) and a net sink at intermediate surface salinities (5 < S P < 15). This is consistent with the weak buffer capacity of the freshwater. Given the short residence time of surface waters in the Saguenay Fjord (∼ 1.5 d), the influence of gas exchange across the air-sea interface is negligible on the DIC pool. Likewise, Dinauer and Mucci (2017) reported that the surface waters in the St. Lawrence Estuary near Tadoussac (at the mouth of the fjord) are highly supersaturated in CO 2 with respect to the atmosphere and only the highly productive waters of the Lower St. Lawrence Estuary manage to draw down the surface pCO 2 to near atmospheric values. In other words, degassing of the metabolic CO 2 accumulated in the river and upper estuary is slow. Thus, changes in temperature-normalized pCO 2 primarily reflect changes in DIC by mixing and biological activity. Hence, discrepancies between results of the mixing model and field measurements can be ascribed to microbial respiration and photosynthesis.
In May 2016, the surface waters of the fjord were clearly supersaturated in oxygen (Fig. 9), implying that photosynthesis dominated over respiration. This would explain the rapid seaward (increasing S P ) decrease in pCO 2(SW-SST) , faster than the mixing model predicts (Fig. 8), and the strong negative NDIC (i.e., change in NDIC relative to the saline waters at the head of the fjord) throughout the fjord (Fig. 10), as CO 2 (i.e., DIC) is taken up by photosynthesizing organisms -most likely diatoms (Chassé and Côté, 1991). In May 2018, surface waters were slightly undersaturated in oxygen, between 90 % and 100 % saturation, and NDIC was positive over most of the fjord. Very similar trends were observed in June 2017, with near-saturation oxygen concentrations (between 95 % and 101 % saturation) and mostly positive NDIC values throughout the transect. Thus, during these sampling periods, biological activity was dominated by microbial respiration (Fig. 10), elucidating the minor deviation between the pCO 2(SW-SST) and the model results (Fig. 8), especially near the head of the fjord. Additionally, it is interesting to note that NDIC is chronically negative for all sampling months near the 45 km mark.
The difference between the May 2016 and May 2018 biological responses to the spring freshet can potentially be explained by the difference in total freshwater discharge from the Saguenay River to the fjord. The freshwater discharge in May 2018 was approximately 20 % larger than in May 2016. Whereas the surface salinities recorded in both years throughout most of the fjord were nearly identical (S P = 0.5-4 in 2016, S P = 0.7-5 in 2018), the greater delivery of soil porewater and associated CDOM in May 2018 may have inhibited local productivity due to light absorption by CDOM (Lavoie et al., 2007). Consistent with this interpretation is the fact that the pCO 2(SW-SST) at the head of the fjord (St. Fulgence) in May 2016 was slightly higher than in May 2018 but was drawn down much faster downstream (Fig. 8).
There does not appear to be a clear biological signal in the November 2017 data, as little variation is observed between the measured and modeled pCO 2 . Furthermore, Fig. 10 indicates that neither respiration nor photosynthesis dominated during this period as the NDIC varies between weakly positive and negative values. Likewise, the September 2014 data reveal little biological activity, although a slight dominance of respiration (i.e., positive NDIC with a drop in %O 2 saturation) can be observed near the mouth of the fjord (Fig. 10), hence explaining the slight deviation from the mixing model.
These results highlight the importance of the freshwater plume from the Saguenay River in regulating the pCO 2 dynamics in the fjord. Winds, in addition to regulating the gas exchange coefficient, are also known to have a direct influence on air-sea CO 2 fluxes by driving upwelling of CO 2rich waters along with the entrainment of nutrients in surface waters, thus increasing biological activity (Wanninkhof and Triñanes, 2017). However, wind speeds are relatively low in the studied system (1.89 m s −1 < u < 4.2 m s −1 , Table 2), implying a calm sea state (Frankignoulle, 1998), hence reinforcing that changes in pCO 2(SW-SST) can mainly be attributed to microbial respiration and photosynthesis modulated by water renewals rather than winds.

Summary and conclusions
Results of the OMP analysis reveal that SRW and CIL are the dominant source water types to the fjord and determine the structure of its water column. Mixing of marine waters with SRW at the head of the fjord leads to the formation of a brackish surface layer (wedge), while the CIL replenishes the bottom waters of the fjord. The analysis further unveiled a small contribution of the LSLE bottom water to the bottom waters of the fjord, adding to the complexity of the water column structure. The SLRW has a negligible influence on the water properties in the fjord, except near its mouth -sampling of the very turbulent waters directly over the fjord's first and shallowest sill would help refine the contribution of the SLRW to the fjord's surface waters.
The magnitude and sign of the pCO 2 across the airwater interface in the Saguenay Fjord, mostly determined by the pCO 2(SW) , as pCO 2(air) varied only slightly over the sampling period, are mostly modulated by the freshwater discharge and the salinity of the surface waters. The surface waters of the fjord are a source of CO 2 to the atmosphere at high freshwater discharge and a sink of CO 2 at intermediate surface salinities (5 < S P < 15), especially at near-freezing temperatures. Direct measurements of surface water pCO 2 acquired in May 2018 confirmed this trend. Biological activity alters the surface water pCO 2 , with both photosynthesis and respiration impacting the waters depending on the Figure 10. (a) NDIC (i.e., change in NDIC relative to the saline waters at the head of the fjord) distribution in the Saguenay Fjord surface waters. Data were normalized to a common salinity (average for each sampling period) according to the method of Friis et al. (2003). Error bars show standard error of the mean for NDIC values. (b) Apparent oxygen utilization (AOU) against NDIC. Error bars show standard error of the mean for AOU values. sampling month. This conclusion is supported by the oxygen saturations observed in the surface waters of the fjord, as well as the downstream NDIC trend along the fjord's main axis. Given the short residence time of surface waters in the Saguenay Fjord (∼ 1.5 d), the influence of gas exchange on spatial variations in the pCO 2 across the air-sea interface along the main axis of the fjord is negligible. Overall, the fjord serves as a source of CO 2 to the atmosphere during the ice-free season, with an average yearly outgassing flux of 0.78 ± 0.16 mol m −2 yr −1 .
A study of the dissolved inorganic carbon budget of fjords not only affords information on their carbon dioxide levels (i.e., source or sink of CO 2 with respect to the atmosphere) and surface water chemistry but also provides insights into the magnitude of gas exchange and the amount of biological activity it sustains. In addition to biological production, upwelling, water temperature and the spreading of freshwater plumes all regulate pCO 2 in coastal systems. Wind speed is also critical as it impacts sea state and the efficiency of gas exchange at the air-sea interface (Chen et al., 2013). Hence, the importance of wind on controlling the CO 2 flux needs to be further investigated, especially in high-latitude fjords where strong winds are often focused along narrow channels between steep cliffs and narrow inlets.
Anthropogenic activities, through climate change, are altering the continental water cycle, along with the flows of carbon, nutrients and sediment to the coastal oceans (Borges, 2005) and thus the sequestration of anthropogenic CO 2 by the oceans. In glacial fjords, freshwater discharge will be modified by accelerated glacier melting and their ultimate demise. Knowledge of the spatial and temporal variability of glacial discharge and seasonal ice cover in high-latitude fjords is essential to estimate their influence on circulation patterns and freshwater export, both of which impact local productivity, terrestrial carbon burial and export as well as CO 2 fluxes in these complex ecosystems. Finally, an improved understanding of the coastal carbon cycle will require a more comprehensive spatial and temporal coverage of surface mixed layer pCO 2 data in these environments, ideally from direct measurements.  Data availability. Data presented in this paper are available upon request from one of the authors (louise.delaigue@mail.mcgill.ca).
Author contributions. AM and LD conceived the project. AM acquired and processed the data prior to 2016. LD conducted the data analysis and wrote the first draft of the paper, whereas AM provided editorial and scientific recommendations. HT provided results of alkalinity and dissolved inorganic carbon analyses and scientific recommendations.