Sensitivity of the regional ocean acidification and carbonate system in Puget Sound to ocean and freshwater inputs

Since the industrial revolution, human activities such as fossil fuel combustion, cement production, and deforestation have resulted in the release of greenhouse gases to the atmosphere. Accumulated anthropogenic emissions of carbon dioxide (CO₂) since 1750 add up to 555 ± 85 gigatons of carbon (GtC), of which 240 ± 10 GtC have accumulated in the atmosphere and led to an approximate 40% increase in current atmospheric CO₂ concentrations with respect to pre-industrial levels (IPCC, 2013). The amount of CO₂ sequestrated by the ocean (155 ± 30 GtC) has led to the lowering of global surface ocean pH from 8.2 to 8.1 (about a 30% increase in acidity) since the beginning of the industrial era (Orr et al., 2005; Feely et al., 2009; Rhein et al., 2013). This gradual acidification of seawater due to increasing atmospheric CO₂ concentrations has been referred to as “ocean acidification” (Broecker and Clark, 2001; Caldeira and Wickett, 2003).

While ocean acidification is a global problem driven by the global increase in atmospheric CO₂, regional factors have the potential to exacerbate acidification in coastal areas. Coastal oceans are affected directly not only from local emissions of CO₂ and other greenhouse gases but also from freshwater sources. For instance, nutrients and organic carbon from land runoff can lead to an enhanced cycle of algae-growth/algae-mortality. The increase in photosynthesizing algae intensifies oxygen production and CO₂ consumption in surface coastal waters during the day, but the same algae produce CO₂ during respiration at night. Moreover, as algae die and sink, bacteria decompose the additional organic matter reaching the bottom, releasing CO₂ and consuming oxygen (sometimes, even generating hypoxic conditions). Furthermore, riverine pH and total alkalinity (TA) are typically lower than seawater values, such that rivers can dilute TA and decrease pH in their areas of influence. In the Puget Sound basin, river pH values range from approximately 6.5 to 8.5 and river TA, from about 200 to 1300 μmol kg^–1 or typically about 10% to 40% of seawater TA (Hallock, 2009), due to minerals leached from soils and/or the decomposition of organic matter (such as plant material, freshwater algae, sewage effluents, etc.). Note that rivers typically also have lower concentrations of dissolved inorganic carbon (DIC) than seawater. While low river DIC tends to increase pH, the low river TA decreases pH and overpowers the effect of DIC in the rivers of Puget Sound. In addition, anthropogenic freshwater point sources such as storm water systems, industrial outfalls, and wastewater treatment plants may contribute nutrients and organic matter to near-shore areas or deeper waters, affecting acidification and the carbonate system (DIC, TA, pH, etc.) locally, especially in poorly flushed areas.

Lastly, coastal regions influenced by wind-driven upwelling are also vulnerable to the pH decline in the open ocean, because deeper ocean waters often enter the shelf during upwelling conditions (Feely et al., 2008; Allen and Durrieu de Madron, 2009). This is the case of the Pacific Northwest coastal waters, where several monitoring programs indicate low and declining pH and related carbonate-system parameters (Wootton and Pfister, 2012; Feely et al., 2012a; Murray et al., 2015). Consequences of this coastal acidification were already observed between 2005 and 2009, when significant production failures at the Pacific Northwest oyster larvae hatcheries were linked to acidic (corrosive) seawater arriving at the West Coast (Washington State Blue Ribbon Panel on Ocean Acidification, 2012; Feely, 2012b). Seawater with low pH also has low aragonite saturation state (Ω_A); aragonite is the least stable form of calcium carbonate formed by many shell-building organisms and its saturation state indicates whether aragonite tends to precipitate (Ω_A > 1) or dissolve (Ω_A < 1). Therefore, the upwelling of waters with lower pH and Ω_A negatively affected the ability of oyster larvae to keep their protective carbonate shells, contributing to the drastic collapse in hatchery larvae between 2005 and 2009.

