<figure><img src="image_1.png" /></figure>

HAL
open science

Improving and Harmonizing El Niño Recharge Indices

Takeshi Izumo, Maxime Colin

► To cite this version:

Takeshi Izumo, Maxime Colin. Improving and Harmonizing El Niño Recharge Indices. Geophysical Research Letters, 2022, 49 (23), 10.1029/2022GL101003. hal-04018742

HAL Id: hal-04018742

https://hal.sorbonne-universite.fr/hal-04018742

Submitted on 8 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2022GL101003

**Key Points:**

* To clarify the debate on the Recharge Oscillator index, we develop an objective approach optimizing equation fit to observations

* The recharge index must be based on the slow component of sea level/thermocline depth, taken independently from the fast zonal tilt mode: this reconciles usual metrics

* The optimal index is this independent component averaged in the equatorial and southwestern Pacific. It is better suited for operational diagnostics

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

T. Izumo,
takeshi.izumo@ird.fr

**Citation:**

Izumo, T., & Colin, M. (2022). Improving and harmonizing El Niño recharge indices. *Geophysical Research Letters, 49*, e2022GL101003. https://doi.org/10.1029/2022GL101003

Received 24 AUG 2022
Accepted 14 OCT 2022

**Author Contributions:**

**Conceptualization:** Takeshi Izumo,
Maxime Colin

**Formal analysis:** Takeshi Izumo

**Investigation:** Takeshi Izumo, Maxime
Colin

**Methodology:** Takeshi Izumo, Maxime
Colin

**Validation:** Takeshi Izumo, Maxime
Colin

**Visualization:** Takeshi Izumo

**Writing – original draft:** Takeshi Izumo

**Writing – review & editing:** Takeshi
Izumo, Maxime Colin

© 2022. The Authors.

This is an open access article under
the terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_3.png" /></figure>

## Improving and Harmonizing El Niño Recharge Indices

Takeshi Izumo¹,² ID and Maxime Colin³,⁴,⁵ ID

¹EIO, IRD-Ifremer-UPF-ILM, Institut de Recherche pour le Développement (IRD), Université de la Polynésie Française (UPF), Tahiti, French Polynesia, ²Formerly at LOCEAN-IPSL, IRD, Sorbonne Université-CNRS-IRD-MNHN, Paris, France, ³Climate Change Research Centre, University of New South Wales (UNSW), Sydney, NSW, Australia, ⁴Laboratoire GEPASUD, University of French Polynesia (UPF), Tahiti, French Polynesia, ⁵Now at Leibniz Centre for Tropical Marine Research (ZMT), Bremen, Germany

**Abstract** El Niño Southern Oscillation (ENSO) is the leading mode of interannual climate variability, with large socioeconomic and environmental impacts. The main conceptual model for ENSO, the Recharge Oscillator (RO), considers two independent modes: the fast zonal tilt mode in phase with central-eastern Pacific Temperature ($T_e$), and the slow recharge mode in phase quadrature. However, usual indices (western or equatorial sea level/thermocline depth $h$) do not orthogonally isolate the slow recharge mode, leaving it correlated with $T_e$. Furthermore the optimal index is currently debated. Here, we develop an improved recharge index by objectively optimizing the RO equations fit to observations. (a) $T_e$-variability is regressed out, to build $h_{ind}$, statistically-independent from $T_e$. Capturing the pure recharge, $h_{ind}$ reconciles usual indices. (b) The optimum is equatorial plus southwestern Pacific $h_{ind\_eq+sw}$ (because of ENSO Ekman pumping meridional asymmetry). Using $h_{ind\_eq+sw}$, the RO becomes more consistent with observations. $h_{ind\_eq+sw}$ is more relevant for ENSO operational diagnostics.

**Plain Language Summary** El Niño and La Niña events have important impacts globally. A key element for long-lead forecasts is the recharge state of the tropical Pacific Ocean, as captured in the Recharge Oscillator (RO) conceptual model. The RO considers two independent modes of oceanic variability, a fast adjustment process and a slow recharge/discharge process. However, usual recharge indices mix these two modes of variability, and can thus lead to ambiguous operational diagnostics of the actual oceanic recharge state. Here we develop a better recharge index, independent of the fast mode, which reconciles typical indices and allows us to go beyond the current geographical debate on the optimal metrics. We use an objective approach optimizing the RO resemblance to observations to find the optimal index: the independent sea level (or equivalently thermocline depth) averaged over the equatorial and southwestern tropical Pacific. We recommend this simple and unambiguous index for El Niño operational forecasts diagnostics.

## 1. Introduction

The El Niño Southern Oscillation (ENSO) is the leading mode of climate interannual variability, with large socioeconomic and environmental impacts (e.g., Clarke, 2008; Jin et al., 2020; Neelin et al., 1998; Timmermann et al., 2018; Wang and Picaut, 2004 reviews). While ENSO predictability skill is now rather good at short leads, it has to be improved at longer leads (e.g., Barnston et al., 2012, 2019). One key element for long-lead seasonal forecasts is the recharge state of the tropical Pacific (be it in terms of Oceanic Heat Content (OHC), thermocline depth $h$ or sea level anomaly (SLA)), as it brings long oceanic memory across ENSO phases. Its role is formalized in the Recharge Oscillator (RO) conceptual model of ENSO (e.g., Clarke, 2010; Clarke et al., 2007; Jin, 1997a, 1997b; Wyrtki, 1985). For example, during a La Niña, easterlies favor a slow accumulation of OHC in the western and equatorial Pacific. This recharge will progressively favor positive Sea Surface Temperature (SST) anomalies in the central-eastern equatorial Pacific ($T_E$) and thus El Niño onset through the Bjerkness positive feedback. The El Niño event will in turn lead to a discharge favoring a reversal to La Niña conditions. The RO can thus well explain ENSO cyclic nature (cf. Section 3 for RO equations).

Yet there is a debate on the best recharge metric, for example, western or equatorial Pacific, sea level or thermocline depth. In the RO, the equatorial Pacific basin adjustment is separated into two independent modes: (a) the fast adjustment mode associated to a zonal tilt of the thermocline, in phase with zonal equatorial wind stress $\tau_x$ and $T_E$ (Figure 1a); (b) the slow adjustment recharge mode in phase quadrature with the fast mode (and thus with $\tau_x$ and $T_E$; Figure 1b). Several indices have been developed for this slow recharge mode that all bring predictability
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Schematics of the optimal orthonormal phase space for the Recharge Oscillator (RO). Panels a and b show respectively the fast adjustment zonal tilt mode dependent of $T_E$ (h regressed onto $T_E$, i.e., K) and the pure recharge mode ($h_{ind}$ regressed onto our default best recharge index $h_{ind\_ew+sw}$), independent of (uncorrelated to) $T_E$ (h unit: satellite SLA in cm; c.f. Sections 1 and 3 for details). Black arrows illustrate other used indices based on usual h. They are hence partly correlated to $T_E$.</figcaption></figure>

skill. Wyrtki (1985) and Jin (1997a, 1997b) originally focused on the western equatorial Pacific $h_w$, using equivalently SLA, OHC or the actual thermocline depth for h. Jin's sketch, Meinen and McPhaden (2000), and Burgers et al. (2005) focused on the mean equatorial band: $h_{eq}$ (e.g., Warm Water Volume (WWV) based on the 20°C isotherm depth Z20, or SLA, as usual thermocline depth proxies). $h_{eq}$ has become a widely-used ENSO recharge index. Theoretically, in the RO model, $h_{eq}$ is independent (i.e., uncorrelated at lag 0) from the fast tilt mode (and thus from $T_E$) within Jin (1997a, 1997b) approximations.

This independence/orthogonality property is essential, so that the recharge metric captures solely the pure recharge mode and does not mix it with the fast adjustment mode. But in observations, $h_{eq}$ (and its proxies) is not independent (Figure 1; upper-right panel of Figure S2 in Supporting Information S1). It is physically ambiguous, mixing the fast adjustment with the slow recharge mode, and thus potentially misleading: $h_{eq}$ is dominated by a short-term Ekman convergence leading to a temporary fast “increase” with El Niño westerlies (thus not a true long-term recharge in the RO sense), about two times larger (for SLA-based $h_{eq}$) than the long-term slow discharge expected from the RO theory. $h_{eq}$ is hence first an index of the fast equatorial Kelvin wave rather than of the RO long-term recharge/discharge process (Izumo, Lengaigne, et al., 2018; see also Neske and McGregor, 2018). Thus $h_{eq}$ is strongly positively correlated to the fast mode (and thus to $\tau_x$ and $T_E$). Hence, Izumo, Lengaigne, et al. (2018) recommended the use of a western Pacific index, $h_w$ (in agreement with Ballester, Petrova, et al., 2016; Boschat et al., 2013; Graham et al., 2015; Jin et al., 2020; Lai et al., 2015; Ramesh and Murtugudde, 2013; Petrova et al., 2017). Yet, such western Pacific OHC index is also not independent of $T_E$. It is partly negatively correlated to the fast zonal tilt mode (Figure 1a). $h_w$ could thus also lead to ambiguous diagnostics of the actual oceanic state.

Therefore the fast adjustment mode needs to be removed in order to obtain an operationally-useful indicator of the fully-isolated recharge process. Here we will hence develop an improved recharge index $h_{ind}$ independent of the fast mode. We will show that previously-defined indices converge when using $h_{ind}$. We will develop an objective approach optimizing the skill of the RO differential equations fit to observations for the ($T_E,h_{ind}$) pair to find the optimal averaging region for $h_{ind}$: $h_{ind\_eq+sw}$.

## 2. Data and Methods

Here we use classical observations and reanalyzes datasets, indices and statistical methods, detailed in Text S1 in Supporting Information S1 (statistics robust thanks to a sufficiently-large number of effective degrees of freedom, ~85 to ~140). The ENSO index $T_E$ is Niño3.4 relative Sea Surface Temperature (RSST, i.e.,
SST minus its 20°N–20°S tropical mean), as recommended by Izumo et al. (2020) and Van Oldenborgh et al. (2021) because atmospheric tropical deep convection interannual anomalies are rather related to RSST than to SST. Theoretically, SST is the variable directly involved for the term $F_1 * h$ in the $dT_E/dt$ equation (cf. Section 3). Yet, the recharge process is driven by windstress, itself directly driven by atmospheric deep convection and thus by RSST. Therefore, Niño3.4 RSST is better for the term $F_2 * T_E$ in the $dh/dt$ equation (Text S3 in Supporting Information S1). Using SST instead of RSST, or Niño3 instead of Niño3.4, makes $T_E$ slightly less correlated to equatorial Pacific $\tau_x$ (i.e., ocean-atmosphere coupling), but leads to very similar results (Table S2 in Supporting Information S1).

# 3. Improving the Recharge Oscillator Recharge Index

## 3.1. Traditional RO Framework Revisited

To derive RO $dT_E/dt$ tendency equation (see Jin et al., 2020 review for derivation), some physically-reasonable assumptions are used. (a) $\tau_x$ is proportional to $T_E$, that is, tropical deep convection and related $\tau_x$ respond quickly to $T_E$. (b) The fast oceanic response (i.e., quasi-instantaneous, timescales faster than ~2–3 months) to $\tau_x$ leads to a positive Bjerkness feedback term in the $dT_E/dt$ equation (through both the zonal advective and thermocline feedbacks) that is proportional to $\tau_x$ and thus to $T_E$: $R_{BJ\_o} T_E$. (c) Atmospheric fluxes are approximated as a weak Newtonian damping proportional to $T_E$: $-r_{damp\_o} T_E$. (d) A deepening of the thermocline depth $h$ related to a recharge favors positive $T_E$ on time scales longer than 2–3 months: $F_1 h$ (see Text S2 in Supporting Information S1 on mechanisms). Therefore:

