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

Continuing emissions of fossil fuels and increasing land use are changing Earth's climate in rapid and unprecedented ways. In recent decades, anthropogenic climate change has had increasing effects on the ocean; ocean warming, ocean acidification, and ocean deoxygenation are three primary stressors that are impacting ocean biogeochemistry globally (Gruber 2011; Bopp et al. 2013). Increasing ocean temperature can change the rate of biological processes (Bopp et al. 2013), decrease species fitness (Mora et al. 2013), and cause species distribution shifts (Grieve et al. 2016; Petatán-Ramírez et al. 2019). Circulation changes are expected as well, since warmer surface oceans will result in a more stratified water column in the open ocean (Gruber 2011), though nearshore and estuarine dynamics may be more complex. Additionally, higher temperatures result in decreased ocean solubility and the creation of more hypoxic areas (oxygen-minimum zones) (Keeling et al. 2010; Gruber 2011; Mora et al. 2013). Marine absorption of CO₂ is lowering overall ocean pH and carbonate availability (acidification), making it more difficult for calcifying organisms to build their shells (Doney et al. 2009; Kroeker et al. 2013). The combination of impacts on ocean biogeochemistry and additional climate change-related stressors has synergistic effects on marine ecosystems that vary by region, including altered food-web dynamics, community composition, and energy flows, as well as a potential reduction in the ecosystem services provided by marine systems (Gruber 2011; Doney et al. 2012). In order to develop effective mitigation and conservation strategies, climate change impacts must be investigated further on a regional and local basis (Marshall et al. 2017).

Puget Sound, a large fjordal estuarine system within the Salish Sea on the North American Pacific coast, faces many climate change-related impacts. Here, increasing temperatures are projected to result in earlier and reduced spring snowmelt, shifts in timing and intensity of upwelling, species distribution changes, decreased salmon migration and survival, and increased bloom seasons for harmful algal bloom species (Reum et al. 2011; Khangaonkar et al. 2019; Sobocinski 2021). 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 well-resolved 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 finer-scale 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.

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 Sound – a 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.

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

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

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.

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

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 short-term 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 long-term 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.

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 low-emissions 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²), 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 ssp370 – and 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³s⁻¹) (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³s⁻¹ 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).

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 well-simulated 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 Sound – including 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-term – in 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.

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

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.

S.W., H.M-L., and I.K. contributed to conception and design; S.W. and H.M-L. contributed to acquisition of data; S.W., H.M-L. and D.P.-R. contributed to analysis and interpretation of data; S.W., H.M-L., I.K. and D.P.-R. drafted and/or revised the article; S.W., H.M-L., I.K. and D.P.-R. approved the submitted version for publication.

We would like to thank Chris Harvey for reviewing a previous version of the document. Manuel Dorantes assisted with correcting the references.

This work was funded by the WA SeaGrant project R/HCE-25 ‘Evaluating the effects of Southern Resident orcas recovery actions and external threats in the marine ecosystem of Puget Sound’ (Hem Nalini Morzaria-Luna and Isaac Kaplan) and by the NOAA Ernest F. Hollings Undergraduate Scholarship Program (Stevie Walker). Computing resources were funded by the Microsoft Azure AI for Earth grant (Hem Nalini Morzaria-Luna).

All relevant data are available from an online repository or repositories. All code is available in Github: https://github.com/steviewalker/Puget_Sound_Downscal-ing. Downscaled projections are archived in the Knowledge Network for Biocomplexity data repository: https://knb.ecoinformatics.org/view/urn:uuid:4e23e90b-87d8-42d4-aec5-80159d2c813d.

The authors declare there is no conflict.