The relative importance of the different drivers of acidification in the coastal ocean (e.g., open ocean, freshwater sources, local atmospheric inputs) varies by location and is not completely understood. For instance, while upwelling-driven acidification is the key driver of coastal acidification along the outer coast of the Pacific Northwest, estuaries and inland marine waters may also be strongly influenced by freshwater sources. Here, we aim to improve the understanding of the drivers of the carbonate system in Puget Sound, a fjord system that lies within Washington State and is subject to both ocean inputs from upwelled seawater and freshwater sources. Puget Sound is the southern part of the Salish Sea, which also comprises the Strait of Juan de Fuca and the Strait of Georgia (Figure 1a). Following an estuarine circulation, Pacific Ocean waters over the continental shelf enter the Strait of Juan de Fuca through the bottom layers (in particular, upwelled ocean waters of deeper origin entrain in summer and fall) and, while undergoing some mixing, flow both north and south towards the deeper waters of the Strait of Georgia and Puget Sound, respectively. The state of the carbonate system and ocean acidification in Puget Sound is of concern, as more than 30% of its marine species depend on the ability to calcify in order to make shells, skeletons, and other hard body parts (Washington State Blue Ribbon Panel on Ocean Acidification, 2012). Farmed and wild shellfish fisheries are a significant component of the regional economy; e.g., annual sales of farmed shellfish in Washington State (the country’s leading producer of farmed oysters, clams, and mussels) accrued almost $150 million in 2013 (U.S. Department of Agriculture, 2014). Tribal communities of the region depend on shellfish for commercial, cultural, ceremonial, and subsistence purposes. While observations of the carbonate system variables in Puget Sound started quite recently (approximately last decade), we know that acidification conditions vary strongly from place to place and across seasons. In winter, observations showed well mixed, corrosive waters (low pH and Ω_A), while stratification in summer confined corrosive waters to deeper subsurface areas, with the lowest pH and Ω_A values found in Hood Canal (Feely et al., 2010; Reum et al., 2014).

Our specific focus in this study was to investigate the relative role of ocean and freshwater flows in the carbonate system and regional ocean acidification in Puget Sound. For this purpose, we used a coupled physical-biogeochemical numerical model of the Salish Sea that represents the carbonate system as well as other water quality variables. The model allowed us to perform sensitivity analyses to study the importance of oceanic and freshwater forcing in the carbonate system (DIC, TA, pH, Ω_A, etc.). We found that the open ocean plays an important role in defining the carbonate conditions in Puget Sound, particularly at depth, which was expected, given a previous study on dissolved oxygen in the Salish Sea (Roberts et al., 2014). However, the open ocean was not the only key driver; we also found that the rivers exert a major control in the region at all depths, due to the dilution of DIC and TA they cause in this estuarine system.

While our current focus is on Puget Sound, the model extends to the Straits of Juan de Fuca and Georgia (Figure 1b); therefore, we refer to the model as the Salish Sea Model (SSM). The hydrodynamic component of the SSM is an application of the unstructured grid Finite-Volume Community Ocean Model (FVCOM; Chen et al., 2003). The unstructured grid framework allows for the representation of complex shoreline geometry, waterways, and islands in the Salish Sea. The biogeochemical component is an adaptation of the Integrated Compartment Model (CE-QUAL-ICM; Cerco and Cole, 1993, 1994), which we refer to as FVCOM-ICM. Crucial improvements from the original CE-QUAL-ICM are the addition of a two-layer sediment module (Di Toro, 2001) and the inclusion of carbonate system variables (TA, DIC) in the water column. In total, FVCOM-ICM currently tracks 16 state variables in the water column: diatoms, dinoflagellates, dissolved oxygen (DO), nitrate (NO₃^–), ammonium (NH₄⁺), phosphate (PO₄^3–), labile and refractory dissolved organic carbon (LDOC, RDOC), labile and refractory particulate organic carbon (LPOC, RPOC), labile and refractory dissolved organic nitrogen (LDON, RDON), labile and refractory particulate organic nitrogen (LPON, RPON), total alkalinity (TA) and dissolved inorganic carbon (DIC). A schematic diagram represents the processes included in the model (Figure 2). FVCOM-ICM is run offline and uses hydrodynamic fields previously computed by FVCOM using the same model grid (the offline approach is described in detail in Kim and Khangaonkar, 2012). The SSM has been described in detail (Yang et al., 2010; Khangaonkar et al., 2011, 2012a; Kim and Khangaonkar, 2012), so here we describe the improvements to the previously published version of the model; namely, the addition of a sediment module and the carbonate system.

