Future temperature and salinity in Puget Sound, Washington State, under CMIP6 climate change scenarios

Climate change will reshape estuarine ecosystems through bottom-up and top-down processes, directly affecting species at all trophic levels. To better understand future regional climate change effects on sea surface temperature and salinity, we used empirical downscaling to derive high-resolution time series of future sea surface temperature and salinity in Puget Sound (Washington State, USA). Downscaling was based on scenario outputs of two coarse-resolution Coupled Model Intercomparison Project Phase 6 (CMIP6) general circulation models, GFDL-CM4 and CNRM-CM6-1-HR. We calculated 30-year climatologies for historical and future simulations, calculated the anomalies between historical and future projections, interpolated to a ﬁ ne-resolution, and applied these downscaled anomalies to a Regional Ocean Modeling System (ROMS) time series, yielding short-term and long-term delta-downscaled forecasts. Downscaled output for Puget Sound showed temperature and salinity variability between scenarios and models, but overall, there was a strong model agreement of future surface warming and freshening in Puget Sound. Spatially, we found regional differences for both temperature and salinity, including higher temperatures in South Puget Sound and lower salinity in Whidbey Basin. Interpreting and applying downscaled CMIP6 projections of temperature and salinity will help assess climate change vulnerability and inform future ecosystem-based management decisions in Puget Sound and other coastal and estuarine systems.

