Causal dependences between the coupled ocean–atmosphere dynamics over the tropical Pacific, the North Pacific and the North Atlantic

The causal dependences (in a dynamical sense) between the dynamics of three different coupled ocean–atmosphere basins, the North Atlantic, the North Pacific and the tropical Pacific region (Nino3.4), have been explored using data from three reanalysis datasets, namely ORA-20C, ORAS4 and ERA-20C. The approach is based on convergent cross mapping (CCM) developed by Sugihara et al. (2012) that allows for evaluating the dependences between variables beyond the classical teleconnection patterns based on correlations.

The use of CCM on these data mostly reveals that (i) the tropical Pacific (Nino3.4 region) only influences the dynamics of the North Atlantic region through its annual climatological cycle; (ii) the atmosphere over the North Pacific is dynamically forcing the North Atlantic on a monthly basis; (iii) on longer timescales (interannual), the dynamics of the North Pacific and the North Atlantic are influencing each other through the ocean dynamics, suggesting a connection through the thermohaline circulation.

These findings shed a new light on the coupling between these three different regions of the globe. In particular, they call for a deep reassessment of the way teleconnections are interpreted and for a more rigorous way to evaluate dynamical dependences between the different components of the climate system.

1 Introduction

In environmental sciences, statistical quantities are essential tools to characterize the properties of a system, the most familiar of which are the mean, the variance and the correlation in space or time. Correlations are very often associated with the notion of dependences in climate sciences, assuming that a certain variable is influencing the one to which it is correlated. Although it is true when dealing with a Gaussian linear system, the picture becomes far more complicated when dealing with a nonlinear system. In particular, for nonlinear deterministic dynamical systems, the correlation between variables is neither sufficient nor necessary for dependence between these variables, e.g., Granger (1969, 2004) and Sugihara et al. (2012).

Teleconnections in the form of correlations between distant points in space are used in climate sciences in order to evaluate the link of a dynamical process in some part of the world with a distant target. In particular, an important question that has attracted a lot of attention in the past decades is to know whether the tropical Pacific system forces the dynamics of the climate system in the extratropics. In this context, important teleconnections are found between events like El Niño or La Niña and the temperature and precipitation patterns all over the world, e.g., Fraedrich and Müller (1992), Brönnimann (2007) and Lau (2016). The origin of these teleconnections for the North Atlantic and North Pacific is currently explained through the concept of the atmospheric bridge that allows for the transfer of information from one basin to another, e.g., Sardeshmukh and Hoskins (1988), Alexander et al. (2002), Yu and Lin (2016) and Lau (2016). This explanation assumes that there is a causality principle leading to these teleconnections, mostly going from the tropical Pacific to the remote regions. This view originating from teleconnection patterns should however be taken with care, since co-variability does not imply interdependences as already mentioned above. Another possible explanation of these teleconnections is the influence of an external driver on both variables that are correlated, even if they are dynamically independent.

How can we then measure this dependence? Answering this question is difficult, as discussed by Clive Granger in his Nobel lecture in 2003 (Granger, 2004). One way proposed by Granger (1969) is to use the information on the predictability of the system with or without the influence of the variable expected to be the cause. In this context, two forecasting models should be developed, one with and the other without the variable investigated as predictor (Mokhov et al., 2011; Mosedale et al., 2006; Tirabassi et al., 2015). A drawback of the approach is precisely the necessity to build such a forecasting model. Moreover, as discussed in detail in the Supplement of Sugihara et al. (2012), the approach can lead to ambiguous results when applied to nonlinear deterministic dynamical systems.

A powerful method has been recently proposed by Sugihara et al. (2012), known as convergent cross mapping (CCM), which is a method suitable for nonlinear deterministic dynamical systems as it is based on analogs of the current state. This method has also been tested with success when nonlinear dynamical systems are affected by noise (Mønster et al., 2017) and in coupled dynamical systems in order to identify the leading element of the coupling (BozorgMagham et al., 2015). It has also been recently used to disentangle the link between galactic cosmic rays and the variations of the global temperature (Luo et al., 2015; Tsonis et al., 2015; Ye et al., 2015) or between environmental drivers and influenza (Deyle et al., 2016). This is the method that will be used in the present study. Alternative approaches based on the transfer of information are very appealing and a lot of progress have been made in that direction (Liang, 2014, 2015; Liang and Kleeman, 2005; Runge et al., 2012). We however do not pursue that direction and leave the use of these techniques for a follow-up study.

In the present work, we address the question of causality dependence (in a dynamical sense) between the tropical Pacific, the North Atlantic and the North Pacific coupled ocean–atmosphere dynamics at monthly to interannual timescales, in order to clarify the remote role of these different climate subsystems on the others. This will be done by first constructing low-order systems based on projections on a few Fourier modes that are assumed to dominate the dynamics in each of the basins. This projection has already been applied successfully in the context of the analysis of the coupling between the ocean and the atmosphere over the North Atlantic (Vannitsem and Ghil, 2017). Once these projections are identified, the time series associated with each of these modes can be analyzed using the CCM approach. The specific choice of regions and Fourier modes, although reasonable, is however a little arbitrary. The current analysis is therefore a proof of concept on investigating the question of dependences in the Earth system dynamics using the CCM approach, keeping in mind that an extension of the analysis to other regions and state (or phase) space representations should be explored.

Section 2 will introduce the technique of CCM. In Sect. 3, the datasets and the projections used are described. The datasets are coming from reanalyses performed at the European Centre for Medium-Range Weather Forecasts (ECMWF). The results of the application of CCM on these data are then presented in Sect. 4. The main conclusions and future works are outlined in Sect. 5.

2 Convergent cross mapping

In the present work, as in Sugihara et al. (2012), two variables recorded as a function of time, say X(t) and Y(t), are said causally linked if they are coming from the same dynamical system. In this case, a two-way dependence relation is present between them and the information gathered from one of the variables should, in principle, provide information on the other one. This type of dependence should not be confused with the case where a dynamical system is forced by an external driver, in which case there is a one-way dependence from the driver to the dynamical system. Disentangling the nature of this coupling is crucial in science when one is interested in describing the dynamics of the system under investigation. This is particularly true when the system is very complicated, as the Earth system is.