Sediment flux models range from simple empirical relationships (Soetaert et al., 2000; Fennel et al., 2006; Hetland and DiMarco, 2008) or constant fluxes (Soetaert et al., 2000; Scully, 2010) to complex process representations within one or more sediment layers, each representing a particular chemical environment (Emerson et al., 1984; Di Toro, 2001; Gypens et al., 2008). Two-layer models are often used as a compromise between computational efficiency and depth-resolution, while still providing acceptable level of accuracy (Brady et al., 2013; Testa et al., 2013). We coupled the SSM with the two-layer model of Di Toro (2001), since it is well documented, has been found to accommodate sufficiently complex processes at an acceptable level of accuracy, and has been applied to a wide range of freshwater and marine water systems (e.g., Tillman et al., 2004; Brady et al., 2013; Testa et al., 2013). The two sediment layers represent a surface thin aerobic layer of variable thickness followed below by a thicker (~10 cm) anaerobic layer. Particulate organic carbon and nitrogen (including all forms of particulate organic matter, POM, from phytoplankton and detritus) settling from the water column into the sediments are deposited in the deeper, thicker layer. There, this POM is partitioned into three labilities, from labile to basically inert. Decomposition of POM in the sediment (i.e., diagenesis) produces dissolved forms of carbon and nitrogen in the sediment pore water; these solutes react and are transported between the two layers, returned to the overlying water or released as gases (e.g., methane and nitrogen gas). Furthermore, some of the reactions require DO (e.g., aerobic remineralization, nitrification), which is transferred from the overlying water into the top sediment layer. Other reactions included in the sediment model are denitrification, methane oxidation, and sulfide oxidation. For equations and more details on this model, please refer to Di Toro (2001) and Roberts et al. (2015).

To represent the carbonate system in the water column, TA and DIC were included in the model. In contrast to other state variables, which have units of grams per cubic meter (i.e., g m^–3, which is equivalent to mg L^–1), TA and DIC are modeled in terms of millimole per cubic meter (mmol m^–3). In the case of TA, it increases due to processes that produce NH₄⁺ or consume NO₃^–, while consumption of NH₄⁺ and production of NO₃^– decrease TA. Therefore, TA in the model increases due to new primary production, remineralization of LDON and RDON, water column denitrification, and sediment fluxes of NH₄⁺; it decreases due to regenerated primary production, water column nitrification, and sediment fluxes of NO₃^– (see Appendix A.1).

The processes that consume DIC are primary production and water column nitrification. Production of DIC occurs through remineralization of LDOC and RDOC, water column denitrification, phytoplankton losses by predation, basal metabolism and photorespiration, and sediment fluxes (see Appendix A.2). The model also includes air-sea exchange of CO₂, which can either increase or decrease DIC depending on whether the surface partial pressure of CO₂ (pCO₂) is lower or higher than the atmospheric pCO₂ (the latter is an input to the model). Therefore, we also introduced the calculation of the seawater pCO₂ and carbonate system constants such as the solubility of CO₂ (K₀) using the same equations and approach as in the CO₂SYS software by Lewis and Wallace (1998) and van Heuven et al. (2011). The user needs to provide as input to the model the choice of K₁, K₂ and K_SO4 dissociation constants for carbonic acid and bisulfate ion and the desired formulation of the borate-to-salinity ratio to be used. In our current setup, we used the formulations by Lueker et al. (2000) for K₁ and K₂, the K_SO4 of Dickson (1990) and the borate-to-salinity ratio of Uppström (1974), as used in other recent observational studies (e.g., Takahashi et al., 2014; Fassbender et al., 2017). The gas transfer velocity currently used for both air-sea exchange of CO₂ and DO is from Wanninkhof (2014), but the user is allowed to choose among several options (e.g., Wanninkhof, 1992; Nightingale et al., 2000; Ho et al., 2006).

The model was set up to represent year 2008 based on the DIC and TA observations available to us at the time of developing the pH module. This year had DIC and TA data from two University of Washington (UW) Puget Sound Regional Synthesis Model (PRISM) cruises in two different months (February and August), covering all of Puget Sound and the Strait of Juan de Fuca (Feely et al., 2010; Reum et al., 2014; Figure 1a). Previous to our 2008 simulation, the model had been calibrated to year 2006 (all except carbonate system variables) and all model parameters remained unchanged for the 2008 simulation (see Khangaonkar et al., 2012a) for information on the 2006 calibration before adding the carbonate-system and sediment modules and Table A1 for a subset of model parameters in the current FVCOM-ICM). The model was run for two years, both representing 2008 conditions. We refer to the first year as a “spin-up” time, which is the period needed to lessen the influence of initial conditions on model behavior; the second year is a stable simulation that can be compared against observations.

