Marine CO₂ system variability along the northeast Pacific Inside Passage determined from an Alaskan ferry

Information on marine CO₂ system variability has been limited along the northeast Pacific Inside Passage despite the region's rich biodiversity, abundant fisheries, and developing aquaculture industry. Beginning in 2017, the Alaska Marine Highway System M/V Columbia has served as a platform for surface underway data collection while conducting twice weekly ∼1600 km transits between Bellingham, Washington, and Skagway, Alaska. Marine CO₂ system patterns were evaluated using measurements made over a 2-year period, which revealed the seasonal cycle as the dominant mode of temporal variability. The amplitude of this signal varied spatially and was modulated by the relative influences of tidal mixing, net community production, and the magnitude and character of freshwater input. Surface water pH_T (total hydrogen ion scale) and aragonite saturation state (Ω_arag) were determined using carbon dioxide partial pressure (pCO₂) data with alkalinity derived from a regional salinity-based relationship, which was evaluated using intervals of discrete seawater samples and underway pH measurements. High-pCO₂, low-pH_T, and corrosive Ω_arag conditions (Ω_arag<1) were seen during winter and within persistent tidal mixing zones, and corrosive Ω_arag values were also seen in areas that receive significant glacial melt in summer. Biophysical drivers are shown to dominate pCO₂ variability over most of the Inside Passage except in areas highly impacted by glacial melt. pH_T and Ω_arag extremes were also characterized based on degrees of variability and severity, and regional differences were evident. Computations of the time of detection identified tidal mixing zones as strategic observing sites with relatively short time spans required to capture secular trends in seawater pCO₂ equivalent to the contemporary rise in atmospheric CO₂. Finally, estimates of anthropogenic CO₂ showed notable spatiotemporal variability. Changes in total hydrogen ion content ([H⁺]_T), pH_T, and Ω_arag over the industrial era and to an atmospheric pCO₂ level consistent with a 1.5 ^∘C warmer climate were theoretically evaluated. These calculations revealed greater absolute changes in [H⁺]_T and pH_T in winter as opposed to larger Ω_arag change in summer. The contemporary acidification signal everywhere along the Inside Passage exceeded the global average, with specific areas, namely Johnstone Strait and the Salish Sea, standing out as potential bellwethers for the emergence of biological ocean acidification (OA) impacts. Nearly half of the contemporary acidification signal is expected over the coming 15 years, with an atmospheric CO₂ trajectory that continues to be shaped by fossil–fuel development.

1 Introduction

Atmospheric carbon dioxide (CO₂) has increased over the industrial era from 278 ppm in 1765 to 414 ppm in 2020 due to the emissions of CO₂ from fossil fuel combustion and land use change, which combined have mobilized a total of 690 ± 80 Gt of carbon (Friedlingstein et al., 2021). So far over the industrial era, an estimated 180 ± 35 Gt of this carbon pool has transferred into the ocean (Friedlingstein et al., 2021), known as the oceanic anthropogenic CO₂ component (Sabine et al., 2004), and led to changes in the marine CO₂ system, including reduced carbonate ion content ([CO${}_{\mathrm{3}}^{\mathrm{2}-}$]) and pH_T (total hydrogen ion scale) and increased total hydrogen ion content ([H⁺]_T) and CO₂ partial pressure (pCO₂). These marine CO₂ system changes are collectively referred to as “ocean acidification” (Caldeira and Wickett, 2003; Doney et al., 2009; Feely et al., 2004a, 2009), and two recent assessments estimate an average pH_T decline for the global surface ocean on the order of 0.1 units over the industrial era (Jiang et al., 2019; Lauvset et al., 2020). In conjunction with this pH_T decline, reductions in [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] have simultaneously decreased the saturation states (Ω) of carbonate biominerals, with aragonite as the most soluble carbonate biomineral typically targeted in biological studies investigating the effects of ocean acidification (OA). Ω_arag is a ratio of the product of [CO${}_{\mathrm{3}}^{\mathrm{2}-}$] and calcium content over the solubility product for aragonite, and this ratio dictates the thermodynamic favorability of aragonite precipitation. If Ω_arag is >1, precipitation is favored over dissolution. Globally, average surface Ω_arag is estimated to have declined by 0.53 units (Lauvset et al., 2020). These assessments of global average pH_T and Ω_arag decline over the industrial era are based on calculations of anthropogenic CO₂ content; however, long-term change in both pH_T and Ω_arag resulting from anthropogenic CO₂ input is captured in multidecadal open-ocean time series datasets (Bates et al., 2014; Doney et al., 2020; Franco et al., 2021).

Along the continental margins, seawater conditions may not track the global average surface ocean pH_T decline, particularly in the northeast Pacific where seawater is less buffered than in some other ocean regions, thereby making it more sensitive to increasing anthropogenic CO₂ (Feely et al., 2008, 2018; Lauvset et al., 2020; Cai et al., 2020; Jiang et al., 2015). One estimate of pH_T decline on this margin suggests twice the global average based on fossil foraminifera shells preserved in marine sediments (Osborne et al., 2020). Even with the potential for larger pH_T decline along the northeast Pacific margin, putting this change into context can be challenging. Given that pH_T is negative log₁₀ of [H⁺]_T, the absolute change in [H⁺]_T varies based on the initial pH_T for the same degree of pH_T decline (Fassbender et al., 2021, 2017). For example, a 0.1-unit pH_T decrease with an initial pH_T of 7.6 will result in an absolute [H⁺]_T change of 6.3 nmol kg⁻¹, whereas the same degree of pH_T decrease with an initial pH_T of 8.4 will drive a 1 nmol kg⁻¹ [H⁺]_T change. This clarification is important because the absolute change in acidity can be different despite the same relative change in pH_T, and confusion may be enhanced when considering that some continental margins likely have experienced different relative pH_T change compared to the global surface ocean average (Osborne et al., 2020; Evans et al., 2019; Pacella et al., 2018; Salisbury and Jönsson, 2018), including in some coastal ecosystems currently being evaluated for their OA mitigation potential (Ricart et al., 2021; Kroeker et al., 2021).

The magnitude of OA-driven marine CO₂ system changes and an ecosystem's mitigation potential are both critical areas of research because negative impacts are already being felt by some vulnerable marine species. Along the northeast Pacific continental margin, larval shellfish mortality within hatcheries has been tied directly to low Ω_arag (Barton et al., 2012), and some adaptation measures to avoid such conditions have been developed (Barton et al., 2015). Other shell-forming marine species in this region are also exhibiting impacts from OA, including Dungeness crab (Bednarsek et al., 2020; Berger et al., 2021) and pteropods (Bednarsek et al., 2017, 2021; Mekkes et al., 2021). The general consensus is that calcifying species may be the most directly impacted (Kroeker et al., 2013; Haigh et al., 2015; Marshall et al., 2017), although sensitivity to OA appears to be very species, life stage, and population specific (Doney et al., 2020) with the potential for compensatory mechanisms to help sustain populations (Peck et al., 2018; Bednarsek et al., 2021). However, there is a high likelihood of enhanced vulnerabilities by other co-occurring stressors like warming (Kroeker et al., 2013) and reduced oxygen content (Gobler and Baumann, 2016). Biological stressors, such as viral pathogens and harmful algal species, may also become more prevalent or virulent in association with changes in marine CO₂ chemistry and warming (Raven et al., 2020; Asplund et al., 2013). The sum of both the direct and indirect effects from OA and other co-stressors threatens marine food webs (Jin et al., 2020), harvested species (Ekstrom et al., 2015), and dependent coastal communities (Mathis et al., 2015); understanding this threat demands assessing how the marine CO₂ system has and will evolve through time.

Determining long-term trends in coastal settings is difficult because of inherent high variability resulting from a number of processes unique to the land–ocean interface. Physical forcing from upwelling-favorable winds or tide-induced vertical mixing can result in surface water pCO₂ that is super-saturated with respect to the atmosphere, whereas high rates of primary production draw down surface water pCO₂ to well below atmospheric levels. Additionally, freshwater input from land can act to dilute total dissolved inorganic carbon (TCO₂) and total alkalinity and reduce pCO₂ (Meire et al., 2015), or, alternatively, increase pCO₂ through the respiration of riverine organic matter (Ward et al., 2017). These processes all occur on different timescales, are not uniformly important across coastal settings, and collectively act to make resolving the relatively small anthropogenic CO₂ signal difficult to disentangle from environmental variability. Resolving environmental variability, even to the point of capturing seasonal cycles, remains a challenge in many settings due to a lack of measurements (Hales et al., 2008). In the northeast Pacific between British Columbia (BC) and southeast Alaska (AK), modeling efforts have aided in addressing this knowledge gap and have indicated the relative significance of freshwater input (Siedlecki et al., 2017; Hauri et al., 2020) and its source character (Pilcher et al., 2016), as well as projected warming, deoxygenation, and acidification on multidecadal timescales (Holdsworth et al., 2021). However, observations remain essential to evaluate model output and confirm our understanding of the governing processes that shape marine CO₂ system variability, particularly in nearshore settings.

2 Study region

