A Dynamically Downscaled Ensemble of Future Projections for the California Current System

Given the ecological and economic importance of eastern boundary upwelling systems like the California Current System (CCS), their evolution under climate change is of considerable interest for resource management. However, the spatial resolution of global earth system models (ESMs) is typically too coarse to properly resolve coastal winds and upwelling dynamics that are key to structuring these ecosystems. Here we use a high-resolution (0.1°) regional ocean circulation model coupled with a biogeochemical model to dynamically downscale ESMs and produce climate projections for the CCS under the high emission scenario, Representative Concentration Pathway 8.5. To capture model uncertainty in the projections, we downscale three ESMs: GFDL-ESM2M, HadGEM2-ES, and IPSL-CM5A-MR, which span the CMIP5 range for future changes in both the mean and variance of physical and biogeochemical CCS properties. The forcing of the regional ocean model is constructed with a “time-varying delta” method, which removes the mean bias of the ESM forcing and resolves the full transient ocean response from 1980 to 2100. We found that all models agree in the direction of the future change in offshore waters: an intensification of upwelling favorable winds in the northern CCS, an overall surface warming, and an enrichment of nitrate and corresponding decrease in dissolved oxygen below the surface mixed layer. However, differences in projections of these properties arise in the coastal region, producing different responses of the future biogeochemical variables. Two of the models display an increase of surface chlorophyll in the northern CCS, consistent with a combination of higher nitrate content in source waters and an intensification of upwelling favorable winds. All three models display a decrease of chlorophyll in the southern CCS, which appears to be driven by decreased upwelling favorable winds and enhanced stratification, and, for the HadGEM2-ES forced run, decreased nitrate content in upwelling source waters in nearshore regions. While trends in the downscaled models reflect those in the ESMs that force them, the ESM and downscaled solutions differ more for biogeochemical than for physical variables.

Given the ecological and economic importance of eastern boundary upwelling systems like the California Current System (CCS), their evolution under climate change is of considerable interest for resource management. However, the spatial resolution of global earth system models (ESMs) is typically too coarse to properly resolve coastal winds and upwelling dynamics that are key to structuring these ecosystems. Here we use a highresolution (0.1 • ) regional ocean circulation model coupled with a biogeochemical model to dynamically downscale ESMs and produce climate projections for the CCS under the high emission scenario, Representative Concentration Pathway 8.5. To capture model uncertainty in the projections, we downscale three ESMs: GFDL-ESM2M, HadGEM2-ES, and IPSL-CM5A-MR, which span the CMIP5 range for future changes in both the mean and variance of physical and biogeochemical CCS properties. The forcing of the regional ocean model is constructed with a "time-varying delta" method, which removes the mean bias of the ESM forcing and resolves the full transient ocean response from 1980 to 2100. We found that all models agree in the direction of the future change in offshore waters: an intensification of upwelling favorable winds in the northern CCS, an overall surface warming, and an enrichment of nitrate and corresponding decrease in dissolved oxygen below the surface mixed layer. However, differences in projections of these properties arise in the coastal region, producing different responses of the future biogeochemical variables. Two of the models display an increase of surface chlorophyll in the northern CCS, consistent with a combination of higher nitrate content in source waters and an intensification of upwelling favorable winds. All three models display a decrease of chlorophyll in the southern CCS, which appears to be driven by decreased upwelling favorable winds and enhanced stratification, and, for the HadGEM2-ES forced run, decreased nitrate content in upwelling source waters in nearshore regions. While trends in the downscaled models reflect those in the ESMs that force them, the ESM and downscaled solutions differ more for biogeochemical than for physical variables.