The initial and boundary conditions for variables other than TA and DIC have been described elsewhere (Khangaonkar et al., 2012a, 2012b), so here we focus on the setup of the new carbonate variables and sediment component. The latter was initialized by assuming the sediment is at steady state with the initial depositional fluxes of POM to the sediment layer (based on initial settling fluxes). Initial conditions for DIC were taken as the mean of the observations for all the available PRISM cruises from 2008 to 2014; i.e., a uniform value was applied to the whole model domain, as was the case for all physical and biological variables. Initial conditions for TA were calculated as a function of initial salinity (S), following the regression TA = 646.7 + 47.7 * S, where TA is in μmol kg^–1 (Fassbender et al., 2017). The open boundary conditions at the Juan de Fuca and Johnstone Straits were also calculated as regressions of S, temperature (T) and DO (D Ianson, personal comm.): DIC = –17.51 * T – 0.95 * DO + 2440.62; TA = 470.13 + 52.85 * S if S < 33.85 and TA = –4932.6 + 212.44 * S if 33.85 ≤ S ≤ 34.65 (where DIC, TA, and DO are given in μmol kg^–1 and T in °C). Rivers and other freshwater sources are also treated as boundary conditions (Khangaonkar et al., 2012a), with both concentrations and flows prescribed. For the rivers, observations of pH and TA from US and Canadian gauged rivers were used to derive regressions as a function of river flow (see Appendix B.1 and B.2); then, DIC was calculated using both TA and pH values. In the case of the other freshwater point sources (e.g., wastewater treatment plants and industrial outfalls, referred to as point sources from now on), we used mean monthly time series of pH and median values of TA based on data from US and Canadian plants (Appendix B.3) and then calculated the corresponding time series for DIC. For atmospheric pCO₂ concentrations, at this time we used 400 μatm for lack of 2008 data in the Salish Sea.

To understand how well the model represents the observations, we used both qualitative and quantitative analyses. For the latter, we calculated statistical metrics that allow measuring model performance in terms of correlation (Pearson’s R²), accuracy (root mean square errors or RMSE) and bias with respect to observed values. In order to conclude that the model is able to represent the observations, the R² should be high (a perfect match of model and observations would provide R² = 1) and RMSE and bias should be low (these metrics would be zero for a perfect match). In addition, we overlapped observations to simulated time series and vertical profiles of different variables and at different locations to assess the ability of the model to simulate the observations both temporally and spatially. The selected stations (Gordon Point, West Point, Hood Canal and Saratoga Passage) were located in four different regions of Puget Sound (South Sound, Central Basin, Hood Canal and Saratoga Passage, respectively; see Figure 1b).

The observations used to evaluate model performance came from different sources. As mentioned at the beginning of Section 2.2, the TA and DIC observations belong to two PRISM cruises in February and August (Figure 1a, magenta diamonds). These cruises also provided other data, such as T, S, DO, NO₃^–, NH₄⁺, pH and saturation state of aragonite, Ω_A (the last two were derived from the TA and DIC observations using CO₂SYS). Furthermore, monthly monitoring efforts by the Washington State Department of Ecology Marine Monitoring Unit (ECY-MMU (Keyzers, 2014; Bos et al., 2015)), the King County Marine and Sediment Assessment Group (Smith, 2003), and the Padilla Bay NOAA National Estuarine Research Reserve System (NERRS, 2015) provided observations of T, S, DO, NO₃^–, and NH₄⁺ for 2008 (Figure 1a, green squares, red circles, and black triangles, respectively).

We compared the baseline simulation against three sensitivity experiments (Table 1). Each of the three experiments had a specific change in the setup to analyze its effect on the model results. In particular, we used the first year (spin-up) of the simulations for the comparisons of model evolution (i.e., time series plots), such that all comparisons would start from the same initial conditions.