Coastal acidification is being exacerbated in Puget Sound (Bianucci et al. 2018), and hypoxic water area is projected to increase from ,1 to 16% (Khangaonkar et al. 2019).These changes are expected to have negative impacts at multiple trophic levels, including for Chinook salmon (Oncorhynchus tshawytscha), coho salmon (Oncorhynchus kisutch), and Southern Resident orcas (Orcinus orca; Southern Residents hereafter)three culturally and ecologically important species in the region (Southern Resident Orca Task Force 2019).In particular, the anticipated decline in Chinook salmon could be detrimental for Southern Residents, which rely on Chinook salmon as their primary food source (Ford et al. 2010).Puget Sound faces a myriad of additional interacting anthropogenic impacts that will cause ecological, physical, and biogeochemical changes (Puget Sound Partnership 2019).
The complexity of these changes reflects the need for improved climate modeling techniques that can help describe potential local-scale impacts of climate change and develop informed climate policy (Tommasi et al. 2017).Rapidly advancing climate modeling approaches are allowing for more detailed predictions about the future of Earth's climate.The latest projections on climate change come from general circulation models (GCMs) from the Coupled Model Intercomparison Project Phase 6 (CMIP6), a major internationally coordinated effort to provide climate projections based on a variety of emissions and land use scenarios (O'Neill et al. 2016).The results from CMIP6 are used by the Intergovernmental Panel on Climate Change (IPCC) to inform global synthesis reports and create mitigation and adaptation strategies (O'Neill et al. 2016).GCMs are based on fundamental equations of physics and are dynamically coupled to include atmosphere-ocean interactions, biogeochemical cycling, land and sea ice, and soil and vegetation (Hayhoe et al. 2017;Daron et al. 2018).
Because GCMs have a coarse resolution (25-300 km), downscaling techniques must be used to obtain results at finer resolutions that are applicable at regional and local scales (Hayhoe et al. 2017;Drenkard et al. 2021).Downscaling involves translating information from GCMs to achieve finer spatial-and temporal-scale dynamics, and there are three main downscaling techniques used (Stock et al. 2011;Ekström et al. 2015).In dynamical downscaling, initial and boundary conditions for regional climate models (RCMs), which explicitly model physical processes but with a smaller, more accurate grid box size (1-50 km), are supplied by a GCM (Hayhoe et al. 2017;Tommasi et al. 2017).Dynamical downscaling is useful for modeling unprecedented and rapidly changing climate conditions, but can sometimes inherit regional biases from GCMs, and it is computationally expensive (Stock et al. 2011;Tommasi et al. 2017).Another approach is statistical downscaling, which combines historical observational data with large-scale patterns observed in the GCM output to create high-resolution results at the scale of the observational data (Hayhoe et al. 2017).One limitation of statistical downscaling is its dependence on large observational datasets, which are not available in every region, and its assumed statistical relationships (Hayhoe et al. 2017;Tommasi et al. 2017).Finally, empirical downscaling involves using a scaling factor, which is defined either as a linear dependence to a local pattern or as the difference between and historical and future climatology, to apply a projected change from a GCM to historical climate data time series or regional model output (Ekström et al. 2015).Empirical downscaling is computationally inexpensive compared to dynamical and statistical downscaling (Tommasi et al. 2017), but it cannot represent interactions between large-scale and local changes and assumes the relationships between variables will not change in the future (Ramirez-Villegas & Jarvis 2010; Ekström et al. 2015).It is likely least accurate in coastal systems that are not wellresolved by GCMs, however the low computational intensity of empirical downscaling makes it a viable method for rapid assessment, which is becoming especially important due to the urgency of understanding climate change impacts in coastal areas with large populations.These three downscaling methodologies can be effectively replicated for coastal ecosystems in which recent oceanography is well characterized and future projections of temperature and salinity can be informed by GCMs (Grieve et al. 2016;Holdsworth et al. 2021;Pozo Buil et al. 2021).Overall, downscaling techniques produce finerscale models that are relevant for making local resource management decisions (Drenkard et al. 2021).Downscaled climate change projections are particularly useful in region-specific studies regarding fisheries, protected species, and ecosystem modeling (Grieve et al. 2016;Marshall et al. 2017;Khangaonkar et al. 2019;Hollowed et al. 2020).
This study aims to quantify projected changes in Puget Sound ocean temperature and salinity under future climate change.As with many estuaries and coastal systems globally, the physics of Puget Sound and the Salish Sea over recent decades are well studied and modeled (Sutherland et al. 2011;Soontiens et al. 2016;Khangaonkar et al. 2017;MacCready et al. 2021), but high-resolution future predictions under climate change are a pressing need (Puget Sound Partnership 2019).We applied empirical downscaling techniques to derive high-resolution projections of future conditions in Puget Sound, based on low and high emission scenario outputs of coarse-resolution GCMs from CMIP6.Our focus was on deriving downscaled projections of sea surface temperature (SST) and salinity for years 2020-2050 and 2070-2100 in Puget Sound and analyzing how projections vary across Puget Sound basins.Data from this study will be used to drive an end-to-end Atlantis ecosystem model of oceanography, food webs, and human activities in Puget Sound (Morzaria-Luna et al. 2022), which will allow for a rapid assessment of vulnerability of the Puget Sound ecosystem to climate change.Empirical downscaling approaches are commonly used as a creative solution to complement dynamical downscaling studies (Drenkard et al. 2021).By incorporating the latest CMIP6 results, this study augments past oceanographic projections for Puget Sound using CMIP5, specifically Khangaonkar et al. (2019), which estimated an average ocean temperature increase of 1.51 °C in the Salish Sea under a high-emissions scenario (RCP8.5).Additionally, our study demonstrates the uses and limitations of empirically downscaling GCM output for coastal ocean ecosystems and regional modeling analysis.We offer an important example of how empirical downscaling can be efficiently developed to answer questions of ecosystem-based management (EBM) and climate-vulnerability of a marine system.In the context of Puget Sound, an urban estuary where assessment of climate change is urgently needed to protect endangered species, including salmon and the Southern Residents, our study is a first-step sensitivity analysis of Puget Sound vulnerability to SST and salinity changes.

METHODS
We developed downscaled projections for SST and sea surface salinity (practical salinity, called salinity hereafter) in Puget Sound using the Delta Method (Ramirez-Villegas & Jarvis 2010).As detailed in the following subsections, SST and salinity anomalies were calculated from between the hindcast (1984-2014) and short-term (2020-2050) and long-term (2070-2100) forecast climatologies using GFDL-CM4 (Held et al. 2019) and CNRM-CM6-1-HR (Voldoire et al. 2019) GCM model outputs.Next, we interpolated the anomalies using thin-plate spline (TPS) interpolation and applied them to a 2-year (2005)(2006) temperature and salinity time series generated by the high-resolution MoSSea Regional Ocean Modeling System (ROMS) (Modeling the Salish Sea; Sutherland et al. 2011).The adjusted time series will be used to drive the Atlantis Ecosystem Model for Puget Sounda spatially explicit, three dimensional model that integrates physical, chemical, ecological, and anthropogenic processes using dynamic, two-way coupling (Morzaria-Luna et al. 2022).All analyses were carried out using the R statistical framework (R Core Team 2021).All code is available in Github: https://github.com/steviewalker/Puget_Sound_Downscaling. Downscaled projections are archived in the Knowledge Network for Biocomplexity data repository: https://knb.ecoinformatics.org/view/urn:uuid:4e23e90b-87d8-42d4-aec5-80159d2c813d.

Study area
Puget Sound is a large, complex estuarine system located in Washington State, USA, in the Northeast Pacific Ocean.It is part of the broader Salish Sea region, which also encompasses the Strait of Georgia and the Strait of Juan de Fuca.Within Puget Sound, we distinguish four major basins for the purpose of analysis: Hood Canal, Main Basin (Admiralty Inlet and Central Puget Sound), South Puget Sound, and Whidbey Basin (Figure 1).Many different habitats, home to a wide diversity of species, are found in Puget Sound, including saltmarshes, tidal flats, and eelgrass beds (Ruckelhaus & McClure 2007;Reum et al. 2011).Puget Sound has the characteristics of a partially mixed estuary, with river input mixing with strong, three to four meter tides (MacCready 2017).The Sound has an average water residence time of 90-180 days due to strong turbulent mixing by the tides, however, because of its relatively shallow sills, circulation and dispersal of water, sediment, organisms, and contaminants can be limited in some areas (Sutherland et al. 2011;MacCready et al. 2021).Each basin has distinct physical differences due to bathymetry and differing amounts of saltwater and freshwater influence, and this complexity suggests that future conditions will involve an interplay of global trends superimposed on local, basin-specific characteristics.For example, waters in Hood Canal are highly stratified and marked by strong temperature and dissolved oxygen variations across depth, area, and season.It has a larger SST range (∼8-19 °C) than other basins (MacCready 2020).In contrast, the Main Basin experiences seasonal temperature stratification with warmer surface waters in the summer (∼14 °C), due to river input and solar radiation, and a well-mixed, cooler water column in the winter (∼8 °C) due to reduced solar radiation and wind (Ruckelhaus & McClure 2007;MacCready 2020).Across Puget Sound, there are seasonal variations in salinity, which are characteristic of coastal waters.During spring the water freshens and then salinity gradually increases as river runoff is mixed by the action of waves, winds, and tides; there are also early fall freshets rushes of fresh water flowing into the sea (Megia 1956).Average salinity in the Strait of San Juan de Fuca in the upper layer ranges from 30.7 in January to 31.6 in October, and in the bottom layer from 33.4 in February to 33.95 in August (Megia 1956).The Main Basin has an average surface salinity range of 26.5-29.0,while Hood Canal has an average surface salinity range of 21.0-26.5 (MacCready 2020).Salinities in South Puget Sound are higher than in the Main Basin; seasonal variation in salinity follows precipitation closely, with an average minimum of 28.4 in April to a maximum of 29.8 in November (Megia 1956).emission reduction and increased environmental degradation; and ssp585, business-as-usual fossil fueled world with little to no climate mitigation (O'Neill et al. 2016;Riahi et al. 2017).

ScenarioMIP and GCM selection
Two GCMs were selected for data analysis: the NOAA Geophysical Fluid Dynamics Laboratory's CM4.0 physical climate model (GFDL-CM4; Held et al. 2019) and the Centre National de Recherches Météorologiques and Cerfacs' CNRM-CM6-1-HR physical climate model (Voldoire et al. 2019).Selection was based on CMIP6 models that participated in ScenarioMIP, had the finest ocean grid resolution available (∼25 km), and spatially resolved some or all of the Salish Sea and Puget Sound (Supplementary Material, Figure S1).Fine ocean grid resolution was an important priority for the purpose of this project because it reduces uncertainty from interpolation.All of the models evaluated are summarized in Supplementary Material, Table S1, however, only two met our criteria.Thus, results in this study should be viewed more as a sensitivity test of Puget Sound to different climate change scenarios rather than projections, as a full model ensemble would be needed to make projections reliable (Brekke et al. 2008).GFDL-CM4 ran two scenarios from ScenarioMIP, ssp245 and ssp585, whereas CNRM-CM6-1 ran scenarios ssp126, ssp245, ssp370, and ssp585.GFDL-CM4 coarsely resolves all parts of the Salish Sea (25 km resolution), including Puget Sound (Supplementary Material, Figure S1).The GFDL-CM4 model simulations show a larger warming trend near the end of the historical runs as compared to observations; although the model shows low mean ocean temperature bias along the WA coast and nearby Eastern Pacific (Supplementary Material, Figure S2; Rayner et al. 2003).The CNRM-CM6-1 model shows a low ocean warming bias along the WA coast and the Pacific Ocean (Supplementary Material, Figure S2; Rayner et al. 2003), but exhibits a warm bias between summer and fall (Voldoire et al. 2019).This model has a 25 km ocean grid resolution and spatially resolves much of the Salish Sea (e.g., the Strait of Juan de Fuca and Strait of Georgia), but not Puget Sound (Supplementary Material, Figure S1).

Empirical downscalingthe Delta method
Empirical downscaling, also called the Delta Method, (Ramirez-Villegas & Jarvis 2010; Jones 2013; Ekström et al. 2015) is a relatively straightforward approach to GCM downscaling that is used to produce high-resolution seasonal climatologies for a region.It involves calculating interpolated anomalies of time series climatologies from the GCM relative to a model hindcast ( Jones 2013;Ekström et al. 2015).In order to do this, we first obtained hindcasts and forecasts from two CMIP-endorsed GCMs, CNRM-CM6-1-HR and GFDL-CM4, for our two focus variables: SST and salinity.These focus variables were chosen because they are two key physical variables that impact circulation and biological processes, they will both be impacted by climate change, and they will be used to force the Atlantis ecosystem model for Puget Sound (Morzaria-Luna et al. 2022).Available data on each CMIP6 model and scenario were obtained from the Earth System Grid Federation data repository (https://esgf-node.llnl.gov/search/cmip6/),using the ESGF Search RESTful API (https://earthsystemcog. org/projects/cog/esgf_search_restful_api).We standardized the data and calculated single delta SST and salinity climatologies for a 30-year hindcast period , short-term forecast period (2020-2050), and long-term forecast period (2070-2100) based on the different scenarios from CMIP6 (O'Neill et al. 2016).Forward simulations of GCMs models part of the CMIP6 simulations extend from 2015 to 2100, and these simulations are commonly divided in near-future and far-future time periods (see for example, Li et al. 2021 andHamed et al. 2022), while the historical period extends prior to 2014 until 1850.Anomalies were then calculated from the difference between forecasts and hindcast climatologies for both variables in each scenario from each GCM, resulting in one short-term and one long-term change factor (delta) for each model grid cell.We interpolated the anomalies using the centroids of GCM cells as points of interpolation and applied TPS interpolation across these points to obtain anomalies at finer resolutions: 3 km for visualizing the downscaled anomalies (Figure 2) and 0.1 km for mapping onto the fine-resolution ROMS model, as detailed below.

Salish Sea regional oceanographic model and the Atlantis ecosystem model
After interpolation, the 0.1-km resolution anomalies were applied to the daily-resolved, regional ocean temperature and salinity time series that drive an end-to-end ecosystem model for Puget Sound (Morzaria-Luna et al. 2022).This model is built on the Atlantis modeling framework (Audzijonyte et al. 2019) and its spatial and temporal fields of circulation, temperature, and salinity are forced with output from a high-resolution realistic numerical simulation of the Salish Sea referred to as MoSSea (Modeling the Salish Sea; Sutherland et al. 2011).MoSSea is based on the ROMS numerical framework, which solves hydrostatic, incompressible, Reynolds-averaged momentum and tracer conservation equations with a terrain-following vertical coordinate and a free surface (Shchepetkin & McWilliams 2005).The MosSSea model was forced with measured flow from 16 rivers, tides, atmospheric forcing by wind stress and heat flux, and open ocean boundary conditions.When validating this model against observed temperature and salinity CTD casts, Sutherland et al. (2011) found that the model captured the seasonal cycles of SST and salinity well, showing positive skill scores and significant positive correlations with 7 out of 8 CTD stations' observations of surface salinity and surface temperature.MoSSea showed the most discrepancies between modeled and observed data in Hood Canal, which is difficult to model because it is the narrowest and deepest part of Puget Sound; overall MoSSea was about 0.5 units saltier than observations (Sutherland et al. 2011).MoSSea was run for the years 2005 and 2006, with the initial fields in Puget Sound for each year extrapolated from available CTD cast observations.These years were chosen for the model run because of optimal overlap with observations on the coastal shelf.The choice of years is mostly arbitrary, as the 2-year period is meant to capture the seasonal cycle in Puget Sound, and the Delta Method we applied is adding a single time-step increase to the MoSSea 2005-2006 baseline output.The MoSSea model's horizontal domain is a spherical, stretched Cartesian grid extending from longitude À127°to À122°, and latitude 45°to 50°N.The grid resolution is as fine as 280 m and stretches to 3.1 km at the boundaries; there are 20 vertical layers.We chose to downscale the anomalies to a 0.1 km resolution so they would fit within the bounds of each Atlantis model polygon derived from MoSSea.In order to adapt the regular-shaped MoSSea grid to the irregular Atlantis spatial polygons, fields were averaged at 12-hour time steps, and oceanography for each Atlantis polygon was interpolated to the nearest velocity grid point on the MoSSea grid.For application in this paper, we superimposed the anomalies onto the averaged time series for each Atlantis Ecosystem Model for Puget Sound model polygon.

RESULTS
We first calculated interpolated anomalies downscaled at 3 km resolution for the entire Salish Sea and surrounding regional waters (Figure 2, Supplementary Material, Figures S3 and S4).There was an increase in SST between the short-term Results were downscaled to a resolution of 3 km.See Supplementary Material, Figure S3 and S4 for interpolated anomaly projections for other scenarios and the CNRM-CM6-1-HR model.Salinity anomalies are reported in salinity units, and SST anomalies are reported in units of °C.
(2020-2050) and long-term (2070-2100) projections, under every scenario and for both models.The CNRM-CM6-1-HR model predicted a lower SST anomaly overall than the GFDL-CM4 model.The spatial patterns of SST change were relatively consistent between the two models; slightly higher anomalies were found within the Strait of Georgia and Puget Sound for any given scenario (Supplementary Material, Figure S3).These spatial differences were more pronounced in the long-term high-emissions scenario ssp585, with significant increases in ocean heat content.For the SST interpolated anomaly projections, the GFDL-CM4 model had slightly warmer projections than the CNRM-CM6-1-HR model, with a maximum SST anomaly of 4.67 °C (long-term ssp585 projection scenario) and a minimum SST anomaly of 0.85 °C (short-term ssp245 projection scenario).The SST anomaly for CNRM-CM6-1-HR ranged between a maximum of 3.96 °C (long-term ssp585 projection scenario), and a minimum SST anomaly of 0.52 °C (short-term ssp245 projection scenario).
Salinity projections, across all models and scenarios, showed increased freshening of the Salish Sea in the long-term (Figure 2; Supplementary Material, Figure S4).However, the magnitude of freshening differs between models; the CNRM-CM6-1-HR model predicts greater freshening than the GFDL-CM4 model.Furthermore, there were differences in the spatial pattern between the models as well.The GFDL-CM4 model showed a less negative salinity anomaly, and in some cases, a positive salinity anomaly around the mouth of the Columbia River, in South Puget Sound, and in the Northern Strait of Georgia (Figure 2).In contrast, for CNRM-CM6-1-HR, the Salish Sea as a whole had a more negative salinity anomaly than the Pacific Ocean, with slightly less negative anomalies in South Puget Sound only in some projections (Supplementary Material, Figure S4).Across the spatial domain, the salinity anomaly for the GFDL-CM4 model ranged between a maximum of 0.83 salinity units (short-term ssp585 projection) and a minimum of À0.60 (long-term ssp585 projection).For CNRM-CM6-1-HR, there was a wider range in the salinity anomaly across the spatial domain: a maximum salinity anomaly of 0.17 for the short-term ssp370 projection, and a minimum salinity anomaly of À2.16 for the long-term ssp585 projection.The interpolated regional anomalies were then applied to the ROMS time series for Puget Sound (MoSSea, Sutherland et al. 2011), and these high-resolution (0.1 km) downscaled results showed variations between scenarios for each model.The seasonal patterns come from the 2005 to 2006 ROMS time series, shown as a black baseline in Figures 3 and 4. The final downscaled results reflect the high-resolution temporal patterns in the ROMS and the single time-step increase from the CMIP6 interpolated anomaly.For SST, the temporal patterns from ROMS reflect seasonal changes in solar insolation, with the warmest temperatures found in the summer and the coldest in the winter (Figure 3).For salinity, the jagged pattern is due to seasonal and annual changes in snowmelt and precipitation, with rapid decreases in salinity during the rainy winter months and during the spring snowmelt, and a slow steady increase during the drier summer months (Figure 4).
Applying the interpolated anomalies over the ROMS time series data allows us to assess temperature and salinity values at a finer spatial scale, as opposed to the anomaly results discussed previously.In both short-term and long-term projections and for every downscaled scenario, SST is projected to increase (Figure 3).SST changes are higher in the long-term projections (maximum average SST ¼ 19.6 °C, minimum average SST ¼ 8.7 °C, Supplementary Material, Table S2) compared to the shortterm projections (maximum average SST ¼ 16.5 °C, minimum average ¼ 8.0 °C, Supplementary Material, Table S2).Here, the maxima and minima refer to averages over all Atlantis Ecosystem Model for Puget Sound spatial polygons.The average downscaled SST increase from the ROMS baseline for the long-term ssp585 projection is 3.5 °C for the CNRM-CM6-1-HR model and 3.9 °C for the GFDL-CM4 model (Supplementary Material, Table S2).The intensity of SST warming increases from the lowest emissions scenario, ssp126, to the highest emissions scenario, ssp585.
Both models predict a decrease in salinity in the long-term for every scenario (Figure 4).In the short-term, the CNRM-CM6-1-HR model predicted a decrease in salinity, while the GFDL-CM4 model predicted a slight increase in salinity.Freshening was greater in the long-term for both models (maximum average salinity ¼ 30.9, minimum average salinity ¼ 25.2, Supplementary Material, Table S2) compared to the short-term (maximum average salinity ¼ 31.4,minimum average salinity ¼ 26.5, Supplementary Material, Table S2).As with SST maxima and minima, the minimum and maximum average salinity refers to the average over all spatial polygons.The average downscaled salinity decrease from the ROMS baseline for the longterm ssp585 projection is À1.9 for the CNRM-CM6-1-HR model and À0.5 for the GFDL-CM4 model (Supplementary Material, Table S2).Intensity of freshening is more pronounced at higher emission scenarios than at lower emissions scenarios.Variability in both SST and salinity change is greater in the long-term projections, which is due to greater uncertainty among SSPs, whereas the uncertainty in the GCMs is likely about the same for each time period.
Only scenarios ssp245 and ssp585 can be compared between both models, because as previously noted, GFDL-CM4 did not run scenarios ssp126 and ssp370.Overall, the GFDL-CM4 model results were slightly warmer and saltier than the CNRM-CM6-1-HR model results (Figure 5).Salinity results showed more inter-model variation than SST.In other words, the two models had very high agreement for SST results, but salinity, while highly correlated between models, was greater in the GFDL-CM4 model than in CNRM-CM6-1-HR.There is greater range in salinity in the long-term compared with the short-term (Figure 5), which is likely an outcome of model configuration differences between the GFDL-CM4 model and the CNRM-CM6-1-HR model.
Final downscaled results overlaid onto the ROMS data showed basin-specific spatial differences within Puget Sound.Overall, the highest average projected SSTs were found within Hood Canal and South Puget Sound, while the lowest projected SSTs were found within the Eastern Strait of Juan de Fuca (Figure 6).The projected freshest waters were found near areas of major river input, mainly the Snohomish and Skagit Rivers in the Whidbey basin (Figure 6).Projected salinity was higher near the Strait of Juan de Fuca and the Strait of Georgia (Figure 6).The degree of freshening and warming increased from the short-term to the long-term projections for each scenario, and the most severe projected freshening and warming occurred at higher emission scenarios in the long-term (Supplementary Material, Figures S6 and S7).

Range of climate sensitivity in empirically downscaled projections
We developed short-term (2020-2050) and long-term (2070-2100) downscaled projections for SST and salinity in Puget Sound, using two GCM models, CNRM-CM6-1-HR model and GFDL-CM4.Over the 21st century, both models predict warming and freshening of Puget Sound.This general result is in-line with the Puget Sound-specific findings in Khangaonkar et al. (2019) and the broader findings in the California Current (Xiu et al. 2018).Warming and freshening is greater in the long-term than in the short-term.The severity of these effects is dependent on the CMIP6 emissions scenario.In ssp585, representing a fossil-fuel dependent, high-emissions world, there will be greater warming and freshening than in ssp126, which represents a sustainable, low-emissions world (Riahi et al. 2017).The CNRM-CM6-1-HR model and the GFDL-CM4 model have strong agreement on SST increases and less agreement on salinity.Spatially, the projected SST increase in Puget Sound is the greatest in Hood Canal and South Puget Sound, and the projected salinity decrease is the greatest in Main Basin.
The models differed in their predictions of future anomalies, though more in terms of salinity than SST.The GFDL-CM4 model had slightly warmer predictions compared to the CNRM-CM6-1-HR model, up to 0.4 °C greater (Supplementary Material, Table S2).Salinity patterns were more varied than SST, with the GFDL-CM4 model showing a slight increase in salinity in the short-term and the CNRM-CM6-1-HR showing greater long-term freshening, up to 1.4 units fresher (Supplementary Material, Table S2).This distinction between models does not mean one is necessarily more accurate than the other; it is likely a result of different ocean model component configurations and differences in the terrestrial model component configurations that could affect projected precipitation and snowmelt changes.The CNRM-CM6-1-HR model uses the NOCS-ORCA1 configuration of the ocean model NEMO version 3.6 (Voldoire et al. 2019), and the GFDL-CM4 model uses the MOM6 configuration of OM4.0 (Held et al. 2019).The different parameterizations used in the GCMs could account for different model responses to the same forcing, as in the case of salinity; internal variability and potential biases in GCMs can be exacerbated at regional scales given their coarse resolution (Brunner et al. 2020).Furthermore, both regional models and GCMs tend to reproduce surface and subsurface salinity worse than temperature, with a typical bias toward fresher than observed conditions (Renault et al. 2021;Liu et al. 2022).CNRM-CM6-1-HR has a larger global average surface salinity bias compared to GFDL-CM4 (-0.445 vs. À0.121s;Liu et al. 2022), which could be one additional reason CNRM-CM6-1-HR predicts greater freshening in Puget Sound (Figure 4).Unlike salinity, both of the models exhibited a low temperature bias along the Washington coast and Pacific Ocean (Supplementary Material, Figure S2), which could explain greater agreement in SST results.Including additional models to expand a multi-model comparison would have helped understand how choice of GCM influences model predictions due to uncertainty in different model parameterizations (Tommasi et al. 2017).However, choosing the GCM that resolved Puget Sound with the best spatial resolution was more important for the purpose of this project in order to reduce uncertainty from interpolating.The GDFL-CM4 model and the CNRM-CM6-1-HR model were the only two models we found that spatially resolved part or all of the Salish Sea on their horizontal grid, participated in ScenarioMIP, and had the data needed available for download (Supplementary Material, Table S1 and Figure S1).
Each scenario resulted in different SST and salinity anomalies for Puget Sound, and the anomalies increased from the lowemissions scenario ssp126 up to the high-emissions scenario ssp585.Within the ScenarioMIP framework, scenarios have been updated relative to CMIP5 to include both the representative concentration pathways (RCPs) of radiative forcing (W/m 2 ), and new shared socioeconomic pathways (SSPs) that identify different societal outcomes over the 21st century (O'Neill et al. 2014;Riahi et al. 2017).For the four different scenarios used in this study, ssp126 represents the most sustainable development, lowest emissions, reduced socioeconomic and political inequality, and the most reduced climate change impacts (Riahi et al. 2017).This scenario also coincides with the least amount of freshening and warming in Puget Sound.Warming and freshening increases in the intermediate scenarios -ssp245 and ssp370and is greatest in scenario ssp585, the business-as-usual scenario in which fossil-fuel reliant development occurs and energy intensive lifestyles become globally widespread (Riahi et al. 2017).The range of outcomes for Puget Sound are dependent upon which scenario the world follows over the 21st century, which we see in Figures 3 and 4 as the downscaled time-steps for each scenario increase from ssp126 to ssp585.
The spatial differences in Figure 6 are primarily generated from the MoSSea ROMS model (Sutherland et al. 2011), which simulates circulation in the Salish Sea.The warmest regions of downscaled projected SST, Hood Canal, and South Puget Sound (Figures 1 and 6) are both connected to the Main Basin by narrow passages that are constrained in part by shallow sills.SST in these regions is likely greater due to a more constricted water flow and a narrower topography.We see lower SSTs in the Northern Basin, where there is strong tidal inflow from the Strait of Juan de Fuca into Puget Sound (∼16,000 m 3 s À1 ) (Khangaonkar et al. 2019).The freshest waters in Figure 6 are seen near areas of major river input, such as the Whidbey Basin, where 50% of freshwater input to Puget Sound comes from the Skagit river (Sobocinski 2021) and both the Skagit and Snohomish rivers have an average discharge of ∼1,000 m 3 s À1 to this basin (Sutherland et al. 2011).However, because in Puget Sound the tidal influence is stronger than river input, the effects of river discharge on salinity are only noticeable in regions nearest to river mouths.The average downscaled results echo the interpolated anomalies from the CMIP6 results because the Delta Method is adding a step increase to the ROMS time series.Imposing the anomalies to the ROMS time series preserves seasonal (Figures 3 and 4) and regional (Figure 6) differences from the ROMS, and provides high-resolution values for SST and salinity that are more useful than anomalies alone for EBM applications, such as the Atlantis Ecosystem Model for Puget Sound (Morzaria-Luna et al. 2022).

Potential uses and limitations of empirical downscaling for estuaries and coastal systems
The empirical downscaling methodology used in this paper is an important first step to developing CMIP6 climate projections for Puget Sound until a more detailed dynamical downscaling, similar to the methodology used for the Salish Sea by Khangaonkar et al. (2019) with CMIP5, can be carried out.However, since only two CMIP6 GCMs spatially resolved part or all of the Salish Sea, we did not conduct downscaling with a full model ensemble.Thus, these results should be viewed as a sensitivity test of Puget Sound to different climate change scenarios rather than projections (Brekke et al. 2008).For the purpose of ecosystem assessment, conservation planning, and other regional applications, empirical downscaling provides a method of rapid and computationally inexpensive climate assessment with high climate realism (Ramirez-Villegas & Jarvis 2010; Ekström et al. 2015;Drenkard et al. 2021), which is useful given the urgency of developing climate change adaptation strategies for highly populated coastal regions, like Puget Sound.It is important to remember that empirical downscaling does not make GCM output more accurate or reliable ( Jones 2013;Drenkard et al. 2021).Empirical downscaling cannot account for changes in variability on the local scale because the climate signal is coming from the coarse-resolution GCM (Ekström et al. 2015).Regional feedbacks, such as local precipitation changes or circulation changes, are also not wellsimulated by the GCM (Mearns et al. 2003;Ekström et al. 2015).Furthermore, empirical downscaling also assumes that the anomaly patterns calculated from the GCMs will hold true in the future (Ramirez-Villegas & Jarvis 2010), which we cannot assume will be the case for Puget Sound as complex regional processes cannot be modeled with the coarse resolution of the GCMs.Drenkard et al. (2021) suggest that empirical methods, similar to our delta approach here, are not a panacea but can help us frame key uncertainties and trends in data-rich systems such as Puget Sound.Moreover, the review by Drenkard et al. (2021) identifies other instances of empirical and statistical downscaling, in the context of managing living marine resources.Though empirical methods may be less skillful (see Ekström et al. 2015), we view them as a rapid stopgap tool to begin exploring CMIP6 climate change implications for Puget Sound marine species, and those in other regions.In particular, these methods are a tool to begin to transition from qualitative, expert-based climate risk assessment for stocks such as endangered salmon (Crozier et al. 2019) to quantitative projections informed by ecosystem models (Morzaria-Luna et al. 2022).Such quantitative models are needed to conserve and manage coastal species such as salmon, which face cumulative impacts from climate change, urbanization, contaminants, fisheries, and other actions (Sobocinski 2021).Thus, even with the limitations of empirical downscaling, it is the best method available for efficient incorporation into ecosystem modeling.
We chose to use GCMs that resolved Puget Sound with the highest initial resolution available (25 km; Supplementary Material, Figure S1) in order to make the downscaled anomalies as sensitive to the region as possible.Even still, GFDL-CM4 simulates Puget Sound with a small number of grid cells, and CNRM-CM6-1-HR only simulates the Strait of Georgia and the Strait of Juan de Fuca (Supplementary Material, Figure S1).A 25 km resolution is too coarse to represent many key circulation patterns that influence temperature and salinity in Puget Soundincluding explicit tidal forcing and horizontal shear across the Strait of Juan de Fuca, which often has opposing flows on each side of the Strait (Thomson et al. 2007).GCM resolution is also too coarse to model alterations in streamflow and the amount of precipitation, which are two factors that will drive salinity changes in Puget Sound (Sobocinski 2021).These processes will not be captured in the downscaled anomalies.However, Puget Sound circulation is well-simulated at a high resolution in the MoSSea regional model (Sutherland et al. 2011), to which we applied our downscaled anomalies.Thus, for the purpose of developing a first estimate of the sensitivity of Puget Sound to a range of future climate conditions under CMIP6, the choice of downscaling method may be less critical at this timescale.
Though we provide both short-term and long-term projections, we expect uncertainty to be higher in the long-term projections.If we categorize uncertainty in terms of scenario uncertainty, model uncertainty, and internal variability (Cheung et al. 2016), we expect scenario uncertainty to increase and dominate over the long-termin other words, the global decisions about emissions become stronger drivers of predicted ocean conditions than decisions about GCMs and downscaling technique.Overall, empirical downscaling was chosen for the purpose of rapid, straightforward assessment, but dynamical downscaling for Puget Sound using CMIP6 models will be an important next step in developing more accurate SST and salinity predictions.

Management and ecosystem modeling applications of downscaled projections
Improving accuracy in SST and salinity projections will be important in the coming decades because it will allow better evaluation of the ecosystem impacts of climate change in Puget Sound.In particular, understanding how climate change will affect all trophic levels and culturally significant species such as Chinook salmon, Pacific herring (Clupea pallasii), and Southern Residents is of key relevance to Puget Sound.Both Southern Residents and Puget Sound Chinook salmon are listed under the Endangered Species Act (as endangered and threatened, respectively), and Pacific herring are a critical forage species for salmon and protected seabirds and marine mammals.Increased ocean temperatures can be harmful to salmon during multiple stages of their life cycle, impacting spawning and migration, increasing mortality, and the risk of pathogens (Battin et al. 2007;Mauger et al. 2015).Higher ocean temperatures could also lead to higher Pacific herring embryo mortality (Villalobos et al. 2020).Furthermore, the loss in snowpack will reduce salmon spawning area in rivers in the Puget Sound watershed, leading to an expected decline in salmon population (Battin et al. 2007).This has been identified as the highest Puget Sound salmon exposure risk (Crozier et al. 2019).Freshening and warming has the potential to lead to more stratification in Puget Sound, which directly affects primary production through a changing nutrient supply (Mauger et al. 2015;Xiu et al. 2018).Warming will also continue to change phytoplankton dynamics in Puget Sound, leading to increased harmful algal blooms, which can be toxic to fish (Southern Resident Orca Task Force 2019; Sobocinski 2021).The fate of Southern Resident orcas is fundamentally tied to Chinook salmon, the orca's primary food source (Ford et al. 2010).Compounding ecosystem stressors may have rippling effects at every trophic level, potentially leading to largely reduced food resources for the Southern Residents (Southern Resident Orca Task Force 2019).Climate change will increase ecosystem vulnerability to other anthropogenic impacts, such as population growth projected in coastal areas.
We produced downscaled anomalies for the whole Salish Sea and downscaled time series for Puget Sound.Downscaled anomalies can be used to assess changes in species distributions (Grieve et al. 2016;Petatán-Ramírez et al. 2019), and to guide resource-specific, basin-specific, or ecosystem-scale climate adaptation and resilience strategies, which is an emerging priority for the Puget Sound Partnership (2019), the State of Washington agency tasked with recovery of Puget Sound habitats, resources and services.The downscaled time series also will be used to drive scenarios using the Atlantis model for Puget Sound, a deterministic simulation model designed to support strategic decision making for marine resource management (Audzijonyte et al. 2019).Within Atlantis, temperature directly affects primary production, respiration, and other metabolic processes, and both temperature and salinity can dictate habitat use (Audzijonyte et al. 2019).Therefore, downscaled ocean projections will directly influence simulated growth of individuals, population dynamics, and spatial and trophic relationships of species ranging from phytoplankton to fish and marine mammals.Results from Atlantis ecosystem model simulations will help inform recommendations for EBM in Puget Sound (Puget Sound Partnership 2019), particularly when it comes to evaluating how climate change, management decisions, ecological changes, and other influences will effect Chinook salmon populations, as well as endangered Southern Resident orcas (Morzaria-Luna et al. 2019).Other Atlantis ecosystem modeling projects have used downscaled projections from GCMs to drive their models, such as the Benguela and Agulhas Currents Atlantis model (Ortega-Cisneros et al. 2018), the California Current Atlantis model (Marshall et al. 2017), and the Nordic and Barents Sea Atlantis model (Hansen et al. 2019).Empirical downscaling is an effective way to obtain finer spatial scale resolutions needed for making ecosystem level analysis.
As discussed in Section 4.2, our downscaling approach has limitations, and these must be considered before and during applications related to management and decision making.For instance, our downscaled projections will drive Atlantis ecosystem models, which are intended as a strategic tool to evaluate cumulative impacts and categories of policy options (Fulton et al. 2011), not to dictate year-by-year tactical decisions.Thus, the goal is to allow the models to include scenarios (not necessarily precise predictions) of realistic climate change and other drivers such as nutrients and contaminants, with the goal of identifying management strategies robust to these changes (Fulton 2011).In this context, empirical downscaling of Puget Sound oceanography is a useful first step.For complex ecosystem models, such as Atlantis, to move from science into policy, they are often complemented by an ensemble of simpler models (Kaplan et al. 2019), and are evaluated in terms of fits to data and model performance (Kaplan & Marshall 2016).

CONCLUSIONS
Overall, warming SST and freshening are expected to intensify the most in the long-term high-emissions scenarios (ssp585), where climate change impacts from a fossil-fuel intensive economy are the greatest.The downscaled projections show agreement between the GFDL-CM4 model and CNRM-CM6-1-HR model, but there is more variability for salinity.These SST and salinity changes are driven by warming atmospheric temperatures and shifts in precipitation patterns in Washington State.Warming ocean temperatures are already causing dynamic shifts at all levels of the food web, particularly an increase in harmful algal blooms (Sobocinski 2021), reduction in survival for threatened Chinook salmon (Crozier et al. 2019), and increasing starvation risk for the endangered Southern Residents (Southern Resident Orca Task Force 2019).Because density is dominated by salinity relative to temperature, circulation in Puget Sound does not appear to be affected by warming ocean temperatures (Khangaonkar et al. 2019).Understanding how ocean conditions will vary based on CMIP6 climate change scenarios is critical for evaluating the intensity of change in Puget Sound, within individual basins and in the ecosystem as a whole.With high-resolution downscaled SST and salinity projections, we will be able to more accurately model ecosystem conditions using the Atlantis Ecosystem Model for Puget Sound (Morzaria-Luna et al. 2022) and complete rapid assessments of EBM decisions.Incremental improvements in climate downscaling will support informed policy decisions and conservation recommendations that will be critical for informing the vital signs of Puget Sound (Puget Sound Partnership 2019) and understanding climate change threats to this ecosystem.
CMIP6 builds upon climate modeling advances made in CMIP5 by introducing a new set of scenarios from the Scenario Model Intercomparison Project (ScenarioMIP), which integrates future emission projections with societal concerns (O'Neill et al. 2016).Four different SSPs were used in this study: ssp126, rapid emission reductions and intensive climate mitigation; ssp245, some emission reductions and intermediate climate mitigation; ssp370, international disparities in

Figure 1 |
Figure 1 | Map of the Salish Sea.The four main basins in Puget Sound include Hood Canal, South Puget Sound, Main Basin (Admiralty Inlet and Central Puget Sound), and Whidbey Basin.

Figure 2 |
Figure 2 | Interpolated anomalies for the GFDL-CM4 model scenario ssp585.The figure shows both short-term and long-term projections.Results were downscaled to a resolution of 3 km.See Supplementary Material, FigureS3and S4 for interpolated anomaly projections for other scenarios and the CNRM-CM6-1-HR model.Salinity anomalies are reported in salinity units, and SST anomalies are reported in units of °C.

Figure 3 |
Figure 3 | Average SST change for each scenario as forced by the ROMS time series.Results are shown for both models and the long-term and short-term projections.The black line shows the temporal patterns that come from the ROMS baseline, which ran over a 2-year period from 2005 to 2006.Scenarios are color coded and range from low (ssp126) to high (ssp585) emissions.ssp126 and ssp370 are not included in the GFDL-CM4 model because it only ran two scenarios, ssp245 and ssp585.

Figure 4 |
Figure 4 | Average sea surface salinity change for each scenario as forced by the ROMS time series.Results are shown for both models and the long-term and short-term projections.The black line shows the temporal patterns that come from the ROMS baseline, which ran over a 2-year period from 2005 to 2006.Scenarios are color coded and range from low (ssp126) to high (ssp585) emissions.ssp126 and ssp370 are not included in the GFDL-CM4 model because it only ran two scenarios, ssp245 and ssp585.

Figure 5 |
Figure 5 | Variability between models for average SST and salinity change for ssp585 short-term and long-term projections.Temporal patterns come from the ROMS time series, which ran from 2005 to 2006.See Supplementary Material, Figure S5 for the ssp245 model variability.

Figure 6 |
Figure 6 | Spatial variations in downscaled projected SST and salinity for GFDL-CM4 ssp585 high-emissions scenario.Results are shown for both the short-term and long-term.Choropleth map polygons come from the Atlantis Ecosystem Model for Puget Sound polygons (Morzaria-Luna et al. 2022).The value within each polygon is the calculated average over the course of the 2-year ROMS time series with the downscaled anomaly overlaid.See Supplementary Material, Figures S6 and S7 for additional choropleth maps.