Several approaches have been developed in the recent past that allow for analyzing the dependences between time series. Granger causality (GC) analysis is a celebrated approach based on the evaluation of the predictability of a variable in the absence or the presence of an hypothetical driver (Granger, 1969). As indicated in Sugihara et al. (2012), its application is restricted to separable systems for which the driver can be effectively removed. In nonlinear deterministic dynamical systems in which all the variables are interconnected, the GC approach does not provide the desired answer on the effective link between the variables; see the examples given in Sugihara et al. (2012). These authors therefore propose to approach the problem of causality in systems governed by deterministic dynamical systems by considering that the variables are indeed sharing the same attractor and that these variables can therefore provide information on each other. In addition, if an external driver is forcing the dynamical system under interest, the knowledge on the dynamics of the system can provide information on the driver, but not the opposite. It must be stressed here that the notion of causality is used in the very specific context of dynamical systems theory discussed above, that should not be confused with the traditional cause–effect relationship. From now on, the words “causality” and “dependences” will be used equivalently.

The original method proposed by Sugihara et al. (2012) is based on the Takens' reconstruction theorem: given a time delay τ, an embedding dimension E of an Euclidean space and the time series of a variable X, a reconstructed attractor, M_x, can be built. Each variable of that state space corresponds to a given delay, and each point of the reconstructed attractor is obtained from the time series as follows: $\mathit{X}\left(t\right)=X\left(t\right),X(t-\mathit{\tau }),X(t-\mathrm{2}\mathit{\tau }),\mathrm{\dots },X(t-(E-\mathrm{1}\left)\mathit{\tau }\right)$. For each state space point, X(t), of the reconstructed attractor, a set of close points are selected, called analogs, based on a distance, typically the Euclidean distance $d_i=\sqrt{{\sum }_j\left(X_j\right(t)-X_{j,i}(t){)}^{\mathrm{2}}}$, where X_j(t) and X_j,i(t) are the (delay) coordinates of the reference point and the ith analog, respectively. This distance can be used provided the solutions of the dynamical system display sufficiently smooth properties (e.g., continuity) as it is often assumed for large-scale geophysical fluid dynamics. Other distances could also be used. Since these are close to X(t), they should share some dynamical properties that can be exploited to make predictions starting from the current situation at time t. This set of analogs can also be used to recover the value of another variable, say Y(t), contemporary to X(t). The idea is to use the analogs found around X(t) to predict the expected variable Y(t), denoted $\widehat{Y}\left(t\right)$, as

$$\begin{array}{}\text{(1)}& \widehat{Y}\left(t\right)=\sum _{i=\mathrm{1}}^{E+\mathrm{1}}{Y}_i\times w_i,\end{array}$$

using weights defined as

$$\begin{array}{}\text{(2)}& w_i={\displaystyle \frac{\exp \left(-\frac{d_i}{min{d}_j}\right)}{{\sum }_i\exp \left(-\frac{d_i}{min{d}_j}\right)}},\end{array}$$

where the distances d_i are associated with the E+1 analogs obtained for the variable X (Tsonis et al., 2015), and Y_i is the values of Y contemporary to the ith analog on M_x. The number of analogs, E+1, is chosen such that one can form a simplex around the point X(t). The quantity mind_j denotes the minimum of d_j of the $j=\mathrm{1},..,E+\mathrm{1}$ analogs around the reference point X(t). The development of this type of nonlinear forecasting traces back to several seminal papers; see, for instance, Casdagli (1991) and Elsner and Tsonis (1992) for a detailed discussion. In their original versions, more general weights w_i were proposed that should be fitted through a least squares approach (Casdagli, 1991). Sugihara and May (1990) proposed to simplify the approach by using a simpler variant based on exponential functions depending on the distance between the analogs and the reference point. This weighting penalizes analogs that are far from the reference point, and the normalization by the minimum distance allows for having weights based only on the relative distance. This technique works well as discussed in Sugihara and May (1990). Moreover, it does not need any additional parameter, implying that the approach is parsimonious.

The dynamical relation between the variables X and Y can then be studied by comparing the actual value Y(t) to the inferred value $\widehat{Y}\left(t\right)$ obtained using analogs of X(t) on the M_x attractor. Repeating this method starting from different times t, the correlation coefficient between Y(t) and $\widehat{Y}\left(t\right)$ can be computed¹:

$$\begin{array}{}\text{(3)}& \mathit{\rho }={\displaystyle \frac{\mathrm{cov}(Y,\widehat{Y})}{{\mathit{\sigma }}_Y{\mathit{\sigma }}_{\widehat{Y}}}},\end{array}$$

where $\mathrm{cov}(.,.)$ denotes the covariance between the variables Y(t) and $\widehat{Y}\left(t\right)$, and σ_Y and ${\mathit{\sigma }}_{\widehat{Y}}$ are their standard deviations. The values of ρ are thus lying between −1 and 1, and it is also known as the Pearson correlation coefficient. High values of ρ indicate that the estimation of Y is good. However, this correlation does not necessarily mean that there is causality, as already emphasized in the introduction. For instance, there could be a confounding factor Z that influences both X and Y in the same manner. In that case, X and Y would behave similarly, and therefore there will be a correlation between Y and $\widehat{Y}$ (Sugihara et al., 2012).

To solve that problem, the method as described above can be repeated for increasing sizes L of the samples. For each sample of size L, an attractor is built, from which analogs can be isolated and the correlation coefficient can be computed. Note that different ways to build these L-size libraries were proposed without substantial differences (Ye et al., 2015). We therefore adopt the approach proposed in Tsonis et al. (2015) by randomly selecting a set of L events. If Y influences X, the effect of Y will be present in the reconstructed attractor M_x. By increasing the length L, more information on the time series of X is gathered, and therefore the selection of the analogs on the attractor is better. If there is a causality relation of Y on X, ρ will increase with L.

On the contrary, if there is no causality, the added information on the variable X will not give any information regarding Y, and the correlation coefficient will not increase with L. For instance, if a confounding factor Z affects both X and Y (that are otherwise dynamically independent of each other), they will contain similar information, and the inference of $\widehat{Y}$ based on X will display a correlation with Y which is independent of L (Sugihara et al., 2012). This provides a criterion on the role of a variable Y on X. Another important behavior of ρ(L) as a function of L is that the rate of increase is related to the strength of coupling. Note also that in an ideal context where the attractor can be reconstructed with precision and for L going to infinity, the correlation should converge to 1. In practical situations, this precision and the asymptotic limit are never reached. The convergence is then limited to a certain level by the presence of observational error, the approximation of the dynamics (like when a low-dimensional approximation is made of the full system) and the length of the series, L.

Figure 1 Monthly-averaged time series of the projections of the Atlantic, Pacific and tropical fields on the dominant modes of the dynamics, as obtained from the ORA-20C ocean reanalysis and the ERA-20C atmosphere reanalysis. (a–c) The geopotential at 500 hPa projected on F₁, the ocean potential temperature at 5 m deep projected on ϕ₂ and the sea surface height projected on ϕ₂ for the Atlantic. (d–f) As for the top row but for the Pacific. (g–i) Zonal velocity at 200 and 500 hPa, and the ocean potential temperature at 5 m deep averaged over the Nino3.4 region. All time series are standardized (with a mean equal to 0 and a variance equal to 1).

Figure 2 Monthly time series of the projections of the Atlantic, Pacific and tropical fields on the dominant modes of the dynamics, as obtained from the ORAS4 ocean reanalysis and the ERA-20C atmosphere reanalysis. (a–c) The geopotential at 500 hPa projected on F₁, the ocean potential temperature at 5 m deep projected on ϕ₂ and the sea surface height projected on ϕ₂ for the Atlantic. (d–f) As for the top row but for the Pacific. (g–i) Zonal velocity at 200 and 500 hPa, and the ocean potential temperature at 5 m deep averaged over the Nino3.4 region. All time series are standardized.

The CCM method requires the knowledge of the embedding dimension E and the time delay necessary for reconstructing the attractor M_x from the time series. Estimating the embedding dimension based, for instance, on the estimates of the correlation dimension of the attractor is very challenging when the expected embedding dimension is high since the approach needs to select close analogs to work properly (Kantz and Schreiber, 1995). It therefore needs very long time series that are usually not affordable (Nicolis, 1998; Van den Dool, 1994). Thus, a way to overcome this problem is to increase progressively the embedding dimension and see whether the results are robust or not. For the delay τ, one usually uses a time period for which successive situations become sufficiently decorrelated but not too much. Different methods are usually proposed to evaluate this delay, for instance, based on decorrelation times, or simply by trial and error (Casdagli, 1991; Parker and Chua, 1989). In the present cases, these delays should be relatively short for the atmosphere but much longer for the ocean as it can be guessed by inspecting the time series of the right panels of Figs. 1 and 2. For the latter, we are therefore facing an important problem since the decorrelation time (or delay) is not substantially smaller than the length of the time series.

A practical alternative is to build an attractor from a set of contemporary variables that are relevant to the dynamics from an expert evaluation. In such a case, a set of E=N variables at the same time t are used as entries of X to represent the attractor (in fact a projection of the full attractor in a subspace of N variables), and the analogs around a specific state space point X(t) can then be found in the same way as above. These analogs can be used to define the weights (Eq. 2) that are in turn used to predict Y(t). The influence of Y on X can then be inferred by computing Eq. (3). In order to see what the impact of the modification of the approach is, it has been applied in the context of a low-order system, a coupled ocean–atmosphere model of 36 ordinary differential equations developed in Vannitsem (2015), for which some results have been reported in Appendix B. One important result is the ability of the method to isolate dominant links between the projected attractor (the target of the analysis) and specific variables. The nature of the link can sometimes be directly related to terms present in the dynamical equations but not always due to the multivariate construction of the analogs on the projected attractor. Likewise, the absence of relationship inferred from the CCM in the present framework does not imply that there could not be some dependences when other projections of the full state space are used. The conclusions reached are therefore dependent on the specific configuration used and other experimental designs are necessary to corroborate the conclusions. This is planned for a future investigation.

Overall, this analysis demonstrates that the approach should be able to isolate important dependences between variables. It provides some confidence in the CCM algorithm, but we should keep in mind that the system explored in Appendix B is relatively simple and the application of CCM on more sophisticated climate models is worth performing. This is left for a future study whose results will be compared with the ones of the present analysis. Note also that when the series are much shorter than discussed in Appendix B, correlation can also be negative, indicating that the total length of the series has an important impact. This will also be further discussed in Sect. 4.1.

The CCM method with this modification will be applied on the data presented in the next section. Note that in order to evaluate the impact of the random sampling of L events in the datasets, we can repeat the sampling a certain number of times and infer a mean (or a median when strong asymmetries are present) and a standard deviation (here using the Fisher Z test). In the experiments that will be described below, this approach is adopted and each correlation value is estimated over a large number of samples. The algorithm is sketched in Appendix A.

3 Time series based on reanalysis datasets

The dynamics of the coupled ocean–atmosphere system has been recently investigated by adopting a novel approach which finds its roots in the low-order modeling of such dynamical systems (Vannitsem, 2015; Vannitsem et al., 2015). It consists in projecting key fields of the large-scale dynamics of the system on a few sets of modes that are dominating its dynamics. In Vannitsem and Ghil (2017), the coupling between the ocean and the atmosphere over the Atlantic has been investigated by projecting the geopotential at 500 hPa on the mode $F_{\mathrm{1}}=\sqrt{\mathrm{2}}\cos (\mathit{\pi }y/L_y)$, and the ocean potential temperature field at a certain depth (close to the surface) and the sea surface height on the mode ${\mathit{\varphi }}_{\mathrm{2}}=\mathrm{2}\sin (\mathit{\pi }x/L_x)\sin (\mathrm{2}\mathit{\pi }y/L_y)$. Note the sea surface height is a proxy for the upper-layer ocean stream-function field. F₁ is one of the largest-scale Fourier modes of the atmospheric field that is confined in an x-periodic β channel with free-slip boundary conditions in the y direction, while ϕ₂ is one of the dominant Fourier modes compatible with free-slip boundary conditions in a rectangular, L_x×L_y closed basin. The latter mode corresponds to the typical structure of a double gyre in such a closed ocean basin. We thus expect the projection of the geopotential on F₁(x,y) to provide information on the intensity of the large-scale eastward zonal transport in the atmosphere, while the projection of the temperature and stream-function field in the ocean on ϕ₂(x,y) will allow us to evaluate the strength of the dominant component of the meridional gradient of temperature and the intensity of the double-gyre dynamics in the ocean, respectively. The domain chosen in terms of the spherical coordinates is 55^∘ W $\le \mathit{\lambda }\le$ 15^∘ W, 25^∘ N $\le \mathit{\varphi }\le$ 60^∘ N, with $(x=\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{0}},y=\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}})$, here (λ₀=305^∘, ϕ₀=25^∘) and $(L_x={\mathrm{40}}^{\circ },L_y={\mathrm{35}}^{\circ })$. Note that the domain used here is the same as in Vannitsem and Ghil (2017) but a typographical error on the domain of projection is reported in the Supplement of Vannitsem and Ghil (2017). The time series obtained for the North Atlantic will be denoted NA${}_{{\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}}}$, NA${}_{{\mathit{\theta }}_{\mathrm{o},\mathrm{2}}}$, NA${}_{{\mathit{\eta }}_{\mathrm{o},\mathrm{2}}}$ for the projection the geopotential at 500 hPa on the mode F₁, the ocean potential temperature field at 5 m deep and the sea surface height on the mode ϕ₂, respectively.

A similar approach can be performed for the North Pacific, except that the domain is now larger in the zonal direction. In this case, the spherical-rectangle domain is 165–225^∘ E, 25–60^∘ N. The series obtained for this domain will be denoted NP${}_{{\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}}}$, NP${}_{{\mathit{\theta }}_{\mathrm{o},\mathrm{2}}}$, NP${}_{{\mathit{\eta }}_{\mathrm{o},\mathrm{2}}}$ as for the North Atlantic. Note that for both basins the projected time series contains the dominant part of the variability. It however does not preclude that other important processes are missing in the description here. Further analyses with more modes are certainly worth doing in the future.

For the tropical Pacific, one can wonder what kind of variables should be considered. First, a dominant variable in the tropical Pacific is the mean temperature of the upper ocean layer, known to be associated with the dynamics of El Niño. It is also known that the Walker circulation is considerably affected by the upper-layer ocean temperature, and vice versa (Philander, 1990). Let us therefore consider for now the mean ocean potential temperature in the Nino3.4 region, known to have strong correlation with the variability in the North Atlantic region (Brönnimann, 2007). For the characterization of the Walker circulation, the zonal winds at 500 and 200 hPa over the same domain are chosen. They provide some information on the position and the strength of the Walker circulation over the Nino3.4 region. The series obtained will be denoted as NI_U200, NI_U500, NI${}_{{\mathit{\theta }}_{\mathrm{o},\mathrm{av}}}$.

This approach of reducing the dynamics of the ocean and the atmosphere to a few spectral large-scale components may at first sight look arbitrary. However, for the two midlatitude basins, these modes possess the largest amplitudes (Vannitsem and Ghil, 2017), and for the tropical Pacific, it is known that these large-scale flows are strongly affected by the interaction between the ocean and the atmosphere. Moreover, we are interested in the basin-scale interaction between midlatitudes and the tropics. If such an interaction exists, we expect that these should be visible through the analysis of these large-scale fields. It is clear that these specific variables do not represent the full dynamics, and additional analyses with more modes are worthwhile, in particular to see what is the role of the main currents present in the ocean like the Gulf Stream or the Kuroshio Current.

Three different reanalysis datasets from ECMWF are used. The ERA-20C dataset provides a continuous reanalysis for the atmosphere of the 20th century which assimilates observations from surface pressure and surface marine winds only. It is produced using the Integrated Forecasting System (IFS) model cycle Cy38r1, and detailed information can be found on the website of ECMWF at https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-20c (last access: 21 August 2018); see also the report on the quality of this reanalysis dataset (Poli et al., 2015). It covers the period 1900–2010.

The second dataset is the Ocean reanalysis ORAS4 obtained using the Nucleus for European Modelling of the Ocean (NEMO) model. The ocean model is forced by the heat, momentum and freshwater fluxes at the upper-surface and ocean observations. For the upper-surface fluxes, the ERA40 reanalysis dataset is used from September 1957 to December 1988, then ERA-Interim from January 1989 to December 2009. In 2010, these fluxes are provided by the ECMWF operational analyses. For the sea surface temperature (SST) and ice products, ERA40 (until December 1981) and the Reynolds dataset are used. Finally, observational data within the ocean are the temperature and salinity profiles from September 1957 to December 2010, and the sea level anomalies from November 1992 onward. More information on the datasets used and the model configuration can be found at https://www.ecmwf.int/en/research/climate-reanalysis/ocean-reanalysis (last access: 21 August 2018), and more information on the quality of this product can be found in Balmaseda et al. (2013). The period covered by the dataset used in the present study is fixed from January 1958 to December 2010.

Finally, the third reanalysis dataset used is ORA-20C which is a 10-member ensemble of ocean reanalysis covering the 20th century using atmospheric forcing from ERA-20C. This dataset is more homogeneous than ORAS4 since the atmospheric forcing is consistent during the whole 20th century. As pointed out in de Boisséson and Balmaseda (2016), the uncertainty is large during the first part of the century before the assimilation process constrains all the members of the ensemble to a state more consistent with other reanalysis products. One can suspect this dataset to be better than ORAS4 for the dynamics of the ocean during the second part of the century since the state of the ocean has gotten some time to adjust toward a representative climatology during the first half. We will therefore use data from January 1958 consistently with the ORAS4 dataset to December 2009, the final date of data availability at the time when the present work has been conducted.

Thus, the atmospheric data from ERA-20C that will be used in the present work will cover the same periods as the ones fixed for ORAS4 and ORA-20C, respectively. All data used are monthly values.

The different monthly-averaged time series obtained by projecting the fields on the two Fourier modes are grouped by zones: three for the North Atlantic (containing one series for the atmosphere and two series for the ocean), three for the North Pacific (as for the North Atlantic) and three for the tropical Pacific (two series for the atmosphere and one series for the ocean). The nine time series based on the ERA-20C and ORA-20C reanalyses are displayed in Fig. 1 for the three regions. The three series in each zone will constitute a three-dimensional projection of the local coupled ocean–atmosphere dynamics. The same projections but using ERA-20C and ORAS4 are displayed in Fig. 2.

Table 1 Correlation coefficients between the different variables for the two reanalysis datasets. On the left are ORAS4 and ERA-20C; on the right are ORA-20C and ERA-20C.

Let us first briefly investigate the covariance structure of these time series. Table 1 displays the covariances between the different time series of Figs. 1 and 2, with the covariances for the series of Fig. 1 on the left side of each column, and the ones corresponding to series of Fig. 2 on the right side. There are a few remarkable correlations. First, there are the ones between the atmospheric fields NA${}_{{\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}}}$, NI_U200, NI_U500 and NP${}_{{\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}}}$, suggesting that some key variables of the global dynamics have been selected. The ocean temperature modes, NA${}_{{\mathit{\theta }}_{\mathrm{o},\mathrm{2}}}$, NP${}_{{\mathit{\theta }}_{\mathrm{o},\mathrm{2}}}$ and NI${}_{{\mathit{\theta }}_{\mathrm{o},\mathrm{av}}}$, are also highly correlated. Interestingly, the NP${}_{{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}}}$ is anti-correlated with NI${}_{{\mathit{\theta }}_{\mathrm{o},\mathrm{av}}}$. Another remarkable result is that the transport in the North Pacific, NP${}_{{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}}}$, is highly correlated with the transport in the North Atlantic, NA${}_{{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}}}$, although different amplitudes are found for the two ocean reanalysis datasets. Some other correlations are much less robust when one looks at the two different ocean datasets, in particular associated with the transport and the ocean temperature over the North Atlantic (second and third columns in Table 1). These differences should be associated with the different approaches to force the ocean model and reflect important uncertainties in reconstructing the past evolution of the Earth system.

Figure 3 CCM as a function of the length L of the samples, as obtained from the monthly time series displayed in Fig. 1 for the ERA-20C and ORA-20C reanalyses. Each line with symbols corresponds to the influence of one variable on a specific coupled ocean–atmosphere basin. The specific variables are denoted in the caption corresponding to each line in each panel.

Important correlations appear in the datasets explored, suggesting that common information is present in the different coupled ocean–atmosphere basins discussed here. These correlations are presumably highly dependent on the seasonal cycle affecting the Earth system. Further analysis by removing the seasonal signal could be done to clarify the correlations between the anomalies found in each basins. We will not, however, go in that direction in the present work since there are several ways to do it and since the seasonal signal is part of the dynamics itself. It suffices here to recognize that links exist between these basins, whose nature will be clarified by using the CCM algorithm discussed in Sect. 2. It will be shown that the annual cycle is also affecting the CCM results and two different ways to disentangle its role in the causality analysis will be proposed.

4 Results of the application of convergent cross mapping

4.1 Reanalyses: ERA-20C/ORA-20C

Let us start the analysis by investigating the CCM for the series displayed in Fig. 1. In Fig. 3a, the correlation ρ(L) of the inferred variable $\widehat{Y}$ with the actual solution Y is shown for the three variables of the tropical Pacific. These inferred variables are obtained by building analogs in the North Atlantic as described in Sect. 2. Note that a 95 % confidence interval is provided based on the Fisher Z-transform test based on resampling a certain number of times the samples of length L, indicating that the influence of the two atmospheric tropical series is significant. This confidence interval for large L is larger than for smaller values due to the fact that we are reaching the limit of the number of data points. Note also that for this analysis, the analogs have been selected to be at least separated by a period of 12 months. Longer exclusion periods have been tested without substantial differences.

An increase of the correlation is found as a function of L, suggesting that the North Atlantic depends on the two large-scale atmospheric variables selected for the tropical Pacific. The correlation for the ocean temperature of the tropical region is however negative. As already mentioned previously, this feature may occur when the time series is too short and when there is no significant coupling between the variables. This suggests that the prediction of the tropical temperature based on analogs over the Atlantic is very poor and therefore indicates the absence of influence of the tropical ocean temperature. Thus, the dominant influence is from the zonal flows in the atmosphere, associated with the dynamics of the Walker circulation.

The influence of the three variables of the North Pacific on the North Atlantic is very important, as shown in Fig. 3b with the three CCM increasing and significantly positive. The North Pacific ocean dynamics has a larger influence than the North Pacific ocean temperature on the Atlantic, as reflected by the larger amplitude of the CCM values. Note that the reduction of the state space coordinates associated with X from three to two also provides interesting results with a dominating influence from the atmospheric tropical Pacific variables on the two-dimensional projection (NA${}_{{\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}}}$, NA${}_{{\mathit{\theta }}_{\mathrm{o},\mathrm{2}}}$). However, the increase before saturation of ρ(L) is much more limited than when using the three variables to build the North Atlantic projection of the attractor (not shown). The latter analysis suggests that the dependences between these regions are better elucidated based on the three-dimensional space.

The impacts of the North Atlantic and the North Pacific regions on the tropical Pacific are displayed in Fig. 3c and d. All CCM values of Fig. 3c are almost flat as a function of L, suggesting that even when it is positive (very significant for the geopotential at 500 hPa), there is no dependence of the tropical Pacific on the dynamics over the Atlantic. A slightly different picture emerges for the influence of the North Pacific on the tropical Pacific, with a slightly increasing CCM for the geopotential at 500 hPa over the North Pacific suggesting an influence of the upper-air dynamics over the North Pacific on the tropical Pacific. Note that sometimes it is difficult to have a clear-cut answer on the increase or not of a correlation as a function of L. There is a degree of arbitrariness that should be alleviated using other approaches as the ones that will be discussed later based on the temporal averaging or based on surrogates.

In Fig. 3e and f, the CCM values characterizing the influence of the North Atlantic (Fig. 3e) and tropical Pacific (Fig. 3f) on the North Pacific are displayed. A clear increase of CCM associated with the influence of the North Atlantic ocean dynamics on the North Pacific is found (Fig. 3e). The CCM values for the two other variables are also slightly increasing with L, suggesting a dependence of these two variables on the North Pacific. Interestingly, the influence of the geopotential of the North Pacific on the North Atlantic (Fig. 3b) is more important than the one from the North Atlantic to the North Pacific since the correlation is higher. This asymmetry looks reasonable since the dominant flow is eastward in the northern extratropics. Finally, Fig. 3f suggests that the atmospheric zonal flow over the tropical Pacific influences the North Pacific dynamics but not the ocean temperature.

Figure 4 CCM as a function of the length L of the samples, as obtained from the series of Fig. 1 after averaging over a sliding window of 12 months for the ERA-20C and ORA-20C reanalyses. Each line with symbols corresponds to the influence of one variable on a specific coupled ocean–atmosphere basin. The specific variables are denoted in the caption corresponding to each line at each panel.

To summarize, the analysis reveals that (i) the upper-air tropical Pacific dynamics and the North Pacific ocean and atmosphere dynamics influence the dynamics over the North Atlantic; (ii) the upper-air North Pacific dynamics influences the tropical Pacific; and (iii) the North Atlantic ocean dynamics and the upper-air tropical Pacific dynamics influence the dynamics over the North Pacific. These results seem to support the view of several authors, e.g., Sardeshmukh and Hoskins (1988), Alexander et al. (2002) and Lau (2016), that there is an atmospheric bridge dependence from the tropical Pacific to the North Atlantic, via the North Pacific. Furthermore, a dependence between the North Pacific and the North Atlantic emerges, which is mostly oriented from the North Pacific to the North Atlantic since the correlations are larger in Fig. 3b than in Fig. 3e for the atmospheric variable, in line with the findings of Drouard et al. (2015). However, another very interesting dynamical signature emerges, suggesting that the North Pacific and the North Atlantic ocean dynamics are dependent on each other. A possible candidate for this coupling is the thermohaline planetary circulation that affects both regions of the globe.

To disentangle the importance of the thermohaline circulation which displays a variability on very long timescales, one can investigate the dependence properties when longer timescales are considered. An average of the dataset has been performed using a sliding window of 12 months. This approach allows for removing most of the impact of the annual signal, while keeping a number of data points large enough to perform the CCM analysis. Larger windows could be used but it will introduce very long time correlations that could penalize the selection of analogs since one needs analogs that are sufficiently uncorrelated in time. The results are displayed in Fig. 4. A first remarkable result is the fact that for a 1-year average, the mutual dependence of the dynamics over the North Atlantic and Pacific is dominated by the ocean dynamics (see Fig. 4b and e). The CCM values of the other variables in these two panels are flat and close to or below 0. At the same time, new dependences emerge between the North Pacific and the tropical Pacific (Fig. 4d), in particular for the ocean transport. This further supports the conjecture that the three regions are coupled via the large-scale global ocean dynamics. This is presumably linked with the thermohaline circulation, but further analyses using additional variables and tropical regions are needed in order to disentangle this point. This will be the subject of a future work.

Figure 5 CCM as a function of the length L of the samples, as obtained from surrogate three-dimensional attractors built by superimposing random anomalies to the annual climatological cycle, based on the data from ERA-20C and ORA-20C reanalyses. Two different surrogates have been used. Each line with symbols corresponds to the influence of one actual variable (not a surrogate series) on a specific coupled ocean–atmosphere surrogate attractor. The specific variables are denoted in the caption corresponding to each line in each panel.

Another very important finding in Fig. 4 is the fact that CCM values are now close to 0 (and do not display any dependence in L) for the influence of the tropical Pacific on the North Atlantic, and vice versa (Fig. 4a and c). It suggests that the dependence between these two regions is confined to timescales smaller than a year. One can therefore wonder whether this dependence is associated purely with the annual cycle or with some specific tropical events like El Niño or La Niña, or in other words, if it is mostly associated with the climatology of the tropical Pacific or not.

To clarify this point, surrogate three-dimensional attractors have been built by using the monthly means (averaged over the different years) on which random anomalies with the appropriate variance are superimposed. These random anomalies were simulated assuming a Gaussian distribution around each monthly mean. The variance of the distribution is estimated using the anomalies of the corresponding month for all years in the datasets.

These new attractors are then used to predict the true variable of interest Y. Figure 5 displays the corresponding panels that should be compared to Fig. 3. Two different random surrogates have been built, implying that two curves are displayed in each panel for each variable, Y. The first remarkable result is that the CCM values of the ocean dynamics variables of Fig. 5b and e are now flat and close to 0, suggesting that the CCM values for the transport in the two northern ocean basins found in Fig. 3 are indeed indicating a dynamical coupling between the two basins beyond the annual climatological variations. Also in Fig. 5b, the CCM values for the influence of the North Pacific series (ocean temperature and geopotential at 500 hPa) on the surrogate attractor of the Atlantic considerably decrease, suggesting the importance of the influence of the North Pacific on the North Atlantic beyond the annual climatological influence.

However, when looking at the results in Fig. 5a, the CCM values based on the use of the surrogate attractors are very close to the one obtained with the actual attractors. This surprising feature suggests that there is no influence between the tropical Pacific and the North Atlantic beyond the annual climatological variations, or in other words, that the tropical Pacific variability does not influence the anomalies over the North Pacific. This has a very strong implication in the sense that there is no dynamical link between an event like El Niño or La Niña and the anomalies over the North Atlantic. A similar picture is found for the influence of the tropical Pacific on the North Pacific as illustrated in Fig. 5f.

Finally, to further test the robustness of these results, a complementary way to clarify whether monthly anomalies between different basins are indeed related to each other is to directly apply CCM to these anomalies. The results are displayed in Appendix C (Fig. C1). The same conclusions are reached with the absence of dependences between the variables over the North Atlantic basin and the tropical Pacific (Fig. C1a and c), and a strong mutual dependence between the ocean dynamics over the North Atlantic and the North Pacific (Fig. C1b and e).

Figure 6 CCM as a function of the length L of the samples, as obtained from the monthly time series displayed in Fig. 2 for the ERA-20C and ORAS4 reanalyses. Each line with symbols corresponds to the influence of one variable on a specific coupled ocean–atmosphere basin. The specific variables are denoted in the caption corresponding to each line in each panel.

In summary, in the limit of the data at our disposal, the analysis suggests that the anomalies associated with the dynamics over the North Atlantic and North Pacific cannot be inferred based on the variability of the variables we have used so far in the tropical Pacific. Note that the previous analysis is made for all seasons and without distinctions between certain types of events, say strong El Niño events, as it is usually done when analyzing the effect of the El Niño–Southern Oscillation (ENSO) over other regions of the globe, e.g., Brönnimann (2007). Such a split between seasons and/or events is worth performing, but the time series are already short and the selection of certain events will considerably reduce the statistics. The absence of a link between the tropical Pacific and the North Atlantic coupled dynamics could also reflect the non-stationary properties of the teleconnections between the North Atlantic and the tropical Pacific as documented in López-Parages et al. (2016), Johnson and Kosaka (2016) and Goss and Feldstein (2017) and references therein. The analysis of long climate runs of state-of-the-art models with the approach used here would be very useful in that respect.

4.2 ERA-20C/ORAS4

Let us now consider the second ocean reanalysis, ORAS4, while keeping the same atmospheric reanalysis. This second investigation should allow us to clarify the robustness of our findings. Figure 6 shows the results of the computation of CCM that should be compared with the results presented in Fig. 3. One remarkable feature is the absence of dependences between the North Atlantic and the North Pacific for the ocean dynamics; see black full circles in Fig. 6b and e. This result considerably differs from the one obtained with ORA-20C. Another important difference is a larger amplitude of the influence of the tropical Pacific ocean temperature on the North Pacific and of the North Pacific variables on the North Atlantic. The other dependences are more robust.

Figure 7 CCM as a function of the length L of the samples, as obtained from the series of Fig. 2 after averaging over a sliding window of 12 months for the ERA-20C and ORAS4 reanalyses. Each line with symbols corresponds to the influence of one variable on a specific coupled ocean–atmosphere basin. The specific variables are denoted in the caption corresponding to each line in each panel.

These results suggest that the dynamics within the ocean differs considerably between ORAS4 and ORA-20C as already suggested by the covariances displayed in Table 1. These differences are probably due to (i) the fact that these reanalyses are obtained with different atmospheric forcing, specifically ORA-20C with ERA-20C and ORAS4 with different atmospheric reanalysis products, and to (ii) the fact that for ORA-20C the ocean model started at the beginning of 1900, while for ORAS4 it started at the end of 1957. We may conjecture that the ORA-20C reanalysis dataset is more reliable since a more consistent atmospheric forcing has been applied and the ocean model has more time to equilibrate around its climate. However, care should be taken here in drawing definitive conclusions on that. A better approach to disentangle which of these reanalyses provides the correct answer is to investigate a full coupled ocean–atmosphere reanalysis obtained for the whole 20th century.

The investigation of the dependences for sliding averages over 12 months displayed in Fig. 7 also suggests a suppression of the dependences for most of the variables, at the exception of the one associated with the influence of the North Pacific ocean temperature on the tropical Pacific (Fig. 7d). It is also remarkable that a dependence emerges from the North Atlantic ocean dynamics to the North Pacific but much weaker than for the other ocean reanalysis dataset.

Finally, for the sake of completeness, the application of the CCM to the monthly anomalies is displayed in Fig. C2. Here, as for the ERA-20C/ORA-20C dataset, no link between the anomalies of the tropical Pacific and the North Atlantic is found (Fig. C2a and c). An important difference is however visible in the influence of the North Pacific ocean temperature on the North Atlantic and on the tropical Pacific (Fig. C2b and d). Again, this contrasts with the results obtained with the other reanalysis dataset.

5 Conclusions

The causality between the dynamics of three different coupled ocean–atmosphere basins, the North Atlantic, the North Pacific and the tropical Pacific region (Nino3.4), has been explored using data from three different reanalysis datasets, ORA-20C, ORAS4 and ERA-20C. The approach used is convergent cross mapping, developed by Sugihara et al. (2012), which allows to go beyond the classical teleconnection patterns and which provides a clear signature of the interdependences between series or regions. The analysis reveals a few very important facts that should help improving our understanding of the remote influence of large-scale dynamical processes and in particular the impact of the tropical Pacific coupled dynamics on the extratropics:

• The tropical Pacific coupled ocean–atmosphere dynamics does not seem to have an impact on the extratropics beyond the annual climatological cycle. This very surprising result suggests that there is little hope to improve predictability in the extratropics based on information on the variability in the tropical Pacific. This result needs however more attention and a thorough inspection of other datasets and long climate model runs in order to be confirmed. It is possible in particular that a more detailed analysis based on the selection of specific events, like strong El Niño or La Niña events, will provide new information on the interactions between the tropical Pacific and the rest of the world, beyond the climatological annual signal. However, such an investigation needs considerably more data than the ones used in the present work and therefore calls for the development of even longer coupled reanalyses or the use of long climate model runs.

• The atmosphere over the North Pacific considerably influences the North Atlantic (beyond the climatological annual signal), in agreement with the results found, for instance, in Drouard et al. (2015).

• The results presented here for the two ocean reanalysis datasets disagree on the nature of the dependence between the North Pacific basin and the other ones. ORA-20C indicates a dynamical influence (transport), while ORAS4 suggests more an influence dominated by the ocean temperature. This difference is probably related to the specific data assimilation setup used in each case. To clarify what is really at work here, new reanalysis datasets are needed. We can conjecture that a more reliable one could be built by using a coupled ocean–atmosphere data assimilation system running continuously on the longest period possible.

• The interdependence between the North Atlantic and the North Pacific on longer timescales than a year seems to be important and is probably related to the coupled ocean–atmosphere dynamics on long timescales. One could conjecture that the thermohaline circulation should play a role on this link. Additional analyses with longer datasets, with climate model runs, but also with the analysis of additional basins like the tropical Atlantic or the Indian monsoon region, are necessary to clarify this role.

The present work has demonstrated the urgent necessity to go beyond the teleconnection patterns for the investigation of the interaction between the different components of the climate system, using tools recently developed in the context of nonlinear sciences. Teleconnection patterns do provide information on a co-variability (which is an interesting piece of information per se) but not on the influence of a region on another through a dynamical coupling. A common forcing can, for instance, induce a correlation between two variables even if these are perfectly independent in a dynamical sense.

New analyses will be performed along the lines drawn above, in particular in climate model runs. Several long control climate runs of CMIP5 models are available and can be analyzed in the same perspective as in the present work, in order in particular to evaluate the impact of ENSO on the dynamics of the North Pacific and the North Atlantic.

Finally, it should be stressed that the specific design of the CCM approach adopted here based on a low-dimensional projection of the full ocean–atmosphere attractor does not allow to have a one-to-one correspondence of the influence of one variable on another. Moreover, other subspaces could be used that can provide different results. More variables should then be considered, such as projections on additional Fourier modes, or using projections on a few empirical orthogonal functions. These analyses will allow to evaluate the robustness of the present results.

Appendix A: The CCM algorithm

The different steps of the CCM algorithm used to compute the correlation coefficient are detailed in Fig. A1.

Figure A1 Algorithm used to compute the CCM presented in Sect. 2, based on the series discussed in Sect. 3.

Appendix B: CCM applied to an idealized model

To test the CCM technique described in Sect. 2, we use a dynamical system recently developed in Vannitsem (2015). It consists of a set of 36 ordinary differential equations representing the large-scale dynamics of a coupled ocean–atmosphere system at midlatitudes. The equations are described in Vannitsem (2015) and in its Supplement. The atmospheric model is based on the vorticity equations of a two-layer, quasi-geostrophic flow defined on a β plane. The ocean dynamics is based on the reduced-gravity, quasi-geostrophic shallow-water model on a β plane. For the ocean, it is assumed that temperature is a passive scalar transported by the ocean currents, but the oceanic temperature field displays strong interactions with the atmospheric temperature through radiative and heat exchanges.

All fields are developed in Fourier series on a β plane as

$$\begin{array}{}\text{(B1)}& {\displaystyle }& {\displaystyle }\mathit{\delta }{T}_{\mathrm{o}}=\sum _{i=\mathrm{2}(\mathrm{and}\ne \mathrm{5})}^{\mathrm{8}}{\mathrm{\Theta }}_{\mathrm{o},i}{\mathit{\varphi }}_i,{\mathrm{\Psi }}_{\mathrm{o}}=\sum _{i=\mathrm{1}}^{\mathrm{8}}{\mathrm{\Psi }}_{\mathrm{o},i}{\mathit{\varphi }}_i,\text{(B2)}& {\displaystyle }& {\displaystyle }{\mathit{\psi }}_{\mathrm{a}}=\sum _{i=\mathrm{1}}^{\mathrm{10}}{\mathrm{\Psi }}_{\mathrm{a},i}{F}_i,\mathit{\delta }{T}_{\mathrm{a}}=\sum _{i=\mathrm{1}}^{\mathrm{10}}{\mathrm{\Theta }}_{\mathrm{a},i}{F}_i,\end{array}$$

where ${\mathrm{\Theta }}_{\mathrm{a},i}=({\mathit{\psi }}_{\mathrm{a},i}^{\mathrm{1}}-{\mathit{\psi }}_{\mathrm{a},i}^{\mathrm{3}})/\mathrm{2}$ and ${\mathrm{\Psi }}_{\mathrm{a},i}=({\mathit{\psi }}_{\mathrm{a},i}^{\mathrm{1}}+{\mathit{\psi }}_{\mathrm{a},i}^{\mathrm{3}})/\mathrm{2}$, with ψ¹ and ψ³ the stream functions in the upper and lower layers of the atmosphere. Ψ_o is the stream-function field in the ocean. δT_o and δT_a are temperature anomaly fields with respect to spatially averaged reference temperatures. The modes ϕ_i used for the ocean are compatible with free-slip boundary conditions in a closed basin, while F_i are modes used for the atmospheric fields compatible with free-slip boundaries in the meridional direction and periodic boundaries in the zonal direction; see also the paper of Vannitsem et al. (2015).

The CCM analysis is performed on the solutions generated by the model with the same parameters as in Fig. 3 of Vannitsem (2015) with a surface friction coefficient C=0.006 kg m⁻² s⁻¹. The model is forced with seasonal variations of the solar input as discussed in Vannitsem (2015), and the solutions are averaged over 1 month (1∕12 of the 365 days of the model year). The three variables used for building the attractor and the analogs are $({\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}},{\mathrm{\Theta }}_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}})$ as for the North Atlantic and North Pacific datasets discussed in Sect. 3. Then several other variables are used to see if they have causality relations with the three variables used to build the attractor. The length of the time series is L=6000 months.

The results are shown in Fig. B1 for different variables and also for a sliding window of 12 months in Fig. B1e and f. As it can be seen in Fig. B1a and b, the only atmospheric variable (explored so far) influencing the dynamics of $({\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}},{\mathrm{\Theta }}_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}})$ is Θ_a,1, since the correlation is high and increases as a function of L. This variable is strongly linked to the dynamics of Ψ_a,1 in the equations of the model, since it is the only one influencing (linearly) the evolution of Ψ_a,1. The others do not seem to have a strong influence. In Fig. B1c and d, the impact of some ocean variables is illustrated with no apparent influence of Ψ_o,5 and Θ_o,3. However, all other variables have different levels of influence with a clear increase of the correlation as a function of L.

Finally, in Fig. B1e–h, the impact of using a sliding average is illustrated. First, the CCM plotted in Fig. B1f indicates that the influence of Θ_a,1 is not removed, suggesting its essential role on the dynamics at different timescales of motion. Second, the influence of the ocean modes which were found to play a role at monthly timescales is further enhanced.

This brief analysis of a low-order system based on the CCM algorithm described in Sect. 2 indicates that it is a powerful tool to isolate the influence of certain variables on others in the system. Note that the angle of approach adopted in Sect. 2 by considering a low-order projection of the full state space as a target does not allow to have detailed information on the nature of the coupling between the variables since what is inferred is a global influence of a variable on a subset of other variables. The analysis also opens new questions on the role of the different variables in this low-order model and the dependences on the specific subspace selected as the target. This problem is out of the scope of the present work but is worth pursuing in the future. We can now proceed with this approach in the context of the datasets presented in Sect. 3.

Figure B1 CCM as a function of the length L of the samples, as obtained from monthly time series of the low-order coupled ocean–atmosphere model integration. Panels (a) to (d) display the values for the influence of a set of model variables on $({\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}},{\mathrm{\Theta }}_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}})$ at a monthly timescale. Panels (e) to (h) display the values for the influence of a set of model variables on $({\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}},{\mathrm{\Theta }}_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}})$ after the application of a sliding average over 12 months. Each line with symbols corresponds to the influence of one variable on $({\mathrm{\Psi }}_{\mathrm{a},\mathrm{1}},{\mathrm{\Theta }}_{\mathrm{o},\mathrm{2}},{\mathrm{\Psi }}_{\mathrm{o},\mathrm{2}})$. The specific variables are denoted in the caption corresponding to each line in each panel.

Appendix C: Application of CCM to anomalies

The CCM is applied on the anomalies of the different time series displayed in Figs. 1 and 2. The anomalies are defined as $X(t,\mathit{\tau })-\mathit{\mu }\left(\mathit{\tau }\right)$, where μ(τ) is the average value over all years for month τ, and the argument t is the year for which the anomalies are computed. The analysis of these results is done in the core of the text.

Figure C1 CCM as a function of the size L of the samples, as obtained from the anomaly of the monthly time series displayed in Fig. 1 for the ERA-20C and ORA-20C reanalyses. Each line with symbols corresponds to the influence of one variable on a specific coupled ocean–atmosphere basin. The specific variables are denoted in the caption corresponding to each line in each panel.

Figure C2 CCM as a function of the size L of the samples, as obtained from the anomaly of the monthly time series displayed in Fig. 2 for the ERA-20C and ORAS4 reanalyses. Each line with symbols corresponds to the influence of one variable on a specific coupled ocean–atmosphere basin. The specific variables are denoted in the caption corresponding to each line in each panel.