In the first experiment, all freshwater sources (rivers and point sources) were allowed to flow with their temperature (T) and salinity (S = 0) signatures in the hydrodynamic model. The biogeochemical model was then run offline after the hydrodynamic model, but without any freshwater sources. The result is a simulation where freshwater sources affected circulation, T, and S, but did not create any dilution of water-quality variables; i.e., it is a hypothetical scenario where freshwater enters the model with nutrients, DO, TA, and DIC at concentrations equal to those of the ambient seawater near the discharge areas. We called this simulation the “Freshwaters at Ambient Concentrations” experiment (or “FW@AmbientC”); it allowed studying the effect of the loadings from rivers and point sources on the acidification state of Puget Sound.

The second experiment studied the effect of changes in the ocean boundary conditions in Puget Sound (“High DIC_obc”). In this simulation, DIC open boundary concentrations at Juan de Fuca Strait were increased by 40 mmol m^–3 (~2% increase), which is approximately equivalent to decreasing the pH by ~0.1 units (consistent with the overall expected global decrease of 0.06–0.32 pH over the 21st century; IPCC, 2013).

In order to compare the role of the open ocean in regulating pH in Puget Sound with that of the freshwater sources, we compared the changes produced by the High DIC_obc simulation with those of a third experiment. This time, DIC in freshwater sources was increased by 2% (analogous to the increase in the open boundaries in the High DIC_obc). This “High DIC_freshwater” experiment led to decreases in river pH ranging from 0.02 to 0.21 (annual mean decrease for all rivers combined was 0.08) and a decrease of 0.22 in point sources. The comparison of the changes by both the High DIC_obc and High DIC_freshwater allowed for the analysis of the response to a similar relative change in DIC in the open boundary vs. local freshwater sources.

We compared observations from Puget Sound and Strait of Juan de Fuca from 2008 against model results in several different ways. First, we compared the whole observational dataset against the corresponding model values in the form of scatter plots of different variables (Figure 3). The statistical metrics resulting from this comparison indicated an acceptable overall performance of the model in these two regions (Table 2). For instance, R² values ranged between ~0.3 and 0.8, and they were all significant at more than 99% confidence. The lowest correlations were for S and chlorophyll (Chl) with R² = 0.37 and 0.25, respectively (both still significant at 99%). One of the reasons that led to the low R² for S was the still-insufficient model resolution in near-shore regions affected by river plumes; sometimes a station belongs to the river plume in the observations but not in the model and vice versa. Similarly, Chl is patchy and highly variable, making it difficult for the model to accurately represent the spatial and temporal variability. The RMSE values indicate that on average the model was within 1.5°C, salinity of 1.3, 2.78 μg Chl L^–1, 1.8 mg DO L^–1, 0.08 mg NO₃ L^–1, 70.3 μmol DIC kg^–1, 60.9 μmol TA kg^–1, and 0.14 pH units of the observations. The biases were slightly negative (i.e., the mean of model predictions underestimated the mean of observations), except for the case of T, which presented a warm bias. Overall, while the differences between model and observations (particularly in S and Chl) suggest that there is room for model improvements (e.g., increase resolution in narrow inlets, intertidal regions, and around islands), the statistical metrics are significant and within reasonable values.

We selected four stations (see Section 2.3 and Figure 1b) to compare the modeled time series at surface and bottom against the observations (Figures 4 and 5). The model successfully represented the seasonal cycle of the observations, also capturing features related to high frequency variability (e.g., tidal mixing, storm systems, variability in freshwater discharge). In particular, both model and observations showed larger temporal variability at surface (Figure 4) than at bottom (Figure 5), due to the higher frequency of surface forcing (e.g., air-sea exchange, storms, runoff, photosynthetic processes) compared with that of bottom forcing (incoming waters from the open ocean, sediment fluxes, remineralization processes). For the carbonate system variables, observations for 2008 were only available in February and August, so observations for other years were added to the time series plots (blue diamonds and pink crosses). While we do not expect a perfect match between our 2008 model and observations belonging to other years, the good agreement indicates that both the seasonal cycle and the magnitudes of the 2008 simulation are reasonable.

Vertical profiles for the same four stations (Figure 6) for February (circles) and August (triangles) 2008 showed an overall agreement between model and observations for the carbonate variables. The model reproduced the observed DIC vertical stratification at all stations and both months, while modeled TA tended to have slightly more vertical structure than the observations in Gordon Point and West Point. The errors in both DIC and TA propagate into the calculations of pH and Ω_A, resulting in larger deviations between observed and modeled pH and Ω_A. Nevertheless, in most cases these modeled variables showed a good representation of the observed vertical structure. The largest relative errors between modeled and observation-based pH and Ω_A were found in summer in Gordon Point, due to the combination of errors in TA (model higher than observations) and DIC (model lower than observations).