$$ \frac{dT_E}{dt} = R_o T_E + F_1 h \qquad (1) $$

the first term on the right hand side representing the net effect of Bjerkness positive feedback and Newtonian damping ($R_o = R_{BJ\_o} - r_{damp\_o}$) and the second one, $F_1$, the recharge/discharge influence on $T_E$ (following Jin et al., 2020 notations).

Concerning $dh/dt$ equation, in the RO, the slow recharge mode response is the temporal integral of $\tau_x$ (e.g., Izumo et al., 2014). Negative $T_E$ associated with easterly anomalies will progressively recharge the equatorial Pacific (see Section 3.3 and Text S2 in Supporting Information S1). This is formalized as a term $-F_{2\_o} T_E$ in the $dh/dt$ equation:

$$ \frac{dh}{dt} = -F_{2\_o} T_E - \epsilon_o h \qquad (2) $$

the second term $-\epsilon_o h$ being formally a Newtonian damping, expected to be weak.

That is, the tendency equation for the vector $\begin{pmatrix} T_E \\ h \end{pmatrix}$ in matrix form is:

$$ \frac{d}{dt} \begin{pmatrix} T_E \\ h \end{pmatrix} = \begin{pmatrix} R_o & F_1 \\ -F_{2\_o} & -\epsilon_o \end{pmatrix} \begin{pmatrix} T_E \\ h \end{pmatrix} \qquad (3) $$

This RO linear equation is simple, but it remains unclear which geographical box to use for $h$. Furthermore, a standard multivariate linear regression fit minimizing rms error (second method in Burgers et al., 2005) gives significantly different coefficients for the classical normalized metrics $h_w$ and $h_{eq}$: $R_o = \epsilon_o = +0.06 \pm 0.04$ and $-0.18 \pm 0.05$ months$^{-1}$ for $h_w$ and $h_{eq}$ respectively, $F_1 = F_2 = 0.15 \pm 0.04$ and $0.25 \pm 0.05$ months$^{-1}$ respectively. The fact that $R_o$ and $\epsilon_o$ can change sign so easily for different classical recharge metrics makes it hard to interpret them physically. It evidences some physical inconsistency (Figure 1).

To resolve these geographical debate and physical issues objectively, an empirical way to find the best $h$ is to have no a priori on the best region and index. Thus, one may intuitively think of searching for the $h(x, y, t)$ region that statistically optimizes the skill of Equation 3. We will actually show later that we need two steps to construct
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Searching for the optimal recharge index h, representing the actual slow recharge mode independent of the fast mode in RO equations. Usual non-orthogonal basis ($T_E$, h) in the left set of panels, and orthogonal basis ($T_E$, $h_{ind}$) in the right one. Panels a and b show respectively the coefficients $R_o(x, y)$ and $F_1(x, y)$ for the multivariate regression of $dT_E/dt(t)$ onto $T_E(t)$ and $h(x, y, t)$ (Equation 4), multiplied by the STDs of $T_E$ (=1) and of $h(x, y, t)$ respectively, to measure their respective contributions to $dT_E/dt(t)$. Panels e and f are as a and c, but for coefficients for $dh(x, y, t)/dt$ regression. Panel i shows $r_{Te,h}$ skill using $h(x, y, t)$. Black line in panel j shows this skill for h averaged over: the equatorial band (5°N–5°S) with its western edge fixed to 120°E and its eastern edge varying, given by the x-axis. Red line is for a two-rectangle region, the classical equatorial box (5°N–5°S, 120°E–80°W) plus a southern 5°S–15°S box with the same 120°E fixed western edge and its eastern edge varying. The right set of panels is the same for the here-developed orthogonal basis ($T_E$, $h_{ind}$). The boxes overlaid represents the two-rectangle equatorial + SouthWest (eq + sw) optimal region finally chosen to average $h_{ind}$. This defines the suggested improved index: $h_{ind\_eq+sw}$. See text for details.</figcaption></figure>

a better index: (a) independence to $T_E$, (b) geographical optimization of the averaging region. Let us first try the geographical optimization alone:

$$ \frac{d}{dt} \begin{pmatrix} T_E(t) \\ h(x, y, t) \end{pmatrix} = \begin{pmatrix} R_o(x, y) & F_1(x, y) \\ -F_{2_o}(x, y) & -\epsilon_o(x, y) \end{pmatrix} \begin{pmatrix} T_E(t) \\ h(x, y, t) \end{pmatrix} + \begin{pmatrix} \text{residual}_T(x, y, t) \\ \text{residual}_h(x, y, t) \end{pmatrix} \quad (4) $$

In practice we do at each spatial point (x, y) a fit with two multivariate linear regressions, one for Equation 1 and one for Equation 2, with Pearson correlation skill $r_{Te}$ and $r_h$ respectively. The squareroot of the RV-coefficient (Rho-Vectoriel; Escoufier & Robert, 1976), "$r_{Te,h}$," is the equivalent of the Pearson correlation for matrix form, that is, a measure of the skill for the 2D Equation 4. Where the skill is the highest (residuals variance minimized) should inform us of the best h region to capture the RO processes. The highest skill is in the equatorial and southwestern tropical Pacific (Figure 2i), suggesting that the best region for averaging $h(x, y, t)$ could combine these two regions. However, the regression coefficients contributions have a puzzling spatial distribution (Figures 2a, 2b, 2e, and 2f). They are highly spatially-correlated with large opposing signs. In the central Equatorial Pacific, $T_E$ and h influences on $dT_E/dt$ would be very large, but negative for $R_o$ and positive for $F_1$ respectively, which is unphysical. This is actually a statistical artifact because $T_E(t)$ and $h(x, y, t)$ are not statistically independent (i.e., not orthogonal), which makes it difficult to interpret physically the regression coefficients. It therefore suggests that the ($T_E$, h) basis is not ideal.

## 3.2. Toward New RO Tendency Equations for $T_E$ and Independent $h_{ind}$

A large component of h responds rapidly to $\tau_x$ and $T_E$. This fast response component of h is correlated to (i.e., linearly dependent of) $T_E$ at timescales longer than ~2–3 months in the equatorial Pacific: Figure 1a shows the downwelling pattern along the central to eastern equatorial Pacific and the zonal seesaw pattern in the western Pacific, that is, the fast adjustment mode. We suggest to regress out from $h(x, y, t)$ this rapidly-responding component correlated to $T_E(t)$, so as to only focus on the slowly-responding independent component of h, hereafter $h_{ind}$.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Observed trajectories of the RO in the usual ($T_E$, $h$) space (upper panels) and new orthonormal ($T_E$, $h_{ind}$) space (lower panels), for classical $h_{eq}$ and $h_w$ indices (1st and 2nd columns), and optimal $h_{eq+sw}$ index (3rd column). (a) Trajectory of the system for the pair of coordinates ($T_E$, $h_{eq}$), $T_E$ in horizontal axis, $h$ in vertical axis, each month being a point, with year indicated in color, for ORAS5 SLA. Other panels are similar, but for other recharge indices. Indices are all normalised. The shared variance (square of correlation) between recharge indices is shown in parenthesis. (d) schematised elliptical trajectory if the recharge index is positively correlated to $T_E$, such as $h_{eq}$. Such index is ambiguous. For example, a positive anomaly does not necessarily correspond to an actual recharge state (blue dot). (h) circular trajectory for the recharge index independent component, consistent with ideal RO trajectory, and unambiguous. Color dots indicate different positions in the cycle.</figcaption></figure>

(Figure 3b of Izumo, Lengaigne, et al., 2018). $h_{ind}$ corresponds to the slow recharge mode related to the slower basin adjustment in disequilibrium with wind stress (e.g., Alory & Delcroix, 2002; Clarke, 2010; Fedorov, 2010; Jin, 1997a, 1997b; Masuda et al., 2009; Thual et al., 2013; Zhu et al., 2017; Izumo, Lengaigne, et al., 2018). Thus, we can formally decompose $h$:

$$h(x, y, t) = h_{\text{fast\_mode(dependent,corr.to Te)}} + h_{\text{slow\_recharge\_mode(independent,uncorr.to Te)}}$$

$$h(x, y, t) = K(x, y) T_E(t) + h_{\text{ind}}(x, y, t) \qquad (5)$$

$K(x, y)$ being the regression coefficient of $h$ onto $T_E$ (Figure 1a).

$h_{ind}$ is the pure recharge component, that allows us to describe the system with an orthonormal basis: ($T_E$, $h_{ind}$) (Figure 1). An orthonormal basis is more acceptable physically and mathematically.

With this linear transform, Equation 4 can then be rewritten as:

$$\frac{d}{dt} \begin{pmatrix} T_E(t) \\ h_{\text{ind}}(x, y, t) \end{pmatrix} = \begin{pmatrix} R(x, y) & F_1(x, y) \\ -F_2(x, y) & -\epsilon(x, y) \end{pmatrix} \begin{pmatrix} T_E(t) \\ h_{\text{ind}}(x, y, t) \end{pmatrix} + \begin{pmatrix} \text{residual}_T(x, y, t) \\ \text{residual}_{h_{\text{ind}}}(x, y, t) \end{pmatrix} \qquad (6)$$

with $R = R_0 + K F_1$, $F_2 = F_{2_o} + \epsilon_0 K + K R_0 + K^2 F_1$, $\epsilon = \epsilon_0 + K F_1$ and $\text{residual}_{h_{ind}} = \text{residual}_h - K \text{residual}_T$.

So the $F_1$ term representing the recharge influence on $T_E$ remains the same for $h_{ind}$, and the physics behind is also the same. $F_2$ represents the influence of easterly anomalies related to $T_E$ on the $h_{ind}$ recharge. It differs from $F_{2_o}$.

Figures 2c, 2d, 2g, and 2h show the coefficients obtained from the multivariate regression onto ($T_E$, $h_{ind}$). $R$ and $\epsilon$ are uniformly negligible (~0.00 ± 0.03 months⁻¹). $F_1$ and $F_2$ vary spatially similarly, with highest values in the equatorial and southwest Pacific corresponding to the highest skill (Figure 2k and Figure S1 in Supporting Information S1). Hence the best $h_{ind}$ region should combine these two regions. The picture is also clearer and more consistent with RO theory: $\epsilon$ is negligible, consistent with a weak damping (due mainly to oceanic mixing,
e.g., Fedorov, 2010). $r_{Te,h\_ind}$ is larger than $r_{Te,h}$ (Figure 2i), further confirming that ($T_E$, $h_{ind}$) is a better basis of RO phase space.

Here we explain the spatial patterns of Figure 2's various panels (full explanation in Text S2 in Supporting Information S1), robust among datasets and periods (cf. Figure S2 in Supporting Information S1), and why there are similarities among some of them. $F_2$ map shows us how $h_{ind}(x, y, t)$ would look like if ENSO windstress anomalies would blow for a long time, for example, because of long-lasting La Niña conditions. The slow recharge is as expected in the western and central equatorial Pacific, through: (a) downwelling equatorial Rossby waves to the west (Jin, 1997a, 1997b; Wyrtki, 1985), and off-equatorial ones in the southwest; (b) upwelling equatorial Kelvin waves to the east forcing coastal Kelvin waves propagating poleward along the eastern boundary and thus leading to a leakage of negative OHC anomalies toward the poles along the eastern boundary (Izumo, Lengaigne, et al., 2018; Wyrtki, 1985).

