1High Predictability of Tropical Pacific Decadal Variability Dominated by

2

Oceanic Rossby Waves

3
4

Xian Wu1*, Stephen G. Yeager2, Clara Deser2, Antonietta Capotondi3,4, Andrew T. Wittenberg5,5and Michael J. McPhaden6

6
7

1 Atmospheric and Oceanic Sciences Program, Princeton University, Princeton, NJ

8

2 National Center for Atmospheric Research, Boulder, CO

9

3 Cooperative Institute for Research in Environmental Sciences, University of Colorado

10

Boulder, Boulder, CO

11

4 National Oceanic and Atmospheric Administration /Physical Sciences Laboratory, Boulder,

12

CO

13

5 National Oceanic and Atmospheric Administration/Geophysical Fluid Dynamics Laboratory,

14

Princeton, NJ

15

6 National Oceanic and Atmospheric Administration /Pacific Marine Environmental Laboratory,

16

Seattle, WA

17
18
19
20
21
22
23
24
25
26
27

*Corresponding author: Xian Wu, xw2794@princeton.edu

128

ABSTRACT

29

Despite its pronounced global impacts, tropical Pacific decadal variability (TPDV) is poorly

30

predicted by current climate models due to model deficiencies and a limited understanding of its

31

underlying mechanisms. Using observational data and a hierarchy of model simulations including

32

decadal hindcasts, we find that decadal isopycnal depth variability driven by oceanic Rossby

33

waves in the tropical Pacific provides the most important source of predictability for TPDV. The

34

predictability arising from initial isopycnal depth conditions is further amplified throughout

35

decadal predictions by tropical ocean-atmosphere coupling and variations in the strength of36

subtropical-tropical cells in the Pacific. Regional initialization experiments that effectively isolate

37

the impact of different ocean basins on TPDV predictability highlight the essential role of the38

tropical Pacific. This study enhances our understanding of the mechanisms governing TPDV

39

predictability, offering crucial insights for improving the accuracy of decadal predictions.

240

Introduction

41

Decadal variations of sea surface temperature (SST) in the tropical Pacific can affect global

42hydroclimate and marine ecosystems 1,2, modulate global mean surface temperature changes 3–5,

43

and interact with the El Niño-Southern Oscillation (ENSO) phenomenon, the leading mode of

44interannual climate variability6–10. However, tropical Pacific decadal SST variations are poorly

45

predicted by the Coupled Model Intercomparison Project Phase 5/6 (CMIP5/6) decadal

46

retrospective forecasts, especially the internal tropical Pacific decadal variability (TPDV)

47

associated with ocean initialization (after removing the effect from external forcings) 11–13. This48

low skill in the Pacific sector contrasts with the high skill for SST in most regions of the Indian49

and Atlantic Oceans, which has been attributed to the response to external forcing and/or ocean

50

initialization11–14.

51

The accuracy of decadal predictions of TPDV relies on the potential predictability provided by52

oceanic processes or external forcings, model representations of these mechanisms, and the realism

53

of oceanic state estimates used to initialize the decadal forecasts. Securing these conditions is54

challenging due to the complex processes that could affect TPDV10,15–18, and systematic model

55

biases in simulating the climatology, variability, and forced changes in the tropical Pacific, as well56

as their interactions with other ocean basins 19–23. Uncertainties in mechanistic understanding and57

model biases of TPDV are challenging to reduce, given limited observational data, particularly for

58

oceanic fields. Therefore, it remains difficult to determine to what extent the low decadal

59

prediction skill for Pacific SSTs is due to intrinsic limits or deficiencies of forecast systems.

60

While numerous studies have investigated the complex origins and mechanisms of TPDV, less61

attention has been devoted to understanding whether and how these mechanisms provide sources62

of prediction skill in retrospective forecasts. The null hypothesis for TPDV is that it is a residual63

of ENSO decadal changes. The leading Empirical Orthogonal Function (EOF) mode of TPDV

64

shows a basin-wide ENSO-like spatial anomaly pattern (Fig. 1a) and is related to random changes

65

in the relative number of warm (El Niño) and cold (La Niña) events over different epochs8,24. The

66

second EOF mode of TPDV displays a zonal dipole pattern in the tropical Pacific (Fig. S1a) and67

is associated with decadal modulation of ENSO amplitude or asymmetries6,25. In contrast to these68

ENSO residual explanations, other studies suggest an extratropical contribution to TPDV from the

69

North or South Pacific26–30. In the extratropics, stochastic atmospheric variability can be integrated3

70

by the ocean due to its large thermal inertia, producing low-frequency SST variability31–34. The

71

resulting extratropical low-frequency SST variability can then influence the tropical Pacific via

72

thermodynamic and dynamical processes, particularly wind-evaporation-SST (WES) and low

73

cloud-SST feedbacks, which propagate wind stress and SST anomalies associated with the Pacific74

meridional mode into the equatorial western-to-central Pacific29,30,35,36.

75

Although the ENSO residual effect and stochastic atmospheric variability lack preferred

76

timescales and are inherently unpredictable on decadal timescales37, these random processes can77

initiate slow oceanic processes which likely determine the timescale of TPDV and provide a source78

of predictability10,17,18,38. Based on observational and/or modeling studies, several oceanic

79

mechanisms in the Pacific have been proposed to contribute to TPDV, including off-equatorial

80

oceanic Rossby wave activity, spiciness advection, and variations in the strength of the subtropical-

81

tropical cells (STCs). Decadal-scale off-equatorial oceanic Rossby wave reflections at the western82

boundary of the Pacific serve as a delayed negative feedback for TPDV by affecting equatorial

83

Pacific thermocline depth, similar to the dynamics driving ENSO phase transitions on interannual

84

timescales39–44. STCs are the upper-ocean overturning circulations connecting the subtropical and85

equatorial Pacific oceans, which influence tropical Pacific SSTs through mean advection of86

temperature anomalies (𝑣̅ 𝑇′) or variations of STC strength (𝑣′𝑇% ). Surface water masses in the87

subtropics that are subducted into the pycnocline may move equatorward and upwell to the surface

88

upon reaching the equator, affecting the equatorial Pacific SSTs45. However, subsequent studies

89

suggest that this subtropical thermal subduction cannot efficiently reach the equator due to energy

90

dissipation, dispersion in the form of planetary-scale oceanic waves, and perturbation from winds

91

at lower latitudes46. Density-compensated temperature anomalies, known as ocean “spiciness”47,92

can propagate more effectively along isopycnal surfaces from the subtropics to the tropics48–51.

93

Alternatively, variations in STC strength can affect the rate of transport of the relatively constant94

water masses, influencing the equatorial Pacific upwelling, with enhanced upwelling bringing

95

colder subsurface water to the surface and reduced upwelling having the opposite effect 46,52,53.

96

The extent to which these various slow oceanic processes contribute to the predictability of TPDV

97

and their relative importance remains unclear.

98

Other ocean basins may also influence tropical Pacific decadal climate through both99

atmospheric and oceanic pathways10,54. For example, decadal-scale cooling in the central-eastern

100

tropical Pacific, a main factor driving global surface warming slowdowns, has been linked to4

101

tropical Atlantic SST warming23,55–59. However, some studies suggest that the tropical Atlantic’s

102

impact on the tropical Pacific might be overstated in regional SST-restoring experiments, which103

can overestimate the upward net surface heat fluxes in the tropical Atlantic compared to104

observations60,61. Further, the Atlantic-Pacific connection may be artificially amplified by the way

105

in which internal variability is defined, leading to spurious linkages62. Investigations using an

106empirical model of the tropical Pacific-Atlantic systems, where the inter-basin feedbacks can

107

cleanly be removed, indicate that Atlantic-Pacific coupling damps TPDV63. Other studies also

108

argue for the role of the tropical Indian Ocean in affecting tropical Pacific decadal climate64,65. It

109

is plausible that inter-basin interactions could play a role in TPDV predictions, but these linkages

110and causalities need to be verified with appropriate sensitivity experiments.

111

Despite the low real-world prediction skill, earlier studies based on perfect model experiments

112

suggest potential multiyear predictability in the Pacific43,66. Recent studies underscore the113

degradation of tropical Pacific decadal prediction skill by volcanic eruptions due to inadequate

114

model representation of volcanic forcing and response, while forecast systems that exclude

115

volcanic eruptions show high skill67–70. Case studies indicate that particular phase transitions of116

TPDV can be retrospectively predicted when a model is properly initialized from strong ENSO

117

event71,72. Other studies indicate that spurious ENSO conditions due to initialization shock in the118

first forecast year will degrade the long-range skill of decadal hindcasts73,74. Increasing the spatial

119

resolution of decadal forecasts to resolve ocean eddies also improves the prediction skill in the120

eastern tropical Pacific because of a more realistic prediction of the SST trend in the Southern

121

Ocean and associated teleconnections75. These findings collectively imply the presence of122

potential predictability for TPDV, which could arise from internal climate processes both within123

and beyond the Pacific.

124

Several important questions remain. What are the external (greenhouse gas, aerosols, and125

volcanic eruptions) and internal (oceanic initialization) factors affecting decadal prediction skill in126

the tropical Pacific? For the skill arising from internal oceanic initial conditions, what is the relative

127

importance of the various mechanisms discussed above? Do other ocean basins play a role in128

affecting TPDV predictions? To answer these questions, we make use of observational and129

reanalysis data and a hierarchy of model simulations conducted with the Community Earth System

130

Model version 1 (CESM1) at 1° horizontal resolution. In particular, we first evaluate the CESM1

131