The Freshwaters at Ambient Concentrations experiment (red lines in Figures 7 and 8) showed a strong response of the carbonate-system variables during the spin-up year. Time series of DIC, TA, pH, and Ω_A at the four different locations showed marked deviations with respect to the baseline (black lines in Figures 7 and 8). With the loading of TA and DIC from freshwater sources at ambient seawater concentrations, the concentrations of these variables became higher than in the baseline simulation (close to ocean values). Therefore, the dilution of DIC and TA by freshwater sources in Puget Sound is an important process for all carbonate variables, as can also be seen in the pH and Ω_A time series plots (lower panels in Figures 7 and 8). Furthermore, the horizontal distribution of the annual mean differences of bottom pH between the second (stable) year of the Freshwaters at Ambient Concentrations and the baseline simulations (Figure 9a) emphasizes that the influence of freshwater loading is seen far away from discharge areas. Mean annual changes in bottom pH extended all over Puget Sound, with increases from 0.05 at Admiralty Inlet to >2 at the mouth of large rivers (mean increase over Puget Sound ~0.3). Note that DIC and TA have opposing effects on pH and Ω_A; i.e., higher DIC decreases pH and Ω_A, while higher TA increases them. Therefore, while the lack of low DIC and TA from the rivers in the Freshwaters at Ambient Concentrations simulation increased both the ambient DIC and TA, the latter led to the resulting increase in pH and Ω_A. Changes in ambient TA in this experiment were more prominent than changes in DIC because other processes (e.g., photosynthesis, air-sea exchange) still acted to modulate DIC, while they did not control TA as much.

Other biogeochemical variables were also affected in the Freshwaters at Ambient Concentrations experiment. The relative change in DO with respect to the baseline simulation was smaller than the change in TA at surface during the spin-up year, but still of a considerable magnitude (Figure 10). While TA increased up to 70% at the surface of the four stations, DO decreased up to 40%. Furthermore, DO tended to decrease at the surface (blue lines in Figure 10) due to the reduction in primary production (resulting from less nutrients supplied by freshwater sources) and generally increased at the bottom (gray lines in Figure 10) because of the decrease in organic matter flux and remineralization at depth. In contrast, TA increased both at the bottom and the surface, because biological fluxes of TA are more than one order of magnitude smaller than advective fluxes, leading to a relatively more important role of freshwater loadings. Bottom changes in both DO and TA were of comparable magnitude, especially at the end of the first year of the simulation when the circulation in the model was already established. The relative changes in the bottom layers were smaller than at surface layers due to the controlling effect of incoming ocean waters at bottom.

The time series of the High DIC_obc experiment are shown only for the bottom of the water column because that is the region of maximum influence of incoming ocean waters (gray lines in Figure 8). As only DIC was changed in the open boundary, only carbonate-system variables were affected (DIC does not limit primary production in the model). By the end of the first year of the simulation, bottom DIC increased in the four stations by ~18–21 mmol m^–3 (similar values in μmol kg^–1), pH decreased by ~0.07–0.08 units and Ω_A, by up to 0.09. Therefore, changes in the open boundary conditions propagate to most of Puget Sound within a year. The annual mean of bottom pH differences between the second year of both simulations (High DIC_obc minus Baseline; Figure 9b) further emphasize that most of the Central Basin showed an annual mean pH decrease of ~0.1, while Hood Canal and parts of South Sound showed decreases of 0.06 or lower. The shallower areas influenced by river discharges (e.g., inlets in South Sound, Skagit Bay) exhibited smaller decreases in pH due to the controlling effect of freshwater flows.

The annual mean changes in bottom pH for the High DIC_obc with respect to the baseline simulation can be compared with those of the High DIC_freshwater (Figure 9b and c, respectively). Increasing DIC at the freshwater sources by 2% (High DIC_freshwater) had a negligible effect on pH in the deeper, main basins of Puget Sound compared with the 2% increase at the open boundary (High DIC_obc). In the former case (Figure 9c), pH decreased by 0.16 close to the Fraser River and by < 0.1 near the largest Puget Sound rivers; the mean decrease in the bottom waters of Puget Sound was ~0.006. The decreases in Puget Sound main basins (areas in red/orange in Figure 9c) are around one to two orders of magnitude smaller than the changes due to the High DIC_obc scenario (Figure 9b).

