Nonlinear Forced Change and Nonergodicity: The Case of ENSO-Indian Monsoon and Global Precipitation Teleconnections

We study the forced response of the teleconnection between the El Niño–Southern Oscillation (ENSO) and the Indian summer monsoon (IM) in the Max Planck Institute Grand Ensemble, a set of Earth system ensemble simulations under historical and Representative Concentration Pathway (RCP) forcing. The forced response of the teleconnection, or a characteristic of it, is defined as the time dependence of a correlation coefficient evaluated over the ensemble. We consider the temporal variability of spatial averages and that with respect to dominant spatial modes in the sense of Maximal Covariance Analysis, Canonical Correlation Analysis and Empirical Orthogonal Function analysis across the ensemble. A further representation of the teleconnection that we define here takes the point of view of the predictability of the spatiotemporal variability of the Indian summer monsoon. We find that the strengthening of the ENSO-IM teleconnection is robustly or consistently featured in view of various teleconnection representations, whether sea surface temperature (SST) or sea level pressure (SLP) is used to characterize ENSO, and both in the historical period and under the RCP8.5 forcing scenario. It is found to be associated dominantly with the principal mode of ENSO variability. Concerning representations that involve an autonomous characterisation of the Pacific, in terms of a linear regression model, the main contributor to the strengthening is the regression coefficient, which can outcompete even a declining ENSO variability when it is represented by SLP. We also find that the forced change of the teleconnection is typically nonlinear by 1) formally rejecting the hypothesis that ergodicity holds, i.e., that expected values of temporal correlation coefficients with respect to the ensemble equal the ensemble-wise correlation coefficient itself, and also showing that 2) the trivial contributions of the forced changes in means and standard deviations are insignificant here. We also provide, in terms of the test statistics, global maps of the degree of nonlinearity/nonergodicity of the forced change of the teleconnection between local precipitation and ENSO.