simulation of TPDV in the free-running coupled preindustrial control simulation and in historical5

132

ocean-only simulations forced with observed surface conditions. We investigate both external and133

internal influences on TPDV prediction skill by comparing uninitialized historical simulations and134

initialized decadal retrospective forecasts. We focus on the initialized hindcasts conducted without

135

volcanic forcing as these exhibit high prediction skill in the tropical Pacific69, and attribute the136

predictability arising from initialization to specific oceanic processes. Finally, the impacts of137

different ocean basins on TPDV are established using a set of regional initialization experiments.

138

Results139

a. Predictability and prediction skill of TPDV140

We evaluate how observed TPDV during 1955–2022 is reproduced by a series of CESM1

141

simulations with varying levels of observational constraints (Fig. 1). We define TPDV as the142

leading EOF of quadratically-detrended 10-year running mean SST anomalies in the tropical143

Pacific (20°S–20°N; 120°E–80°W). This 10-year running mean is used to facilitate comparisons

144

between observations and initialized forecasts as explained below, and it yields results similar to145

that using an 8–40-year band-pass filter10). During 1955–2022, observed TPDV exhibits basin-

146

wide SST anomalies in the tropical Pacific (Fig. 1a), fluctuating between warm and cold phases

147

on decadal timescales (Fig. 1b). During the positive phase of TPDV, the thermocline becomes

148

deeper in the eastern equatorial Pacific but shallower in the western part (Fig. 1c). The observed

149

temporal and spatial features of TPDV are well reproduced by the ocean-sea ice simulation forced

150

with observed surface forcing and fluxes (FOSI; Methods). Decadal Prediction Large Ensemble

151(DPLE) and DPLE without volcanic forcing (DPLE_NoVolc; Methods) are initialized from the152

FOSI oceanic and sea ice states on Nov 1 of each year during 1954–2015, and the ensemble-mean

153

hindcasts are averaged over forecast year (FY) 1–10 for each initialization date, in order to examine154

the predictability of decadal anomalies. At FY1–10, DPLE_NoVolc shows a correlation skill

155[anomaly correlation coefficient (ACC) = 0.57) in predicting the observed (standardized) first156

principal component (PC1) timeseries of TPDV during 1955–2022, which is substantially higher

157

than that for DPLE (ACC = 0.22). The 0.57 PC1 correlation between DPLE_NoVolc and158

observations indicates that over 30% of the observed PC1 variance is predictable with ocean

159

initialization. The lower skill of DPLE compared to DPLE_NoVolc is related to an excessive

160

tropical Pacific cooling response to large volcanic eruptions in the 1960s and 1980s69. In contrast

161

to the initialized decadal forecasts, the ensemble mean of the uninitialized CESM1 large ensemble
6162(LE; Methods), which represents the model’s response to external forcings, shows a negative

163

correlation (-0.33) with the observed TPDV PC1 timeseries, and the magnitude of the externally164

forced TPDV PC1 is much weaker than TPDV PC1 in both the observations and initialized165

forecasts. This contrast between initialized and uninitialized simulations suggests that the observed

166

TPDV during the last ~70 years is largely driven by internal climate variability rather than external167

radiative forcing. We also present the same analysis for the second EOF mode in Fig. S1, which

168

exhibits a zonal-dipole SST pattern over the tropical Pacific in observations and accounts for only

169

16% of the total decadal variance. The second PC (PC2) is linked to the decadal modulation of

170

ENSO characteristics and is poorly predicted by the decadal forecasts and the uninitialized171

simulations (Fig. S1h–j).

172
173
174
175

Fig. 1 Tropical Pacific decadal variability from 1955 to 2022 in observations and CESM1
simulations/hindcasts. (a–e) Global SST and (k–o) equatorial Pacific subsurface ocean
temperature patterns associated with TPDV in (a, f, and k) observations (SST in ERSSTv5 and

7

176
177
178
179
180
181
182
183
184
185
186
187
188

ocean temperature in EN4), (b, g, and l) the forced ocean-sea ice simulation (FOSI), and theensemble-means of (c, h, and m) CESM1 DPLE, (d, I, and n) DPLE_NoVolc, and (e, j, and o) LE.
The EOF analysis is conducted for quadratically detrended and 10-yr running-mean SST
anomalies in the tropical Pacific (20°S–20°N, 120°E–80°W; outlined by the black box in a. to
calculate (f–j) the timeseries of the standardized first principal component (PC1). The numbers in
the top-right corner in a–e denote the percentage of total decadal variance explained by the leading
EOF mode in each dataset. The associated pattern is displayed using the regression of quadraticallydetrended and 10-yr running-mean (a–e) SST anomalies and (k–o) equatorial Pacific (3°S–3°N)
subsurface ocean temperature anomalies onto PC1. The abscissa in f–j represents the start year of
the 10-yr averaging window (e.g., 1955 represents 10-yr average anomalies spanning from 1955
to 1964 for observation, FOSI, and LE, and corresponds to the hindcasts averaged across FY1–10
initialized in Nov 1954 for DPLE and DPLE_NoVolc). The r value indicates the correlations ofeach model simulation with observed PC1.

189

The ensemble-mean decadal hindcasts show weaker amplitude and excessive westward

190

extension of TPDV SST anomalies compared to observations (Fig. 1a–d), due to the predictability

191

limit and/or inherent model bias. The magnitude of the global SST regression onto the standardized

192

PC1 in the decadal hindcasts is weaker than observed for most of the ocean basins, mostly because193

the ensemble averaging reduces variance by isolating the predictable component of variance76. The

194

free-running simulation of CESM1 under preindustrial conditions shows an amplitude of TPDV

195

natural variability more comparable to that in observations (Fig. S2). The pattern bias of TPDV,

196

however, exists even in the preindustrial control simulation (Fig. S2), suggesting that it arises from197

inherent biases in CESM177. We also evaluate the SST and subsurface ocean temperature

198

variability associated with TPDV in other observational or reanalysis datasets. ERSSTv5 and199

HadISST give similar results for SST, while EN4 and ORAS4 show differences in the amplitude

200

and structure of equatorial Pacific subsurface temperature anomalies (cf. Fig. 1 and S2). In the201

following sections, we will present the results based on FOSI in the main text, and those based on

202

observational datasets in the supplementary materials.

203

Consistent with the spatiotemporal features of TPDV captured by EOF analysis in the different

204

simulations (Fig. 1), the anomaly correlation skill (ACC) for detrended decadal SST anomalies in205

the tropical Pacific is increased by including ocean initialization (DPLE minus LE; Fig. 2e and j)

206

or excluding volcanic forcing (DPLE_NoVolc minus DPLE; Fig. 2d and i). The ACC results are

207

overall similar whether verifying against FOSI (Fig. 2) or observations (ERSSTv5 and EN4; Fig.

208

S3). The increased skill in DPLE_NoVolc compared to DPLE is most pronounced over the central

209

tropical Pacific, where the solar radiation reduction associated with volcanic aerosol forcing can

8

210

most effectively influence the ocean mixed layer heat budget and SST69. The skill enhancement211

by ocean initialization in DPLE relative to LE is confined to the eastern Pacific, due to intrinsic

212

model bias in simulating excessive westward extension of TPDV SST anomalies. The local

213

correlation skill for SST in the central-eastern Pacific reaches ~0.7 in DPLE_NoVolc, higher than214

that for EOF PC1 (0.57; Fig. 1).

215

We further investigate the ACC skill of ocean temperature in the equatorial Pacific as a216function of depth and longitude (Fig. 2). The significant impact of volcanic forcing is limited to217

the upper ~100 meters in the central equatorial Pacific (Fig. 2i), which is generally above the218thermocline depth (black curve in Supplemental Fig. S4a). The influence of volcanic forcing on

219

ocean temperatures is largely confined above the thermocline depth and acts primarily through220

shortwave radiation and ENSO dynamical processes69. This suggests that the same subsurface

221

oceanic process may operate in both DPLE and DPLE_NoVolc, but cannot provide a source of222

predictability for upper layer temperature and SST when the tropical Pacific system is perturbed

223

by volcanic forcing. Related to the western Pacific pattern bias in subsurface temperatures (Fig.

2241), ACC is insignificant (or even negative) in the western Pacific for the upper 250 meters. In

225subsequent sections, we will explore the oceanic processes contributing to the predictable

226

component of TPDV.

9

227
228
229
230
231
232
233
234
235
236

Fig. 2 Decadal prediction skill is improved by ocean initialization and excluding volcanic
forcing. Anomaly correlation coefficient (ACC) of quadratically detrended 10-yr running-mean(a–e) SST and (f–j) equatorial Pacific (3°S–3°N) subsurface ocean temperature during 1955–2022
in the ensemble-mean forecasts averaged over FY1–10 for (a and f) DPLE and (b and g)
DPLE_NoVolc, and 10-yr running-mean ensemble mean of (c and h) LE, and their differences [(dand i) DPLE_NoVolc minus DPLE, and (e and j) DPLE minus LE]. The ACC skill is verified
against FOSI (see Fig. S3 for ACC skill verified against observations). Stippling indicates values
that are not significant at the 90% confidence level, based on bootstrapping across both time andensemble members (see Methods).

237

b. Pacific oceanic mechanisms contributing to TPDV predictability

238

In this section, we will explore the influence of three leading oceanic mechanisms (Rossby

239

waves, spiciness, and STCs) on TPDV in both FOSI and hindcasts. We illustrate the characteristics

240of the Rossby wave reflection mechanism by showing variability of the depth of the σ_θ= 25.5 kg

241

m-3 isopycnal (Fig. S4). The climatological depth of this isopycnal aligns closely with the depth242

where large variability of subsurface ocean temperature is observed for both equatorial (3°S–3°N;

243

Fig. S4b) and off-equatorial Pacific regions (15°S–20°N; Fig. S4f). Its vertical displacement
10

244

mostly reflects adiabatic temperature variability associated with oceanic wave propagation, which245

is most pronounced at 10°–15° latitude as shown in Fig. S3j and consistent with earlier studies78,79.

246

Additionally, this isopycnal (σ_θ= 25.5 kg m-3) lies below the equatorial Pacific thermocline and247

can remain unaffected by volcanic forcing during a 10-year forecast period (Fig. 2; Fig. S4a). To

248

isolate the low-frequency variability associated with TPDV (Fig. 3a), we use a 10-year running249

mean filter to smooth the annual mean fields that are dominated by ENSO variance (Fig. S5).

250

However, in contrast to the equatorial Pacific where ENSO variability dominates, off-equatorial

251

Pacific subsurface temperature shows more pronounced decadal variability (Fig. S4h). In addition,252

off-equatorial Pacific isopycnal depth shows robust decadal variations even in unfiltered data (Fig.

253

S5).

254

During the positive phase of TPDV, decadal SST warming in the central-to-eastern equatorial

255

Pacific (3°S–3°N; ~160°E–80°W) is concurrent with isopycnal deepening and vice versa during256

the negative phase of TPDV in FOSI (Fig. 3a). Phase changes of TPDV SST anomalies occurred

257

around 1962, 1977, and 1997 (note that the year in the y-axis in Fig. 3a denotes the start year of258

the 10-year average). These SST phase changes were preceded by isopycnal depth anomalies

259

propagating along the equator from the western to the eastern Pacific. These precursor isopycnal

260

depth anomalies in the western equatorial Pacific are further linked to off-equatorial (10°–15°)

261

isopycnal depth anomalies. The results suggest that there might be a delayed oceanic feedback262

contributing to TPDV phase changes. During the positive phase of TPDV, positive (10°–15°N)

263

and negative (15°–10°S) wind stress curl anomalies generate off-equatorial upwelling Rossby

264

waves. These waves propagate toward the western boundary of the Pacific, then travel equatorward

265

along the western boundary (130°–135°E or 155°–160°E zonally averaged) and reflect as

266

upwelling equatorial Kelvin waves (3°S–3°N), which may cause thermocline shoaling and SST

267

cooling in the eastern equatorial Pacific.

268

In FOSI, equatorial Pacific decadal SST anomalies are significantly correlated with the state269

of western equatorial Pacific isopycnal depth anomalies occurring more than 7 years earlier and270

with even 13 years earlier precursory isopycnal conditions in the off-equatorial Pacific (Fig. 3b271

and c). Similar results are found for decadal SST anomalies in the tropical Pacific (20°S–20°N;

272

not

273

HadISST/ORAS4/NOAA20CR in Fig. S7) shows similar results regarding preceding isopycnal274

depth conditions in the western equatorial Pacific influencing TPDV equatorial Pacific SST

shown).

