Strong Oceanic Forcing on Decadal Surface Temperature Variability Over Global Ocean

Sea surface temperature (SST) variability on decadal timescales has been associated with global and regional climate variability and impacts. The mechanisms that drive decadal SST variability, however, remain highly uncertain. Many previous studies have examined the role of atmospheric variability in driving decadal SST variations. Here we assess the strength of oceanic forcing in driving decadal SST variability in observations and state‐of‐the‐art climate models by analyzing the relationship between surface heat flux and SST. We find a largely similar pattern of decadal oceanic forcing across all ocean basins, characterized by oceanic forcing about twice the strength of the atmospheric forcing in the mid‐ and high latitude regions, but comparable or weaker than the atmospheric forcing in the subtropics. The decadal oceanic forcing is hypothesized to be associated with the wind‐driven oceanic circulation, which is common across all ocean basins.

Part of the reason for the uncertainty of the forcing mechanism of decadal SST variability stems from the lack of comprenhensive and long term subsurface oceanic observations.One strategy that could help resolving this problem was proposed by Bjerknes: inferring oceanic forcing from the surface heat flux (SHF) (downward positive) (Bjerknes, 1964).He hypothesized that the SST-SHF correlation should change from positive for interannual variability to negative for decadal variability, reflecting the change of the dominant driving from atmospheric to oceanic forcing.Gulev et al. (2013) and later work (O'Reilly & Zanna, 2018;R. Zhang, 2017;R. Zhang et al., 2016) provided strong observational evidence supporting the Bjerknes hypothesis in the mid-latitude North Atlantic where the SST-SHF correlation indeed changes from positive to negative with low-pass filtering, implying an active role for oceanic processes (hereafter referred to as "oceanic forcing") in generating decadal SST variability there.In the theoretical framework of stochastic climate model (Cane et al., 2017;Frankignoul & Hasselmann, 1977;Hasselmann, 1976;R. Zhang, 2017), Z. Liu et al. (2023) (hereafter LGD23) recently demonstrated that this correlation sign reversal further requires oceanic forcing with a long persistence time (relative to the SST persistence time), effectively functioning as a red noise forcing.Moreover, they developed a scheme (hereafter LGD scheme) to directly assess the separate magnitudes of decadal oceanic and atmospheric forcing.An application of the LGD scheme to observations and a state-of-the-art climate model showed that decadal oceanic forcing is stronger than atmospheric forcing in the mid-latitude North Atlantic.
Over the global ocean, however, there has been no systematic assessment of the relative strengths of decadal atmospheric and oceanic forcing.Here, we present the first comprehensive assessment of the strenghth of oceanic forcing on decadal SST variability over global oceans in observations and state-of-the-art climate models.We show that decadal oceanic forcing tends to be stronger than atmospheric forcing over the mid-and high latitude ocean (>30°), but is weaker than atmospheric forcing over the subtropical ocean (10-30°).

CMIP6 Model Simulations
We use the model output of a total of 12 CMIP6 models: ACCESS-CM2, CanESM5, CESM2, EC-Earth3, FGOALS-g3, GFDL-CM4, GISS-E2-1-G, IPSL-CM6A-LR, MIROC6, MPI-ESM1-2-LR, MRI-ESM2-0 and NorESM2-LM.These models are selected based on several criterions, mainly, the availability of required variables, sufficient length of control run, latest model version for those present in both CMIP5 and CMIP6.We analyze the historical simulations of the time period 1900-2014 and piControl simulations longer than 500 years.The multi-model mean (MMM) is derived by averaging the specified variable (e.g., correlations, forcing ratio) across CMIP6 models, calculated with equal weight on each simulation.Model data are first interpolated to the 2°× 2°grid points before analyses.All CMIP6 model output are publicly available at https://esgf-node.llnl.gov/search/cmip6/.

Observational Data
The sea surface temperatures are derived from the Hadley Center's Sea Ice and Sea Surface Temperature (HadISST) data (Rayner, 2003) (https://www.metoffice.gov.uk/hadobs/hadisst).Two surface heat flux data sets are used to check the robustness of the analysis.The first heat flux data set is the net surface heat flux (including sensible heat flux, latent heat flux, net shortwave radiation, and net longwave radiation) and the 10 m wind speed, from the NOAA-CIRES twentieth Century Reanalysis V3 (20CR) (Slivinski et al., 2021) (https://www.psl.noaa.gov/data/gridded/data.20thC_ReanV3.html).The second heat flux data is the observational turbulent heat flux (including sensible heat flux and latent heat flux) from OAFlux data spanning the years 1958-2010 (Yu et al., 2008).The ocean surface salinities from the ISHII data version 6.13 (Ishii et al., 2006) (https://rda.ucar.edu/datasets/ds285.3) are also used for our estimation.All data are interpolated to 2°× 2°grid points before analyses.

Qualitative Assessment of Oceanic Forcing in SST-SHF Correlation
We first study the role of ocean forcing across the global ocean qualitatively by analyzing the SST-SHF correlation from 1900 to 2014 using the quasi-observational twentieth Century Reanalysis (20CR) data sets and the historical simulations of the climate models of the latest Coupled Model Intercomparison Project Phase six (CMIP6).All model and observation data are monthly anomalies after the removal of the seasonal cycle, the linear trend, and the ENSO signal (Text S1 in the Supporting Information).
Figures 1a and 1c present the global map of the SST-SHF correlations for the monthly (unfiltered) and decadal (low-pass filtered with 10-year moving average) variability, respectively, in the 20CR reanalysis.In the North Atlantic, the correlation changes from positive monthly correlation to significantly negative decadal correlation over most of the mid-and high latitudes (>30 o N), and this feature also appears robust in the correlations made using the direct observational heat flux (OAFlux) (Figures S1a and S1c in Supporting Information S1).These correlations are consistent with previous studies (Gulev et al., 2013;O'Reilly et al., 2016), implying a dominant role of oceanic forcing for decadal variability, but of atmospheric forcing for interannual variability.The dominant oceanic role is also robust across the CMIP6 models, as shown for MMM (Figure 1d) and each model member (Figure S2 in Supporting Information S1, second column).It is noted that the decadal correlation remains positive at high latitudes after 1958 in both sets of heat flux observations (Figures S1c and S1d in Supporting Information S1) but this positive correlation becomes much weaker if the analysis period is extended to 1900 in 20CR, which is the only heat flux data set available before 1958 (Figure 1c).This change of decadal correlation after 1958 may be suggestive of a more dominant atmospheric forcing on decadal variability post-1958 from external forcing (He et al., 2023), but strong oceanic influence pre-1958 from internal climate variability.One interesting feature of particular notice here is that the positive correlation remains present in the decadal correlation over the subtropics (10 o -30 o N) in observations (Figure 1c; Figure S1c in Supporting Information S1) and, particularly, in CMIP6 models (Figure 1d), implying dominant atmospheric driving there for  S3g in Supporting Information S1) (Mignot, 2004;Mignot & Frankignoul, 2003) and positive heat flux feedback (Figure S3j in Supporting Information S1) (Park et al., 2005).
decadal variability, as for interannual variability.The correlation over the North Atlantic is thus suggestive of a pattern of strong decadal oceanic forcing in the mid-and high latitudes, but weak decadal oceanic forcing in the subtropics.
The implied strong decadal oceanic forcing, most interestingly, appears in observations over the mid-and high latitudes across other ocean basins, including North Pacific, South Pacific, western and central tropical Pacific, and Western India Ocean, for the periods post- 1900 (Figure 1c) and post-1958 (Figures S1c and S1d in Supporting Information S1).This implies a dominance of decadal oceanic forcing over many other regions over the globe.Meanwhile, positive decadal correlation is still present in the central subtropic North Pacific, eastern subtropic South Pacific, eastern tropical Pacific, South Atlantic and Southern Ocean, implying the important role of atmospheric forcing there.The monthly correlation remains positive across most global ocean, reflecting the dominant atmospheric forcing on interannual variability, with some exceptions of negative correlation in midlatitudes of the North Pacific, Southern Ocean and equatorial ocean.The exception is the equatorial ocean region, where the correlation is negative from monthly to decadal variability, consistent with our expectation of dominant oceanic forcing associated with upwelling process there across time scales.In CMIP6, the decadal correlations are predominantly negative in the mid-and high latitudes, but positive in the subtropics, across the globe in the North Pacific, South Pacific, South Indian Ocean, South Atlantic and Southern Ocean, as well as in the North Atlantic (Figure 1d; Figure S2 in Supporting Information S1, second column).There are some discrepancies in decadal correlation patterns between observations (Figure 1c) and models (Figure 1d).The model tends to display a zonally uniform latitudinal band structure, whereas observations show some longitudinal variation.This discrepancy could be attributed, partly, to the noise removal of multi-model ensemble mean in the model.
The largely zonal structure of the correlations in models makes this characteristic pattern more distinct in the zonal mean (Figures 2a and 2b, red curves), with negative decadal correlation in the mid-and high latitudes, but positive decadal correlation in the subtropics, largely symmetric to the equator; relatively, the monthly correlation tends to be positive (or less negative to close to zero) over the entire ocean, except in the equatorial ocean.This characteristic correlation pattern also appears detectable, qualitatively, in observations of 20CR and OAFlux (Figures 2a and 2b, black curves and gray markers).The zonal mean decadal correlations in observations appear negative in the subtropics, instead of positive in models, but this negative decadal correlation tends to be less negative than in the mid-latitude, at least in the Northern Hemisphere, qualitatively consistent with the models.There is also a clear discrepancy between observations and CMIP6 models in the equatorial Pacific: the correlation is more negative across the models than in observations for both interannual and decadal correlations, a point to be returned later.

Quantifying Decadal Ocean Forcing Using LGD Scheme
The correlation analysis is useful for inferring the dominant driving, but only qualitatively.Indeed, the sign reversal of low-passed correlation depends not only on the relative strength of oceanic to atmospheric forcing (Cane et al., 2017), but also on the associated SST feedbacks (R. Zhang, 2017).LGD23 shows that the key forcing parameter for the correlation structure is the effective forcing ratio R = σ o 2 / λ o σ a 2 / λ a , where σ o 2 and σ a 2 are the variances of decadal oceanic and atmospheric forcing, λ o and λ a are the SST feedback damping rates, of the oceanic and atmospheric forcing, respectively (Z.Liu et al., 2023) (Appendix A).The correlation map can be largely understood in terms of R in the red noise stochastic climate theory (LGD23).In the mid-and high latitudes, the effective oceanic forcing is significantly stronger than the atmospheric forcing (R > 1) leading to the negative correlation on decadal time scale in both observations and models (Figures 1e and 1f); In the subtropics, however, effective oceanic forcing is less strong (Figure 1e) or even weaker (Figure 1f, r < 1) than the atmospheric forcing such that the decadal correlation tends to be less negative or slight positive in models and, to some extent, in the 20CR observation.This relation also appears detectable in observations of OAFlux (Figure S1e in Supporting Information S1) and each model member (Figure S2 in Supporting Information S1, third column).
As a further step, we use the LGD scheme to quantify the magnitudes of decadal oceanic and atmospheric forcing σ o and σ a over the global ocean from the monthly data of SST, SHF, and sea surface salinity (SSS) (Appendix A).
The LGD scheme involves the estimation of the SST feedbacks from lagged covariances and the reconstruction of atmospheric/oceanic forcing as a residual from the stochastic climate model.The estimated decadal forcing is rather robust, because of the suppression of sampling error and noises by the low-pass filtering (LGD23).
In the observation of 20CR, the ratio of decadal ocean over atmosphere forcing σ o /σ a is estimated over the global ocean in Figure 1g, with the zonal mean shown in Figure 2c.Over the North Atlantic and North Pacific, decadal ocean forcing is stronger than atmospheric forcing (σ o /σ a > 1) over most of the mid-latitude, by a factor of 1.5-2 in the zonal mean.Regionally, decadal oceanic forcing maximizes in the interior basins, with the amplitude over 3 times stronger than the atmospheric forcing in the regions of the central North Atlantic and central North Pacific (σ o /σ a > 3).The forcing ratio further decreases toward the subtropics, implying an relatively weaker oceanic forcing in the subtropics.In the subtropical North Atlantic and North Pacific (Figure 1g), there are regions where the atmospheric forcing is modestly stronger than the oceanic forcing (0.5 < σ o /σ a < 1), such that, in the zonal mean (Figure 2c), the oceanic forcing becomes weaker, approximately 1-2 times the atmospheric forcing.The relatively weaker ocean forcing in the subtropics than mid-latitude is qualitatively consistent with that using the OAFlux observation post-1958 (Figure 2, circle markers) and 20CR post-1958 (Figure 2, crossing markers).
The estimated forcing ratio σ o /σ a in MMM (Figure 1h) and each model member (Figure S2 in Supporting Information S1, forth column) show a pattern somewhat similar to observations (Figure 1g), characterized by a larger forcing ratio in the mid-to high latitudes compared to the subtropics in the Northern Hemisphere.The alignment between CMIP6 models and observation of 20CR becomes more apparent in the zonal mean (Figure 2c), with the observed forcing ratio falls within the model's spread across most latitudes.Additionally, we note that decadal ocean forcing is about 1.5-2 times greater than the atmospheric forcing in the mid-latitude North Atlantic and North Pacific (Figure 2c), and less than the atmospheric forcing in the subtropics (σ o /σ a < 1).Over the Southern Hemisphere oceans, the pattern of forcing ratio in the CMIP6 models is largely similar to the Northern Hemisphere, symmetric to the equator.In observations, however, the forcing ratio appears less coherent with no clear latitudinal structure.Overall, the stronger oceanic forcing in the mid-latitude than in the subtropics is consistent qualitatively with the expectation from the correlation analysis above, which shows a more pronounced negative correlation in the mid-latitude than in the subtropics on decadal timescale.This result is robust regardless of the filtering of ENSO, the inclusion of radiative forcing, the length of the low-pass window and the spatial scale of interest (not shown).
Quantitatively, there are discrepancy between observations and models.Overall, models tend to show a smaller decadal forcing ratio σ o /σ a than in observations in the subtropics, but a higher forcing ratio in the high-latitude and equatorial regions (Figure 2c).These discrepancies may be caused partly by uncertainties in observations, especially in the Southern Hemisphere where historical observations are sparse.More likely, these discrepancies are caused by model biases in the forcing mechanism of decadal climate variability.Models may have overestimated oceanic forcing or underestimated the atmospheric forcing in the subtropics, but the opposite in the high latitude and equatorial region.For example, the systematic overestimation of decadal oceanic forcing in the equatorial region may be related to the well-known cold tongue bias across almost all climate models so far (Li & Xie, 2014;Lin, 2007;T. Liu et al., 2023), which is consistent with an excessive oceanic forcing, and excessive equatorial cold tongue, in the equatorial Pacific.Much further studies are needed to fully understand these model biaes.

Mechanism of Decadal Ocean Forcing
Our assessment of strong decadal oceanic forcing over the global ocean has an important implication to the mechanism of decadal oceanic variability.Previous debate on the role of decadal oceanic forcing has focused on the North Atlantic region, notably for the AMV (Cane et al., 2017;Clement et al., 2015Clement et al., , 2016; R. Zhang, 2017; R. Zhang et al., 2016).Those studies that favor oceanic forcing of the AMV often associated the decadal oceanic forcing with the Atlantic Meridional Overturing Circulation (AMOC) (Gulev et al., 2013;O'Reilly et al., 2016;O'Reilly & Zanna, 2018;R. Zhang et al., 2016R. Zhang et al., , 2019)).Since the AMOC is unique in the Atlantic, the associated oceanic forcing should exhibit some distinct features from other basins, if it is indeed the major decadal forcing in the Atlantic.However, our assessed decadal oceanic forcing shows a similar pattern across all ocean basins: strong ocean forcing in the mid-latitude, but weak in the subtropics.This common forcing pattern implies a common mechanism for decadal SST variability across ocean basins and in both hemispheres, excluding AMOC as the major mechanism.
Given the roughly zonal pattern of the decadal ocean forcing in all ocean basins, most clearly in climate models, we speculate this decadal oceanic forcing is associated with the wind-driven circulation, which is characterized by subtropical and subpolar gyres that are indeed largely similar across different ocean basins.This decadal oceanic forcing, however, is unlikely to be simply associated with the temperature advection by the surface Ekman flow Q ek (Appendix C).This follows that the sign reversal of SST-SHF correlation must be driven by a red noise ocean forcing of persistence time comparable or longer than that of SST (LGD23).The Q ek , however, is a fast forcing such that it is unable to drive the sign reversal of SST-SHF correlation.Indeed, the Ekman flow is an instantaneous response to surface atmospheric wind (after geostrophic adjustment in hours) such that it should have a time scale comparable with atmospheric variability and SHF.This is confirmed in the autocorrelation analysis of SST, SHF, and Q ek over the globe in 20CR and models (Figure S4 in Supporting Information S1).The damping rate of Ekman advection ranges from 1 to 2 (month) 1 , comparable with that of the SHF of about 2 (month) 1 , but significantly greater than the SST damping rate of approximately 0.3 (month) 1 .
Given the comparable time scale of SHF and Q ek , the two forcing can be combined as the effective (fast) surface forcing SHF + Q ek , whose covariance with SST <T,H + Q ek > can be compared with that of SHF covariance <T, H>.In both the 20CR and CMIP6, there is little difference of low-passed covariance between <T,H> and <T, H + Q ek > in the zonal mean (Figures 3c and 3d) and regionally over the globe (Figures S5e-S5h in Supporting Information S1).This confirms the negligible role of Ekman advection on SST decadal variability.Meanwhile, the Ekman advection does contribute systematically to the unfiltered covariance in both the 20CR and CMIP6 models: <T, H + Q ek > is smaller than <T,H> in the subtropical region, but larger than <T, H> in the mid-latitude region (Figures 3a and 3b; Figure S5a-S5d in Supporting Information S1).This suggests that the Ekman advection opposes the SHF in the easterly wind region, but reinforces the SHF in the mid-latitude westerly region.

Geophysical Research Letters
10.1029/2023GL107401 The different impacts of Ekman advection on SHF between the easterly and westerly regions is caused by the cancellation/reinforcement of the meridional Ekman advection with the total wind speed and, in turn, turbulent heat fluxes, as studied previously for the North Pacific (Lee et al., 2008) and North Atlantic (Buckley et al., 2014).
This implied decadal ocean forcing is unlikely to be caused by external forcing either.Recent studies have also suggested that multidecadal SST variability can be influenced by external forcing factors such as the volcanic forcing and anthropogenic aerosol (Bellomo et al., 2018;Booth et al., 2012;Mann et al., 2021;Murphy et al., 2021;Otterå et al., 2010).These external forcing factors, however, are not the major contributors to the decadal oceanic forcing, at least, in the CMIP6 models.This can be seen by comparing the correlations and forcing ratio between the historical simulations and the preindustrial control simulations (piControl) of the same model.Unlike the historical simulations that are forced by variable external forcing, the piControl is forced by constant external forcing such that all variability is caused by internal coupled ocean-atmosphere variability.The SST-SHF correlation remains almost identical for monthly correlation and only differs modestly on the decadal correlation, leaving the decadal forcing ratio σ o /σ a almost the same between the historical and piControl simulations (Figure S6 in Supporting Information S1).we speculate, the most likely mechanism for the decadal oceanic forcing is the temperature advection associated with geostrophic oceanic current that is forced by variable wind-driven Ekman pumping (T.Delworth et al., 1993;Frankignoul et al., 2002;Jin, 1997;Latif & Barnett, 1994;Marshall et al., 2001).A further study along this direction will be presented elsewhere.

Conclusion and Discussions
This study presents the first quantitative estimation of decadal oceanic and atmospheric forcing on SST variability across the global ocean in observations and state-of-the-art climate models.In the mid-and high latitude, decadal ocean forcing is stronger than the atmospheric forcing by about 2-3 times.In the subtropics, however, decadal oceanic forcing tends to be comparable or even weaker than atmospheric forcing.Importantly, this decadal ocean forcing tends to exhibit a common pattern in all ocean basins, suggestive of the role of wind-driven circulation.
Ocean dynamics is known to affect SST variability strongly in the region of strong wind-driven ocean currents from interannual to longer time scales, notably the equatorial upwelling zone (Philander, 1989) and the region of intensified western boundary currents (Bishop et al., 2017).However, it has remained unclear how the winddriven circulation affect decadal SST variability across most of the interior ocean.We speculate that the most likely mechanism in the interior ocean could be attributed to the non-Ekman component of wind-driven circulation, that is, the geostrophic current, via entrainment to affect the surface variability.Furthermore, given the nature of multiple decadal and interdecadal variability in the Atlantic (Kushnir, 1994;Msadek & Frankignoul, 2009;Sutton & Hodson, 2005;R. Zhang et al., 2019;R. Zhang & Delworth, 2006) and Pacific (Y.Zhang et al., 1997;Mantua et al., 1997;Deser et al., 2004;Newman et al., 2016;Stolpe et al., 2017;Di Lorenzo et al., 2008), our study does not exclude the significant role of additional mechanisms of ocean forcing on decadal variability, notably the AMOC.Overall, our study highlights the key role of oceanic processes in the generation of decadal variability across the globe ocean.Dynamically, the specific nature of the implied decadal oceanic forcing remains to be further explored in the future.

Appendix A: LGD23 Theory and Scheme for Forcing Estimation
For convenience, we briefly describe some relevant results of the stochastic climate theory of LGD23 (Z.Liu et al., 2023), and the LGD scheme of forcing estimation.The stochastic climate model takes the form of (Cane et al., 2017;R. Zhang, 2017) where the net atmospheric forcing H(t) and net oceanic forcing O(t) are: The λ a and λ o are the damping rates associated with SST feedback, combined giving the total SST damping rate b = λ a +λ o .The f a (t) is a white noise forcing representing synoptic atmospheric variability (Hasselmann, 1976).The f o (t) is introduced in LGD23 as a red noise forcing of damping rate d, reflecting slow subsurface oceanic processes as simulated in where k(t) is a white noise forcing, which is assumed to reflect faster oceanic processes such as convection, mixing, wave braking, entrainment and small scale eddies.
LGD23 showed that one key forcing element that determines the structure of the SST-SHF correlation r TH is the effective forcing ratio = σ o 2 / λ o σ a 2 / λ a , where σ 2 o = < f o , f o > and σ 2 a = < f a , f a > are respectively the variances of the lowpass oceanic forcing f o and atmospheric forcing f a in the limit of long low-pass window length L → ∞.The other element is the damping rate of the red noise oceanic forcing d.These two forcing parameters determine the correlation in three regimes.For strong atmospheric forcing of R ≤ 1, the correlation is positive r TH > 0 across all time scales, regardless of low-pass filtering.On the other end, for strong and slow oceanic forcing such that R > 1 + d 1 , the correlation is always negative r TH < 0 regardless of time filtering.For intermediate strength of oceanic forcing, such that 1 + d 1 > R > 1, the correlation r TH reverses the sign from positive to negative with increasing low-pass window .It should be noted that the effective forcing ratio R depends not only on the relative magnitude of the forcing, but also the relative damping rate of the forcing.
LGD23 further developed a residual-low pass scheme (LGD scheme) to estimate the magnitudes of decadal oceanic and atmospheric forcing from monthly heat flux, SST and SSS in observations or climate models.The LGD scheme proceeds in five steps Step 1: Estimate damping rates b, λ o , λ a (Figures S3a-S3i in Supporting Information S1): The total SST damping rate b is estimated from the SST lag-1 autocorrelation coefficient, while the ocean forcing damping on SST λ o is estimated from the sea surface salinity lag-1 autocorrelation coefficient (Hall & Manabe, 1997;R. Zhang, 2017) such that

Geophysical Research Letters
Step 2: Estimate the heat flux H in the dimension of temperature tendency ( o C s 1 ): H is estimated from the dimensional heat flux H (in W m 2 ) as H = H/ (ρC p h e ) .Here, h e is an effective mixed layer depth estimated indirectly from the air-sea feedback as follows.In the original heat flux dimension (W m 2 ), the equivalent form of Equation A2 is H = μ A T + f a , where the feedback coefficient μ a is related to the atmospheric damping rate on SST as λ a = μ A /(ρC p h e ), with ρ = 1024 kg m 3 and C p = 3,850 Jkg 1 K 1 as the density and specific heat of seawater, respectively.Thus, h e is estimated as (Figures S3m-S3o in Supporting Information S1): Here, μ A is estimated using the Equilibrium Feedback Analysis (EFA) (Frankignoul et al., 1998;Z. Liu et al., 2008) as the ratio of the covariance of H lagging SST as (Figures S3j-S3l in Supporting Information S1): Step 3: Estimate atmospheric forcing f a : the atmospheric forcing is derived as the sum Step 4: Estimate oceanic forcing f o : The ocean forcing is estimated from the SST equation 1 as the residual with T t being estimated using the monthly SST in finite difference.
Step 5: Estimating the variance ratio of decadal ocean/atmospheric forcing in Equations A8 and A9 at the lowpass limit σ 2 o / σ 2 a = < f o , f o >/ < f a , f a > at L → ∞.In this work, the low-pass window is 10 years, leaving σ 2 o and σ 2 a as the variance of decadal forcing.
There is one more step in which the damping rate of the oceanic forcing d is estimated from the reversal of lowpassed correlation.This step is not shown here because it is not used here.Note that, in the LGD scheme, the variance ratio σ 2 o /σ 2 a is estimated before the estimation of the damping rate of the red noise oceanic forcing d.The former estimation is generally more robust than the latter estimation (Z.Liu et al., 2023).
It should also be pointed out that the LGD scheme is invalid in some regions, especially in observations.The failure can be induced by a greater estimated λ o than b, leading to negative λ a (Figures S3g and S3h in Supporting Information S1), or EFA estimation of local positive heat flux feedback (Figures S3j and S3k in Supporting Information S1).These may involve dominant non-local process in changing SSS and SST (Mignot, 2004;Mignot & Frankignoul, 2003).

Appendix C: Estimation of Ekman Advection
According to the upper ocean heat budget, we have where ρ = 1024 kg m 3 is the density of seawater, C p = 3,850 Jkg 1 K 1 is the specific heat of seawater, h is the mixed layer depth (here we use the effective mixed layer depth in Equation A6, H is the dimensional heat flux (in W m 2 ), and Q ocn is oceanic heat transport.
We decompose oceanic heat transport Q ocn into Ekman part Q Ek and non-Ekman part Q nonEk .
In CMIP6 models, the Ekman advection is calculated by monthly wind stress: where ρ = 1024 kg m 3 is the density of seawater, C p = 3,850 Jkg 1 K 1 is the specific heat of seawater, h ek is the Ekman layer depth, u E and v E are the meridional Ekman velocity in each layer, U E and V E are the Ekman transport.
Here we make an approximation that Ekman layer depth is shallower than the mixed layer depth, so that the SST gradient is independent of depth.
The Ekman transport is given by: where τ x and τ y are zonal and meridional wind stress, f is the Coriolis parameter.
Due to unavailability of wind stress output in the 20CR data set, we utilize wind speed at 10 m to calculate the wind stress by the bulk formula (Large & Pond, 1981).where u 10 and v 10 are the wind speed at 10 m above the sea surface, ρ air = 1.23 kg/m 3 is the air density, C D = 0.0013 is the dimensionless drag coefficient.

Figure 1 .
Figure 1.Global SST-SHF correlations and forcing ratio for 1900-2014 (a), (b) Unfiltered SST-SHF correlations (rTH) (c), (d) Low-pass SST-SHF correlations (e), (f) Effective forcing ratio (g), (h) Low-pass forcing ratio.The left and right columns respectively are 20CR and the multi-model mean of CMIP6 historical simulations.In (a) and (c), yellow hatched area and purple dots indicate the statistical significance at the 85% and 95% confidence levels, respectively, based on Student's t-test using the effective degrees of freedom (Confidence levels for individual models are shown in Figure S2 in Supporting Information S1).Dashed black lines represent the boundary of tropical regions which are not the focus of this study.The gray shading in (e) and (g) represent missing values associated with nonlocal effects (FigureS3gin Supporting Information S1)(Mignot, 2004;Mignot & Frankignoul, 2003) and positive heat flux feedback (FigureS3jin Supporting Information S1)(Park et al., 2005).

Figure 2 .
Figure 2. Global zonal mean of SST-SHF correlations and forcing ratio.(a) Unfiltered SST-SHF correlations (rTH).(b) Lowpass SST-SHF correlations.(c) Low-pass forcing ratio.The black solid lines represent the 20CR from 1900 to 2014.The thick red lines and thin pink lines represent the multi-model mean of CMIP6 historical simulations and individual model.The gray cross and circle markers respectively represent the 20CR and OAFlux from 1958 to 2010.The dashed gray lines mark the tropic regions, which are not the focus of this study.

Figure 3 .
Figure 3. Global zonal mean covariance between SST, SHF and Ekman advection Q ek (a), (b) Unfiltered covariances (c), (d) Low-pass covariances.The left and right columns respectively are 20CR and the multi-model mean of CMIP6 historical simulations.The blue solid lines represent the covariance between SST and SHF.The red solid lines represent the covariance between SST and SHF + Q ek .The thin blue and pink lines in (b) and (d) represent the individual models in CMIP6.The dashed gray lines represent the boundary of tropic regions which are not the focus of this study.