Interestingly the recharge is also in the southwest, because of asymmetric Ekman pumping (poleward shift of the South Pacific Convergence Zone; SPCZ) forcing locally downwelling (in the La Niña case) and thus slow off-equatorial downwelling Rossby waves progressively recharging the southwestern Pacific. Note that the western boundary coastline meridional asymmetry would conversely favor a larger northwest recharge, as shown by sensitivity experiments with the LCS model (Linear Continuously Stratified model; McCreary, 1980; configuration of Izumo, Lengaigne, et al., 2018; Figure S3 in Supporting Information S1). $F_1$ physically represents the slow recharge mode influence on $T_E$ through several mechanisms (Text S2 in Supporting Information S1).

$R_o$, $\epsilon_o$ and $F_{2\_o}$ spatial patterns can be explained as follow: $R_o \approx \epsilon_o \approx -KF_1$ and $F_{2\_o} \approx F_2 + K^2F_1$ (as $\epsilon \approx 0$ and $R \approx 0$). These relationships explain the paradox (pointed out in Section 3.1) with large positive and negative $R_o$ and $\epsilon_o$ values found for classical $h_w$ and $h_{eq}$ respectively: $R_o \approx \epsilon_o > 0$ (because $K < 0$ in the west) for $h_w$ and $R_o \approx \epsilon_o < 0$ (because $K > 0$) for $h_{eq}$, even if the actual damping $\epsilon$ and net feedback $R$ are weak. So the “damping” term ($\epsilon_o$) and the “positive feedback” ($R_o$) would be artificially large. This further confirms that we should use $h_{ind}$ rather than the full $h$ polluted by the fast mode that biases the RO model.

Such transforms from observable variables to new variables that are more relevant physically are often done in physics and geosciences, for example, to create potential temperature, relative SST, rotating PCs (Takahashi et al., 2011), decomposition into baroclinic modes, spherical harmonics... We just want to clearly isolate the independent influence of the slow recharge mode from the fast tilt mode, and to optimize the RO understanding and metrics (we are not trying to add new physical processes).

Considering $h_{ind}$ instead of $h$ allows us to go beyond the RO metrics debate, by reconciling indices: they become much closer when considering their independent component. For example, $h_{ind\_w}$ and $h_{ind\_eq}$ are much better correlated ($r^2$ increases from 10% to 69%) with now similar coefficients (Figures S4, S5 and Table S2 in Supporting Information S1). And Z20 and SLA-based indices become even closer (Figure S6 in Supporting Information S1). Trajectories in the ($T_E$, $h_{ind}$) phase diagram (Dommenget & Al-Ansari, 2022; Kessler, 2002) also become closer (Figure 3; shown here for longer ORAS5 SLA data set to have a larger density; see Figures S7 and S8 in Supporting Information S1 for satellite SLA and ORAS5 Z20). The transform is geometrically simply a linear transform of coordinates from the ($T_E$, $h$) to the ($T_E$, $h_{ind}$) orthogonal coordinate space. Trajectories in the latter space are more circular, closer to idealized RO model circular trajectory.

To conclude, $h_{ind}$ reconciles various recharge metrics such as $h_{ind\_w}$ and $h_{ind\_eq}$. Here below we would like to further improve them by objectively finding $h_{ind}$ optimal averaging region.

### 3.3. Objectively Finding the Optimal Region: Equatorial + South-West ($h_{ind\_eq+sw}$)

The optimal region for $h_{ind}(x, y, t)$ will be the one optimizing the skill $r_{Te,h\_ind}$ of the RO tendency Equation 6

for $\begin{pmatrix} T_E \\ h_{ind} \end{pmatrix}$:

$$ \frac{d}{dt} \begin{pmatrix} T_E \\ h_{ind} \end{pmatrix} = \begin{pmatrix} R & F_1 \\ -F_2 & -\epsilon \end{pmatrix} \begin{pmatrix} T_E \\ h_{ind} \end{pmatrix} \quad (7) $$
$r_{Te,h\_ind}$ (Figure 2k) is the largest in the equatorial (~5°N–5°S) and southwest (~5°S–15°S) Pacific. The combination of these two boxes should be our optimal region, and $r_{Te,h\_ind}$ should be even better thanks to the spatial averaging (reducing noise). To choose objectively the optimal averaging region, we have tested several options. We have first tested rectangular boxes, for example, by averaging $h_{ind}(x, y, t)$ from the usual fixed 120°E western end to a varying eastern edge, for various latitudinal bands. For the classical 5°N–5°S band, the best skill ($r_{Te,h\_ind}$ = 0.61) is found for an eastern edge around ~90–80°W (Figure 2l, black line), consistent with classical $h_{eq}$ (5°N–5°S, 120°E–80°W). Another relevant band is the 5°N–15°S, 120°E to ~150°W, leading to a similarly high skill (not shown). A “hybrid” choice with two rectangles is even better. Adding to the optimal classical equatorial band (5°N–5°S, 120°E–80°W) a second box in the southwest, along the 5°S–15°S band, from the same 120°E western edge to a varying eastern edge, further improves the skill to $r_{Te,h\_ind}$ = 0.69 for an optimal eastern edge around ~170°W (Figure 2l, red line). We have similarly tested all the other edges of this two-rectangle region, also testing ORAS5 SLA and Z20 (Figures S9a, S9b, and S10 in Supporting Information S1). This is our best allround and sufficiently-simple index, hereafter $h_{ind\_eq+sw}$ (note that $h_{ind\_eq+sw}$ = $h_{eq+sw\_ind}$, the regression being linear). Interestingly, this choice also conveniently minimizes the correlation between $h_{eq+sw}$ and $T_E$ (Figures S9c and S9d in Supporting Information S1), making $h_{eq+sw}$ closer to $h_{ind\_eq+sw}$. The skill increase from $h_{ind\_eq}$ to $h_{ind\_eq+sw}$ is statistically significant, more thanks to $r_{h\_ind}$ ($p$ = 0.01|0.01|0.01) than to $r_{Te}$ ($p$ = 0.09|0.01|0.23) for satSLA/ORAS5_SLA/ORAS5_Z20 respectively). And trajectories are even smoother (Figure 3g). Including the southwest agrees with Izumo et al. (2018a; see also Santoso et al., 2017; Ramesh and Murtugudde, 2013). $h_{ind\_eq+sw}$ is now our default best recharge index.

Using normalized $h_{ind\_eq+sw}$, we obtain: $F_1 \approx F_2 \approx 0.17 \pm 0.03$, $R \approx \varepsilon \approx 0.00 \pm 0.03$ months⁻¹. These coefficients are robust among datasets and periods (Table S1 in Supporting Information S1). $R$ and $\varepsilon$ are negligible. This reduces the parameter space. And the RO equations system has the form of a harmonic oscillator excited by noise (Burgers et al., 2005), with $d²X(t)/dt² \approx -F_1F_2X(t), X(t)$ being $T_E$ or $h_{ind}$, and angular frequency being the Wyrtki index $W = (F_1F_2)^{1/2} \approx F_1 \approx F_2$ (i.e., reasonable eigenperiod of ~3.1 year). Even if theoretically undamped, the oscillator is actually still damped because residual is not a pure red noise but includes non-linear terms neglected in our first order linear approximation (Table S2 in Supporting Information S1).

# 4. Conclusion

## 4.1. Summary

Here we have defined a simple Pacific recharge index $h_{ind}$ independent of the fast mode, by regressing out $T_E$-related variability. It unambiguously represents the slow recharge mode, with more physically-consistent RO parameters, conversely to classical recharge indices. $h_{ind}$ harmonizes recharge indices: by taking $h_{ind}$, averages over the usual western/equatorial Pacific regions based on Z20/OHC/SLA have much more similar time series, equation parameters and phase trajectories, so that all indices converge to a single one. We have also objectively searched for the optimal averaging region to have the most realistic RO tendency equations. The optimal index is $h_{ind\_eq+sw}$ averaged over the classical equatorial band (5°N–5°S, 120°E–80°W), extended in the southwest until 15°S (5°S–15°S, 120°E–170°W).

In practice, obtaining $h_{ind\_eq+sw}$ is straightforward: (a) average SLA over “eq + SW” box and then normalize, (b) remove its dependent part $KT_E: h_{ind\_eq+sw} = h_{eq+sw} - K_{eq+sw} T_E$, with regression coefficient $K_{eq+sw} \approx 0.26$ (details in Text S1 in Supporting Information S1).

## 4.2. Discussion

For ENSO operational diagnostics, the ($T_E, h_{ind\_eq+sw}$) basis is more relevant to describe the system trajectory than the usual ($T_E, h_{eq}$) and ($T_E, h_w$) pairs (Figure 3). If the latter are used, a situation with anomalous $h$ is ambiguous. Independent $h_{ind\_eq+sw}$ better represents actual precursory recharge anomalies. Only with $h_{ind}$ sign (recharged/discharged state) can one get $dT_E/dt$ sign directly (while $h$ sign is not sufficient, as the fast tilt mode signal could blur a weaker long-term build-up; see yellow/blue dots in idealized schematics Figures 3d and 3h).

To compare $h_{ind\_eq}$ and $h_{ind\_eq+sw}$ forecasting skills, we have done a preliminary assessment by using the simple multivariate linear regression model combining the recharge index, Indian Ocean Dipole (IOD) index and $T_E$, all in September-November, to hindcast $T_E$ peak in November-January 14 months later. Adding the Southwest
improves the skill for all datasets/periods, more clearly for SLA than for Z20 (Table S3 in Supporting Information S1; since we have statistically-significant contributions from $h_{ind\_eq + sw}$ and IOD, but not from $T_E$ itself, this new recharge index and updated datasets confirm earlier studies of Izumo et al., 2010, 2014, 2016; Dayan et al., 2014; Jourdain et al., 2016).

Results are robust among datasets. Only some subtle second order differences remain between SLA and Z20 (OHC being inbetween) for $h_{ind}$, likely related to different weighting of the first baroclinic modes (Figures S2, S6, S8, and S10 in Supporting Information S1). Note that we are not trying to address “how to best estimate the thermocline depth” (Vijayeta et al., 2020), but rather “how to best isolate the recharge mode.” SLA and Z20 each have their own advantages with presently-available datasets: Z20 seems better in a perfectly-observed and non-warming ocean (Tables S1 and S3 in Supporting Information S1), but SLA is observed globally by satellite in near real-time and is more available in climate models outputs (e.g., Coupled Model Intercomparison Project, CMIP).

Here we have neglected possible seasonal cycles of parameters and asymmetries/non-linearities for the sake of simplicity. Knowing ENSO seasonal phase-locking, taking into account such seasonal cycles could be one next step. We could include asymmetrical/non-linear terms (Table S2 in Supporting Information S1). To have a simple $h_{ind}$ definition, we have also neglected a weak non-linearity in the fast response of $h$ to $T_E$ (Figure 3c; Figures S7c and S8c in Supporting Information S1) likely due to $\tau_x$ anomalies longitudinal position, further to the east during El Niño than during La Niña, favoring a larger discharge. This non-linearity would lead to a quadratic term in the $dh_{ind}/dt$ equation.

Our method is based on a statistical optimisation of regression coefficients from observations to fit the RO model, by projecting on an orthogonal basis, which is conceptually satisfying, and practically necessary to determine eigenmodes. But correlations do not mean causality, and this study therefore does not aim to challenge RO physical validity. Although using the physical $h$ variable is necessary to infer causality and derive RO equations, we still argue that having orthogonal coordinates thanks to $h_{ind}$ allows for clearer recharge diagnostics, and a properly isolated recharge mode.