Analysisofobservational

data(ERSSTv5/EN4/ERA5

in

Fig.

S6;

11

275

anomalies. However, in both EN4 and ORAS4, isopycnal variability south of the equator makes a276

more important contribution to western equatorial Pacific isopycnal depth variability (Figs. S6 and277

S7), in contrast to FOSI which shows a stronger contribution from the north (Fig. 3). This

278

discrepancy in the relative importance of oceanic conditions in the south vs. the north between

279

FOSI and observations is likely related to uncertainties in the wind stress data used to force the280

ocean model component of CESM114,80 (cf. Fig. 3 and Figs. S6 and 7), and differences in the ocean

281

solutions, assimilation, and statistical correction methods used in the observational datasets.

282

Considerable differences also exist in subsurface temperature and salinity fields between EN4 and283

ORAS4, surface wind stress between ERA5 and 20CR, and SST fields between ERSSTv5 and284

HadISST (Figs. S6–7). These observational uncertainties emphasize the need for improved ocean

285

observations and data assimilation methods to better understand the mechanisms of TPDV.

286
287
288

Fig. 3 Relation of oceanic Rossby wave reflection to TPDV. a. Longitude-time sections ofquadratically detrended and 10-yr running-mean isopycnal depth anomalies (m; shading) along the

12

289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309

off-equatorial (10°–15°N and 15°S–10°S), western Pacific boundary (130°–135°E and 155°–
160°E), and equatorial (3°S–3°N) waveguides during 1948–2015 (The year value in the y-axis
represents the start year of any 10-yr averaging window). In the equatorial segment, SST anomalies(°C, contours at intervals of 0.1; positive contours in red solid, negative contours in blue dashed,and zero contours in thick gray) are overlaid. In the off-equatorial segment, wind stress curl
anomalies (N m-3; contours at intervals of 0.5*10-8; positive contours in black solid, negative
contours in black dashed, zero contours in thick gray) are overlaid and smoothed with a nine-point
local smoothing. Note that the longitude axis is reversed for the off-equatorial segment to show
Rossby wave reflection at the western boundary. b. Timeseries of quadratically detrended and 10yr running-mean filtered SST anomalies (°C) in the equatorial Pacific (3°S–3°N, 180°–80°W;
thickened black curve; right y-axis), isopycnal depth anomalies (m) in the western equatorial
Pacific (3°S–3°N, 120°E–160°W; blue curve; left y-axis), in the northern off-equatorial western
Pacific (10°–15°N, 120°E–160°W; orange curve) and southern off-equatorial western Pacific(15°–10°S, 155°E–160°W; green curve). c. Lead-lag correlation of the 10-yr running-meanequatorial Pacific SST anomalies with the 10-yr running-mean equatorial Pacific SST anomalies(black), isopycnal depth anomalies in the western equatorial Pacific (blue curve), in the northern
off-equatorial western Pacific (orange curve) and southern off-equatorial western Pacific (green
curve) during a range of 20 lead to lag years. Negative lags correspond to isopycnal depth
anomalies leading the TPDV event peak; positive lags correspond to isopycnal depth anomalies
lagging the TPDV event peak. The filled circles indicate correlations that are statisticallysignificant at the 90% confidence level, based on a bootstrapping method (Methods).

310

Next, we compare TPDV SST anomalies with STC strength, estimated using the maximum of311

the Pacific zonally-averaged meridional overturning streamfunction as a function of latitude in

312

FOSI (Fig. 4a; positive/negative overturning streamfunction anomalies north/south of the Equator

313

indicate poleward anomalous near-surface transport and intensification of STC). The vertical314

structure of STC climatology and anomalies during different TPDV phases are shown in Fig. S8.

315

During the positive phase of TPDV, the STC slows down, with the near-surface flow exhibiting316

anomalous equatorward convergence, and vice versa during the negative phase of TPDV. The total

317

near-surface meridional transport convergence or divergence (near-surface meridional transport

318

values at 9°S minus 9°N; Fig. 4b) tends to significantly lag the equatorial Pacific decadal SST

319

anomalies by 2–3 years (Fig. 4b,c) and may act to strengthen TPDV SST anomalies via modulation

320

of equatorial upwelling once a TPDV phase is initiated. Changes in the STC strength are not321

independent of the westward propagation of Rossby waves, which alter the zonal slope of the322

pycnocline and create meridional geostrophic current anomalies. Consequently, the collective323

effect of Rossby wave propagation leads to a lag in zonally integrated meridional streamfunction

324

variations relative to SST variations53. Here, we use the Pacific zonally integrated transport, which325

may better capture the equilibrium STC response to surface wind stress forcing, while the zonally

13

326

interior transport (excluding the western boundary current) reflects the transient STC response and327

may shorten the lag between STC and SST decadal variability.

328

To illustrate the role of spiciness advection on TPDV, we present a latitude-time diagram of329

decadal temperature anomalies on the time-varying isopycnal depth (σ_θ= 25.5 kg m-3; Fig. 4d).

330

Spiciness anomalies originate in the subtropical eastern Pacific and propagate to the western

331

tropical Pacific (the pathway is denoted in Fig. S4l). The spiciness of the western equatorial Pacific332(5°S–5°N) is largely controlled by the advection from the South Pacific (5°S–0°; Fig. 4e) and tends

333

to be out-of-phase with the equatorial Pacific SST anomalies (Fig. 4f). These subsurface spiciness

334

anomalies thus damp TPDV SST anomalies induced by anomalous equatorial upwelling. Based335

on the diagnostic analyses shown in Figs 3 and 4, isopycnal depth variability and associated336

Rossby wave activity appear to be the most important precursors for TPDV during 1948–2015 in

337

FOSI.

338
14

339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359

Fig. 4 Relation of Pacific STC and spiciness anomalies to TPDV. a. Latitude-time sections ofquadratically detrended and 10-yr running-mean filtered vertical maximum of the Pacific zonal
average meridional overturning streamfunction (Sv; shading) and tropical Pacific (160°E–120°W)
SST (°C; contours at intervals at 0.1; positive contours in red solid, negative contours in blue
dashed, and zero contours in gray thickened) during 1948–2015 (The abscissa represents the startyear of each 10-yr averaging window). b. Timeseries of quadratically detrended and 10-yr runningmean filtered SST (°C) anomalies in the equatorial Pacific (3°S–3°N, 180°–80°W; thickened black
curve; right y-axis), STC (Sv) strength at 9°N (orange curve; sign flipped so that positive/negativevalues denote equatorward/poleward near-surface transport; left y-axis), 9°S (green curve), andthe total convergence/divergence between 9°N and 9°S (blue curve). c. Lead-lag correlation of theequatorial Pacific SST anomalies in year 0 with the equatorial Pacific SST anomalies (black), STC
strength at 9°N (orange curve), 9°S (green curve), and in total (blue curve) during a range of 20
lead to lag years. Positive (negative) lags indicate that TPDV leads (lags) the STC strength
anomalies. The filled circles indicate correlations that are statistically significant at the 90%
confidence level based on a bootstrapping method. d. As in a., but the shading indicates the
spiciness anomalies (°C; left y-axis) along the path denoted by dots in Fig. S4l. e. Timeseries ofquadratically detrended and 10-year running-mean filtered SST (°C) anomalies in the equatorial
Pacific as in panel b, and spiciness anomalies (°C) along the advection pathway from panel d,meridionally averaged in the northern equatorial Pacific (0°–5°N; orange curve), southernequatorial Pacific (5°S–0°; green curve), and the total equatorial Pacific (5°S–5°N; blue curve). f.
As in c., but spiciness anomalies in e.

360

To confirm the role of different oceanic processes in the decadal hindcasts, we correlate the361

predicted equatorial Pacific SST in FY1–10 during 1955–2016 with concurrent oceanic fields in

362

FY1–10 of DPLE_NoVolc (Fig. 5a–c) and the corresponding initial conditions in November

363(Nov0) in FOSI during 1954–2015 (Fig. 5d–f). In FY1–10, positive decadal SST anomalies in the364

tropical Pacific are associated with isopycnal deepening in the eastern equatorial Pacific and365

shoaling in the west (Fig. 5a–c). The predicted equatorial Pacific SST index in FY1–10 does not366

show significant correlations with the initial SST anomalies over most regions in the tropical

367

Pacific. In contrast, it shows significant correlations with subsurface temperature anomalies in the368western equatorial Pacific and with isopycnal depth anomalies not only over the equatorial western

369

Pacific but also in the off-equatorial Pacific region (10°–20° latitude), suggesting the critical role370

of Rossby wave initialization.

371

To show the propagation and influence of subsurface processes through the 10-year forecast372

period in DPLE_NoVolc, we correlate the predicted equatorial Pacific SST index in FY5–10 with373

several fields at the concurrent (FY5–10) and earlier 6-year forecast periods (FY4–9, 3–8, 2–7,

374

and 1–6). This analysis reveals that equatorial Pacific SST anomalies in FY5–10 (Fig. 5g5) can be375

traced back to subsurface temperature anomalies in the equatorial Pacific (Fig. 5h1), along with

15

376

off-equatorial and equatorial isopycnal depth anomalies in FY1–6 (Fig. 5i1). As significant SST

377

correlations intensify in FY2–7 (Fig. 5g2), surface westerly anomalies start to develop over the378

western equatorial Pacific in FY3–8 (Fig. 5i2), deepening the isopycnal in the eastern equatorial

379

Pacific (Fig. 5i2). The coupling of SST, wind, and isopycnal depth can enhance the predictability

380

provided by the initial subsurface temperature anomalies. The equatorial Pacific SST index in381

FY5–10 also shows significant correlations with the Pacific STC overturning streamfunction,

382

originating in the South Pacific (20°Sº–10°S) in FY1–6 (Fig. 5j1). When a positive TPDV starts

383

to develop at FY1–6, the positive overturning streamfunction anomalies in the subsurface South

384

Pacific weaken the south Pacific branch of anticlockwise STC (with negative values in dashed385

contours), amplifying TPDV SST anomalies in the later lead times. In DPLE, however, the tropical

386

Pacific SST anomalies at FY5–10 show much weaker correlations with oceanic conditions at387

earlier lead times than in DPLE_NoVolc, suggesting that volcanic forcing perturbs the linkage

388

between the initial condition memory and SST variability at late lead times (Fig. S9).

16

389
390
391
392
393
394
395
396
397
398
399
400
401
402

Fig. 5 Source and persistence of predictability in initialized decadal forecasts. (a–f)
Correlation maps with detrended SST anomalies in the central-eastern equatorial Pacific (3°S–3°N,
180°E–120°W; denoted by the white box in a) averaged in FY1–10 during 1955–2016 in
DPLE_NoVolc. Correlations (color shading) are calculated with the quadratically detrended (a andd) global SST, (b and e) ocean temperature in the equatorial Pacific (3°S–3°N), and (c and f)
tropical Pacific isopycnal depth (σ_θ= 25.5 kg m-3) anomalies at (a–c) FY1–10 during 1955–2016,
and (d–f) the corresponding FOSI initial conditions in Nov0 during 1954–2015. (g–j) Correlation
maps with detrended SST anomalies in the central-eastern equatorial Pacific (3°S–3°N, 180°E–
120°W; denoted by the white box in g5) averaged in FY5–10 during 1955–2016 in DPLE_NoVolc.
Correlations (color shading) are calculated with the detrended (g1–5) global SST, (h1–5) ocean
temperature in the equatorial Pacific (3°S–3°N), (i1–5) tropical Pacific isopycnal depth (σ_θ= 25.5
kg m-3), surface winds, and (j1–5) STC overturning streamfunction at FY1–6, FY2–7, FY3–8,
FY4–9, and FY5–10 during 1955–2016. STC climatology during 1964–2015 in Fig. S8 are
17

403
404
405
406

overlaid in j1–5 using black contours at intervals of 8 Sv [positive (negative) contours in solid(dashed) contours and zero contours thickened]. The stippling indicates shaded values that are notsignificant at the 90% confidence level, based on a two-tailed Student’s t-test. Only significant
correlations for surface wind vectors are shown.

407

Insignificant correlations between predicted TPDV SST anomalies at FY1–10 and initial408

tropical Pacific SST anomalies in Nov0 (Fig. 5) suggest that the initial ENSO state does not affect409

the overall TPDV predictability at 10-year lead. The minor role of ENSO is further examined using

410

a Singular Value Decomposition (SVD) analysis81 between predicted tropical Pacific SST and the411

global SST initial conditions in Nov0 (Fig. S10). The leading SVD mode (SVD1) explains 59%

412

of the total squared covariance between predicted tropical Pacific SST and the global SST initial413

conditions. The timeseries of the two expansion coefficients of SVD1 show decadal variations,414

and the associated heterogeneous correlation maps onto the SVD timeseries closely mirror the415

linear correlation analysis in Fig. 5. In contrast, the two expansion coefficients of SVD2 show

416

strong interannual variability associated with ENSO, accounting for 15% of the covariance. This

417

result suggests that there is only a weak interannual component in the predicted tropical Pacific418

SST anomalies in FY1–10, manifesting as an equatorial mode affected by the initial ENSO states419

in Nov0.

420

c. Role of Interbasin interactions in affecting TPDV prediction during 1999–2008

421

In Fig. 5d, we also identify significant correlations between the predicted TPDV index and the422

initial SST conditions in several remote ocean areas, including negative correlations over the North

423

Atlantic and positive correlations over the extratropical oceans of the Southern Hemisphere. To

424

test if these correlations indicate causality, and to separate the effects of initial conditions in

425

different ocean basins on the tropical Pacific predictions, we conduct a set of regional initialization

426

experiments for the period 1999–2008 (Fig. 6b–e; Methods). During this period, the tropical

427

Pacific shows negative decadal SST anomalies in observations (Fig. S11a) and FOSI (Fig. 6f),

428

which are well predicted by the DPLE_NoVolc forecasts initialized on Nov 1, 1998 in Fig. S11b

429

& Fig. 6g with different drift correction methods (See Methods). The comparisons among the430

sensitivity forecasts with regional ocean initialization suggest that this skillful prediction of431

tropical Pacific decadal cooling is mainly associated with tropical Pacific initialization (Fig. 6h),432

which shows negative isopycnal depth anomalies in the tropical Pacific in the initial conditions

433(Fig. S11c). The North Atlantic (20°N–60°N) initialization experiment produces positive SST
18

434

anomalies over both the North and northern tropical Atlantic (0°–20°N; Fig. 6d and i). In contrast

435

to previous SST-restoring experimentse.g.,59, the tropical Atlantic induced warming does not436

generate tropical Pacific cooling. This is because the net surface heat flux over the tropical Atlantic437

is downward (Fig. 6n), and so SST warming there is not an active forcing that can affect438

atmospheric teleconnections, consistent with previous studies61,62. The tropical Atlantic shows439

insignificant SST anomalies in the Full initialization experiment (Fig. 6g), possibly due to440

competing effects from the tropical Pacific (Fig. 6h) and North Atlantic (Fig. 6i). The Southern

441

Hemisphere ocean initialization (Southern Ocean, Indian Ocean, and Atlantic Ocean sectors) is442

responsible for the predicted warming over the Southern Ocean in the Full Initialization

443

experiment (cf. Figs. 6g and j). However, the sign of the predicted ensemble-mean Southern Ocean

444

warming is opposite to that observed and seen in FOSI (Fig. S11a and Fig. 6f) and might contribute445

to the tropical Pacific warming in the SH initialization experiment (Fig. 6j). This suggests that the446

predicted tropical Pacific cooling in the Full initialization experiment might be underestimated due

447

to prediction errors in the Southern Ocean.

448

Given that interbasin interactions may vary from case to case, additional regional initialization

449

experiments with more start dates are required to robustly isolate the influence of different ocean

450

basins on TPDV prediction skill in general. A more refined analysis of the role of different ocean

451

regions is also necessary. For example, further experiments are needed to distinguish the impacts452

of isopycnal initial conditions in the equatorial vs. off-equatorial Pacific on TPDV prediction skill.

453

Although the correlation analysis does not suggest an influence from initial North Pacific SST

454

anomalies on TPDV prediction skill (Fig. 5d), this might be due to inherent model underestimation455

of the coupling between the extratropics and tropics82. Similarly, additional experiments are456

required to isolate the role of the Indian Ocean. The significant correlations between the predicted

457

TPDV index and the initial Indian SST conditions are confined to the extratropical South Indian

458

Ocean, rather than the tropical region, which has been suggested as an important factor for459

TPDV64,65. Notably, most areas of the Indian Ocean exhibit significant prediction skill in detrended

460

SST in both decadal forecasts and LE, suggesting that the Indian Ocean decadal prediction skill is461

largely attributable to external forcing (Figs. S13 and S14). In the selected case study of 1999–

462

2008, the cooling in the Indian Ocean is not well predicted by any of the experiments (Fig. 6g–j),463

potentially diminishing its influence on the tropical Pacific.

19

464
465
466
467
468
469
470
471
472
473
474
475

Fig. 6 The influence of different ocean basins on TPDV predictability during 1999–2008. SST(°C) anomalies in November 1998 in (a) ERSSTv5 and (b) FOSI. (c–e) The areas with SST values
on November 1 of 1998 denote the ocean regions that have the initialization of full-depth ocean
temperature and salinity anomalies added on the climatology from FOSI; the November 1 FOSI
climatology during 1954–2015 is used everywhere else. Quadratically detrended SST (°C)
anomalies during 1999–2008 in (f) FOSI and 10-member ensemble mean forecasts initialized on
November 1st, 1998, including (g) Full initialization experiment, (h) Tropical Pacific initialization,(i) North Atlantic initialization, and (j) Southern Ocean initialization. (k–o) as in (f–j) but for thesurface net heat flux (SHF; W m-2; positive values heat the ocean). Stippling indicates values thatare not significant at the 90% confidence level, based on bootstrapping across ensemble members(see Methods).

476

Summary and Discussion

477

Our study investigates a set of oceanic mechanisms that may provide a source of predictability

478

for TPDV based on observations and a hierarchy of model simulations and hindcasts. The CESM1

20

479

decadal prediction system without volcanic forcing (DPLE_NoVolc) shows high skill in predicting

480

observed TPDV during the 1950s to the present. The high predictability of TPDV arises from the481

subsurface ocean initial conditions, particularly decadal isopycnal depth variability associated with

482

oceanic Rossby wave adjustments in the tropical Pacific. The predictability inherent in the initial

483

isopycnal conditions is further reinforced by the response of the subtropical cells and ocean-484

atmosphere coupling processes in the tropical Pacific throughout the 10-year forecasts. Although485

we can rely on the initial isopycnal conditions for useful TPDV prediction skill, the origins of486

decadal isopycnal depth variability and the associated surface wind stress fluctuations driving

487

isopycnal variability remain less clear, as also suggested by Capotondi et al.18. We also find that488

the initial ENSO SST state in a particular year plays a very minor role in affecting the overall

489

TPDV predictability, but ENSO decadal modulation may contribute to establishing decadal wind490

and isopycnal variability associated with TPDV.

491

To investigate the potential role of inter-basin interactions in TPDV prediction skill, we492

conducted regional initialization experiments which suggest that the predictability of TPDV SST

493

cooling during 1999–2008 arises mainly from initial conditions in the tropical Pacific as opposed494

to those in the North Atlantic or tropical Atlantic. This result adds to previous work highlighting

495

potential issues in using SST-restoring experiments to study the causality of interbasin ocean

496interactions61. During 1999–2008, CESM1 decadal forecasts erroneously predict positive decadal

497

SST anomalies in the Southern Ocean, which are opposite to those observed and may contribute498

to an erroneous tropical Pacific SST warming. The Southern Ocean errors could be reduced by499

increasing the model resolution used for decadal hindcasts, which enhances prediction skill in the500

tropical Pacific75. On the other hand, some studies suggest that models may underestimate

501

Atlantic-Pacific Ocean interactions due to climatological biases, which displace atmospheric502

convective regimes and sensitivities22,83. The influence of interbasin ocean interactions on decadal503

predictability needs to be further analyzed using regional experiments with more initialization504

dates and other properly designed experiments.

505

Our results are subject to biases in CESM1 at 1° resolution and to considerable observational

506

uncertainties. To investigate the dependence of the results on the choice of model and507

observational datasets, it will be important to examine other CMIP6 decadal forecasts84 as well as

508

Ocean Model Intercomparison Project simulations85. Enhancing observational and reanalysis

509

datasets is crucial for studying decadal ocean variability. Addressing the low CMIP5/6 decadal
21

510

prediction skill in the Pacific Ocean is key for improving the skill of current decadal climate

511

prediction systems, which bridge the gap between seasonal forecasts and centennial climate

512

projections and provide useful information on climate adaptation and resilience for decision513

makers in many sectors of the economy.

514

Methods515

a. CESM1 simulations and forecasts

516

We analyze the dynamics and predictability of TPDV using the CESM1, a global Earth system

517

model consisting of atmosphere, ocean, land, and ice components linked by a flux coupler86. All

518

experiments are conducted using the same model version as in the 40-member CESM1 Large

519

Ensemble (CESM1 LE87) at nominal 1° latitude-longitude resolution, including the Community

520

Atmosphere Model, version 5 (CAM588) with 30 vertical levels; the Parallel Ocean Program,521

version 2 (POP289) with 60 vertical levels; the Community Land Model, version 4 (CLM490); and522

the Los Alamos National Laboratory Community Ice Code, version 4 (CICE491). To compensate

523

for the lack of comprehensive and consistent observations and reanalysis of subsurface oceanic

524

processes that could influence TPDV, we analyze a forced ocean-sea ice simulation (FOSI), in

525

which the ocean and sea ice components are forced with observed atmospheric and surface flux

526

fields. The surface fluxes are derived using bulk formulae based on observed atmospheric fields

527

from the Coordinated Ocean-Ice Reference Experiment forcing dataset92. FOSI provides a realistic528

simulation of SST variability during 1948–2015, despite some inconsistencies in the Southern

529

Ocean93 (Fig. S12). Estimates for subsurface ocean temperatures and salinities are subject to larger

530

observational uncertainties than for SST, as discussed in the following result sections.

531

FOSI provides the ocean and sea ice conditions needed to initialize the Decadal Prediction

532

Large Ensemble (DPLE14) and a parallel decadal forecast set that excludes historical volcanic

533

aerosol forcing (DPLE_NoVolc69). In DPLE, 40-member forecasts are initialized from identical534

oceanic and sea ice conditions from FOSI on November 1st of each year during 1954–2015 and535

integrated for 122 months. The atmosphere and land initial conditions are expected to play a very536

minor role in contributing to decadal-scale climate predictability, and are obtained from a random537

member of the CESM1 LE without any observational constraints except for the historical radiative

538

forcing. The ensemble spread among individual members is created by adding round-off level

22

539

perturbations to the initial atmospheric temperatures. The CMIP5 historical forcings for 1954–

540

2005, and representative concentration pathway (RCP) 8.5 forcings for 2006 and onwards, are541

used as the external forcings for the DPLE forecasts. DPLE_NoVolc follows the DPLE protocol,

542

except that it excludes historical volcanic aerosol forcing during 1954–2005 and has a smaller543

ensemble size of 10. The comparison between DPLE and DPLE_NoVolc isolates the effect of544

historical volcanic forcing on decadal prediction skill and predictability69.

545

To isolate the role of initialization in affecting predictability and prediction skill, we compare546

the DPLE with the uninitialized CESM1 LE14. The CESM1 LE is comprised of 40-member

547

historical simulations subject to the CMIP5 forcings during 1920–2100, including the historical

548

volcanic forcing87. We quadratically detrend all data sets (observed and simulated) to remove the549

forced climate change signal during 1954–2015. Results are very similar if we use other methods

550

of estimating and removing the forced climate signal, including subtracting a linear trend or the551

ensemble mean of the CESM1 LE69. We also evaluate the performance of CESM1 in simulating552

intrinsic TPDV variability in an 1801-yr CESM1 control simulation under preindustrial

553

atmospheric greenhouse gas concentrations.

554

b. Observational datasets

555

We assess the simulation realism and retrospective forecast skill of TPDV using multiple

556

observational and reanalysis datasets available from the 1950s to 2022. The SSTs are taken from557

the National Oceanic and Atmospheric Administration (NOAA) Extended Reconstruction Sea

558

Surface Temperature version 594 (ERSSTv5) dataset at 2° spatial resolution, and the Hadley Centre

559

Sea Ice and SST dataset95 (HadISST) at 1° spatial resolution, during 1954-2022. The ocean

560

temperature and salinity are taken from the Met Office Hadley Center EN496 (EN4) with bias

561

correction97 at 1° horizontal resolution with 42 vertical levels during 1954-2022, and the European

562

Centre for Medium-Range Weather Forecasts (ECMWF) Ocean Reanalysis System 498 (ORAS4)563

at 1° horizontal resolution with 42 vertical levels during 1958–2017. The ocean potential density564

is calculated using the temperature and salinity based on the equation of state for seawater99.

565

Surface wind stress is taken from the fifth generation ECMWF atmospheric reanalysis of global

566

climate100 (ERA5) for 1954-2022. Surface wind stress is also calculated based on surface wind567

components from the NOAA Twentieth Century Reanalysis version 3101 (NOAA20CRv3) during

23

568

1954–2022 for comparisons. All observational data are re-gridded to the model grid at ~1° using

569

a bilinear interpolation before comparison.

570

c. Analysis methods

571

Decadal forecasts are initialized from full-field oceanic states constrained by observations and572

will drift toward the model’s biased climatology as the forecasts progress. To obtain forecast573

anomalies, forecast lead-dependent climatologies are removed from each ensemble member. For

574

each forecast year (FY) 1–10, a climatology is determined by averaging the ensemble mean annual

575

mean forecasts at that lead across 1964–2015. The choice of the 1964–2015 period as the baseline576

climatology ensures a consistent sample size (52 years) for each lead time. Anthropogenic climate

577

change is estimated as a quadratic fit of the ensemble mean forecast anomalies across 1954–2015,578

as a function of lead time (i.e., FY1–10). This estimated climate change signal is then subtracted

579

from individual drift-corrected ensemble members, to obtain detrended forecast anomalies.

580

We evaluate the model performance in capturing observed TPDV using EOF analysis, and by581

calculating the anomaly correlation coefficient (ACC) between the ensemble-mean forecasts and582

observations. Statistical significance of the results is tested using a bootstrapping method14,102. At

583

each spatial location, we determine a nonparametric bootstrap distribution of the forecast ACC,584

ACC difference, correlation, or composite anomalies by resampling (with replacement) the585

forecast ensembles across both the time and/or ensemble member dimensions, and then586

recalculating the statistics for each of the 5,000 samples. The calculation of 5,000 values is587

performed using 10-member ensembles for all model analyses, considering the varied ensemble

588

sizes of DPLE (40 ensemble members), DPLE_NoVolc (10), and LE (40). A positive value is589

deemed significant at the 90% confidence level if its bootstrapped distribution contains fewer than590

500 values below zero (p<500/5,000=0.1), and vice versa for a negative value.

591

Several metrics are used to examine oceanic processes relevant to TPDV. Mixed layer depth592

in CESM1-FOSI is a monthly model output defined using a maximum buoyancy gradient593

criterion103. Monthly thermocline depth is calculated as the depth of the maximum vertical594

temperature gradient. We use the monthly isopycnal depth where the potential density is equal to595

25.5 kg m-3 (σ_θ=25.5 kg m-3) to capture the evolution of both off-equatorial and equatorial oceanic

596

wave characteristics. Spiciness is defined using the monthly temperature along the time-dependent597

monthly isopycnal surface where σ_θ=25.5 kg m-3.
24

598

d. Regional initialization experiments599

To examine the impact of initial ocean conditions in different ocean basins on TPDV

600

predictions, we conducted a set of sensitivity experiments with regional initialization. We focus

601

on the period of 1999–2008, during which the tropical Pacific exhibited negative decadal SST

602

anomalies in observations. This tropical Pacific decadal cooling is well predicted by the603

DPLE_NoVolc forecasts initialized from the global oceanic states estimated by FOSI on

604

November 1, 1998 with a ~10-year lead. We refer to this DPLE_NoVolc forecast as the Global

605

Initialization experiment (Fig. 6). To separate which part of the ocean initialization is more

606

important to this prediction case, we conducted the following sensitivity experiments (Fig. 6): 1)607