Introduction
ENSO teleconnections are widely studied but their forced changes remain to be further explored and understood.Power and Delage [1] provide a multi-model assessment of ENSO-precipitation teleconnection changes based on the CMIP5 archive.They consider, in particular, ENSO-driven precipitation anomalies in tropical regions around the globe, and assess them jointly with changes of mean precipitation.For example, they find that, as a combination, wet DJF anomalies are "strengthened" in the Tibetan plateau both by increasing mean precipitation and wet ENSO teleconnection strength; or, as another type of combination, an increasing JJA dry ENSO teleconnection strength is counteracted by a wetter mean state in Southeast Asia.Haszpra et al. [2] have evaluated only the trend in ENSO-precipitation teleconnection strength changes, however, not in a model ensemble but the so-called "single model initial condition (large) ensemble" (SMILE) CESM1-LE [3].Working with a SMILE has the advantages that the forced response is correctly represented [4,5,6] in that model at least, and seeking a physical interpretation of changes is not faced with confusion at the outset, even if the physics is misrepresented by that one model.The forced response is the time evolution of some ensemble-wise statistics, or, most generically, that of the probability measure carried by the climate snapshot attractor [4,5,6].The ensemble-wise statistics evaluated over the converged ensemble [7] is in a one-to-one correspondence with the external forcing of the nonautonomous dynamical system, hence the term 'forced response'.Without convergence, or considering a single realisation, which may even correspond to observations, the state depends also on initial conditions, so that forced changes cannot be precisely disentangled from changes brought about by internal variability.In an attempt of doing this nevertheless, concerning single realisations, the standard practice is that a trend, linear or not [8], is simply identified as a forced change, before possibly "anomalies", i.e., (what is assumed to be) the internal variability is analysed.Identifying the principal component (PC) of the first empirical orthogonal function (EOF1) [9] obtained without detrending with the forced response [10] is still somewhat arbitrary.It certainly leads to biases [11].Yet, even with SMILEs available, it is very common to see that authors evaluate temporal statistics first -unnecessarily involving a subjective factor -and ensemble statistics afterwards [12,13].Otherwise, the biases are certainly controlled by the magnitude of the climate change signal relative to the intensity of internal variability -a kind of signal-to-noise ratio.Computing the PC1 of EOF1 without detrending, as mentioned just above, is in fact meant to maximise the signal at least, and the EOF1 is referred to [14] as a "fingerprint" of the forced change -the spatial pattern that is supposed to be distorted by internal variability, in any single realisation, the least.See [15] that is concerned rather with the minimisation of the "noise".
It is not pertinent to talk about "advantages" of temporal or ensemble methods over one another, because no situation arises when we need to decide between them.When an ensemble is available, the correct, conceptually sound ensemble method is to be applied if the forced response (including that of the internal variability) is concerned.However, certainly we can make statements only about the given model due to model errors.Therefore, we should speak only about limitations in the respective situations of analysing observational or model ensemble data.
Concerning teleconnections, in particular, the forced response can be identified by the time evolution of the ensemble-wise Pearson correlation coefficient [16,17] (in a simple linear approach) between some anomaly quantities.The anomaly can correspond to simply a spatial (mean of a temporal) mean [18], but also the PC of an EOF concerning variability across the ensemble, called a snapshot EOF (SEOF) [19].Haszpra et al. [2] take the latter approach; however, it can be extended to obtaining anomalies observing the "mutual variability", such as in the sense of Maximal Covariance Analysis (MCA) [9] or Canonical Correlation Analysis (CCA) [20,9].With an interest in a teleconnection and its forced response, MCA and CCA -just like EOF analysis -can also be pursued concerning the variability across the ensemble, whereby we can refer to these methods as SMCA and SCCA.This is perhaps best suited to teleconnection analyses concerning two extended, possibly disconnected, regions.See an application of SMCA to studying the forced change of the coupling of JJA 200-hPa geopotential height and September sea ice concentration in [21].
An important example of this would be the relationship of the summer monsoon precipitation on the Indian subcontinent with the ENSO phenomenon which extends over the whole of the Equatorial Pacific [22,23,24].In a previous publication [18] we examined the forced response of the ENSO-Indian monsoon teleconnection in the Max Planck Institute Grand Ensemble (MPI-GE) [25] representing the monsoon by the average JJAS precipitation over India and the ENSO by either the gridpoint-and SLP-based SOI, or the areal mean of the SST in some extended area in the Equatorial Pacific.Standard choices for the latter are the Niño3, Niño4, and Niño3.4 regions.For the first time we established, via a formal hypothesis test, that standard representations of the teleconnection in the sense of a correlation coefficient between some of the above-mentioned indices strengthened in this model, both in the historical period and under the strong RCP8.5 forcing scenario, although not necessarily monotonically.In fact, using the SOI it was the historical period only when an increase could be detected, and, using the Niño3 index, it was rather the RCP8.5 scenario under which an increase could be detected.This raises the question whether the discrepancy is down more to (1) the physical difference between the SOI and Niño3 indices, or (2) rather that -one based on two distant gridpoints while the other on a limited region of the Equatorial Pacific -they "take somewhat different slices" of the whole ENSO phenomenon, which can matter when even spatial characteristics of ENSO have a forced response beside possibly the "temporal".In fact, using the Niño3.4index, instead of Niño3, the result is closer to that using the SOI, or, even more so using a box-SOI, even to the extent of detecting the same nonmonotonicity: a temporary weakening around the turn of the millennium.In addition, both with Niño3.4 and the box-SOI we can detect changes both under the historical and scenario forcing.On the one hand (I), this seems to rule out that using SST versus SLP makes much difference.We will see below that this is actually not precisely the case.On the other hand (II), it is still a question whether spatial characteristics of ENSO, in terms of SEOFs, or the whole teleconnection, in terms of SMCA or SCCA, would also change, or only the variance as reflected in the magnitudes of the PCs belonging to either of the said instantaneous/snapshot spatial modes.These are the two basic questions that we set out here to investigate.
Furthermore, this study identifies, in terms of a regression model, controlling factors driving the forced change of the ENSO-IM teleconnection strength as being: the ENSO variability, the coupling (regression coefficient), and the intensity of other influences -which can be viewed as noise.We find that the changes of the ENSO-IM teleconnection are not driven only or dominantly by the change of ENSO variability.Nevertheless, it imprints its nonlinearity onto the teleconnection change, which could be an important source of biases as follows.We observe in the MPI-GE a nonmonotonic change in the ENSO variability: after a seemingly monotonic change, a decline follows in the second half of the 21.c. under RCP8.5.Such a feature is absent in the CESM1-LE [2].This nonmonotonicity and possibly others, as already mentioned above, constituting a nonlinear change, prompt that the system should be nonergodic also with respect to correlations, i.e., biases should exist [11] in the temporal correlation coefficient evaluated in, say, multi-decadal time windows.Such a bias is something to bear in mind beside the ample statistical fluctuations of finite-size statistics even under a stationary climate [18].Given that the ensemble size is finite and not so large concerning teleconnections, we develop here a statistical test whereby we can detect nonergodicity, and, subsequently, map out regions of the world where such a nonergodicity can be detected in the MPI-GE in the context of the relationship of local precipitation with ENSO.
The rest of the paper is organised as follows.In Sec. 2 we provide details about our methodology of analysing the forced response of teleconnections, such as the way we pursue EOF, MCA, CCA, as well as the idea of decomposing the change of the teleconnection strength in terms of a simple regression model.In Sec. 3 we provide results both on spatial and temporal aspects of the forced change of the ENSO-IM teleconnections.In Sec. 4 we outline our method of detecting nonergodicity and map out the world with respect to the degree of nonergodicity concerning the synoptic relationship of the local JJA precipitation with ENSO.In Sec. 5 we discuss and summarise our results.
2 Data and methods

Data
The analysis is restricted to the use of the MPI-GE data [25], that is, to a particular Earth system model (ESM), the MPI-ESM.The quoted presentation paper of the data set outlines main features of the state-of-the-art MPI-ESM.
To characterise the ENSO we make use of the 2D fields of either the Sea Level Pressure (SLP) or the Sea Surface Temperature (SST).The Indian summer monsoon rain (ISMR) was calculated from the total (convective and large scale) precipitation rate.The variable codes for SLP, SST and total precipitation rate are 134, 12 and 4, respectively.Monthly mean SLP is accessible publicly at https://esgf-data.dkrz.de/projects/mpi-ge/under variable name 'psl'.The SST can be derived from the top layer of the 3D potential temperature field, whose variable is 'thetao'.Instead of the total precipitation rate [km/s], one can use the precipitation flux [kg/s/m 2 ], whose variable name is 'pr'.