INTRODUCTION
The California Current System (CCS) is one of the four global eastern boundary upwelling systems (EBUSs) characterized by extraordinary biological productivity that supports a variety of human uses including tourism, fisheries, and recreation (e.g., Checkley and Barth, 2009). As in other EBUS, the high level of primary productivity in the CCS is primarily attributed to wind-driven coastal upwelling, which delivers nutrient-rich deep waters to the surface. However, upwelled waters also have low pH, contain high concentrations of respired carbon, and are moderately oxygen-poor, making this region prone to hypoxia and acidification, conditions that can threaten both benthic and pelagic marine life (e.g., Gruber et al., 2012;Turi et al., 2014;Feely et al., 2016). Simulations of changes in the CCS that affect those processes are therefore key to evaluate the direct ecological and economic impact of the future climate on this dynamic ecosystem (Vecchi et al., 2006;Bakun et al., 2010;Sydeman et al., 2014;Rykaczewski et al., 2015;Howard et al., 2020b).
Under future climate scenarios, projected physical changes in the CCS include enhanced stratification, shifts in the timing and intensity of upwelling, and alteration of the properties and relative contributions of source waters advected into the region (e.g., Rykaczewski and Dunne, 2010;Doney et al., 2012;Rykaczewski et al., 2015;Bograd et al., 2019). Ocean warming generally results in increased water column stability (Capotondi et al., 2012), reducing vertical mixing, nutrient supply, and primary productivity in the euphotic zone of subtropical regions (Hoegh-Guldberg and Bruno, 2010). In EBUS regions, however, it has been proposed that heterogeneous warming of the land and ocean could intensify the surface atmospheric pressure gradient and consequently enhance equatorward winds and coastal upwelling (i.e., the Bakun hypothesis; Bakun, 1990;Bakun et al., 2015). Global climate models suggest that the Bakun hypothesis is overly simplified; projected changes in upwelling are likely to be season-and latitude-dependent, with intensification of CCS upwelling during spring, especially in the southern CCS, followed by a significant weakening of upwelling during the summer months, especially in the northern region of the CCS (e.g., García-Reyes et al., 2015;Rykaczewski et al., 2015;Wang et al., 2015). These anthropogenic trends in upwelling are not likely to be emergent until mid-century (∼2050) in the southern CCS and even later (∼2080) in the northern CCS (Brady et al., 2017), as strong decadal variability remains the dominant signal on shorter time horizons.
In the biogeochemical realm, reduced ventilation and circulation of the North Pacific may increase nutrient concentrations, lower pH, and decrease oxygen concentrations in the deep source waters of the CCS (Rykaczewski and Dunne, 2010;Van Oostende et al., 2018;Xiu et al., 2018). These source water changes could result in increased primary productivity over the continental shelf, which, when the increased organic matter is subsequently remineralized and combined with reduced oxygen in source waters, would result in a significant expansion of hypoxic areas (Dussin et al., 2019). Similarly, Hauri et al. (2013) projected that, by 2050, the nearshore mean surface pH of the CCS will move outside the envelope of present-day variability and the aragonite saturation horizon of the central CCS will shoal into the upper 75 m, causing near-permanent undersaturation in subsurface waters. However, many of these biogeochemical projections describe results from individual climate models; biogeochemical responses to climate change in the CCS vary considerably across models  and the robustness of projected changes must be determined.
Earth system models (ESMs) are essential tools for future climate studies. ESMs are atmosphere-ocean-land-sea ice general circulation models (GCMs) that have been coupled to biogeochemical models (Taylor et al., 2012). While the initial motivation for global ESM development was resolving the partitioning of anthropogenic carbon emissions between the land, ocean, and atmosphere (e.g., Friedlingstein et al., 2006;Arora et al., 2013Arora et al., , 2020, the inclusion of ocean biogeochemical components provided new insights into ocean acidification, deoxygenation, and productivity changes (Bopp et al., 2013;Kwiatkowski et al., 2020). However, the utility of ESMs on regional scales and especially in upwelling systems is limited by their coarse resolution (Boville and Gent, 1998;Mote and Mantua, 2002;Palmer, 2014). The ocean component in most GCMs is too coarse to resolve fine-scale processes including coastal upwelling, mesoscale eddy activity, and coastal trapped waves (Stock et al., 2011). Their atmosphere component is also too coarse to resolve nearshore structure in the winds, which can exert a strong influence on the physical and biogeochemical signature of coastal upwelling (e.g., Mote and Mantua, 2002;Capet et al., 2004;Jacox and Edwards, 2012;Small et al., 2015;Renault et al., 2016;Franco et al., 2018). Dynamical downscaling (i.e., using coarse global fields to force a high-resolution regional model) is one commonly used approach to produce highresolution ocean projections at regional scale (Drenkard et al., in review). This approach has been used in other EBUS (Benguela Upwelling, e.g., Machu et al., 2015;Humboldt Upwelling, e.g., Echevin et al., 2012Iberian Upwelling, Miranda et al., 2013), and in the CCS (Auad et al., 2006;Li et al., 2014;Xiu et al., 2018;Arellano and Rivas, 2019;Dussin et al., 2019;Howard et al., 2020a). For the CCS, these previous high-resolution projections were conducted by downscaling single climate projections forced by just one ESM, analyzed short future periods (10-30-year time-slices), or applied idealized perturbations in some variables that did not account for all physical and biogeochemical forcing. They often do not consider inherent present-climate bias in ESMs or capture the full envelope of future climate uncertainties, which in the case of biogeochemistry projections significant contribution comes from model uncertainties Frölicher et al., 2016).
Dynamically downscaled projections allow us to recover important regional-scale features that are missing or poorly represented in their coarse-resolution parents models. But, the high-resolution projection is tied to the large-scale forcing from the ESMs, including their biases, and performing a bias correction on the ESM forcing has the potential to significantly improve dynamical downscaling simulations (Bruyère et al., 2014;Machu et al., 2015;Xu et al., 2019;Drenkard et al., in review). Here we use a high-resolution regional oceanbiogeochemical coupled model and apply a "time-varying" delta approach to dynamically downscale three different ESMs. The selected ESMs span a wide range of potential future changes in both the mean and variance of physical and biogeochemical properties for the CCS. Our objectives are to produce and present downscaled projections of climate change in the CCS at 0.1 degree. (i.e., mesoscale resolving) horizontal resolution from 1980 until 2100 under the Representative Concentration Pathway (RCP) 8.5 scenario. We examine a range of ecosystem-relevant variables including sea surface temperature (SST), chlorophyll, and subsurface nitrate and oxygen. We analyze changes in those variables for the middle and end-of century and the plausible mechanisms associated with those changes. The remainder of the paper is organized as follows: section "Materials and Methods" describes the regional modeling approach, the ESM output used and the bias correcting method applied to construct the forcing of the regional model. Section "Results" describes the projected changes in ecosystem-relevant variables, and discussion of the results and conclusion are presented in sections "Discussion" and "Concluding Remarks, " respectively.

Coupled Physical-Biogeochemical Model: ROMS-NEMUCSC
To produce the high-resolution future projections in the CCS, we use the Regional Ocean Modeling System (ROMS) coupled with a biogeochemical model (NEMUCSC) based on the North Pacific Ecosystem Model for Understanding Regional Oceanography (NEMURO). ROMS is a free-surface, hydrostatic, primitive equation ocean model that uses stretched, terrainfollowing coordinates in the vertical and orthogonal curvilinear coordinates in the horizontal (Shchepetkin and McWilliams, 2005;Haidvogel et al., 2008). The model configuration, developed by the University of California Santa Cruz, covers the region 30-48 • N and 115.5-134 • W (midway down Baja California to just south of Vancouver Island) with 0.1 • (∼10 km) horizontal resolution and 42 vertical sigma levels (Veneziani et al., 2009 1 ). NEMURO is a medium complexity nutrient-phytoplankton-zooplankton (NPZ) model specifically developed and parameterized for the North Pacific (Kishi et al., 2007). The model includes three limiting macro-nutrients (i.e., nitrate, ammonium, and silicate), two phytoplankton groups (nanophytoplankton and diatoms), three zooplankton groups (micro-, meso-, and predatory zooplankton) and three detritus pools (DON, PON, and opal). Here, we use a customized version of NEMURO, called NEMUCSC, which is specifically parameterized for the CCS and augmented with oxygen and carbon cycling (Cheresh and Fiechter, 2020;Fiechter et al., 2018Fiechter et al., , 2020. The coupled ocean-biogeochemistry regional model will be hereafter referred to as "ROMS-NEMUCSC." We first use ROMS-NEMUCSC to perform a historical control simulation (CTRL) for 1980-2010. Initial and ocean lateral boundary conditions are derived from the Simple Ocean Data Assimilation version 2.1.6 (SODA; Carton and Giese, 2008) at monthly resolution. The atmospheric surface forcing is derived from the global atmospheric reanalysis from European Centre for Medium-Range Weather Forecasts version 5 (ERA-5; Hersbach et al., 2020) at 1-h and ∼30 km resolution, with the exception of surface winds, which were obtained from ERA-5 for 1980-1987 and from the 0.25 • (∼25 km) Cross-Calibrated Multi-Platform wind product (CCMP1; Atlas et al., 2011) for 1988-2010 at 6h resolution. In the nearshore region of the CCS, especially off northern California and Oregon, we have found that CCMP1 more closely reproduces the observed magnitude of summertime winds and leads to better representation of biogeochemical processes near the coast. Air-sea fluxes are computed in ROMS internally using bulk formulae (Fairall et al., 1996a,b;Liu et al., 1979). Physical variables are stored daily and used to force the biogeochemical component "offline" (i.e., NEMUCSC is run independently and driven by daily surface atmospheric fluxes and oceanic fields from ROMS). Initial and boundary conditions are derived from the 2009 World Ocean Atlas (WOA) climatology (Garcia et al., 2010a,b) for nutrients and oxygen, from the Global Ocean Data Analysis Project (Key et al., 2004) for total alkalinity, and from the empirical relationship of Alin et al. (2012) for DIC using monthly temperature and oxygen. Initial and boundary conditions for ammonium, phytoplankton, zooplankton, and detritus are set to a small value (0.1 mmol N m −3 ), noting that the evolution of these quantities is dominated by surface ocean dynamics that adjust rapidly to simulated macronutrients in the interior of the model domain.  Table 1. We downscaled each ESM for the period from 1980 to 2100 using historical forcing  and the RCP8.5 climate change scenario . We focus on the RCP8.5 scenario since model uncertainty dominates scenario uncertainty for biogeochemical (and to a lesser degree, physical) changes in the CCS. Indeed, for biogeochemical variables in the CCS the range of projections under the RCP2.6 and RCP4.5 scenarios is fully contained within the model uncertainty under the RCP8. The question of which ESMs to include in a multi-model ensemble is an active research area. The models that best match the observed historical climate are not necessarily the ones that will most faithfully represent future climate sensitivity, and more process-based model selection methods (e.g., with emerging constraints) are being developed (Eyring et al., 2019;Hall et al., 2019). Given the lack of definitive selection criteria, we chose three ESMs that capture the CMIP5 range of projected future changes in both the mean and variance of physical and biogeochemical CCS properties. While they show agreement in the sign of the projected SST change in the CCS, they differ in the magnitude of warming and disagree on the sign of the primary production (PP) change (Figure 1). HAD projects the warmest SST anomalies by 2100 and a sharp decline in PP by around 2050. Relative to the CMIP5 ensemble mean, GFDL and IPSL project weak and moderate increases in SST, respectively, and both project increased PP by 2100 (Figure 1), counter to the declining PP trend in the CMIP5 ensemble mean (Bopp et al., 2013). Similarly, these three models have different projections for changes in the magnitude of interannual variability under climate change. The three models span the range of potential changes in physical and biogeochemical variance, which may increase, decrease, or remain unchanged depending on the variable and model (Supplementary Figure 2).

Downscaling Approach
Output from the selected ESMs provide the surface and ocean boundary conditions for ROMS-NEMUCSC in a oneway downscaling approach. To reduce biases exhibited in the ESMs historical simulations, we apply a "time-varying delta" bias-correction method prior to downscaling the ESMs output. We first estimate time-varying deltas (DELTA) for a forcing atmospheric variable (ATM) by subtracting the ESM's historical monthly climatology (ESM CLM ; years 1980-2010) from the whole period of interest (i.e., 1980-2100).
DELTA 1980−2100 = ESM 1980−2100 − ESM CLM,1980CLM, −2010 (1) These monthly time-varying deltas are first bilinearly interpolated in space and time to the resolution of the reanalysis data used to force the control run (REAN), and then added to the monthly reanalysis climatology (REAN CLM ). Finally, we add high-frequency variability from the reanalysis (REAN HF ), computed as the residual after removing a 30-day running mean, as the absence of this high-frequency variability and associated damping of upwelling can lead to a biased ecosystem state (Gruber et al., 2006(Gruber et al., , 2011 Bias correction for the ocean and biogeochemical variables is handled similarly, but with several adjustments for the lower frequency of available ESM output. The high-frequency component (REAN HF ) is excluded, since the available ESM output for 3D physical and biogeochemical variables have monthly and annual temporal resolution, respectively. In addition, the time-varying deltas for the biogeochemical variables are calculated annually, again because the ESM biogeochemical outputs are available at annual resolution ( Table 1). Model initialization is handled similarly to the ocean boundaries by applying deltas from the year 1980 to our CTRL initial conditions. The time-varying delta method for downscaling retains the observed historical climatology and high-frequency (submonthly) variability while inheriting the long-term change and interannual variability of the parent ESM. Relative to a "fixed delta" method that compares a historical period to a future one (e.g., Alexander et al., 2019;Shin and Alexander, 2020), the time-varying delta method has advantages of capturing projected changes in interannual variability, and resolving the full climate change evolution, including potentially non-linear impacts that would be missed when the transient response is excluded. An alternative method for prescribing the high-frequency variability would be to use the daily ESM forcing, which may be modulated by low-frequency variability (e.g., El Niño Southern Oscillation events alter storm tracks in the CCS). We retain the observed historical high-frequency variability from the reanalysis for several reasons. First, it allows us to force the model with highresolution temporal forcing (6-hourly for the winds and hourly for the rest of ATMs), which resolves the daily cycle and is consistent with the historical forcing. Second, high frequency variability is likely to differ substantially within the space of a single ESM grid cell, particularly near the coast, and will not be resolved by the ESMs. In other EBUS, previous works applied statistical downscaling methods in the wind forcing to overcome the limitation of the coarse resolution models (e.g., Goubanova et al., 2011;Machu et al., 2015;Bonino et al., 2019). While future changes in high frequency atmospheric forcing is not the focus of this study, it should be a topic of further research.
We refer to the downscaled simulations from ROMS-NEMUCSC, driven by GFDL, IPSL, and HAD as ROMS-GFDL, ROMS-IPSL, and ROMS-HAD, respectively, in the following sections and figures.

Observation Data for Model Evaluation
To evaluate the CTRL and the historical period of the downscaled simulations (i.e., ROMS-GFDL, ROMS-IPSL, and ROMS-HAD), we use the Optimum Interpolation SST data set, which combines in situ and satellite-based observations, from the National Oceanic and Atmospheric Administration (NOAA OISST v2; Reynolds et al., 2007) Garcia et al., 2013a;nitrate, Garcia et al., 2013b;temperature,Locarnini et al., 2013).

RESULTS
In the following sections, we describe projected changes in the future physical and biogeochemical states of the CCS, first for the long-term mean changes for the 2070-2100 period with respect to the historical period (1980-2010) (i.e., Future -Historical) and then the time-dependent changes for the full transient runs .

Model Evaluation
The historical CTRL reproduces many aspects of the CCS physical and biogeochemical variability (Supplementary  Figures 3, 4). The simulated annual mean SST has a weak warm bias over most of the domain and a weak cold bias along the northern boundary. The Southern California Bight (SCB) shows the largest SST bias (∼1.4 • C), which is comparable with other numerical simulations of the CCS (e.g., Veneziani et al., 2009;Renault et al., 2016) and is likely due to the poor representation of cloud cover and the cross-shore gradient in alongshore winds resulting from coarse spatial resolution in the forcing reanalysis data (Renault et al., 2020). The annual mean subsurface patterns of nitrate and dissolved oxygen in the CTRL simulation are also well represented compared to climatological values from the WOA. Although the observed and simulated values of nitrate and oxygen differ near the coast, which cannot be resolved by the coarse horizontal resolution of the WOA grid (1 • × 1 • ), the model reasonably reproduces the offshore spatial gradient and variability of climatological subsurface nitrate and dissolved oxygen (Supplementary Figures 3, 4). The observed and simulated annual mean chlorophyll exhibits similar onshore-offshore gradients. However, the main differences are located in the nearshore regions, where the model overestimates chlorophyll between the northern San Francisco Bay (SFB) and Cape Blanco, and underestimates it in the northern coastal domain (>44 • N), SFB, and the SCB (Supplementary Figure 3).

Alongshore Wind Stress and Vertical Velocity
We use meridional wind stress as a proxy for the alongshore or upwelling favorable winds. In general, the three ESMdriven ROMS simulations exhibit an overall intensification (i.e., more negative) and a northern displacement (>39 • N) of the meridional wind stress (Figure 2), consistent with prior analysis Xiu et al., 2018). These equatorward wind stress anomalies oppose the mean wind stress (poleward) at the northern latitudes. However, there are some notable intermodel differences in the spatial variability. ROMS-GFDL and ROMS-IPSL show a strong increase centered in the northern nearshore regions (<200 km from the coast) of the CCS, while ROMS-HAD shows the center of this intensification offshore and in the northwest corner of the domain. While ROMS-GFDL exhibits the strongest intensification across the CCS domain, ROMS-IPSL exhibits a weaker intensification that is confined along the coastal domain, and ROMS-HAD shows this intensification confined to the northern part of the CCS. A strong decrease in alongshore wind stress over the southern part of the domain is projected in both ROMS-IPSL and ROMS-HAD, offshore in the former and nearshore in the latter.
We use vertical velocity at 50 m as a measure for upwelling (Figure 2). The projected increase in the equatorward meridional wind stress over the northern coastal region of the CCS enhances coastal upwelling (upward and positive vertical velocity changes) in the three ESM-driven downscaled ROMS simulations. While ROMS-GFDL and ROMS-IPSL projections show an overall upward velocity all along the coast of the domain, the intense weakening of the alongshore-wind stress projected by ROMS-HAD reduces coastal upwelling in the southern part of the domain.

Sea Surface Temperature
The SST response to projected climate change includes warming over the entire domain in all three ESM-driven ROMS simulations (Figure 2). Although the three ROMS simulations agree on the warming of the CCS region, the magnitude of change differs considerably among simulations. ROMS-GFDL exhibits the weakest projected warming (<2.5 • C), ROMS-IPSL a moderate warming (>3 • C), and the ROMS-HAD the strongest warming (>3.5 • C). In addition to the differences in overall climate sensitivity, there are differences among the spatial warming patterns: in the ROMS-IPSL the warming offshore of the SCB is relatively smaller, while ROMS-HAD presents a more homogeneous warming pattern. In ROMS-IPSL, this weaker warming extends north along the coast up to ∼43 • N. In ROMS-GFDL, the warming is lessened in the northern coastal waters (>39 • N), suggesting that increased wind-driven upwelling may partially mitigate the warming locally, while the strongest warming occurs in the offshore waters of the SCB.

Subsurface Nitrate
The ecological effects of upwelling critically depend on the biogeochemical properties and composition of upwelled source waters. To explore the subsurface biogeochemical properties, we select a depth of 150 m (Chhak and Di Lorenzo, 2007), which is deeper than the seasonal historical mixed layer depth (MLD) which ranges from 30 to 100 m in the whole domain. All three downscaled simulations project an increase of the subsurface nitrate concentration in the offshore waters (>100 km from the coast), with a larger increase by the end of the 21st century in ROMS-GFDL (∼5 mmol/m 3 ) than in the other two simulations (∼2.5-3 mmol/m 3 ) (Figure 3). However, the projections differ in the coastal waters (<100 km from the coast). While ROMS-GFDL shows a coastwide nitrate increase, the nitrate increase in ROMS-IPSL is much weaker and limited to the southern and northern CCS, and ROMS-HAD projects a decrease in subsurface nitrate along the coast, which is strongest in the southern part of the domain (Figure 3).

Subsurface Oxygen
Below the seasonal pycnocline, upwelling systems are characterized by both high nutrient concentrations and low oxygen concentrations. The projected dissolved oxygen changes at 150 m depth mirror the nitrate changes described above -all three models exhibit large-scale oxygen declines. ROMS-HAD diverges from the others with a weak increase in dissolved oxygen that extends in the coastal waters from the southern boundary up to ∼40 • N (Figure 3). Other ecologically important measures of the oxygen environment -bottom oxygen concentration and depth of the hypoxic boundary -show similar changes. While ROMS-GFDL and ROMS-IPSL projections feature a decline of bottom dissolved oxygen on the shelf, which is more intense FIGURE 2 | Control run (1980-2010) mean and future changes (2070-2100 relative to historical) for the physical variables. Rows show from top to bottom: meridional wind stress, vertical velocity at 50 m depth (w), and sea surface temperature (SST). Columns show data for the historical period in the control simulation (CTRL, first column) and the future changes from the high-resolution downscaled projections: ROMS-GFDL (second column), ROMS-IPSL (third column), and ROMS-HAD (fourth column). Negative values of wind stress changes correspond to an increase in equatorward wind stress. Black contour marks 100 km from shore. Vertical velocity field has been spatially smoothed out for better representation.
(∼50 mmol m 3 ) in the northern and southern coastal waters, ROMS-HAD projects a weak increase (∼25 mmol m 3 ) of the bottom dissolved oxygen that extends from the southern coastal waters up to 40 • N (Figure 3). Future changes in depth of the hypoxic boundary layer have a similar spatial pattern as the O 2 at 150 m. The hypoxic boundary, defined here as the depth at which dissolved oxygen = 1.43 ml/l (i.e., 63.86 mmol/m 3 ), is projected to shoal by up to ∼150 m in the offshore waters. Along the coast there is a weak shoaling, up to ∼60-70 m, in ROMS-GFDL and ROMS-IPSL, while ROMS-HAD projects changes of similar magnitude but opposite sign off southern and central California (Figure 3).

Phytoplankton Biomass
We present total phytoplankton biomass -the sum of NEMUCSC's small (nanophytoplankton) and large (diatoms) phytoplankton groups -in units of chlorophyll concentration, which is estimated from NEMUCSC's nitrogen units using fixed ratios for C:N (106:16) and C:Chl (50:1 and 100:1 for small and large phytoplankton, respectively, Goebel et al., 2010). By the end of the century, ROMS-HAD projects a decrease of surface chlorophyll along the entire coast; ROMS-IPSL projects a weaker decrease along most of the coast but a slight increase in the northern regions; and in contrast, ROMS-GFDL projects an increase centered in in the northern CCS (>39 • -45 • N), coincident with the region of most pronounced increases in upwelling favorable winds (Figure 4). Offshore, ROMS-IPSL and ROMS-HAD project a decrease in surface chlorophyll everywhere, while in ROMS-GFDL surface chlorophyll change is mostly positive north of 35 • N. For ROMS-GFDL and ROMS-HAD, projected changes in vertically integrated chlorophyll over the upper 50 m exhibit similar patterns to surface chlorophyll changes. However, for ROMS-IPSL the slight increase in surface chlorophyll in the nearshore northern CCS is offset by subsurface chlorophyll decrease such that the upper 50 m integrated chlorophyll is projected to decline (Figure 4). In sum, phytoplankton biomass is projected to decline throughout the study domain with the exception of the northern CCS, where enhanced biomass is projected in ROMS-GFDL and ROMS-IPSL, though only for a shallow nearshore layer in the latter.

Vertical Sections
All three downscaled simulations show an overall increase in temperature in the subsurface, with the warming being surfaceintensified such that there is also an increase in stratification (i.e., more tightly packed potential density contours or isopycnals) by the end of the century (Figure 5). ROMS-HAD and ROMS-IPSL project warming changes up to ∼2 • C at 300 m depth, and warming of ∼1 • C extends to 500 m or more. As for SST, ROMS-GFDL indicates less subsurface warming than the other two projections, with weak warming in the northern latitudes, almost no change in the southern latitudes, and even a slight cooling of ∼0.5 • C around ∼150 m depth in the nearshore southern CCS (32 • N; Figure 5).
The inter-model consistency of the projected nitrate changes is spatially dependent. All three projections show subsurface increases in nitrate concentration offshore, though ROMS-GFDL shows a stronger increase with greater horizontal and vertical  extent compared to the other models (Figure 6). However, the changes diverge near the coast. At 32 • N and ∼200 m depth, in the core of the California Undercurrent (CU; i.e., ∼26.5 kg m −3 ), there is an increase of nitrate concentration in ROMS-GFDL but a decrease in ROMS-HAD. In both models, the inshore change in concentration extends north of 32 • N, but with reduced magnitude (Figure 6). In contrast, ROMS-IPSL exhibits no discernible change in nearshore nitrate concentrations.
The future changes of spatial patterns of dissolved oxygen generally follow those of nitrate concentration with opposite signs (Figure 7), though the oxygen declines are more pronounced and widespread than the nitrate increases. There is  a strong decrease in dissolved oxygen nearly everywhere in all of the projections, with the most pronounced decreases occurring slightly deeper in the south than in the north. Accompanying the nitrate concentration changes around the core of the CU, ROMS-GFDL and ROMS-IPSL project a decrease in oxygen, which is again stronger in ROMS-GFDL, while ROMS-HAD projects an increase of oxygen in the CU core at 32 • N that vanishes quickly to the north. The association of enhanced nitrate and decreased oxygen is consistent with "older" subsurface waters adjacent to the California Current upwelling (e.g., Rykaczewski and Dunne, 2010; waters adjacent to the California Current upwelling have been isolated from the ocean surface for a longer period, allowing more nitrate to accumulate and more oxygen to be lost through aerobic remineralization of sinking organic material).

Alongshore Wind Stress
Projections of wind stress across the three simulations share a common meridional pattern with intensification of equatorward (upwelling-favorable) wind stress in the northern CCS and unchanged or slightly weakened winds in the southern and central CCS (Figure 8). This pattern is consistent with previous analysis of CMIP5 models , but there are noticeable differences between simulations in the magnitude and timing of change throughout the 21st century. In ROMS-GFDL, the strongest intensification in upwelling favorable winds, centered between 40 and 45 • N, becomes apparent by 2030 and persists through the end of the century. ROMS-IPSL shows low frequency variability without strong linear trends; in particular there is a period of considerably weakened upwelling favorable winds north of 35 • N from approximately 2060-2080 followed by an intensification from 2080 through the end of the century. The prominence of these signals throughout the projections suggests the influence of decadal variability in this system, which can delay the emergence of anthropogenic upwelling trends until the late 21st century (Brady et al., 2017). ROMS-HAD has a consistent trend of equatorward wind stress intensification in northern latitudes and weakening in southern latitudes, with the division between increasing and decreasing upwelling favorable winds at ∼40 • N.

Sea Surface Temperature
Surface ocean warming is one of the most robust ocean responses to climate change, and all three simulations project increases in SST across all latitudes (Figure 8). SST continues to increase with time in all simulations, with ROMS-GFDL projecting the weakest increase (∼2.5 • C by 2090), ROMS-IPSL projecting a stronger increase (∼4 • C by 2090), and ROMS-HAD projecting the strongest increase (∼5 • C by 2090). ROMS-GFDL projects slightly stronger warming at southern latitudes, while ROMS-IPSL and ROMS-HAD project nearly uniform warming across the meridional extent of the coastal domain.

Subsurface Nitrate Concentration
Regional Ocean Modeling System-Geophysical Fluid Dynamics Laboratory projects a strong increase in subsurface nitrate concentration, with a trend that emerges from the natural variability (i.e., exceeds 1 standard deviation of the historical variability) by ∼2050 and is strongest from 30 • to 35 • N (Figure 8). ROMS-IPSL exhibits a similar spatial pattern, but with a reduced trend that exceeds decadal variability only in the southern part of the domain and that is inherited from the parent ESM. The nearshore decrease in nitrate concentration in ROMS-HAD emerges clearly in roughly 2060 in the southern latitudes and spreads northward across most of the domain by the end of the century (Figure 8).

Subsurface Dissolved Oxygen
Projected dissolved oxygen concentrations generally follow those of the nitrate concentrations, with opposite signs, across all three models (Figure 8). The tight coupling between nitrate and oxygen is evident both in the decadal variability and in the long-term trends, though the oxygen trends are even more pronounced. While all three models project a strong decrease in dissolved oxygen, particularly in the northern latitudes, the time when the dissolved oxygen changes becomes prominent differs among models. Trends emerge from the natural variability across all latitudes by ∼2060 and ∼2040 in ROMS-GFDL and ROMS-IPSL, respectively. In ROMS-HAD, an oxygen decrease has already emerged in the northern latitudes (>43 • N), while decadal variability remains the dominant signal in the rest of the domain, and beyond 2070, a weak increase in dissolved oxygen emerges in the southern part of the domain, consistent with the decrease of nitrate there.

Vertically Integrated Phytoplankton Biomass
All three models project the main changes in chlorophyll, vertically integrated over 0-50 m, along a latitude band between 35 and 45 • N, the most productive portion of the CCS (Figure 8). ROMS-HAD and ROMS-IPSL project a chlorophyll decrease starting in ∼2040 and ∼2050, respectively, near 37 • N, which strengthens and extends across most of the coastal domain by the end of the century in ROMS-HAD. In contrast, ROMS-GFDL projects an increase of vertically integrated chlorophyll (∼8 mg m 2 that becomes distinct ∼2050), off Northern California and Oregon (∼40-45 • N) and a decrease off Central California (∼35-37 • N). These changes in ROMS-GFDL are likely driven by both the intensification of the alongshore winds, which increase the upwelling between 40 • and 45 • N, and the strong increase of the subsurface nitrate concentration. In ROMS-HAD, the negative trends in vertically integrated chlorophyll are largely consistent with trends in the weakening of the upwellingfavorable wind stress and the decrease of subsurface nitrate concentration along the southern (<39 • N) CCS coast. However, along the northern (>39 • N) CCS coast, these negative trends in chlorophyll indicate that the intensification of alongshore winds cannot offset the depletion of subsurface nitrate. However, the dynamics driving chlorophyll changes are not necessarily straightforward, as phytoplankton biomass can be influenced by changes in upwelling strength, subsurface nutrients, upper ocean stratification, and top-down control by grazers. The relative influence of these potentially competing drivers under historic and future forcing continues to be a topic of interest, which we discuss more in the following section.

Coastal Trends ESMs Comparisons
To highlight the impact of the downscaling process, we compare coastal trends of key ecosystem variables with those of the ESMs (Figure 9). All projections show a strong warming of the coastal CCS. SST increases rapidly in the ROMS-IPSL and in ROMS-HAD from the late 2020s, and late 2030s. By the end of the century, SST increases in ROMS-IPSL reach ∼4 • C and >3.5 • C in the northern and southern CCS, respectively, while in ROMS-HAD warming reaches >4.5 • C in both the northern and southern regions. ROMS-GFDL shows a moderate increase in SST, starting later (∼2040), and reaching ∼2 • C and >2.5 • C, in the northern and southern CCS, respectively, by the end of the century. SST trends in the downscaled projections strongly followed those from the ESMs, resulting in no significant spread between the downscaled projections and their ESMs counterparts. Compared to the evolution of the ROMS-ESM SST means, the SST ESMs exhibit an offset of the trends due to their different bias with respect to the historical period (∼1-2 • C, Supplementary Figure 5).
Biogeochemical variables in ROMS-ESMs show different signs in the trends along the CCS coast. These changes and their meridional differences in magnitude are consistent with those from ESMs (Figures 9, 10). In ROMS-GFDL and ROMS-HAD, the evolution of subsurface nitrate and oxygen agree more closely with their ESM counterparts. However, the evolution of these anomalies in ROMS-IPSL separate from those in the IPSL, showing a significant intra-model spread especially in the southern coastal region.
Regional Ocean Modeling System-Institut Pierre Simon Laplace and ROMS-HAD depict a moderate and a strong decline of coastal chlorophyll along the coast, reaching −0.5 and −10 mg m −2 by the end of the century, respectively. Meanwhile, ROMS-GFDL depicts a moderate increase of chlorophyll in the northern CCS coast. The evolution of the ROMS-GFDL chlorophyll trends agrees more closely with those from GFDL, showing a pronounced decadal variability. However, ROMS-IPSL and ROMS-HAD show a larger spread in chlorophyll with respect to their counterparts. While anomalies of chlorophyll in IPSL remain unchanged during the century, the chlorophyll in ROMS-ISPL declines, showing larger variability than the chlorophyll projection in IPSL (Figure 9).
Because of the inherent ESM bias with respect to the historical period, the evolution of the mean values of the biogeochemical variables also exhibit an offset compared to those in the ROMS-ESMs. Offset magnitudes are larger in subsurface nitrate and oxygen in GFDL and IPSL, while chlorophyll offset is larger in HAD (Supplementary Figure 5).

DISCUSSION
We used a coupled physical-biogeochemical model (ROMS-NEMUCSC) at 0.1-degree (∼10 km) resolution to produce downscaled projections of climate change in CCS from 1980 to 2100 under the high emissions RCP8.5 scenario. To capture the spread of projections, we selected three ESMs, GFDL-ESM2M, HadGEM2-ES, and IPSL-CM5A-MR, that span the CMIP5 range for future changes in both the mean and variance of physical and biogeochemical CCS properties. The downscaled runs provide information on potential future physical and biogeochemical states in the CCS at higher resolution than those available from current generation ESMs. More importantly, they provide insight into which future changes are robust (i.e., with different projections agreeing on the sign of change, if not the magnitude) and which are uncertain even in a qualitative sense (i.e., when models diverge on the sign of changes). They also illustrate  benefits and limitations of dynamical downscaling in the context of climate change projection. We discuss these points below.

Similarities and Differences Among ESM-Driven ROMS Projections
All three downscaled models show an overall intensification of the meridional wind stress (i.e., equatorward wind stress or upwelling-favorable winds) in the northern part of the domain (>40 • N), consistent with previous work Wang et al., 2015;Xiu et al., 2018), with ROMS-IPSL showing the weakest intensification and ROMS-GFDL the strongest. In the southern region, ROMS-GFDL and ROMS-IPSL show a weaker intensification of the upwelling-favorable wind stress, while ROMS-HAD projects a weakening wind stress. Although, ROMS-GFDL and ROMS-IPSL simulations project an intensification of the upwelling-favorable wind stress for the last thirty years of the end of the century, ROMS-IPSL projects a weakening of the meridional wind stress in coastal waters from approximately 2060-2080. Since low frequency variability in the dynamically downscaled surface winds remains dependent on the overlying larger-scale wind, this period of winds stress relaxation along the CCS coast in ROMS-IPSL may reflect variability in the North Pacific represented by the ESM model parent IPSL. Further analysis demonstrates that decadalscale variability in upwelling-favorable wind stress in all three models is correlated with basin-scale climate oscillations (i.e., the Pacific Decadal Oscillation, PDO; Mantua et al., 1997). For example, weakened upwelling-favorable wind stress in ROMS-IPSL during 2050-2080 occurs during a positive phase of the PDO (Supplementary Figures 6, 7, and Supplementary Table 1). Without careful assessment, an analysis of the projections from ROMS-IPSL for this particular and short period of time (i.e., 2060-2080) could lead to misattributing this signal to anthropogenic climate change instead of natural climate variability (e.g., Deser et al., 2012). Thus, it is important to consider the simulation's full transient period, or longer time-scales projections, as recommended in the downscaling protocols by Drenkard et al. (in review). It also justifies our approach of using a time-varying delta method to generate the downscaled projections.
In the subsurface offshore, the three ESM-driven ROMS simulations project similar enrichment of nitrate and a decrease in dissolved oxygen as well as a shoaling of the hypoxic boundary layer, all of which follow their ESM counterparts. In the coastal regions, ROMS-GFDL and ROMS-IPSL biogeochemical responses are similar in sign and aligned to those off the coast, featuring an average of ∼30-40% increase/decrease of subsurface nitrate/dissolved oxygen, ∼25% reduction of bottom oxygen on the shelf, and a ∼50 m shallower hypoxic boundary layer. ROMS-HAD projects a different response featuring a decrease of subsurface nitrate (∼20%) that by 2070 extends along the entire coast, an increase of subsurface oxygen (30%) in the southern region that extends to the bottom of the shelf, and a deepening of the hypoxic boundary layer by ∼30 m. Changes in subsurface nitrate and oxygen in our ROMS projections are largely consistent with recent downscaling studies that include these ESMs (Xiu et al., 2018;Howard et al., 2020a). One exception is subsurface oxygen in ROMS-HAD along the coast, for which we find an increase in oxygen coupled to the decrease in subsurface nitrate. Howard et al. (2020a) find that subsurface oxygen decreases in their downscaled HAD projection, though that decrease is much weaker than in their GFDL and IPSL projections. A plausible explanation for this may be associated with the downscaling method and/or the specific ensemble member from HAD used to force the highresolution projections.
Changes in nutrient concentrations and nutrient ratios have already been observed in upwelling source waters of the CCS (e.g., Bograd et al., 2015). ROMS-GFDL and ROMS-IPSL projected trends in nitrate enrichment and deoxygenation are consistent with previous studies that suggest that biogeochemical changes in upwelling source waters could be a first order effect in changes in the biogeochemistry of the CCS, with the decrease in ventilation of North Pacific interior waters, increase in stratification and the subsequent deepening of the isopycnals, being the proposed driving mechanisms (Rykaczewski and Dunne, 2010;Xiu et al., 2018;Dussin et al., 2019;Howard et al., 2020a). Moreover, deoxygenation of the upwelling source waters is the main driver of changes in hypoxia and expansion of hypoxic areas on the shelf (Dussin et al., 2019). Under the future RCP8.5 scenario, HAD and IPSL project relatively small changes of Equatorial Undercurrent (EUC) transport compared to GFDL, accompanied by a slight compression of the oxygen minimum zones (OMZ) in the former and expansion in the latter (Shigemitsu et al., 2017). These models' biases (i.e., representation of the EUC) and discrepancies in the future projections of the complex equatorial dynamics and the OMZ (Cabré et al., 2015;Busecke et al., 2019) could explain the diverse trends in the downscaled projections. The strong decrease of subsurface oxygen near the southern boundary projected by ROMS-GFDL might be associated with the weakening of the EUC and transport of older waters (i.e., less oxygen and higher nutrient concentrations) into the CCS from the south. In contrast, the increase of subsurface oxygen in the projected ROMS-HAD might be related to the compression of the OMZ and the transport of younger waters to the CCS. Future changes in the equatorial dynamics and their connections to the CCS source waters will be explored in future work.
The projected enrichment of nutrients in subsurface waters, combined with increased upwelling caused by the intensification of the meridional wind stress, could be the driving mechanisms for increases in the phytoplankton biomass of ∼50% in ROMS-GFDL. Enrichment of nutrients near the depth of upwelled source waters may increase primary productivity in the CCS (Rykaczewski and Dunne, 2010), but our results suggest that the increase in productivity is reinforced by the intensification of upwelling favorable winds projected to occur in the northern coastal CCS in ROMS-GFDL. Conversely, nutrient-poor source waters combined with the weakening of upwelling winds (15-18%) likely drive the projected strong decrease of chlorophyll (∼50%) in the coastal CCS by ROMS-HAD. This differing response of ROMS-HAD could come from inheriting the large-scale biogeochemical changes in HAD that propagate through the CCS boundaries.
All the models show that future ocean warming in the CCS is surface intensified, indicating enhanced stratification (Figure 5 and Supplementary Figure 8). In the coastal domain, ROMS-GFDL and ROMS-IPSL project the larger stratification and reduced (i.e., shallower) MLD around ∼34-40 • N (between Point Conception and Cape Mendocino), while ROMS-HAD projects a weaker stratification and less shoaling of the MLD in that region (Supplementary Figure 8). Increased stratification can potentially counteract the effect of intensifying winds and source water nutrient enrichment by limiting the nutrient flux to the mixed layer and decreasing PP (Di Lorenzo et al., 2005;García-Reyes et al., 2015;Jacox et al., 2015). Thus, enhanced stratification, in combination with weakened upwelling-favorable winds, could explain why phytoplankton biomass decreases are projected in the southern CCS by ROMS-GFDL and along most of the coast by ROMS-IPSL despite increasingly nutrient-rich subsurface waters. The uncertainty related to the interplay between the thermal stratification and intensification of upwelling winds should be explored further.

Impact of Regional Downscaling
Given the added effort, computational cost, and storage required associated with dynamical downscaling, it is worth considering the added benefit relative to the uncertainties in the parent model forcing. We compared the ROMS projections with those from their ESM counterparts and showed that both SST and chlorophyll trends are qualitatively consistent between ROMS simulations and their ESM parents (Figure 9). ROMS-ESM SST changes follow those of their ESM parents closely in magnitude, producing a similar range of the uncertainties. From the ESM parent models, GFDL and IPSL agree on the sign of the NPP and surface chlorophyll change (i.e., increasing through the century), while the CMIP5 ensemble mean projects declines in NPP in this region with the same sign trend as HadGEM2-ES, but with much smaller magnitude (Figure 1). Along the CCS coast, ROMS-GFDL and ROMS-HAD also agree on the trend sign and the latitudinal expressions of the chlorophyll changes compared to their ESM parents, but not in the magnitude in the latter. On the other hand, ROMS-IPSL projects a clear negative trend in chlorophyll while the trend in IPSL is weaker (Figure 9 and Supplementary Figure 9). Differences in the sign and magnitude of the trends might be related to nonlinear interactions in the biogeochemistry as a response to the upwelling dynamics that are better resolved in the downscaled models. Compared to the SST changes, some ROMS-ESM biogeochemical changes largely differ in magnitude from their ESM parents, producing larger intra-model spread in particular between ROMS-HAD and HAD in chlorophyll, and ROMS-IPSL and IPSL in subsurface nitrate and oxygen. However, there is an overall reduction of the uncertainty of the biogeochemistry projections in the ROMS-ESMs, likely due to using the same biogeochemical model and parameterization for all three downscaled simulations. A summary of these main physical and biogeochemical changes in the CCS is shown in Figure 10.
The larger biogeochemical differences between ROMS simulations and their ESM parent models illustrate the impact of regional downscaling's ability to resolve local processes impacting productivity (Echevin et al., 2020). Various factors may explain why biogeochemical variables are more sensitive to regional downscaling than physical variables: (1) highspatial resolution for the circulation is needed to simulate upwelling dynamics and their ecosystem responses; (2) complex nonlinear interactions occur between the physical changes and ecosystem responses; and (3) the formulation of the biogeochemical component varies among ESMs and is different than NEMUCSC, which has been specifically calibrated for the CCS.

CONCLUDING REMARKS
The high-resolution downscaled projections presented here belong to a broader end-to-end modeling project (Future Seas 2 ) that integrates climate, biogeochemical, ecosystem, and socioeconomic models to evaluate the impacts of climate change in the CCS and evaluate fisheries management strategies under projected future conditions. The regional ocean models provide the physical and biogeochemical foundation for higher trophic level models of various types, including individual based models and species distribution models (Smith et al., 2021), highlighting the utility of high-resolution projections for regional marine resource applications. However, given that differences between downscaled models and their ESM parents can be small relative to the differences between ESMs, any specific application should consider the value of dynamical downscaling in the context of other sources of uncertainty. Until large ensembles of eddyresolving global or regional models are computationally feasible, we suggest that a fruitful approach is to combine coarser resolution large ensembles with dynamical downscaling of select runs informed by analyses like the one presented here to assess how representative basin-scale changes are translated to shelfscale responses.

DATA AVAILABILITY STATEMENT
The regional model output is available on request from the corresponding author. CMIP5 model output is available through a distributed data archive developed and operated by the Earth System Grid Federation data portal (ESGF, https://esgf-node. llnl.gov/search/cmip5/). ERA-5, the fifth generation of ECMWF atmospheric reanalyses of the global climate, is available from the Copernicus Climate Change Service Climate Data Store website (https://cds.climate.copernicus.eu/cdsapp#!/home). The Cross-Calibrated Multi-Platform (CCMP) wind vector analysis product version 1.1 is available from (http://www.remss.com/ measurements/ccmp/). NOAA High Resolution SST data by the NOAA/OAR/ESRL PSL, Boulder, Colorado, USA is available from their website https://psl.noaa.gov/. Sea-viewing Wide Fieldof-view Sensor (SeaWiFS) chlorophyll data from the NASA Goddard Space Flight Center, Ocean Ecology Laboratory, Ocean Biology Processing Group is available from their website https: //oceancolor.gsfc.nasa.gov/data/seawifs/. Finally, World Ocean Atlas database series are derived from the Word Ocean Database, and are available from https://www.nodc.noaa.gov/ OC5/indprod.html.

AUTHOR CONTRIBUTIONS
MPB and MGJ designed the research. MPB and MAA downscaled the ESM models. MPB conducted the physical variables projections, and analyzed the data. JF produced the biogeochemical variables projections. MPB, MGJ, JF, MAA, SJB, ENC, CAE, RRR, and CAS wrote the manuscript.

FUNDING
This work was supported by a grant from the National Atmospheric and Oceanic Administration (NOAA) CPO Coastal and Climate Applications (COCA) program and the NOAA Fisheries Office of Science and Technology (NA17OAR4310268).

ACKNOWLEDGMENTS
Many thanks to Isaac Schroeder and Liz Drenkard for helpful comments on an earlier version of the manuscript.