A control ensemble initialized with the global FOSI climatology during 1954–2015 (Climatology

608

initialization); 2) An ensemble initialized with the global FOSI climatology and full-depth ocean609

temperature and salinity anomalies in the North Atlantic (North Atlantic initialization; 20°–60°N);

6103) Similar to 2), but focusing on Tropical Pacific initialization (20°S–20°N); and 4) Southern

611

Hemisphere Ocean initialization, including the Southern Ocean, tropical Indian Ocean and tropical

612

South Atlantic Ocean. To avoid discontinuities in oceanic forcing at the boundaries, we apply a613

linear interpolation over 5° latitude bands. Each experiment includes 10-member forecasts subject

614

to external radiative forcings from November 1998 to December 2008, as was done for the Global

615

Initialization experiment. Forecast anomalies are computed by subtracting the control simulation

616(Experiment 1) from Experiments 2–4 and the Global Initialization experiment. This approach to617

calculating forecast anomalies assumes a consistent time-dependent model drift across all

618

experiments and is used to prevent the extensive computing resources required to reproduce619

initialized forecasts for every year from 1954 to 2015 for each experiment. The forecast anomalies

620

for the Global Initialization experiment, obtained by removing the control simulation, are very

621

similar to the anomalies calculated using the traditional drift correction (Fig. S11b) and quadratic