Representations of the ENSO-IM teleconnection
The scalar quantity of the ISMR averaged in time over the JJAS monsoon season and in space over India is the same as that calculated in [18], in order to secure a correspondence with the so-called AISMR rain gauge data set [26] (excluding some states of India; see their Fig. 1).Scalar quantities to represent ENSO variability are also the same as in [18], as follows.E).We will refer to the average SSTs under 1. and 2. as the Niño3 and Niño3.4 indices, respectively, and to the pressure differences under 3. and 4. as the SOI and box-SOI indices (SOI being short for the Southern Oscillation Index), respectively.See a discussion on this in Sec.4.a of [18].
We correlate the "near-synoptic" JJA mean of any of Niño3, Niño3.4,SOI, box-SOI with the JJAS (model) AISMR.The correlation coefficient r is evaluated across the ensemble, as done in [18], which results in annual time series.These time series provide representations of the forced response of the ENSO-IM teleconnection.However, the finite ensemble size entails large statistical error fluctuations, which badly mask the forced response signal.That is, the estimate r(t) rather poorly represents the unknown true r(t) signal.Therefore, we aspire only to detect a change in certain time periods.For this, we employ the Mann-Kendall test [27] in the same way as done in [18], aiming to reject the hypothesis of stationarity against the alternative hypothesis of a monotonic trend masked by serially uncorrelated noise constituted by realisations of iid rv's (independent identically distributed random variables).
We will evaluate correlations of the AISMR also with the PCs (which are scalars) belonging to EOFs.PCs and EOFs can be also composed with respect to (wrt.) the ensemble, as a "snapshot" (never mind the seasonal means), hence the name SEOF, as proposed by [2,19] 1 .The SEOFs are evaluated here wrt. the Equatorial Pacific box (30 • S, 30 • N, 150 • E, 295 • E).Traditional EOFs are decomposing spatio-temporal variability as a sum of standing waves, i.e., time-independent spatial patterns modulated by arbitrary (but uncorrelated) temporal signals, whereby the spatial patterns are given by the EOFs, and the temporal signals by the PCs.Snapshot EOFs are the same but different time steps are replaced by different ensemble members (i.e., realizations of the dynamics).EOFs are orthogonal, and each new EOF is defined such that the corresponding PC has the largest possible variance, leading to an eigen-problem.The variances of the PCs are singular values belonging to the EOFs being the singular vectors.We find them by using Matlab's svd as done in [31].The ordering of the singular values wrt.magnitude provides a natural ordering of the EOFs, which are denoted as EOF1, EOF2, etc., and like-wise we write PC1, PC2, etc. Correlating e.g. the PC1s of either the SLP or the SST field (in the same box) with AISMR have more in common than the representations given by e.g. the pairs of Niño3-AISMR and SOI-AISMR, because Niño3 and SOI do not derive from the same mathematical definition of dominance2 ; and one could think that more similar representations yield more similar results.Therefore, firstly, we wish to see if we can establish a robust representation of the ENSO-IM teleconnection by having possibly a better match between the forced response signals belonging to the SLP and SST than between those using the classic indices 1.-4.Secondly, the teleconnection strength can change possibly because of a shift in the "centre of action" of the ENSO phenomenon, something that can impact on the teleconnection representations by 1.-4.differently.Thirdly, including higher EOF modes can reveal more predictive power of the full e.g.SST field compared with e.g.Niño3 alone 3 .In this regard we are interested in whether the correlations to do with the higher modes are more susceptible to the anthropogenic forcing.We will only show results for the second modes.
We can pursue the three points of inquiry above also in terms of MCA and CCA.MCA and CCA are similar to the EOF analysis, but consider two fields, and find the orthogonal modes for which the covariance and the correlation, respectively, is maximal between the corresponding PCs of the two fields.We carry out these analyses by applying Matlab's svd and canoncorr, respectively.These methods will allow us to take into account a change in spatial characteristics on the side of the IM too.We will use the box of (5 • S, 40 • N, 65 • E, 100 • E) on the side of the IM.We evaluate the correlation coefficient between the PC1s of the respective 1st SMCA/SCCA modes as well as the PC2s.That is, the correlations of these PC2s are not analogous to those of the AISMR and PC2 of EOF2.We anticipate that the SMCA and SCCA lead to higher correlations than the SEOF analysis, stemming from their very definition.In this regard, as the fourth point of inquiry, we want to see if higher correlations are more susceptible to anthropogenic forcing.Instead of the complete field, the SCCA is performed on the PCs corresponding to the first 4 SEOFs, since it is ill-defined on the full fields.We will speak about SCCA modes as the linear combination of the SEOF modes in terms of the coefficients defining the so-called canonical variables.Indeed, thanks to the orthogonality of the SEOFs, the projection of the full fields onto these modes yield the SCCA PCs whose correlations were maximised by the SCCA.

Decomposition of the forced change of the teleconnection strength
Instead of trying to find physical mechanisms responsible for a particular forced change signal in the teleconnection strength, here we pursue only a statistical attribution of that.In terms of a linear regression model, Ψ = aΦ + ξ, underpinning the correlation coefficient, r, we can attribute changes in the teleconnection strength to three factors.We can identify these factors by considering the formula: These factors are: • ENSO-IM "coupling" a being the regression coefficient; Note that, like r, all σ Φ and σ ξ are defined in terms of the variability wrt. the ensemble, for every year separately.That is, like r(t), we have time series σ Φ (t) and σ ξ (t), and also a(t).Based on these, in a simple way, we can say that e.g. a change in r(t) in a time period is due to a change in σ Φ (t) if in that period no change is seen wrt.a(t) and σ ξ (t).Or, perhaps we can also say that if σ ξ /a shows no change.In fact, it has been suggested to us that perhaps the strengthening of the ENSO-IM teleconnection is due in this model to an increasing ENSO variability.
Although apriori this possibility is perhaps hard to exclude, we note that it would be rather specific, and could possibly be seen accidental, to the ENSO-IM teleconnection given that ENSO-precipitation teleconnections can either strengthen or weaken in different places on Earth [1].This is shown in Fig. 5 of [2] for the CESM1-LE; and the MPI-GE data shows a rather similar picture as seen in our Fig.S1 in the Supplementary Material.
In order to check the hypothesis of the dominant role of ENSO variability, one can evaluate σ Φ directly (i.e., compute the standard deviation of Φ over the ensemble), on the one hand, and σ ξ /a by simply inverting Eq. ( 1), on the other.(The only issue with the latter is that the sign of a remains indeterminate.)Anticipating that σ ξ /a is not constant in time, we can evaluate a by, first, directly evaluating the IM variability σ Ψ , and, subsequently, using the textbook formula: (We can, of course, evaluate a directly by linear regression, but calculating the standard deviation is easier, and we already have r on hand.)In turn, having now also a on hand, σ ξ can be obtained by inverting Eq. ( 1).From eq. ( 1) we see the following: if in a time period [t 1 , t 2 ] the ENSO variability increases as σ Φ (t 2 ) = βσ Φ (t 1 ), β > 1, and we have also a decrease σ ξ (t 2 )/a(t 2 ) = ασ ξ (t 1 )/a(t 1 ), α < 1, then αβ > 1 would mean that the ENSO variability increase, in comparison with the σ ξ /a decrease, has a larger effect on the increase of r(t) in that period.Given the very noisy time series, the appropriate approach would be a statistical test attempting to reject the null-hypothesis of αβ = 1.However, it is not clear to us how this can be done.As for a preliminary analysis, we simply evaluate αβ in all possible time periods [t 1 , t 2 ] assuming linear changes of σ Φ (t) and σ ξ (t)/a(t).That is, α and β are determined from the respective linear regression model fits.
3 The forced response of the ENSO-IM teleconnection