Some studies, using usual $h$ indices, have put in question the RO and suggested that the delayed oscillator (DO, Suarez and Schopf, 1988; Battisti & Hirst, 1989) is more realistic (e.g., Graham et al., 2015; Linz et al., 2014). One reason for this “RO versus DO” debate could be the misleading character of usual $h$ indices because of fast tilt mode (dependent) component, which tends to artificially reduce RO skill. It would be interesting to compare the RO to the other ENSO oscillators, using $h_{ind\_eq + sw}$ instead, to clarify whether theories, climate models and observations would better agree.

## Acknowledgments

MC would like to thank Steven Sherwood for his essential support, supervision, and mentoring, in particular throughout the pandemic, which made this collaboration possible. We would like to thank the two reviewers, notably Dr. Dietmar Dommenget for his critical thinking and interesting discussions, who have pushed us to improve, clarify and strengthen the results, as well as Dr. Matthieu Lengaigne, Dr. Jérôme Vialard, and Pr. Fei-Fei Jin for discussions. T.I. and M.C. would like to thank UPF for hosting M.C. because of COVID-related closed Australian borders; this latter fact was paradoxically essential for allowing us to make this work mature. We also thank UNSW Sydney for allowing a flexible work agreement for MC during the pandemic. M. Colin was first supported by the Australian ARC Grant FL150100035 and kindly hosted by UPF, and then directly funded by the ZMT's core budget. T. Izumo was funded mainly by IRD core budget, with some funding from ANR ARISE project.

## Data Availability Statement

