ITAT 2016 Proceedings, CEUR Workshop Proceedings Vol. 1649, pp. 102–109 http://ceur-ws.org/Vol-1649, Series ISSN 1613-0073, c 2016 N. Jajcay, M. Paluš Statistical modelling in climate science Nikola Jajcay1,2 and Milan Paluš1 1 Dept. of Nonlinear Dynamics and Complex Systems, Institute of Computer Science, Academy of Sciences of the Czech Republic 2 Dept. of Atmospheric Physics, Faculty of Mathematics and Physics, Charles University in Prague Abstract: When it comes to modelling in atmospheric and the model errors [3] and these are intrinsic. The problem climate science, the two main types of models are taken with initial errors is usually tackled by considering an en- into account – dynamical and statistical models. The for- semble of model forecasts (instead of just one realization - mer ones have a physical basis: they utilize discretized dif- integration from single initial state), starting with slightly ferential equations with a set of conditions (boundary con- different initial conditions. The model errors are intrinsi- ditions + present state as an initial condition) and model cally connected with the exponential error growth emerg- the system’s state by integrating the equations forward in ing from the chaotic behaviour related to nonlinearities in time. Models of this type are currently used e.g. as a nu- discretized equations [4]. This limits the predictability of merical weather prediction models. The statistical models such GCMs to 6-10 days maximum (e.g. [5]). are considerably different: they are not based on physical mechanisms underlying the dynamics of the modelled sys- tem, but rather derived from the analysis of past weather patterns. An example of such a statistical model based 1.1 Statistical models on the idea of linear inverse modelling, is examined for modelling the El Niño – Southern Oscillation phenomenon with a focus on modelling cross-scale interactions in the The second kind of models used in climate science are sta- temporal sense. Various noise parameterizations and the tistical models. In their design, they are considerably dif- possibility of using a multi-variable model is discussed ferent than the dynamical models in the sense that they among other characteristics of the statistical model. The are not based on physical mechanism underlying the dy- prospect of using statistical models with low complexity namics of the modelled system, but rather derived from as a surrogate model for statistical testing of null hypothe- the analysis of past weather patterns. Probably the most ses is also discussed. used concept is that of inverse stochastic model [6], where the model is designed, then estimated using past data and, finally, stochastically integrated forward in time to obtain 1 Modelling in climate science the prediction. The disadvantages connected to this type of models consist of the selection of variables that capture the Climate models, which rely on the use of quantitative system we are trying to model. Other possible issue could methods to simulate interactions in the climate system, are be the non-stationarity of the modelled system - since the one of the most important tools to predict and asses future statistical model does not involve the underlying physi- climate projections or to study the climate of the past. In cal mechanisms, just the interaction between subsystems general, two types of models are mainly used: dynamical (ignoring hidden variables), the model estimated on some models and statistical models. The base for a dynamical subset of the past data may not correctly capture all possi- model is a set of discretized differential equations which ble states of the system. In other words, the training period are integrated forward in time from the present state, pos- of the past data used to estimate the statistical model may ing as an initial condition. The most prominent example of not capture the full phase space of the modelled system. the usage of dynamical models is without doubt a general circulation model (GCM hereafter). It employs a mathe- The motivation for building a statistical model for par- matical model of circulation of the planetary atmosphere ticular phenomenon, apart from its forecasting, would be and oceans, therefore it uses the Navier-Stokes equations to scale down the complexity of the problem. When we on a rotating sphere (describing a motion of viscous fluid) find some e.g. nonlinear interactions in the observed data, with thermodynamic terms for energy sources and sinks. and we are interested in uncovering the mechanisms, con- The above described model is used in numerical weather structing a models of different complexity and seeking prediction, to infer the reanalysis datasets of the past cli- such interactions in them would help to expose the mech- mate and for future climate projections in climate model anisms and shed some light on the problem. intercomparison projects CMIP3 [1] and CMIP5 [2]. In the following sections, the inverse stochastic model The uncertainties of the forecast arisen from the GCM for forecasting the El Niño - Southern Oscillation (ENSO models are usually classified into two types: the first one is hereafter) phenomenon is built following [7], with the fo- related to the initial errors (errors in determining the “true” cus on various noise parametrizations and possible use of present state of the climate), while the second one is due to multiple variables. Statistical Modelling in Climate Science 103 Figure 1: ENSO phenomenon, its phases and mechanisms: (left) neutral, (center) positive and (right) negative. Figures taken from [10]. 2 Data-based ENSO model (Fig. 1 right) is acting to strengthen the Walker circula- tion, to move the area of persistent precipitation even more The ENSO phenomenon exhibits strong interannual cli- westward, the differences in surface pressure is now larger mate signal and has a great economic and societal impact. and the thermocline is even more tilted. The ENSO tends It originates from the coupled ocean-atmosphere dynam- to naturally oscillate between these three phases without a ics of the tropical Pacific [8], but has a strong influence on distinct period (there is no distinct peak in ENSO signal’s circulation and air-sea interaction also outside the tropical spectrum) and the reasons why are still largely unknown. belt through teleconnections associated with it [9]. The important aspect of ENSO is that its positive phase The ENSO phenomenon expresses itself as a sea surface - El Niño is generally characterized by a larger magnitude temperature (SST hereafter) anomaly and exists in three than its negative phase - La Niña. This statistical skewness distinct phases - the neutral, positive (El Niño) and nega- is one of the indicators that, at least to some extent, the dy- tive (La Niña). The basic physical mechanisms for each of namics of ENSO involves nonlinear processes [11]. At the the phases are depicted in Fig. 1. The normal state of the same time, the most detailed numerical dynamical mod- equatorial Pacific (Fig. 1 left) is warm SST in the western els seem to severely underestimate this nonlinearity [12], basin, near Australia and cold SST in the eastern basin, hence the quality of the forecast is not satisfactory. near the coast of Peru. Above the warm water in the west, From the reviews of statistical models for ENSO fore- the deep convection takes place, where warm and moist air casting before 2000 [13] it was clear, that majority of mod- is ascending to the border of troposphere, creating an area els were still linear, but lately the nonlinear models are of low atmospheric pressure and area of persistent precip- getting more attention (e.g. [14]). In the following, we de- itation. From the upper part of the troposphere, the air is scribe easy-to-interpret nonlinear model for ENSO fore- moving eastward and then it descends already as cold and casting. dry, creating an area of high atmospheric pressure above the eastern equatorial Pacific. From this basin, the air is 2.1 Inverse models blowing westward on the surface, in agreement with the The concept of inverse stochastic models are used as the trade winds, finishing the circulation loop known as the starting point in developing the ENSO model. Let x(t) be Walker circulation. The easterly surface air flow triggers the state vector of anomalies, so x(t) = X(t) − X, where the oceanic surface current to flow poleward, effectively X(t) is the climate state vector (could be multi- or univari- removing water from the surface, thus the water needs to ate climate observations e.g. temperature, pressure etc. be replaced and this is due to the upwelling, where in the or a PCA time series from eigen-decomposition of some equatorial area, the water is upwelled from roughly 50 me- climate field) and X is its time-mean. The evolution of ters depth to the surface. Since the thermocline (a border anomalies could be expressed as between cold deep ocean and warm surface ocean) is lo- cated below 50 meters in the west, the upwelled water is ẋ = Lx + N(x) (1) warm, but in the eastern Pacific the thermocline level is where L is a linear operator, N represents the nonlinear above the 50 meters, thus the upwelled water is cold, cre- terms and dot denotes time derivative. ating the cold SST in the east and warm SST in the west. The simplest type of inverse models is linear inverse The warm phase of ENSO (Fig. 1 center) creates a warm models (LIM, [6]). By assuming, in eq. (1), that N(x)dx ≈ SST anomaly in the eastern Pacific, acting to weaken the Tx dt + dr(0) , where T is the matrix describing linear feed- Walker circulation, to move the area of persistent precip- backs of unresolved (hidden) processes on x and dr(0) is a itation eastward, to diminish the difference between east- white-noise process, eq. (1) could be written as ern and western Pacific surface pressures and to level off the thermocline. Reversely, the negative phase of ENSO dx = B(0) xdt + dr(0) , B(0) = L + T. (2) 104 N. Jajcay, M. Paluš The matrix B(0) and the covariance matrix of the noise The eqs. (4) and (5) describe a wide variety of processes Q ≡ hr(0) r(0)T i can be directly estimated from the ob- in a fashion that explicitly accounts for the modeled pro- served statistics of x by multiple linear regression [15]. cess x feeding back on the noise statistics. The linear mul- The state vector x, or predictor-variable vector, consists of tilevel model is obtained by assuming Ai ≡ 0 and c(0) ≡ 0 amplitudes of corresponding principal components (PCA in eq. (4). Details of the methodology and further discus- analysis [16] yields spatial patterns - empirical orthogonal sion could be found in [7]. functions and its respective time series - principal com- It is well known, that the extreme ENSO events tend to ponents), while the vector of response variables contains occur in boreal winter. From several ways to include this their tendencies ẋ. phase locking to the annual cycle, the alternative approach used here is to include seasonal dependence in the dynam- ical part of the first level. Namely, we assume the matrix 2.2 Nonlinear multilevel model B(0) and vector c(0) to be periodic, with period T = 12 The assumptions of linear, stable dynamics and of additive months: white-noise used to construct LIMs are only valid to cer- B(0) = B0 + Bs sin(2πt/T ) + Bc cos(2πt/T ), tain degree of approximation. In particular, the stochastic (0) forcing dr(0) typically involves serial correlations, and, in c = c0 + cs sin(2πt/T ) + cc cos(2πt/T ) (6) addition, the matrices B(0) and Q obtained from the data In this case, the whole record is used to estimate four exhibit substantial dependence on the lag, that was used to seasonal-dependent coefficients. The model is trained fit them [17]. The two modifications of the basic inverse in the leading EOF (empirical orthogonal function) model, that address both nonlinearity and serial correla- space [16] of tropical Pacific SST anomalies. The opti- tions are taken into account, as in [18]. mal number of state-vector components and the degree of The first modification is obtained by assuming polyno- nonlinearity has to be assessed by cross-validation. The mial, rather than linear form of N(x) in eq. (1), in par- parameters in this paper were used as in [7]. ticular, a quadratic dependence. The ith component Ni (x) could be written as   3 Results (0) (0) Ni (x) ≈ xT Ai x + ti x + ci dt + dri (3) In this section, the brief results are presented of how the The matrices Ai represent the blocks of a third-order ten- statistical model is able to simulate the ENSO signal. The (0) sors, while the vectors bi = li + ti are the rows of the skill of the model is determined in the sense of basic linear matrix B = L + T (as in eq. (2)). These objects, as well (0) ENSO metrics such as the amplitude of the ENSO signal, as components of the vector c(0) , are estimated by multiple the seasonality (since the seasonality is important aspect polynomial regression [19]. of ENSO dynamics) and finally, the power spectrum of ENSO signal. The model is employed as described in the The second modification, considering the serial correla- previous section, the matrices and vectors are estimated tions in residual forcing, is due to the multilevel structure from the previous data and then the model is integrated to of our model. In particular, consider the ith component of obtain the time series of same length as the training data. the first, main level of the inverse stochastic model Since the model is stochastic (forced by a white noise),   (0) (0) (0) we employed an ensemble of 20 members. Each member dxi = xT Ai x + bi + ci dt + dri , (4) is integrated with slightly different initial conditions and these members are referred to as realizations. where x = {xi } is the state vector and matrices Ai , vec- (0) (0) The basic ENSO metric is its amplitude, which could tors bi and the components ci of the vector c(0) as well be characterized by the standard deviation of SST anoma- (0) as the components ri of the residual forcing vector r(0) lies averaged over Nino3.4 box (bounded by 5◦ S - 5◦ N are determined by the least squares. The additional model and 120◦ W - 170◦ W). In Fig. 2 we can see the ENSO level is added to express the known increments dr(0) as a amplitude as derived from the Nino3.4 index [20] (thick linear function of an extended state vector [x, r(0) ]. We black line), along with 20 realizations from the data-based estimate this level’s residual forcing again by the least ENSO model, both linear and quadratic (gold for linear, squares. More levels are added the same way, until the red for quadratic). Lth level’s residual, r(L+1) , becomes white in time, and its As can be seen, the linear model slightly overestimates lag-0 correlation matrix converges to constant, hence the ENSO amplitude, while the quadratic model slightly underestimate the ENSO amplitude. From the spread of (0) (1) (1) dri = bi [x, r(0) ]dt + ri dt, the ensemble members we could infer that the model is (1) (2) (2) sensitive to initial conditions and the forcing. Still, the dri = bi [x, r(0) , r(1) ]dt + ri dt, ensemble averages for both models are within reasonable ... distance from the data borderline, therefore in this aspect (L) (L+1) (L+1) dri = bi [x, r(0) , . . . , r(L) ]dt + ri dt (5) the model performs adequately. Statistical Modelling in Climate Science 105 0.95 ENSO amplitude tive models are more flat in this area of frequencies. In 0.90 the higher frequencies (around annual frequency and less) the power spectra are in agreement. In general, the spectra 0.85 of modelled time series could be said to copy the actual STD of SSTA -- Nino3.4 0.80 Nino3.4 time series. The power spectra were computed 0.75 using the Welch method [21]. 0.70 10 8 0.65 0.60 10 7 linearan linear linear model 1 linear model 2 linear model 3 linear model 4 linear model 5 linear model 6 linear model 7 linearar model 98 linear model 10 linear model 11 linear model 12 linear model 13 linear model 14 linear model 15 linear model 16 linear model 17 linear model 18 linear model 19 l 20 quadratic mod . quadratic model 1 quadratic model 2 quadratic model 3 quadratic model 4 quadratic model 5 quadratic model 6 quadratic model 7 tic mo el 8 quadratic moddel 9 quadratic model 10 quadratic model 11 quadratic model 12 quadratic model 13 quadratic model 14 quadratic model 15 quadratic model 16 quadratic model 17 quadratic model 18 tic mo el 19 del 20 quadra ean quad line model mode me m quadra 10 6 power Figure 2: ENSO amplitude as standard deviation of SST 10 5 anomalies in data (black line) and in 20 realizations of lin- ear (gold) and quadratic (red) model. 10 4 Other metric connected with ENSO amplitude is its sea- sonality. As written above, the ENSO phenomenon ex- 10 3 10 -1 10 0 hibits seasonal changes in variance, with elevated variance period [1/yr] in winter months and lower variance in spring and summer months. This can be also seen in Fig. 3, where the monthly Figure 4: ENSO power spectra estimated using the Welch variance is plotted for the data and for both models. Both method in data (black curve) and in 20 realizations of lin- models are capable of modelling higher variance in winter ear (gold) and quadratic (red) model. Thicker lines repre- months and drop in variance through spring and summer, sent the mean over 20 realizations in the respective model. although the difference in variance is higher in data than in both models. Still, the ensemble averages are reasonably close to the data. 4 Noise parametrization in the model 1.2 The statistical model, once estimated, is integrated for- 1.1 ward in time and forced by a noise - usually a realization 1.0 of spatially correlated random process. In the most in- tuitive and basic case, the last level residuals’ covariance 0.9 matrix is estimated and decomposed using Cholesky fac- STD of SSTA 0.8 torization yielding a lower triangular matrix R. When the model is integrated, the random realization of white noise 0.7 is multiplied by the matrix R, yielding spatially corre- lated white noise which is used as a random forcing in the 0.6 model. The results for quadratic and linear ENSO models 0.5 from the previous section were obtained using this simple noise parametrization, and the question is whether looking 0.4 deeper into the residuals’ structure could aid the model’s Jan b r r y Jun Jul g p t v c performance. Oc Ma Ap De Ma No Fe Au Se month Figure 3: ENSO seasonality as standard deviation per 4.1 Dependence on the system’s state month in data (black curve) and in 20 realizations of linear First refinement for the noise parametrization arises from (gold) and quadratic (red) model. Thicker lines represent the concept of modelling climate processes which exhibit the mean over 20 realizations in the respective model. low-frequency variability (LFV). In this method, we find and select noise samples, snippets, from the past noise The last metric taken into account was the power spec- (residuals) which have forced the system during short time trum of Nino3.4 time series. The spectrum for the Nino3.4 intervals that resemble the LFV phase just preceding the data and both linear and quadratic model realizations can currently observed state, and then use these snippets (or in- be seen in Fig. 4. The main peak in data occurs at roughly formation contained in them) to drive the current state into 5 year period, but still the ensemble averages for respec- the future. For full methodology and discussion, see [22]. 106 N. Jajcay, M. Paluš The found past noise snippets can be used in two differ- 0.90 ENSO amplitude ent ways. The first one (as used in [22]) seeks various snip- pets from the past observations and then directly uses them 0.85 to force the model as an ensemble. When e.g. we find 4 STD of SSTA -- Nino3.4 intervals which resemble the LFV phase, we integrate the model 4 times using all 4 noise snippets directly and than 0.80 average over them. The second version (as used in our study) is to find, say, 100 samples of the past noise clos- 0.75 est to the current state of the system, cluster them together and create covariance matrix from them. Afterwards, the linlienaerar linear 1 linear 2 linear 3 linear 4 linear 5 linear 6 linear 7 line 8 linearar 9 linear 10 linear 11 linear 12 linear 13 linear 14 linear 15 linear 16 linear 17 linear 18 linear 19 20 conditan cond. conditioned 1 conditioned 2 conditioned 3 conditioned 4 conditioned 5 conditioned 6 conditioned 7 con ioned 8 conditditioned 9 conditioned 10 conditioned 11 conditioned 12 conditioned 13 conditioned 14 conditioned 15 conditioned 16 conditioned 17 conditioned 18 conditioned 19 20 Cholesky decomposition is used to obtain the matrix R and ioned mean finally, the random white noise realization is multiplied by me the matrix R. Using this matrix, the spatial covariance of 1.2 the forcing is dependent on the current state of the sys- tem. In both noise parametrizations, the current system 1.1 state could be estimated in multiple ways: either using correlation of the SSA time series, or using the Euclidean 1.0 distance in the subspace spanned by first few EOFs. As can be seen in Fig. 5, although the amplitude statis- STD of SSTA 0.9 tics are not substantially shifted, the transient from high- variance winter period to low-variance spring and summer 0.8 are better captured by the later model, with noise forcing 0.7 conditioned on system’s state. The power spectra for both models are practically the same (not shown). 0.6 4.2 Seasonal dependence of the forcing 0.5 Jan b r r y Jun Jul g p t v c Oc Ma Ap De Ma No Fe Au Se Although the seasonal dependence of the model is cap- month tured in model’s dynamics by fitting the seasonally depen- dent matrices B(0) and c(0) (recall eq. (6)), our analysis Figure 5: ENSO amplitude (upper) and seasonality (bot- showed, that the last level’s residuals still exhibit season- tom) in data (black curve) and in 20 realizations of linear ally dependent amplitude. To address this issue, we com- (gold) and linear with conditioned noise on the system’s puted the standard deviations for each month from the last state (red) model. Thicker lines represent the mean over level’s residuals, then fitted the 5 harmonics of the annual 20 realizations in the respective model. cycle to capture the seasonal dependence, removed this dependence from the residuals, then estimated covariance matrix and subsequently the matrix R and finally gener- The two latter modifications bring just a slight improve- ated spatially correlated white noise realization which was ments into ENSO metrics (not shown), but could have multiplied back by the requisite seasonal amplitude to ac- more substantial advancements in modelling different at- count for the seasonally dependent amplitude of the forc- mospheric phenomena. ing. The fitted harmonics of the annual cycle were selected as 5 Synchronization and causality in the Pi = cos(2πit/T ) + sin(2πit/T ), i = 1, . . . , 5 (7) observed and modelled data and then regressed on the seasonally varying standard de- viation of the last level’s residuals. Better understanding of the complex dynamics of the at- mosphere and climate is one of the challenges for con- temporary science. Considering the climate system as 4.3 Using extended covariance matrix a complex network of interacting subsystems [23] is a The last modification to the noise is to use the extended co- new paradigm bringing new data analysis methods help- variance matrix instead of lag-0 covariance matrix. When ing to detect, describe and predict atmospheric phenom- evaluating system’s state we do not take just the state clos- ena [24]. A crucial step in constructing climate networks est to the current state of the model, but, say 5 consecu- is inference of network links between climate subsys- tive months and construct the extended matrix out of this tem [25]. Directed links determine which subsystems in- snippet. Then the matrix is decomposed using Cholesky fluence other subsystems, i.e. uncover the drivers of at- factorization and used as a spatial correlation matrix R is mospheric phenomena. Inference of causal relationships random forcing generation. from climate data is an intensively developing research Statistical Modelling in Climate Science 107 field, e.g. [26, 27]. Typically, a causal relation is sought between different variables or modes of atmospheric vari- ability. Paluš [28] has open another view at the complexity of atmospheric dynamics by uncovering causal relations or information flow between dynamics on different time scales in the same variable. Recently, phase-phase and also phase-amplitude interactions between dynamics on different temporal scales were observed in the ENSO dy- namics (captured by the Nino3.4 index) using the ap- proach as in [28]. Shortly, we use the continuous wavelet transform to the time series for particular time scales to Figure 6: Phase synchronization (left) and phase- obtain the instantaneous phase and amplitude of the oscil- amplitude causality (right) in Nino3.4 time series. Shown latory mode as is the significance (over 95th percentile against 500 Fourier ψ(t) = s(t) + iŝ(t) = A(t)e(iφ (t)) , (8) transform surrogates) of k-nearest neighbours estimate of mutual information and conditional mutual information. ŝ(t) φ (t) = arctan , (9) s(t) q to uncover the mechanisms of these interactions and shed A(t) = s2 (t) + ŝ2 (t). (10) more light onto the dynamics of ENSO in general. We constructed the ENSO model and repeated the above anal- Then the time series of phase and / or amplitude are used ysis to modeled ensemble of the Nino3.4 time series. to study the interactions. We adopt measures from infor- mation theory, namely mutual information and conditional mutual information, where the mutual information could be expressed as p(x, y) I(X;Y ) = ∑ ∑ p(x, y) log , (11) x∈X y∈Y p(x)p(y) where p(·) is the probability distribution or joint probabil- ity distribution and X and Y are our time series of either phase or amplitude derived from the ENSO SST data. Fi- nally, the measures we are interested in could be written as: Figure 7: Phase synchronization (left) and phase- amplitude causality (right) in modeled Nino3.4 time series • phase synchronization – I(φ1 (t); φ2 (t)), by the data-based model. Shown is the aggregate of 5 real- • phase-phase causality – I(φ1 (t); φ2 (t + τ) − izations of k-nearest neighbours estimate of mutual infor- φ2 (t)|φ2 (t)), mation and conditional mutual information. Significance against 500 Fourier transform surrogate data. • phase-amplitude causality – I(φ1 (t); A2 (t + τ)|A2 (t), A2 (t − η), A2 (t − 2η)), As seen from the analysis of modelled data (Fig. 7), the main phase synchronization bands (annual cycle with quasi-biennal cycle and combination frequencies) are also 5.1 Interactions in the data captured by the modelled data, while the phase - ampli- As can be seen from Fig. 6, in ENSO dynamics captured tude interactions are not very well captured. This might by the Nino3.4 index, the synchronization of annual cycle arise from the low complexity of the model, or the absence with quasi-biennal and combination frequencies (frequen- of some nonlinear interactions in the model design (apart cies that arise from the interactions between annual and from quadratic). the most prominent ENSO period) is observed. Also, the 4-6 year cycle of phase in ENSO dynamics influence the quasi-biennal range of the amplitude time series. 6 Modelling surrogate data with statistical model 5.2 Interactions in the model Surrogate data (or analogous data) is a method to generate Our goal was to simulate the nonlinear cross-scale inter- synthetic data set (time series) that preserve some of the actions in the model. This is important since it might help statistical properties, while omitting the others. One way 108 N. Jajcay, M. Paluš of using them, is to test statistical significance by contra- diction. This involves posing a null hypothesis describing some kind of a process and then generating an ensemble of surrogate data according to null hypothesis using Monte Carlo methods. One of the most used technique for gener- ating surrogate data is the Fourier transform surrogate [29] (FT surrogates), which preserve the linear correlations in the data (periodogram or spectrogram, including autocor- relation) of the time series, but omits any other interactions in them. As an example, consider two intertwined Lorenz sys- tems, where one of them drives the other. Now, using the Figure 8: Phase synchronization (left) and phase- time series in one dimension, say the x dimension from amplitude causality (right) in modelled Nino3.4 time se- both Lorenz systems, we can use some method for detect- ries by the data-based model. Shown is the aggregate of ing causality, e.g. conditional mutual information between 5 realizations of k-nearest neighbours estimate of mutual the two time series of two Lorenz systems. We get the information and conditional mutual information. Signifi- value of conditional mutual information, but this is still cance against 500 surrogate time series created with data- not enough to interpret it in the means of whether there based model. is a causal relationship between them or the result arose by chance. For this purpose, we construct an ensemble of Fourier transform surrogate data (which qualitatively pre- minimum. This way, we can get better idea of the statisti- serves properties of the time series, but allows no causal cal significance of the interactions between subsystems, in relationship between them) and repeat the analysis using particular the nonlinear ones, since we are testing against the very same method on this ensemble and finally com- the model with just linear interactions. pare the value for actual data with the histogram of values obtained from the ensemble of surrogate data. When the 7 Conclusions value from the data exceeds some percentile (e.g. 95th ) of the surrogate data distribution, we say that the causal re- Statistical modelling in climate science is continuously lationship is significant in comparison with e.g. 500 FT getting more attention, since their usage is not limited to surrogates. forecast some of the phenomena of interest (like ENSO), When studying nonlinear cross-scale interactions in but could also be used to infer some of the statistical prop- time series using the above method, the statistical test in- erties and relationships among different subsystems. Since volves creating an ensemble of surrogate, synthetic time the statistical models live in phase space of particularly re- series and repeat the analysis for the whole ensemble. duced dimensionality, when we could observe the inter- Then we computed the percentile, where the observed in- actions of interest, the identification of their sources will teractions could not arose by random chance. Of course, become more feasible. one could use Fourier transform method to generate the We showed that the statistical model with the right set- surrogate time series, effectively posing a null hypothe- tings, which were selected based on careful inspection of sis of a linear process which has the same spectrum to the modelled system, could generate synthetic time series that of an observed data. On the other hand, one can cre- of interest, copying the desired properties of the system - ate a more sophisticated null hypothesis by exploiting the both linear and nonlinear statistics. Since the stochastic- options of a data-based model: when one consider just a ity is the important aspect of the data-based model, var- linear model, omit the dynamical seasonal dependence in ious parametrization techniques exist to correctly model B(0) and c(0) terms (as in eq. (6)) and use the simplest noise the system’s external forcing. Finally, the possibility of parametrisation (just consider the spatial covariance struc- usage of the low complexity model as surrogate data was ture), the model will omit the nonlinear interactions and discussed, showing advantages of usage of such technique could pose as a surrogate data model copying the basic to infer statistical significance. statistical properties of a modelled time series. This way, The outlook for future work combines various differ- the analysis would show whether the cross-scale interac- ent paths which appeared. One direction would be fo- tions are arising from the seasonal dependent dynamics, or cusing on statistical modelling itself, experimenting with from nonlinear (e.g. quadratic) interactions between sub- various variable model, with input time series and their systems and so on. preprocessing and so on and so forth. Other direction When comparing Fig. 6 (testing against 500 Fourier would be connecting the statistical models with dynami- transform surrogates) and Fig. 8 (testing against 500 data- cal ones, in the sense, that statistical models could be used based model surrogates), the significant interactions are for parametrization of e.g. sub-grid phenomena (micro- virtually the same, expect in the latter, the “fluctuations” physics of clouds, local convection etc.) in large coupled (or they might be false positives as well) are attenuated to atmospheric-oceanic models. Statistical Modelling in Climate Science 109 References [19] McCullagh, P., and J. A. Nelder: Generalized Linear Mod- els. Chapman and Hall (1989) [20] Rayner N. A., D. E. Parker, E. B. Horton, C. K. Folland, [1] Meehl, G. A., C. Covey, T. Delworth, M. Latif, B. McA- L. V. Alexander, D. P. Rowell, E. C. Kent and A. Kaplan: vaney, J. F. B. Mitchell, R. J. Stouffer, and K. E. Taylor: Global analyses of sea surface temperature, sea ice, and The WCRP CMIP3 multi-model dataset: A new era in cli- night marine air temperature since the late nineteenth cen- mate change research. A Bull. Amer. Meteor. Soc. 88 (2007) tury. J. Geophys. Res. 108 (2003) 4407 1383–1394 [21] Welch, P. D.: The use of Fast Fourier Transform for the [2] Taylor, K.E., R.J. Stouffer and G.A. Meehl: An Overview estimation of power spectra: A method based on time aver- of CMIP5 and the experiment design. A Bull. Amer. Meteor. aging over short, modified periodograms. IEEE Trans. Audio Soc. 93 (2012) 485–498 AU-15 (1967) 70–73 [3] Bjerknes, V.: Dynamic meteorology and hydrology, Part [22] Chekroun, M. D., D. Kondrashov and M. Ghil: Predicting II. Kinematics. Gibson Bros., Carnegie Institute, New York. stochastic systems by noise sampling, and application to the (1911) El Niño-Southern Oscillation. P. Natl. Acad. Sci. USA 108 [4] Lorenz, E. N.: Deterministic nonperiodic flow. J. Atmos. Sci. (2011) 11766–11771 20 (1963) 130–141 [23] A. A. Tsonis and P. J. Roebber: The architecture of the [5] Van den Dool, H. M.: Long-range weather forecasts through climate network. Physica A: Statistical Mechanics and its numerical and empirical methods. Dyn. Atmos. Oceans 20 Applications, 333 (2004) 497–504 (1994) 247–270 [24] S. Havlin, D. Y. Kenett, E. Ben-Jacob, A. Bunde, R. Cohen, [6] Penland, C.: Random forcing and forecasting using princi- H. Hermann, J. Kantelhardt, J. Kertész, S. Kirkpatrick, J. pal oscillation pattern analysis. Mon. Weat. Rev. 117 (1989) Kurths, et al.: Challenges in network science: Applications 2165–2185 to infrastructures, climate, social systems and economics. [7] Kravtsov, S., D. Kondrashov, and M. Ghil: Multilevel regres- The European Physical Journal Special Topics, 214 (2012) sion modeling of nonlinear processes: Derivation and appli- 273–293 cations to climate variability. J. Climate 18 (2005) 4404— [25] M. Paluš, D. Hartman, J. Hlinka, and M. Vejmelka: Dis- 4424 cerning connectivity from dynamics in climate networks. [8] Philander, S. G. H.: El Niño, La Niña, and the Southern Nonlinear Processes in Geophysics 18 (2011) 751–763 Oscillation. Academic Press (1990) [26] Ebert-Uphoff and Y. Deng: Causal discovery for climate re- [9] Alexander, M. A., I. Bladé, M. Newman, J. R. Lanzante, N.- search using graphical models. Journal of Climate 25 (2012) C. Lau, and J. D. Scott: The atmospheric bridge: The in- 5648–5665 fluence of ENSO teleconnections on air–sea interaction over [27] Y. Deng and I. Ebert-Uphoff: Weakening of atmospheric the global oceans. J. Climate 15 (2002) 2205–2231 information flow in a warming climate in the community cli- [10] WIKIPEDIA . ORG: El Niño–Southern Oscillation mate system model. Geophysical Research Letters 41 (2014) https://en.wikipedia.org/wiki/El_Ni%C3%B1o% 193–200 E2%80%93Southern_Oscillation, downloaded June 27, [28] Paluš, M.: Multiscale atmospheric dynamics: Cross- 2016. frequency phase-amplitude coupling in the air temperature. [11] Ghil, M. and A. W. Robertson: Solving problems with Phys. Rev. Lett. 112 (2014) 1–5 GCMs: General circulation models and their role in the cli- [29] Theiler, J., S. Eubank, A. Longtin, B. Galdrikian, and J. mate modeling hierarchy. Academic Press, (2000) 285–325 Doyne Farmer: Testing for nonlinearity in time series: The [12] Hannachi, A., D. B. Stephenson, and K. R. Sperber: method of surrogate data. Physica D 58 (1992) 77–94 Probability-based methods for quanfifying nonlinearity in the ENSO. Clim. Dyn. 20 (2003) 241–256 [13] Ghil, M. and N. Jiang: Recent forecast skill for the El Niño / Southern Oscillation. Geophys. Res. Lett. 25 (1998) 171- –174 [14] Timmermann, A., H. U. Voss, and R. Pasmanter: Empir- ical dynamical system modeling of ENSO using nonlinear inverse techniques. J. Phys. Oceanogr. 31 (2001) 1579–1598 [15] Wetherill, G. B.: Regression Analysis with Applications. Chapman and Hall (1986) [16] Hannachi, A., I. T. Jolliffe and D. B. Stephenson: Empirical orthogonal functions and related techniques in atmospheric science: A review. Int. J. Climatol. 27 (2007) 1119–1152 [17] Penland, C. and M. Ghil: Forecasting Northern Hemi- sphere 700-mb geopotential height anomalies using empir- ical normal modes. Mon. Wea. Rev. 121 (1993) 2355—2372 [18] Kondrashov, D., S. Kravtsov, A. W. Robertson and M. Ghil: A Hierarchy of Data-Based ENSO Models. J. Climate 18 (2005) 4425–4444