Dominant patterns of interaction between the tropics and mid-latitudes in boreal summer: causal relationships and the role of timescales

Tropical convective activity represents a source of predictability for mid-latitude weather in the Northern Hemisphere. In winter, the El Niño–Southern Oscillation (ENSO) is the dominant source of predictability in the tropics and extratropics, but its role in summer is much less pronounced and the exact teleconnection pathways are not well understood. Here, we assess how tropical convection interacts with mid-latitude summer circulation at different intra-seasonal timescales and how ENSO affects these interactions. First, we apply maximum covariance analysis (MCA) between tropical convective activity and mid-latitude geopotential height fields to identify the dominant modes of interaction. The first MCA mode connects the South Asian monsoon with the mid-latitude circumglobal teleconnection pattern. The second MCA mode connects the western North Pacific summer monsoon in the tropics with a wave-5 pattern centred over the North Pacific High in the mid-latitudes. We show that the MCA patterns are fairly insensitive to the selected intra-seasonal timescale from weekly to 4-weekly data. To study the potential causal interdependencies between these modes and with other atmospheric fields, we apply the causal discovery method PCMCI at different timescales. PCMCI extends standard correlation analysis by removing the confounding effects of autocorrelation, indirect links and common drivers. In general, there is a two-way causal interaction between the tropics and mid-latitudes, but the strength and sometimes sign of the causal link are timescale dependent. We introduce causal maps that show the regionally specific causal effect from each MCA mode. Those maps confirm the dominant patterns of interaction and in addition highlight specific mid-latitude regions that are most strongly connected to tropical convection. In general, the identified causal teleconnection patterns are only mildly affected by ENSO and the tropical mid-latitude linkages remain similar. Still, La Niña strengthens the South Asian monsoon generating a stronger response in the mid-latitudes, while during El Niño years the Pacific pattern is reinforced. This study paves the way for process-based validation of boreal summer teleconnections in (sub-)seasonal forecast models and climate models and therefore works towards improved sub-seasonal predictions and climate projections.

1 Introduction

Tropical–mid-latitude teleconnections in boreal summer can have a great impact on surface weather conditions in the northern mid-latitudes (Ding and Wang, 2005; O'Reilly et al., 2018; Wang et al., 2001). Still, the direct influence of the El Niño–Southern Oscillation (ENSO) on the mid-latitude circulation is weaker in summer than in winter (Branstator, 2002; Schubert et al., 2011; Thomson and Vallis, 2018). Instead, in summer convective activity related to the Northern Hemisphere tropical monsoon systems can profoundly influence surface weather conditions in the mid-latitudes (Branstator, 2014; Ding and Wang, 2005; O'Reilly et al., 2018; Rodwell and Hoskins, 1996). In contrast, mid-latitude wave trains and cyclonic activity at intra-seasonal timescales can modulate the tropical monsoons and have been linked to extreme rainfall events in the Indian region (Lau and Kim, 2011; Vellore et al., 2014, 2016). Therefore, tropical and mid-latitude regions are likely connected in complex two-way teleconnection patterns operating at a range of sub-seasonal timescales (Di Capua et al., 2020; Ding and Wang, 2005, 2007).

During boreal summer, the South Asian monsoon (SAM), represents one of the most powerful features of the tropical/subtropical circulation. Characterized by heavy rainfall over central India and the Bay of Bengal, the SAM has strong intra-seasonal variability associated with alternating active and break phases, linked to the boreal summer intra-seasonal oscillation (BSISO, Choudhury and Krishnan, 2011; Gadgil and Joseph, 2003; Goswami et al., 1998; Krishnamurti and Surgi, 1987; Krishnan et al., 2000; Rao, 1976; Saha et al., 2012; Suhas et al., 2012). The western North Pacific summer monsoon (WNPSM) represents the Pacific counterpart to the SAM and is identified by strong rainfall over the sub-tropical western North Pacific (Li and Wang, 2005). Similar to the SAM, the WNPSM also exhibits strong intra-seasonal oscillations (Wang and Xu, 1997).

Latent heat release due to strong monsoonal rainfall can influence subtropical and mid-latitude regions via Rossby wave teleconnections. The SAM has been connected to subtropical arid conditions in the North African region via the so-called monsoon–desert mechanism, creating reinforced descending motions over the Sahara during strong SAM phases (Rodwell and Hoskins, 1996; Stephan et al., 2019). This mechanism is fairly well captured by both climate (Cherchi et al., 2014) and seasonal forecast models (Beverley et al., 2019).The SAM is also connected to mid-latitude circulation via its interaction with the circumglobal teleconnection pattern (CGT), a wave pattern with five centres of action encircling the northern mid-latitudes and affecting temperature and precipitation there (Ding and Wang, 2005; Kripalani et al., 1997). This wave-5 like CGT pattern can be identified through inter-annual to intra-seasonal (weekly) timescales, and it is likely connected with the SAM via two-way causal links (Di Capua et al., 2020). Seasonal forecast models are biased in their representation of the CGT, with typically a too weak CGT signal (Beverley et al., 2019). Therefore, seasonal forecasts miss an important source of predictability on intra-seasonal timescales, primarily in summer (Weisheimer and Palmer, 2014).

The WNPSM has been shown to influence precipitation anomalies over North America via its relation to the western Pacific–North America (WPNA) pattern (Wang et al., 2001). The WNPSM is shown to be related to surface pressure conditions over East Asia, with high pressure anomalies during years characterized by stronger WNPSM activity (Nitta, 1987). The WNPSM area also represents a genesis region for tropical cyclones in the North Pacific (Briegel and Frank, 1997). The WNPSM is weaker during the decaying phase of El Niño (Wang et al., 2001) and its related circulation anomalies provide a link from ENSO to the East Asian summer monsoon (Yim et al., 2008). Thus, in summary, the SAM appears particularly important for sub-seasonal variability over Eurasia, while the WNPSM is important for the Pacific–North American sector.

ENSO, operating at inter-annual timescales, might primarily influence the mid-latitude circulation via its effect on the SAM strength (Ding et al., 2011). During boreal summers preceding La Niña phases, a strengthening of the Walker circulation can enhance SAM rainfall, while El Niño phases often have an opposite effect (Joseph et al., 2011; Ju and Slingo, 1995; Terray et al., 2003; Wu et al., 2012). However, this relationship depends on the longitudinal position of the strongest El Niño-related sea surface temperature (SST) anomalies (Krishna Kumar et al., 2006) and potentially has weakened over recent decades (Chakraborty and Krishnamurti, 2003; Krishna Kumar et al., 1999; Srivastava et al., 2017; Xavier et al., 2007). At inter-annual timescales, anomalous tropical convection in the central–eastern Pacific related to ENSO has also been shown to affect mid-latitude circulation over the Euro-Atlantic sector as well as temperature and precipitation anomalies over Europe during boreal summer (O'Reilly et al., 2018). Trends in tropical SSTs play a crucial role in the interdecadal changes of this tropical–extratropical teleconnection (O'Reilly et al., 2019).

In general, a major challenge faced by teleconnection research is to understand the underlying physical processes and associated cause–effect relationships. Past observational studies have typically employed correlation analysis or linear regression techniques. Such analyses can however be dominated by spurious correlations and therefore can give only limited insight into cause–effect relationships. On the other hand, model-based studies can be affected by biases in the representation of circulation and precipitation characteristics (Beverley et al., 2019; Schubert et al., 2011; Weisheimer and Palmer, 2014), which can feed back on each other. Also, although perturbation and sensitivity experiments can point towards potential causal relationships, they do not necessarily reveal the causal links between tropical and mid-latitudinal features, since the uncovered relationship may not be the dominant one.

In recent years, several approaches have been applied to identify causal relationships in climate and atmospheric sciences (Runge et al., 2019b), ranging from Granger causality (McGraw and Barnes, 2018, 2020; Samarasinghe et al., 2019) to causal (Bayesian) graphical models (Ebert-Uphoff and Deng, 2012a, b; Horenko et al., 2017; Pearl, 2000) and conditional independence-based network discovery methods for time series (Runge et al., 2019a). These studies have shown the ability of causal discovery tools to improve our understanding of several atmospheric circulation interactions such as Arctic–mid-latitudes connections (McGraw and Barnes, 2020; Samarasinghe et al., 2019), synoptic-scale disturbances between boreal summer and boreal winter (Ebert-Uphoff and Deng, 2012a), and the relationship between ENSO and surface temperature in the American continent (McGraw and Barnes, 2018).