All data used is open data: OISSTv2 (https://psl.noaa.gov/data/gridded/data.noaa.oisst.v2.html), CMAP1 (https://psl.noaa.gov/data/gridded/data.cmap.html), ERA5-ORAS5 (https://www.ecmwf.int), Copernicus SLA (https://doi.org/10.48670/moi-00148), and HadiSST (https://www.metoffice.gov.uk/hadobs/hadisst/). We acknowledge use of NOAA pyferret open source software for analyses and figures.

## References

Alory, G., & Delcroix, T. (2002). Interannual sea level changes and associated mass transports in the tropical Pacific from TOPEX/Poseidon data and linear model results (1964–1999). *Journal of Geophysical Research*, 107(C10), 17–21. https://doi.org/10.1029/2001jc001067

Ballester, J., Petrova, D., Bordoni, S., Cash, B., García-Díez, M., & Rodó, X. (2016). Sensitivity of El Niño intensity and timing to preceding subsurface heat magnitude. *Scientific Reports*, 6(1), 1–9. https://doi.org/10.1038/srep36344

Barnston, A. G., Tippett, M. K., L'Heureux, M. L., Li, S., & DeWitt, D. G. (2012). Skill of real-time seasonal ENSO model predictions during 2002–11: Is our capability increasing? *Bulletin of the American Meteorological Society*, 93(5), 631–651. https://doi.org/10.1175/bams-d-11-00111.2

Barnston, A. G., Tippett, M. K., Ranganathan, M., & L'Heureux, M. L. (2019). Deterministic skill of ENSO predictions from the North American multimodel ensemble. *Climate Dynamics*, 53(12), 7215–7234. https://doi.org/10.1007/s00382-017-3603-3

Battisti, D. S., & Hirst, A. C. (1989). Interannual variability in a tropical atmosphere-ocean model: Influence of the basic state, ocean geometry and nonlinearity. *Journal of the Atmospheric Sciences*, 46(12), 1687–1712. https://doi.org/10.1175/1520-0469(1989)046<1687:iviata>2.0.co;2

Boschat, G., Terray, P., & Masson, S. (2013). Extratropical forcing of ENSO. *Geophysical Research Letters*, 40(8), 1605–1611. https://doi.org/10.1002/grl.50229

Burgers, G., Jin, F.-F., & van Oldenborgh, G. J. (2005). The simplest ENSO recharge oscillator. *Geophysical Research Letters*, 32(13), L13706. https://doi.org/10.1029/2005GL022951
Clarke, A. J. (2008). *An introduction to the dynamics of El Niño and the southern oscillation*. Elsevier.

Clarke, A. J. (2010). Analytical theory for the quasi-steady and low-frequency equatorial ocean response to wind forcing: The “tilt” and “warm water volume” modes. *Journal of Physical Oceanography*, 40(1), 121–137. https://doi.org/10.1175/2009jpo4263.1

Clarke, A. J., Van Gorder, S., & Colantuono, G. (2007). Wind stress curl and ENSO discharge/recharge in the equatorial Pacific. *Journal of Physical Oceanography*, 37(4), 1077–1091. https://doi.org/10.1175/jpo3035.1

Dayan, H., Vialard, J., Izumo, T., & Lengaigne, M. (2014). Does sea surface temperature outside the tropical Pacific contribute to enhanced ENSO predictability? *Climate Dynamics*, 43(5–6), 1311–1325. https://doi.org/10.1007/s00382-013-1946-y

Dommenget, B. D., & Al-Ansari, M. (2022). Asymmetries in the ENSO phase space. *Climate Dynamics*, 1–20. https://doi.org/10.1007/s00382-022-06392-0

Escoufier, R., & Robert, P. (1976). A unifying tool for linear multivariate statistical methods: The Rv-coefficient. *Applied Statistics*, 25(3), 257–265. https://doi.org/10.2307/2347233

Fedorov, A. V. (2010). Ocean response to wind variations, warm water volume, and simple models of ENSO in the low-frequency approximation. *Journal of Climate*, 23(14), 3855–3873. https://doi.org/10.1175/2010jcli3044.1

Graham, F. S., Brown, J. N., Wittenberg, A. T., & Holbrook, N. J. (2015). Reassessing conceptual models of ENSO. *Journal of Climate*, 28(23), 9121–9142. https://doi.org/10.1175/jcli-d-14-00812.1

Izumo, T., Lengaigne, M., Vialard, J., Luo, J. J., Yamagata, T., & Madec, G. (2014). Influence of Indian Ocean Dipole and Pacific recharge on following year's El Niño: Interdecadal robustness. *Climate Dynamics*, 42(1–2), 291–310. https://doi.org/10.1007/s00382-012-1628-1

Izumo, T., Lengaigne, M., Vialard, J., Suresh, I., & Planton, Y. (2018). On the physical interpretation of the lead relation between warm water volume and the El Niño southern oscillation. *Climate Dynamics*, 52(5), 2923–2942. https://doi.org/10.1007/s00382-018-4313-1

Izumo, T., Vialard, J., Dayan, H., Lengaigne, M., & Suresh, I. (2016). A simple estimation of equatorial Pacific response from windstress to untangle Indian Ocean Dipole and Basin influences on El Niño. *Climate Dynamics*, 46(7–8), 2247–2268. https://doi.org/10.1007/s00382-015-2700-4

Izumo, T., Vialard, J., Lengaigne, M., de Boyer Montegut, C., Behera, S. K., Luo, J. J., et al. (2010). Influence of the Indian Ocean Dipole on following year's El Niño. *Nature Geoscience*, 3, 168–172. https://doi.org/10.1038/ngeo760

Izumo, T., Vialard, J., Lengaigne, M., & Suresh, I. (2020). Relevance of relative sea surface temperature for tropical rainfall interannual variability. *Geophysical Research Letters*, 47(3), e2019GL086182. https://doi.org/10.1029/2019gl086182

Jin, F. F. (1997a). An equatorial ocean recharge paradigm for ENSO. Part I: Conceptual model. *Journal of the Atmospheric Sciences*, 54(7), 811–829. https://doi.org/10.1175/1520-0469(1997)054<0811:aeorpf>2.0.co;2

Jin, F. F. (1997b). An equatorial ocean recharge paradigm for ENSO. Part II: A stripped-down coupled model. *Journal of the Atmospheric Sciences*, 54(7), 830–847. https://doi.org/10.1175/1520-0469(1997)054<0830:aeorpf>2.0.co;2

Jin, F. F., Chen, H. C., Zhao, S., Hayashi, M., Karamperidou, C., Stuecker, M. F., et al. (2020). Simple ENSO models. In *El Niño Southern Oscillation in a changing climate* (pp. 119–151).

Jourdain, N., Lengaigne, M., Vialard, J., Izumo, T., & Sen Gupta, A. (2016). Further insights on the influence of the Indian Ocean Dipole on the following year's ENSO from observations and CMIP5 models. *Journal of Climate*, 29(2), 637–658. https://doi.org/10.1175/JCLI-D-15-0481.1

Kessler, W. S. (2002). Is ENSO a cycle or a series of events? *Geophysical Research Letters*, 29(23), 40–41. https://doi.org/10.1029/2002gl015924

Lai, A. W. C., Herzog, M., & Graf, H. F. (2015). Two key parameters for the El Niño continuum: Zonal wind anomalies and Western Pacific sub-surface potential temperature. *Climate Dynamics*, 45(11), 3461–3480. https://doi.org/10.1007/s00382-015-2550-0

Linz, M., Tziperman, E., & MacMartin, D. G. (2014). Process-based analysis of climate model ENSO simulations: Intermodel consistency and compensating errors. *Journal of Geophysical Research*, 119(12), 7396–7409. https://doi.org/10.1002/2013JD021415

Masuda, S., Awaji, T., Toyoda, T., Shikama, Y., & Ishikawa, Y. (2009). Temporal evolution of the equatorial thermocline associated with the 1991–2006 ENSO. *Journal of Geophysical Research*, 114(C3), C03015. https://doi.org/10.1029/2008jc004953

McCreary, J. P. (1980). Modelling wind-driven ocean circulation. *JIMAR 80-0029, HIG 80-3*. p. 64. University of Hawaii. Honolulu.

Meinen, C. S., & McPhaden, M. J. (2000). Observations of warm water volume changes in the equatorial Pacific and their relationship to El Niño and La Niña. *Journal of Climate*, 13(20), 3551–3559. https://doi.org/10.1175/1520-0442(2000)013<3551:oowwvc>2.0.co;2

Neelin, J. D., Battisti, D. S., Hirst, A. C., Jin, F. F., Wakata, Y., Yamagata, T., & Zebiak, S. E. (1998). ENSO theory. *Journal of Geophysical Research*, 103(C7), 14261–14290. https://doi.org/10.1029/97jc03424

Neske, S., & McGregor, S. (2018). Understanding the warm water volume precursor for ENSO events and its interdecadal variation. *Geophysical Research Letters*, 45(3), 1577–1585. https://doi.org/10.1002/2017GL076439

Petrova, D., Koopman, S. J., Ballester, J., & Rodó, X. (2017). Improving the long-lead predictability of El Niño using a novel forecasting scheme based on a dynamic components model. *Climate Dynamics*, 48(3–4), 1249–1276. https://doi.org/10.1007/s00382-016-3139-y

Ramesh, N., & Murtugudde, R. (2013). All flavours of El Niño have similar early subsurface origins. *Nature Climate Change*, 3(1), 42–46. https://doi.org/10.1038/nclimate1600

Santoso, A., Mcphaden, M. J., & Cai, W. (2017). The defining characteristics of ENSO extremes and the strong 2015/2016 El Niño. *Reviews of Geophysics*, 55(4), 1079–1129. https://doi.org/10.1002/2017rg000560

Suarez, M. J., & Schopf, P. S. (1988). A delayed action oscillator for ENSO. *Journal of the Atmospheric Sciences*, 45(21), 3283–3287. https://doi.org/10.1175/1520-0469(1988)045<3283:adaoafe>2.0.co;2

Takahashi, K., Montecinos, A., Goubanova, K., & Dewitte, B. (2011). ENSO regimes: Reinterpreting the canonical and modoki El Niño. *Geophysical Research Letters*, 38(10), 1–5. https://doi.org/10.1029/2011GL047364

Thual, S., Dewitte, B., Ayoub, N., & Thual, O. (2013). An asymptotic expansion for the recharge–discharge model of ENSO. *Journal of Physical Oceanography*, 43(7), 1407–1416. https://doi.org/10.1175/JPO-D-12-0161.1

Timmermann, A., An, S. I., Kug, J. S., Jin, F. F., Cai, W., Capotondi, A., et al. (2018). El Niño–southern oscillation complexity. *Nature*, 559.

Van Oldenborgh, G. J., Hendon, H., Stockdale, T., L'Heureux, M., De Perez, E. C., Singh, R., & Van Aalst, M. (2021). Defining El Niño indices in a warming climate. *Environmental Research Letters*, 16(4), 044003. https://doi.org/10.1088/1748-9326/abe9ed

Vijayeta, ., Dommenget, D., & McGregor, S. (2020). *ENSO in a changing climate simulations*. PhD thesis. Monash University.

Wang, C., & Picaut, J. (2004). Understanding ENSO physics—A review. *Earth's Climate: The Ocean–Atmosphere Interaction, Geophys. Monogr.*, 147, 21–48.

Wyrtki, K. (1985). Water displacements in the Pacific and the Genesis of El Niño cycles. *Journal of Geophysical Research*, 90(C4), 7129–7132. https://doi.org/10.1029/jc090ic04p07129

Zhu, X., Greatbatch, R. J., & Claus, M. (2017). Interannual variability of tropical Pacific sea level from 1993 to 2014. *Journal of Geophysical Research: Oceans*, 122(1), 602–616. https://doi.org/10.1002/2016jc012347
## References From the Supporting Information

An, S. I., Tziperman, E., Okumura, Y. M., & Li, T. (2020). ENSO irregularity and asymmetry. In *El Niño Southern Oscillation in a changing climate* (pp. 153–172).

Ballester, J., Bordoni, S., Petrova, D., & Rodó, X. (2015). On the dynamical mechanisms explaining the western Pacific subsurface temperature buildup leading to ENSO events. *Geophysical Research Letters*, *42*(8), 2961–2967. https://doi.org/10.1002/2015gl063701

Ballester, J., Bordoni, S., Petrova, D., & Rodó, X. (2016). Heat advection processes leading to El Niño events as depicted by an ensemble of ocean assimilation products. *Journal of Geophysical Research: Oceans*, *121*(6), 3710–3729. https://doi.org/10.1002/2016jc011718

Bretherton, C. S., Widmann, M., Dymnikov, V. P., Wallace, J. M., & Bladé, I. (1999). The effective number of spatial degrees of freedom of a time-varying field. *Journal of Climate*, *12*(7), 1990–2009. https://doi.org/10.1175/1520-0442(1999)012<1990:tenosd>2.0.co;2

Capotondi, A., Wittenberg, A. T., Newman, M., Di Lorenzo, E., Yu, J. Y., Braconnot, P., et al. (2015). Understanding ENSO diversity. *Bulletin of the American Meteorological Society*, *96*(6), 921–938. https://doi.org/10.1175/bams-d-13-00117.1

Cibot, C., Maisonnave, E., Terray, L., & Dewitte, B. (2005). Mechanisms of tropical Pacific interannual-to-decadal variability in the ARPEGE/ORCA global coupled model. *Climate Dynamics*, *24*(7), 823–842. https://doi.org/10.1007/s00382-004-0513-y

Gadgil, S., Joseph, P. V., & Joshi, N. V. (1984). Ocean–atmosphere coupling over monsoon regions. *Nature*, *312*(5990), 141–143. https://doi.org/10.1038/312141a0

Gasparin, F., & Roemmich, D. (2017). The seasonal march of the equatorial Pacific upper-ocean and its El Niño variability. *Progress in Oceanography*, *156*, 1–16. https://doi.org/10.1016/j.pocology.2017.05.010

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., et al. (2020). The ERA5 global reanalysis. *Quarterly Journal of the Royal Meteorological Society*, *146*(730), 1999–2049. https://doi.org/10.1002/qj.3803

Izumo, T. (2005). The equatorial undercurrent, meridional overturning circulation, and their roles in mass and heat exchanges during El Niño events in the tropical Pacific Ocean. *Ocean Dynamics*, *55*(2), 110–123. https://doi.org/10.1007/s10236-005-0115-1

Izumo, T., Khodri, M., Lengaigne, M., & Suresh, I. (2018). A subsurface Indian Ocean dipole response to tropical volcanic eruptions. *Geophysical Research Letters*, *45*(17), 9150–9159. https://doi.org/10.1029/2018gl078515

Jin, F.-F., Lin, L., Timmermann, A., & Zhao, J. (2007). Ensemble- mean dynamics of the ENSO recharge oscillator under state- dependent stochastic forcing. *Geophysical Research Letters*, *34*(3), L03807. https://doi.org/10.1029/2006GL027372

Johnson, N. C., & Kosaka, Y. (2016). The impact of eastern equatorial Pacific convection on the diversity of boreal winter El Niño teleconnection patterns. *Climate Dynamics*, *47*(12), 3737–3765. https://doi.org/10.1007/s00382-016-3039-1

Johnson, N. C., & Xie, S.-P. (2010). Changes in the sea surface temperature threshold for tropical convection. *Nature Geoscience*, *3*(12), 842–845. https://doi.org/10.1038/ngeo1008

Khodri, M., Izumo, T., Vialard, J., Janicot, S., Cassou, C., Lengaigne, M., et al. (2017). Tropical explosive volcanic eruptions can trigger El Niño by cooling tropical Africa. *Nature Communications*, *8*(1), 1–13. https://doi.org/10.1038/s41467-017-00755-6

Kumar, A., & Hu, Z. Z. (2014). Interannual and interdecadal variability of ocean temperature along the equatorial Pacific in conjunction with ENSO. *Climate Dynamics*, *42*(5), 1243–1258. https://doi.org/10.1007/s00382-013-1721-0

McGregor, S., Ramesh, N., Spence, P., England, M. H., McPhaden, M. J., & Santoso, A. (2013). Meridional movement of wind anomalies during ENSO events and their role in event termination. *Geophysical Research Letters*, *40*(4), 749–754. https://doi.org/10.1002/grl.50136

Okumura, Y. M. (2019). ENSO diversity from an atmospheric perspective. *Current Climate Change Reports*, *5*(3), 245–257. https://doi.org/10.1007/s40641-019-00138-7

Palanisamy, H., Cazenave, A., Delcroix, T., & Meyssignac, B. (2015). Spatial trend patterns in the Pacific ocean sea level during the altimetry era: The contribution of thermocline depth change and internal climate variability. *Ocean Dynamics*, *65*(3), 341–356. https://doi.org/10.1007/s10236-014-0805-7

Rayner, N. A., Parker, D. E., Horton, E. B., Folland, C. K., Alexander, L. V., Rowell, D. P., et al. (2003). Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late 19th century. *Journal of Geophysical Research*, *108*(D14), 4407. https://doi.org/10.1029/2002JD002670

Reynolds, R. W., Rayner, N. A., Smith, T. M., Stokes, D. C., & Wang, W. (2002). An improved in situ and satellite SST analysis for climate. *Journal of Climate*, *15*(13), 1609–1625. https://doi.org/10.1175/1520-0442(2002)015<1609:aiisas>2.0.co;2

Rebert, J. P., Donguy, J. R., Eldin, G., & Wyrtki, K. (1985). Relations between sea level, thermocline depth, heat content, and dynamic height in the tropical Pacific Ocean. *Journal of Geophysical Research*, *90*(C6), 11719–11725. https://doi.org/10.1029/jc090ic06p11719

Williams, I. N., & Patricola, C. M. (2018). Diversity of ENSO events unified by convective threshold Sea Surface temperature: A nonlinear ENSO index. *Geophysical Research Letters*, *45*(17), 9236–9244. https://doi.org/10.1029/2018gl079203

Xie, P., & Arkin, P. A. (1997). Global precipitation: A 17-year monthly analysis based on gauge observations, satellite estimates, and numerical model outputs. *Bulletin of the American Meteorological Society*, *78*(11), 2539–2558. https://doi.org/10.1175/1520-0477(1997)078<2539:gpayma>2.0.co;2

Yokoi, T., Tozuka, T., & Yamagata, T. (2008). Seasonal variation of the Seychelles dome. *Journal of Climate*, *21*(15), 3740–3754. https://doi.org/10.1175/2008jcli1957.1

Zuo, H., Balmaseda, M. A., Tietsche, S., Mogensen, K., & Mayer, M. (2019). The ECMWF operational ensemble reanalysis–analysis system for ocean and sea ice: A description of the system and assessment. *Ocean Science*, *15*(3), 779–808. https://doi.org/10.5194/os-15-779-2019
# Improving and harmonizing El Niño recharge indices

Takeshi Izumo<fnref n="1" />,<fnref n="2" /> and Maxime Colin<fnref n="3" />,<fnref n="4" />,<fnref n="5" />

Submitted to GRL on February 21st, 2022

Revised, August 24th, 2022

Accepted, October 14, 2022

## Contents of this file

Texts S1 to S3; Tables S1 to S3; Figures S1 to S10

## Introduction

The supplementary information contains 3 supplementary texts, 3 tables and 10 figures:

* Supplementary Text S1. *Data and methods*

* Supplementary Text S2. *Discussion on the mechanisms behind $T_E$ and $h_{ind}$ tendency equations*

* Supplementary Text S3. *Discussion on the choice of relative SST*

* Supplementary Table S1. Comparison of the $r_{Te,h_ind}$ skill (using optimal $h_{ind_eq+sw}$) between various data, periods, fields, filtering and ENSO indices.

* Supplementary Table S2. Robustness of the ($T_E$, $h_{ind}$) tendency equation coefficients.

<footnote n="1">¹ Institut de Recherche pour le Développement (IRD), EIO laboratory, Université de la Polynésie Française (UPF), Tahiti, French Polynesia</footnote>

<footnote n="2">² formerly at IRD, Sorbonne Université - CNRS-IRD-MNHN, LOCEAN Laboratory, IPSL, Paris, France</footnote>

<footnote n="3">³Climate Change Research Centre, University of New South Wales, Sydney, New South Wales, Australia</footnote>

<footnote n="4">⁴Laboratoire GEPASUD, University of French Polynesia, Tahiti, French Polynesia</footnote>

<footnote n="5">⁵now at Leibniz Centre for Tropical Marine Research (ZMT), Bremen, Germany</footnote>
* Suppl. Table S3. 14-month lead ENSO hindcasts based on various recharge indices and IOD.

* Suppl. Fig. S1. Respective contributions of the $r_{Te}$ and $r_{h\_ind}$ skills to the $r_{Te,h\_ind}$ skill.

* Suppl. Fig. S2. Equivalent (partly) to Fig. 2 but for ORAS5 reanalysis (1959-2018), comparing SLA, OHC and Z20.

* Suppl. Fig. S3. ENSO asymmetrical Ekman pumping.

* Suppl. Fig. S4. Regression maps of $h(x,y,t)$ or $h_{ind}(x,yt)$ onto $h(t)$ or $h_{ind}(t)$ indices.

* Suppl. Fig. S5. Harmonizing the various recharge indices by using $h_{ind}$ index.

* Suppl. Fig. S6. Comparing ORAS5 SLA and satellite SLA, and comparing SLA, Z20 and OHC in ORAS5, for $h_{ind\_eq+sw}$ index (shown from 1993).

* Suppl. Fig. S7: as Fig. 3, but for satellite SLA (1993-2021).

* Suppl. Fig. S8. As Fig. 3, but for ORAS5 Z20.

* Suppl. Fig. S9. Finding the best averaging box for $h_{ind}$ to optimize both $T_E$ and $h_{ind}$ tendency equations of RO.

* Suppl. Fig. S10. Same as S9ab for ORAS5 SLA and Z20.

## Supplementary Text S1. Data and methods

Here we use classical monthly datasets: Optimum Interpolation SST OISSTv2 based on *in situ* observations and satellite measurements for the recent period (November 1981-Mar2022; Reynolds et al. 2002), HadiSSTv1.1 SST (1870-Mar2022; Rayner et al. 2003) when using longer SLA/Z20/OHC from ECMWF ORAS5 oceanic reanalysis extended version (1959-2018; Zuo et al. 2019; similar results with ORAS5 SST instead of HadiSST; not shown), CMAP1 precipitation (from 1979; Xie and Arkin 1997), and windstress from the latest ECMWF ERA5 atmospheric reanalysis (from 1979 also; Hersbach et al. 2020; using former ERA-Interim reanalysis giving similar results; not shown).

We mostly use SLA as an accurate proxy of thermocline depth/OHC anomalies (e.g. Rebert et al. 1985; Gasparin and Roemmich 2017; Palanisamy et al. 2015), measured from satellites (1993-mid2021; Copernicus product), allowing a better and more homogeneous spatio-temporal sampling than the usual 20°C isotherm depth Z20. We also verified that the results are robust with sensitivity tests shown in Supplementary Tables S1, S2 and S3 and Supplementary Figures S2, S6, S7, S8 and S10. Results are similar at 1st order with SLA, OHC and Z20 in ORAS5. For SLA, we subtract the 60°S-60°N global average at each time step to remove sea level rise global trend due to global warming, and we also remove any additional regional trend through a linear regression (results are similar without this regional detrending).

For the ENSO index, $T_E$, we use the usual Niño3.4 region (170°W-120°W, 5°N-5°S), a reasonable compromise to capture central Pacific as well as eastern Pacific ENSO events (e.g. Takahashi et al. 2011, Capotondi et al. 2015; the present study focusses on typical ENSO events and neglect 2nd order spatial diversity). We define $T_E$ as Niño3.4 relative SST (RSST, i.e. SST minus its 20°N-20°S tropical mean), as recommended by Izumo et al. (2020) and Van Oldenborgh et al. (2021) because atmospheric tropical deep convection interannual anomalies are rather related to RSST than to SST, notably in the
presence of external forcing (see also Johnson and Kosaka 2016; Khodri et al. 2017; Izumo, Khodri et al. 2018; Williams and Patricola 2018, and Okumura 2019). The deep convection threshold SST>~27-28°C (e.g. Gadgil et al. 1984) translates into RSST>~1°C, a threshold that remains valid with global warming (e.g. Johnson and Xie 2010). See Supplementary Text S3 for an extensive discussion justifying the relevance of RSST for RO equations.

$T_E$ is normalized (divided by its standard deviation; STD), so that a 1-value represents typical ENSO amplitude. $h$ is normalized when a regional average is done (e.g. usual 5°N-5°S, 120°E-80°W for $h_{eq}$ and 5°N-5°S, 120°E-155°W for $h_w$, and for $h_{ind\_eq+sw}$), but not when taken at each (x,y) point, so that related regression maps represent typical anomalous amplitudes (cm for SLA, m for Z20). Using SST instead of RSST, or Niño3 instead of Niño3.4 region, makes $T_E$ slightly less correlated to equatorial Pacific $\tau_x$ (i.e. ocean-atmosphere coupling), confirming that Niño3.4 RSST is a relevant choice. They anyway lead to very similar results (Suppl. Text S3), with slightly weaker skills when using Niño3 (Suppl. Table S1).

Here we use typical statistical methods. The monthly seasonal cycle (computed by averaging each month of the year over the whole period) is removed and intraseasonal variations are filtered out by a 4-month Hanning filter (3-point Hanning filter in pyferret software for monthly time series), so that periods lower than ~2-3 months are removed (considering a fast adjustment to be with timescales less than ~2-3 months). For ORAS5 analysis over its whole extended period 1959-2018, we use in addition to the intraseasonnal low-pass filter a highpass Hanning filter (14-year window; i.e. 7 years cut on each side) to remove interdecadal variability with cutting frequency at periods of ~10-year (results are robust without this highpass filter, with just some increased interdecadal noise). Coefficients obtained from multivariate linear regression fits are shown multiplied by the STD of their multiplying variable ($h$ or $h_{ind}$, cf. section 3; STD($T_E$) is already 1), so as to represent typical amplitudes.

For statistical significance, we use Steigers Z-method for difference between two dependent correlations from a single sample (Hotelling William method giving similar results; using https://www.quantitativeskills.com/sisa/statistics/correl.htm) and typical two-tailed Student t-tests for 90% confidence intervals. Using formulae (30) of Bretherton et al. (1999), we have about one effective degree of freedom every 4 months, i.e. ~85 for 1993-mid2021 satellite SLA period (~140 for filtered ORAS5 data), a sufficiently large number leading usually to strong statistical significances of the results.

Here we explain how to obtain $h_{ind\_eq+sw}$ in details (and give the $K$ dependance to data chosen for $h$) : 1) average over eq+sw box (equatorial box 5°N-5°S, 120°E-80°W + southwest box 5°S-15°S, 120°E-170°W) the SLA (using $h(x,y,t) \sim$ SLA(x,y,t), detrended, smoothed, with seasonal cycle removed; cf. start of this text S1) and then normalize, 2) remove its dependent part $KT_E: h_{ind\_eq+sw} = h_{eq+sw} - K_{eq+sw} T_E$, $T_E$ being Niño3.4 relative SST anomaly (smoothed and normalized) with regression coefficient $K_{eq+sw} \approx 0.26[0.30/0.21]$ for satellite SLA[ORAS5_SLA/Z20] ($K_{eq+sw}$ is conveniently already small thanks to sw addition to eq region, conversely to usual eq and w regions, for which $K_{eq}=+0.70[+0.68/+0.38]$ and $K_w=-0.40[-0.36/-0.51]$ respectively).
**Supplementary Text S2. Discussion on the mechanisms behind $T_E$ and $h_{ind}$ tendency equations**