The Inside Passage is a roughly 1600 km network of coastal waterways that spans the semi-enclosed Salish Sea, the central and northern BC coast, and southeast AK. This nearshore region is a key interface between the Pacific coastal temperate rainforest (O'Neel et al., 2015; Bidlack et al., 2021) and a highly productive continental shelf ecosystem (Ware and Thomson, 2005; Jackson et al., 2015). The area experiences a wide array of physical and biogeochemical drivers including intense tidal currents within narrow passages that induce persistent vertical mixing (Whitney et al., 2005; Dosser et al., 2021), strong autumn and winter storms (Stabeno et al., 2004), high runoff from rainfall and snowmelt and glacial-melt sources (Morrison et al., 2012; Beamer et al., 2016; Edwards et al., 2020; Neal et al., 2010), high terrestrial organic carbon input (Edwards et al., 2020; Oliver et al., 2017; St. Pierre et al., 2021), and remotely forced influences such as El Niño events and marine heat waves (Bond et al., 2015; Jackson et al., 2018).

The coastal ocean from BC to southeast AK has large under-sampled areas (Hales et al., 2008; Evans and Mathis, 2013) and coarse temporal information on marine CO₂ system variability based on direct measurements (Evans and Mathis, 2013; Evans et al., 2012; Tortell et al., 2012) except within the Salish Sea where seasonal and spatial patterns are more constrained (Evans et al., 2019; Cai et al., 2021; Ianson et al., 2016; Fassbender et al., 2018a; Feely et al., 2010; Lowe et al., 2019). We reduced this information gap by outfitting a passenger ferry within the Alaska Marine Highway System (AMHS) fleet, the M/V Columbia, with instrumentation to monitor surface ocean conditions along the Inside Passage (Fig. 1). We report here on surface underway measurements made from November 2017 to October 2019, and we use this dataset to describe marine CO₂ system patterns and quantify the relative importance of key drivers in shaping the observed variability. We also evaluate marine CO₂ system extremes and their timing along the ferry transit, which likely has implications for the exposure histories of vulnerable species. Finally, we estimate the anthropogenic CO₂ content accrued over the industrial era; assess the impact this perturbation has had on [H⁺]_T, pH_T, and Ω_arag; and theoretically gauge the extent of acidification implied by the Paris Agreement (UNFCC, 2015) to limit global warming to preferably 1.5 ^∘C relative to pre-industrial levels, what we refer to as the “1.5 ^∘C acidification level”. The so-called “remaining 1.5 ^∘C carbon budget” translates to an atmospheric CO₂ level which would be reached with all potential mitigation pathways (Rogelj et al., 2018) and therefore can be viewed as the best case scenario for the maximum acidification owing to anthropogenic CO₂ input.

Figure 1 Surface salinity along the Inside Passage expressed as the coefficient of variation (CV; %; a) computed for $\mathrm{0.03}{}^{\circ }\times \mathrm{0.03}{}^{\circ }$ monthly grid cells from underway measurements made from the M/V Columbia between 3 November 2017 and 2 October 2019. Areas of the highest salinity CV are due to large freshwater input, and black circles with labels mark the Alaska Marine Highway System terminals: Bellingham, WA (B); Ketchikan, AK (K); Wrangell, AK (W); Petersburg, AK (P); Juneau, AK (J); Haines, AK (H); Skagway, AK (Sk); and Sitka, AK (Si). Also shown is the surface chlorophyll CV (%; b) from Moderate Resolution Imaging Spectroradiometer (MODIS) Level 3, 4 km mapped data from February to October 2018 and 2019. Areas of high chlorophyll CV reflect instances of biomass accumulation presumably owing to periods of high primary productivity. Areas labeled in the figure panel are the states of Washington (WA) and Alaska (AK), the province of British Columbia (BC), Lynn Canal (LC), Sergius Narrows (SN), Wrangell Narrows (WN), Dixon Entrance (DE), Johnstone Strait (JS), Vancouver Island (VI), and the Salish Sea (SS).

3 Methods

3.1 Underway instrumentation

The AMHS M/V Columbia transited the 1600 km Inside Passage on a weekly basis (Sect. S1). Surface (∼2 m) seawater pCO₂ data were obtained from measurements of CO₂ mixing ratio (xCO₂) made using a General Oceanics 8050 (GO8050) pCO₂ measuring system following recommended protocols (Pierrot et al., 2009). Seawater was drawn into the M/V Columbia through an intake located in the bow thruster engine room and supplied to the GO8050 and ancillary sensors using a $\frac{\mathrm{1}}{\mathrm{2}}$ HP self-priming centrifugal pump (AMT 429A-98 or similar) located ∼2 m from the seawater intake. Temperature measurements were made at the seawater intake using an SBE 38 digital oceanographic thermometer with an accuracy reported by Sea-Bird Electronics of 0.001 ^∘C. Seawater was then circulated from the bow thruster room up one deck to the car deck and then aft approximately 30 m along the starboard side to where the GO8050 and ancillary sensors were installed. The seawater circulation loop was split at this location between the GO8050 wet box and an ancillary sensor loop consisting of an SBE 45 MicroTSG thermosalinograph and an Aanderaa 4330F oxygen optode. The accuracy of the temperature measurement from the SBE 45 was reported as 0.001 ^∘C when interfaced with the SBE 38, and the accuracy of the conductivity measurement was 0.0003 S m⁻¹. Salinity, computed from conductivity and temperature, is reported here on the Practical Salinity Scale (PSS-78). The accuracy of the Aanderaa 4330F oxygen optode reported by the manufacturer Xylem was <1.5 %. All ancillary sensors were serviced annually, and a multipoint calibration was conducted on the oxygen optode at the Aanderaa facility in Norway. Oxygen data from the Aanderaa 4330F were output in µmol L⁻¹, salinity-corrected using the approach described in Bittig et al. (2018), and then density-corrected to µmol kg⁻¹. Oxygen data are reported here as the difference from saturated values (ΔO₂).

Seawater entered the GO8050 wet box at ∼2.9 L min⁻¹ and ∼10 psi and then was circulated into a water-jacketed (to minimize warming) primary showerhead equilibrator with a liquid volume of ∼0.5 L and to a smaller secondary equilibrator with a liquid volume of ∼0.1 L. The primary equilibrator was maintained at ambient pressure on the car deck by a vent that was plumbed to the secondary equilibrator and then to the primary equilibrator. The pressure difference between inside the primary equilibrator and the car deck was monitored using a Setra pressure transducer (model 239) with a 0.15 hPa uncertainty. The secondary equilibrator serves to pre-equilibrate (make-up) air entering the primary equilibrator through the vent due to any loss through the headspace gas recirculation loop. A flow meter was present at the opening of the vent in order to monitor make-up air flow into the equilibrator. The headspace gas volume of the primary equilibrator was ∼0.8 L, and seawater temperature was monitored within the primary equilibrator using a Fluke thermometer (model 1523) and thermistor probe (model 5610) with an uncertainty of 0.01 ^∘C.

Atmospheric air was drawn from an intake on the foredeck to the GO8050 wet box. Both the equilibrator headspace gas and atmospheric air were dried using a condenser (Peltier thermoelectric cooling device) and Permapure Nafion drying tubes in order to minimize the correction for water vapor content associated with band-broadening within the infrared gas analyzer located in the GO8050 dry box. The analyzer housed in the GO8050 dry box was a LI-COR LI840A CO₂–H₂O gas analyzer with a root-mean-square noise level for 1 Hz measurements of 1 ppm reported by the manufacturer. Dried equilibrator headspace and atmospheric gases were supplied to the analyzer from the wet box at ∼0.1 L min⁻¹. In addition to the analysis of equilibrator headspace and atmospheric gases, four standard gases of known mixing ratio (150, 349, 449, and 850 ppm; Praxair) were also plumbed to provide gas flow to the GO8050. Praxair standard gases were evaluated by calibrating a LI840A using the Praxair gases, and then the calibrated analyzer was used to measure the CO₂ content of a World Meteorological Organization (WMO) traceable standard gas cylinder from the National Oceanic and Atmospheric Administration Earth System Research Laboratories (ESRL) Greenhouse Gas Global Reference Network. The Praxair gas standard calibrated LI840 was able to reproduce the certified ESRL standard to within 0.1 %.

The GO8050 was controlled using National Instruments LabVIEW software run on a PC laptop computer. The software controls data acquisition from the GO8050, an interface box connecting the SBE 38 and the SBE 45, the Aanderaa 4330F, the primary equilibrator temperature and pressure sensors, a Vaisala digital barometer (0.07 mbar accuracy) with a model 61002 Gill pressure port and GPS antenna positioned adjacent to the atmospheric air intake, and the LI840A, while also controlling a Valco Instruments Co. Inc. (VICI) multi-port actuator that cycles between the gas streams plumbed to the dry box. The software captured measurements from all ancillary sensors as well as analyses of the four gas standards of known CO₂ content, 12 measurements of atmospheric CO₂, and 240 seawater CO₂ measurements in a cycle that was repeated every 8.5 h with a 2 min measurement frequency. The seawater and atmospheric CO₂ analyses were run in a sequence of three atmospheric measurements and 60 seawater measurements that was repeated four times between standardization. Analyses of each gas standard were interpolated over the time record of the dataset and used to create calibration functions needed to correct the raw LI840A xCO₂ data. Calibrated seawater xCO₂ data in dry air were quality-controlled and then converted to CO₂ partial pressure (pCO₂) in wet air at the equilibrator temperature by using atmospheric pressure measured by the LI840A plus the differential pressure recorded in the equilibrator corrected for the removal of water vapor. Finally, seawater pCO₂ in wet air was adjusted to sea surface temperature using the offset between SBE 45 temperature recorded at the GO8050 and intake temperature from the SBE 38 located at the seawater intake with zero lag (0.3 ^∘C ± 0.17 ^∘C), as the lag between the temperature measurements at these two locations was determined to be less than the measurement frequency.

Total uncertainty in our pCO₂ measurements is the combined uncertainties from calibration, equilibrator temperature, equilibrator pressure, and the warming correction added in quadrature. At contemporary atmospheric pCO₂ levels near 400 µatm, these component uncertainties would equate to 0.4, 0.06, 0.17, and 3.9 µatm. Typically, underway pCO₂ measurement uncertainties are reported as a function of uncertainties in the equilibrator temperature and pressure and the water vapor pressure (Wanninkhof et al., 2013). Considering our dried gas stream that minimizes uncertainty from water vapor pressure, the uncertainty from just the equilibrator temperature and pressure, with inclusion of the calibration uncertainty, would equal 0.44 µatm. However, taking into account uncertainty in the warming correction (while still not addressing deviations from a constant pCO₂ temperature sensitivity) increases the pCO₂ uncertainty to 3.92 µatm. While we prefer the more typical assessment that points to a lower pCO₂ uncertainty, we use a conservative ∼1 % pCO₂ uncertainty below to estimate the uncertainties in derived marine CO₂ system parameters.

In June 2019, a BioGeoChemical SUrface MOnitoring-system (BGC-SUMO) was configured with the GO8050 to provide underway pH measurements on the total hydrogen ion scale (pH_T). The BGC-SUMO measures pH_T, temperature, and nitrate concentration, although the latter measurement was not successful on this vessel. The pH_T was measured using a Deep-Sea DuraFET, consisting of an ion-sensitive field effect transistor (ISFET) and a chloride ion-selective electrode as the reference (Johnson et al., 2016). Thus, seawater is unmodified, and no chemicals are added as it flows through the BGC-SUMO. The pH_T sensor was calibrated prior to deployment on the M/V Columbia, and its performance was verified based on discrete samples taken alongside the sensor (n=9) throughout the deployment (Takeshita et al., 2018). Based on this comparison, we assume an uncertainty in pH_T of 0.01. Maintenance on all instrumentation configured aboard the M/V Columbia was conducted during service stops in Ketchikan to prevent biofouling.

3.2 Discrete sample collection

Discrete seawater samples were collected on two ferry trips in November 2017 and August 2018. Samples were drawn from the seawater supply line immediately upstream of the GO8050 into rinsed 350 mL amber soda-lime glass bottles and analyzed for TCO₂ and pCO₂ within a month of collection following methods described elsewhere (Evans et al., 2019). Briefly, TCO₂ and pCO₂ were analyzed from the same sample bottle in this order at the Hakai Institute's Quadra Island Field Station using a Burke-o-Lator pCO₂–TCO₂ analyzer. The TCO₂ measurement was achieved by acidification and gas stripping followed by non-dispersive infrared detection using a LI-COR LI840A and consumed ∼60 mL of sample. TCO₂ measurements were adjusted using correction factors developed through the analysis of certified reference materials (CRMs) from Andrew Dickson (Scripps Institute of Oceanography), with typical correction factors between 0.99 and 1.01. Uncertainty in the discrete TCO₂ measurement was determined to be 0.3 % (Evans et al., 2019). The pCO₂ measurement was achieved by headspace gas recirculation between the LI840A and the sample bottle in a closed loop until equilibration of the headspace gas with the seawater sample pCO₂ was obtained (roughly 6 min). Uncertainty in the discrete pCO₂ measurement was determined to be 1 %. The TCO₂ measurement was subsequently headspace gas corrected (Wanninkhof and Thoning, 1993), and then alkalinity (Alk) was computed using the pCO₂ and head space gas-corrected TCO₂ data with a MATLAB version of CO2SYS (Sharp et al., 2021) and the carbonic acid dissociation constants of Waters et al. (2014), bisulfate dissociation constant of Dickson et al. (1990), fluoride and hydrogen association constants from Perez and Fraga (1987), and boron $/$ chlorinity ratio of Uppström (1974). Alk computed in this way excluded contributions from organic acids, phosphate, and silicate.

3.3 Calculations

3.3.1 Gap filling, marine CO₂ calculations, and gridding

The record of underway measurements from the M/V Columbia contained a number of data gaps related to service interruptions, the largest of which was between October 2018 and March 2019 when the ferry went into dry dock for the winter. However, there was a period from 25 August to 2 October 2019 when only the direct pCO₂ measurements were compromised due to an issue with the LI-COR. Subsequently, pCO₂ was estimated indirectly using pH_T measurements and a regional Alk–salinity relationship (Evans et al., 2015). To fill these missing data, pH_T measurements were interpolated to the measurement time of the GO8050. pCO₂ was then computed using the time-matched pH_T data with the relationships described above and derived Alk. Missing measured pCO₂ observations in late 2019 were filled with the computed values.

Seawater pH_T, [H⁺]_T, and Ω_arag were computed for the entire dataset using the salinity, intake temperature, the gap-filled pCO₂ record, and Alk derived from salinity (Evans et al., 2015) with the dissociation constants and relationships described above using a MATLAB version of CO2SYS (Sharp et al., 2021). Uncertainty in pH_T, [H⁺]_T, and Ω_arag derived from our pCO₂ record coupled with salinity-based Alk determinations was assessed using the error propagation routine from Orr et al. (2018) updated in the most recent MATLAB version of CO2SYS (Sharp et al., 2021). Combined standard uncertainties for pH_T, [H⁺]_T, and Ω_arag were computed using the previously described 1 % pCO₂ uncertainty, the reported 17.21 µmol kg⁻¹ uncertainty in the regional Alk–salinity relationship (Evans et al., 2015), and the default uncertainties for the dissociation constants within the error propagation routine. pH_T, [H⁺]_T, and Ω_arag uncertainties were computed across the range of observed Alk values and with pCO₂ computed across a range TCO₂:Alk ratios spanning 0.85 to 1 for each corresponding Alk value (Fig. S1). These calculations were done at a constant temperature and with salinity ranging from 10 to 32 corresponding with the range of Alk values. Importantly, marine CO₂ system data quality falls into two objectives as defined by the Global Ocean Acidification Observing Network (Newton et al., 2015; Tilbrook et al., 2019): (1) weather and (2) climate. The weather data quality objective is thought sufficient for identifying spatial and temporal patterns excluding long-term trends, which is considered more appropriate for data reaching the stringent climate quality objective to assess. The mean pH_T, [H⁺]_T, and Ω_arag uncertainties from our calculations are 0.01, 0.23 nmol kg⁻¹, and 0.07 respectively. These values meet the Global Ocean Acidification Observing Network weather data quality objective. However, we note that uncertainties vary across the range of values considered. For instance, Ω_arag uncertainty is higher at higher Ω_arag values, whereas pH_T uncertainty is higher at lower salinity and Alk (Fig. S1).

To evaluate basic statistics along the M/V Columbia transit, including means, coefficients of variation (CV), and lower 5th percentiles, as well as assess seasonal drivers and the time of detection that are both described below, observations were gridded by isolating and averaging data within 0.03^∘ by 0.03^∘ grid cells. It is important to note that due to the 2018–2019 winter data gap, gridded averages likely over-represent spring and summer relative to the autumn and winter conditions. This grid size equalled roughly 6 km² across the latitudinal range of the M/V Columbia transit. Analyses using gridded data were only conducted on grid cells containing more than 40 measurements.

3.3.2 Seasonal drivers

Assessing the drivers of seasonal pCO₂ variations requires isolating the thermodynamic and biophysical components of the variability. The process to achieve this is described in Sect. S2, and results in isolating the pCO₂ temperature component (pCO₂ T component), the pCO₂ salinity component (pCO₂ S component), and remaining variability from biophysical drivers. Seasonal amplitudes of each component of pCO₂ variability were assessed, and the ratio of the amplitude of thermodynamic (T, S, or combined TS) to biophysical drivers (B_T, B_S, B_TS; where subscript denotes the removed terms) defines which is more important for determining pCO₂ variability on an annual basis (Takahashi et al., 2002; Fassbender et al., 2018b).

3.3.3 Severity and time of detection

We determined the severity of derived pH_T and Ω_arag in each grid cell based on the lower 5th percentile as in Chan et al. (2017) and the timing of severe conditions as the mode of all months of observations less than or equal to the lower 5th percentile of each grid cell. We also assessed the time of detection (ToD) within each grid cell of the M/V Columbia transit in order to guide future observational efforts targeting the identification of long-term change. ToD is similar to the time of emergence used in climate studies (Henson et al., 2017) with the exception that it includes measurement uncertainty (Carter et al., 2019b). Both of these terms represent the time required for a secular trend, in our case increasing seawater pCO₂ from anthropogenic CO₂ uptake, to emerge from the “noise” in an environmental dataset. Monthly mean pCO₂ is computed from the observations occurring within each grid cell, and then the observations are differenced from the monthly mean in order to compute de-seasonalized anomalies (i.e., removing the large-amplitude seasonal cycle from the noise). The standard deviation of the de-seasonalized anomalies was combined in quadrature with the pCO₂ measurement uncertainty to represent the remaining environmental noise and compute ToD as

$$\begin{array}{}\text{(1)}& \text{ToD}={\displaystyle \frac{\mathrm{2}\times \text{noise}}{p{\mathrm{CO}}_{\mathrm{2}}\text{growth rate}}},\end{array}$$

where the pCO₂ growth rate used here was 2.5 µatm yr⁻¹ and is approximately the average of annual values from the National Oceanic and Atmospheric Administration ESRL over the 2014–2019 period. Importantly, we consider ToD as a guiding metric. The growth of seawater pCO₂ can vary across coastal settings and may or may not be entirely driven by anthropogenic CO₂ input (Laruelle et al., 2018; Salisbury and Jönsson, 2018). For example, changes in nutrient input from runoff can alter the pCO₂ growth rate from an expected anthropogenic CO₂-driven signal (Turk et al., 2019). Therefore, we present ToD only to discuss how these data might be used to target observing efforts and not as absolute values.

3.3.4 Anthropogenic CO₂

Anthropogenic CO₂ content was estimated using the ΔTCO₂ approach (Takeshita et al., 2015; Pacella et al., 2018; Evans et al., 2019), which is a simplification of the ΔC^∗ method (Sabine et al., 2002; Gruber et al., 1996), and assumes a constant TCO₂ disequilibrium with the atmosphere defined as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathrm{\Delta }{\mathrm{TCO}}_{\mathrm{2},\mathrm{diseq}}={\mathrm{TCO}}_{\mathrm{2},\mathrm{obs}}-{\mathrm{TCO}}_{\mathrm{2}}\\ & \left(\mathrm{atm}p{\mathrm{CO}}_{\mathrm{2}},\mathrm{current}\mathrm{year}-\text{age},{\mathrm{Alk}}_{\mathrm{der}},T_{\mathrm{obs}},S_{\mathrm{obs}}\right),\end{array}\end{array}$$

where TCO₂,obs is the observed TCO₂ and ${\mathrm{TCO}}_{\mathrm{2}}\left(\mathrm{atm}p{\mathrm{CO}}_{\mathrm{2}},\text{current year}-\text{age},{\text{Alk}}_{\mathrm{der}},T_{\mathrm{obs}},S_{\mathrm{obs}}\right)$ is the TCO₂ content that would result from equilibration with the atmospheric pCO₂ at the time of last contact with the atmosphere (current year minus the age of the water mass), at the derived Alk and at the observed temperature and salinity. The time of last contact with the atmosphere represents the age of a water mass in years and is zero for most surface measurements except in areas where deep water is mixed to the surface. Water mass age was estimated by dividing the measured apparent oxygen utilization (AOU, or the inverse of ΔO₂) by the oxygen utilization rate (OUR). A value of 4.1 µmol kg⁻¹ yr⁻¹ for OUR was taken from the literature for Pacific subarctic upper water (Feely et al., 2004b) and used in this calculation. Using the ΔTCO_2,diseq term and assuming patterns in derived Alk and observed temperature and salinity are largely invariant, the TCO₂ for a given year can be estimated by

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\mathrm{TCO}}_{\mathrm{2}},\mathrm{year}={\mathrm{TCO}}_{\mathrm{2}}\\ & \left(\mathrm{atm}p{\mathrm{CO}}_{\mathrm{2},\mathrm{year}}-\text{age},{\mathrm{Alk}}_{\mathrm{der}},T_{\mathrm{obs}},S_{\mathrm{obs}}\right)\\ & +\mathrm{\Delta }{\mathrm{TCO}}_{\mathrm{2},\mathrm{diseq}},\end{array}\end{array}$$

