Nonlinear forced change and nonergodicity: The case of ENSO-Indian monsoon and global precipitation teleconnections
Tamás Bódai, Gábor Drótos, Kyung-Ja Ha, June-Yi Lee, Tímea Haszpra, Eui-Seok Chung
NNonlinear forced change and nonergodicity: The case ofENSO-Indian monsoon and global precipitation teleconnections
Tam´as B´odai, G´abor Dr´otos, Kyung-Ja Ha, June-Yi Lee,T´ımea Haszpra and Eui-Seok ChungSeptember 7, 2020
Abstract
We study the forced response of the teleconnection between the El NioSouthern Oscillation (ENSO) andglobal precipitation in general and the Indian summer monsoon (IM) in particular in the Max Planck InstituteGrand Ensemble. The forced response of the teleconnection is defined as the time-dependence of a correlationcoefficient evaluated over the ensemble. The ensemble-wise variability is taken either wrt. spatial averages ordominant spatial modes in the sense of Maximal Covariance Analysis or Canonical Correlation Analysis or EOFanalysis. We find that the strengthening of the ENSO-IM teleconnection is robustly or consistently featured inview of all four teleconnection representations, whether sea surface temperature (SST) or sea level pressure (SLP)is used to characterise ENSO, and both in the historical period and under the RCP8.5 forcing scenario. The maincontributor to this strengthening in terms of a linear regression model is the regression coefficient, which canoutcompete even a declining ENSO variability in view of using the SLP. We also find that the forced changeof 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-wisecorrelation coefficient itself, and also showing that (2) the trivial contributions of the forced changes of e.g. themean SST and/or precipitation to temporal correlations are insignificant here. We also provide, in terms of thetest statistics, global maps of the degree of nonlinearity/nonergodicity of the forced change of the teleconnectionbetween local precipitation and ENSO.
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 onthe CMIP5 archive. They consider, in particular, ENSO-driven precipitation anomalies in tropical regions aroundthe globe, and assess them jointly with changes of mean precipitation. For example, they find that, as a combina-tion, wet DJF anomalies are “strengthened” in the Tibetan plateau both by increasing mean precipitation and wetENSO teleconnection strength; or, as another type of combination, an increasing JJA dry ENSO teleconnectionstrength is counteracted by a wetter mean state in Southeast Asia. Haszpra et al. [2] have evaluated only the trendin ENSO-precipitation teleconnection strength changes, however, not in a model ensemble but the so-called “singlemodel initial condition (large) ensemble” (SMILE) CESM1-LE [3]. Working with a SMILE has the advantagesthat the forced response is correctly represented [4, 5, 6] in that model at least, and seeking a physical interpretationof 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 theprobability measure carried by the climate snapshot attractor [4, 5, 6]. The ensemble-wise statistics evaluated overthe converged ensemble [7] is in a one-to-one correspondence with the external forcing of the nonautonomous dy-namical system, hence the term ‘forced response’. Without convergence, or considering a single realisation, whichmay even correspond to observations, the state depends also on initial conditions, so that forced changes cannot beprecisely 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 forcedchange, before possibly “anomalies”, i.e., (what is assumed to be) the internal variability is analysed. Identifyingthe 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 SMILEsavailable, it is very common to see that authors evaluate temporal statistics first – unnecessarily involving a sub-jective factor – and ensemble statistics afterwards [12, 13]. Otherwise, the biases are certainly controlled by the1 a r X i v : . [ phy s i c s . a o - ph ] S e p agnitude of the climate change signal relative to the intensity of internal variability – a kind of signal-to-noiseratio. Computing the PC1 of EOF1 without detrending, as mentioned just above, is in fact meant to maximise thesignal at least, and the EOF1 is referred to [14] as a “fingerprint” of the forced change – the spatial pattern that issupposed to be distorted by internal variability, in any single realisation, the least. See [15] that is concerned ratherwith the minimisation of the “noise”.It is not pertinent to talk about “advantages” of temporal or ensemble methods over one another, because nosituation arises when we need to decide between them. When an ensemble is available , the correct, conceptuallysound ensemble method is to be applied if the forced response (including that of the internal variability) is con-cerned. 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 ensembledata.Concerning teleconnections, in particular, the forced response can be identified by the time evolution of theensemble-wise Pearson correlation coefficient [16, 17] (in a simple linear approach) between some anomaly quan-tities. The anomaly can correspond to simply a spatial (mean of a temporal) mean [18], but also the PC of anEOF concerning variability across the ensemble , called a snapshot EOF (SEOF) [19]. Haszpra et al. [2] take thelatter approach; however, it can be extended to obtaining anomalies observing the “mutual variability”, such as inthe sense of Maximal Covariance Analysis (MCA) [9] or Canonical Correlation Analysis (CCA) [20, 9]. With aninterest in a teleconnection and its forced response, MCA and CCA – just like EOF analysis – can also be pursuedconcerning the variability across the ensemble, whereby we can refer to these methods as SMCA and SCCA. Thisis perhaps best suited to teleconnection analyses concerning two extended, possibly disconnected, regions. Seean application of SMCA to studying the forced change of the coupling of JJA 200-hPa geopotential height andSeptember sea ice concentration in [21].An important example of this would be the relationship of the summer monsoon precipitation on the Indiansubcontinent with the ENSO phenomenon which extends over the whole of the Equatorial Pacific [22, 23, 24]. Ina previous publication [18] we examined the forced response of the ENSO-Indian monsoon teleconnection in theMax Planck Institute Grand Ensemble (MPI-GE) [25] representing the monsoon by the average JJAS precipitationover India and the ENSO by either the gridpoint- and SLP-based SOI, or the areal mean of the SST in someextended area in the Equatorial Pacific. Standard choices for the latter are the Ni˜no3, Ni˜no4, and Ni˜no3.4 regions.For the first time we established, via a formal hypothesis test, that standard representations of the teleconnection inthe sense of a correlation coefficient between some of the above-mentioned indices strengthened in this model, bothin the historical period and under the strong RCP8.5 forcing scenario, although not necessarily monotonically. Infact, using the SOI it was the historical period only when an increase could be detected, and, using the Ni˜no3 index,it was rather the RCP8.5 scenario under which an increase could be detected. This raises the question whether thediscrepancy is down more to (1) the physical difference between the SOI and Ni˜no3 indices, or (2) rather that – onebased on two distant gridpoints while the other on a limited region of the Equatorial Pacific – they “take somewhatdifferent slices” of the whole ENSO phenomenon, which can matter when even spatial characteristics of ENSOhave a forced response beside possibly the “temporal”. In fact, using the Ni˜no3.4 index, instead of Ni˜no3, theresult 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˜no3.4 andthe box-SOI we can detect changes both under the historical and scenario forcing. On the one hand (I), this seemsto rule out that using SST versus SLP makes much difference. We will see below that this is actually not preciselythe 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 inthe magnitudes of the PCs belonging to either of the said instantaneous/snapshot spatial modes. These are the twobasic questions that we set out here to investigate.Furthermore, this study identifies, in terms of a regression model, controlling factors driving the forced changeof 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-IMteleconnection are not driven only or dominantly by the change of ENSO variability. Nevertheless, it imprints itsnonlinearity onto the teleconnection change, which could be an important source of biases as follows. We observein the MPI-GE a nonmonotonic change in the ENSO variability: after a seemingly monotonic change, a declinefollows in the second half of the 21. c. under RCP8.5. Such a feature is absent in the CESM1-LE [2]. Thisnonmonotonicity and possibly others, as already mentioned above, constituting a nonlinear change, prompt thatthe system should be nonergodic also with respect to correlations, i.e., biases should exist [11] in the temporalcorrelation coefficient evaluated in, say, multi-decadal time windows. Such a bias is something to bear in mindbeside the ample statistical fluctuations of finite-size statistics even under a stationary climate [18]. Given thatthe ensemble size is finite and not so large concerning teleconnections, we develop here a statistical test whereby2e can detect nonergodicity, and, subsequently, map out regions of the world where such a nonergodicity can bedetected 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 analysingthe forced response of teleconnections, such as the way we pursue EOF, MCA, CCA, as well as the idea ofdecomposing the change of the teleconnection strength in terms of a simple regression model. In Sec. 3 we provideresults both on spatial and temporal aspects of the forced change of the ENSO-IM teleconnections. In Sec. 4 weoutline our method of detecting nonergodicity and map out the world with respect to the degree of nonergodicityconcerning the synoptic relationship of the local JJA precipitation with ENSO. In Sec. 5 we discuss and summariseour results. The analysis is restricted to the use of the MPI-GE data [25], that is, to a particular Earth system model (ESM), theMPI-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 SurfaceTemperature (SST). The Indian summer monsoon rain (ISMR) was calculated from the total (convective and largescale) 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/ undervariable name ‘psl’. The SST can be derived from the top layer of the 3D potential temperature field, whosevariable is ‘thetao’. Instead of the total precipitation rate [km/s], one can use the precipitation flux [kg/s/m ],whose variable name is ‘pr’. The scalar quantity of the ISMR averaged in time over the JJAS monsoon season and in space over India is thesame as that calculated in [18], in order to secure a correspondence with the so-called AISMR rain gauge dataset [26] (excluding some states of India; see their Fig. 1). Scalar quantities to represent ENSO variability are alsothe same as in [18], as follows. 1. Average SST in the standard Ni˜no3 region represented by the box (5 ◦ S, 5 ◦ N,210 ◦ E, 270 ◦ E); 2. Average SST in the standard Ni˜no3.4 region (5 ◦ S, 5 ◦ N, 190 ◦ E, 240 ◦ E); 3. SLP difference, p diff = p Tahiti − p Darwin , between the closest gridpoints to Tahiti and Darwin. 4. The difference of the averageSLPs in the boxes (5 ◦ S, 5 ◦ N, 80 ◦ E, 160 ◦ E) and (5 ◦ S, 5 ◦ N, 200 ◦ E, 280 ◦ E). We will refer to the average SSTsunder 1. and 2. as the Ni˜no3 and Ni˜no3.4 indices, respectively, and to the pressure differences under 3. and 4. asthe SOI and box-SOI indices (SOI being short for the Southern Oscillation Index), respectively. See a discussionon this in Sec. 4.a of [18].We correlate the “near-synoptic” JJA mean of any of Ni˜no3, Ni˜no3.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 annualtime 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 responsesignal. That is, the estimate ˆ r ( t ) rather poorly represents the unknown true r ( t ) signal. Therefore, we aspire onlyto detect a change in certain time periods. For this, we employ the Mann-Kendall test [27] in the same way asdone in [18], aiming to reject the hypothesis of stationarity against the alternative hypothesis of a monotonic trendmasked by serially uncorrelated noise constituted by realisations of iid rv’s (independent identically distributedrandom variables).We will evaluate correlations of the AISMR also with the PCs (which are scalars) belonging to EOFs. PCs andEOFs 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] . 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 standingwaves, i.e., time-independent spatial patterns modulated by arbitrary (but uncorrelated) temporal signals, wherebythe spatial patterns are given by the EOFs, and the temporal signals by the PCs. Snapshot EOFs are the samebut different time steps are replaced by different ensemble members (i.e., realizations of the dynamics). EOFs areorthogonal, and each new EOF is defined such that the corresponding PC has the largest possible variance, leadingto an eigen-problem. The variances of the PCs are singular values belonging to the EOFs being the singular Ensemble-wise EOFs, SEOFs, were used earlier in the context of probabilistic ensemble forecasting in e.g. [28, 29, 30], in which case theconcern was certainly not the forced response. svd as done in [31]. The ordering of the singular values wrt. magnitudeprovides 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 havemore in common than the representations given by e.g. the pairs of Ni˜no3-AISMR and SOI-AISMR, becauseNi˜no3 and SOI do not derive from the same mathematical definition of dominance ; and one could think thatmore similar representations yield more similar results. Therefore, firstly, we wish to see if we can establisha robust representation of the ENSO-IM teleconnection by having possibly a better match between the forcedresponse signals belonging to the SLP and SST than between those using the classic indices 1.-4. Secondly, theteleconnection 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 higherEOF modes can reveal more predictive power of the full e.g. SST field compared with e.g. Ni˜no3 alone . Inthis regard we are interested in whether the correlations to do with the higher modes are more susceptible to theanthropogenic 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 similarto the EOF analysis, but consider two fields, and find the orthogonal modes for which the covariance and thecorrelation, respectively, is maximal between the corresponding PCs of the two fields. We carry out these analysesby applying Matlab’s svd and canoncorr , respectively. These methods will allow us to take into account achange 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 theside of the IM. We evaluate the correlation coefficient between the PC1s of the respective 1st SMCA/SCCA modesas well as the PC2s. That is, the correlations of these PC2s are not analogous to those of the AISMR and PC2of EOF2. We anticipate that the SMCA and SCCA lead to higher correlations than the SEOF analysis, stemmingfrom their very definition. In this regard, as the fourth point of inquiry, we want to see if higher correlationsare more susceptible to anthropogenic forcing. Instead of the complete field, the SCCA is performed on the PCscorresponding to the first 4 SEOFs, since it is ill-defined on the full fields. We will speak about SCCA modes asthe 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 SCCAPCs whose correlations were maximised by the SCCA.
Instead of trying to find physical mechanisms responsible for a particular forced change signal in the teleconnectionstrength, 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: r = 1 (cid:113) σ ξ /aσ Φ ) . (1)These factors are: • ENSO variability σ Φ = std [Φ] ; • ENSO-IM “coupling” a being the regression coefficient; • Noise strength σ ξ = std [ ξ ] .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 cansay 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 thatperhaps 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 couldpossibly be seen accidental, to the ENSO-IM teleconnection given that ENSO-precipitation teleconnections can Although, the rationale behind defining simple indices is such that we want to work with a simple-to-compute quantity (e.g. an arealaverage or difference between two locations) that correlates very strongly with the dominant variability in terms of EOFs, which latter can beseen to represent a more natural method of decomposing the variability of a field. Given that the PCs are uncorrelated, one can calculate the so-called coefficient of multiple correlation or determination [9] as the squareroot of the squared correlation coefficients between the “predictand” and the individual PCs, which provides a framework for identifying the contribution of the different PCs to predictability in terms of a multi-dimensional linear regression model. σ Φ 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 isnot constant in time, we can evaluate a by, first, directly evaluating the IM variability σ Ψ , and, subsequently, usingthe textbook formula: a = r σ Ψ σ Φ . (We can, of course, evaluate a directly by linear regression, but calculating the standard deviation is easier, and wealready 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 , t ] the ENSO variability increases as σ Φ ( t ) = βσ Φ ( t ) , β > , and we have also a decrease σ ξ ( t ) /a ( t ) = ασ ξ ( t ) /a ( t ) , α < , then αβ > would mean thatthe ENSO variability increase, in comparison with the σ ξ /a decrease, has a larger effect on the increase of r ( t ) inthat period. Given the very noisy time series, the appropriate approach would be a statistical test attempting to rejectthe 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 , t ] assuming linear changes of σ Φ ( t ) and σ ξ ( t ) /a ( t ) . Thatis, α and β are determined from the respective linear regression model fits. We start by inspecting the forced response of spatial characteristics of the ENSO and IM variability and thoseof their teleconnection. Fig. 1 shows the first modes of the SEOF, SMCA, SCCA analyses on the side of theENSO, based on the SST, smoothed by averaging within four subsequent 50-year periods starting from 1900. Thetemporal averaging brings it really out that hardly any change is visible in the EOF1 and the MCA1 mode, evenupon the strong RCP8.5 forcing. Carreric et al. [13] suggests that the maximum of the EOF1 mode shifts to theeast in the CESM1-LE by showing that some centre (C) mode [32] shifts to the east more in comparison with thewestward shift of some east (E) mode . This does not appear to be so here, although the error of the results (withrespect to any year within a given period) would be difficult to be estimated. Unlike the EOF1 and the MCA1mode (which are very similar to each other), the CCA1 mode does undergo a considerable change by the late 21stcentury. It is even more interesting to observe this change to have a different sign in the middle of the ENSOdomain (highlighted by the Ni˜no3 and Ni˜no3.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 intoaccount the patterns in the reference period) that important areas for ENSO itself and for its teleconnection withthe Indian monsoon get more decoupled. In particular, the middle of the basin gains importance for ENSO butloses 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 thechanges in consecutive rows suggests that we can see meaningful signals. Otherwise, the supplementary video https://youtu.be/RdHjbjZxaZQ shows ample year-to-year change in the form of a flickering. Given thatthe 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), arelargely unchanged, too. We note, however, that the obtained temporal-mean second modes are less reliable thanthe 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 wepick a reference year, and calculate the spatial correlation of the modes in all years with the mode in the referenceyear. When the sign of the correlation is negative, then we flip the sign/orientation of the mode. However, whenthe statistical error corrupting the modes is too large, the pattern matching might not be possible. It seems to be thecase that the pattern matching is feasible for the first modes, but fails in certain years for the second modes. See thesupplementary 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) anomaliesfrom a temporal window for the evaluation of the modes. That is, we mix temporal and ensemble-wise variability. This situation is not shown conclusively in [2] because only three snapshots are shown belonging to single years, and it is not clear howmuch statistical error is associated with them, leaving an uncertainty about the locations of the peaks.
5e emphasize, however, that the anomalies are obtained the same way as done without smoothing: simply bysubtracting the ensemble-mean belonging to single years, so that the only subjective element is the choice of thewindow length. The supplementary videos ( https://youtu.be/bmbPXAtq_p4 https://youtu.be/R78ycuvOI6c https://youtu.be/4pWPf2oUE-I https://youtu.be/RdHjbjZxaZQ ) show thatthe 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 beincurred – whether by averaging 50 annual “snapshot” modes or by pooling anomaly data from the full 50-yearperiod – because there is no reason to think that the finite-ensemble-size estimate is unbiased . It is certainlycomputationally far less intensive by pooling the data, because then the expensive SVD needs to be performedonly once.On the side of the IM, Fig. 2 evidences a change in spatial characteristics, namely, a shift of some “telecon-nection hotspots” in Northern India and over the sea to North-Northwest and to North, respectively, at least in thelast 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 Westof the land, but one may reasonably argue that the teleconnection with ENSO (captured by CCA1 and MCA1) isstrongly bound to the strongest Indian variability (captured by EOF1 and MCA1). In view of this, the coincidenceof 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 andits changes concerning the relationship of the gridpoint-wise precipitation with the PC1 of EOF1 of ENSO, andfind similar patterns to those of the first SEOF/SMCA/SCCA modes, which further supports the close relation-ship between the Indian precipitation variability and the ENSO-influence. Note that many areas of the strongestchanges fall outside the AISMR region. In particular, most of the sea loses very much of its relevance for theIM variability and for its teleconnection with ENSO (a positive change results in getting closer to zero in theseareas 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. UsingSLP instead of SST on the side of the ENSO, the spatial patterns of precipitation over India are largely unchangedfor 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 forthe first ENSO modes, CCA1 gives prominence to the variability in the westernmost part of the selected box incomparison with EOF1 and MCA1. Presumably this is in large part for the intuitive reason that for a maximalcorrelation (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 ofEOF2 and MCA2. The SLP shows a similar behaviour. As for the IM, MCA2 is an outlier, but matching is verygood in all other cases. The time evolution of the correlation coefficients is given in the top row of Fig. 3 for the SST-based representationsof ENSO. This time evolution is basically a long-term increase for the Ni˜no indices (reproduced from [18]) andfor the first modes. The second modes behave differently: their correlation coefficients are practically constant, sothat they really cannot have higher than secondary importance for changes in the teleconnection strength.For the Ni˜no indices and the first modes, the details of the detectable time evolution are presented in Fig. 4by temporal maps of the Mann–Kendall test statistic for stationarity, following [18]. The general increase inteleconnection strength is detected by every quantifier of ENSO and the IM. Although a change is hardly detectedwithin the 21st c. part alone (upper right triangle), the gradual darkening of the shades in the upper left box suggestsa continued strengthening of the teleconnection. Note that the radiative forcing is much stronger in the 21st c. thanin the 20th c., so that any strengthening in the 21st c. must be nonlinearly slowed down in comparison with the20th c.The further details also mostly match in all representations utilizing the AISMR (Figs. 4(a)-(c)), although anearly-21st-c. drop does not surpass the threshold for 5% significance for Ni˜no3, and the late-21st-c. increase israther moderate for Ni˜no3.4. This feature is even more pronounced in the EOF1-EOF1 and the MCA1 represen-tation (Figs. 4(d)-(e)), whereas the CCA1 representation comes with a very prominent gradual increase extendingto all of the 21st c. (Fig. 4(f)). Interestingly, all representations based on principal modes on the IM side place theslight drop at the turn of the 20th-21st centuries considerably earlier than the AISMR-based representations. As an example, the sample correlation coefficient for normally distributed rv.s is negatively biased [33] similar away from the middle (see again in Fig. 1), this is why thestrengthening 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 regressionmodel introduced in Sec. 2.3. We emphasize that we do not perform statistical tests regarding nonstationarity otherthan what is to do with Figs. 4, 7, S8 and S9, concerning r only. Therefore, the attribution of a change in r todifferent factors is not rigorous but rather tentative.Fig. 3 suggests that ENSO variability, σ Φ , first stagnates then increases with time in the early 21st c. irrespectiveof the choice of the ENSO representation. This seems to be followed by a decrease, which, however, remains minorrelative to the previous increase for Ni˜no3.4, EOF1 and MCA1. We propose this to be the reason why the overallincrease in ENSO variability can dominate changes in the correlation coefficient by the late 21st c. for theserepresentations of ENSO, see in Figs. 5(b), (d) and (e). According to Fig. 5, however, all other changes andall changes for other ENSO representations are more closely linked to a decrease in a/σ ξ . Whether the generalincrease in teleconnection strength is dominated by changes in a/σ ξ or σ Φ , the 21st-c. slowdown of this increaseis 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. Thisjump is related to a drop in the time evolution of a and serves as a presumable explanation for the drop in thecorrelation 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 theresult of a “winning”, but only slightly. The approximate balancing between these two factors may presumablybe associated with both the left and right MCA1 modes being very closely resembling the corresponding EOF1modes. Among others, we have thus seen that changes in the strength of ENSO and its coupling with the IM undergoindependent changes, which may exert opposite influence on the ENSO-IM teleconnection (as represented by thePearson correlation coefficient). We note that there should be no apriori universal connection between an increasingENSO variability and coupling a of the ENSO and IM variability. Such a coupling to do with the precipitation invarious locations can certainly have opposite trends (bottom row of Fig. S1). The coupling itself should emerge asan 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 butimportant exception is a typically decreasing ENSO variability. We note that the same domain as for SST is obvi-ously 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 shouldnot explain the universal decrease of σ Φ , however, because it is decreasing also when ENSO is represented bythe SOI and the box-SOI, which are defined by data outside of the used domain of the modes. Nevertheless, wehave 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. Further-more, the principal modes are not concentrated on the edges of the domain, but are largely uniform, and slightlyconcentrated in the middle.Concerning the next-to-dominant modes of variability, the correlations do not seem to feature much of a forcedresponse (Figs. 3, 6, S8, S9). Using the SST, the ENSO east-west dipole variability is weakening, which is theopposite of the increasing variability projected onto the dominant modes. The weakening dipole variability is more Although, the levels of correlation are noticeably different, which should be because the “western hump” of the MCA1 mode is somewhatmore emphasized than that of the EOF1.
7r less countered by the other factors. On the other hand, using the SLP, there seems to be no much of a changeeven wrt. the drivers, except for a short period at the end of the 21st c.
Whenever the forced change is nonlinear, whether due to a nonlinear progression of the forcing or a nonlinearresponse characteristic, the estimation of the momentary climatology in terms of an ensemble-mean by a temporalmean in a finite window should be biased [11]. That is, the ensemble-mean of the temporal mean would not beequal to the ensemble-mean (say, at the middle of the time window), which is termed as nonergodicity. Higher-order statistics should be biased even in the unlikely case that the forced change of the ensemble-mean is linearand the internal variability would not feature a forced change. The linear Pearson correlation coefficient is noexception, i.e., its (traditional) time-based evaluation will be biased even if the ensemble-means of the correlatedquantities exhibit a temporally linear forced response. Notwithstanding, a linear time evolution of the Pearsoncorrelation 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 anonlinear change of r ( t ) , whether it is a monotonic but degressive change (represented by a concave graph) or anonmonotonic one, as seen in Figs. 3, 6, 4, 7. However, the evaluation of the bias or the degree of nonergodicityis not straightforward in this case since the “signal-to-noise ratio” is rather poor thanks to the relatively smallensemble size. Nevertheless, we can attempt to at least detect nonergodicity by performing a statistical test. Thequantity 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 ˆ z = arctanh(ˆ r ) (e.g. ˆ r distinguishes the finite- N sample estimate of r ) provides approximately normally distributed iid rv’s (wrt. thedifferent years) of standard deviation approximately / √ N − for large enough N and any true value r of thecorrelation coefficient. The 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. Insuch a case, the only reason that the expectation (cid:104) (ˆ z τ − z ) (cid:105) (2)at any time does not match the theoretical value ( τ − − is the existence of a bias: (cid:104) ˆ z τ (cid:105) (cid:54) = z, (3)i.e., nonergodicity. We have four issues to consider. First, z is not available because of the finite N ensemblesize. Therefore, we cannot evaluate nonergodicity at any time (every year). Instead, we can concern nonergodicityoverall in an interval of length T calculating v = T (cid:88) t =1 N (cid:88) n =1 (ˆ z τ,t,n − ˆ z t ) . (4)Second, had ˆ z τ and ˆ z in the above been computed from combinations (without repetition) from the same sampleof size T = N , Cochran’s theorem [34] would dictate that cv (where c is an appropriate constant factor) would bedistributed according to a χ -distribution. This would allow for calculating the p-value the usual way by evaluatingthe χ -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 correlationrealisations 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 whichdeviations 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 χ -distribution even in our setting. Third, ENSO isa 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 is0.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, the0.95 quantile of v . The latter is about 670 using the parameters T = 220 , τ = 31 , N = 63 ; while having a serially This approximation assumes the correlated variables to follow a Gaussian distribution, but non-Gaussianity has been checked to have afarly negligible effect for some basic ENSO-quantifiers and AISMR in [18]. 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 sampletemporal correlation coefficient, i.e., it should be larger than / √ N − . We tested the effect of this by introducingan additive perturbation, a t = a + δa t , where δa t is modeled again as an AR(1) process, setting φ = 0 . andsuch a noise strength that std [ δa t ] = 0 . . With this the 0.95 quantile is 674; that is, the effect is rather small. Thiscorroborates well with the findings of [40, 41], namely, that even if there was a low-frequency modulation of theENSO-IM teleconnection strength, it is too weak to detect from century-long observations. The cancellation effectdescribed 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 representationgiven by the Ni˜no3-ISMR pair. That is, it is well above the 0.95 quantile, corresponding to a minuscule p-value, sothat the hypothesis of ergodicity can be rejected with extremely high confidence. We also evaluate the test statisticscorresponding to the global correlation maps seen in Fig. S1. The result of this can be presented as a global maptoo (Fig. 8), in which a contour for the 0.95 quantile 670 encloses regions where ergodicity can be rejected at theusual significance level. We can see such regions not just in the tropics or in the Equatorial Pacific, but all over theworld. Nevertheless, the highest levels of bias/nonergodicity is found indeed at the centre of ENSO variability andelsewhere 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, achanging ensemble-mean or standard deviation of Φ or Ψ . We redo the calculation of the test statistics in order toeliminate 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 τ . Theresult (Fig. S10) is hardly distinguishable from the original one (Fig. 8) with respect to the patterns, only the valuesof v are slightly off. That is, in this case the change of the mean state hardly contributes to the bias. In fact, theforced change of the standard deviation can also be a source of bias. We believe that just as the ensemble-meanchange, this is also a negligible effect, considering especially that in τ = 30 yr windows the changes of the σ Φ ’s inFigs. 3 and 6 seem unlike to be detectable. For this reason, nonergodicity should robustly imply the nonlinearityof r ( t ) . We have re-examined the forced response of the ENSO-Indian monsoon (IM) teleconnection as conveyed by theMPI-GE data. One main increment taken by the new analysis is the consideration of spatial aspects of any forcedchange. This was achieved by determining the dominant as well as the next-to-dominant modes of variabilityin terms of two leading EOFs and modes of Maximum Covariance Analysis (MCA) and Canonical CorrelationAnalysis (CCA). In turn, the current teleconnection strength was defined by the ensemble-wise Pearson correlationcoefficient between the scalar data series obtained by projecting SST/SLP and precipitation fields onto any of themodes. The plurality of the representation of the ENSO-IM teleconnection allowed us to build a robust picture offorced change.We have found that the dominant modes of variability in all representations, similarly to the classical ENSOindices, convey a picture of a strengthening teleconnection in terms of a statistical test. However, a slight droparound 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 curiousnonmonotonicity in the change of the ENSO variability in most representations: about midway in the 21st c. rathersuddenly 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 anENSO-IM coupling strength – also plays a strong role. Although the model’s noise strength undergoes a similarincrease, which has an opposing role, the change of the coupling is found to dominate. The latter turns out to bethe central piece of the robust or consistent picture of the strengthening ENSO-IM teleconnection. We conjecturethat a temporary drop in the coupling strength at the turn of the 20th-21st centuries also has an important effect: aslight 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 ofboth the coupling and noise strength should involve both thermodynamic and dynamic factors. We emphasize thatall these findings are identified rather robustly across several representations of the ENSO-IM teleconnection. Wemust also admit, however, that some of these features are not pronounced in all representations and may even notpass a strict detectability threshold in some of them.These differences between different representations are due to changes in the relevant spatial patterns. Themost interesting finding is that the pattern of maximal ENSO variability and that of maximal correlation with the9M precipitation undergo opposing trends in the middle of the ENSO domain: the middle becomes more importantfor ENSO itself, while it loses importance for the teleconnection. On the side of the IM, however, the intrinsic andteleconnective patterns change very similarly. In particular, most areas over the sea lose much of their importancefor both aspects.It is essential to note that based on observations the ENSO-IM teleconnection appears to have weakened sinceabout 1980 [42, 18]. In [18] we argued that several reasons could be responsible for this mismatch. Among themone is that the MPI-ESM associated with the MPI-GE is not consistent with the observations of the 20th c. IndianOcean (IO) warming and IM precipitation decline: in terms of the long-term trend, the ranges of simulated reali-sations, 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 declineof ISMR, concerning La Ni˜na years, in particular. However, it is a question yet to be answered if the decline ofprecipitation 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 fluc-tuations of temporal correlation coefficients – the larger in magnitude the shorter the time window – even understationary unforced conditions, and even without decadal variability of the correlated signals. Evaluating tempo-ral correlations is troubled by further problems when there is a forced change, namely, that (i) when the forcedchange 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, removinga linear slope in a time window, does not seem to be a good fix, as we do not know whether there is any forcedchange at all, being the very question that we are asking, and even under stationary conditions one would findspurious 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 bedetected and, to a certain extent, quantified. We have developed an ad hoc statistical test to do just that, and appliedit to the ENSO-IM teleconnection as well as the relationship of global precipitation and ENSO. We have foundmany regions of the world where a bias could be detected. The reasons for the nonlinearities that such biases implycan be manifold, not only driven by the nonlinearity of the change of ENSO variability, and it should be carefullyexamined in the future.Our methodology indeed ensures that the detected biases imply nonlinearity (i.e., not only nonlinearity impliesbiases, but also the other way round). In our context of the historical and scenario change of ENSO-precipitationteleconnections, 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). Thiscan be seen to complement the Mann-Kendall test (in the special case of the nonstationarity of the correlationcoefficient), because the latter takes the alternative hypothesis of a monotonic change, while our test includesnonmonotonic changes too.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationshipsthat could be construed as a potential conflict of interest.
Author Contributions
TB developed any new code, including those implementing SMCA and SCCA, carried out all the numerical anal-yses and created all the figures for the paper. GD and subsequently TB originated the ideas of SMCA and SCCA,respectively. TB proposed the breaking down of the change of the correlation coefficients into three constituentsaccording to the linear regression model. TB devised the statistical test for nonergodicity of correlations motivatedby Axel Timmermann’s enquiry. GD, J-Y L, E-S C suggested to calculate the differences as e.g. in Fig. 1; andto evaluate the relative importance of ENSO variability as seen in Fig. 5 developed and implemented by TB. TBand GD lead the writing of the manuscript and the interpretation of results. K-J H, J-Y L, E-S C contributed to themanuscript. 10 unding
TB was supported by the Institute for Basic Science (IBS), under IBS-R028-D1. G.D. acknowledges financialsupport from the Spanish State Research Agency through the Mar´ıa de Maeztu Program for Units of Excellence inR&D under grant no. MDM-2017-0711, from the National Research, Development and Innovation Office (NKFIH,Hungary) under grant no. K125171, and from the European Social Fund through the fellowship no. PD/020/2018under the postdoctoral program of the Government of the Balearic Islands (CAIB, Spain).
Acknowledgments
TB thanks Malte Stuecker for useful discussions and Axel Timmermann for 1. his interest in nonergodicity and2. discussing whether nonergodicity implies the nonlinear change of correlations and vice versa. TB thanks KeithRodgers for informing him about the volume of
Variations that includes [21]. GD is thankful to B. Stevens, T.Mauritsen, Y. Takano, and N. Maher for providing access to the output of the MPI-ESM ensembles.
Supplemental Data
Supplementary Figures are available online, on the webpage of the article. Supplementary Videos are availableonline on YouTube, as cited in the text.
Data Availability Statement
The data generated and analysed for this study are available upon reasonable request. For reproducing purposes,one can make use of the data available at https://esgf-data.dkrz.de/projects/mpi-ge/ , as sug-gested in Sec. 2.
References [1] Scott B. Power and Franois P. D. Delage. El NioSouthern Oscillation and Associated Climatic Conditionsaround the World during the Latter Half of the Twenty-First Century.
Journal of Climate , 31(15):6189–6207,2018.[2] T. Haszpra, M. Herein, and T. B´odai. Investigating ENSO and its teleconnections under climate change in anensemble view – a new perspective.
Earth System Dynamics , 11(1):267–280, 2020.[3] J. E. Kay, C. Deser, A. Phillips, A. Mai, C. Hannay, G. Strand, J. M. Arblaster, S. C. Bates, G. Danabasoglu,J. Edwards, M. Holland, P. Kushner, J.-F. Lamarque, D. Lawrence, K. Lindsay, A. Middleton, E. Munoz,R. Neale, K. Oleson, L. Polvani, and M. Vertenstein. The Community Earth System Model (CESM) LargeEnsemble Project: A Community Resource for Studying Climate Change in the Presence of Internal ClimateVariability.
Bulletin of the American Meteorological Society , 96(8):1333–1349, 2015.[4] Tams Bdai and Tams Tl. Annual variability in a conceptual climate model: Snapshot attractors, hystere-sis in extreme events, and climate sensitivity.
Chaos: An Interdisciplinary Journal of Nonlinear Science ,22(2):023110, 2012.[5] Gbor Drtos, Tams Bdai, and Tams Tl. Probabilistic Concepts in a Changing Climate: A Snapshot AttractorPicture.
Journal of Climate , 28(8):3275–3288, 2015.[6] T. T´el, T. B´odai, G. Dr´otos, T. Haszpra, M. Herein, B. Kasz´as, and M. Vincze. The Theory of Parallel ClimateRealizations.
Journal of Statistical Physics , 2019.[7] G´abor Dr´otos, Tam´as B´odai, and Tam´as T´el. On the importance of the convergence to climate attractors.
TheEuropean Physical Journal Special Topics , 226(9):2031–2038, 2017.[8] Christian L. E. Franzke. Nonlinear climate change.
Nature Climate Change , 4(6):423–424, 2014.119] Hans von Storch and Francis W. Zwiers.
Statistical Analysis in Climate Research . Cambridge UniversityPress, 1999.[10] Byeong-Hee Kim and Kyung-Ja Ha. Observed changes of global and western Pacific precipitation associatedwith global warming SST mode and mega-ENSO SST mode.
Climate Dynamics , 45(11):3067–3075, 2015.[11] G´abor Dr´otos, Tam´as B´odai, and Tam´as T´el. Quantifying nonergodicity in nonautonomous dissipative dy-namical systems: An application to climate change.
Phys. Rev. E , 94:022214, Aug 2016.[12] Benjamin Vega-Westhoff and Ryan L. Sriver. Analysis of ENSO’s response to unforced variability andanthropogenic forcing using CESM.
Scientific Reports , 7(1):18047, 2017.[13] Aude Carr´eric, Boris Dewitte, Wenju Cai, Antonietta Capotondi, Ken Takahashi, Sang-Wook Yeh, GuojianWang, and Virginie Gu´emas. Change in strong Eastern Pacific El Ni˜no events dynamics in the warmingclimate.
Climate Dynamics , 54(1):901–918, 2020.[14] Benjamin D. Santer, John C. Fyfe, Susan Solomon, Jeffrey F. Painter, C´eline Bonfils, Giuliana Pallotta, andMark D. Zelinka. Quantifying stochastic uncertainty in detection time of human-caused climate signals.
Proceedings of the National Academy of Sciences , 116(40):19821–19827, 2019.[15] A. Timmermann. Detecting the Nonstationary Response of ENSO to Greenhouse Warming.
Journal of theAtmospheric Sciences , 56(14):2313–2325, 1999.[16] Mtys Herein, Jnos Mrfy, Gbor Drtos, and Tams Tl. Probabilistic Concepts in Intermediate-Complexity Cli-mate Models: A Snapshot Attractor Picture.
Journal of Climate , 29(1):259–272, 2016.[17] Vineel Yettella, Jeffrey B. Weiss, Jennifer E. Kay, and Angeline G. Pendergrass. An Ensemble CovarianceFramework for Quantifying Forced Climate Variability and Its Time of Emergence.
Journal of Climate ,31(10):4117–4133, 2018.[18] Tam´as B´odai, G´abor Dr´otos, M´aty´as Herein, Frank Lunkeit, and Valerio Lucarini. The Forced Responseof the El Ni˜noSouthern OscillationIndian Monsoon Teleconnection in Ensembles of Earth System Models.
Journal of Climate , 33(6):2163–2182, 2020.[19] Tmea Haszpra, Dniel Topl, and Mtys Herein. On the Time Evolution of the Arctic Oscillation and RelatedWintertime Phenomena under Different Forcing Scenarios in an Ensemble Approach.
Journal of Climate ,33(8):3107–3124, 2020.[20] Wolfgang H¨ardle and L´eopold Simar.
Canonical Correlation Analysis , pages 321–330. Springer BerlinHeidelberg, Berlin, Heidelberg, 2007.[21] T. Haszpra, D. Top´al, and M. Herein. Detecting forced changes in internal variability using large ensembles:On the use of methods based on the “snapshot view”. In C. Deser and K. Rodgers, editors,
New researchon climate variability and change using initial-condition Large Ensembles. Variations. , volume 18, pages36–43. UCAR, 2 2020.[22] Bin Wang.
The Asian Monsoon . Springer-Verlag, 2006.[23] Vimal Mishra, Brian V. Smoliak, Dennis P. Lettenmaier, and John M. Wallace. A prominent pattern of year-to-year variability in indian summer monsoon rainfall.
Proceedings of the National Academy of Sciences ,109(19):7213–7217, 2012.[24] Bin Wang, Baoqiang Xiang, Juan Li, Peter J. Webster, Madhavan N. Rajeevan, Jian Liu, and Kyung-Ja Ha.Rethinking indian monsoon rainfall prediction in the context of recent global warming.
Nature Communica-tions , 6(1):7154, May 2015.[25] Nicola Maher, Sebastian Milinski, Laura Suarez-Gutierrez, Michael Botzet, Mikhail Dobrynin, Luis Ko-rnblueh, Jrgen Krger, Yohei Takano, Rohit Ghosh, Christopher Hedemann, Chao Li, Hongmei Li, ElisaManzini, Dirk Notz, Dian Putrasahan, Lena Boysen, Martin Claussen, Tatiana Ilyina, Dirk Olonscheck,Thomas Raddatz, Bjorn Stevens, and Jochem Marotzke. The Max Planck Institute Grand Ensemble: En-abling the Exploration of Climate System Variability.
Journal of Advances in Modeling Earth Systems ,11(7):2050–2069, 2019. 1226] B. Parthasarathy, A. A. Munot, and D. R. Kothawale. All-India monthly and seasonal rainfall series: 1871–1993.
Theoretical and Applied Climatology , 49(4):217–224, Dec 1994.[27] Henry B. Mann. Nonparametric Tests Against Trend.
Econometrica , 13(3):245–259, 1945.[28] Minghua Zheng, Edmund K. M. Chang, and Brian A. Colle. Ensemble Sensitivity Tools for AssessingExtratropical Cyclone Intensity and Track Predictability.
Weather and Forecasting , 28(5):1133–1156, 2013.[29] Minghua Zheng, Edmund K. M. Chang, Brian A. Colle, Yan Luo, and YueJian Zhu. Applying Fuzzy Clus-tering to a Multimodel Ensemble for U.S. East Coast Winter Storms: Scenario Identification and ForecastVerification.
Weather and Forecasting , 32(3):881–903, 2017.[30] Minghua Zheng, Edmund K. M. Chang, and Brian A. Colle. Evaluating U.S. East Coast Winter Storms in aMultimodel Ensemble Using EOF and Clustering Approaches.
Monthly Weather Review , 147(6):1967–1987,2019.[31] H. Bj¨ornsson and S. A. Venegas. A Manual for EOF and SVD Analyses of Climatic Data. Technical report,n.a.[32] K. Takahashi, A. Montecinos, K. Goubanova, and B. Dewitte. ENSO regimes: Reinterpreting the canonicaland Modoki El Nio.
Geophysical Research Letters , 38(10), 2011.[33] R. A. Fisher. Frequency distribution of the values of the correlation coefficient in samples from an indefinitelylarge population.
Biometrika , 10(4):507–521, 1915.[34] W. G. Cochran. The distribution of quadratic forms in a normal system, with applications to the analysis ofcovariance.
Mathematical Proceedings of the Cambridge Philosophical Society , 30(2):178191, 1934.[35] Alexander Gershunov and Tim P. Barnett. Interdecadal Modulation of ENSO Teleconnections.
Bulletin ofthe American Meteorological Society , 79(12):2715–2726, 12 1998.[36] Christopher Torrence and Peter J. Webster. Interdecadal Changes in the ENSOMonsoon System.
Journal ofClimate , 12(8):2679–2690, 08 1999.[37] V. Krishnamurthy and B. N. Goswami. Indian MonsoonENSO Relationship on Interdecadal Timescale.
Journal of Climate , 13(3):579–595, 02 2000.[38] Lakshmi Krishnamurthy and V. Krishnamurthy. Influence of pdo on south asian summer monsoon andmonsoon–enso relation.
Climate Dynamics , 42(9):2397–2410, May 2014.[39] Takeshi Watanabe and Koji Yamazaki. Decadal-Scale Variation of South Asian Summer Monsoon Onset andIts Relationship with the Pacific Decadal Oscillation.
Journal of Climate , 27(13):5163–5173, 06 2014.[40] Alexander Gershunov, Niklas Schneider, and Tim Barnett. Low-Frequency Modulation of the ENSOIndianMonsoon Rainfall Relationship: Signal or Noise?
Journal of Climate , 14(11):2486–2492, 06 2001.[41] Kyung-Sook Yun and Axel Timmermann. Decadal monsoon-enso relationships reexamined.
GeophysicalResearch Letters , 45(4):2014–2021, 2018.[42] K. Krishna Kumar, Balaji Rajagopalan, and Mark A. Cane. On the Weakening Relationship Between theIndian Monsoon and ENSO.
Science , 284(5423):2156–2159, 1999.[43] S. Aneesh and S. Sijikumar. Changes in the la nia teleconnection to the indian summer monsoon duringrecent period.
Journal of Atmospheric and Solar-Terrestrial Physics , 167:74 – 79, 2018.13igure 1: Forced change of the first modes of JJA-mean SST variability in the Equatorial Pacific. Temporal meansare taken in four consecutive 50-year periods starting from 1900. The top row displays the temporal mean in thefirst 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 Indianmonsoon are shown in Fig. 2. The Ni˜no3 and Ni˜no3.4 boxes are marked in gray and black, respectively.
Supplementary Material r ( t ) and the drivers of change; see the main text. Thedifferent columns correspond to different representations of the ENSO-IM teleconnections. For columns 1-3, 7the IM is represented by the average monsoon rain ISMR. Columns 1-6 concern the respective dominant modes ofvariability, and columns 7-10 concern the next-to-dominant modes of variability.16igure 4: Test statistics of the Mann-Kendall test for the stationarity of r ( t ) or − r ( t ) , depending on which of thesegives an overall positive value. The diagrams correspond to columns 1-6 of Fig. 3 such as: (a) Ni˜no3-AISMR; (b)Ni˜no3.4-AISMR; (c) EOF1-AISMR; (d) EOF1-EOF1; (e) MCA1; (f) CCA1. Red and blue shades correspond to p < . , i.e., a detection of nonstationarity at that significance level.17igure 5: Relative importance of σ Φ versus a/σ ξ to the change of r ( t ) . The diagrams correspond to those ofFig. 4 showing the Mann-Kendall test statistics for r ( t ) . Color is applied only where r ( t ) changes significantly.The color saturates where αβαβ