Here we explain physically and mathematically the spatial patterns of Fig. 2’s various panels, very robust among datasets and periods (cf. Suppl. Fig. S2). And why there are similarities among some of them. Along these explanations we discuss the mechanisms behind the various tendency equations terms (descriptions done here by convention in the recharge case related to previous La Niña conditions; vice versa in the discharge case related to previous El Niño conditions).

$F_2$ map (Fig. 2h; and its equivalents in Suppl. Fig. S2) shows us how $h_{ind}(x,y,t)$ would look like if ENSO windstress anomalies would blow for a long time, e.g. after long-lasting La Niña conditions. This physical meaning of $F_2$ is mathematically because the temporal integral of $dh_{ind}/dt$ equation (6) corresponds to the regression of $h_{ind}(x,y,t)$ onto the temporal integral of $T_E$, as $\varepsilon$ is negligible. The slow recharge in the western and central equatorial Pacific shown by $F_2$ map is consistent with former EOF analyses (e.g. Meinen and McPhaden 2000, Alory and Delcroix 2002, Clarke et al. 2007, Clarke 2010, Kumar and Hu 2014). As mentioned in the main manuscript, this western and equatorial recharge is through: 1) downwelling equatorial Rossby waves to the west (Wyrtki 1985, Jin 1997a, 1997b), and off-equatorial ones in the southwest; 2) upwelling equatorial Kelvin waves to the east forcing coastal Kelvin waves propagating poleward along the eastern boundary and thus a leakage of negative $h$ anomalies towards the poles along the eastern boundary (Wyrtki 1985, Izumo, Lengaigne et al. 2018; vice versa for El Niño case). The recharge has also been suggested to be through Sverdrup transport towards the equator (Jin 1997a,b) but this is now debated (Clarke 2010; Zhu et al. 2018; Izumo, Lengaigne et al. 2018). Whatever the mechanisms, they are all formally well represented at 1st order by the term $-F_2T_E$.

Concerning the recharge meridional asymmetry with a larger recharge in the southwest shown by $F_2$ map, it is caused by Ekman pumping asymmetry, but not by the western boundary coastline meridional asymmetry, as mentioned in the main manuscript. This asymmetrical Ekman pumping is mainly due south of ~7°S to windstress curl related to a poleward shift of the SPCZ (Suppl. Fig. S3). It forces locally downwelling (in the La Niña case) and thus slow off-equatorial downwelling Rossby waves progressively recharging the southwestern Pacific (see also Alory and Delcroix 2002; Cibot et al. 2005).