Spatial aspects
We start by inspecting the forced response of spatial characteristics of the ENSO and IM variability and those of their teleconnection.Fig. 1 shows the first modes of the SEOF, SMCA, SCCA analyses on the side of the ENSO, based on the SST, smoothed by averaging within four subsequent 50-year periods starting from 1900.The temporal averaging brings it really out that hardly any change is visible in the EOF1 and the MCA1 mode, even upon the strong RCP8.5 forcing.Carreric et al. [13] suggests that the maximum of the EOF1 mode shifts to the east in the CESM1-LE by showing that some centre (C) mode [32] shifts to the east more in comparison with the westward shift of some east (E) mode 4 .This does not appear to be so here, although the error of the results (with respect to any year within a given period) would be difficult to be estimated.Unlike the EOF1 and the MCA1 mode (which are very similar to each other), the CCA1 mode does undergo a considerable change by the late 21st century.It is even more interesting to observe this change to have a different sign in the middle of the ENSO domain (highlighted by the Niño3 and Niño3.4 boxes) than the minor change in the EOF1 and the MCA1 mode.Since the MCA1 mode is supposed to mainly reflect the ENSO variability (see later), unlike the CCA1 mode, which truly captures the variability relevant for the teleconnection, this discrepancy presumably means (taking into account the patterns in the reference period) that important areas for ENSO itself and for its teleconnection with the Indian monsoon get more decoupled.In particular, the middle of the basin gains importance for ENSO but loses importance for the teleconnection.The supplementary figure Fig. S2 suggests the same picture for the SLP.
Although each panel of Fig. 1 washes together 50 separate entities, the mostly monotonic emergence of the changes in consecutive rows suggests that we can see meaningful signals.Otherwise, the supplementary video https://youtu.be/RdHjbjZxaZQshows ample year-to-year change in the form of a flickering.Given that the forcing is rather smooth and gradual, this is clearly just a finite-ensemble-size statistical error fluctuation.
Interestingly, the second modes of ENSO variability, both in terms of SST (Fig. S3) and SLP (Fig. S4), are largely unchanged, too.We note, however, that the obtained temporal-mean second modes are less reliable than the first modes for the following reason.The sign of a mode is arbitrary, and we need to choose a sign convention, according to which we flip the modes returned by the svd Matlab routine in years when needed.To this end we pick a reference year, and calculate the spatial correlation of the modes in all years with the mode in the reference year.When the sign of the correlation is negative, then we flip the sign/orientation of the mode.However, when the statistical error corrupting the modes is too large, the pattern matching might not be possible.It seems to be the case that the pattern matching is feasible for the first modes, but fails in certain years for the second modes.See the supplementary videos https://youtu.be/bmbPXAtq_p4 and https://youtu.be/4pWPf2oUE-I.
We can achieve a smoothing of the temporal evolution of the modes by pooling snapshot (annual) anomalies from a temporal window for the evaluation of the modes.That is, we mix temporal and ensemble-wise variability.
We emphasize, however, that the anomalies are obtained the same way as done without smoothing: simply by subtracting the ensemble-mean belonging to single years, so that the only subjective element is the choice of the window length.The supplementary videos (https://youtu.be/bmbPXAtq_p4https://youtu.be/R78ycuvOI6c https://youtu.be/4pWPf2oUE-Ihttps://youtu.be/RdHjbjZxaZQ)show that the effect of smoothing by a 5-year window is well visible, but would still allow for modes to be far "out of shape".Considering the comparison of 50-year periods, however, it is not clear to us which way would smaller errors be incurred -whether by averaging 50 annual "snapshot" modes or by pooling anomaly data from the full 50-year period -because there is no reason to think that the finite-ensemble-size estimate is unbiased 5 .It is certainly computationally far less intensive by pooling the data, because then the expensive SVD needs to be performed only once.
On the side of the IM, Fig. 2 evidences a change in spatial characteristics, namely, a shift of some "teleconnection hotspots" in Northern India and over the sea to North-Northwest and to North, respectively, at least in the last 50 years.Note that the EOFs, the MCA modes and the CCA modes, including their changes, are very similar.The biggest dissimilarity is again shown by the CCA mode as the relative weakness of a peak over the sea West of the land, but one may reasonably argue that the teleconnection with ENSO (captured by CCA1 and MCA1) is strongly bound to the strongest Indian variability (captured by EOF1 and MCA1).In view of this, the coincidence of EOF1 and the MCA1 mode for ENSO may be attributed to the dominating variability of ENSO in covariance (as opposed to correlations, which determine the CCA1 mode).In Fig. 2, we show also the correlation map and its changes concerning the relationship of the gridpoint-wise precipitation with the PC1 of EOF1 of ENSO, and find similar patterns to those of the first SEOF/SMCA/SCCA modes, which further supports the close relationship between the Indian precipitation variability and the ENSO-influence.Note that many areas of the strongest changes fall outside the AISMR region.In particular, most of the sea loses very much of its relevance for the IM variability and for its teleconnection with ENSO (a positive change results in getting closer to zero in these areas of negative sign).Like the first modes, so do the second modes show changes in the last 50 years (Fig. S6).This further confirms that the teleconnection with the ENSO region is geographically rearranged over India.Using SLP instead of SST on the side of the ENSO, the spatial patterns of precipitation over India are largely unchanged for both the first and second modes (Figs.S5, S7).These features can be seen also in the supplementary videos https://youtu.be/bmbPXAtq_p4 and https://youtu.be/R78ycuvOI6c.
As a curiosity, we further elaborate on the similarity and the differences between different modes.As for the first ENSO modes, CCA1 gives prominence to the variability in the westernmost part of the selected box in comparison with EOF1 and MCA1.Presumably this is in large part for the intuitive reason that for a maximal correlation (not covariance) it is what happens in the Pacific the closest to India that matters the most.For the SST, this comparison persists to an extent for the second EOF, MCA, CCA modes, with somewhat less similarity of EOF2 and MCA2.The SLP shows a similar behaviour.As for the IM, MCA2 is an outlier, but matching is very good in all other cases.

The forced evolution of correlation coefficients and the drivers of their change
The time evolution of the correlation coefficients is given in the top row of Fig. 3 for the SST-based representations of ENSO.This time evolution is basically a long-term increase for the Niño indices (reproduced from [18]) and for the first modes.The second modes behave differently: their correlation coefficients are practically constant, so that they really cannot have higher than secondary importance for changes in the teleconnection strength.
For the Niño indices and the first modes, the details of the detectable time evolution are presented in Fig. 4 by temporal maps of the Mann-Kendall test statistic for stationarity, following [18].The general increase in teleconnection strength is detected by every quantifier of ENSO and the IM.Although a change is hardly detected within the 21st c. part alone (upper right triangle), the gradual darkening of the shades in the upper left box suggests a continued strengthening of the teleconnection.Note that the radiative forcing is much stronger in the 21st c. than in the 20th c., so that any strengthening in the 21st c. must be nonlinearly slowed down in comparison with the 20th c.
The further details also mostly match in all representations utilizing the AISMR (Figs. 4(a)-(c)), although an early-21st-c.drop does not surpass the threshold for 5% significance for Niño3, and the late-21st-c.increase is rather moderate for Niño3.4.This feature is even more pronounced in the EOF1-EOF1 and the MCA1 representation (Figs.4(d)-(e)), whereas the CCA1 representation comes with a very prominent gradual increase extending to all of the 21st c. (Fig. 4(f)).Interestingly, all representations based on principal modes on the IM side place the slight drop at the turn of the 20th-21st centuries considerably earlier than the AISMR-based representations.These differences may be linked to the shape and rearrangement of principal modes as described in Section 3.1.The slight dichotomy between Figs. 4(a)-(c) and (d)-(f) may be explained by pattern rearrangements on the IM side dominantly falling outside the AISMR region but being similar for the different principal modes.Understanding the details would require better temporal resolution than in Fig. 2, but using AISMR clearly favours the detection of the long-term increase in teleconnection strength, which is naturally explained by non-AISMR areas, especially the sea, losing importance.On the other hand, the discrepancy between Figs. 4(d)-(e) and (f) is presumably linked to the slightly increasing importance of the middle part of the ENSO basin from the point of view of ENSO variability (reflected by EOF1, and also MCA1, as mentioned) together with the strong loss of importance of this area for the teleconnection (as CCA1 shows; see Fig. 1).It is thus not a surprise that using EOF1 or MCA1 result in lower susceptibility for detecting an increasing teleconnection strength (in comparison with Fig. 4(f)).By chance, the middle area of importance loss as seen by CCA1 coincides well with the Niño3.4 and the Niño3 box in the late 21st c. and around the turn of the 20th-21st centuries, respectively (see in Fig. 1), hence the respective insensitivities of these boxes to the increase and the drop in teleconnection strength in these periods (Figs.4(b) and (a)).Note that changes in EOF1 and CCA1 are rather similar away from the middle (see again in Fig. 1), this is why the strengthening in the late 21st c. is much better caught in Fig. 4(c) than in Fig. 4(b).
Having established the observable changes, we can now ask about their drivers in terms of the linear regression model introduced in Sec.2.3.We emphasize that we do not perform statistical tests regarding nonstationarity other than what is to do with Figs. 4, 7, S8 and S9, concerning r only.Therefore, the attribution of a change in r to different factors is not rigorous but rather tentative.
Fig. 3 suggests that ENSO variability, σ Φ , first stagnates then increases with time in the early 21st c. irrespective of the choice of the ENSO representation.This seems to be followed by a decrease, which, however, remains minor relative to the previous increase for Niño3.4,EOF1 and MCA1.We propose this to be the reason why the overall increase in ENSO variability can dominate changes in the correlation coefficient by the late 21st c. for these representations of ENSO, see in Figs.5(b), (d) and (e).According to Fig. 5, however, all other changes and all changes for other ENSO representations are more closely linked to a decrease in a/σ ξ .Whether the general increase in teleconnection strength is dominated by changes in a/σ ξ or σ Φ , the 21st-c.slowdown of this increase is undoubtedly due to the decrease in σ Φ in the same period.This must be so, since Fig. 3 shows that the ratio a/σ ξ is not just decreasing for almost all representations, but there is no evident slowdown in this decrease either.An important nontriviality in the time evolution of a/σ ξ is a jump at the turn of the 20th-21st centuries.This jump is related to a drop in the time evolution of a and serves as a presumable explanation for the drop in the correlation coefficient at the same time.It turns out from Fig. 3 that both a (the inherent coupling strength) and σ ξ (the ENSO-independent part of the IM variability) typically undergo a considerable increase otherwise, with the result of a "winning", but only slightly.The approximate balancing between these two factors may presumably be associated with both the left and right MCA1 modes being very closely resembling the corresponding EOF1 modes. 6mong others, we have thus seen that changes in the strength of ENSO and its coupling with the IM undergo independent changes, which may exert opposite influence on the ENSO-IM teleconnection (as represented by the Pearson correlation coefficient).We note that there should be no apriori universal connection between an increasing ENSO variability and coupling a of the ENSO and IM variability.Such a coupling to do with the precipitation in various locations can certainly have opposite trends (bottom row of Fig. S1).The coupling itself should emerge as an interplay of various physical processes; both thermodynamic and dynamic processes.
Making use of the SLP mostly agrees with the findings obtained with the SST, see Figs. 6-7.The only but important exception is a typically decreasing ENSO variability.We note that the same domain as for SST is obviously not suited to SLP: the corresponding principal modes do not capture the essence of the ENSO phenomenon, as the concentration of the high-amplitude areas to the edges of the domain indicates (see Fig. S2).This should not explain the universal decrease of σ Φ , however, because it is decreasing also when ENSO is represented by the SOI and the box-SOI, which are defined by data outside of the used domain of the modes.Nevertheless, we have checked (not shown) that a narrower box defining the EOF/MCA/CCA modes, extending to 10 • N and 10 • S, still leads to a mostly decreasing ENSO variability, with respective correlation levels largely unchanged.Furthermore, the principal modes are not concentrated on the edges of the domain, but are largely uniform, and slightly concentrated in the middle.
Concerning the next-to-dominant modes of variability, the correlations do not seem to feature much of a forced response (Figs. 3, 6, S8, S9).Using the SST, the ENSO east-west dipole variability is weakening, which is the opposite of the increasing variability projected onto the dominant modes.The weakening dipole variability is more or less countered by the other factors.On the other hand, using the SLP, there seems to be no much of a change even wrt.the drivers, except for a short period at the end of the 21st c.

Nonergodicity
Whenever the forced change is nonlinear, whether due to a nonlinear progression of the forcing or a nonlinear response characteristic, the estimation of the momentary climatology in terms of an ensemble-mean by a temporal mean in a finite window should be biased [11].That is, the ensemble-mean of the temporal mean would not be equal to the ensemble-mean (say, at the middle of the time window), which is termed as nonergodicity.Higherorder statistics should be biased even in the unlikely case that the forced change of the ensemble-mean is linear and the internal variability would not feature a forced change.The linear Pearson correlation coefficient is no exception, i.e., its (traditional) time-based evaluation will be biased even if the ensemble-means of the correlated quantities exhibit a temporally linear forced response.Notwithstanding, a linear time evolution of the Pearson correlation coefficient r(t) itself may lead to a vanishing bias, and this is what we will elaborate on in this section.
Concerning the ENSO-IM teleconnection in the MPI-GE, its various representations, in fact, clearly feature a nonlinear change of r(t), whether it is a monotonic but degressive change (represented by a concave graph) or a nonmonotonic one, as seen in Figs. 3, 6, 4, 7.However, the evaluation of the bias or the degree of nonergodicity is not straightforward in this case since the "signal-to-noise ratio" is rather poor thanks to the relatively small ensemble size.Nevertheless, we can attempt to at least detect nonergodicity by performing a statistical test.The quantity of the so-called test statistics can in turn serve as some quantifier of the degree of nonergodicity.
To the end of constructing a suitable test, we observe that the Fisher transform ẑ = arctanh(r) (e.g.r distinguishes the finite-N sample estimate of r) provides approximately normally distributed iid rv's (wrt.the different years) of standard deviation approximately 1/ √ N − 3 for large enough N and any true value r of the correlation coefficient. 7The same applies, of course, to the finite-τ -window sample correlation coefficient r τ when Φ(t), Ψ(t) are not auto-correlated and their ensemble-means and standard deviations remain constant in time.In such a case, the only reason that the expectation at any time does not match the theoretical value (τ − 3) −1 is the existence of a bias: i.e., nonergodicity.We have four issues to consider.First, z is not available because of the finite N ensemble size.Therefore, we cannot evaluate nonergodicity at any time (every year).Instead, we can concern nonergodicity overall in an interval of length T calculating Second, had ẑτ and ẑ in the above been computed from combinations (without repetition) from the same sample of size T = N , Cochran's theorem [34] would dictate that cv (where c is an appropriate constant factor) would be distributed according to a χ 2 -distribution.This would allow for calculating the p-value the usual way by evaluating the χ 2 -distribution at the level given by the calculated test statistics.However, our setting is somewhat different, which might not be possible to tackle analytically.Nevertheless, one can sample the test statistics and, therefore, determine its quantiles to arbitrary precision.We have done this sampling simply by generating sample correlation realisations by generating realisations of correlated random variables X and Y , where Y = aX + ξ, and X and ξ are normally distributed independent rv's.This is a further assumption for the nature of our variables, from which deviations certainly exist but which are presumably moderate enough to have a secondary effect for the results.Note that the moving-window temporal correlation coefficient rτ,t is calculated upon pregenerating x t and y t , t = 1, . . ., T .We have checked that cv does seem to follow a χ 2 -distribution even in our setting.Third, ENSO is a dynamical phenomenon and Φ, in any representation, is in fact auto-correlated.We have checked that taking X to be governed by an auto-regressive process of order 1, X t+1 = φX t + ξ x , such that the lag-1 autocorrelation is 0.3, shifts the distribution of v little wrt.its std.We use this value of 0.3 in order to calculate the p-value, or, the 0.95 quantile of v.The latter is about 670 using the parameters T = 220, τ = 31, N = 63; while having a serially uncorrelated X t , it is 667.As a fourth issue, there might be a low-frequency influence on the teleconnection, say, via a time-dependent a [35,36,37,38,39].This should increase the width of the distribution of the sample temporal correlation coefficient, i.e., it should be larger than 1/ √ N − 3. We tested the effect of this by introducing an additive perturbation, a t = a 0 + δa t , where δa t is modeled again as an AR(1) process, setting φ = 0.8 and such a noise strength that std[δa t ] = 0.05.With this the 0.95 quantile is 674; that is, the effect is rather small.This corroborates well with the findings of [40,41], namely, that even if there was a low-frequency modulation of the ENSO-IM teleconnection strength, it is too weak to detect from century-long observations.The cancellation effect described by [37] may be an alternative view of this.
Turning to our application, the test statistics is found to be 756 for the ENSO-IM teleconnection representation given by the Niño3-ISMR pair.That is, it is well above the 0.95 quantile, corresponding to a minuscule p-value, so that the hypothesis of ergodicity can be rejected with extremely high confidence.We also evaluate the test statistics corresponding to the global correlation maps seen in Fig. S1.The result of this can be presented as a global map too (Fig. 8), in which a contour for the 0.95 quantile 670 encloses regions where ergodicity can be rejected at the usual significance level.We can see such regions not just in the tropics or in the Equatorial Pacific, but all over the world.Nevertheless, the highest levels of bias/nonergodicity is found indeed at the centre of ENSO variability and elsewhere in the Equatorial Pacific.
We note that even if r(t) changes linearly, or not at all, there could be a trivial source of the bias, namely, a changing ensemble-mean or standard deviation of Φ or Ψ.We redo the calculation of the test statistics in order to eliminate this trivial source.It is done by subtracting the respective ensemble-mean time series from the Φ and Ψ time series, separately for each ensemble member, before calculating the temporal correlation coefficient rτ .The result (Fig. S10) is hardly distinguishable from the original one (Fig. 8) with respect to the patterns, only the values of v are slightly off.That is, in this case the change of the mean state hardly contributes to the bias.In fact, the forced change of the standard deviation can also be a source of bias.We believe that just as the ensemble-mean change, this is also a negligible effect, considering especially that in τ = 30 yr windows the changes of the σ Φ 's in Figs. 3 and 6 seem unlike to be detectable.For this reason, nonergodicity should robustly imply the nonlinearity of r(t).

Discussion and conclusions
We have re-examined the forced response of the ENSO-Indian monsoon (IM) teleconnection as conveyed by the MPI-GE data.One main increment taken by the new analysis is the consideration of spatial aspects of any forced change.This was achieved by determining the dominant as well as the next-to-dominant modes of variability in terms of two leading EOFs and modes of Maximum Covariance Analysis (MCA) and Canonical Correlation Analysis (CCA).In turn, the current teleconnection strength was defined by the ensemble-wise Pearson correlation coefficient between the scalar data series obtained by projecting SST/SLP and precipitation fields onto any of the modes.The plurality of the representation of the ENSO-IM teleconnection allowed us to build a robust picture of forced change.
We have found that the dominant modes of variability in all representations, similarly to the classical ENSO indices, convey a picture of a strengthening teleconnection in terms of a statistical test.However, a slight drop around the turn of the 20th-21st centuries also turns out to be present irrespective of the particular representation, and a late-21st-c.slowdown in the increase of the strength as well.The latter mostly has to do with a curious nonmonotonicity in the change of the ENSO variability in most representations: about midway in the 21st c. rather suddenly it starts to decline.In terms of a linear regression model that can be associated with evaluating the (linear) Pearson correlation coefficient, the typically increasing regression coefficient -which can be viewed as an ENSO-IM coupling strength -also plays a strong role.Although the model's noise strength undergoes a similar increase, which has an opposing role, the change of the coupling is found to dominate.The latter turns out to be the central piece of the robust or consistent picture of the strengthening ENSO-IM teleconnection.We conjecture that a temporary drop in the coupling strength at the turn of the 20th-21st centuries also has an important effect: a slight decrease in the Pearson correlation coefficient at that time.
We leave it for future work to attribute these statistical features to physical effects.However, the change of both the coupling and noise strength should involve both thermodynamic and dynamic factors.We emphasize that all these findings are identified rather robustly across several representations of the ENSO-IM teleconnection.We must also admit, however, that some of these features are not pronounced in all representations and may even not pass a strict detectability threshold in some of them.
These differences between different representations are due to changes in the relevant spatial patterns.The most interesting finding is that the pattern of maximal ENSO variability and that of maximal correlation with the IM precipitation undergo opposing trends in the middle of the ENSO domain: the middle becomes more important for ENSO itself, while it loses importance for the teleconnection.On the side of the IM, however, the intrinsic and teleconnective patterns change very similarly.In particular, most areas over the sea lose much of their importance for both aspects.
It is essential to note that based on observations the ENSO-IM teleconnection appears to have weakened since about 1980 [42,18].In [18] we argued that several reasons could be responsible for this mismatch.Among them one is that the MPI-ESM associated with the MPI-GE is not consistent with the observations of the 20th c.Indian Ocean (IO) warming and IM precipitation decline: in terms of the long-term trend, the ranges of simulated realisations, wrt.both the IO temperature and ISMR, do not contain the observations as respective single realisations -by far (see Figs. 11 and 13 of [18]).See [43], for example, for the role of the IO warming regarding the decline of ISMR, concerning La Niña years, in particular.However, it is a question yet to be answered if the decline of precipitation does actually translate into a weakening of the ENSO-IM teleconnection in the precise sense of a forced change.
Another possibility for the reason for the said mismatch was posited in [18] to be the very considerable fluctuations of temporal correlation coefficients -the larger in magnitude the shorter the time window -even under stationary unforced conditions, and even without decadal variability of the correlated signals.Evaluating temporal correlations is troubled by further problems when there is a forced change, namely, that (i) when the forced change is nonlinear, the temporal correlations are biased (Sec.4), and that (ii) the forced signal is not known, and, therefore, the anomalies that are meant to be correlated cannot even be constructed.A detrending, say, removing a linear slope in a time window, does not seem to be a good fix, as we do not know whether there is any forced change at all, being the very question that we are asking, and even under stationary conditions one would find spurious trends due merely to internal variability.
When an ensemble is available concerning a single model at least, the biases to do with point (i) above can be detected and, to a certain extent, quantified.We have developed an ad hoc statistical test to do just that, and applied it to the ENSO-IM teleconnection as well as the relationship of global precipitation and ENSO.We have found many regions of the world where a bias could be detected.The reasons for the nonlinearities that such biases imply can be manifold, not only driven by the nonlinearity of the change of ENSO variability, and it should be carefully examined in the future.
Our methodology indeed ensures that the detected biases imply nonlinearity (i.e., not only nonlinearity implies biases, but also the other way round).In our context of the historical and scenario change of ENSO-precipitation teleconnections, robust nonlinearity has been found.Given that nonlinearity implies (presupposes) a forced change, we turn out to have on hand a statistical test for the nonstationarity of teleconnections (linear correlations).This can be seen to complement the Mann-Kendall test (in the special case of the nonstationarity of the correlation coefficient), because the latter takes the alternative hypothesis of a monotonic change, while our test includes nonmonotonic changes too.4a-c of [2] concerning the CESM1-LE) and corresponding linear trends (top right; the same is shown in Fig. 5a of [2]) and temporal EOF1 for the full evolution of the correlation maps (middle right).We note that if the response (characteristic) is linear in every gridpoint and the forcing is quasi-statically slow, then the single eigenvalue that is nonzero is the one belonging to EOF1; the forced response of the field is a perfect standing wave.JJA-mean SST and precipitation data is used.

