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

While ocean acidification was first investigated as a global phenomenon, coastal acidification has received significant attention in recent years, as its impacts have been felt by different socio-economic sectors (e.g., high mortality of shellfish larvae in aquaculture farms). As a region that connects land and ocean, the Salish Sea (consisting of Puget Sound and the Straits of Juan de Fuca and Georgia) receives inputs from many different sources (rivers, wastewater treatment plants, industrial waste treatment facilities, etc.), making these coastal waters vulnerable to acidification. Moreover, the lowering of pH in the Northeast Pacific Ocean also affects the Salish Sea, as more acidic waters get transported into the bottom waters of the straits and estuaries. Here, we use a numerical ocean model of the Salish Sea to improve our understanding of the carbonate system in Puget Sound; in particular, we studied the sensitivity of carbonate variables (e.g., dissolved inorganic carbon, total alkalinity, pH, saturation state of aragonite) to ocean and freshwater inputs. The model is an updated version of our FVCOM-ICM framework, with new carbonate-system and sediment modules. Sensitivity experiments altering concentrations at the open boundaries and freshwater sources indicate that not only ocean conditions entering the Strait of Juan de Fuca, but also the dilution of carbonate variables by freshwater sources, are key drivers of the carbonate system in Puget Sound.


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

RESEARCH ARTICLE
Sensitivity of the regional ocean acidification and carbonate system in Puget Sound to ocean and freshwater inputs Laura Bianucci *, ‡ , Wen Long * , Tarang Khangaonkar * , Gregory Pelletier † , Anise Ahmed † , Teizeen Mohamedali † , Mindy Roberts †, § and Cristiana Figueroa-Kaminsky † While ocean acidification was first investigated as a global phenomenon, coastal acidification has received significant attention in recent years, as its impacts have been felt by different socio-economic sectors (e.g., high mortality of shellfish larvae in aquaculture farms).As a region that connects land and ocean, the Salish Sea (consisting of Puget Sound and the Straits of Juan de Fuca and Georgia) receives inputs from many different sources (rivers, wastewater treatment plants, industrial waste treatment facilities, etc.), making these coastal waters vulnerable to acidification.Moreover, the lowering of pH in the Northeast Pacific Ocean also affects the Salish Sea, as more acidic waters get transported into the bottom waters of the straits and estuaries.Here, we use a numerical ocean model of the Salish Sea to improve our understanding of the carbonate system in Puget Sound; in particular, we studied the sensitivity of carbonate variables (e.g., dissolved inorganic carbon, total alkalinity, pH, saturation state of aragonite) to ocean and freshwater inputs.The model is an updated version of our FVCOM-ICM framework, with new carbonate-system and sediment modules.Sensitivity experiments altering concentrations at the open boundaries and freshwater sources indicate that not only ocean conditions entering the Strait of Juan de Fuca, but also the dilution of carbonate variables by freshwater sources, are key drivers of the carbonate system in Puget Sound.

Keywords:
Ocean acidification; Carbonate system dynamics; Physical-biogeochemical model; Ocean numerical model; Salish Sea; Puget Sound 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 nearshore 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 upwellingdriven 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.

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 andCole, 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 3 -), ammonium (NH 4 + ), phosphate (PO 4 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., 2011Khangaonkar et al., , 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 4 + or consume NO 3 -, while consumption of NH 4 + and production of NO 3 -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 4 + ; it decreases due to regenerated primary production, water column nitrification, and sediment fluxes of NO 3 -(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 2 , which can either increase or decrease DIC depending on whether the surface partial pressure of CO 2 (pCO 2 ) is lower or higher than the atmospheric pCO 2 (the latter is an input to the model).Therefore, we also introduced the calculation of the seawater pCO 2 and carbonate system constants such as the solubility of CO 2 (K 0 ) using the same equations and approach as in the CO 2 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 1 , K 2 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 1 and K 2 , 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 2 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).

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.Bianucci et al: Sensitivity of the regional ocean acidification and carbonate system in Puget Sound to ocean and freshwater inputs Art. 22, page 5 of 16 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(Khangaonkar et al., , 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 3) and then calculated the corresponding time series for DIC.For atmospheric pCO 2 concentrations, at this time we used 400 µatm for lack of 2008 data in the Salish Sea.

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

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.
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.32pH 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.

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

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, airsea 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.

Sensitivity to ocean open boundary conditions
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  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

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

Figure 1 :
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: https:// doi.org/10.1525/elementa.151.f1

Figure 2 :
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: https://doi.org/10.1525/elementa.151.f2 (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.

Figure 4 :
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: https://doi.org/10.1525/elementa.151.f4

Figure 5 :
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: https://doi.org/10.1525/elementa.151.f5

Figure 6 :
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: https://doi.org/10.1525/elementa.151.f6

Figure 7 :
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: https://doi.org/10.1525/elementa.151.f7

Figure 8 :Figure 9 :Figure 10 :
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 DIC obc (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 DIC obc 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: https://doi.org/10.1525/elementa.151.f8

Table 1 :
Description of model experiments.DOI: https://doi.org/10.1525/elementa.151.t1 obc DIC at Strait of Juan de Fuca open boundary condition increased by 40 mmol m -3 ( ~2% increase with respect to Baseline) High DIC freshwater DIC in all freshwater sources increased by 2% with respect to Baseline

Table 2 :
Statistical metrics a of model global performance for year 2008.DOI: https://doi.org/10.1525/elementa.151.t2While R 2 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. a