$F_1$ physically represents the slow recharge mode influence on $T_E$. $F_1$ tends to have a similar pattern to $F_2$, statistically because $h_{ind}$ tends to co-vary in the west and southwest Pacific (as the spatial pattern of ENSO-related windstress remains at first order similar among ENSO events). Physically, the thermocline depth deepening related to a recharge (La Niña case) favors positive $T_E$ through the thermocline feedback in the central Pacific (where coupling with convection is large; Clarke et al. 2007) and eastern Pacific (Wyrtki 1985, Jin 1997ab). Southwest and equatorial thermocline deepening can also increase $T_E$ by increasing the transport-weighted temperatures $T_{conv}$ and $T_{EUC}$ of STCs (Shallow Subtropical/Tropical meridional overturning Cells) lower branches. These are the meridional pycnocline convergence $T_{conv}$ (knowing that pycnocline convergence is larger in the southwest than in the northwest) and the equatorial undercurrent (EUC) $T_{EUC}$ (e.g. Izumo 2005, Ballester et al. 2015, 2016ab). These influences could be formalized as two terms, $F_{conv}T_{conv}$ and $F_{EUC}T_{EUC}$, in the $dT_E/dt$ equation (where $F_{conv}$ and $F_{EUC}$ are constant
<table><thead><tr><th>r<sub>Te,h_ind</sub> skill<br>Comparing<br>satellite/ORAS5</th><th>Sat. SLA<br>(1993-2021)</th><th>Sat. SLA<br>(1993-2018)</th><th>ORAS5 SLA<br>(1993-2018)</th><th>ORAS5 Z20<br>(1993-2018)</th><th>ORAS5 OHC<br>(1993-2018)</th></tr></thead><tbody><tr><td>from 1993</td><td>0.69</td><td>0.68</td><td>0.66</td><td>0.71</td><td>0.65</td></tr><tr><td>Comparing<br>SLA/Z20/OHC<br>in ORAS5</td><td>SLAfilt</td><td>Z20filt</td><td>SLA</td><td>Z20</td><td>OHC</td></tr><tr><td>from 1959</td><td>0.66</td><td>0.72</td><td>0.63</td><td>0.65</td><td>0.60</td></tr><tr><td>Testing ENSO<br>indices</td><td>Sat. SLA<br>(1993-2021)</td><td>ORAS5 SLA</td><td>ORAS5 Z20</td><td>ORAS5<br>SLAfilt</td><td>ORAS5<br>Z20filt</td></tr><tr><td>Niño3.4 RSST</td><td>0.69</td><td>0.64</td><td>0.67</td><td>0.66</td><td>0.72</td></tr><tr><td>Niño3.4 SST</td><td>0.70</td><td>0.65</td><td>0.66</td><td>0.68</td><td>0.72</td></tr><tr><td>Niño3 RSST</td><td>0.65</td><td>0.62</td><td>0.66</td><td>0.66</td><td>0.72</td></tr><tr><td>Niño3 SST</td><td>0.66</td><td>0.63</td><td>0.65</td><td>0.64</td><td>0.72</td></tr></tbody></table>

Suppl. Table S1. Comparison of the *r*<sub>Te,h_ind</sub> skill (using optimal *h*<sub>ind_eq+sw</sub>) between different data, periods, fields, filtering and ENSO indices. *Upper rows:* comparison between satellite and ORAS5 on the same period, from 1993 (also verifying the unsensitivity to the change of the final year from 2021 to 2018). *Central rows:* skill stable over the long period 1959-2018 using ORAS5, with interdecadal variations filtered out ('filt') or not. *Bottom rows:* verifying equivalence between ENSO indices. We have also verified the non-sensitivity to SST datasets: e.g. ORAS5 SST leads to almost no change (here we use OISST for recent periods and HadiSST for long periods from 1959).
<table><thead><tr><th colspan="5">Comparison for satellite SLA (1993-2021)</th></tr><tr><th></th><th>R<sub>o</sub></th><th>F<sub>1</sub></th><th>F<sub>2</sub><sub>o</sub></th><th>ε<sub>o</sub></th></tr></thead><tbody><tr><td>h<sub>w</sub></td><td>0.06±0.04</td><td>0.15±0.04</td><td>0.15±0.04</td><td>0.06±0.04</td></tr><tr><td>h<sub>eq</sub></td><td>-0.18±0.05</td><td>0.25±0.05</td><td>0.25±0.05</td><td>-0.18±0.05</td></tr><tr><td>h<sub>eq+sw</sub></td><td>-0.05±0.03</td><td>0.18±0.03</td><td>0.18±0.04</td><td>-0.05±0.04</td></tr><tr><td></td><th>R</th><th>F<sub>1</sub></th><th>F<sub>2</sub></th><th>ε</th></tr><tr><td>h<sub>ind w</sub></td><td>0.00±0.04</td><td>0.13±0.04</td><td>0.13±0.04</td><td>0.00±0.04</td></tr><tr><td>h<sub>ind eq</sub></td><td>0.00±0.04</td><td>0.16±0.04</td><td>0.17±0.04</td><td>0.00±0.04</td></tr><tr><td>h<sub>ind eq+sw</sub></td><td>0.00±0.04</td><td>0.18±0.04</td><td>0.18±0.04</td><td>0.00±0.04</td></tr><tr><td></td><th>R</th><th>F<sub>1</sub></th><th>F<sub>2</sub></th><th>ε</th></tr><tr><td>ORAS5 SLA (1993-2018)</td><td>0.00±0.02</td><td>0.17±0.02</td><td>0.17±0.02</td><td>0.00±0.02</td></tr><tr><td>ORAS5 SLAfilt (1959-2018)</td><td>0.00±0.03</td><td>0.16±0.03</td><td>0.16±0.03</td><td>0.00±0.02</td></tr><tr><td>ORAS5 Z20filt (1959-2018)</td><td>0.00±0.03</td><td>0.16±0.03</td><td>0.16±0.02</td><td>0.00±0.02</td></tr></tbody></table>

**Suppl. Table S2. Robustness of the ($T_E$, $h_{ind}$) tendency equation coefficients.** The sensitivity tests here show that, when using independent $h_{ind}$ instead of usual $h$, values of the regression coefficients become stable for different recharge indices, data, fields (SLA and Z20) and periods: $F_1 \sim F_2 \sim 0.16$ to $0.18 \pm 0.03$, $R \sim \varepsilon \sim 0.0 \pm 0.03$ (for normalized indices, i.e. unit in month$^{-1}$; 90% confidence intervals shown). Conversely, for usual $h$-based indices, the coefficients are not stable (cf. section 3).

Note that the RO, even if theoretically undamped as $R$ and $\varepsilon$ are negligible, is actually still damped because residual$_T$ (equation 6) is not a pure red noise but includes non-linear terms neglected in our 1st order linear approximation. Indeed, if we add to the $T_E$ tendency equation a simple non-linear term $cT_E^2$ (quadratic term notably related to the non-linear response of convection to $T_E$; e.g. Jin et al. 2020; see An et al. 2020 review on more complicate possibilities, e.g. a multiplicative noise, see also e.g. Jin et al. 2007, Graham et al. 2015, their equation (23)), the fit gives $c=+0.05 \pm 0.03$ and $R= -0.01 \pm 0.03$ month$^{-1}$. $R$ becomes weakly negative (not at the 90% level but robust for longer ORAS5 SLA and Z20; not shown), suggesting that the oscillator is weakly damped.
<table><thead><tr><th></th><th>r<sub>Te,h ind</sub> skill</th><th>a (IOD)</th><th>b (h<sub>ind</sub>)</th><th>c (T<sub>E</sub>)</th></tr></thead><tbody><tr><td colspan="5"><strong>For satellite SLA (1993-2022):</strong></td></tr><tr><td><strong>h<sub>ind eq</sub></strong></td><td><strong>0.50</strong></td><td><strong>-0.55±0.36</strong></td><td><strong>0.24±0.28</strong></td><td><strong>0.30±0.36</strong></td></tr><tr><td><strong>h<sub>ind eq+sw</sub></strong></td><td><strong>0.58</strong></td><td><strong>-0.48±0.34</strong></td><td><strong>0.39±0.26</strong></td><td><strong>0.26±0.34</strong></td></tr><tr><td colspan="5"><strong>For ORAS5 SLAfilt (1966-2011):</strong></td></tr><tr><td><strong>h<sub>ind eq</sub></strong></td><td><strong>0.54</strong></td><td><strong>-0.50±0.28</strong></td><td><strong>0.32±0.21</strong></td><td><strong>0.10±0.28</strong></td></tr><tr><td><strong>h<sub>ind eq+sw</sub></strong></td><td><strong>0.60</strong></td><td><strong>-0.43±0.27</strong></td><td><strong>0.42±0.20</strong></td><td><strong>0.05±0.26</strong></td></tr><tr><td colspan="5"><strong>For ORAS5 Z20filt (1966-2011):</strong></td></tr><tr><td><strong>h<sub>ind eq</sub></strong></td><td><strong>0.65</strong></td><td><strong>-0.43±0.25</strong></td><td><strong>0.48±0.19</strong></td><td><strong>0.05±0.25</strong></td></tr><tr><td><strong>h<sub>ind eq+sw</sub></strong></td><td><strong>0.68</strong></td><td><strong>-0.40±0.25</strong></td><td><strong>0.52±0.18</strong></td><td><strong>0.03±0.24</strong></td></tr></tbody></table>

Suppl. Table S3. 14-month lead ENSO hindcasts based on IOD (Indian Ocean Dipole) and various recharge indices. 14-month lead hindcasts of $T_E$ in NDJyr1 (November of year 1 to January of year 2) using a IOD + b $h_{ind}$ + c $T_E$ in $SON_{yr0}$ (September to November of year 0), with coefficients a, b, c of the multivariate linear regression estimated from least-square fit. The 90% confidence intervals are shown, with the regression coefficients in bold when significantly different from zero at the 90% level.

Main messages of Table S3 are:

1) adding the SouthWest box (5°S-15°S, 120°E-170°W) improves the skill for all datasets/periods, more clearly for SLA than for Z20.

2) ORAS5 Z20 skill is better than ORAS5 SLA one, suggesting that thermocline depth metrics could perform better in a perfectly-observed ocean. Yet SLA has the advantage of being observed by satellite in near real-time, and the additional advantage of being more available in climate models outputs (e.g. CMIP) than Z20.

3) we have statistically-significant contributions from $h_{ind\_eq+sw}$ and IOD, but not from $T_E$ itself (the IOD contribution a is always significantly negative, while the coefficient c for ENSO itself is never significant at 90% level). Hence this improved recharge index with updated datasets confirms earlier studies emphasizing the IOD influence on following year's ENSO phase (Izumo et al. 2010, 2014, 2016, Dayan et al. 2014, Jourdain et al. 2016). Regression coefficients have qualitatively similar values among datasets/periods.
a) skill of the regression of $dT_E/dt$ onto ($T_e, h_{ind}$) : $r_{Te}$

b) skill of the regression of $dh_{ind}/dt$ onto ($T_e, h_{ind}$) : $r_{h_{ind}}$

<figure><img src="image_7.png" /><figcaption>Suppl. Fig. S1. Respective contributions of the $r_{Te}$ skill and $r_{h_{ind}}$ skill to the $r_{Te,h_{ind}}$ skill. a) Pearson correlation skill $r_{Te}$ for the regression of $dT_E/dt$ onto ($T_e, h_{ind}(x,y,t)$) (equation 4). b) same but for $r_{h_{ind}}$ for the regression of $dh_{ind}/dt$.</figcaption></figure>

**Main messages of Fig. S1:**

Both $r_{Te}$ and $r_{h_{ind}}$ contribute to the skill $r_{Te,h_{ind}}$ (square-root of the RV coefficient) shown in Fig. 2k for the tendency equation of the vector ($T_e, h_{ind}$), since we can mathematically show that $r_{Te,h_{ind}}^2 \approx (r_{Te}^2 + r_{h_{ind}}^2)/2$ for normalized vectors, as $T_E$ and $h_{ind}$ are by definition uncorrelated, and as their tendencies are also almost uncorrelated (<0.1). Note that $r_{Te}$ tends to have a similar spatial pattern to $r_{h_{ind}}$, but of stronger amplitude (though with also relatively larger central-equatorial values). The similarity of the spatial patterns, notably with the strong skill in the west and southwest Pacific, is likely because $h_{ind}$ tends to vary in phase in the west-central equatorial Pacific and southwest tropical Pacific (cf. Supplementary Text S2). $r_{Te}$ tends to be of larger amplitude than $r_{h_{ind}}$, likely because $dT_E(t)/dt$ is less noisy ($T_E$ being a box average) than $dh_{ind}(x,y,t)/dt$, which is for each (x,y) point. The pair ($T_E(t), h_{ind}(x,y,t)$) can thus more easily explain statistically $dT_E(t)/dt$ variance than $dh_{ind}/dt$ variance.
<figure><img src="image_8.png" /></figure>

**Suppl. Fig. S2. Equivalent (partly) to Fig. 2 (which is based on 1993-2021 satellite SLA) but for ORAS5 (60-yr long period 1959-2018), comparing SLA (left), OHC (middle) and Z20 (right).**
First upper-left panel is the skill map, as Fig. 2k (with the same color scale as Fig. 2k). 2^nd^ upper panel is the regression of *h* onto *T*<sub>E</sub> (the fast adjustment mode), as Fig. 1a. Panels of 2^nd^ and 3^rd^ row panels are for usual non-orthogonal (*T*<sub>E</sub>,*h*<sub>reg</sub>) basis, as Fig. 2a,b,e,f. Panels of 4^th^ and 5^th^ rows are for orthogonal (*T*<sub>E</sub>,*h*<sub>ind</sub>) basis, as Fig. 2c,d,g,h. SLA is as in Fig. 2 in cm. OHC is the oceanic heat content from ocean surface to bottom in J (color scale from -1.8*10<sup>9</sup> to +1.8*10<sup>9</sup> J in 1^st^ row, and from -1.6*10<sup>8</sup> to +1.6*10<sup>8</sup> J in 3^rd^ and 5^th^ rows). Z20 is the 20°C isotherm depth in m (no interdecadal high-pass filter applied here; only weak differences if applied; not shown).

