When it comes to modelling in atmospheric and climate science, the two main types of models are taken into account - dynamical and statistical models. The former ones have a physical basis: they utilize discretized differential equations with a set of conditions (boundary conditions + present state as an initial condition) and model the system's state by integrating the equations forward in time. Models of this type are currently used e.g. as a numerical weather prediction models. The statistical models are considerably different: they are not based on physical mechanisms underlying the dynamics of the modelled system, but rather derived from the analysis of past weather patterns. An example of such a statistical model based 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 temporal sense. Various noise parameterizations and the possibility of using a multi-variable model is discussed among other characteristics of the statistical model. The prospect of using statistical models with low complexity as a surrogate model for statistical testing of null hypotheses is also discussed.
Climate models, which rely on the use of quantitative
methods to simulate interactions in the climate system, are
one of the most important tools to predict and asses future
climate projections or to study the climate of the past. In
general, two types of models are mainly used: dynamical
models and statistical models. The base for a dynamical
model is a set of discretized differential equations which
are integrated forward in time from the present state,
posing as an initial condition. The most prominent example of
the usage of dynamical models is without doubt a general
circulation model (GCM hereafter). It employs a
mathematical model of circulation of the planetary atmosphere
and oceans, therefore it uses the Navier-Stokes equations
on a rotating sphere (describing a motion of viscous fluid)
with thermodynamic terms for energy sources and sinks.
The above described model is used in numerical weather
prediction, to infer the reanalysis datasets of the past
climate and for future climate projections in climate model
intercomparison projects CMIP3 [
The uncertainties of the forecast arisen from the GCM
models are usually classified into two types: the first one is
related to the initial errors (errors in determining the “true”
present state of the climate), while the second one is due to
the model errors [
The second kind of models used in climate science are
statistical models. In their design, they are considerably
different than the dynamical models in the sense that they
are not based on physical mechanism underlying the
dynamics of the modelled system, but rather derived from
the analysis of past weather patterns. Probably the most
used concept is that of inverse stochastic model [
The motivation for building a statistical model for particular phenomenon, apart from its forecasting, would be to scale down the complexity of the problem. When we find some e.g. nonlinear interactions in the observed data, and we are interested in uncovering the mechanisms, constructing a models of different complexity and seeking such interactions in them would help to expose the mechanisms and shed some light on the problem.
In the following sections, the inverse stochastic model
for forecasting the El Niño - Southern Oscillation (ENSO
hereafter) phenomenon is built following [
The ENSO phenomenon exhibits strong interannual
climate signal and has a great economic and societal impact.
It originates from the coupled ocean-atmosphere
dynamics of the tropical Pacific [
The ENSO phenomenon expresses itself as a sea surface temperature (SST hereafter) anomaly and exists in three distinct phases - the neutral, positive (El Niño) and negative (La Niña). The basic physical mechanisms for each of the phases are depicted in Fig. 1. The normal state of the equatorial Pacific (Fig. 1 left) is warm SST in the western basin, near Australia and cold SST in the eastern basin, near the coast of Peru. Above the warm water in the west, the deep convection takes place, where warm and moist air is ascending to the border of troposphere, creating an area of low atmospheric pressure and area of persistent precipitation. From the upper part of the troposphere, the air is moving eastward and then it descends already as cold and dry, creating an area of high atmospheric pressure above the eastern equatorial Pacific. From this basin, the air is blowing westward on the surface, in agreement with the trade winds, finishing the circulation loop known as the Walker circulation. The easterly surface air flow triggers the oceanic surface current to flow poleward, effectively removing water from the surface, thus the water needs to be replaced and this is due to the upwelling, where in the equatorial area, the water is upwelled from roughly 50 meters depth to the surface. Since the thermocline (a border between cold deep ocean and warm surface ocean) is located below 50 meters in the west, the upwelled water is warm, but in the eastern Pacific the thermocline level is above the 50 meters, thus the upwelled water is cold, creating the cold SST in the east and warm SST in the west.
The warm phase of ENSO (Fig. 1 center) creates a warm SST anomaly in the eastern Pacific, acting to weaken the Walker circulation, to move the area of persistent precipitation eastward, to diminish the difference between eastern and western Pacific surface pressures and to level off the thermocline. Reversely, the negative phase of ENSO (Fig. 1 right) is acting to strengthen the Walker circulation, to move the area of persistent precipitation even more westward, the differences in surface pressure is now larger and the thermocline is even more tilted. The ENSO tends to naturally oscillate between these three phases without a distinct period (there is no distinct peak in ENSO signal’s spectrum) and the reasons why are still largely unknown.
The important aspect of ENSO is that its positive phase
- El Niño is generally characterized by a larger magnitude
than its negative phase - La Niña. This statistical skewness
is one of the indicators that, at least to some extent, the
dynamics of ENSO involves nonlinear processes [
From the reviews of statistical models for ENSO
forecasting before 2000 [
The concept of inverse stochastic models are used as the starting point in developing the ENSO model. Let x(t) be the state vector of anomalies, so x(t) = X(t) − X, where X(t) is the climate state vector (could be multi- or univariate climate observations e.g. temperature, pressure etc. or a PCA time series from eigen-decomposition of some climate field) and X is its time-mean. The evolution of anomalies could be expressed as x˙ = Lx + N(x) (1) where L is a linear operator, N represents the nonlinear terms and dot denotes time derivative.
The simplest type of inverse models is linear inverse
models (LIM, [
The assumptions of linear, stable dynamics and of additive
white-noise used to construct LIMs are only valid to
certain degree of approximation. In particular, the stochastic
forcing dr(0) typically involves serial correlations, and, in
addition, the matrices B(0) and Q obtained from the data
exhibit substantial dependence on the lag, that was used to
fit them [
The first modification is obtained by assuming polynomial, rather than linear form of N(x) in eq. (1), in particular, a quadratic dependence. The ith component Ni(x) could be written as
Ni(x) ≈
xT Aix + tix + ci(0) dt + dri(0)
(3)
The matrices Ai represent the blocks of a third-order
tensors, while the vectors bi(0) = li + ti are the rows of the
matrix B(0) = L + T (as in eq. (2)). These objects, as well
as components of the vector c(0), are estimated by multiple
polynomial regression [
The second modification, considering the serial correlations in residual forcing, is due to the multilevel structure of our model. In particular, consider the ith component of the first, main level of the inverse stochastic model dxi = xT Aix + bi(0) + ci(0) dt + dri(0), (4) where x = {xi} is the state vector and matrices Ai, vectors bi(0) and the components ci(0) of the vector c(0) as well as the components ri(0) of the residual forcing vector r(0) are determined by the least squares. The additional model level is added to express the known increments dr(0) as a linear function of an extended state vector [x, r(0)]. We estimate this level’s residual forcing again by the least squares. More levels are added the same way, until the Lth level’s residual, r(L+1), becomes white in time, and its lag-0 correlation matrix converges to constant, hence dri(0) dri(1)
. . . dri(L) = = = bi(1)[x, r(0)]dt + ri(1)dt, bi(2)[x, r(0), r(1)]dt + ri(2)dt, bi(L+1)[x, r(0), . . . , r(L)]dt + ri(L+1)dt (5)
The eqs. (4) and (5) describe a wide variety of processes
in a fashion that explicitly accounts for the modeled
process x feeding back on the noise statistics. The linear
multilevel model is obtained by assuming Ai ≡ 0 and c(0) ≡ 0
in eq. (4). Details of the methodology and further
discussion could be found in [
It is well known, that the extreme ENSO events tend to occur in boreal winter. From several ways to include this phase locking to the annual cycle, the alternative approach used here is to include seasonal dependence in the dynamical part of the first level. Namely, we assume the matrix B(0) and vector c(0) to be periodic, with period T = 12 months:
B(0) c(0) = =
B0 + Bs sin(2πt/T ) + Bc cos(2πt/T ),
c0 + cs sin(2πt/T ) + cc cos(2πt/T )
(6)
In this case, the whole record is used to estimate four
seasonal-dependent coefficients. The model is trained
in the leading EOF (empirical orthogonal function)
space [
In this section, the brief results are presented of how the statistical model is able to simulate the ENSO signal. The skill of the model is determined in the sense of basic linear ENSO metrics such as the amplitude of the ENSO signal, the seasonality (since the seasonality is important aspect of ENSO dynamics) and finally, the power spectrum of ENSO signal. The model is employed as described in the previous section, the matrices and vectors are estimated from the previous data and then the model is integrated to obtain the time series of same length as the training data. Since the model is stochastic (forced by a white noise), we employed an ensemble of 20 members. Each member is integrated with slightly different initial conditions and these members are referred to as realizations.
The basic ENSO metric is its amplitude, which could
be characterized by the standard deviation of SST
anomalies averaged over Nino3.4 box (bounded by 5◦S - 5◦N
and 120◦W - 170◦W). In Fig. 2 we can see the ENSO
amplitude as derived from the Nino3.4 index [
As can be seen, the linear model slightly overestimates
the ENSO amplitude, while the quadratic model slightly
underestimate the ENSO amplitude. From the spread of
the ensemble members we could infer that the model is
sensitive to initial conditions and the forcing. Still, the
ensemble averages for both models are within reasonable
distance from the data borderline, therefore in this aspect
the model performs adequately.
tive models are more flat in this area of frequencies. In
the higher frequencies (around annual frequency and less)
the power spectra are in agreement. In general, the spectra
of modelled time series could be said to copy the actual
Nino3.4 time series. The power spectra were computed
using the Welch method [
Other metric connected with ENSO amplitude is its seasonality. As written above, the ENSO phenomenon exhibits seasonal changes in variance, with elevated variance in winter months and lower variance in spring and summer months. This can be also seen in Fig. 3, where the monthly variance is plotted for the data and for both models. Both models are capable of modelling higher variance in winter months and drop in variance through spring and summer, although the difference in variance is higher in data than in both models. Still, the ensemble averages are reasonably close to the data.
The last metric taken into account was the power spectrum of Nino3.4 time series. The spectrum for the Nino3.4 data and both linear and quadratic model realizations can be seen in Fig. 4. The main peak in data occurs at roughly 5 year period, but still the ensemble averages for respec1.1 1.0 0.9 A T S foS0.8 D T S0.7 0.6 0.5 0.4 108 107 106 r e w o p 105 104 103 10-1 period [1/yr] 100
Noise parametrization in the model The statistical model, once estimated, is integrated forward in time and forced by a noise - usually a realization of spatially correlated random process. In the most intuitive and basic case, the last level residuals’ covariance matrix is estimated and decomposed using Cholesky factorization yielding a lower triangular matrix R. When the model is integrated, the random realization of white noise is multiplied by the matrix R, yielding spatially correlated white noise which is used as a random forcing in the model. The results for quadratic and linear ENSO models from the previous section were obtained using this simple noise parametrization, and the question is whether looking deeper into the residuals’ structure could aid the model’s performance. 4.1
First refinement for the noise parametrization arises from
the concept of modelling climate processes which exhibit
low-frequency variability (LFV). In this method, we find
and select noise samples, snippets, from the past noise
(residuals) which have forced the system during short time
intervals that resemble the LFV phase just preceding the
currently observed state, and then use these snippets (or
information contained in them) to drive the current state into
the future. For full methodology and discussion, see [
The found past noise snippets can be used in two
different ways. The first one (as used in [
As can be seen in Fig. 5, although the amplitude statistics are not substantially shifted, the transient from highvariance winter period to low-variance spring and summer are better captured by the later model, with noise forcing conditioned on system’s state. The power spectra for both models are practically the same (not shown). 4.2
Although the seasonal dependence of the model is captured in model’s dynamics by fitting the seasonally dependent matrices B(0) and c(0) (recall eq. (6)), our analysis showed, that the last level’s residuals still exhibit seasonally dependent amplitude. To address this issue, we computed the standard deviations for each month from the last level’s residuals, then fitted the 5 harmonics of the annual cycle to capture the seasonal dependence, removed this dependence from the residuals, then estimated covariance matrix and subsequently the matrix R and finally generated spatially correlated white noise realization which was multiplied back by the requisite seasonal amplitude to account for the seasonally dependent amplitude of the forcing. The fitted harmonics of the annual cycle were selected as
Pi = cos(2πit/T ) + sin(2πit/T ), i = 1, . . . , 5 (7) and then regressed on the seasonally varying standard deviation of the last level’s residuals. 4.3
The last modification to the noise is to use the extended covariance matrix instead of lag-0 covariance matrix. When evaluating system’s state we do not take just the state closest to the current state of the model, but, say 5 consecutive months and construct the extended matrix out of this snippet. Then the matrix is decomposed using Cholesky factorization and used as a spatial correlation matrix R is random forcing generation. 0.90 40.85 . 3 o n i-N A T fSS0.80 o TD S inearear1ear2ear3ear4ear5ear6ear7ear8ear9ear10ear11ear12ear13ear14ear15ear16ear17ear18ear19ear20 eanl lin lin lin lin lin lin lin lin lin lin lin lin lin lin lin lin lin lin lin lin m d. d1 d2 d3 d4 d5 d6 d7 d8 d9d10d11d12d13d14d15d16d17d18d19d20 eanconditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditioneditione mcon con con con con con con con concon con con con con con con con con con con Jan
Feb Mar
Apr
5 Synchronization and causality in the
observed and modelled data
Better understanding of the complex dynamics of the
atmosphere and climate is one of the challenges for
contemporary science. Considering the climate system as
a complex network of interacting subsystems [
Paluš [
, arctan qs2(t) + sˆ2(t).
Then the time series of phase and / or amplitude are used to study the interactions. We adopt measures from information theory, namely mutual information and conditional mutual information, where the mutual information could be expressed as
I(X ;Y ) = ∑ ∑ p(x, y) log x∈X y∈Y p(x, y) p(x)p(y) , where p(·) is the probability distribution or joint probability distribution and X and Y are our time series of either phase or amplitude derived from the ENSO SST data. Finally, the measures we are interested in could be written as: (8) (9) (10) (11) • phase synchronization – I(φ1(t); φ2(t)), • phase-phase
φ2(t)|φ2(t)), • phase-amplitude causality τ)|A2(t), A2(t − η), A2(t − 2η)), – causality –
I(φ1(t); φ2(t + τ) −
I(φ1(t); A2(t +
As can be seen from Fig. 6, in ENSO dynamics captured by the Nino3.4 index, the synchronization of annual cycle with quasi-biennal and combination frequencies (frequencies that arise from the interactions between annual and 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.
Our goal was to simulate the nonlinear cross-scale interactions in the model. This is important since it might help to uncover the mechanisms of these interactions and shed more light onto the dynamics of ENSO in general. We constructed the ENSO model and repeated the above analysis to modeled ensemble of the Nino3.4 time series.
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
captured by the modelled data, while the phase -
amplitude interactions are not very well captured. This might
arise from the low complexity of the model, or the absence
of some nonlinear interactions in the model design (apart
from quadratic).
Surrogate data (or analogous data) is a method to generate
synthetic data set (time series) that preserve some of the
statistical properties, while omitting the others. One way
of using them, is to test statistical significance by
contradiction. 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
generating surrogate data is the Fourier transform surrogate [
As an example, consider two intertwined Lorenz systems, where one of them drives the other. Now, using the time series in one dimension, say the x dimension from both Lorenz systems, we can use some method for detecting causality, e.g. conditional mutual information between the two time series of two Lorenz systems. We get the value of conditional mutual information, but this is still not enough to interpret it in the means of whether there 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 preserves properties of the time series, but allows no causal relationship between them) and repeat the analysis using the very same method on this ensemble and finally compare the value for actual data with the histogram of values obtained from the ensemble of surrogate data. When the value from the data exceeds some percentile (e.g. 95th) of the surrogate data distribution, we say that the causal relationship is significant in comparison with e.g. 500 FT surrogates.
When studying nonlinear cross-scale interactions in time series using the above method, the statistical test involves creating an ensemble of surrogate, synthetic time series and repeat the analysis for the whole ensemble. Then we computed the percentile, where the observed interactions could not arose by random chance. Of course, one could use Fourier transform method to generate the surrogate time series, effectively posing a null hypothesis of a linear process which has the same spectrum to that of an observed data. On the other hand, one can create a more sophisticated null hypothesis by exploiting the options of a data-based model: when one consider just a linear model, omit the dynamical seasonal dependence in B(0) and c(0) terms (as in eq. (6)) and use the simplest noise parametrisation (just consider the spatial covariance structure), the model will omit the nonlinear interactions and could pose as a surrogate data model copying the basic statistical properties of a modelled time series. This way, the analysis would show whether the cross-scale interactions are arising from the seasonal dependent dynamics, or from nonlinear (e.g. quadratic) interactions between subsystems and so on.
When comparing Fig. 6 (testing against 500 Fourier transform surrogates) and Fig. 8 (testing against 500 databased model surrogates), the significant interactions are virtually the same, expect in the latter, the “fluctuations” (or they might be false positives as well) are attenuated to minimum. This way, we can get better idea of the statistical significance of the interactions between subsystems, in particular the nonlinear ones, since we are testing against the model with just linear interactions. 7
Statistical modelling in climate science is continuously getting more attention, since their usage is not limited to forecast some of the phenomena of interest (like ENSO), but could also be used to infer some of the statistical properties and relationships among different subsystems. Since the statistical models live in phase space of particularly reduced dimensionality, when we could observe the interactions of interest, the identification of their sources will become more feasible.
We showed that the statistical model with the right settings, which were selected based on careful inspection of the modelled system, could generate synthetic time series of interest, copying the desired properties of the system both linear and nonlinear statistics. Since the stochasticity is the important aspect of the data-based model, various parametrization techniques exist to correctly model the system’s external forcing. Finally, the possibility of usage of the low complexity model as surrogate data was discussed, showing advantages of usage of such technique to infer statistical significance.
The outlook for future work combines various different paths which appeared. One direction would be focusing on statistical modelling itself, experimenting with various variable model, with input time series and their preprocessing and so on and so forth. Other direction would be connecting the statistical models with dynamical ones, in the sense, that statistical models could be used for parametrization of e.g. sub-grid phenomena (microphysics of clouds, local convection etc.) in large coupled atmospheric-oceanic models.