Domain Editor-in-Chief: Jody W. Deming; University of Washington, WA, US
Associate Editor: Lisa A. Miller; Fisheries & Oceans Canada, CA

1. Introduction

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 (CO2) 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 CO2 concentrations with respect to pre-industrial levels (IPCC, 2013). The amount of CO2 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 CO2 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 CO2, regional factors have the potential to exacerbate acidification in coastal areas. Coastal oceans are affected directly not only from local emissions of CO2 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 CO2 consumption in surface coastal waters during the day, but the same algae produce CO2 during respiration at night. Moreover, as algae die and sink, bacteria decompose the additional organic matter reaching the bottom, releasing CO2 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).

Figure 1 

Maps of Salish Sea with location of observations, stations used for analysis, and model domain. (a) Location of observations from Puget Sound Regional Synthesis Model cruises (PRISM, pink diamonds), Washington State Department of Ecology Marine Monitoring Unit (ECY-MMU, green squares), King County Marine and Sediment Assessment Group (brown circles), and Padilla Bay NOAA National Estuarine Research Reserve System (black triangles). Some regions mentioned in the text are also shown. (b) Salish Sea Model (SSM) domain. The resolution of this unstructured grid ranges from 60 m in narrow inlets to ~3 km near the mouth of the Juan de Fuca Strait. Stars indicate the four stations used for analysis: Gordon Point, West Point, Hood Canal, and Saratoga Passage. DOI:

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.

2. Salish Sea model with FVCOM-ICM

2.1. Model description

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 (NO3), ammonium (NH4+), phosphate (PO43–), 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.

Figure 2 

Biogeochemical model diagram. Schematic diagram of variables and processes represented in the model (total alkalinity and dissolved oxygen not shown). LDOM and RDOM stand for labile and refractory dissolved organic matter, respectively, and LPOM and RPOM, for labile and refractory particulate organic matter. DOI:

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 NH4+ or consume NO3, while consumption of NH4+ and production of NO3 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 NH4+; it decreases due to regenerated primary production, water column nitrification, and sediment fluxes of NO3 (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 CO2, which can either increase or decrease DIC depending on whether the surface partial pressure of CO2 (pCO2) is lower or higher than the atmospheric pCO2 (the latter is an input to the model). Therefore, we also introduced the calculation of the seawater pCO2 and carbonate system constants such as the solubility of CO2 (K0) using the same equations and approach as in the CO2SYS 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 K1, K2 and KSO4 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 K1 and K2, the KSO4 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 CO2 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).

2.2. Model setup

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 pCO2 concentrations, at this time we used 400 μatm for lack of 2008 data in the Salish Sea.

2.3. Approach to assess model performance

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 R2), 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 R2 should be high (a perfect match of model and observations would provide R2 = 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, NO3, NH4+, pH and saturation state of aragonite, ΩA (the last two were derived from the TA and DIC observations using CO2SYS). 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, NO3, and NH4+ for 2008 (Figure 1a, green squares, red circles, and black triangles, respectively).

2.4. Description of sensitivity experiments

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.

Table 1

Description of model experiments. DOI:

Experiment name Description

Baseline 2008 simulation (setup described in Section 2.2)
Freshwater at Ambient Concentrations (or FW@AmbientC) Concentrations of biogeochemical variables at freshwater sources (rivers and point sources) equivalent to ambient seawater concentrations near the discharge areas; performed by including freshwater sources in the hydrodynamics but not in the offline biogeochemical run
High DICobc DIC at Strait of Juan de Fuca open boundary condition increased by 40 mmol m–3 (~2% increase with respect to Baseline)
High DICfreshwater DIC in all freshwater sources increased by 2% with respect to Baseline

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 DICobc”). 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 DICobc 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 DICobc). This “High DICfreshwater” 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 DICobc and High DICfreshwater allowed for the analysis of the response to a similar relative change in DIC in the open boundary vs. local freshwater sources.

3. Results and discussion

3.1. Model performance

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, R2 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 R2 = 0.37 and 0.25, respectively (both still significant at 99%). One of the reasons that led to the low R2 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 NO3 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.

Figure 3 

Model vs. observations for the whole Salish Sea Model domain for 2008. Scatter plots of model vs. observations for the whole water column and all stations shown in Figure 1a for a) Temperature (°C), b) Salinity, c) dissolved oxygen (mg L–1), d) nitrate (mg L–1), e) Chlorophyll (μg L–1), f) Dissolved inorganic carbon (μmol kg–1), g) Total alkalinity (μmol kg–1), and h) pH (total scale). Dark solid lines indicate the least-mean-squares regression slope (see R2 values and other metrics in Table 2) and thin dashed lines indicate the 1:1 slope. DOI:

Table 2

Statistical metricsa of model global performance for year 2008. DOI:

Variable (unit) R2 b RMSE Bias N