**Main messages of Fig. S2 are:**

1) ORAS5 SLA gives similar results to satellite SLA (we have also tested the sensitivity to the periods chosen: ORAS5-based maps are similar for 1959-2018 and shorter 1993-2018 period; not shown).

2) SLA, OHC and Z20 are qualitatively similar, with some differences in patterns and in the correlation between $T_E$ and $h$. OHC spatial patterns tend to be in between SLA and Z20 patterns, but the skill is weaker. Note that the lag-0 correlation along the central-eastern equatorial Pacific (i.e. for $h_{eq}$ and $h_{eq+sw}$) with $T_E$ is weaker for Z20 than for SLA, possibly because of 1st baroclinic mode playing more role for Z20 than for SLA, which will be more influenced by higher modes, notably the 2nd baroclinic mode (we have compared each baroclinic mode contribution in the LCS simulation forced by ERA-I windstress of Izumo, Lengaigne et al. 2018; not shown). These different baroclinic mode contributions could partly explain some differences in results between Izumo, Lengaigne et al. (2018) study, based on the LCS model that take into account several baroclinic modes, and Neske and McGregor (2018) study based on a shallow water model, which only simulates the 1st baroclinic mode.
<figure><img src="image_9.png" /></figure>

**Suppl. Fig. S3.** ENSO asymmetrical Ekman pumping. **a)** Zonal wind stress $\tau_x$ (color, N.m$^{-2}$), precipitation (contours, mm.day$^{-1}$) and wind stress $\tau$ (vectors) regressed on $T_E$. **b)** Regression of the 150°E-140°W zonal mean of Ekman pumping (black, 10$^{-5}$ m/s) on $T_E$, and Ekman pumping decomposition into its wind stress curl (red) and beta (green) terms, only defined out of the 5°N-5°S equatorial band (shown by vertical lines).

**Main messages of Fig. S3:**

In the southern hemisphere, the decomposition shows that south of ~7°S, the wind stress curl term is the main contributor to Ekman pumping, while north of ~7°S, the beta term is also important. Hence, Ekman pumping asymmetry is mainly due to wind stress curl (curl itself mostly due to $d\tau_x/dy$; not shown), with some added asymmetry from the beta term, proportional to $\tau_x$ (Yokoi et al. 2008), large between ~7°S and 5°S.

Note that the western boundary coastline meridional asymmetry would conversely favor a larger northwest recharge (in the La Niña case discussed here below; discharge in the El Niño case), as shown by an LCS (Linear Continuously Stratified model; McCreary 1980) idealized experiment with a symmetric zonal wind stress forcing (Izumo, Lengaigne et al. 2018, their Fig. 5). We have done a similar experiment, but with a more realistic ENSO-like asymmetric windstress (like the one of panel a), with easterlies shifted to the south during La Niña, plus a reversal to westerlies south of ~15°S, related to a poleward shift of the SPCZ (vice versa for the El Niño case; see also Alory and Delcroix 2002; Cibot et al. 2005; McGregor et al. 2013). This asymmetry indeed favors a larger southwest recharge (not shown).
<figure><img src="image_10.png" /></figure>

Suppl. Fig. S4. Regression maps of *h(x,y,t)* or *h<sub>ind</sub>(x,y,t)* onto various *h(t)* or *h<sub>ind</sub>(t)* recharge indices. Spatial patterns associated with *h*<sub>eq</sub> (left column), *h*<sub>w</sub> (middle column) and *h*<sub>eq+sw</sub> (right column) look significantly closer for *h*<sub>ind</sub> (2<sup>nd</sup> and 4<sup>th</sup> rows) than for usual full *h* (1<sup>st</sup> and 3<sup>rd</sup> rows). I.e. *h*<sub>ind</sub> harmonizes recharge metrics, making indices converge. First two rows are based on satellite SLA, the last two on ORAS5 Z20filt. Recharge indices are normalized, so that these maps show typical amplitudes (units: cm for SLA and m for Z20).
<figure><img src="image_11.png" /><figcaption>Suppl. Fig. S5. Using $h_{ind}$ harmonizes the various recharge indices. a) time series of $h_{eq}$ (green), $h_{w}$ (red) and $h_{eq+sw}$ (black) from satellite SLA. b) as a but for their independent components $h_{ind\_eq}$ (green), $h_{ind\_w}$ (red) and $h_{ind\_eq+sw}$ (black).</figcaption></figure>

Main message of Fig. S5:

Indices are much better correlated for independent component than for full signal:

$r(h_{ind\_eq}, h_{ind\_w}) = 0.79$ ($r^2=62\%$) instead of 0.18 (3%)

$r(h_{ind\_eq}, h_{ind\_eq+sw}) = 0.91$ (82%) instead of 0.79 (62%)

$r(h_{ind\_w}, h_{ind\_eq+sw}) = 0.88$ (78%) instead of 0.67 (44%)

This confirms that using $h_{ind}$, i.e. removing the dependent component, helps in reconciliating and harmonizing the various recharge indices used in the literature. Furthermore, $h_{ind\_eq+sw}$ is close both to $h_{ind\_eq}$ and $h_{ind\_w}$.
a) Comparison between sat. SLA and ORAS5 SLA for $h_{ind\_eq+sw}$

<figure><img src="image_12.png" /></figure>

b) Comparison in ORAS5 between SLA, Z20 and OHC

<figure><img src="image_13.png" /></figure>

Suppl. Fig. S6. Comparing a) ORAS5 SLA (red) with satellite SLA (black), and b) SLA (red), Z20 (green) and OHC (blue) in ORAS5, for $h_{ind\_eq+sw}$ index (shown from 1993). Panel a validates ORAS5 SLA, highly correlated (r=0.97; r²=95%) to observed satellite SLA for $h_{ind\_eq+sw}$ (and also for $h_{ind\_eq}$ and $h_{ind\_w}$; not shown; ORAS5 shown here without high-pass decadal filter here for comparison with satellite SLA; note that ORAS5 may not be as good before 1993 before which the reanalysis does not assimilate satellite sea level observations). Panel b shows that, in ORAS5, SLA is highly correlated to Z20 (r=0.94; r²=88%) and OHC (r=0.96; r²=92%), i.e. SLA, Z20 and OHC are almost equivalent for $h_{ind\_eq+sw}$ (and also for $h_{ind\_eq}$ and $h_{ind\_w}$; not shown).

Note that using $h_{ind}$ instead of usual $h$ also harmonizes SLA-based and Z20-based metrics: e.g. the squared correlation between SLA-based $h_{eq}$ and Z20-based $h_{eq}$ is r²=81% (88%) instead of 72% (81%). Adding the Southwest further increases r² to 85% (90%) for full $h_{eq+sw}$ and even to 88% (92%) for independent $h_{ind\_eq+sw}$. The r² value is given over the 1993-2019 period (in parenthesis for whole ORAS5 1959-2018 period with interdecadal high-pass filter applied).
<figure><img src="image_14.png" /><figcaption>Suppl. Fig. S7: as Fig. 3, but for satellite SLA over its shorter available period (1993-2021).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Suppl. Fig. S8. As Fig. 3, but for ORAS5 Z20filt. Results are robust, with only 2nd order differences: e.g. extreme El Niño events cause a weaker discharge in Z20 than in SLA (and thus a weaker El Niño/La Niña asymmetry). Hence, trajectories in bottom panels look even more circular and closer to the ideal RO circular trajectories when using Z20 instead of SLA (and when removing interdecadal variability, as done here).</figcaption></figure>
<figure><img src="image_16.png" /></figure>

Suppl. Fig. S9. Finding the best averaging box for *h*<sub>ind</sub> to optimize both *T*<sub>E</sub> and *h*<sub>ind</sub> tendency equations of RO. Panel a is almost similar to Fig. 2j, showing the optimization of the eastern edge longitude. The black line shows the skill *r*<sub>Te,h ind</sub> for *h*<sub>ind</sub> averaged over: the equatorial band (5°N-5°S) with its western edge fixed to 120°E and its eastern edge varying, given by the x-axis. The black line shows that the optimal eastern edge for the equatorial band is at equatorial Pacific eastern boundary, 80°W. Red line is for a two-rectangle region, this classical equatorial box (5°N-5°S, 120°E-80°W) plus a southern 5°S-15°S box with the same western edge fixed at 120°E and its eastern edge varying. Its optimal is around 170°W (indicated by a dashed black vertical line, i.e. the optimal eq+sw region). The green and blue lines are similar to the red line, but specifically for *r*<sub>Te</sub> and for *r*<sub>h ind</sub> respectively. Panel c shows the associated lag-0 correlation between *T*<sub>E</sub> and *h*, for the same regions as the black and red lines in panel a. The lag-0 correlation interestingly is almost the lowest for the optimal eq+sw region. This is an additional interest of this region, e.g. the difference between *h*<sub>ind eq+sw</sub> and *h*<sub>eq+sw</sub> is much weaker than between *h*<sub>ind eq</sub> and *h*<sub>eq</sub>. Concerning the western edge longitude, the optimal is around 100°E-120°E, so we have kept 120°E to stay in Pacific, as Meinen and McPhaden (2000) choice. Right panels are showing the optimization of the latitudinal edges of the *h*<sub>ind</sub> two-rectangle averaging region. The black curve is for the southern edge of the southwest rectangle (120°E-170°W, 5°S to y<sub>sw</sub> indicated in horizontal axis, the other region being fixed to the classical equatorial box 5°N-5°S, 120°E-80°W). y<sub>sw</sub> optimum is around 16-12°S (cf. black dashed vertical line at 14°S), confirming that 15°S is a good simple choice. The green dashed curve is for the southern edge of the equatorial band rectangle (120°E-80°W; 5°N-
$y_{eqS}$, the southwest region being 120°E-170°W, 15°S to $y_{eqS}$). The optimum is around 6-2°S, confirming that 5°S is a good simple choice. The red curve is for equatorial band northern edge (120°E-80°W; 5°S-$y_{eqN}$, the southwest region being 120°E-170°W, 15°S to 5°S). The optimum is around 2-6°N, confirming that 5°N is a good simple choice. Note that theoretically, we could even do a multivariate optimization. We have done some tests that confirm our choice (e.g. if $y_{sw}=14°S$, $y_{eqS}=4°S$ and $y_{eqN}=3°N$, skill is 0.70, only slightly larger than 0.69 for 15°S, 5°S, 5°N, both skills being much higher than 0.61 for classical $h_{eq}$), within the range of uncertainties (e.g. differences between satellite SLA shown here and ORAS5 Z20filt in next figure). All these analyses confirm that the optimal (and still simple) region combines the classical equatorial band (5°S-5°N, 120°E-80°W) and the southern band 15°S-5°S, from 120°E to about 170°W, i.e. $h_{ind\_eq+sw}$.

<figure><img src="image_17.png" /><figcaption>Suppl. Fig. S10. Same as S9ab (based on satellite SLA) but using ORAS5 SLAfilt (upper) or ORAS5 Z20filt (lower). Upper panels show that ORAS5 SLA results are consistent with satellite ones, with similar optimal longitude and latitudes (skill overall slightly weaker). Lower panels show that Z20 is overall consistent with SLA, with however less sharp latitudinal edges and an optimal eastern edge of southwest rectangle at ~155°W instead of ~170°W. But the difference between these two longitude choices is *in fine* not significant anyway (almost similar skill, and correlation of 0.99 between the two definitions). Therefore, SLA and Z20 are almost similar for $h_{ind\_eq+sw}$ (correlation of 0.96, cf. Fig. S6b).</figcaption></figure>