Here, we use a causal inference approach to study the relationships between the Northern Hemisphere mid-latitudes and the tropical belt during boreal summer at different intra-seasonal timescales. We apply a causal discovery approach making use of the so-called PCMCI (Peter and Clark algorithm combined with the Momentary Conditional Independence approach; see Sect, 2.3) method (Runge et al., 2019a) and then estimate physically interpretable causal links weights by (standardized) multivariate regression. The resulting weighted network representation of causal interdependencies is referred to as a Causal Effect Network (CEN) (Kretschmer et al., 2016). Expanding our understanding of the corresponding physical mechanisms has the potential to improve seasonal and sub-seasonal forecasts in boreal summer. The main advantage of causal discovery tools is that they can identify and remove spurious correlations (Runge et al., 2015b, 2019a; Runge, 2018) and thus provide insight into the potential causal relationships (McGraw and Barnes, 2018; Runge et al., 2014). Building upon this advanced methodology, we introduce a new concept called causal maps, visually highlighting causally related spatial structures. Finally, we assess the role of the ENSO background state on the identified causal relationships between the tropical belt and the mid-latitude circulation. The remainder of this paper is organized as follows: Sect. 2 presents the data and methods used in this analysis. Section 3 describes the results obtained by applying CEN and causal maps to the identified research questions. Section 4 provides a discussion of the obtained results in the context of the existing literature and finally, Sect. 5 presents a short summary and conclusions.

2 Data and methods

2.1 Data

In our analysis, we diagnose monsoon characteristics and Northern Hemisphere circulation features using outgoing longwave radiation (OLR) at the top of the atmosphere, geopotential height at 200 hPa (Z200) and 2 m surface temperature (T2m) data from the ERA-Interim Reanalysis (Dee et al., 2011) for the period 1979–2018 ($\mathrm{1.5}{}^{\circ }\times \mathrm{1.5}{}^{\circ }$). Strong tropical convection is characterized by high cloud tops and thus by low emission temperatures, which in turn correspond to low OLR values (Krishnan et al., 2000). In the tropical belt, OLR can be used as a proxy of convective activity and therefore rainfall. To select different ENSO phases, we use the monthly Niño3.4 index from NOAA (data are available at https://www.esrl.noaa.gov/psd/gcos_wgsp/Timeseries/Nino34/, last access: 12 October 2020), representing the central Pacific SST anomalies. El Niño and La Niña events are discerned by periods of December-to-February Niño3.4 index values larger than 0.5 ^∘C or smaller than −0.5 ^∘C, respectively. Then, we identify El Niño summers as those preceding the El Niño peak in winter and La Niña summers as those preceding the La Niña peak in winter. We use the Niño3.4 index since it has been shown to have a relatively strong connection to Indian monsoon rainfall (Krishna Kumar et al., 2006).

We also use the BSISO index as defined by Kikuchi et al. (2012) and Kikuchi and Wang (2010) (data are available at http://iprc.soest.hawaii.edu/users/kazuyosh/ISO_index/data/BSISO_25-90bpfil_pc.txt, last access: 12 October 2020) in order to describe the phase and amplitude of the BSISO characterizing the large-scale driver of active and break events over India. Causal discovery tool techniques require detrended anomalies centred at zero. Therefore, all data are linearly detrended and anomalies are calculated relative to an individual year's mean seasonal state by removing both the mean seasonal cycle and the year's mean seasonal state (i.e. the seasonal average from May to September, MJJAS) (Di Capua et al., 2020; Ding and Wang, 2007). Removing the year's mean seasonal state, and thus excluding the influence of inter-annual variations of the involved mechanisms, is essential to analyse intra-seasonal variability of atmospheric components that present a strong inter-annual variability, such as the SAM.

2.2 Maximum covariance analysis

To extract the dominant co-variability patterns reflecting interactions between mid-latitude circulation in the Northern Hemisphere and tropical convection at intra-seasonal timescales, we follow Ding et al. (2011) and apply maximum covariance analysis (MCA) to OLR fields (used as a proxy for convective activity) in the tropical belt (15^∘ S–30^∘ N, 0–360^∘ E) paired with Z200 fields in the northern mid-latitudes (25–75^∘ N, 0–360^∘ E).

MCA identifies the patterns that explain the greatest squared covariance between two different fields (Ding et al., 2011; Wiedermann et al., 2017) and ranks them according to their explained squared covariance fraction (SCF) (Wilks, 2011). Among the available correlation-based methods to highlight strong co-variability and reduce the dimensionality of a spatio-temporal dataset, MCA allows identification of patterns in pairs of variables that evolve simultaneously and may be causally related (via, e.g. dynamical coupling between multiple climatological fields). MCA detects patterns that can explain shared covariance, which cannot be achieved using other dimensionality reduction methods that consider individual variables separately, such as empirical orthogonal function (EOF) analysis. However, to provide a complete picture, we will also discuss the corresponding EOF patterns and the fraction of variance explained for comparison with our MCA results.

Each MCA mode provides two coupled (2D) spatial patterns (one for tropical OLR and one for mid-latitude Z200) and two associated (1D) time series (the time-dependent MCA scores or pattern amplitudes for both fields), describing the magnitude (prominence) and phase (sign) of those patterns for each time step. These (1D) time series are obtained by calculating the scalar product between each MCA spatial pattern (2D field) and the original spatial field of the associated variable at each time step as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }A={\mathit{u}}^T\mathbf{X},\text{(2)}& {\displaystyle }& {\displaystyle }B={\mathit{v}}^T\mathbf{Y},\end{array}$$

where A and B represent the two MCA scores for Z200 and OLR, X and Y are two matrices representing the Z200 and OLR fields, and u and v are the coupled patterns that maximize their covariance c, defined as follows:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}c& =\mathrm{cov}[A,B]=\mathrm{cov}\left[{\mathit{u}}^T\mathbf{X},{\mathit{v}}^T\mathbf{Y}\right]\\ & ={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}\left[{\mathit{u}}^T\mathbf{X}{\left({\mathit{v}}^T\mathbf{Y}\right)}^T\right]={\mathit{u}}^T{\mathbf{C}}_{xy}\mathit{v}\end{array},\end{array}$$

and

$$\begin{array}{}\text{(4)}& {\mathbf{C}}_{xy}={\displaystyle \frac{\mathrm{1}}{n-\mathrm{1}}}{\mathbf{XY}}^T,\end{array}$$

with n denoting the number of observation times.

Here, we select the first two MCA modes that represent the dominant patterns of co-variability between tropical convection and mid-latitude circulation, and we calculate time series for each MCA mode. These time series will be used as inputs for the causal discovery algorithm (see Sect. 2.3 and 2.4).

2.3 PCMCI and Causal Effect Networks

PCMCI is a causal discovery method based on the PC algorithm (named after its inventors Peter and Clark; see Spirtes et al., 2000) combined with the Momentary Conditional Independence approach (MCI, Runge et al., 2019). Given a set of univariate time series (called “actors”), PCMCI estimates their time series graph representing the conditional independencies among the time-lagged actors. In the context of the present work, actors are user-selected based on theoretical knowledge to represent either a specific component of the atmospheric circulation or surface conditions estimated with MCA (A, B) or an individual grid point time series C(lat, long). Assuming linear dependencies, PCMCI uses partial correlations to iteratively test conditional independencies and remove spurious links arising from autocorrelation effects, indirect links or common drivers. For example, if an actor Z drives X at lag −1 and Y at lag −2, then X and Y will be correlated, but the partial correlation ρ(X_τ−1, $Y_t\left|Z_{\mathit{\tau }-\mathrm{2}}\right)$ will be zero. PCMCI efficiently conducts partial correlation tests to identify which links cannot be explained by other time-lagged actors. Compared to the standard PC algorithm, PCMCI better deals with autocorrelation and high-dimensional sets of actors (Runge et al., 2019a). The output of PCMCI is a p value for each time-lagged causal link.

It is important to note that the term causal rests on specific assumptions (Runge, 2018; Spirtes et al., 2000), most importantly that it should be understood as “causal relative to the set of analysed actors”. Therefore, adding (or removing) an actor can alter the result of PCMCI, highlighting the importance of having an expert-guided hypothesis underlying the choice of the selected set of actors. In addition, using partial correlation for a conditional independence test implies further assumptions, such as the stationarity and linearity of the relationships. To control for multiple testing among the multiple grid locations in causal maps, we apply a false discovery rate (FDR) correction (Benjamini and Hochberg, 1995).

Based on the reconstructed network among the actor variables (at some significance level, α), we determine the causal parents as the incoming links to each actor (C(lat, long), A, B), which can come from the pasts of A, B or C(lat, long), i.e. $\mathit{\{}{A}_{\mathit{\tau }=-\mathrm{1}}$, $B_{\mathit{\tau }=-\mathrm{1}}$, C(lat, long)${}_{\mathit{\tau }=-\mathrm{1}}$, …, $A_{\mathit{\tau }=-{\mathit{\tau }}_{\max }}$, $B_{\mathit{\tau }=-{\mathit{\tau }}_{\max }}$, C(lat, long)${}_{\mathit{\tau }=-{\mathit{\tau }}_{\max }}\mathit{\}}$. In this analysis, A and B represent the two MCA scores obtained for a selected MCA mode, while C(lat, long) represents the grid point time series of a 2D field, e.g. T2m or Z200. In its first step, PCMCI iterates through partial correlations with increasing cardinality of conditions to remove the influence of common drivers and indirect links and estimate a preliminary set of parents. The first iteration of PC (cardinality 0) calculates the correlation between a selected time series, e.g. A_τ=0, and the past of any other available time series, $\mathit{\{}{A}_{\mathit{\tau }=-\mathrm{1}}$, $B_{\mathit{\tau }=-\mathrm{1}}$, C(lat, long)${}_{\mathit{\tau }=-\mathrm{1}}$, …, $A_{\mathit{\tau }=-{\mathit{\tau }}_{\max }}$, $B_{\mathit{\tau }=-{\mathit{\tau }}_{\max }}$, C(lat, long)${}_{\mathit{\tau }=-{\mathit{\tau }}_{\max }}\mathit{\}}$, including its own past $A_{\mathit{\tau }=-\mathrm{1},\mathrm{\dots },-{\mathit{\tau }}_{\max }}$. For illustration purposes, here we provide an example for C(lat, long), where ρ denotes the correlation and τ is the lag that is being used in the network (in this example, ${\mathit{\tau }}_{\max }=-\mathrm{2}$):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},A_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.32},p=\mathrm{0.01};\\ & \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},A_{\mathit{\tau }=-\mathrm{2}}\right)=\mathrm{0.13},p=\mathrm{0.1};\\ & \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},B_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.35},p=\mathrm{0.005};\\ & \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},B_{\mathit{\tau }=-\mathrm{2}}\right)=\mathrm{0.23},p=\mathrm{0.058};\\ & \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.41},p=\mathrm{0.01};\\ & \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{2}}\right)=-\mathrm{0.16},p=\mathrm{0.06}.\end{array}\end{array}$$

Applying a significance level α=0.05, only three actors are significantly correlated with C(lat, long) at the chosen time lag. These form the initial preliminary set of parents for C(lat, long) and are ordered by the strength of their correlation:

$$\begin{array}{}\text{(6)}& P_{C(\mathrm{lat},\mathrm{long})}^{\mathrm{0}}=\left\{C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}},B_{\mathit{\tau }=-\mathrm{1}},A_{\mathit{\tau }=-\mathrm{1}}\right\}.\end{array}$$

Next, partial correlations between C(lat, long) and each actor in $P_{C(\mathrm{lat},\mathrm{lon})}^{\mathrm{0}}$ are calculated by conditioning on the strongest preliminary parent:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}}\mathrm{|}{B}_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.35},p=\mathrm{0.02};\\ & \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},B_{\mathit{\tau }=-\mathrm{1}}\mathrm{|}C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.28},p=\mathrm{0.03};\\ & \mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},A_{\mathit{\tau }=-\mathrm{1}}\mathrm{|}C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.25},p=\mathrm{0.04}.\end{array}\end{array}$$

Parents with significant partial correlations will enter the second set of preliminary parents:

$$\begin{array}{}\text{(8)}& P_{C(\mathrm{lat},\mathrm{long})}^{\mathrm{1}}=\left\{C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}},B_{\mathit{\tau }=-\mathrm{1}},A_{\mathit{\tau }=-\mathrm{1}}\right\}.\end{array}$$

Next, the partial correlation is calculated conditioning on the two strongest parents:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathit{\rho }& \left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}}\mathrm{|}{B}_{\mathit{\tau }=-\mathrm{1}},A_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.31},\\ & p=\mathrm{0.03};\\ \mathit{\rho }& \left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},B_{\mathit{\tau }=-\mathrm{1}}\mathrm{|}C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}},A_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.23},\\ & p=\mathrm{0.04};\\ \mathit{\rho }& \left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},A_{\mathit{\tau }=-\mathrm{1}}\mathrm{|}C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}},B_{\mathit{\tau }=-\mathrm{1}}\right)=\mathrm{0.12},\\ & p=\mathrm{0.08}.\end{array}\end{array}$$

Since it is not possible to further increase the dimension of the condition set from the PC step, the preliminary parents converge as follows:

$$\begin{array}{}\text{(10)}& P_{C(\mathrm{lat},\mathrm{long})}^{\mathrm{2}}=\left\{C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}},B_{\mathit{\tau }=-\mathrm{1}}\right\}.\end{array}$$

By repeating this step for each variable (and for each longitude and latitude position), preliminary sets of parents are estimated. Let us assume that in our example we also obtain the following set:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& P_A^{\mathrm{3}}=\left\{B_{\mathit{\tau }=-\mathrm{1}},A_{\mathit{\tau }=-\mathrm{2}}\right\},\\ & P_B^{\mathrm{2}}=\left\{B_{\mathit{\tau }=-\mathrm{1}}\right\}.\end{array}\end{array}$$

In the MCI step, partial correlation is calculated again between each pair of actors (at different time lags) conditional on the above-estimated sets of preliminary parents, whereby both sets of parents are conditioned upon. To give one example, this would lead to the following result:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\mathit{\rho }& \left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},A_{\mathit{\tau }=-\mathrm{1}}\mathrm{|}{P}_{C(\mathrm{lat},\mathrm{long})}^{\mathrm{2}},P_A^{\mathrm{3}}\right)\\ & =\mathit{\rho }\left(C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=\mathrm{0}},A_{\mathit{\tau }=-\mathrm{1}}\mathrm{|}C(\mathrm{lat},\mathrm{long}{)}_{\mathit{\tau }=-\mathrm{1}}{B}_{\mathit{\tau }=-\mathrm{1}},\right.\\ & \left.B_{\mathit{\tau }=-\mathrm{2}},A_{\mathit{\tau }=-\mathrm{3}}\right)=\mathrm{0.1},p=\mathrm{0.3}.\end{array}\end{array}$$

Note that the parents of $A_{\mathit{\tau }=-\mathrm{1}}$ are shifted in time by $\mathit{\tau }=-\mathrm{1}$. After repeating Eq. (12) for each possible pair of actors and for time lags from 0 to τ_max, those parents that are significant in the MCI test will then form the final set of causal parents for each actor. We refer to Runge et al. (2019a) for a more detailed discussion and explanation of the algorithm design and extensive numerical experiments.

Finally, we estimate the Causal Effect Network (CEN) (Kretschmer et al., 2016; Runge et al., 2015a) among A, B and C(lat, long) by applying standardized multiple regression of each actor onto its causal parents identified via PCMCI, i.e. for Y∈ in A_t, B_t, C(lat, long)_t and the parents P:

$$\begin{array}{}\text{(13)}& Y_t=\sum _i{\mathit{\beta }}_i{X}_i+{\mathit{\eta }}_Y,\end{array}$$

where X_i∈P{Y}, i=1, …, N; i.e. the set of N parents of Y_t and η_Y is the residual of Y_t (the difference between the observed value Y_t and the value obtained by the linear regression on the causal parents $\sum _i{\mathit{\beta }}_i{X}_i$). Note that there can be different numbers N of parents for each actor.

Finally, the strength of a causal link $X_{t-\mathit{\tau }}\to Y_t$ is expressed in terms of the path coefficient β, which can be interpreted as the change in the expectation of Y_t (in units of its standard deviation – SD) induced by raising X_t−τ by 1 SD, while keeping all other parents of Y_t constant. Thus, for β=0.5, a change in a causal parent of 1 SD at lag −1 corresponds to a change 0.5 SD in the analysed actor at lag 0 (Runge et al., 2015a). The influence of an actor on itself is referred to as the autocorrelation path coefficient, which must not be confused with the Pearson autocorrelation. A detailed description of the PCMCI algorithm is available in Runge et al. (2019), while recent applications can be found in Kretschmer et al. (2016, 2018) and Di Capua et al. (2019, 2020).

Figure 1 Schematic explanation of causal maps. Panel (a) shows the correlation maps for A with C(lat, long) (left) and B with C(lat, long) (right). Panel (b) shows an example of a CEN constructed with A, B and C(lat, long) for a specific geographical position (identified with a diamond in the 2D maps). Panel (c) shows the corresponding causal maps, showing the path coefficients β from A to C, conditioned on B and all autocorrelations (left) and from B to C, conditioned on A and all autocorrelations (right). The “|” denotes the conditioned-out actor: A for the right-hand side and B for the left-hand side. See the text for further description.

2.4 Causal maps

To explore the causal effects that a specific actor has on a selected 3D (lat, long, time) atmospheric field, we introduce the concept of causal maps. Conceptually, causal maps are similar to correlation maps, as they show the spatial pattern of the relationship between a 3D climate dataset (covering two spatial dimensions plus time) and a 1D time series. However, instead of computing correlations between the time variations at each grid point and one additional time series, here we apply the PCMCI + CEN approach with actors consisting of the two MCA score time series (A, B) and each individual grid point time series (C(lat, long)). The causal map then shows the path coefficient β from one of the MCA scores (as one actor) to this grid point, conditioned on all remaining actors. For a set of two actor time series (A and B in Fig. 1) and one time-varying atmospheric field C, we can thus derive two causal maps: one from A to C(lat, long), conditioned on B and on the autocorrelation in all actors, and one from B to C(lat, long), conditioned on A and on all autocorrelation effects. Figure 1 provides an illustrative example of this type of analysis. Both correlation maps (Fig. 1a) indicate a positive value for the specific geographical location highlighted with the black diamond. The CEN constructed for A, B and C(lat, long) at this grid point is plotted in Fig. 1b and shows that only B is causally connected to C. The correlation between A and C is thus due to an indirect link via B (or to a common driver not included in the CEN). This is also seen in the causal maps showing the path coefficient β, which for the B→C(lat, long) link is positive (right-hand side) but is non-significant for the A→C(lat, long) link (left-hand side). In causal map visualization we can directly illustrate the effect of a specific actor on a global field (taking into account the influence of autocorrelation), indirect links and common driver effects due to other competing variables.

Here, we will derive causal maps using the time series obtained with MCA for modes 1 and 2 and Z200, OLR and T2m fields, both for the entire time period (1979–2018) and for two subsets depicting different ENSO phases, to assess how the ENSO background state influences the causal relationships. El Niño (La Niña) summers are defined as summers preceding the El Niño (La Niña) peak in boreal winter. We thus obtain 14 La Niña years and 13 El Niño years (see Table S1 in the Supplement for a list of corresponding years and Fig. S1 in the Supplement for the associated SST anomaly composites). Although the strongest SST anomalies related to the ENSO phase are found in winter, warm (cold) SST patterns related to El Niño (La Niña) phases are already clearly developed during the preceding summers.

Finally, to test the robustness of our causal maps to the choice of time period, we calculate causal maps for a range of sub-periods. In 10 trials we removed 10 % of the record (4 years). For ENSO-phase-dependent causal maps, we have shorter time series, and we thus remove 1 year in each trial, leaving a set of 14 causal maps for La Niña events and 13 causal maps for El Niño events. As a result, we obtain an ensemble of causal maps and apply the false discovery rate correction to p values of each single map. Then, both for the full period (1979–2018) and for El Niño and La Niña years separately, we masked out areas where less than 70 % of the trials indicated a significant causal link, giving an indication of the robustness of our findings and at the same time suppressing noise.

A summary of the abbreviation and variable used in this analysis can be found in Table A1, while the parameters used for the PCMCI algorithm are reported in the Supplement.

3 Results

3.1 Tropical–mid-latitude interactions: maximum covariance analysis

Figure 2a–d show the coupled patterns for the first two MCA modes between weekly tropical OLR and mid-latitude Z200. Figure 2e shows the associated time series of MCA scores for all four patterns (two for each MCA mode), which were obtained as explained in Sect. 2.2.

Figure 2 MCA of mid-latitude Z200 and tropical OLR at intra-seasonal timescales. Panels (a, b) show the first MCA mode for mid-latitude Z200 (25–75^∘ N) and tropical OLR (15^∘ S–30^∘ N), respectively, at the weekly timescale. The first MCA highlights the circumglobal teleconnection (CGT) pattern in the mid-latitudes and the South Asian monsoon (SAM) in the tropical belt. Panels (c, d) are the same as panels (a, b) but for the second MCA mode. This mode depicts the North Pacific High (NPH) in the mid-latitudes and the western North Pacific summer monsoon (WNPSM) in the tropical belt. The squared covariance fraction (SCF) of each MCA mode is given on top of the panels. Panel (e) shows the time series of MCA scores for the two MCA modes at a weekly timescale. Each MCA pattern has its own time series, i.e. one for tropical OLR and one for mid-latitude Z200 (note that different y axes are used).

The first two MCA modes highlight the two key patterns of boreal summer monsoonal activity in the tropics, along with the co-varying mid-latitude Z200 patterns. In both modes, the mid-latitudes are characterized by a zonally oriented circumglobal wave pattern with a wavenumber close to 5 (i.e. roughly five centres of action). However, the two wave patterns are phase-shifted and aligned with the longitudinal position of the strongest monsoonal convection in the tropics.

The first MCA mode explains 18 % of the squared covariance (squared covariance fraction, SCF) and shows a CGT-like wave-5 pattern in mid-latitude Z200. The Pearson correlation between the two time series of MCA scores for the first mode is r∼0.5. The spatial correlation with the weekly CGT pattern, as defined by Di Capua et al. (2020), is 0.52 (Fig. 2a). The CGT pattern also represents the second most important pattern in boreal summer mid-latitude circulation (Di Capua et al., 2020; Ding and Wang, 2005). This wave-5 pattern is linked to the South Asian monsoon (SAM) activity via its positive centre of action east of the Caspian Sea (see Fig. 2a). Applying MCA, we find that the CGT pattern co-varies with a band of enhanced tropical convective activity that extends from the Arabian Sea towards Southeast Asia, with a peak of convective activity over the Bay of Bengal (Fig. 2b) (Kang et al., 1999). Using OLR composites and the Kikuchi Boreal Summer Intra-seasonal Oscillation (BSISO) index, we explicitly show that the temporal evolution of SAM convective activity as defined in Fig. 2b at weekly timescales closely resembles the evolution of the BSISO (Goswami and Ajaya Mohan, 2001; Saha et al., 2012) (see Figs. S2 and S3 and further discussion in the Supplement). Therefore, we explicitly link the region of low OLR identified in Fig. 2b over the northern Indian Ocean and the Indian subcontinent to SAM activity as described in the literature. Note that we name each MCA pattern after a characteristic regional feature, but the analysis is applied to the larger geographical domains, as shown in Fig. 2.

The second mode of co-variability explains a SCF of 14 % and is characterized by a region of strong positive Z200 anomalies located at ∼45^∘ N, over the western North Pacific, directly to the west of the dateline (i.e. the most prominent centre of action of the mid-latitude wave). The Pearson correlation between the two time series of MCA scores for the first mode is r∼0.6. We will refer to this pattern as the North Pacific High (NPH) (Fig. 2c). The NPH is the summer counterpart of the North Pacific subtropical high, which characterizes boreal winter. During summer, this high-pressure region is displaced northward by the start of the monsoon season over the western Pacific Ocean and replaces the Aleutian Low (Lu, 2001; Riyu, 2002). The NPH is associated with a region of enhanced convection over the sub-tropical western North Pacific, related to the western North Pacific summer monsoon (WNPSM) convective activity (Fig. 2d) (Li and Wang, 2005; Nitta, 1987; Wang et al., 2001). The WNPSM core domain extends from 110–160^∘ E and 10–20^∘ N, while the boundary with the SAM is located over the South China Sea (Murakami and Matsumoto, 1994). The WNPSM is characterized by a late sudden onset (end of July) and a peak in rainfall activity during August and September, which is different from the SAM that features an earlier onset (in June) and peak rainfall activity during July–August.

We compare the patterns obtained with MCA with those obtained with EOF analysis of Z200 and OLR fields (see Fig. S4). We find that the closest match of the Z200 MCA mode 1 pattern is with Z200 EOF 2 (spatial correlation ∼ 0.8), while the closest match of Z200 MCA mode 2 is with EOF 1 (spatial correlation ∼ 0.6). OLR MCA mode 1 has the closest match with EOF 2 (spatial correlation ∼ 0.5), while OLR MCA mode 2 has the closest match with EOF 5 (spatial correlation ∼ 0.4). Thus, in general our MCA patterns also reflect the first two EOFs of Z200 and OLR, indicating that they explain an important fraction of the regional variability. Nevertheless, here we are interested in those patterns that can explain shared covariance, which cannot be achieved by using EOF analysis alone. Therefore, we use the MCA-defined patterns for the following part of the analysis.

We also investigate whether the obtained MCA patterns are sensitive to the choice of OLR in representing tropical convective activity. Using vertical velocity, another proxy for tropical convection where strong convective activity is represented by enhanced upward motion, shows qualitatively the same patterns as those in Fig. 2a–d (see Fig. S5). When velocity potential is used instead of OLR, the first MCA mode still closely resembles the OLR/Z200 MCA mode 1, while the second MCA mode only partly captures features in the western Indian Ocean (see Fig. S6).

Application of MCA to 4-weekly data gives nearly identical MCA patterns but with somewhat lower magnitude of the Z200 and OLR anomalies (see Fig. S7). In this case, we define both 4-weekly and weekly MCA scores by projecting 4-weekly MCA patterns onto 4-weekly and weekly data, respectively (see Fig. S7e and f). In this way, we check whether the analysis is robust given different definitions of the MCA patterns.

3.2 Influence of tropical–mid-latitude MCA modes on Northern Hemisphere circulation

To show how each MCA mode affects the circulation and surface conditions in the Northern Hemisphere, we calculate causal maps for the influence of SAM, CGT, WNPSM and NPH time series (as defined in Fig. 2e) on selected atmospheric fields in the Northern Hemisphere (15^∘ S–75^∘ N, 0–360^∘ E). Although we use ${\mathit{\tau }}_{\max }=-\mathrm{2}$ and τ_min=0, we plot only β values for lag −1 (week), as β values for longer time lags are mostly non-significant. This way the past behaviour of each actor, with potential confounding effects, is also accounted for in the corresponding grid point CEN. Note that we only show robust links as defined in Sect. 2.4, and the masks used to plot the results are shown in Figs. S8 and S9.

Figure 3 Influence of MCA mode 1 on Northern Hemisphere circulation. (a) Correlation map between the weekly SAM time series and the Z200 field. (b) Same as (a) but for the correlation between weekly CGT time series and the Z200 field. (c) Path coefficient β for the link SAM${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$ for a three-actor CEN built with SAM, CGT and Z200. Here, the “|” denotes the conditioned-out actor: CGT. (d) Same as (c) but for the link CGT${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$. The “|” denotes the conditioned-out actor: SAM. (e, g) Same as (c) but for the influence of SAM on OLR and T2m fields, respectively. (f, h) Same as (d) but for the influence of CGT on OLR and T2m fields, respectively. Only path coefficients β with p<0.05 (accounting for the effect of serial correlations) and the robustness mask (see Fig. S8) are shown. The dashed black line located at 30^∘ N shows the border between the tropical and the mid-latitude belt that separates OLR and Z200 analyses.

Figure 3 shows the causal maps for weekly Z200, OLR and T2m fields with SAM and CGT time series, including correlation maps for weekly Z200 fields with SAM and CGT time series. Referring to the schematic illustrated in Fig. 1 and following the PCMCI algorithm explanation (Sect. 2.3), here the A and B time series are represented by the SAM and CGT time series, respectively, while C(lat, long) is represented by Z200, OLR and T2m fields. In the mid-latitudes, the correlation map between Z200 and SAM (Fig. 3a) shows some similarities with the correlation map between Z200 and CGT (Fig. 3b), with negative correlation regions over central Europe and Scandinavia (Region 2 and Region 4) and over the eastern North Pacific and eastern Canada visible in both plots (regions 3 and 6). Both correlation maps also display a positive correlation over North Africa (Region 1), though with smaller values in the CGT plot. The causal map for the link SAM${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$ (after removing the effects of the CGT and of the past of both SAM and Z200) shows that the path coefficient β remains pronounced over northern India and North Africa (Region 1 in Fig. 3c), with values of β∼0.3–0.4. Interestingly, those regions disappear completely in the causal map for the link CGT${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$ (after removing the effects from SAM and of the past of both CGT and Z200). Thus, in this way, we can separate the signal coming from SAM convective activity from signals originating from the CGT pattern. Also, the causal maps in Fig. 3c and d indicate that the influence of SAM on Europe (negative path coefficients shown by Region 2 in Fig. 3c) and the North Pacific (seesaw of positive and negative path coefficients shown by Region 3 in Fig. 3c) is not mediated via the CGT. This influence is weaker than that over North Africa, with values of β∼0.1–0.3. In turn, the influence of SAM on other mid-latitude regions (over the North Atlantic, Region 4; over some parts of East Asia, Region 5; and Canada and the USA, Region 6; see Fig. 3d) is mediated via the CGT, with values of β∼0.1–0.3.

Figure 4 Influence of MCA mode 2 on Northern Hemisphere circulation. (a) Correlation between the weekly WNPSM time series and the Z200 field. (b) Same as (a) but for the correlation between weekly NPH time series and the Z200 field. (c) Path coefficient β for the link WNPSM${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$ in a three-actor CEN built with WNPSM, NPH and Z200. Here, the “|” denotes the conditioned-out actor: NPH. (d) Same as (c) but for the link NPH${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$. Here, the “|” denotes the conditioned-out actor: WNPSM. (e, g) Same as (c) but for the influence of WNPSM on OLR and T2m fields, respectively. (f, h) Same as (d) but for the influence of NPH on OLR and T2m fields, respectively. Only path coefficients β with p<0.05 (accounting for the effect of serial correlations) and the robustness mask (see Fig. S9) are shown. The dashed black line located at 30^∘ N shows the border between the tropical and the mid-latitude belt that separates OLR and Z200 analysis.

In the mid-latitudes, the causal maps for OLR and T2m (Fig. 3e–h) are largely consistent with those obtained for Z200, with negative β values for OLR, representing wet anomalies, overlapping with negative β values for T2m, representing colder anomalies, and dry anomalies (positive β values for OLR) overlapping with warm T2m (positive β values for T2m). The CGT influence is mostly concentrated in the mid-latitude regions, where the same Regions 4 to 6 identified in Fig. 3d can also be found in Fig. 3f and h. A significant and consistent negative causal effect of the CGT pattern on OLR values in the tropical regions can only be seen in a small area in the western Indian Ocean (Region 1 in Fig. 3f). Again, the OLR and T2m causal maps indicate that the SAM has a direct influence on North Africa and Europe as well as tropical Africa (Region 1 in Fig. 3e and g). Over the Indian peninsula and Indochina, strong convective motions (negative β values for OLR in Fig. 3e) are accompanied by colder temperature (negative β values for T2m in Fig. 3g), related to increased precipitation and, consequently, decreased surface temperatures during active SAM activity. The influence of SAM on the western North Pacific identified by Region 3 in Fig. 3c is also detected in OLR (Region 3 in Fig. 3e). Negative β values found over Region 2 in Fig. 3c are only slightly visible in OLR and T2m. However, we should also remind the reader that our causal maps show only the most robust links (see Sect. 2.4).

Figure 4 shows the same set of results but for the second MCA mode consisting of WNPSM- and NPH-pattern-related time series. Here, our A and B time series are represented by the WNPSM and NPH time series while C(lat, long) is again represented by either Z200, OLR or T2m fields. As expected, both correlation maps resemble the Z200 field of MCA 2 (Fig. 2c and d) with two characteristic features: an arch-shaped wave pattern in the North Pacific (Regions 7–9 in Fig. 4a and b), with a prominent positive correlation over the NPH region and over western North America, and two weaker centres of action over the Eurasian continent (Region 10 and 11 in Fig. 4b and Region 11 in Fig. 4a). The corresponding causal maps based on CENs are given in Fig. 4c (path coefficient β for the link WNPSM${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$) and Fig. 4d (for the link NPH${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$). If we compare the correlation maps (Fig. 4a and b) with the causal maps (Fig. 4c and d), we find great similarity in the spatial structures of the Z200 patterns over the North Pacific in both figures (Region 7 in Fig. 4c and d) with β∼0.1–0.3 for the influence of NPH on Z200 and β∼0.1–0.2 for the influence of WNPSM on Z200. These causal maps show that the influence of the WNPSM on Z200 (after removing the effect of the NPH) is confined to the North Pacific alone (Regions 7 and 8 in Fig. 4c). The causal effect of the NPH pattern on Z200 (after removing effects of WNPSM) shows the two most prominent regions, displaying a significant positive path coefficient β∼0.2–0.4 over the NPH region (Region 7 in Fig. 4d) and over the US West Coast (Region 9 in Fig. 4d). Region 7 in the Pacific sector coincides with that found for the WNPSM causal map (Fig. 4c). This suggests that the NPH is reinforced both by convective activity of the WNPSM and by the mid-latitude wave pattern localized in the North Pacific. Regions 10 and 11 found in Fig. 4a and b disappear in both Fig. 4c and d, showing that the correlation in these regions is mostly explained by Z200 activity in the mid-latitude itself (note that in the CEN we also condition on the past of Z200) or by other factors not considered in this analysis.

Next, we compute the causal maps for the influence of WNPSM and NPH on weekly OLR and T2m fields (Fig. 4e, f and g, h, respectively). The impact of the WNPSM on OLR and T2m fields is very weak, though it is possible to recognize some negative β∼0.1–0.2 over the Arabian Sea and over the WNPSM area (Region 8 in Fig. 4e). The impact of NPH on OLR and T2m causal maps shows some similarities with the correlation map shown in Fig. 4b. T2m and OLR show the strongest effect over North America (Region 9 in Fig. 4f and h) with β∼0.2–0.4, and this result is largely consistent with that obtained for Z200, with positive Z200 anomalies hinting at warm and dry weather in the mid-latitudes related to an active WNPSM. Over Eurasia, it is possible to recognize Regions 10 and 11 in both OLR and T2m causal maps but with smaller regions and lower β∼0.1–0.2 (Fig. 4f and h). Thus, these results highlight the importance of the NPH pattern in shaping surface temperatures across the northern mid-latitudes.

Using weekly MCA scores obtained from 4-weekly MCA patterns (Fig. S7) gives consistent results, showing that the analysis is robust when a different definition of the MCA pattern is chosen (see Figs. S10 and S11). Causal maps calculated for 4-weekly Z200 for both MCA 1 and MCA 2 show less significant results, likely due to the limited time series length (not shown).

3.3 The influence of ENSO on tropical–mid-latitude causal interactions

Next, we assess how the ENSO background state influences the causal relationships between mid-latitude and tropical patterns in boreal summer. We recalculate the causal maps for both MCA scores and Z200 fields, where A and B time series are represented by the SAM and CGT time series for the first MCA mode and by the WNPSM and NPH time series for the second MCA mode, while C(lat, long) is represented by Z200. As for Figs. 3 and 4, the robustness mask used to plot the results shown in Fig. 5 is shown in Fig. S12.

Figure 5 Causal maps and ENSO influence. Panel (a) shows the β for the link SAM${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$, a three-actor CEN built with SAM, CGT and Z200 during El Niño years. Here, the “|” denotes the conditioned-out actor: CGT. Panel (b) is the same as (a) but for the link CGT${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$. The “|” denotes the conditioned-out actor: SAM. Panel (c) is the same as (a) but for La Niña years. Panel (d) is the same as (c) but for the link CGT${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$. Panels (e, g) are the same as (a, c) but for the link WNPSM${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$ from a three-actor CEN built with WNPSM, NPH and Z200. Panels (f, h) are the same as (e, g) but for the link NPH${}_{\mathit{\tau }=-\mathrm{1}}\to \mathrm{Z}{\mathrm{200}}_{\mathit{\tau }=\mathrm{0}}$. Only path coefficients β with p<0.05 (accounting for the effect of serial correlations) and the robustness mask (see Fig. S12) are shown. The dashed black line located at 30^∘ N shows the border between the tropical and the mid-latitude belt which separates OLR and Z200 analysis.

In general, the strength and the sign of the patterns seen in the causal maps (Figs. 3 and 4) are not markedly affected by ENSO, though we can see higher absolute values of β∼0.1 in Fig. 5 with respect to Figs. 3 and 4. During La Niña years, the effect of SAM on the Sahara Desert intensifies, and its effects on the Tibetan Plateau and in the mid-latitudes are also more pronounced (Region 1, 2 and 6 in Fig. 5c). This is likely related to stronger SAM convective activity during La Niña summers. The region of negative causal effect of SAM on central Europe, is present only during La Niña summers (with a β∼0.2–0.3 and larger in area when compared to the 1979–2018 causal map in Fig. 3c), and it disappears during El Niño summers. This signal is possibly linked to the strong positive causal effect over the Sahara (Region 1). The region of positive causal effect over the Tibetan Plateau shows the same intensity as for the 1979–2018 period (β∼0.2–0.3) though it remains more confined over the Indian subcontinent and the Tibetan Plateau with respect to the 1979–2018 period (Region 1 in Fig. 5c). During La Niña, we also see an area of positive β values over North America (Region 6 in Fig. 5c) that is not present during El Niño.

During El Niño summers, the influence of the SAM is almost completely absent in the mid-latitude regions, with only one region of low β values over the eastern North Pacific still being present (Region 3 in Fig. 5a). The positive β values over the tropical Pacific found in Fig. 3c disappear during La Niña and only some residues of this region are seen during El Niño years (Region 8 in Fig. 5a).

In the western North Pacific, the most notable feature is the presence of both the WNPSM and NPH on the North Pacific only during El Niño summers (Fig. 5e and f). During those summers, the positive causal effect of the WNPSM over the western North Pacific (Region 7 and 8 in Fig. 5e) intensifies in magnitude (absolute β∼0.3–0.4) with respect to the 1979–2018 mean pattern (Fig. 4c), although the geographical extent of Region 7 shrinks. Over the western tropical Pacific, in correspondence with the La Niña warm pool, a region of positive causal effect is shown (Region 8 in Fig. 5e). These features disappear during La Niña summers. Thus, the second MCA mode (the WNPSM–NPH pair) has its strongest effect during El Niño summers, whereas the first MCA mode (SAM–CGT pair) is important during both La Niña and El Niño summers but with different spatial characteristics.

Calculating MCA pattern during different ENSO phases does not change the results in a qualitative way, although the order of the patterns is reversed in La Niña summers (see Fig. S13). Moreover, we have checked whether the distribution of the spatial correlation between each MCA mode and the respective Z200 and OLR fields changes during different ENSO phases and found that ENSO does not affect the frequency of each pattern in a significant way (see Fig. S14).

As for the 1979–2018 causal maps, when weekly MCA scores obtained from 4-weekly MCA patterns (Fig. S7) are used to provide ENSO-dependent causal maps, consistent results are obtained (see Fig. S15). Further analysis of possible physical explanations is provided in Sect. 4.

3.4 MCA causal interactions

Finally, we study the role of timescales on the causal interaction patterns presented above. We create CENs between the two time series of scores for each MCA mode, as identified in Figs. 2 and S7, and do so for weekly and 4-weekly data for the 1979–2018 period. Figure 6 plots the path coefficient β for two separate sets of CENs built for MCA mode 1 (SAM with CGT, Fig. 6a) and MCA mode 2 (WNPSM with NPH, Fig. 6b) for both 4-weekly and weekly timescales. As for the causal maps, we use ${\mathit{\tau }}_{\max }=-\mathrm{2}$ for weekly data and ${\mathit{\tau }}_{\max }=-\mathrm{1}$ for 4-weekly data. In both cases, τ_min=0.

Figure 6 Two-way causal link between tropical OLR and mid-latitude Z200. Shown is the path coefficient for pairs of MCA time series. The CGT is studied along with the SAM, while the NPH is analysed together with the WNPSM. Panel (a) shows the path coefficient β for the link SAM${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{CGT}}_{\mathit{\tau }=\mathrm{0}}$ over the 1979–2018 period (first row) and path coefficient β for the link CGT${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{SAM}}_{\mathit{\tau }=\mathrm{0}}$ (second row). The 4-weekly β values are shown in the left column, and weekly β values are shown in the right column. Panel (b) is the same as (a) but for WNPSM${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{NPH}}_{\mathit{\tau }=\mathrm{0}}$ and NPH${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{WNPSM}}_{\mathit{\tau }=\mathrm{0}}$ links, respectively. The β values with p<0.1 (0.05) are identified with one (two) asterisks.

At the 4-weekly timescale, the WNPSM–NPH pair does not show significant causal links (Fig. 6b). The SAM–CGT pair shows two fairly strong causal links with absolute values β∼0.3–0.4, though with different signs (Fig. 6a). The northward link, i.e. SAM${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{CGT}}_{\mathit{\tau }=\mathrm{0}}$, shows a positive β∼0.4: a 1 SD shift in the SAM leads to a ∼0.4 SD positive shift in the CGT 4 weeks later (Fig. 6a). The southward link, i.e. CGT${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{SAM}}_{\mathit{\tau }=\mathrm{0}}$, shows $\mathit{\beta }\sim -\mathrm{0.3}$, meaning that at this timescale a more intense CGT pattern leads to a weakening of the SAM pattern 4 weeks later (Fig. 6a).

At the weekly timescale, both the WNPSM${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{NPH}}_{\mathit{\tau }=\mathrm{0}}$ and the NPH${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{WNPSM}}_{\mathit{\tau }=\mathrm{0}}$ links show a β∼0.1–0.2, indicating that the two-way link has a similar magnitude in both southward and northward directions (Fig. 6b). At this timescale, the path coefficient β for the SAM${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{CGT}}_{\mathit{\tau }=\mathrm{0}}$ link is about a factor of 4 smaller than that for the 4-weekly timescale (Fig. 6a). The southward link, CGT${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{SAM}}_{\mathit{\tau }=\mathrm{0}}$, shows a positive β∼0.2 that is about twice as strong as the northward link (Fig. 6a). Thus, the influence of the SAM on the CGT pattern is weak (but present) at shorter (weekly) timescales but much stronger at longer (4-weekly) timescales.

Finally, we tested how the CENs change when the 4-weekly signal is removed from the weekly time series: each 4-weekly mean is removed from the four values of the corresponding weekly data (Fig. S16). This way, we attempt to isolate the dominant timescales of physical processes behind the different MCA patterns. This is similar in rationale to removing the effects of inter-annual variability before quantifying intra-seasonal variability. Results for the first MCA (SAM and CGT) show that the path coefficient β for the CGT${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{SAM}}_{\mathit{\tau }=\mathrm{0}}$ link remains almost unaffected (see Fig. S16a). This suggests that this southward link typically operates at weekly timescales (rather than 4-weekly), which is rather intuitive since mid-latitude variability dominates at synoptic timescales. In contrast, the path coefficient β for the northward link (SAM${}_{\mathit{\tau }=-\mathrm{1}}\to {\mathrm{CGT}}_{\mathit{\tau }=\mathrm{0}}$) becomes insignificant when the 4-weekly signal is removed from the weekly time series, suggesting that the influence from the tropics via the SAM pattern operates at longer, 4-weekly timescales. Removing the 4-weekly signal from the weekly time series for the second MCA (WNPSM and NPH) roughly halves the path coefficient β for both the northward and southward link (see Fig. S16b).

4 Discussion

In our analysis, we have found that the dominant patterns of interaction between the tropics and mid-latitudes remain qualitatively similar across different sub-seasonal timescales (weekly and 4-weekly averages) (Figs. 2 and S7). Two pairs of co-varying patterns are identified: (a) convective activity of the South Asian monsoon (SAM) paired with a mid-latitude circumglobal wave-5 resembling the circumglobal teleconnection (CGT) pattern and (b) convective activity of the western North Pacific summer monsoon (WNPSM) paired with a second wave-5 circumglobal wave pattern with its strongest action centre represented by the North Pacific High (NPH). This second mid-latitude wave pattern is phase-shifted with respect to the CGT pattern, to the longitudinal position of WNPSM monsoonal convection in the tropics. These patterns of sub-seasonal co-variability between the mid-latitudes and tropics during boreal summer are remarkably similar to those identified by Ding et al. (2011) for inter-annual timescales. This consistency across timescales (from weekly over monthly to inter-annual) suggests that the inter-annual patterns originate from a summation of the same patterns at sub-seasonal timescales. Still, the strength and sign of the causal interactions are timescale dependent. At longer timescales (from monthly to seasonal), slowly varying components such as tropical SST and associated regions of convective activity dominate tropical–mid-latitude interactions. Therefore, on these timescales the causal links from the tropics to mid-latitudes tend to be stronger. At shorter (weekly) timescales, in general a two-way positive feedback between the tropics and mid-latitudes is found, although strong variability in the mid-latitudes dominates over the tropical convection and thus the reversed southward pathways become stronger (Fig. 6). Moreover, we have introduced a novel visualization approach – termed causal maps – that can provide regionally specific information on the causal influence of a specific teleconnection source, and how this signal gets mediated. In this way, we identify mid-latitude regions and surface weather conditions that are influenced by tropical drivers by taking into account the linear influence of both MCA patterns together (for each MCA mode). The strongest causal effect of SAM convection is found over the Saharan region and depicts the monsoon–desert mechanism (Rodwell and Hoskins, 1996). Also important is the effect of SAM on the central Asian CGT centre of action (see Di Capua et al., 2020). The SAM also appears to directly influence geopotential heights in the North Pacific and central European surface temperatures 1 week ahead (Fig. 3c and e). The influence of the CGT pattern is stronger over the mid-latitude regions, nevertheless some influence on the Indian Ocean is detected (Fig. 3d and f), further supporting the results shown in Fig. 6. In the North Pacific there is a clear positive influence from the WNPSM towards the NPH patterns reflecting a Hadley cell-like circulation (Fig. 4d). The direct causal effect from the NPH on surface weather conditions is particularly strong in central North America, while its direct effect back on the tropics is weak (Fig. 4d, f and h). Thus, for MCA mode 2, the signal from the WNPSM towards the NPH is consistent both in simple CEN (Fig. 6) and causal maps (Fig. 4), while the direct influence of the NPH on the tropical belt is present but weaker and less robust (see Fig. S9).

In the tropical belt, the processes behind the identified MCA patterns are linked to strong convection related to monsoon activity. Though tropical convection is characterized by heavy precipitation with a typical duration of less than a day, the latent heat release can act as a Rossby wave source for up to 2 weeks after the initial forcing is removed (Branstator, 2014). Moreover, while individual convective events are short-lived, the regions of dominant convective activity in the tropics change on much longer timescales, such as in response to the BSISO (30–60 d). Thus, this finding could serve as a possible explanation for why the patterns identified at weekly and 4-weekly timescales show great similarity (see also the discussion surrounding Figs. S2 and S3). It appears reasonable to assume that the tropics operate at longer timescales, providing potential sources of predictability at seasonal-to-sub-seasonal (S2S) timescales. In contrast, mid-latitude circulation in summer is weaker than in winter and is characterized by circumglobal wave trains with typical timescales of about 1 or 2 weeks (Di Capua et al., 2020; Ding and Wang, 2007; Kornhuber et al., 2016).

On the western North Pacific side, our findings linking the WNPSM convective activity to the NPH, and in turn to a wave-5 circumglobal wave train that affects surface weather condition in the mid-latitudes, further support results from previous studies suggesting that convective activity related to this oceanic monsoon system can enhance the high pressure found in the North Pacific mid-latitudes and that this affects weather conditions in western North America (Chou et al., 2003; Wang et al., 2001). Four centres of action over northern and eastern Europe, central Asia, central North America, and the central North Pacific are identified in the T2m causal map (Fig. 4h). Eastern Europe and central North America were also identified to be teleconnected regions associated with a boreal summer wave-5 by Kornhuber et al. (2020), who highlighted the risk for potential bread-basket failures. Recent evidence indicates that land–atmosphere interactions and increased aridity in mid-latitude regions such as North America and Europe may constitute an enhancing mechanism for the amplification of circumglobal quasi-stationary Rossby wave events during boreal summer (Teng and Branstator, 2019). Therefore, our results support the importance of the role of Pacific forcing for this wave-5 circumglobal wave pattern. Although other mechanisms could also be relevant in exciting and maintaining this pattern, the link to the WNPSM convection may hold the potential to affect seasonal forecasts and climate risks, such as heat waves and simultaneous crop failures.

We have applied our causal map analyses to specific ENSO phases to assess the role of El Niño and La Niña in modulating the interactions between mid-latitude circulation and tropical convection in boreal summer. These analyses suggest that in general the ENSO phase does not change the qualitative nature of the causal relationships between different MCA patterns: the signs and strengths of the causal links are largely unaffected (see Fig. 5). Moreover, the same MCA patterns occur both in La Niña and El Niño summers, and their frequency is hardly affected (Figs. S13 and S14). Nevertheless, during La Niña summers, the effect of the SAM–CGT mode is reinforced over Europe, North Africa and the Indian subcontinent and reaches northward towards Canada, while during El Niño summers the effect of the SAM is mainly confined to the tropical belt. For the WNPSM-NPH pattern, a clear asymmetry between El Niño and La Niña summers is shown, with a stronger signal during El Niño years (Fig. 5e and f) that is absent during La Niña years. Although this effect is not very large, it is still significant. During La Niña summers, SAM exerts a stronger causal effect on the Tibetan High than over the entire 1979–2018 period, along with a reinforced monsoon–desert mechanism and a stronger effect on European circulation. This could be because under La Niña conditions the SAM circulation is supported by a favourable Walker circulation (Ju and Slingo, 1995; Terray et al., 2003). The same applies to the southward link: although ENSO does not alter the standardized causal effect from the CGT to SAM, a stronger CGT pattern in the mid-latitudes would have a stronger absolute effect on SAM activity at the weekly timescale. At inter-annual timescales, Ding et al. (2011) show that the SAM–CGT pair is strongest during the developing phase of ENSO. Therefore, our results further support the hypothesis that ENSO acts on the CGT pattern via its influence on SAM activity, in agreement with Ding et al. (2011). This finding is also in agreement with previous work showing that, at decadal timescales, the CGT and Silk Road patterns intensify under Pacific decadal oscillation (PDO) negative (i.e. La Niña-like) forcing (Stephan et al., 2019). During El Niño summers, the SAM shows a more prominent effect on the tropical Pacific. Nevertheless, since we condition on the effect of the CGT, we cannot exclude that this strong signal over the Niño-3.4 region may be due to an element outside our CEN. In the North Pacific, causal maps for different ENSO phases show stronger activity of both WNPSM and NPH links during El Niño summers, consistent with previous literature (Chou et al., 2003; Liu et al., 2016). During El Niño events, tropical convection shifts together with SST anomalies towards the central-eastern Pacific, which may favour WNPSM convective activity. In contrast, during La Niña summers convection is shifted towards the Maritime Continent and the western tropical Indian Ocean, reducing convective activity over the central Pacific and WNPSM region. A weaker WNPSM system may in turn be more prone to the influence of mid-latitude variability on the NPH.

Future projections describe an increase in monsoon precipitation associated with increasing global mean temperature and thermodynamic arguments (Menon et al., 2013; Turner and Annamalai, 2012). Quantifying teleconnections between the tropics and mid-latitudes is important in order to better understand and constrain future changes in boreal summer circulation, as uncertainty may arise due to changing connections to remote regions. While simulations show great uncertainty in the ENSO response to global warming (Cai et al., 2015; Chen et al., 2017a, b, 2015), observations show a La Niña-like warming trend in central-western Pacific SST (Kohyama et al., 2017; Mujumdar et al., 2012). A better understating of these teleconnections in observations can help to improve S2S forecasts. Verifying the existence and strength of causal teleconnections in forecast models could help diagnose the origin of model biases. For example, one could disentangle whether lower forecast skill (such as in the mid-latitude regions in summer) is related to local processes or to a misrepresentation of remote drivers. Beverley et al. (2019) showed that the CGT representation in seasonal forecasts is too weak. The CGT is important for predictability of summer extremes, and its relationship with the SAM may provide some information to improve predictability. Therefore, these methods could help answering the question “where do model biases come from?” and help in developing a physics-based bias correction. At the same time, CENs provide an encoded statistical predictive model, which can be used for actual forecasting (Di Capua et al., 2019; Kretschmer et al., 2017; Lehmann et al., 2020). Our analysis shows that at a 4-weekly timescale, the effect of SAM on the CGT pattern has a path coefficient β∼0.4, thus indicating some potential for predictability. Future studies should analyse how the causal links between these teleconnection patterns are reproduced in corresponding state-of-the-art S2S forecast and climate models, respectively.

Finally, it should not be forgotten that in the context of the present work, causal interpretation rests upon several assumptions, such as the causal Markov condition, faithfulness, causal sufficiency, stationarity of the causal links and assumptions about the dependence-type (Runge, 2018). These assumptions can be violated in a real system and it is important to be aware of the associated typical challenges for causal discovery in Earth system sciences (Runge et al., 2019a). Causal sufficiency requires that all relevant actors in a specific system are accounted for. Here, given the limited set of actors analysed, we cannot rule out that other excluded actors may act as important (common) drivers. Therefore, the obtained links can be considered causal only with respect to the specific set of actors used here. However, the absence of a link can still be interpreted as a likely indication that no direct physical connection among the respective variables exists (Runge, 2018). Moreover, we assume linear dependencies and stationarity for the detection of causal links. While linearity has been shown to be a useful assumption in previous work (Di Capua et al., 2020), monsoon dynamics behave partly non-linearly, and therefore our causal networks only capture some part of the underlying mechanisms by construction. Additionally, the SAM teleconnections might well behave in an non-stationary manner on decadal timescales (Di Capua et al., 2019; Robock et al., 2003). We therefore cannot exclude the possibility that (multi-)decadal oscillations such as the Pacific Decadal Oscillation may influence our results. However, the amount of reliable data is limited, and this prohibits the application of non-linear measures or the study of effects due to non-stationarity.

5 Conclusions

We have analysed the interdependencies and spatial effects of the two main MCA modes of co-variability between tropical convection and Northern Hemisphere mid-latitude circulation in boreal summer. The first MCA pair connects the circumglobal teleconnection (CGT) pattern in the mid-latitudes with the South Asian monsoon (SAM) convection, while the second MCA pair connects the western North Pacific summer monsoon (WNPSM) convection with a second circumglobal pattern related to the North Pacific High (NPH). These patterns appear qualitatively independent of the analysed timescales and emerge in weekly, 4-weekly and inter-annual analyses. The strength of the causal links is timescale dependent. In particular, the influence of SAM on CGT is strongest at the 4-weekly timescale, while the reversed link is stronger at weekly timescale. The patterns and sign of the standardized causal effect links are also not strongly affected by ENSO. During La Niña years the effect of the SAM on the mid-latitudes intensifies, while we find statistically significant links for the WNPSM effect on the mid-latitudes only for El Niño years. Moreover, the boreal summer intra-seasonal oscillation exerts strong control on the SAM convection at various lags.

Furthermore, we have introduced causal maps, a new application of the concept of Causal Effect Networks and have highlighted how this method can overcome limitations of correlation maps by removing spurious links. These causal maps further confirm our findings by showing a general positive two-way causal relationship between the dominant modes. Moreover, they highlight specific regions in the mid-latitudes that are particularly affected by the tropical modes (e.g. Eurasia, North America). These findings provide an improved understanding of the interactions between tropical convective activity and circumglobal wave trains that characterize mid-latitude circulation in boreal summer. This may help toward improving sub-seasonal forecasts and constraining future projections of boreal summer circulation. Further work shall assess whether these causal relationships are captured by general circulation models and whether this knowledge can be used to improve seasonal forecasts over the mid-latitudes.

Appendix A

Table A1 Abbreviations.