where TCO₂,year is the TCO₂ content for a specific year and ${\mathrm{TCO}}_{\mathrm{2}}\left(\mathrm{atm}p{\mathrm{CO}}_{\mathrm{2},\mathrm{year}}-\text{age},{\mathrm{Alk}}_{\mathrm{der}},T_{\mathrm{obs}},S_{\mathrm{obs}}\right)$ is the TCO₂ content that would be realized if that surface water mass were in equilibrium with the atmospheric pCO₂ that occurred during a given year, corrected for the age of the water mass, and at the contemporary derived Alk and observed temperature and salinity. The anthropogenic CO₂ content is then determined as the difference between the TCO₂ for a given year and the TCO₂ content for the year 1765. Historical atmospheric CO₂ mol fractions based on observations and projected atmospheric CO₂ for the shared socio-economic pathways (SSPs) were obtained from Meinshausen et al. (2020) using their data portal (http://greenhousegases.science.unimelb.edu.au, last access: May 2020) and converted to pCO₂ assuming standard atmospheric pressure. pH_T, [H⁺]_T, and Ω_arag were computed for each year from 1765 onward using the TCO₂ estimated for a given year with the modern derived Alk and observed temperature and salinity. It is important to acknowledge that uncertainty in estimating anthropogenic CO₂ content using this approach is at least 5 µmol kg⁻¹ based on similarities with the ΔC^∗ method (Sabine et al., 2002). Uncertainty stems from a number of sources, including the key assumptions of constant TCO₂ disequilibria and unchanged variation in the natural carbon cycle, temperature, and salinity. Inadequacies in these assumptions can lead to biases in anthropogenic CO₂ (Matsumoto and Gruber, 2005), which in turn influences estimations of past and future pH_T, [H⁺]_T, and Ω_arag.

4 Results and discussion

4.1 Time and space variability

Over 244 000 seawater temperature, salinity, O₂, and pCO₂ measurements were made on the M/V Columbia during 135 north- and south-bound transits of the Inside Passage over a 2-year period. These data revealed substantial spatiotemporal variability in surface seawater conditions along this 1600 km stretch of coastline. The spatial and temporal mosaic captured by these measurements (Figs. 2 and 3) portrays two key features of the Inside Passage: (1) the dominant mode of temporal variability is the seasonal cycle, and (2) there is regional variability in the seasonal cycle amplitude that is modulated by the relative influences of tidal mixing, net community production, and the magnitude and character of freshwater input.

Between November and March, cold seawater spanned the entire Inside Passage, with the coldest water in southeast AK generally near 4 ^∘C, but 0.5 ^∘C was observed near Juneau. Seasonal warming in most regions began in April and occurred earlier in the Salish Sea, which was consistent with satellite observations that identified earlier seasonal warming in this region relative to coastal areas to the north (Jackson et al., 2015). Surface salinity was fresher throughout the year in the Salish Sea, although variability in salinity was larger in southeast AK (Figs. 1 and 2) where seasonal freshwater delivery to the coastal ocean contributes 41 % of the freshwater input to the Gulf of Alaska (Edwards et al., 2020). The combined estimates of discharge from each major watershed along the Inside Passage from Edwards et al. (2020) and Morrison et al. (2012) indicate that over 570 km³ yr⁻¹ of freshwater enters the northeast Pacific from southeast AK, an amount that exceeds the Mississippi River discharge (Dai and Trenberth, 2002). Along the BC portion of the Inside Passage, discharge is near 390 km³ yr⁻¹ with almost a quarter of this amount originating from the Fraser River. Despite lower runoff from BC, its influence on salinity manifests earlier than the peak freshwater input in southeast AK (Fig. 2) due to the high contribution of snowmelt to the late spring and early summer discharge (Morrison et al., 2012). In southeast AK, seasonal reduction in salinity began in May and reached the summer minima in August. Low-salinity conditions were uniform over an area that encompassed Lynn Canal and the inside waterways around Juneau (Figs. 1 and 2). The late summer minimum in salinity reflects the significant contribution of glacial melt in the most northern portion of the Inside Passage (Neal et al., 2010; Edwards et al., 2020). Seasonal variation in temperature and salinity was reduced in some confined waterways, such as Johnstone Strait (Fig. 1), owing to the influence of intense tidal mixing in these areas that dampens the seasonal cycle amplitude (Dosser et al., 2021; Whitney et al., 2005). Despite the emergence of a marine heat wave in the North Pacific in late 2019 (Amaya et al., 2020), spring and summer patterns in temperature and salinity appeared similar between 2018 and 2019.

Figure 2 Sea surface temperature (SST; ^∘C; a) and salinity (b) measured between 3 November 2017 and 2 October 2019. The x and y axes represent longitude and latitude, respectively, and with the coastline and terminal positions shown as in Fig. 1 and time increasing along the z axis.

Across most of the Inside Passage, ΔO₂ and pCO₂ showed an inverse relationship (Fig. 3). Where ΔO₂ values were positive, pCO₂ was undersaturated with respect to the atmosphere, and this combination likely reflects the influence of primary productivity exceeding rates of organic matter respiration, i.e., positive net community production (NCP). Abiotic changes in ΔO₂ result from changes in temperature and salinity as well as bubble injection and wave breaking (Juranek et al., 2019), although we contend that the latter two drivers may be of lesser importance in the protected Inside Passage waterways. Seasonal warming increases both ΔO₂ and pCO₂. However, we observed an increase in ΔO₂ with a corresponding decrease in pCO₂, which strongly suggests that O₂ supersaturation and pCO₂ drawdown resulted from positive NCP (Tortell et al., 2012; Juranek et al., 2019). In areas outside of the influence of tidal mixing, the signals of O₂ supersaturation and pCO₂ drawdown were initiated in response to the spring phytoplankton bloom, and generally were sustained through summer until the autumn storm season commenced (Evans et al., 2019; Fassbender et al., 2018a). An exception was Lynn Canal in southeast Alaska (Fig. 1) where the relationship between O₂ and pCO₂ diverged in summer (Fig. 3) when the seasonal change in salinity was maximal (Fig. 2). The addition of cold glacial meltwater results in undersaturated surface pCO₂ (Cai et al., 2021; Pilcher et al., 2016; Evans et al., 2014) while also increasing oxygen solubility and subsequently decreasing ΔO₂ (Fig. 3). The diverging character between O₂ and pCO₂ within Lynn Canal in summer dissipated during autumn when salinity increased in response to storm-induced vertical mixing. Autumn marked the transition back to supersaturated pCO₂ with respect to the atmosphere throughout the Inside Passage. Inter-annual variability was apparent in this dataset during the spring and summer months, as 2019 had slightly greater O₂ supersaturation and pCO₂ drawdown during spring in the Salish Sea and on the central BC coast and throughout much of the summer in southeast AK (Fig. 3).

Figure 3 ΔO₂ (µmol kg⁻¹; a) and pCO₂ (µatm; b) measured between 3 November 2017 and 2 October 2019. The x and y axes represent longitude and latitude, respectively, and with the coastline and terminal positions shown as in Fig. 1 and time increasing along the z axis.

Seasonal O₂ supersaturation and pCO₂ drawdown do not occur uniformly along the Inside Passage but in distinct regions separated by areas of tidal mixing that support sustained low-O₂ and high-pCO₂ conditions (Fig. 3) due to the near-continuous ventilation of sub-surface waters (Whitney et al., 2005; Dosser et al., 2021; Evans et al., 2012; Tortell et al., 2012). The most obvious region of intense tidal mixing along this coastline was in Johnstone Strait between Vancouver Island and mainland BC (Fig. 1), but other areas were also evident, including in the narrow waterway north of Sitka known as Sergius Narrows (Fig. 1). As mentioned above, the seasonal amplitude in the tidal mixing zones is reduced because the water column may be completely mixed, and seasonal variation in these areas may more reflect that of sub-surface water entering the mixing zone laterally (Dosser et al., 2021). Since seasonality in these areas is potentially more influenced by sub-surface source waters, the seasonal cycle can be out of phase with adjacent areas outside of the tidal mixing zones. This was most obvious in Johnstone Strait, where high-pCO₂ conditions outside of this area were generally seen during winter, whereas within this region, the highest pCO₂ was in autumn. The highest observed seawater pCO₂ was near 1200 µatm in Johnstone Strait during September. Winter pCO₂ values outside of the tidal mixing zones broadly ranged between 450 and 800 µatm, being higher in regions with less direct connection to the open continental shelf, such as in the semi-enclosed Salish Sea, in areas of the central BC coast, and in southeast AK (Fig. 3). These areas receive high amounts of riverine organic matter (St. Pierre et al., 2021; Oliver et al., 2017; Johannessen et al., 2003) that may be confined to the nearshore zone by winter-time downwelling circulation (Thomson, 1981; Weingartner et al., 2009) and there subsequently remineralized by the microbial community (St. Pierre et al., 2020), leading to elevated surface pCO₂ nearshore that is not seen in the offshore data along this coast (Evans and Mathis, 2013). In tidal mixing zones like Johnstone Strait, the highest pCO₂ in early autumn decreased through winter to a minimum by late spring, albeit with values that were still supersaturated with respect to atmospheric pCO₂. This difference in timing likely reflects the seasonality of sub-surface waters (Dosser et al., 2021), since, without a short residence time (Pawlowicz et al., 2007), these waters would experience a build-up of respiratory CO₂ through the growing season as organic matter rains out of the surface layer and is respired at depth by the microbial community. We suggest this sub-surface respiration signal is ventilated in the tidal mixing zones and is responsible for the early autumn peak in surface pCO₂.

Seawater pH_T and Ω_arag variability was evaluated by employing an Alk–salinity relationship developed from observations spanning a large portion of the region (Evans et al., 2015). Validation of this relationship was done using Alk determined from seawater samples collected during ferry ride-along cruises and processed as described above. These cruises occurred in November 2017 and August 2018 and spanned the dynamic range of observed salinity conditions (Fig. 2). During November, comparison between discrete Alk and salinity-derived Alk was within 2 times the root-mean-square error of the salinity-based relationship (Fig. S2). During August, larger divergence between discrete and salinity-derived Alk occurred in low-salinity water within the northernmost portion of the Inside Passage. Specifically in the area of Lynn Canal, Alk determined from the salinity-based relationship overpredicted bottle-determined Alk by at most 200 µmol kg⁻¹. Salinity-based Alk determination was further evaluated in 2019 by comparing estimated pH_T, computed from directly measured pCO₂ and salinity-based Alk, to directly measured pH_T (Fig. S3). pH_T was measured from June to October 2019 over the period of lowest observed salinities in southeast AK and revealed a similar pattern to the discrete Alk comparison. Divergence between estimated and directly measured pH_T was greatest in seawater with salinity <22 and north of 57^∘ N in the region around Juneau and up Lynn Canal (Fig. S3). In the analysis that follows, we continue to use salinity-based Alk with our gap-filled pCO₂ record to determine components of the marine CO₂ system along the Inside Passage, but we acknowledge that the northernmost region during summer is likely more corrosive for aragonite than our analysis suggests because of our Alk over-predictions in low-salinity water. It is also possible that local deviations from the assumption of proportionality between salinity and calcium made within CO2SYS may counteract a portion of this “missing” corrosive signal in the low-salinity surface water of southeast AK (Beckwith et al., 2019). The presence of proton-binding dissolved organic molecules may cause additional complication by impacting the interpretation of low-salinity Alk measurements used to generate the regional Alk–salinity relationship (Sharp and Byrne, 2020). While the magnitude of how corrosive Ω_arag is during the melt season in the northernmost area of the Inside Passage may be less well-constrained, the drivers and timing of adverse conditions should not deviate from what we describe below. The confounding factors of variable freshwater Alk, interpretations of Alk measurements in the presence of proton-binding dissolved organic molecules, and potential for non-zero calcium end-members in this region all demand further study in order to more accurately assess the magnitude of corrosive summer Ω_arag conditions in these glacial-melt-influenced waters.

Figure 4 pH_T (total scale; a) and Ω_arag (b) derived from measurements made from 3 November 2017 to 2 October 2019. The x and y axes represent longitude and latitude, respectively, and with the coastline and terminal positions shown as in Fig. 1 and time increasing along the z axis.

Patterns in seawater pH_T and Ω_arag were largely the inverse of that for pCO₂ (Fig. 4); areas exhibiting pCO₂ undersaturation with respect to the atmosphere typically co-occurred with high-pH_T and high-Ω_arag conditions, whereas regions with high pCO₂ have low pH_T and Ω_arag. Areas with both high pH_T and Ω_arag have experienced recent positive NCP that would also support O₂ supersaturation and pCO₂ drawdown (Figs. 3 and 4). The evidence of inter-annual variability discussed above for spring and summer O₂ and pCO₂ was apparent for pH_T and Ω_arag, with 2019 exhibiting more frequent occurrences of Ω_arag>3 compared to 2018. pH_T and Ω_arag were lowest in most areas during winter and year-round within tidal mixing zones. Winter Ω_arag values were <1 in all regions that lacked direct connection to the open continental shelf, specifically within the Salish Sea, Johnstone Strait, inside waterways on the central and northern BC coast, and in southeast AK. Corrosive conditions for aragonite persisted throughout the year in Johnstone Strait, and in an area known as Wrangell Narrows between Wrangell and Petersburg (Fig. 1). In the northernmost area of the Inside Passage, a short period of Ω_arag conditions >1 occurred between March and June, resulting from the spring phytoplankton bloom as evidenced by co-occurring O₂ supersaturation and pCO₂ drawdown (Fig. 3). Once the summer melt season commenced, the Inside Passage-wide minimum in Ω_arag was observed in this region despite the over-prediction in Alk in low-salinity water mentioned above. Lynn Canal exhibited the most corrosive conditions for aragonite along the 1600 km M/V Columbia transit due to the large contribution of meltwater in this region (Fig. 2). Such corrosive conditions in glacial-melt-influenced settings have been reported previously in AK (Reisdorph and Mathis, 2013; Evans et al., 2014) as well as in Svalbard (Ericson et al., 2019; Cantoni et al., 2020). Co-occurring corrosive conditions for aragonite (Fig. 4) and undersaturated pCO₂ with respect to the atmosphere (Fig. 3) are unique to cold glacial-melt-influenced coastal regions, which likely enables a positive feedback whereby CO₂ influx from the atmosphere either enhances or prolongs corrosive summer Ω_arag conditions (Evans et al., 2014; Ericson et al., 2019; Cantoni et al., 2020).

4.2 Seasonal drivers

pCO₂ variability is determined by thermodynamic and biophysical forcings, the latter being the sum of the physical and biogeochemical influences of vertical mixing, horizontal transport, NCP, sea–air CO₂ exchange, and calcification. Seasonal variation in pCO₂ reflects the interaction of these terms, which often are competing. For instance, warming and freshwater input have opposing influences on CO₂ solubility such that together they can dampen pCO₂ variability (Cai et al., 2021; Salisbury and Jönsson, 2018). As described (Sect. S2), the ratio of the seasonal amplitude of pCO₂ at T_mean (B; Fig. S4) to the pCO₂ T component (T; Fig. S5) provides information on whether biophysical processes or seasonal warming are more important for shaping pCO₂ variability within a region. Takahashi et al. (2002) describe this as $T/B$ (or as a difference, T–B; see their Fig. 9), where if $T/B$ (or T–B) is greater than 1 (or positive), seasonal temperature change is the dominant process determining pCO₂ variability. The global analysis by Takahashi et al. (2002) suggests that in the area closest to BC and southeast AK, temperature and biophysical processes play equal roles in determining pCO₂. A more recent analysis by Fassbender et al. (2018b) produced similar results for the northeast Pacific with balanced roles of temperature and biophysical processes evident closest to the coast. However, both of these analyses were conducted with large global grids that did not resolve the coastal margin and did not differentiate the role of freshwater given the open-ocean focus.

In the nearshore zone spanning BC and southeast AK, it is essential to account for salinity variation when assessing pCO₂ variability. As pointed out by Sarmiento and Gruber (2006), variations in pCO₂ that result from changes in salinity cannot be evaluated based solely on the salinity sensitivity (Takahashi et al., 1993) because this only accounts for changes in solubility and not the corresponding change in TCO₂ and Alk from a decrease in salinity. Instead, the contribution of changes in salinity to the pCO₂ variability can be evaluated by incorporating TCO₂ and Alk buffer factors into the calculation (Sect. S2). Changes in salinity would result from both freshwater input (decrease) and vertical mixing (increase) and are expressed here as the pCO₂ S component (S; Fig. S6). As mentioned, the pCO₂ S component and pCO₂ T component can be in opposition such that their corresponding influences on pCO₂ are counterbalanced (Fig. S7). However, there are times and locations when these factors are not balanced. Lynn Canal (Fig. 1) during the summer months is an important example of an area and time period when salinity variability exceeds the influence of seasonal warming (Fig. 5). Subtracting both the pCO₂ S component and the pCO₂ T component (TS) from the observed pCO₂ leaves remaining variability associated with NCP, calcification, and gas exchange (Fig. S8). Given that calcification is only episodically important in this region and gas exchange is generally slow (on the order of months), this remaining pCO₂ variability largely reflects the influence of NCP, or CO₂ removal and addition by organic matter production and degradation.

Figure 5 Ratios of the seasonal amplitudes of thermodynamic and biophysical drivers of pCO₂ variability. (a) The ratio of the amplitude of the temperature component (T) to the amplitude of the remaining biophysical components (B_T). (b) The ratio of the amplitude of the salinity component (S) to the amplitude of the remaining biophysical components (B_S). Areas with major river outflows are highlighted in this panel (Taku River (T), Stikine River (St), Nass River (N), Skeena River (Sk), and the Fraser River (F)). (c) The ratio of the amplitude of the combined temperature and salinity components (TS) to the amplitude of the remaining biophysical components (B_TS).

As illustrated in Fig. 5, the biophysical component dominates over the temperature component in shaping pCO₂ variability on an annual basis everywhere along the Inside Passage. Excluding Lynn Canal, the salinity component is also less important than the biophysical component, even in areas adjacent to major river outflows. At the outflows of major rivers, such as the Fraser and Stikine (Fig. 5), the salinity component is an important contributor but still roughly 30 % less than the amplitude of the biophysical component. In Lynn Canal, salinity variance exceeded all other locations along the Inside Passage (Figs. 1 and 2), which resulted in a dominant contribution to the pCO₂ variability (Fig. 5). Temperature counterbalanced some of the salinity component in Lynn Canal, such that this was the only area where there was near equivalence between the combined thermodynamic components and the biophysical drivers in determining pCO₂. These computations show the spatial complexity in the balance between thermodynamic and biophysical drivers in the nearshore zone and that the influence of salinity must be considered with temperature in settings with significant freshwater input. The importance of salinity in shaping marine CO₂ system variability in this region has been discussed previously in modeling studies by Siedlecki et al. (2017) and Hauri et al. (2020), as well as by Pilcher et al. (2016), who evaluated the role of variability in freshwater Alk end-members in enhancing nearshore atmospheric CO₂ uptake. However, the contributions of thermodynamic versus biophysical drivers to the observed variability have not been evaluated to the extent shown here, which indicated the dominance of biophysical drivers over most of the Inside Passage.

The northernmost reach of the Inside Passage is heavily influenced by changes in salinity resulting from the volume of glacial melt water entering this area (Neal et al., 2010; Edwards et al., 2020). Reisdorph and Mathis (2013) first described the influence of meltwater on marine CO₂ chemistry in this region, and subsequent observational and modeling work has considered the de-coupling that can occur between pCO₂ and Ω_arag in locales of significant cold glacial melt discharge (Evans et al., 2014; Ericson et al., 2019; Hauri et al., 2020; Cantoni et al., 2020). Given that atmospheric CO₂ uptake is promoted in glacially influenced regions, these areas may be important amplifiers of OA (Cantoni et al., 2020; Ericson et al., 2019; Evans et al., 2014). Increasing glacial discharge, changes in glacial meltwater Alk as glaciers further recede and the flow path over land to the ocean increases, increasing glacial river temperatures, and increasing organic matter decomposition in glacial rivers are all factors that would modulate the extremely corrosive conditions found within these nearshore environments, as well as the decoupling between pCO₂ and Ω_arag. Given the potential for intensifying positive feedback with further increasing atmospheric pCO₂ and enhanced sea–air CO₂ exchange, thereby amplifying the already extreme Ω_arag conditions, additional research effort should target these areas in order to understand which feedbacks dominate from seasonal to inter-annual timescales.

4.3  Characterizing regional extremes

Identifying regional extremes in the marine CO₂ system is important for characterizing environmental variability, identifying where unfavorable conditions for vulnerable marine species occur more often or more intensely, and pin-pointing areas that may experience faster rates of change from anthropogenic CO₂ input (Feely et al., 2018; Hare et al., 2020). Marine CO₂ system extremes were characterized here based on pH_T and Ω_arag variability and severity (Chan et al., 2017). Importantly, these two descriptors for extremes may not manifest the same way in a region or with the same timing; a region may have severely low pH_T but also low variability and experience severely low pH_T at a different time of year than in adjacent areas. pH_T and Ω_arag extremes can also be temporally mismatched within a region.

Figure 6 Panels (a–c) show the pH_T coefficient of variation (CV), severity, and the timing of severe pH_T conditions, respectively. Panels (d–f) show the same parameters for Ω_arag.

Extremes, as defined by the variability, were regionally similar for pH_T and Ω_arag; both the Salish Sea and select areas in southeast AK exhibited large variability relative to other areas along the Inside Passage (Fig. 6). Contrasting these highly variable areas, Johnstone Strait, Sergius Narrows, and Wrangell Narrows (Fig. 1) all had low variability owing to the influence of persistent tidal mixing. Using severity to portray extremes provided a nearly inverse picture, with Johnstone Strait exhibiting both severe pH_T and Ω_arag, while Lynn Canal had severe Ω_arag values but less severe pH_T. Notably, Ω_arag severity was only above 1 in surface water most exposed to the open northeast Pacific between AK and BC, an area known as Dixon Entrance (Fig. 1), and in Sitka; areas exposed to the open continental shelf generally had less severe and variable pH_T and Ω_arag compared to more confined waters along the Inside Passage.

There were also differences in the timing of severe pH_T and Ω_arag conditions across and within regions (Fig. 6). The majority of Inside Passage waters experienced severe pH_T and Ω_arag between November and February when seawater pCO₂ was highest; however, severe conditions occurred in some areas earlier in the year. In Johnstone Strait, most severe pH_T and Ω_arag occurred in September, whereas most severe pH_T and Ω_arag conditions occurred in June in areas proximal to the Skeena and Stikine outflows due to the influence of the snowmelt freshet (Fig. 5). In Lynn Canal, most severe Ω_arag values were in August coinciding with the peak input of glacial melt and in November for pH_T owing to storm-induced vertical mixing. The variation in timing of most severe conditions along the Inside Passage may have regionally distinct biological implications when these coincide with times when more sensitive life stages of vulnerable species are present. In addition, the different characterization of extremes based on variability and severity has potential implications for adaptation trajectories as, for example, vulnerable organisms in Johnstone Strait would experience a sustained corrosive and moderately stable low-pH_T environment, whereas in much of southeast AK, vulnerable organisms would be subjected to large swings in marine CO₂ system parameters over the year. While some research considers long-term exposure to variable marine CO₂ conditions to be a factor enhancing physiological tolerance to OA (Kapsenberg and Cyronak, 2019), other research suggests that organisms living in persistently low-pH environments might be more locally adapted (Chan et al., 2017). Here we provide information on the locations of both of these types of settings such that future work can move to examine how species fare along this gradient within the Inside Passage.

Figure 7 Standard deviation of pCO₂ anomalies (µatm; a) and the time of detection (years; b) to resolve the secular trend of increasing seawater pCO₂ that tracks the contemporary rise in atmospheric CO₂.

Characterizing extremes is also useful for guiding observing efforts by identifying locations that either minimize the amount of time anticipated to observe an OA-driven change in the marine CO₂ system or capture a key process that may be driving large-amplitude signals in a region (Turk et al., 2019). Areas with low natural variability require fewer years to resolve a secular trend as opposed to regions of high variability (Sutton et al., 2019). Figure 7 shows the standard deviation of de-seasonalized pCO₂ observations (anomalies), representing the environmental noise along the Inside Passage, with the resulting ToD (computed following Eq. 1). The Salish Sea, the area near the Stikine River, and the northern portion of Lynn Canal all have very long ToD, and observing efforts in these regions would be better suited for targeting the processes discussed previously that shape the variability. As these processes may themselves be subject to climate change (Bidlack et al., 2021), tracking their evolving influence on marine CO₂ system parameters will provide valuable information on the dynamics organisms are subjected to in these highly variable environments. On the other hand, Johnstone Strait, the area south of Ketchikan, and Sergius Narrows all have much shorter ToD (Fig. 7). These areas would be ideal for placing observing assets aimed at resolving long-term secular trends. It is also important to understand that the ToD estimates computed here are “forced” values (Turk et al., 2019), in that they are based on seawater pCO₂ increasing at a similar pace to the present atmospheric CO₂ increase (Sutton et al., 2019). In some cases, there may be an observed trend in a time series that differs from the forced trend (Laruelle et al., 2018), and this can reflect either the role of other processes independent of anthropogenic CO₂ increase that modulates the trend in seawater pCO₂ (Turk et al., 2019; Salisbury and Jönsson, 2018) or faster increases in pCO₂ resulting from anthropogenic CO₂ addition in weakly buffered settings (Feely et al., 2018). Areas, such as Johnstone Strait, with low variability that results in short ToD, and with severe pH_T and Ω_arag stemming from higher TCO₂:Alk ratios and weaker buffering, likely will exhibit faster rates of change than estimated by forced ToD values. Identifying and establishing these areas as sentinel sites for tracking OA would optimize coastal observing efforts aimed at resolving long-term secular trends.

4.4 Estimating past and future conditions

We consider below how the marine CO₂ system along the Inside Passage has evolved over the industrial era as well as what additional change might be anticipated if greenhouse gas emissions are reduced to reach the preferable Paris Agreement target of 1.5 ^∘C warming (UNFCC, 2015); referred to here as the 1.5 ^∘C acidification level. It is important to note that our evaluation is theoretical and only considers the role of anthropogenic CO₂ and not the influence of other forcings like increasing temperature or changing freshwater input. Following the approach outlined above, anthropogenic CO₂ for Inside Passage surface waters was determined and showed notable spatiotemporal variability (Fig. 8). This variability was strongly influenced by freshwater input and water mass age, of which the latter ranged from 0 to 35 years (Fig. S9). Maximal age estimates were confined to the areas of persistent tidal mixing and were similar to estimates from other studies for the age of upwelled water present on the northeast Pacific continental shelf (Feely et al., 2008; Murray et al., 2015). This agreement is encouraging considering the potential for inaccurately representing OUR in the calculation due to a likely higher oxygen utilization in the confined nearshore regions (Johannessen et al., 2014; Pawlowicz et al., 2007). It is therefore worth considering how inaccurately estimating water mass age translates to uncertainty in anthropogenic CO₂ and characterization of preindustrial [H⁺]_T, pH_T, and Ω_arag. Inaccurate water mass age estimates are more influential in older water masses at the surface and over time periods when atmospheric pCO₂ is changing faster (i.e., when computing contemporary ΔTCO_2,diseq). Keeping this in mind, we consider how a 50 % uncertainty in the age of surface water in Johnstone Strait impacts our estimation of preindustrial conditions. An overestimate of the age by 50 % results in a lower ΔTCO_2,diseq by 11.2 µmol kg⁻¹ and a higher anthropogenic CO₂ by 11.6 µmol kg⁻¹. This overestimate adjusts the contemporary [H⁺]_T, pH_T, and Ω_arag acidification signals (i.e., the difference between contemporary and preindustrial values) by 0.96 nmol kg⁻¹, −0.03, and −0.08, respectively. An underestimate of the age by 50 % would increase ΔTCO_2,diseq by 9.6 µmol kg⁻¹ and decrease anthropogenic CO₂ by 10 µmol kg⁻¹. The contemporary [H⁺]_T, pH_T, and Ω_arag acidification signals would adjust by −1.4 nmol kg⁻¹, 0.04, and 0.07, respectively. Despite the presence of variability in the age estimate for Johnstone Strait (Fig. S9), we suggest it is unlikely that the age is underestimated. Rather, this variability is more likely a function of seasonal variation in sub-surface water mass oxygen utilization (Johannessen et al., 2014) and the application of a constant OUR; although, given that the presence of older, upwelled water is seasonal along this coastline (Feely et al., 2016), some variability in the age of surface water masses within persistent tidal mixing zones is expected. We therefore use the absolute values of the shifts in [H⁺]_T, pH_T, and Ω_arag resulting from an overestimate of the water mass age by 50 % as uncertainty bounds when considering the contemporary and 1.5 ^∘C acidification levels.

Figure 8 Contemporary anthropogenic CO₂ content (µmol kg⁻¹; a) and the estimated first year when Ω_arag was <1 (b). The x and y axes represent longitude and latitude, respectively, and with the coastline and terminal positions shown as in Fig. 1 and time increasing along the z axis.

Interestingly, areas identified previously as pH_T and Ω_arag extrema based on severity, due to either persistent tidal mixing or glacial melt input, were not locales containing the highest anthropogenic CO₂ content. Instead, the highest values were in regions that experienced the greatest O₂ supersaturation and pCO₂ drawdown during summer (Figs. 3 and 8). The highest estimated values were near 66 µmol kg⁻¹ and similar to other estimates for coastal northeast Pacific surface water (Carter et al., 2019a; Feely et al., 2016); however, much lower values were evident in some locations. The Salish Sea is considered to have more moderate anthropogenic CO₂ levels (Feely et al., 2010; Hare et al., 2020; Evans et al., 2019), but the freshest areas observed in southeast AK exhibited very low anthropogenic CO₂ content. Such weakly buffered areas are likely locations that have been corrosive for aragonite prior to the industrial era. To consider this possibility, we estimated the first year when Ω_arag was <1 along the Inside Passage (Fig. 8) by calculating the anthropogenic CO₂ content accrued each year along with the resulting change in Ω_arag over the industrial era. The most weakly buffered areas, either due to being very fresh or because of ventilation of sub-surface water with high TCO₂ relative to Alk, are shown by this calculation to have been corrosive at least on a seasonal basis since the start of the industrial era (Fig. 8). These naturally corrosive hot spots are being amplified by anthropogenic CO₂ addition such that, for example, Johnstone Strait now experiences under-saturation throughout the year (Fig. 4). The shift to corrosive winter conditions has occurred over more recent decades outside of the mixing zones. Winter surface water in the Salish Sea likely transitioned to Ω_arag<1 beginning around 1950, consistent with the emergence of corrosive winter values found in a previous study (Evans et al., 2019). Similar winter transition timing was evident for Inside Passage waters on the central BC coast, although it has appeared more recently within the last few decades over large portions of southeast AK (Fig. 8).

Figure 9 Contemporary [H⁺]_T values minus the estimated values for 1765 (nmol kg⁻¹; a) and the estimated values for 2035 minus the 1765 values (nmol kg⁻¹; b). The x and y axes represent longitude and latitude, respectively, and with the coastline and terminal positions shown as in Fig. 1 and time increasing along the z axis.

Matthews et al. (2021) determined the allowable future CO₂ emissions that limit global warming to the preferred 1.5 ^∘C level stated in the Paris Agreement (UNFCC, 2015), the so-called “remaining carbon budget”, to be 440 GtCO₂ from 2020 onwards. Using the relationships described in Friedlingstein et al. (2021), full emission of the remaining carbon budget can be equated to a rise in atmospheric CO₂. With 3664 GtCO₂ equalling 1 GtC, and every 2124 GtC of emissions increasing atmospheric CO₂ by 1 ppm, the remaining carbon budget would drive an increase in atmospheric CO₂ of 56 ppm. At 1 atmosphere of pressure, combining the contemporary atmospheric CO₂ mole fraction with the remaining carbon-budget-forced atmospheric CO₂ increase results in an atmospheric pCO₂ of 468 µatm. The time at which this atmospheric pCO₂ would be realized is trajectory-dependent, with both the sustainable development pathway (SSP1) and the fossil-fuel development pathway (SSP5) reaching this atmospheric pCO₂ level at different times (Meinshausen et al., 2020). SSP5 reaches, and surpasses, this value quickly, by roughly 2035. SSP1 takes longer and reaches this level by 2063. Here we use this atmospheric CO₂ target to consider the theoretical 1.5 ^∘C acidification level and report the year 2035 as the fastest trajectory following the central estimate for the year when 1.5 ^∘C warming would be reached if the current rate of warming continues (IPCC, 2018). However, what is important is that this is an anticipated extent of OA reachable with either trajectory. If our society follows SSP1, this is theoretically the most acidification we should expect without considering amplifying processes. However, if we follow a SSP5-type scenario, acidification will surpass what we estimate below.

Figure 10 The contemporary Ω_arag values minus the estimated values for 1765 (a) and the estimated values for 2035 minus the 1765 values (b). The x and y axes represent longitude and latitude, respectively, and with the coastline and terminal positions shown as in Fig. 1 and time increasing along the z axis.

The estimated average percentage increase in [H⁺]_T along the Inside Passage is 40 ± 18 % over the 254 years since the start of the industrial era, although important spatiotemporal variability in acidification was evident (Fig. 9). This average extent of acidification equates to a 0.14 ± 0.02 unit drop in pH_T. The Inside Passage average pH_T change exceeds the global average (Lauvset et al., 2020; Jiang et al., 2019); however, similar to [H⁺]_T, pH_T change varied spatially and temporally with values ranging from −0.06 to −0.20 (Fig. S10). Spatially, [H⁺]_T and pH_T change has been greatest in the more weakly buffered and moderately anthropogenic CO₂-concentrated waters of the Salish Sea and Johnstone Strait. The largest [H⁺]_T change was evident in Johnstone Strait (Fig. 9), with a maximum near 7 nmol kg⁻¹, due to the inherently lower background pH_T level (Fassbender et al., 2021). For comparison, open-ocean surface water has experienced an average [H⁺]_T change of 1.6 nmol kg⁻¹ (Fassbender et al., 2021). The pH_T change within these two settings appeared only marginally different, being 0.17 versus 0.1, but with an over 4-fold increase in [H⁺]_T in Johnstone Strait. On a seasonal basis over most of the Inside Passage, and despite a higher anthropogenic CO₂ signal during summer (Fig. 8), the change in [H⁺]_T and pH_T appeared larger in winter. Note that seasonal variation in the acidification signal exceeded the 0.96 nmol kg⁻¹ and 0.03 [H⁺]_T and pH_T uncertainties, respectively. This seasonality in acidification manifests because of seasonal differences in TCO₂:Alk ratios (Fig. S11) that alter the CO₂ system response to anthropogenic CO₂ increase. During winter, the TCO₂:Alk ratio is closer to unity such that seawater is more weakly buffered and the percent change in pH following a percentage change in TCO₂ has a greater magnitude. Larger percent changes in pCO₂ and Ω_arag are also expected during this season, and this pattern follows the modeled (Fassbender et al., 2018b; Kwiatkowski and Orr, 2018) and observed (Landschützer et al., 2018) changes in the seasonality of surface marine CO₂ parameters at the global scale.

At an atmospheric pCO₂ of 468 µatm, an additional 17 ± 7 % increase in [H⁺]_T, on average, would theoretically be expected for the 1.5 ^∘C acidification level. This implies that nearly half of the acidification experienced thus far over the industrial era will likely occur over the coming 15 years if society maintains the current emissions trajectory. However, a change in emissions trajectory that follows a sustainable development pathway would enable this acidification signal to occur over a longer period of time. It is anticipated that acidification will be further amplified during winter along the Inside Passage and particularly within the semi-enclosed and more weakly buffered waterways (Fig. 8). Johnstone Strait and the Salish Sea will likely continue to experience the largest changes in [H⁺]_T (Fig. 9), and these areas may serve as bellwethers for the emergence of biological OA impacts in a similar manner as how high-latitude settings are viewed (Fabry et al., 2009). Efforts to examine biological impacts in situ should target these regions where we estimate the largest contemporary and 1.5 ^∘C acidification levels. In addition, studies challenging organisms to adverse [H⁺]_T and pH_T levels within experimental settings may benefit from our estimates of the 1.5 ^∘C acidification level, as these could serve as a near-term treatment for diagnosing OA impacts.

Unlike for pH_T, change in Ω_arag over the industrial era has a seasonal maximum in summer (Fig. 10). This characteristic can also be explained by considering the seasonality in seawater TCO₂ and Alk as well as relative versus absolute changes in Ω_arag. Despite the TCO₂:Alk ratio being closer to unity during winter, and the percent change in Ω_arag following a percentage change in TCO₂ being larger during that season (Fig. S11), summer Ω_arag values are much higher than winter values (Fig. 4). A 14 % change in an Ω_arag value near 1 is a smaller absolute change than a 9 % change in an Ω_arag value of 3. Considering the 1.5 ^∘C acidification level, the change in summer Ω_arag values at times may exceed 0.8 units, consistent with an overall reduction in seasonality as anthropogenic CO₂ content continues to increase (Kwiatkowski and Orr, 2018). To our knowledge, differences in the season during which maximum absolute changes in [H⁺]_T and pH_T versus Ω_arag occur have not been widely acknowledged in the literature and point to the need for careful consideration of the specific marine CO₂ system parameter an organism may be most sensitive to (Waldbusser et al., 2014). Seasonally specific changes in the most impactful marine CO₂ system parameter for a sensitive species may or may not align with periods of maximal vulnerability (Hales et al., 2016). Considering how the marine CO₂ system is being modified by increasing anthropogenic CO₂ on a seasonal as well as long-term basis, and which specific variable is most impactful for an organism, are both essential elements for understanding the implications of OA.

5 Conclusions

Through partnership with the Alaska Marine Highway System, we have reduced the information gap on marine CO₂ system variability along the northeast Pacific Inside Passage. This study has documented spatiotemporal variability in surface water along this 1600 km passageway and shown that the dominant mode of temporal variability is the seasonal cycle and that the amplitude of this signal varies spatially and is modulated by the relative influences of tidal mixing, net community production, and the magnitude and character of freshwater input. While this effort advanced our understanding of marine CO₂ system variability in this region, winter observations were limited to a single year, and inter-annual variability was not adequately constrained. Enhancing winter observations and further evaluating the magnitude of inter-annual variations are both important next steps for marine CO₂ system research in this region.

We have highlighted that within the northernmost area of the Inside Passage, deviations in freshwater end-members from the broader regional alkalinity–salinity relationship and the proportionality between salinity and calcium require further study given the potential for positive feedback with atmospheric CO₂ uptake and for modifications in freshwater outflows that can alter the coastal OA signals. The impact of proton-binding organic molecules on the interpretation of alkalinity measurements in this region is also a large unknown. Analysis of seasonal drivers indicated that the biophysical component has a dominant role in shaping variability along most of the Inside Passage but that the combined influences of temperature and salinity balance the biophysical component in glacial-melt-impacted areas of southeast AK where these uncertainties are expected to be the largest. Further research along these lines will be valuable for the research community striving to understand marine CO₂ system patterns in areas exhibiting large-amplitude variation in salinity.

We considered the characterization of pH_T and Ω_arag extremes and recognized that there are regional differences in the manifestation of extremes based on variability versus severity that likely have biological implications. Vulnerable organisms experiencing a sustained corrosive and moderately stable low-pH_T environment may have a differing adaptation trajectory than organisms subjected to large swings in marine CO₂ system parameters over the year. Our diagnosis of these locations should be useful for future studies examining organismal and ecosystem adaptation trajectories within the context of OA. We also used our variability assessment to determine the time of detection and point out that this information can help optimize coastal observing efforts aimed at establishing sentinel sites to resolve long-term secular trends or further evaluating the drivers of large-amplitude variability.

Finally, we estimated the anthropogenic CO₂ content in surface water and considered change over the industrial era and to an atmospheric CO₂ level that corresponds with the exhaustion of the remaining 1.5 ^∘C carbon budget. It was shown that some areas, including the tidally mixed Johnstone Strait and within Lynn Canal, have likely experienced seasonal Ω_arag values <1 over the entire industrial era. Other areas have transitioned to winter Ω_arag values <1 more recently. Seasonal differences were also identified in the absolute changes in [H⁺]_T, pH_T, and Ω_arag. Absolute [H⁺]_T and pH_T changes appeared larger during winter when conditions are more weakly buffered while Ω_arag change was greater during summer. This difference should be a consideration when evaluating biological OA impacts. Looking to the future, the theoretical 1.5 ^∘C acidification level suggests that significant marine CO₂ system changes may develop over the coming 15 years if society continues on a fossil-fuel development emissions trajectory. These estimates of the 1.5 ^∘C acidification level should be useful as a near-term treatment for regional challenge studies aiming to diagnose species responses to OA. Time series observations must be expanded and maintained along the Inside Passage to determine if the 1.5 ^∘C acidification level is realized to the extent suggested by this study.