Figure 1 :
Figure 1: Forced change of the first modes of JJA-mean SST variability in the Equatorial Pacific.Temporal means are taken in four consecutive 50-year periods starting from 1900.The top row displays the temporal mean in the first period, and the subsequent rows display the difference with respect to that in the following periods.Left: EOF, middle: MCA, right: CCA first mode.For the MCA and CCA analysis the boxes on the side of the Indian monsoon are shown in Fig. 2. The Niño3 and Niño3.4 boxes are marked in gray and black, respectively.

Figure 2 :
Figure 2: Same as Fig. 1 wrt. the first three columns, but concerning the Indian summer monsoon.In the fourth column a correlation map and its changes are shown concerning the gridpoint-wise JJAS precipitation and the PC1 of the EOF1 of the SST in the box seen in Fig. 1.

Figure 3 :
Figure 3: The forced evolution of correlation coefficients r(t) and the drivers of change; see the main text.The different columns correspond to different representations of the ENSO-IM teleconnections.For columns 1-3, 7 the IM is represented by the average monsoon rain ISMR.Columns 1-6 concern the respective dominant modes of variability, and columns 7-10 concern the next-to-dominant modes of variability.

Figure 5 :
Figure 5: Relative importance of σ Φ versus a/σ ξ to the change of r(t).The diagrams correspond to those of Fig. 4 showing the Mann-Kendall test statistics for r(t).Color is applied only where r(t) changes significantly.The color saturates where αβ is outside of the range indicated in the colorbars.