T (°C) 0.81 1.48 1.28 67858
Salinity 0.37 1.33 –0.68 66934
Chl (μg L–1) 0.25 2.78 –0.30 66041
DO (mg L–1) 0.64 1.80 –1.56 66538
NO3 (mg L–1) 0.64 0.08 –0.001 1902
DIC (μmol kg–1) 0.59 70.33 –20.13 593
TA (μmol kg–1) 0.66 60.89 –38.75 589
pH (total scale) 0.41 0.14 –0.07 584

a While R2 is dimensionless, RMSE and Bias have the units shown in the first column; N is the total number of comparisons between observations and model.

b All correlations are significant beyond the 1% level.

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.

Figure 4 

Time series of modeled and observed values for the surface layer at four stations. Time series of surface model values (black lines) for seven variables (rows) and four stations (columns). Overlapped observations belong to CTD casts (blue triangles) and bottles (red circles) for 2008; diamonds and crosses correspond to other years in the panels for DIC (μmol kg–1), TA (μmol kg–1), pH (total scale), and ΩA (blue diamonds for PRISM, pink crosses for ECY-MMUdata). DOI:

Figure 5 

Time series of modeled and observed values for the bottom layer at four stations. Time series of bottom model values (black lines) for seven variables (rows) and four stations (columns). Overlapped observations belong to CTD casts (blue triangles) and bottles (red circles) for 2008; blue diamonds correspond to other years of PRISM data in the panels for DIC (μmol kg–1), TA (μmol kg–1), pH (total scale), and ΩA. Scales of the Y-axes are the same as in Figure 4. DOI:

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).

Figure 6 

Depth profiles of modeled and observed values at four stations for 2008. Comparison of modeled values (filled black and gray symbols) and observations from PRISM (open red and blue symbols) in February (circles) and August (triangles). Each column represents a different station and each row, a different variable: DIC (μmol kg–1), TA (μmol kg–1), pH, and ΩA. DOI:

3.2. Sensitivity to freshwater sources

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.

Figure 7 

Time series of sensitivity experiments for the surface layer at four stations. Time series of surface values at four stations (columns) and five variables (rows) for two experiments: Baseline (black) and Freshwater at Ambient Concentrations (FW@AmbientC, red). Experiments are shown for the first year (spin-up time) so that all simulations start with identical values. The five variables shown are chlorophyll (μg L–1), DIC (μmol kg–1), TA (μmol kg–1), pH (total scale) and ΩA. All carbonate-system variables showed a strong response to freshwater loadings at seawater ambient concentrations. DOI:

Figure 8 

Time series of sensitivity experiments for the bottom layer at four stations. Time series of bottom values at four stations (columns) and five variables (rows) for three experiments: Baseline (black), Freshwater at Ambient Concentrations (FW@AmbientC, red), and High DICobc (gray). Experiments are shown for the first year (spin-up time) so that all simulations start with identical values. The five variables shown are chlorophyll (μg L–1), DIC (μmol kg–1), TA (μmol kg–1), pH (total scale) and ΩA. The FW@AmbientC experiment showed the largest deviations with respect to the Baseline simulation. In the High DICobc experiment, changes in DIC by the end of the first year are approximately half of the changes imposed at the open boundary, and pH changes are even closer (~0.08 vs. 0.1 at the open boundary). DOI:

Figure 9 

Mean annual ΔpH (experiment minus baseline) for the bottom layer. Maps of mean annual differences of pH (using the second year of simulation) for the bottom layer of the water column between three experiments and the baseline: (a) Freshwater at Ambient Concentrations (FW@AmbientC), (b) High DICobc, and (c) High DICfreshwater. Color scales of ΔpH are different (up to one order of magnitude) for each panel. In (a), loading from rivers affects all Puget Sound. In (b), the effect of the open boundary is also felt in the whole region, but with more limited influence in shallower areas with river inputs. In (c), changes are only seen in areas close to large rivers. DOI:

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.

Figure 10 

Relative differences between the Freshwater at Ambient Concentrations and Baseline experiments for DO and TA. Relative difference (%) of DO (top row) and TA (bottom row) between the Freshwater at Ambient Concentrations and Baseline simulations at four stations in the first year of the simulation. Positive (negative) values indicate that the former experiment showed higher (lower) values than the baseline. Blue and gray lines indicate the time series at the surface and bottom layers of the model, respectively. The dashed line highlights the 0% change. The relative difference in TA is larger than that of DO in the surface layer, indicating a more important role of freshwater loading for the former variable. DOI:

3.3. Sensitivity to ocean open boundary conditions

The time series of the High DICobc 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 DICobc 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.

3.4. Open boundary vs. freshwater sources

The annual mean changes in bottom pH for the High DICobc with respect to the baseline simulation can be compared with those of the High DICfreshwater (Figure 9b and c, respectively). Increasing DIC at the freshwater sources by 2% (High DICfreshwater) had a negligible effect on pH in the deeper, main basins of Puget Sound compared with the 2% increase at the open boundary (High DICobc). 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 DICobc scenario (Figure 9b).

We note that the comparison of the effects of the High DICobc and High DICfreshwater 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.

4. Conclusions

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 DICfreshwater 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.

Data Accessibility Statement

Model results from this study are available upon request (

Supplemental Files

The supplemental files for this article can be found as follows: