Drivers of Future Extratropical Sea Surface Temperature Variability Changes in the North Paci�c

.


Introduction
Anthropogenic emissions of greenhouse gasses are causing profound changes to the Earth's climate.Changes to the climate mean state have been studied for over half a century (e.g., ref. [1]) and are often used to set targets for reducing greenhouse gas emissions.In contrast, changes to climate variability-characterized statistically by variance and occurrence of extreme events and of importance for regional adaptation strategies-under future warming scenarios are less well understood.
The recent advent of large ensemble climate model simulations offers an opportunity to robustly quantify future variance and extreme event changes [13,[16][17][18].
Conducting a large number of simulations with the same climate model with identical external forcing but perturbed initial conditions allows for a clear identification of the forced signal as it changes over time, leaving only model and scenario uncertainty [19].
In this study, we examined the projected change to sea surface temperature (SST) variability in the North Pacific and its physical drivers using the Community Earth System Model version 2 Large Ensemble (CESM2-LE), which consists of 100 ensemble member simulations [13].Changes to SST variability are of key importance to both physical and biological components of the climate system: SSTs couple the ocean and atmosphere via radiative and turbulent heat fluxes [20] and control many physiological processes of marine organisms [21].The occurrence of marine heatwaves, prolonged periods of anomalously high SST that result in severe ecological and socioeconomic impacts [22], is directly related to SST variability from a moving baseline perspective [23,24].
Strikingly, the projected change in SST variance in CESM2-LE between 1960-2000 and 2060-2100 is not spatially uniform (Fig. 1c, d), and the aim of this study was to identify the drivers responsible for this pattern of variability change.Note that these projected changes in variance directly translate (if the other statistical moments remain constant) to changes of threshold exceedances of upper percentiles (e.g., the 90 th percentile) that are often used to define marine heatwaves (e.g., ref. [25]).
The area-weighted spatial pattern correlation coefficient between the SST standard deviation change (Fig. 1e) and marine heatwave intensity change (Fig. 1f) globally is 0.87.
We used a local linear stochastic SST model to quantify the relative effect of changes to three drivers on the overall change in SST variance: ocean memory, ENSO teleconnections, and stochastic noise forcing.

Data
We used the Community Earth System Model version 2 Large Ensemble in this study.
CESM2 is a coupled Earth system model with active ocean biogeochemistry [26].The model incorporates the CAM6 atmosphere model and POP2 ocean model, both on ∼ 1 • horizontal grids, as well as coupled land, sea ice, wave, marine biogeochemical, and river runoff models.The large ensemble consists of 100 ensemble members run from 1850 to 2100 and forced by CMIP6 historical (1850-2014) and SSP3-7.0protocols (2015-2100) [13].The SSP3-7.0 scenario, which has a high rate of emissions, was selected to investigate climate variability and its projected future changes.Anomalies were calculated by subtracting the ensemble mean from each ensemble member.We excluded SST data from our analysis at grid points where the ensemble-mean sea ice fraction exceeded 15% for any month during the time period considered.We used SSTs from the Hadley Centre Global Sea Ice and Sea Surface Temperature v1.1 dataset (HadISST [27]); sea level pressure and 850-hPa winds from the ECMWF Reanalysis v5 (ERA5 [28]); mixed layer depth from the Ocean Reanalysis System 5 (ORAS5 [29]), which is defined as the depth where the density exceeds the near surface density by 0.01 kg m −3 ; turbulent surface heat fluxes from the 1 • Objectively Analyzed air-sea Fluxes (OAFLUX [30]); and radiative surface heat fluxes from OAFLUX (derived from the ISCCP-D product [31]) and Clouds and Earth's Radiant Energy Systems Energy Balanced and Filled Ed4.2 product (CERES EBAF [32]).Anomalies were calculated by subtracting the climatology for the entire time period used and then detrending with a linear fit.We excluded HadISST data from our analysis at grid points with sea ice cover (i.e., NaN values in the data) during any month from January 1960 to January 2000.
For the radiative heat fluxes, we calculated anomalies separately for OAFLUX (January 1985 to February 2000) and CERES EBAF (March 2000 to December 2022), and then combined the two sets of anomalies.We spatially smoothed this heat flux data using a moving average filter with 3-by-3-grid-cell window size.For computations requiring both heat flux and SST data, we also spatially smoothed the HadISST data in the same manner.Note that the CESM2-LE data was not smoothed.

Marine Heatwave Intensity
Marine heatwaves were defined using a 90 th -percentile threshold for monthly SST anomaly computed for each calendar month using all ensemble members [33].The mean marine heatwave intensity at a given grid point was calculated as the mean SST anomaly of all 90 th -percentile exceedances over time and all ensemble members.

Linear Stochastic-Deterministic Model
To quantify the effect of different drivers on SST variance, we used an extension of the original local linear stochastic climate model [34,35] with seasonally modulated feedback and noise forcing [36,37] and an ENSO teleconnection term [38][39][40].We use the formulation developed in refs.[41][42][43] that includes seasonal modulations in the feedback, noise forcing, and the ENSO teleconnection term: where T ′ is the SST anomaly at a given location, λ is a seasonally modulated feedback coefficient, β is a seasonally modulated ENSO teleconnection coefficient, N is the Niño3.4index (the SST anomaly averaged over 5 • N-5 • S, 170 • W-120 • W), and ξ is stochastic forcing (i.e., "weather noise").Averaged over the annual cycle, λ must be negative so that SST anomalies are damped and do not grow without bound.λ−1 has units of time and represents the decay timescale of SST anomalies, thus we refer to it hereafter to as the "ocean memory" [44].
The parameters λ and β are defined as where ω a is the angular frequency of the annual cycle (2π/12 months −1 ) and λ 1 , λ 2 , β 1 , and β 2 determine the amplitude and phase of the seasonal modulation.Physically, the seasonal modulation of these coefficients reflects seasonal changes of air-sea heat fluxes and the mixed layer heat capacity, the latter which is proportional to the mixed layer depth [41,45].For ease of display we present these coefficients as annual averages in this report (the amplitude and phase of λ and β are shown in Supplementary Fig. 1).
The noise term ξ primarily represents stochastic forcing from the atmosphere.It includes all processes that are uncorrelated with local SST anomalies and remote ENSO forcing, chiefly anomalous air-sea heat fluxes and anomalous Ekman advection of the SST gradient due to weather variability [46].Entrainment and other ocean processes can also contribute to the forcing [47,48].ξ should be nearly white given the fast decorrelation timescale of the atmosphere [34,49].
At each grid point for each ensemble member, equation 1 was fitted to the SST anomaly data using multiple linear regression (see ref, [42]).∂T ′ /∂t was computed using the forward finite difference method.The noise forcing ξ was taken to be the residual from the fit.This residual is well-described by white noise (see Supplementary Fig. 2), supporting the suitability of our choice of theoretical SST model.

SST Feedback Decomposition
The SST feedback coefficient λ is the sum of several different atmospheric and oceanic feedbacks [48,[50][51][52]: where λSH , λLH , λSW , λLW are the feedbacks associated with the sensible, latent, shortwave, and longwave components of the air-sea heat flux, respectively; λent is the feedback due to entrainment as the mixed layer deepens in fall and winter; λdiff is the feedback due to horizontal eddy diffusion, and λother is the feedback due to non-local and other processes not considered here.
We calculate the air-sea heat flux feedbacks given heat flux component x by fitting the following equation using multiple linear regression: where Q ′ x (t) is the heat flux anomaly (defined as positive downward), λ * x is the feedback for that heat flux component (with units Wm −2 K −1 ), and ξ * x (t) is the noise forcing.

λ *
x is related to the feedbacks λx in equation 4 by the following: where ρ is the density of seawater (1024 kg m −3 ), c p is the heat capacity of seawater (4001 J kg −1 K −1 ), and H is the monthly mixed layer depth climatology.To fit this equation to observations, we used the whole time period available for the heat flux data to minimize the error: January 1985 to December 2022 instead of the 1960-2000 period for fitting equation 1.
The feedback due to entrainment is where went is the entrainment velocity climatology, the time derivative of the mixed layer depth climatology H, and T ′ b is the temperature below the mixed layer, with angled brackets denoting the ensemble/time mean (see ref. [50]).If T ′ b is uncorrelated with T ′ , and assuming a mixed layer of average depth 75 meters with an annual cycle amplitude of 100 meters, λent ≈ −0.1 months −1 when averaged over the annual cycle.
Entrainment also leads to the phenomenon of "reemergence": often the SST anomaly from the previous winter persists under the mixed layer during summer and in fall is re-entrained into the mixed layer, leading to the reemergence of SST anomalies [53,54].Reemergence is not modeled in this work.
The feedback due to horizontal eddy diffusion is where κ is the horizontal eddy diffusivity.Assuming SST anomalies with a sinusoidal spatial structure of wavelength L, the feedback can be estimated via scaling analysis as For L ≈ 1000 km (i.e., a length scale of ∼ 160 km) and κ ≈ 500 m 2 s −1 (note that κ is a function of length scale and geographic location; see ref. [55]), λ diff ≈ −0.05 months −1 .
Equation 4 can be rewritten as where ∆ indicates the change between the two time periods, a subscript 0 indicates that the value from the first time period is used and ∆ λH is the change in the air-sea heat flux feedback due to the change in the mixed layer depth climatology.

Applicability of the Linear Stochastic-Deterministic Model
Equation 1 describes SSTs forced solely by the atmosphere: anomalous air-sea heat fluxes and anomalous Ekman advection of the mean SST gradient from stochastic weather processes and remote forcing from ENSO.Contributions to the variance from internal ocean dynamics (e.g., geostrophic advection, mixed layer depth variability, and entrainment) are neglected [56].This simplification is inadequate to explain SST variance in the equatorial oceans, where coupled ocean-atmosphere dynamics in the Pacific give rise to ENSO; in western boundary currents, where ocean dynamics are important [57][58][59]; and in the areas of the North Atlantic and Southern Ocean where the thermohaline circulation contributes to SST variability on long timescales [60,61].
In previous studies, the applicability of a linear stochastic model to SST dynamics was tested by goodness of fit to a theoretical power spectrum [50,57], by establishing a threshold of sea surface height variance over which oceanic processes were assumed to dominate [62], or by comparing advection of SST anomalies with the estimated feedback term [45].
We used an objective criterion based on the lagged covariance of SST anomalies T ′ and net surface heat flux anomalies Kuroshio-Oyashio Extension region.For observations, as with the calculation of the air-sea heat flux feedbacks, this criteria was evaluated using data from January 1985 to December 2022.Supplementary Fig. 3 shows R T Q at several representative locations.

Isolating SST Variance Contribution from Each Driver
Once λ, β, and ξ are determined, the SST variance due to changes in the corresponding drivers-the ocean memory, ENSO teleconnection, and noise forcing-can be isolated.
We used two forward integrations, one isolating the SST anomalies forced only by the ENSO teleconnection T ′ N and the other isolating SST anomalies forced only by noise where k is the time index, m is the month index (k mod 12), and ∆t is the time step (one month).ξ(k) was constructed using a shuffled fit residual (for each ensemble member): for each calendar month, the year was randomly shuffled, producing noise forcing that is temporally uncorrelated (i.e., white) but retains spatial correlations and seasonal variance modulation present in the fit residual.Our results differ little if the original fit residual (that contains both spatial correlations and a slight temporal autocorrelation) or a version in which the time dimension of the noise forcing is shuffled in a different random order at each grid point (and thus is white in both time and space; see Supplementary Fig. 4).
To isolate the change in variance due to the change of each driver, we performed six of these integrations with parameters from different time periods (see Tab. 1).Each integration was initialized with the SST anomaly at the beginning of the specified time period (2060-2100 for case C).We calculated the change in variance due to the  change in each driver using the following expressions: where ∆ x σ 2 (T ′ ) is the change in SST variance due to changes to the driver x, σ 2 (T ′ x,n ) is the variance of the integrated SST time series corresponding to the case letter n (A, B, or C) in Tab. 1.

Statistical Significance Testing
All parameters shown in this report (e.g., σ 2 (T ′ x,n ), λ, β) were calculated for each ensemble member, creating 100 independent samples.Welch's t-test was then used to assess the statistical significance of ensemble-mean changes of these parameters between 1960-2000 and 2060-2100 [64].Except in areas with minimal changes, the null hypothesis of no change between the two time periods is rejected at the 5% level.

Ocean Memory and Its Future Changes
The ocean memory varies considerably across the North Pacific, both in observations and CESM2.Over most of the North Pacific, the ocean memory diagnosed from the observations is between 2-6 months (Fig. 2a).Equatorward of about 20 • N, particularly toward the eastern side the basin, the ocean memory is substantially longer, typically around 9 months.The magnitude of the ocean memory is largely consistent with previous estimations (e.g., refs.[40,56]) and the autocorrelation timescale of largescale modes such as the the Pacific Decadal Oscillation [39].
In the observations, the contribution of the different heat fluxes to the total feedback (Fig. 3a-c The ocean memory in CESM2-LE is similar in magnitude to observations, ranging between about 2 and 9 months, but has a distinct spatial pattern (Fig. 2d, g).The ocean memory is shorter in the western North Pacific than in the east, which can mostly be attributed to strong damping by turbulent heat fluxes (Fig. 3d).As in the observations, the turbulent and radiative feedbacks are dominated by the latent heat and shortwave feedbacks, respectively (see Supplementary Fig. 5).A large area of particularly long ocean memory is present between Hawai ' i and North America, resulting from relatively weak turbulent heat flux damping and positive radiative feedback, likely from the low cloud-SST feedback.
Interestingly, the phases of λ and H differ: λ is most strongly negative between August and December (depending on location) whereas H is deepest between December and March (see Supplementary Fig. 1).That implies that the seasonality of the air-sea heat flux feedbacks play a strong role in the seasonal modulation of λ in addition to that of the mixed layer depth.
In observations, the residual feedback has considerable spatial structure (Fig. 3c), with areas of negative and strongly positive feedbacks.In CESM2-LE, the residual feedback is negative everywhere except for coastal areas off China and Mexico.As related in Section 2.3, entrainment and horizontal eddy diffusion are expected to damp SST anomalies, with a combined feedback on the order of -0.15 months −1 , which corresponds well with the results from CESM2-LE.However, the strong positive feedbacks in observations could be the result of errors in the heat flux and mixed layer depth data.The magnitude of the feedbacks λ * x for different heat flux components are similar between observations and CESM2-LE (see Supplementary Fig. 5).However, the mixed layer depth is typically somewhat deeper in CESM2-LE than in the ORAS5 reanalysis, which would lead to the λrad and λturb being greater in magnitude in observations compared to CESM2-LE.Part of that discrepancy may be due to the different mixed layer definitions used: a density-based definition for ORAS5 (see Section 2.1) and a buoyancy-based definition for CESM2 [68].
In the future climate in CESM2-LE, the ocean memory declines over most of the basin except for a zonally-elongated area in the central North Pacific where it increases (Fig. 2j).The changes to the individual feedbacks are are spatially varied, but it appears that the change in ocean memory is primarily driven by changes to the radiative and residual feedbacks, suggesting that changes in clouds and ocean dynamics are most important for the change in ocean memory.In common with other climate models (e.g., refs.[44,69]), the mixed layer depth in the North Pacific in CESM2-LE is shallower nearly everywhere in the future climate, leading to a reduced heat capacity and correspondingly shorter ocean memory (Fig. 3g).However, the magnitude of the feedback change due to the shallower mixed layer is relatively minor compared to the changes to the other feedbacks, in contrast with the findings of ref. [44], who attributed the projected decline in ocean memory in CMIP6 models primarily to mixed layer depth shallowing.

ENSO Teleconnection and Its Future Changes
The ENSO teleconnection, represented by β multiplied by the standard deviation of Niño3.4, in both observations and CESM2-LE (Fig. 2b, e, h) exhibits the well-known The spatial pattern of the teleconnection in CESM2-LE for 1960-2000 is broadly similar to the observed pattern but is displaced slightly to the west and is somewhat stronger in magnitude.The westward displacement likely is due to the ENSO SST anomaly in CESM2 extending further west than in observations [73].However, in most of the North Pacific the observed teleconnection falls within two cross-ensemble standard deviations.At the center of action in the central North Pacific, the annuallyaveraged teleconnection coefficient β is much stronger in observations than in CESM2-LE for either time period (see Supplementary Fig. 6).However, the ensemble mean Niño3.4 standard deviation in CESM2-LE is about 50% greater than in observations: 1.30 K and 1.26 K for 1960-2000 and 2060-2100, respectively, compared to the observed value of 0.86 K for 1960-2000 in HadISST.Thus, the overall magnitude of forcing of the teleconnection on SST anomalies is comparable between the model and observations.
In CESM2-LE, the ENSO teleconnection pattern shifts to the northeast in the future climate.The teleconnection, both in its effect on atmospheric circulation and SSTs, weakens slightly.That shift likely is caused by the eastward shift of the location of maximum precipitation during ENSO due to the expansion of the western Pacific warm pool (see refs.[74,75]).Changes to the atmospheric waveguide may also contribute to the teleconnection shift.

Noise Forcing and Its Future Changes
The variance of the noise forcing ξ has a broad maximum at 40  These findings have implications for predictability -the generally lower ocean memory and higher noise forcing suggests that predictability of a simple "damped persistence" model will decline in skill in the future climate in most regions.ENSO is a major source of SST predictability on seasonal timescales, hence the shift of its teleconnections results in ENSO-associated changes in predictability in different regions.Our results highlight the importance of studies into future ENSO changes and its regional impacts.
Although this study was focused narrowly on the North Pacific and the CESM2-LE model, our framework should be equally applicable to other extratropical oceans and other climate models.Different large ensemble climate models show considerable diversity in their future ENSO dynamics [5], thus contribution of the various drivers of SST variability may differ greatly between models.This study also did not determine the physical mechanisms responsible for the change in ocean memory and stochastic noise forcing and how they relate to climate mean state changes.We aim to answer these questions in future work.
[50,56,63]).If SST anomalies are both damped and forced by Q ′ , at negative lags (when the ocean leads), R T Q should be negative, corresponding to damping of SST anomalies by Q ′ .At positive lags (when the atmosphere leads), R T Q should be positive, corresponding to forcing of SST anomalies by Q ′ .Thus we considered that any grid point which had R T Q < 0 at negative lags (averaged over lags -3 to -1 months and all ensemble members) and R T Q > 0 at positive lags (averaged over lags 1 to 3 months and all ensemble members) to be well represented by a linear stochastic model forced by the atmosphere.The grid points that did not meet this criterion were excluded from our analysis and are shown as white hashed areas in the Fig.s.As expected these grid points are in areas of high oceanic variability and strong air-sea coupling, such as the equatorial Pacific and

Fig. 2
Fig. 2 Equation 1 parameters fit to HadISST (a)-(c) and CESM2-LE (d)-(i) SST data in shaded contours, with CESM2-LE projected changes on the bottom row (j)-(l).(b), (e), (h) vectors and contours are the 850-hPa winds and sea level pressure regressed onto the Niño3.4index, the latter with 0.25-hPa/K spacing (positive values are solid lines and negative lines are dashed, with a thicker line at the zero contour).Stippling in (d)-(f) indicates that the parameters derived from observations lie outside the 5 th -95 th percentile range of those derived from the CESM2-LE ensemble members.Stippling in (j)-(l) indicates where the changes are not significant at the 5% level.The ocean memory and ENSO teleconnection panels show the mean over the seasonal cycle, and all CESM2-LE panels are the ensemble mean.Locations where the SST data is not well-described by a local linear stochastic model are shown as white hatched areas (see Section 2.5).
) shows strong damping from turbulent heat fluxes (almost entirely the latent heat feedback) particularly in a band at 25 • N in the western North Pacific.Over much of the North Pacific poleward of 20 • N, the radiative heat flux feedback (almost entirely shortwave feedback) is positive, indicative of the low cloud-SST feedback, where negative SST anomalies are associated with increased atmospheric stability, leading to the formation of low clouds which reduce surface shortwave radiation and further cool the ocean [65-67].

Fig. 3
Fig. 3 (a)-(f) Turbulent, radiative, and residual SST feedbacks in HadISST and CESM2-LE for 1960-2000.(g)-(i) Changes to those feedbacks in CESM2-LE between 1960-2000 and 2060-2100, with (j) showing the contribution of the mixed layer depth change.Stippling in (g)-(i) indicates where the changes are not significant at the 5% level.All panels show the feedbacks averaged over the seasonal cycle and the CESM2-LE panels showing the ensemble mean.Locations where the SST data does not meet the criterion described in Section 2.5 are shown as white hatched areas.

3. 4 4 Conclusions
Fig.4hshows the contribution of each driver to the overall variance change by assigning the change due to each driver to a color channel (red=∆ ξ σ 2 (T ′ ), LW ) heat flux feedback, and λres is the residual feedback.λres includes λent , λdiff , λother , and errors in estimating the air-sea feedbacks.From the estimations above, λent + λdiff ≈ −0.15 months −1 , thus we expect λres to have a similar value if there are not substantial errors in the calculation of the feedbacks and contributions from other unmodeled feedbacks.Because the large number of degrees of freedom in CESM2-LE (100 members) allows for robust statistical estimates of the atmospheric feedbacks, we expect λres to primarily reflect damping by entrainment and diffusion.However, for observations/reanalysis, uncertainties in the heat flux, SST, and mixed layer depth data may compound to produce substantial errors in the calculated feedbacks and thus λres may primarily reflect these errors rather than just damping from oceanic processes.

Table 1
Integration Parameters [77]in both the observations and CESM2-LE, stretching from Japan to about 150 • W (Fig.2c, f, i).This coincides with the subarctic SST front and the North Pacific storm track, thus high atmospheric and oceanic variability in this region is expected.The noise in observations has considerably greater variance than in CESM2-LE even though the SST variance is similar.Because SST variance increases with increasing ocean memory (in an AR-1 process; see ref.[76]), the greater noise variance in observations is compensated by the somewhat shorter ocean memory to yield comparable overall SST variance to CESM2-LE.The future change of the noise forcing variance is spatially heterogeneous in CESM2-LE.Although increasing in most areas, particularly in the eastern North Pacific between Hawai ' i and North America, there are areas in the central and southeastern parts of the basin where noise variance decreases.The strong increase in variance north of Japan is potentially due to a poleward shift of the Kuroshio[77].