We note that the comparison of the effects of the High DIC_obc and High DIC_freshwater does not represent a fair evaluation of the future roles of the ocean and freshwater boundaries, particularly because it is unclear how DIC will change in the freshwater sources in the next century. For instance, changes could go in different ways in different river basins and point sources (increase vs. decrease in DIC). Mechanisms responsible for changes in riverine DIC are also likely to lead to changes in TA, further increasing the uncertainty of the pH signature of rivers in the future. River flow discharges and seasonality are also expected to change along with climate, which in turn will affect the carbonate system by modifying the rates of dilution. Nevertheless, the experiments compared in this section allow evaluating the role of a similar relative increase in DIC in both ocean and freshwater boundaries.

The different sensitivity experiments performed in this study helped improve the understanding of the carbonate system in Puget Sound. Our results suggest that both ocean and freshwater sources are key drivers of DIC, TA, and the related variables, pH and Ω_A. Changes in the ocean boundary affected conditions all over Puget Sound, but more strongly in bottom waters due to the characteristic estuarine circulation. Within a year, the change in DIC imposed to the open boundary arrived at regions such as the Central Basin (West Point station) and South Sound (Gordon Point station) after mixing and dilution due to circulation. Shallower bays and inlets with strong riverine influence were more resilient to changes in the open boundary. We also found that freshwater loading has a far-reaching effect, since a simulation with biogeochemical loading at seawater concentrations (Freshwaters at Ambient Concentrations experiment) generated a strong response in the carbonate system throughout all Puget Sound. This experiment is a useful way to determine the important role that the dilution by freshwater sources exerts in the carbonate system; it is not intended to be a realistic scenario or to determine how biogeochemical loading should be altered in order to ameliorate the effects of ocean acidification in Puget Sound. Small changes in freshwater DIC (2% increase, High DIC_freshwater experiment) did not show major changes in the carbonate system away from the major rivers. However, larger changes in freshwater concentrations have the potential to affect carbonate variables in Puget Sound. We emphasize the opposing effects that DIC and TA exert on variables such as pH and Ω_A, as different imbalances between them in rivers, point sources, and/or oceanic forcing could lead to different responses in pH and Ω_A in Puget Sound. This issue highlights the importance of continuing long-term monitoring of these variables in both freshwater sources and incoming ocean waters.

Future work will focus on discerning between the effects of rivers vs. point sources on the carbonate system; these freshwater sources differ both in levels of discharge and loading, as point sources tend to have much lower flows and much higher DIC and TA than rivers. As part of our current efforts to improve model performance, we are working on increasing horizontal resolution in narrow inlets, intertidal regions, and passageways. These changes will allow improving the representation of steep bathymetry in those areas without incurring pressure gradient errors associated with vertical sigma coordinates. Further efforts are aimed at modeling other years, for which richer carbonate system datasets are now available, and extending the model domain by moving the open boundary to the Washington and Vancouver Island shelves. By placing this boundary farther away from the mouth of the Strait of Juan the Fuca, conditions entering the Strait (which are critical for the bottom waters of the Salish Sea) will be generated dynamically by the model (e.g., through wind-driven upwelling), instead of depending on user input. Despite current limitations of the model, the results presented here emphasize the important roles of both open ocean conditions and freshwater inflows in the carbonate system of Puget Sound; the former was an expected result, but we believe that the latter role has been previously overlooked.

Model results from this study are available upon request ([email protected]).

We acknowledge fruitful interactions with Carl Cerco during the early stages of the pH module development.

Funding was provided by EPA Grant Agreements No. PC-00J20101 and PC00J89901, Catalog of Domestic Assistance Number 66.123-Puget Sound Action Agenda: Technical Investigations and Implementation Assistance Program.

The authors have no competing interests to declare.

• Contributed to conception and design: LB

• Contributed to acquisition of data: LB, WL, GP, AA, TM, MR

• Contributed to analysis and interpretation of data: LB, WL, TK, GP, AA

• Drafted and/or revised the article: LB, WL, TK, GP, AA, TM, CFK

• Approved the submitted version for publication: LB, WL, TK, GP, AA, TM, MR, CFK