Figure 6 :
Figure 6: Same as Fig. 3 but the ENSO is characterised based on the SLP.

Figure 7 :
Figure 7: Same as Fig. 4 but the ENSO is characterised based on the SLP.

Figure 8 :
Figure 8: Test statistics of our ad hoc test of nonergodicity regarding the Pearson correlation coefficient.The contour marks the 0.95 quantile.

Figure S 1 :
Figure S 1: Sample snapshot correlation maps at the indicated years (left column; the same is shown in Fig.4a-cof[2] concerning the CESM1-LE) and corresponding linear trends (top right; the same is shown in Fig.5aof[2]) and temporal EOF1 for the full evolution of the correlation maps (middle right).We note that if the response (characteristic) is linear in every gridpoint and the forcing is quasi-statically slow, then the single eigenvalue that is nonzero is the one belonging to EOF1; the forced response of the field is a perfect standing wave.JJA-mean SST and precipitation data is used.

Figure S 2 :
Figure S 2: Same as Fig. 1 (dominant modes) but the ENSO is characterised by the SLP.

Figure S 3 :
Figure S 3: Same as Fig. 1 (SST) but for the next-to-dominant modes.

Figure S 6 :
Figure S 6: Same as Fig. 2 but for the next-to-dominant modes.

Figure S 7 :
Figure S 7: Same as Fig. 6 but the ENSO is characterised by the SLP.

Figure S 10 :
Figure S 10: Same as Fig. 8 but eliminating the contribution from the ensemble-mean changes.