622

detrending (Fig. 6g).

623
624
625
626
627
628
629
630
631

Data availability: NCAR’s Climate Data Gateway provides the output from CESM1 DPLE
(https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.CESM1-CAM5-DP.html), FOSI
(https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.DPLEFOSI.ocn.proc.monthly_ave.html), and LE(https://www.earthsystemgrid.org/dataset/ucar.cgd.ccsm4.cesmLE.html). DPLE_NoVolc isavailable through https://portal.nersc.gov/archive/home/c/ccsm/www/CESM1-CAM5-DP-NoV.
Observational and reanalysis datasets used in this study are available online: ERSSTv5 SSTs
from https://psl.noaa.gov/data/gridded/data.noaa.ersst.v5.html, HadISST from
25

632
633
634
635
636
637
638
639
640
641

http://www.metoffice.gov.uk/hadobs/, ORAS4 from ftp://ftpicdc.cen.unihamburg.de/EASYInit/ORA-S4, EN4 fromhttps://www.metoffice.gov.uk/hadobs/en4/download.html, ERA5 fromhttps://www.ecmwf.int/en/forecasts/dataset/ecmwf-reanalysis-v5, and NOAA20CRv3 fromhttps://www.psl.noaa.gov/data/gridded/data.20thC_ReanV3.html.

642

Acknowledgments: We would like to thank Elizabeth Maroon for helpful discussion on the design643

of the regional initialization experiments, and Nan Rosenbloom for providing the example scripts

644

to run the CESM1 decadal forecasts. We thank Ping Chang and Who Kim for their insightful

645

comments and suggestions on our analysis. We also thank Liping Zhang for her helpful comments646

on the manuscript. The CESM project is supported primarily by the National Science Foundation

647(NSF). The National Center for Atmospheric Research (NCAR) is a major facility sponsored by648

the NSF under Cooperative Agreement 1852977. X.W. was supported by an Advanced Study

649

Program postdoctoral fellowship from NCAR. S.Y acknowledges support from award

650

NA20OAR4310408 of the Climate Variability and Predictability program of NOAA’s Climate

651

Program Office. A.C. was supported by the NOAA Climate Program Office Climate Variability

652

and Predictability Program and by DOE Award # DE-SC0023228. PMEL contribution no.

653

XXXX.

654
655
656
657
658
659

Code availability: Codes generated during this study are available from the corresponding author
upon reasonable request.

Author Contributions: X.W. conceptualized the study, conducted the analysis with detailed
discussions with S.G.Y. and C.D., and wrote the initial draft. All authors contributed to interpreting
the results and editing the manuscript.
Competing interests: The authors declare no competing interests.

660
661

REFERENCES

662
663
664

1. Lorenzo, E. D. et al. Synthesis of Pacific Ocean Climate and Ecosystem Dynamics.
Oceanography 26, 68–81 (2013).

665
666

2. Lorenzo, E. D. et al. Modes and Mechanisms of Pacific Decadal-Scale Variability. Annu Rev
Mar Sci 15, (2022).

26

667
668

3. Kosaka, Y. & Xie, S.-P. Recent global-warming hiatus tied to equatorial Pacific surface
cooling. Nature 501, 403–407 (2013).

669
670
671

4. Meehl, G. A., Hu, A., Arblaster, J. M., Fasullo, J. & Trenberth, K. E. Externally Forced andInternally Generated Decadal Climate Variability Associated with the Interdecadal PacificOscillation. J. Clim. 26, 130408111212000 (2013).

672
673

5. England, M. H. et al. Recent intensification of wind-driven circulation in the Pacific and theongoing warming hiatus. Nature Climate Change 4, 222–227 (2014).

674
675

6. Rodgers, K. B., Friederichs, P. & Latif, M. Tropical Pacific Decadal Variability and Its
Relation to Decadal Modulations of ENSO. J Climate 17, 3761–3774 (2004).

676
677

7. Vimont, D. J. The Contribution of the Interannual ENSO Cycle to the Spatial Pattern ofDecadal ENSO-Like Variability. J Climate 18, 2080–2092 (2005).

678
679

8. Okumura, Y. M., Sun, T. & Wu, X. Asymmetric Modulation of El Niño and La Niña and the
Linkage to Tropical Pacific Decadal Variability. J Climate 30, 4705–4733 (2017).

680
681

9. Sun, T. & Okumura, Y. M. Impact of ENSO‐Like Tropical Pacific Decadal Variability on the
Relative Frequency of El Niño and La Niña Events. Geophys Res Lett 47, (2020).

682
683

10. Power, S. et al. Decadal climate variability in the tropical Pacific: Characteristics, causes,
predictability, and prospects. Science 374, eaay9165 (2021).

684
685

11. Doblas-Reyes, F. et al. Initialized near-term regional climate change prediction. Nat
Commun 4, 1715 (2013).

686
687

12. Kirtman, B., Power, S., Adedoyin, A., Boer, G. & Bojariu, R. Near-term climate change:
projections and predictability. (2013).

688
689

13. Smith, D. et al. Robust skill of decadal climate predictions. Npj Clim Atmospheric Sci 2, 13(2019).

690
691
692

14. Yeager, S. et al. Predicting near-term changes in the Earth System: A large ensemble ofinitialized decadal prediction simulations using the Community Earth System Model. B Am
Meteorol Soc 99, 1867–1886 (2018).

693
694

15. Meehl, G. A. et al. Decadal Climate Prediction: An Update from the Trenches. B Am
Meteorol Soc 95, 243–267 (2014).

695
696

16. Newman, M. et al. The Pacific Decadal Oscillation, Revisited. J Climate 29, 4399–4427(2016).

697
698

17. Liu, Z. & Lorenzo, E. D. Mechanisms and Predictability of Pacific Decadal Variability. Curr
Clim Change Reports 4, 128–144 (2018).

27

699
700

18. Capotondi, A. et al. Mechanisms of tropical Pacific decadal variability. Nat. Rev. EarthEnviron. 1–16 (2023) doi:10.1038/s43017-023-00486-x.

701
702
703

19. Guilyardi, E., Capotondi, A., Lengaigne, M., Thual, S. & Wittenberg, A. T. ENSO Modeling:
History, Progress, and Challenges. Geophys. Monogr. Ser. 199–226 (2020)
doi:10.1002/9781119548164.ch9.

704
705

20. Lee, S. et al. On the future zonal contrasts of equatorial Pacific climate: Perspectives from
Observations, Simulations, and Theories. Npj Clim Atmospheric Sci 5, 82 (2022).

706
707
708

21. Wills, R. C. J., Dong, Y., Proistosecu, C., Armour, K. C. & Battisti, D. S. Systematic Climate
Model Biases in the Large‐Scale Patterns of Recent Sea‐Surface Temperature and Sea‐Level
Pressure Change. Geophys Res Lett 49, (2022).

709
710

22. Li, C., Dommenget, D. & McGregor, S. Trans-basin Atlantic-Pacific connections further
weakened by common model Pacific mean SST biases. Nat. Commun. 11, 5677 (2020).

711
712

23. McGregor, S. et al. Recent Walker circulation strengthening and Pacific cooling amplified by
Atlantic warming. Nat Clim Change 4, 888–892 (2014).

713
714

24. Lin, R., Zheng, F. & Dong, X. ENSO Frequency Asymmetry and the Pacific Decadal
Oscillation in Observations and 19 CMIP5 Models. Adv. Atmos. Sci. 35, 495–506 (2018).

715
716
717

25. Ogata, T., Xie, S.-P., Wittenberg, A. & Sun, D.-Z. Interdecadal Amplitude Modulation of El
Niño–Southern Oscillation and Its Impact on Tropical Pacific Decadal Variability*. J. Clim. 26,
130412124111006 (2013).

718
719

26. Yeh, S.-W. & Kirtman, B. P. ENSO Amplitude Changes due to Climate Change Projections
in Different Coupled Models. J Climate 20, 203–217 (2007).

720
721

27. Chung, C. T. Y., Power, S. B., Sullivan, A. & Delage, F. The role of the South Pacific inmodulating Tropical Pacific variability. Sci. Rep. 9, 18311 (2019).

722
723
724

28. Liguori, G. & Lorenzo, E. D. Separating the North and South Pacific Meridional Modes
Contributions to ENSO and Tropical Decadal Variability. Geophys. Res. Lett. 46, 906–915(2019).

725
726

29. Sun, T. & Okumura, Y. M. Role of Stochastic Atmospheric Forcing from the South andNorth Pacific in Tropical Pacific Decadal Variability. J Climate 32, 4013–4038 (2019).

727
728

30. Zhao, Y. & Lorenzo, E. D. The impacts of Extra-tropical ENSO Precursors on Tropical
Pacific Decadal-scale Variability. Sci Rep-uk 10, 3031 (2020).

729

31. Hasselmann, K. Stochastic climate models Part I. Theory. Tellus 28, 473–485 (1976).

28

730
731

32. Frankignoul, C., Müller, P. & Zorita, E. A Simple Model of the Decadal Response of the
Ocean to Stochastic Wind Forcing*. J Phys Oceanogr 27, 1533–1546 (1997).

732
733

33. Schneider, N., Miller, A. J. & Pierce, D. W. Anatomy of North Pacific Decadal Variability. J
Climate 15, 586–605 (2002).

734
735

34. Kwon, Y.-O. & Deser, C. North Pacific Decadal Variability in the Community Climate
System Model Version 2. J. Clim. 20, 2416–2433 (2007).

736

35. Stuecker, M. F. Revisiting the Pacific Meridional Mode. Sci Rep-uk 8, 3216 (2018).

737
738

36. Kim, H., Kang, S. M., Kay, J. E. & Xie, S.-P. Subtropical clouds key to Southern Ocean
teleconnections to the tropical Pacific. Proc. Natl. Acad. Sci. 119, e2200514119 (2022).

739
740

37. Wittenberg, A. T., Rosati, A., Delworth, T. L., Vecchi, G. A. & Zeng, F. ENSO Modulation:
Is It Decadally Predictable? J Climate 27, 2667–2681 (2014).

741
742

38. Latif, M. & Barnett, T. P. Causes of Decadal Climate Variability over the North Pacific andNorth America. Science 266, 634–637 (1994).

743
744

39. Knutson, T. R. & Manabe, S. Time-Mean Response over the Tropical Pacific to Increased
C02 in a Coupled Ocean-Atmosphere Model. J. Clim. 8, 2181–2199 (1995).

745
746

40. Yukimoto, S. et al. ENSO‐like interdecadal variability in the Pacific Ocean as simulated in acoupled general circulation model. J. Geophys. Res.: Oceans 105, 13945–13963 (2000).

747
748

41. Jin, F.-F. Low-Frequency Modes of Tropical Ocean Dynamics*. J Climate 14, 3874–3881(2001).

749
750
751

42. Meehl, G. A. & Hu, A. Megadroughts in the Indian Monsoon Region and Southwest North
America and a Mechanism for Associated Multidecadal Pacific Sea Surface Temperature
Anomalies. J Climate 19, 1605–1623 (2006).

752
753

43. Power, S. & Colman, R. Multi-year predictability in a coupled general circulation model.
Clim Dyn 26, 247–272 (2006).

754
755
756

44. San, S.-C., Tseng, Y.-H., Ding, R. & Lorenzo, E. D. A key role of off-equatorial subsurface
temperature anomalies in Tropical Pacific Decadal Variability. npj Clim. Atmos. Sci. 7, 109(2024).

757
758

45. Gu, D. & Philander, S. G. H. Interdecadal Climate Fluctuations That Depend on Exchanges
Between the Tropics and Extratropics. Science 275, 805–807 (1997).

759
760

46. Schneider, N. et al. Pacific thermocline bridge revisited. Geophys. Res. Lett. 26, 1329–1332(1999).

29

761
762

47. Munk, W. Internal waves and small-scale processes. Evolution of Physical Oceanography—
Scientific Surveys in Honor of Henry Stommel 264–291 (1981).

763
764

48. Schneider, N. A decadal spiciness mode in the tropics. Geophys. Res. Lett. 27, 257–260(2000).

765
766

49. Yeager, S. G. & Large, W. G. Late-Winter Generation of Spiciness on Subducted Isopycnals.
J Phys Oceanogr 34, 1528–1547 (2004).

767
768
769

50. Zeller, M., McGregor, S., Sebille, E. van, Capotondi, A. & Spence, P. Subtropical-tropicalpathways of spiciness anomalies and their impact on equatorial Pacific temperature. Clim Dynam56, 1131–1144 (2021).

770
771

51. San, S.-C. & Tseng, Y. Aleutian low/PDO forces a decadal subsurface spiciness propagating
mode in the North Pacific. Clim. Dyn. 1–19 (2023) doi:10.1007/s00382-023-06938-w.

772
773

52. McPhaden, M. J. & Zhang, D. Slowdown of the meridional overturning circulation in theupper Pacific Ocean. Nature 415, 603–608 (2002).

774
775

53. Capotondi, A. & Qiu, B. Decadal Variability of the Pacific Shallow Overturning Circulationand the Role of Local Wind Forcing. J. Clim. 36, 1001–1015 (2023).

776

54. Cai, W. et al. Pantropical climate interactions. Science 363, eaav4236 (2019).

777
778

55. Kucharski, F., Kang, -S I, Farneti, R. & Feudale, L. Tropical Pacific response to 20th century
Atlantic warming. Geophys Res Lett 38, n/a-n/a (2011).

779
780

56. Kucharski, F. et al. Atlantic forcing of Pacific decadal variability. Clim Dynam 46, 2337–
2351 (2015).

781
782

57. Li, X., Xie, S.-P., Gille, S. T. & Yoo, C. Atlantic-induced pan-tropical climate change overthe past three decades. Nat Clim Change 6, 275–279 (2015).

783
784

58. Meehl, G. A. et al. Atlantic and Pacific tropics connected by mutually interactive decadaltimescale processes. Nat Geosci 1–7 (2020) doi:10.1038/s41561-020-00669-x.

785
786

59. Ruprich-Robert, Y. et al. Impacts of Atlantic multidecadal variability on the tropical Pacific:
a multi-model study. Npj Clim Atmospheric Sci 4, 33 (2021).

787
788
789
790

60. Kim, W. M., Yeager, S. & Danabasoglu, G. Atlantic Multidecadal Variability and AssociatedClimate Impacts Initiated by Ocean Thermohaline Dynamics Atlantic Multidecadal Variability
and Associated Climate Impacts Initiated by Ocean Thermohaline Dynamics. J Climate 33,
1317–1334 (2019).

791
792

61. O’Reilly, C. H. et al. Challenges with interpreting the impact of Atlantic Multidecadal
Variability using SST-restoring experiments. npj Clim. Atmos. Sci. 6, 14 (2023).

30

793
794
795

62. Deser, C. & Phillips, A. S. Spurious Indo‐Pacific Connections to Internal Atlantic
Multidecadal Variability Introduced by the Global Temperature Residual Method. Geophys. Res.
Lett. 50, (2023).

796
797

63. Zhao, Y. & Capotondi, A. The role of the tropical Atlantic in tropical Pacific climate
variability. npj Clim. Atmos. Sci. 7, 140 (2024).

798
799

64. Luo, J.-J., Sasaki, W. & Masumoto, Y. Indian Ocean warming modulates Pacific climatechange. Proc. Natl. Acad. Sci. 109, 18701–18706 (2012).

800
801

65. Dong, L. & McPhaden, M. J. Why Has the Relationship between Indian and Pacific Ocean
Decadal Variability Changed in Recent Decades? J. Clim. 30, 1971–1983 (2017).

802
803

66. Meehl, G. A., Hu, A. & Tebaldi, C. Decadal Prediction in the Pacific Region. J Climate 23,2959–2973 (2010).

804
805

67. Timmreck, C., Pohlmann, H., Illing, S. & Kadow, C. The impact of stratospheric volcanic
aerosol on decadal‐scale climate predictions. Geophys. Res. Lett. 43, 834–842 (2016).

806
807
808

68. Ménégoz, M., Bilbao, R., Bellprat, O., Guemas, V. & Doblas-Reyes, F. J. Forecasting theclimate response to volcanic eruptions: prediction skill related to stratospheric aerosol forcing.
Environ. Res. Lett. 13, 064022 (2018).

809
810

69. Wu, X., Yeager, S. G., Deser, C., Rosenbloom, N. & Meehl, G. A. Volcanic forcing degrades
multiyear-to-decadal prediction skill in the tropical Pacific. Sci Adv 9, eadd9364 (2023).

811
812

70. Bilbao, R. et al. Impact of volcanic eruptions on CMIP6 decadal predictions: A multi-model
analysis. Earth Syst. Dyn. Discuss. 2023, 1–34 (2023).

813
814

71. Meehl, G. A. & Teng, H. Case studies for initialized decadal hindcasts and predictions for the
Pacific region. Geophys Res Lett 39, n/a-n/a (2012).

815
816

72. Meehl, G. A., Hu, A. & Teng, H. Initialized decadal prediction for transition to positive
phase of the Interdecadal Pacific Oscillation. Nat Commun 7, 11718 (2016).

817
818

73. Pohlmann, H., Kröger, J., Greatbatch, R. J. & Müller, W. A. Initialization shock in decadal
hindcasts due to errors in wind stress over the tropical Pacific. Clim. Dyn. 49, 2685–2693 (2017).

819
820

74. Teng, H., Meehl, G. A., Branstator, G., Yeager, S. & Karspeck, A. Initialization Shock in
CCSM4 Decadal Prediction Experiments. Past Glob. Chang. Mag. 25, 41–46 (2017).

821
822

75. Yeager, S. G. et al. Reduced Southern Ocean warming enhances global skill and signal-tonoise in an eddy-resolving decadal prediction system. npj Clim. Atmos. Sci. 6, 107 (2023).

823
824

76. Boer, G. J., Kharin, V. V. & Merryfield, W. J. Decadal predictability and forecast skill. Clim
Dynam 41, 1817–1833 (2013).

31

825
826
827

77. Wu, X., Okumura, Y. M., DiNezio, P. N., Yeager, S. G. & Deser, C. The Equatorial PacificCold Tongue Bias in CESM1 and Its Influence on ENSO Forecasts. J. Clim. 35, 3261–3277(2022).

828
829

78. Capotondi, A. & Alexander, M. A. Rossby Waves in the Tropical North Pacific and Their
Role in Decadal Thermocline Variability. J. Phys. Oceanogr. 31, 3496–3515 (2001).

830
831

79. Capotondi, A., Alexander, M. A. & Deser, C. Why Are There Rossby Wave Maxima in the
Pacific at 10°S and 13°N? J. Phys. Oceanogr. 33, 1549–1563 (2003).

832
833
834

80. Sun, Z. et al. The impact of wind corrections and ocean-current influence on wind stressforcing on the modeling of Pacific North Equatorial Countercurrent. Ocean Model. 166, 101876(2021).

835
836

81. Bretherton, C. S., Smith, C. & Wallace, J. M. An Intercomparison of Methods for Finding
Coupled Patterns in Climate Data. J. Clim. 5, 541–560 (1992).

837
838

82. Zhao, Y., Lorenzo, E. D., Newman, M., Capotondi, A. & Stevenson, S. A Pacific Tropical
Decadal Variability Challenge for Climate Models. Geophys. Res. Lett. 50, (2023).

839
840
841

83. McGregor, S., Stuecker, M. F., Kajtar, J. B., England, M. H. & Collins, M. Model tropicalAtlantic biases underpin diminished Pacific decadal variability. Nat Clim Change 8, 493–498(2018).

842
843

84. Boer, G. J. et al. The Decadal Climate Prediction Project (DCPP) contribution to CMIP6.
Geosci Model Dev 9, 3751–3777 (2016).

844
845
846

85. Griffies, S. M. et al. OMIP contribution to CMIP6: experimental and diagnostic protocol forthe physical component of the Ocean Model Intercomparison Project. Geosci. Model Dev. 9,3231–3296 (2016).

847
848

86. Hurrell, J. W. et al. The Community Earth System Model: A Framework for Collaborative
Research. B Am Meteorol Soc 94, 1339–1360 (2013).

849
850
851

87. Kay, J. et al. The Community Earth System Model (CESM) Large Ensemble Project: A
Community Resource for Studying Climate Change in the Presence of Internal Climate
Variability. B Am Meteorol Soc 96, 1333–1349 (2015).

852
853

88. Neale, R., Chen, C. & Note …, G. A. Description of the NCAR community atmospheremodel (CAM 5.0). (2010).

854
855

89. Smith, R. et al. The Parallel Ocean Program (POP) reference manual, Ocean component of
the Community Climate System Model (CCSM). LANL Tech. Report LAUR-10-01853, 141.

856
857

90. Lawrence, D. M. et al. Parameterization improvements and functional and structuraladvances in Version 4 of the Community Land Model. J Adv Model Earth Sy 3, (2011).

32

858
859
860

91. Hunke, E. C. & Lipscomb, W. H. CICE: The Los Alamos sea ice model, documentation andsoftware, version 4.0. Los Alamos National Laboratory Tech. Rep. LA-CC-06-012, 76 pp.
(2008).

861
862
863

92. Yeager, S. G., Karspeck, A. R. & Danabasoglu, G. Predicted slowdown in the rate of Atlantic
sea ice loss: PREDICTED RATE OF SEA ICE LOSS. Geophys Res Lett 42, 10,704-10,713(2015).

864
865
866

93. Meehl, G. A. et al. The effects of bias, drift, and trends in calculating anomalies forevaluating skill of seasonal-to-decadal initialized climate predictions. Clim Dyn 59, 3373–3389(2022).

867
868

94. Huang, B. et al. Extended Reconstructed Sea Surface Temperature, Version 5 (ERSSTv5):
Upgrades, Validations, and Intercomparisons. J Climate 30, 8179–8205 (2017).

869
870

95. Rayner, N. et al. Global analyses of sea surface temperature, sea ice, and night marine airtemperature since the late nineteenth century. J Geophys Res 108, (2003).

871
872
873

96. Good, S. A., Martin, M. J. & Rayner, N. A. EN4: Quality controlled ocean temperature andsalinity profiles and monthly objective analyses with uncertainty estimates. J Geophys Res
Oceans 118, 6704–6716 (2013).

874
875
876

97. Gouretski, V. & Reseghetti, F. On depth and temperature biases in bathythermograph data:
Development of a new correction scheme based on analysis of a global ocean database. Deep Sea
Res Part Oceanogr Res Pap 57, 812–833 (2010).

877
878

98. Balmaseda, M., Mogensen, K. & Weaver, A. T. Evaluation of the ECMWF ocean reanalysis
system ORAS4. Quarterly Journal of the Royal Meteorological Society 139, 1132–1161 (2013).

879
880
881

99. McDougall, T. J., Jackett, D. R., Wright, D. G. & Feistel, R. Accurate and Computationally
Efficient Algorithms for Potential Temperature and Density of Seawater. J. Atmos. Ocean.
Technol. 20, 730–741 (2003).

882
883

100. Hersbach, H. et al. The ERA5 global reanalysis. Q. J. R. Meteorol. Soc. 146, 1999–2049(2020).

884
885
886

101. Slivinski, L. C. et al. Towards a more reliable historical reanalysis: Improvements forversion 3 of the Twentieth Century Reanalysis system. Q. J. R. Meteorol. Soc. 145, 2876–2908(2019).

887
888

102. Goddard, L. et al. A verification framework for interannual-to-decadal predictions
experiments. Clim Dynam 40, 245–272 (2013).

889
890
891

103. Large, W. G., Danabasoglu, G., Doney, S. C. & McWilliams, J. C. Sensitivity to Surface
Forcing and Boundary Layer Mixing in a Global Ocean Model: Annual-Mean Climatology. J.
Phys. Oceanogr. 27, 2418–2447 (1997).

33