# Synergistic Dynamic Theory of Complex Coevolutionary Systems: Disentangling Nonlinear Spatiotemporal Controls on Precipitation

###### Abstract

We formulate a nonlinear synergistic theory of coevolutionary systems, disentangling and explaining dynamic complexity in terms of fundamental processes for optimised data analysis and dynamic model design: Dynamic Source Analysis (DSA). DSA provides a nonlinear dynamical basis for spatiotemporal datasets or dynamical models, eliminating redundancies and expressing the system in terms of the smallest number of fundamental processes and interactions without loss of information. This optimises model design in dynamical systems, expressing complex coevolution in simple synergistic terms, yielding physically meaningful spatial and temporal structures. These are extracted by spatiotemporal decomposition of nonlinearly interacting subspaces via the novel concept of a Spatiotemporal Coevolution Manifold. Physical consistency is ensured and mathematical ambiguities are avoided with fundamental principles on energy minimisation and entropy production. The relevance of DSA is illustrated by retrieving a non-redundant, synergistic set of nonlinear geophysical processes exerting control over precipitation in space and time over the Euro-Atlantic region. For that purpose, a nonlinear spatiotemporal basis is extracted from geopotential data fields, yielding two independent dynamic sources dominated respectively by meridional and zonal circulation gradients. These sources are decomposed into spatial and temporal structures corresponding to multiscale climate dynamics. The added value of nonlinear predictability is brought out in the geospatial evaluation and dynamic simulation of evolving precipitation distributions from the geophysical controls, using DSA-driven model building and implementation. The simulated precipitation is found to be in agreement with observational datasets, which they not only describe but also dynamically link and attribute in synergistic terms of the retrieved dynamic sources.

## 1 Introduction

### 1.1 Fundamental Motivation

The Earth System in general, and the Climate system in particular, are inherently coevolutionary and complex [Peixoto and Oort 1992, Nicolis and Nicolis 2007]. Therefore, their fundamental understanding calls for the development and implementation of methodologies that disentangle nonlinear processes and interactions at play, ultimately improving their predictability. These are the core motivations of this paper.

The concepts of coevolution and complexity span a diversity of views depending on the background and context in which they are raised, including in biology, geomorphology and hydro-climatology, e.g. [Ehrlich and Raven 1964, Whipple et al. 1999, Perdigão and Blöschl 2014, Troch et al. 2015]. From all different accounts and perspectives, a general definition can now be formulated with overarching generality, with broad interdisciplinary applications in mind.

A coevolutionary system consists of a multiscale dynamical system wherein the intervening processes are dynamically bound to evolve in an interactive manner across spatiotemporal scales. Whether such systems are complex depend on whether the system as a whole evolves in a manner that cannot be attributed to any individual or linearly combined factor, rather exhibiting features not present or explainable by any such factors alone: emergence.

Complexity thus stems neither from the variables themselves, nor from their linear combinations, but rather from their nonlinear interactions. If the intervening variables are independent from each other (i.e. do not have any dynamical mutual influence) and still interact to produce an emerging process, we are in the presence of a synergy. For instance, secondary waves in triadic wave resonance are synergic processes emerging from nonlinear interacting pairwise independent waves [Pires and Perdigão 2015] or oscillatory processes [Hocke and Kämpfer 2008].

A wise statement from the ancient Greek philosopher Aristotle in his Metaphysica compendia and highlighted in [Nicolis and Nicolis 2007] aptly makes the point on complexity: ”The whole is more than the sum of its parts”.

The essence of the general methodological problem resides in finding the simplest set of non-redundant and dynamically independent controls able to explain and predict the dynamics of a physical process of interest in a general coevolutionary system. Given a set of measurements or signals capturing its behaviour, it is of interest to find and quantitatively characterise dominant features that represent the fundamental behaviour of the system, as represented by statistical or dynamical properties and interactions. This is particularly relevant in high-dimensional cases such as those involving atmospheric variables, where the full description of the system is often beyond reach and the dynamics appear to be random. Albeit the apparent randomness, self-organised and synergistically sustained phenomena may emerge far from thermodynamic equilibrium [Haken 1983], leading to coherent group behaviour that can be described by a simple set of macro-scale variables.

Figure 1 depicts the nature of the problem in visually intuitive terms: (a) a complex coevolutionary system is represented by a broadband spectrum of gray shades with a large degree of redundancy as expressed by shading codependence; (b) a lower-order palette enables the reproduction of the original system using a few key components, whilst retaining some redundancy among shades; (c) the simplest set of non-redundant components consists of just two tones, black and white, the fundamental basis whose members can synergistically interact to generate the entire gray scale representing the original system. This is, in a nutshell, the fundamental set of basis components that we seek.

Bearing in mind the notion of a vector or functional space representing all possible states of the system, our goal is to find a basis for that space, i.e. the fundamental independent ”axes” relative to which the entire dynamical process unfolds. As a functional basis, its components are independent from each other and their synergistic combinations generate the entire state space, i.e. the problem is reduced to the minimum possible dimensions without information loss.

In linear algebra and functional analysis, vector and functional spaces are traditionally generated by linear expansion of their basis elements [Reed and Simon 1980, Roman 2005]. Our interest, however, is in the more general case where the functional space is a spatiotemporal field fully generated by nonlinear combinations of an even lower-dimensional basis without loss of generality. For that purpose, we seek basis elements that are dynamically independent not only in the linear but also in the nonlinear sense. From within the multiplicity of mathematical solutions to this nonlinear problem, we further seek those actually representing fundamental physical processes. For that purpose, the sought basis terms are those for which first principles hold, namely on energy and entropy dynamics.

All in all, we seek the smallest set of physically consistent, non-redundant dynamic processes with full generative power enabled by nonlinear dynamic interactions (i.e. a nonlinear synergistic basis), providing complex systems with an optimal model design in terms of fundamental processes and interactions.

### 1.2 Applied Motivation

While a general methodological undertaking is our main focus, a second, no less important motivation drives our quest: a dynamical understanding of Precipitation and its predictability in terms of underlying fundamental controls.

Precipitation is one of the most widely studied phenomena playing a crucial role in hydro-climate dynamics. Still, its thorough understanding and prediction remains rather elusive due to its dynamical complexity, stemming from nonlinear interactions among intervening processes.

In Statistical Climatology and Hydrology, Precipitation is often treated as a random process, with each event being a particular realisation of a theoretically or empirically inspired statistical distribution, conditioned or not by signatures of assumed controls in the causal chain. Precipitation datasets are then often simulated as stationary or ciclo-stationary stochastic processes, [Wilks and Wilby 1999, Paschalis et al. 2013], or as random variables conditioned to relevant oceanic-atmospheric controls by Bayesian probabilistic models [Kidson and Thompson 1998, Bárdossy and Pegram 2011, Song et al. 2014], including nonlinear statistical downscaling approaches [Pires and Perdigão 2007, Perdigão 2010].

Still in the statistical context, Precipitation datasets can also be studied with complex process network approaches [Donges et al. 2009, Ruddell and Kumar 2009, Boers et al. 2013], wherein intervening processes are interpreted as network nodes, and their relationships as network branches. In practice, the branches are simply statistical codependence measures ranging from traditional linear or nonlinear correlations to information-theoretical diagnostics. Therefore, while containing statistical value, they do not necessarily correspond to interactions in the physical sense.

Albeit the inherent complexity of Precipitation, its stochastic treatment can also be complemented by deterministic approaches of time series analysis, wherein self-similarity and scaling properties are brought out in connection with geometric features of the dynamics. The Fractal-Multifractal method [Puente 1992], while inspired in the stochastic setting of multiscaling cascades in turbulence [Mandelbrot 1974], is a deterministic-geometric approach already used for Precipitation studies [Obregón et al. 2002, Maskey et al. 2015].

Even though data-based analytics are aptly descriptive and quite practical to implement, they may leave out relevant physical information that is crucial to a better understanding and prediction.

### 1.3 Traditional and Information-Theoretical Feature Extraction Methods

While statistical feature extraction has a long tradition in the geosciences, the fundamental dynamical analysis we seek has been more elusive.

Among various statistical techniques with diverse levels of sophistication, one of the most popular approaches to the extraction of spatial and/or temporal patterns is Principal Component Analysis (PCA), which has been widely used to extract dominant modes of spatial and temporal variability in atmospheric and oceanic fields [Horel 1981, Wallace and Gutzler 1981, Karl et al. 1982, Barnston and Livezey 1987, Preisendorfer et al. 1988]. Fundamentally, PCA searches for uncorrelated components under a given inner product, while maximising the amount of variance they explain. The signal is then given as a linear combination of those components.

However, when the processes of interest have non-Gaussian or non-normal distributions, PCA and similar methods will give out uncorrelated principal components (PCs) that are not truly statistically independent. In fact, those components may be further decomposed as a nonlinear combination of a smaller number of underlying independent processes still able to generate the full signal. This problem is common to all factor analysis approaches in which the signal decomposition is done into uncorrelated components. It is worth reiterating that null correlation does not imply independence.

Statistical tools do exist to search for a non-redundant set of statistically independent features within a non-normal signal. One of the most notable examples comes from the Information-Theoretical framework on Independent Component Analysis (ICA) [Comon 1994, Hyvärinen et al. 2001] and equivalent methods. ICA searches for statistically independent components in the general case when the data is non-normally distributed. It does so by transforming a multidimensional dataset into a linear combination of components with minimised statistical dependence or Mutual Information from each other, i.e., as close to statistical independence from each other as possible.

Statistical independence means that the value of any of the components gives no statistical information on the values of other components, and therefore the statistical predictability potential of one variable from another one is null. Still, there will always be a small residual codependence among the components [Cover and Thomas 1991], with the theoretical foundations for a strictly positive lower bound for Mutual Information associated to prescribed cross moments having been derived by [Pires and Perdigão 2012] independently of data size constraints, and by [Pires and Perdigão 2013] for finite-sized data samples.

Within the geophysical sciences, ICA has been applied to the Sea Surface Temperature (SST) field extracting statistically independent modes of variability in [Aires et al. 1999] and an intuitive didactic approach to ICA applied to geophysical-like correlated data was provided in [Aires et al. 2002]. Moreover, the method was also used by [Fodor and Kamath 2003] to separate dominant statistical features in climate data, in [Basak et al. 2004] for weather data mining using the North Atlantic Oscillation as a specific example, and by [Perdigão 2004] to produce statistically independent spatiotemporal teleconnection patterns from non-normally distributed atmospheric fields, using then those components as predictors for meteorological regimes.

For all its merits, ICA still has some limitations: its basic form does not disentangle nonlinearly mixed signals and its results are ambiguous in some situations, such as under rotational degeneracy when multi-dimensional subspaces are jointly Normally distributed, or when nonlinear mixing comes into play.

These limitations may be addressed to some extent by considering the full signal as a mixing of multivariate non-normally distributed sources through the Independent Subspace Analysis [Pires and Ribeiro 2016]. Alternatively, nonlinear generalisations of ICA and PCA can be considered [Pajunen et al. 1996, Oja 1997], along with advanced methods for nonlinear statistical decomposition of kinematic features such as Principal Interacting Patterns (PIP) and Dynamic Mode Decomposition (DMD) [Hasselmann 1988, Tu et al. 2014].

However, the sophistication of these approaches often comes at the expense of the uniqueness and fundamental meaning of the solutions. Moreover, whilst capturing nonlinear statistical dependence among system features on an aggregate level, they do not capture the dynamical dependence within the spatiotemporal domain of analysis, i.e. the small scale behaviour beneath the overall large scale features captured by the statistics. This is a common aspect of all statistical techniques involving blind source separation [Yu et al. 2014].

The core of the problem is that processes that are statistically independent can still be dynamically codependent [Perdigão 2016]. By leaving out the fine print in the dynamics, the aggregate diagnostics have a more limited predictive potential. Even if in the end an aggregation is done on dynamical diagnostics, at least their fundamental structure will have ”learnt” (gathered information) from the dynamics rather than solely from the aggregate statistical features.

### 1.4 Study Outline

One of the key aims of the present study is to disentangle nonlinearly interacting processes in the climate system, in order to bring out fundamental mechanisms that, whilst dynamically independent from each other, cooperatively influence precipitation regimes. For that purpose, this paper introduces a physical source decomposition methodology in section 2, where the fundamentals of dynamic source analysis are formulated. A framework for nonlinear space-time decomposition in general spatiotemporally coevolving fields is then introduced in section 3. The theoretical foundations for evaluation of predictability and dynamic model building are then introduced in section 4. The methodologies introduced in the paper are essentially nonlinear dynamic analysis and model building techniques with broad applications. The application of the novel methodologies to the retrieval of nonlinear atmospheric controls of precipitation is conducted in section 5, followed by an evaluation of dynamic predictability and simulation of evolving precipitation distributions given the retrieved controls, in section 6. The main body of the paper is then completed with an overall discussion and concluding remarks in section 7.

Further to the main body of the paper, three appendices are included. Appendix A addresses the concepts of coevolution and synergy in information-theoretical terms. Appendices B and C introduce new theoretical developments in functional analysis necessary for the spatiotemporal decomposition conducted in the dynamic source analysis framework. These methodologies enable an effective decomposition of nonlinearly interacting functionals (e.g. the spatial and temporal information) even in a coevolutionary setting whereby space and time are nonlinearly entangled by relative celerities.

## 2 Dynamic Source Analysis: a Synergistic Theory of Coevolutionary Systems

### 2.1 Motivation and Fundamentals

#### 2.1.1 Expressing Redundant Coevolution as Non-Redundant Synergy

Coevolutionary complex systems involve codependent observables, i.e. mensurable quantities with some degree of redundancy. This is indeed the case in the current paradigms of coupled dynamical systems. The existence of redundant information suggests that a more efficient formulation could be derived involving a lower number of fundamental non-redundant terms without loss of generality. The fundamental terms are the independent dynamic sources we seek, representing the fundamental ”backbone” of the dynamical system.

While coevolution always entails a degree of redundancy among processes, a synergistic interaction does not. In fact, it rather entails cooperative dynamics wherein the intervening processes jointly produce emerging structures with features that do not exist in any individual or linearly combined source components. These matters are mathematically exemplified in Appendix A.

A purely synergistic interaction is optimally efficient since it leads to the emergence of unprecedented features (innovation), giving the ability of a simple set of processes to synergistically produce an emergent cascade of child processes that ultimately leads to the full observed complexity [Pires and Perdigão 2015]. The polyadic wave resonance in fluids is a simple example of synergistic interactions among parent processes ( generation waves) leading to the emergence of child processes [ generation waves]. This is the case in geophysical phenomena such as nonlinear triadic and quartic interactions among ocean waves [Komen et al. 1996].

Essentially, our idea is to express coevolutionary complex systems involving a large number of codependent observables in terms of synergistic interactions among independent components representing the fundamental underlying processes. For that purpose, we seek a simple nonlinear dynamic basis for complex nonlinear systems.

#### 2.1.2 Unveiling a Simple Nonlinear Basis for Complex Systems

Let denote the set of observables, i.e. the physical variables or fields that are observed, measured and recorded in spatiotemporal datasets, where represents the spatial coordinates and the time. Let denote the set of the hereafter called Dynamic Sources, taken as spatiotemporal functions from a Hilbert space [Courant and Hilbert 1953] where differentials can be defined: a Sobolev space [Sobolev 1938, Sobolev 1963] comprising a class of (square integrable) functions with weak derivatives up to order . Weak differentiability means that differential operators can be defined everywhere in the domain even if traditional differentiability fails in a null-measure set (e.g. a 2D set in a 3D domain). This is particularly useful when evaluating fluid flows over critical transitions, e.g. across weather fronts [Perdigão 2016].

Essentially, the dynamic sources are fundamental spatiotemporal state variables from which depend via an operational set of dynamical relationships in the general form:

(0) |

The operator maps the -dimensional functional domain of the sources, to the -dimensional functional domain of the observables, , both of which are subsets of . In practice, can be seen as an observation operator akin to those involving spatial re-colocation, temporal delaying, differentiation or integration. Without a priori constraints to the spatial and temporal domains of sources and observations, they are not necessarily the same. Moreover, the spatial coordinates may be time-dependent in a coevolutionary setting, leading to a celerity term of the space referential [Perdigão and Blöschl 2014].

The aforementioned variables are generic -dimensional functions living in an -dimensional functional space, with . The dynamical system in equation (2.1.2) then prescribes the dynamics of a spatiotemporally distributed quantity, which can be straightforwardly treated as a multidimensional deterministic field or a multivariate stochastic variable or distribution.

Our goal is to retrieve the dynamic sources from the observational datasets that minimise the dynamic dependence (linear and non-linear) among the components of and maximise their ability to [not necessarily linearly but also nonlinearly] generate the entire state space spun by . In other words, given the observations (the observed dynamics) we seek a functional basis in their state space. That basis shall provide the optimal components underlying the observations, since it represents them with neither redundancies nor information loss.

The retrieval of from would be trivial if their functional relationship would be a known invertible function, by simply taking . However, that is not necessarily the case in most applications. Without further knowledge, the problem would be under-determinate. Therefore, constraints could be put into place, as done statistically with ICA and PCA, so that ambiguities would be overcome. In the present study, the constraints are of physical nature (section 2.5) and the feature extraction performed over the dynamics (sections 2.2 - 2.4). For notational ease, unless otherwise specified, symbols in uppercase bold type shall refer to spatiotemporal variables, i.e. and in the upcoming sections.

The key differences between the proposed Dynamic Source Analysis (DSA) and existing families of methods are summarised in Figure 2. Essentially, traditional feature extraction methods yield uncorrelated but not statistically independent components; information-theoretical methods do provide statistically independent components but not dynamically independent ones; and DSA provides dynamically independent components, thus yielding the lowest level of redundancy in the characterisation of the system.

### 2.2 Nonlinear Diagonalisation of Dynamic Interactions

Dynamic Source Analysis seeks dynamically independent sources, i.e. such that no redundant dynamic interactions exist among them. This brings out a canonic formulation of the dynamical system in terms of non-redundant synergistic interactions, which is naturally more efficient than formulating the dynamics in redundant, coevolutionary terms.

For this purpose, we define Dynamic Interaction of order among the components of as:

(0) |

where is the -order tensor differential operator given by the -times recursive application of the gradient with respect to the spatial components of , and is the time derivative of . In the spectral case where each spatiotemporal component corresponds to a different scale , equation (2.2) provides the Dynamic Scale Interaction of order within .

The dynamically independent sources are thus those for which the following condition holds:

(0) |

where , with , is a functional basis to the dynamical system prescribed by . As such, the off-diagonal elements of vanish, meaning that the cross-dependencies are null not only at bilateral but also multilateral levels. In the linear case, for and corresponds to a diagonal Jacobian matrix of the underlying dynamical system, i.e. null linear couplings. In the general nonlinear case, it means that along the whole trajectory in state space all non-diagonal (i.e. cross) dynamical interactions are null.

The derivatives in equation (2.2) require the knowledge of the phase space trajectory prescribed by . However, that information can be found from the observable dynamics . In fact, since both and come from the same dynamical system, they have a topologically equivalent phase space portrait except if in equation (2.1.2) is under-determined from .

In practice, the dynamical systems of interest shall be generic spatiotemporal fields spanning phase space manifolds. These are topological spaces that, whilst eventually curved, can be mapped, at the vicinity of each point, to and from Euclidean spaces via a functional bijection (invertible one-on-one correspondence or isomorphism). This is the case with fundamental geodesic and relativistic spatiotemporal geometries with broad applications.

One of the most relevant isomorphisms is the transformation between non-Euclidean manifolds and the corresponding Euclidean tangent spaces at the vicinity of each point, which under appropriate differentiability conditions defines a diffeomorphism [Da Silva 2001]. Therefore, albeit the curvature of the spatiotemporal field of study, the topological space spun by its spatiotemporal gradients is conveniently comprised of an Euclidean structure (tangent bundle, union of all local tangent spaces).

A crude motivational image of a tangent bundle can be obtained by depicting a ball covered by flat ”postit” stickers. The ball represents the manifold, each sticker is a local tangent space, and the family of stickers covering the entire ball is the tangent bundle.

The rigorous mathematical treatment of manifolds and associated concepts stems from differential geometry, whose roots date centuries back to the origins of cartography [Spivak 2005].

### 2.3 Beyond Statistical Independence: Differential Geometric Quest for Dynamic Independence

Essentially, DSA can be seen as a dynamic generalisation of ICA wherein the independent components are computed at the tangent space to each point and fundamental physical constraints on energy and entropy dynamics (section 2.5) are taken into consideration. The procedure is conducted throughout the phase space manifold to yield the global dynamic sources. Given the differential geometric isomorphism between local charts and tangent spaces, the independent components computed on the tangent spaces are applicable to the actual manifold chart as well. The global dynamic sources are thus defined in the tangent bundle guiding the underlying phase space manifold.

With this procedure, DSA ensures independence among components at every stage of the dynamics (dynamic independence), rather than only over a statistical aggregate (statistical independence as in ICA and information theory). Moreover, the dynamic information is preserved as the analysis follows the phase space trajectories rather than assigning classes with disregard to the dynamic sequence (as happens when defining a statistical distribution e.g. by assigning bins or classes).

Albeit the technical sophistication of differential geometry, it should be noted that its mathematical essence is actually rather intuitive if we think in geodesic terms. In all quantitative geosciences, local coordinates (e.g. zonal, meridional) are also defined on the local tangent space to the curved geophysical surface of interest at each point. While the local coordinate system evolves from point to point due to the curvature of the surface, a global set of coordinates is then defined with overarching span and simplicity, e.g. the spherical coordinates. Similarly, the dynamic sources are fundamentally global functions that encompass, in a similarly simple manner, the entire phase space manifold. As such, they define a basis grounded on a canonic set of generalised curvilinear coordinates. The practical challenge thus lies in finding such a functional basis for generic dynamical systems.

### 2.4 Dynamically Independent Sources from Coevolving Observables

The problem consists essentially on expressing a complex coevolutionary dynamical system in Taylor expansion terms of a simpler synergistic system of independent dynamic sources :

(0) |

where refers to evaluation at the reference manifold^{2}^{2}2The reference manifold is spun by the invariants of motion and includes not only attractor but also repeller sub manifolds., and

(0) |

where are the spatiotemporal fields characterising the dynamic sources ( as defined in Section 2.1.2).

By applying the condition (2.2) to equation (2.4), the dynamically independent sources sought after are ultimately given by:

(0) |

In practice, a functional diagonalisation is performed to the coevolutionary system towards meeting the condition in (2.2) and thus leading to (2.4). For that purpose, DSA searches for that minimises the off-diagonal functionals defined by

(0) |

with , i.e. the quadratic norm or self inner product of the difference between the left- and right-hand sides of the dynamical independence condition (2.2).

The quadratic norm in equation (2.4) is taken in the Sobolev space with the inner product generically defined by:

(0) | ||||

(0) |

where are generic square-integrable tensor functions, is the order up to which differentials are defined in the Sobolev space, is the order of the differentiation operator , is an unbounded manifold living in , and * denotes complex conjugate transpose. The measure refers in general to any mathematically consistent metric in .

Note that, while coevolutionary dynamics unfold among the components of , no such coevolution occurs among different independent dynamic sources in . In fact, by construction the components interact in a cooperative synergistic manner in the dynamics of whilst retaining their own dynamic independence relative to each other.

### 2.5 Disambiguation by Physical Principles

Eventual ambiguities in the solutions stemming from non-linearities (e.g. solutions of with the same quadratic norm) are overcome by taking the physical solution grounded on first principles, namely the dynamic configuration with maximum entropy production and minimum energy consumption. Entropy and energy are quantified in the mechanistic sense, which enables their characterisation in dynamical systems terms by resorting to topologic properties of the phase space representing the dynamics. In practice, the entropy measure is computed with the topologic entropy [Ott 2002] and the energy as a Morse-Lyapunov function [Grines and Pochinka 2010].

We caution the reader that purely dynamical systems approaches are only kinematic or motion-descriptive, until adequately taking physical principles into consideration. It is important to take into consideration that the fundamental entropy production principles hold at constant energy levels, and energy minimisation principles at constant entropies [Lage 1995]. Therefore, our fundamental physical constraints consist of seeking solutions with minimum isentropic energy and maximum equienergetic entropy production. That is, entropy production is evaluated in the equienergetic subspace and energy dynamics evaluated in the isentropic subspace of the phase space spun by the dynamical system. The joint implementation of these principles brings out the general disambiguated sources over the overall phase space. Moreover, it ensures the physical consistency of the solutions.

Formally, we define the Physically Optimised System among the functional of mathematical possibilities for equienergetic entropy production rates and isentropic energy consumption rates , as the concrete function within for which the equienergetic entropy production rate is maximum and isentropic energy consumption rates is minimum:

(0) |

where the Delta functional extracts the concrete characterisation of the function at .

By applying this definition to the set of mathematical solutions , the physically optimised system becomes:

(0) |

Having retrieved dynamically independent sources from the nonlinearly mixed observables captured by the measurements, it is important to ensure the invariance of the information content within the system, i.e. that no information has been lost or spuriously added in the procedure.

In the dynamic source analysis, the dynamics of the measurements and of their corresponding independent sources are related by a locally smooth homeomorphism, equation (2.4). The physical reason for this is that their respective dynamical systems, and , span the same phase space manifold and invariants of motion, as noted above. Therefore, they also share the same information content, which thus ensures the non-existence of spurious or lost information in the analysis. That is, the referential may change but the underlying physics do not.

## 3 Nonlinear Space-Time Decomposition and the Spatiotemporal Coevolution Manifold

### 3.1 Motus: Expressing Spatiotemporal Processes in Observable Spatial and Temporal Structures

Having retrieved dynamically independent processes, it is also important to contextualise, characterise and depict them in terms of perceptual references such as space and time. In that regard, we will be interested in their spatial structure captured over some time span (e.g. a Spatial Structure over a Climatology); and on the temporal structure captured over some spatial domain (e.g. a Temporal Structure over a Region). For that purpose, observable spatial and temporal structures are sought for the dynamic sources .

A standard procedure in space-time decomposition consists of taking subspace projections of the spatiotemporal process by statistically aggregating it over the orthogonal complement of the subspace of interest. For instance, in a linearly separable space-time coordinate system, spatial patterns can be obtained by averaging the spatiotemporal process over time, and aggregate time series by averaging the process over space. A spatiotemporal process can then be expressed in terms of its temporally and spatially distributed components, e.g. Principal Components (time) and Empirical Orthogonal Functions (space) in Principal Component Analysis. Given their statistical nature, they assist in the study of its temporal and spatial variability.

However, in the aforementioned procedure the dynamics are lost due to the averaging operations. Moreover, spatial and temporal subspaces are not necessarily independent. In fact, the spatial reference can be time-dependent, as expressed via spatiotemporal coevolution [Perdigão and Blöschl 2014].

### 3.2 Decomposition over Interacting Subdomains

In this section, a dynamically preserving methodology is introduced for space-time decomposition of spatiotemporal processes in the general case where spatial and temporal components of the reference frame are not necessarily independent from each other.

Consider a generic functional basis which can be the independent dynamic sources retrieved in the previous section. Whilst independent from each other, the components of are all spatiotemporal dynamical processes. Our aim is to characterise these spatiotemporal components into observable spatial and temporal functionals.

For that purpose, we introduce an Interacting Subspace Decomposition Operator that, once applied to , extracts its spatial and temporal subspace manifolds (sub manifolds) for respectively. The general mathematical definition of the operator is introduced in Appendix B, and its application to is given by:

(0) |

where we introduce a Retrieval Product retrieving the -structure or dimension from . Its general definition is also introduced in the Appendix B. Still in equation (3.2), is the normalised functional basis for the subspace , and is a tensor functional with the rank of the projected subspace. For instance, and for a tridimensional space and unidimensional time.

The spatiotemporal decomposition in (3.2) is performed whilst concentrating all the space-time codependencies into a generalised coevolution term, consisting of the interactions between spatial and temporal submanifolds. In order to mathematically characterise this term, we introduce the concept of a Spatiotemporal Coevolution Manifold, hereby defined as:

(0) |

where represents the inner product defined in equation (2.4), involving multiple-order interactions. By noting that the gradient of a generic coordinate is given by , the first-order gradient differential form of the spatiotemporal coevolution is given in terms of relative celerities: . By definition of gradient, is locally orthogonal to .

In spatiotemporally separable systems, as there is no shared dimension between space and time. In non-separable systems (in which spatial patterns are time dependent), the tensor rank of characterises the dimensionality of the spatiotemporal coevolution manifold.

A simple example is now given for a propagating wave on a coevolutionary spatiotemporal coordinate system , with and as constants. Here, the coevolution manifold is the 1D manifold given by the straight line defined by and the magnitude of the spatiotemporal interaction is given by the relative celerity . The canonic form of that wave in non-redundant coordinates is thus or, converserly, .

By taking the definition of given in the Appendix B [equation (B.1)], the retrieval operation in equation (3.2) can be written in terms of contraction^{3}^{3}3The contraction is an interior product when a tensor operates onto a differential form (such as , which is essentially a gradient). When the aforementioned differential form has rank one, the interior product corresponds to a traditional inner product.
() and outer () products as:

(0) |

where is the orthogonal complement of , depending on alone.

When space and time are independent, and vice-versa. Here, however, space and time are related by the Spatiotemporal Coevolution Manifold , the presence of which means that space and time are not orthogonal complements of each other in the functional sense.

Essentially, the operator performs a contraction or dimensionality reduction in the functional manifold into a submanifold in a defined subspace, say space or time. In practice, and correspond, respectively, to the spatial and temporal structures of the spatiotemporal process .

Rebuilding the spatiotemporal process from its spatial and temporal submanifolds can then be done by introducing an Interacting Subspace Composition Operator and an associated Composition Product :

(0) |

The general mathematical definitions of the composition operator and product are introduced in the Appendix C. Their direct application to the spatiotemporal composition in equation (3.2) is given by:

(0) |

where , the aforementioned coevolution manifold, defines the rank- subspace of spatiotemporal interactions, i.e. where and interact. In practice, it also refers to the relative rate of structural regional and climatological changes, i.e. in space and time, reflecting the dynamic spatiotemporal codependence in the process .

All in all, is retrieved from its spatial and temporal subspaces and by taking their rank-additive outer product and contracting with their interacting subspace spun by . The total rank or dimensionality of is thus that of space plus time minus the rank of the space-time interactions.

Unlike subspace decompositions based solely on projections, in our formulation the dynamics are stored in the coevolution manifold, enabling the recomposition of the spatiotemporal process from its spatial and temporal structures. In a non-coevolutionary setting, wherein dim, the decomposition and composition operators reduce to the usual projection and composition operations involving linearly separable subspaces.

## 4 Evaluation of Dynamical Predictability

### 4.1 Dynamic Interaction Analysis

Having retrieved with DSA a functional basis for the dynamical system of interest, it can now be used to evaluate its predictive power with regards to other processes in which dynamics the system may play an active part.

In this regard, simple linear or nonlinear correlation measures do not suffice, as they provide an aggregate statistical relation without any considerations of dynamic nature. Moreover, they leave out the information content stemming from higher-order nonlinear dynamical interactions.

In order to capture the dynamics at play between processes, we hereby introduce a Dynamic Interaction Analysis (DIA). The procedure basically consists of evaluating the linear and nonlinear dynamical sensitivities between processes, i.e. their dynamic interactions, along with their evolution in time. This will be a fundamental building block for dynamic model design (section 4.2).

Formally, the Dynamic Interaction of order among spatiotemporal multivariate processes , is hereby defined as:

(0) |

This generalises Equation (2.2) for multiple multivariate processes. In particular, the first order dynamic interaction between two processes living in a vector space is given by the Jacobian matrix of their coupled dynamical system.

By evaluating all self and cross sensitivities of various orders within the system, the measure in equation (4.1) captures linear and nonlinear dynamic interactions required to characterise the dynamic feedback structure of the system.

As with the Dynamic Sources, we are also interested in the spatial and temporal structures of the dynamic interactions between processes.

The spatial structure of the interaction of order between processes , captured over a time period can then be expressed by the spatial submanifold of their dynamic interactions, given for each order by:

(0) |

where is as in equation (3.2) and its operand given by equation (4.1).

Analogously for time, the temporal structure captured over a certain spatial region is given by:

(0) |

The Spatiotemporal Dynamic Interaction of order among the processes is then retrievable from their spatial and temporal structures with the -product introduced in equation (3.2):

(0) |

### 4.2 Dynamic Model Design

Having retrieved independent dynamic sources (e.g. atmospheric controls) underneath a coevolutionary dynamical system (e.g. geopotential height fields), the goal of the present section is to establish a framework for the simulation of predictand processes (e.g. precipitation) controlled by . This can be addressed as a simple nonlinear dynamic model, whereby the processes of interest are written as a dynamic function of the sources :

(0) |

where denotes evaluation at the reference manifold, the dynamic sources are governed by equation (2.4), and refers to higher-order dynamics beyond the explicit analytical truncation at order . The order powers account for nonlinear interactions of order within each component of but not between different components or sources due to the disentanglement condition in equation (2.2).

Equation (4.2) provides a general model structure template or recipe to address the dynamical evolution of any system in terms of a functional basis comprised of underlying independent controls. The determination of a functional basis in a dynamical system is an important step to its optimal modelling, as it ensures that the system is expressed in terms of the smallest number of dimensions without loss of generality. By taking into account not only linear elasticities ( for ) but also nonlinear codependencies ( for ), this equation can then address changes in dynamical tendencies in terms of nonlinearly interacting processes.

## 5 Dynamic Source Analysis of Nonlinear Geophysical Controls

### 5.1 Characterisation of the Geophysical Observable and Associated Datasets

Our quest for nonlinearly interacting controls in the hydroclimatic system begins with one of the most popular geophysical observables that embodies the dynamical imprint of such interacting processes: the spatiotemporal Geopotential Height fields at pressure levels , .

In a nutshell, provides information on the geospatial dynamics of atmospheric waves at the pressure level to which refers, e.g. hPa. High and low pressure centres can be detected in geopotential terms as high and low heights of a given pressure level. In particular, the large-scale fields provide information on low-frequency atmospheric waves at the planetary scale.

In their natural form, contains not only information on the state of the atmosphere at the pressure level but also on how it is being affected by underlying forcing mechanisms taking part in its dynamics. In fact, embedded at the fields for each pressure level lie dynamical footprints of processes at other levels, through vertical interactions that ultimately bring in features stemming from other geophysical domains. These include oceanic and land-surface processes at large scales, resulting in geopotential heights actually responding to the difference between the state of the atmosphere over large water and land bodies such as oceans and continents, even considering mid-upper pressure levels such as 500 hPa. Higher pressure levels (corresponding to lower vertical levels) will naturally have a finer sensitivity to the surface properties, at the expense of getting lost in the details and losing the grasp on the bigger picture.

The geophysical mechanisms captured in the geopotential height fields are in turn at the source of the atmospheric circulation regimes that ultimately shape the weather, including moisture advection and thermodynamics, underlying the precipitation dynamics we intend to understand. Comprehensive, introductory treatments on associated atmospheric physical processes are found at [Salby 2012] and [Holton 2004].

One of the most reliable and popular datasets with geopotential height fields stems from the NCEP-NCAR Reanalysis Project [Kalnay et al. 1996]. Therefore, it is the geopotential dataset of choice for our dynamical source analysis. In the present study, the dynamic source analysis methodology presented in section 2 is implemented onto NCEP-NCAR Reanalysis data on the 500 hPa Geopotential Height field for a subset of the Northern Hemisphere spanning [30,70] degrees of latitude North and [-80,40] degrees East of longitude. This area covers most of Europe, the Mediterranean and the North Atlantic.

Our interest is in dominant atmospheric processes at the regional and monthly to multidecadal scales. Therefore, we consider datasets with an adequately coarse spatial resolution ( degree grid-cells), and 30-day moving average of Geopotential Height daily datasets, over a multidecadal period ranging from 1951 until 2012.

### 5.2 Statistical Physics of the Geophysical Controls

In order to evaluate whether traditional statistical tools such as PCA, optimal only for normally distributed linearly mixed data, would be appropriate to the analysis of Geopotential Height fields, we diagnose the information-theoretical measure of Negentropy [Comon 1994, Hyvärinen et al. 2001]. This is a non-negative measure, given by the information entropy deficit of a generic multivariate signal relative to that of a multivariate normally distributed random variable with the same mean and standard deviation as , i.e.:

(0) |

where denotes information entropy and negentropy in the information sense [Shannon 1948]. Higher values of correspond to lower statistical dependence among the components of .

In the present study, Negentropy is computed with the recently developed nonlinear statistical estimator of mutual information by anamorphosis [Pires and Perdigão 2012, Pires and Perdigão 2013] for better robustness against outliers and data sparsity, noting that information entropy is the mutual information between a variable and itself. The procedure consists essentially on performing a bijective, information-invariant, differential-geometric transformation between the original state space and a working space for optimal statistical analysis without loss of information. Recent applications of the aforementioned method include the estimation of nonlinear information codependences among coevolving landscape-climate processes [Perdigão and Blöschl 2014] and the development of an information-theoretic framework for nonlinear scale interactions among chaotic processes in complex systems [Pires and Perdigão 2015].

Figure 3 shows the Negentropy of the 500 hPa Geopotential Height field, with a significance level at 95% of 0.01 nat (10 mnat). Nat is the unit in which information entropy is expressed when natural logarithms are taken in its computation. For most of the covered area Negentropy is higher than 10 mnat, which indicates a significant departure from Normality, especially over the Atlantic, where Negentropy reaches its maxima.

Clearly, the Geopotential Height fields are not normally distributed, which indicates that statistical tools assuming or designed for such distribution are not an optimal way to address these fields.

Naturally though, it might be argued that data could be pre-whitened. However, that would remove precisely the kind of information that is most interesting in the Geopotential Height fields, namely that it is actually not a ”free” variable but rather one under the influence of external processes. In fact, by exhibiting negentropy, their distribution exhibits statistical signatures of underlying controls. This statement stems directly from the definition of Negentropy as entropy reduction with respect to the unbounded entropy maximum corresponding to the Normal distribution with the same mean and variance. The definition of negentropy further dissipates any misconception about information Entropy as being akin to variance: it is not, as distributions with exactly the same variance can still have completely different information entropies, as evaluated in the negentropy measure.

Higher negentropies diagnose higher entropy reductions, which suggest the existence of stronger controls being enforced on the geopotential over the areas where that reduction is strongest. The geopotential thus serves as proxy for such underlying controls. For instance, while it affects wind patterns, it is also affected by them, which can be understood by bearing in mind that just as a gradient in potential energy drives flow, it is progressively depleted by it as the system works towards geopotential balance. Moreover, albeit the evaluation at a certain pressure level, the horizontal heterogeneity of the geopotential height negentropy informs on entropy constraints imposed by vertical dynamics, including vertical moisture and heat exchanges.

Note that negentropy is higher over the subtropical Atlantic than over Europe, suggesting a stronger oceanic contribution towards constraining the Geopotential even at levels where the earth surface would no longer be expected to play a significant role. The aforementioned contrast is more visible particularly at lower latitudes, where the sea is exerting more influence via enthalpy forcing of the atmospheric pressure gradients and subsequent circulation. This may actually bestow more predictability onto the geopotential, as it is endowed with signatures from oceanic processes with longer memory and driving potential than short-term footprints from intrinsic atmospheric dynamics. Further investigation on these matters is deferred to a subsequent study.

The evidence of regions with a significant non-Normal signature and the associated statistical physics elicit the importance of treating the Geopotential height fields with methodologies that take into account the non-normality of the associated atmospheric processes. Failure to do so may result in undesired results, namely the extraction of features that are not fully independent, and a suboptimal dimensionality reduction of the problem.

### 5.3 Dynamic Source Analysis of the Geopotential Height Field

The theory introduced in sections 2 and 3 is now applied to the Geopotential Height Field datasets characterised in section 5.1. A set of spatiotemporal dynamic sources is extracted from the coevolutionary observables (the spatiotemporal geophysical datasets characterised in section 5.1), following section 2.4 and under the physical consistency and nonlinear disambiguation constraints in section 2.5.

In practice, the dynamical interactions are evaluated for , where is the order of the Sobolev space defined in section 2. This corresponds to the maximum order of weak differentiability of the observable fields, which in the illustrative case considered in the present section yields .

In the implementation of the DSA procedure, time derivatives up to fifth order are estimated from the spatiotemporal fields with [Geiser 2007]’s high-order time discretisation scheme based on the application of the Richardson extrapolation [Descombes 2001] to the second-order Crank-Nicolson method [Stoer and Burlish 2002]. In this regard, the time step is the temporal resolution of the spatiotemporal dataset and the field values for each involved cell in the scheme correspond to those from grid point data. An analogous reasoning is devised for spatial derivatives quantifying cross-variable sensitivities: in this case, the discrete spatial derivatives are computed using high-order compact ADI (alternating direction implicit) schemes [Karaa 2006].

While these schemes are originally intended to evaluate discrete derivatives approximating continuous ones (e.g. in computational fluid dynamics), they are also useful for evaluating discrete differentials from datasets, in which case the grid point neighbourhoods provide the scheme inputs natively. The main reason for the aforementioned approach rather than basic discrete derivative schemes stem from the inadequacy of the latter to reliably estimate high-order derivatives from data.

Without a priori restriction in the number of sources, an objective cut-off criterium is established whereby only dynamic sources sharing effective interaction information with the observables are retained. This leads us to two effective dynamic sources, henceforth denoted Meridional Dynamic Source (MDS) and Zonal Dynamic Source (ZDS) due to associated properties unveiled below.

Figure 4 depicts the spatial structure of the Meridional (Figure 4a) and the Zonal (Figure 4b) Dynamic Sources retrieved from the Geopotential Height datasets, i.e. from the observable through Dynamic Source Analysis (section 2). The spatial structure is computed from (3.2) with on , by taking into account all climatology partitions within the period 1951-2012, i.e. of all sub periods obtained by differential partition of the original period.

The spatial structure of the first Dynamic Source shows an improved geospatial signature of the North Atlantic Oscillation, as compared to spatial patterns obtained by PCA (e.g. [Barnston and Livezey 1987]). In fact, PCA-retrieved geospatial patterns in the literature are uncorrelated but not truly independent, due to the non-normality of the Geopotential Height fields as explained in section 5.2. Such PCA-retrieved patterns are actually different combinations of a smaller set of fundamentally independent processes, which leads to an excessive number of Principal Components and Empirical Orthogonal Functions, when a smaller set of fundamental features suffices to fully characterise the datasets.

As for the second Dynamic Source of the Geopotential Heights at 500 hPa, its spatial structure suggests the presence of an oscillation pattern with an East-West (zonal) dipole and centres of action over Greenland and the Baltic. Therefore, this is essentially a Zonal Dynamic Source representing a zonal dipole, whose statistical physics are consistent with the Baltic-Greenland Oscillation (BGO) unveiled by [Perdigão 2004].

Overall, the depiction of both retrieved spatial structures is consistent with the geospatial patterns extracted by [Perdigão 2004] with ICA, which had been performed by information-theoretical minimisation of mutual information and maximisation of information negentropy. Relative to ICA, DSA has the added value that not only the statistical climatological aggregate is independent, but also the dynamics within. That way, dynamic information is retained that can be used not only in downscaling but also in dynamic prediction of evolving precipitation distributions (section 6).

A common misconception associates statistical patterns with physical regimes. By statistically lumping information, the dynamical sequences upon which regimes should be built are lost. Different dynamical processes can lead to the exact same statistical distribution. That way, while the dynamics explain observed statistics, statistics do not elicit underlying dynamics.

In DSA the spatial structure of the underlying dynamics expresses state-spatial energy optima of the Geopotential Height Fields, i.e. its non-permanent, meta-stable states. It should be noted that no perennial stable states exist in the system, as the Earth is permanently far from equilibrium.

The temporal structure of the Dynamic Sources is depicted in Figure 5, by taking the joint temporal structure of all regional partitions within the study area. The Meridional Source is depicted in (a) and the Zonal in (b). These have been obtained by the temporal subspace retrieval of the respective Dynamic Sources, using the retrieval operation introduced in equation (4.1). For visualisation purposes alone, the depicted temporal structures are annually aggregated, which introduces artificial correlations. For that reason, their depiction may appear correlated when in reality the raw dynamic sources on a monthly basis are not.

The values of the temporal structures represent the standardised departure of the respective dynamic sources from their climatological centre of mass. In practice, positive [negative] values in the meridional and zonal cases correspond respectively to northward [southward] and eastward [westward] shifts in the atmospheric wave train and on its associated synoptic systems.

Together with the spatial structure, the temporal structure provides the dynamical evolution of the full geopotential height fields, i.e. their spatiotemporal structure. In mathematical terms, the spatiotemporal structure is the result of the interacting subspace composition (Appendix C) between the temporal and spatial structures.

In the present illustrative case, the manifold relating spatial and temporal structures of the dynamic sources underlying the geopotential height field at 500 hPa is given by the projection of the full spatiotemporal sources onto the coevolution manifold: . From within the multi-order spatiotemporal dependencies, in this particular case the spatiotemporal codependence is only manifested at orders . By integrating in time over the climatological period and evaluating the spatiotemporal interactions of up to fourth order, a spatial pattern emerges practically coinciding with the Negentropy map in Figure 3. Its visual depiction is thus redundant.

By noting that negentropy captures higher-order statistics, it might not be unreasonable to conjecture that beneath those higher-order statistics lie nonlinear spatiotemporal interactions. These in turn might be related to slow-fast climate dynamic interactions bearing in mind the role of spatial structures in representing processes deemed ”slow” relative to the time scale at which the temporal structures are evaluated [Perdigão and Blöschl 2014]. In particular, these signatures might suggest a nonlinear interplay between ”slow” multidecadal and ”fast” monthly dynamics in the climate system, with the former conditioning the regimes of the latter, and the latter feeding back on the structural properties of the former. A detailed investigation of the dynamical links beneath Negentropy beyond these conjectures makes for a fertile discussion and analysis that is outside of the scope of this paper, rather being deferred to a subsequent study.

From a fundamental physical standpoint, we interpret meridional dynamic sources in atmospheric dynamics as being associated to the latitudinal (meridional) curvature of the earth with associated differential surface heating by the sun. Meridional thermal gradients between the equator and the poles then result in poleward heat redistribution, which ultimately results in the meridional components of atmospheric flow. If the earth had cylindrical symmetry around its north-south axis, no such meridional dynamics would take place.

As far as the zonal dynamic sources are concerned, our interpretation consists on the fundamental mechanism underneath being the earth’s rotation. This is what fundamentally introduces the longitudinal (zonal) component in the atmospheric dynamics unleashed by the aforementioned meridional gradient. Zonal thermal gradients and associated baroclinicity between large ocean and land masses, while important, come as modulators of the large-scale zonal dynamics induced by planetary rotation.

Therefore, while observable zonal and meridional flow mutually interfere in the geophysical continuum, the fundamental zonal and meridional processes underneath the dynamics are in essence independent: the meridional gradients come from differential heating due to the shape of the earth, and the zonal from planetary rotation and longitudinal heterogeneities in the energy budget.

## 6 Dynamic Predictability and Simulation of Precipitation from the Geophysical Controls

The quest for geophysical controls on hydro-climatic processes aims at providing a better understanding of fundamental mechanisms and their relevance to the ultimate behaviour of the hydroclimatic system. Having retrieved and discussed such controls, the aim of the present section is to relate them with precipitation regimes: firstly by diagnosing dynamic relations between them, and secondly by dynamically predicting precipitation dynamics from the geophysical controls.

### 6.1 Diagnosing Dynamic Predictability

We begin by analysing measures of dynamic dependence between Precipitation and the Dynamic Sources of the Geophysical Controls for all climatology partitions within the study period.

Formally, the spatial structure of the dynamic interaction of order of the dynamic process (e.g. Precipitation) relative to the sources (e.g. MDS, ZDS) over a climatology can be expressed by taking equation (18) with given by :

(0) |

where is as in equation (3.2) with subspace , denotes, in this application, the joint spatial structure of all climatology partitions within the 1951-2012 study period.

In practice, a normalised dynamic codependence is computed from equation (6.1) as follows:

(0) |

Equation (6.1) yields essentially an adimensional measure of -order predictability among the intervening variables at each point on a map: in particular, Linear Predictability for and added value of Nonlinear Predictability for . The absolute value of the measure ranges from zero in dynamically independent cases to one in fully redundant cases at the respective order.

An integration of this measure over a sub manifold of the dynamical system yields an overall statistical predictability consistent with the Information-Theoretical Correlation introduced by [Pires and Perdigão 2007] as a generalisation of the Equivalent Correlation from [Perdigão 2004]. In particular for the linear case , the integration of the measure yields an aggregate diagnostic corresponding to a traditional correlation ranging from -1 to 1.

In practice, we are interested in effective interactions, i.e. the measured dynamic interactions discounting the spurious interactions between a given number NS of independently Monte Carlo (MC) shuffled dynamic sources and precipitation.

For that purpose, let denote the ensemble mean of the dynamical interactions between NS randomly shuffled realisations of and , with independent shuffling applied to each variable. Then the effective interaction of order is given by

(0) |

In the practical implementation considered in the present section, NS=.

The precipitation data has been retrieved from the Global Precipitation Climatology Centre (NCAR, 2014), and contains global analysis datasets of monthly precipitation on the earth’s land surface based on in situ rain gauge observational data. The data processing procedures conducted by GPCC are detailed in e.g. [Rudolf et al. 1994], [Rudolf and Schneider 2005] and involve statistics exclusively over the observational data.

### 6.2 Dynamic Predictability Results and Discussion

Figures 6 and 7 depict, respectively, the spatial structures of the Linear Predictability () and the added value from Nonlinear Predictability ( for ) of the Monthly Precipitation to the Meridional (a) and Zonal (b) Dynamic Sources, depicted in Figure 4, considering the joint spatial structure of all climatology partitions within the 1951-2012 period. The computations have been performed from equations (23)-(25) with representing the Dynamic Sources and the Monthly Precipitation. A truncation order was established in practice since higher-order terms no longer added significant value to the predictability.

The linear predictability is quite considerable over the centres of action of the underlying atmospheric oscillations, namely the Meridional Dynamic Source (MDS) in (a) and the Zonal Dynamic Source (ZDS) in (b). However, it degrades with the distance to those centres of action. As far as the MDS is concerned, this is consistent with the findings of statistical studies conducting linear sensitivity analysis of monthly precipitation to related atmospheric teleconnection patterns (e.g. [Murphy and Washington 2001] focusing on the Irish and British Isles and [Trigo et al. 2004] on Iberia), and also with the first-order components of the nonlinear information-theoretical analysis conducted by [Pires and Perdigão 2007] at the Euro-Atlantic scale.

Regarding the ZDS, no such studies are known to exist since the oscillation per se is a new finding, providing the dynamical basis underneath the information-theoretical pattern of the BGO unveiled by [Perdigão 2004]. In fact, while the literature features teleconnection patterns with a somewhat related zonal component (e.g. Eurasia-2 pattern in [Barnston and Livezey 1987] known as the Eastern Atlantic / Western Russia teleconnection), the Eurasia patterns are neither statistically independent nor dynamically based, rather relying on hierarchies of explained variance in climatological anomalies of geopotential height fields.

The ”lost” predictability away from main centres of action seen in Figure 6 is essentially ”recovered” when the nonlinear contributions to predictability are taken into account (Figure 7). In fact, the added value is highest at the linear-poorest regions, essentially closing the gap in predictability between the regions with higher linear predictability and those with scarcer one. As far as the MDS is concerned, this is consistent with the findings from [Pires and Perdigão 2007] regarding the nonlinear sensitivity of precipitation to the NAO, to which the statistical physics of MDS are closely related. Here, as there, regions where the linear sensitivity is weak (e.g. Central Europe) actually exhibit significant predictability once the information associated to nonlinear relations between precipitation and NAO is taken into account. In the present work this result is thus revisited with the newly extracted Meridional Dynamic Source, now with a nonlinear dynamical perspective on the NAO complementing the previous nonlinear statistical, information-theoretical one, and further extended to the Zonal Dynamic Source capturing the dynamics underlying the statistical physics of BGO from [Perdigão 2004].

All in all, the influence over precipitation is mostly linear over the centres of action and nonlinear away from them. By taking into account the full nonlinear information, then the overall domain exhibits strong predictability of precipitation from the dynamic sources.

While the linear information refers to the dynamics associated to mean atmospheric flows, the nonlinear stems from higher-order interactions such as turbulent fluxes. These are particularly important along storm tracks and regions dominated by convective behaviour. By bringing out the nonlinear predictability, DSA elicits not only the traditional large-scale synoptic regimes, but also hotspots of convectivity elusive to the linear terms.

The physical reason why smaller convective and turbulent processes can be dynamically inferred from large-scale atmospheric processes resides in the cooperative nature of synergistic interactions formulated in DSA. In fact, primary processes can synergistically cooperate to produce secondary processes at different scales (e.g. large-scale planetary waves resonate into finer-scale secondary perturbations in the atmosphere) via wave resonance across scales, which can be rigorously characterised in both nonlinear dynamic and information-theoretical terms [Pires and Perdigão 2015].

Our dynamical approach complements overall statistics by enabling the extraction of structural information within the domain, namely in its sub-partitions along with its dynamical evolution, i.e. a dynamic rather than lumped approach. By not simply averaging spatiotemporal fields in space and time upon separation, the coevolution manifold is preserved. Therefore, the spatiotemporal interplay is retained, enabling space and time to communicate and accounting for structural changes in space, with the added predictability that ensues under non-invariant non-ergodic climate dynamics.

Moreover, the consideration of all partitions within the spatiotemporal domain of study captures not only the broader Euro-Atlantic (spatial) and 60-year (temporal) scales, but also a whole spectrum within, down to spatiotemporal resolution limits (time step, grid size). Therefore, each pixel in the spatial structures contains not a temporal mean over a single period, but rather an integrated account of dynamic processes at multiple scales.

### 6.3 Dynamic Simulation of Monthly Precipitation from the Geophysical Controls

The dynamical simulation model [equation (4.2)], introduced in section 4.2 as a general simulation framework, is now applied to the simulation of Monthly Precipitation from the Dynamic Sources of the Geophysical Controls disentangled in section 5.3 with the methodologies devised in section 2. Taking and in equation (4.2) bring us to a practical model form:

(0) |

where prescribes the simulated dynamics of monthly precipitation, and is given by equation (2.4) with being the set of dynamic sources (meridional and zonal) of the geopotential height fields at 500 hPa determined in the previous section, and the model is truncated to the order in .

The model addresses the dynamic evolution of precipitation as a function of the dynamic sources representing the dominant modes of spatiotemporal dynamical variability in the geopotential height fields. Moreover, given the reciprocal nature of the Dynamic Interactions, it enables retroaction of Precipitation as it responds to the Dynamic Sources. This is physically important since the weather systems work towards depletion of the baroclinicity at their genesis, a negative feedback stemming from the Second Law of Thermodynamics. While the fundamental basis functionals are independent, their impact on as formulated in is indeed adjusted (internal feedback) as the dynamics unfold in (6.3).

The dynamic model is initialised in a two-step procedure: 1) evaluating the long-term dynamical properties of the reference manifold (), and 2) the short-term spatiotemporal divergence of the phase space flow from . The evaluation in 1) is done via the spatial structure of the dynamic interaction between and the atmospheric controls represented by the dynamic sources . In this regard, the spatial structure plays the role of the spatial legacy of long-term dynamics as in [Perdigão and Blöschl 2014]. The procedure in 2) finds the phase space directions that maximise the equienergetic entropy production, by concentrating the Lyapunov spectrum [Ott 2002] along the highest positive values whilst conserving the spectral integral and thus the energy of the system in that initialisation step. This procedure yields similar effects to the dynamical breeding performed in [Perdigão 2010] in that it projects the dynamics onto the most unstable directions, thus maximising the entropy production of the system in line with the physical principles. However, it differs from [Perdigão 2010] in that the operation is performed directly over the Lyapunov spectrum of the datasets.

Essentially, the initialisation leads to an initial distribution maximising the ability to further produce entropy in a realistic manner. In fact, while the reference manifold in 1) provides the structure of the slow dynamics (the climate or ”personality” of the system), the spatiotemporal divergence in 2) provides the dominant lines of deviation from that core, i.e. the seeds for the faster weather deviations from the climate norm (”mood swings”). The dynamics then unfold as prescribed by (6.3) and (2.4).

### 6.4 Dynamic Simulation Results and Discussion

The dynamic simulation of monthly precipitation from the underlying dynamic sources has been performed for two illustrative regions with predominantly Atlantic and Continental climatic norms: the former roughly encompasses the Loire Region in France, and the latter the Upper Danube Region in Central Europe. In geodesic coordinate terms, the simulations have been performed over the domains and , respectively. The denomination of these climate regions should not be mistaken for those of underlying hydrologic basins.

Overall, the observed Monthly Precipitation () is captured by the distribution produced by the dynamic simulation model (). This is seen in Figures 8 and 9 for the Loire and Upper Danube regions, respectively, where the observations (orange) are found within the modelled distribution (shades of blue), with very few exceptions. In general, the highest values of sit at the highest quantiles of , consistent with those values being upper extremes in the climatology. The converse happens for lower values of : in fact, these sit at the lower tail of the for the corresponding months. These considerations are further supported by Figure 10, showing in which quantile of the observation sits. The shape of the cumulative density function (cdf) is due in part to the fifth-order truncation of the model.

During the drier seasons (Summer in Loire, Winter in Upper Danube), the observed precipitation is captured by lower quantiles of the simulated precipitation, i.e. below their median. Since the outcome of the simulations is not the median but rather the full distribution for each month, and since the distribution does capture the observations at different quantiles, this behaviour does not constitute a modelling bias.

Whether the observations are captured by upper or lower quantiles of the simulated distributions actually comes down to physical arguments. For instance, higher-order turbulent processes associated with storms are captured by upper quantiles of the simulated precipitation. As for the aforementioned dry season behaviour, the observed precipitation naturally sits in lower quantiles of the simulated distributions due to the relative scarcity of precipitable water for precipitation to occur at the levels that would otherwise be expected given favourable low pressure conditions.

While moisture content in the atmosphere is not explicitly characterised, its influence is nonlinearly embedded in the geophysical fields. In fact, since moist air has lower density than drier air, it introduces lower pressure anomalies relative to those of drier air conditions. Therefore, air moisture implicitly contributes to the predictability of precipitation from the geophysical dynamics in the atmosphere - a contribution that would be elusive if only the linear information contained in the geopotential height fields would be taken into account.

By simulating an overall distribution, the model captures observations arising from a diversity of situations with associated quantiles as noted above. In doing so, the model captures not only precipitation outcomes associated to first-order atmospheric dynamics (mean atmospheric flows, large-scale synoptic processes) but also regional extremes arising from higher-order processes. An ensemble of simulated means (the usual paradigm in ensemble prediction) would not aptly capture the behaviour of extremes, as the physics governing the overall distributions entails laws differing from those governing the first moment. This is why our simulations take full distributions rather than ensembles of simulated means.

## 7 Concluding Remarks

Dynamic Source Analysis (DSA) was introduced as a novel nonlinear dynamic analysis and modelling framework for improving the fundamental physical understanding and predictability of complex coevolutionary systems. These were formulated in terms of the smallest set of generative fundamental non-redundant processes, i.e. a functional basis to the dynamical systems. In other words, DSA provided a model building and optimisation approach wherein a canonic structure is brought out for expressing the full system complexity in terms of the minimum set of components without loss of information. This enables data compression beyond the traditional feature extraction and information theoretical methods, by disentangling nonlinearly independent components not only over statistical aggregates but also within each step of the dynamics.

Given the natural mathematical ambiguity in solutions from nonlinear problems, the uniqueness of the DSA basis is ensured with first principles in Physics, namely on energy and entropy production. For that purpose, the physical solutions are dynamically selected from within the mathematically outcomes by automatically taking the energy minima in the solutions’ isentropic subspaces and entropy maxima in their equienergetic subspaces.

The generative ability of the DSA-retrieved functional basis stems from the formulation of synergistic interactions expressing cooperative dynamics among fundamental processes. In doing so, new features are brought as secondary processes that do not exist in any of the intervening primary processes. Further synergistic interactions among primary and secondary processes lead to a third generation of processes and beyond, in a synergistic cascade ultimately leading to the full representation of the system complexity.

While coevolutionary interactions entail mutual influence and thus redundancy, purely synergistic ones entail only innovative outcomes. A model structure grounded on the latter interactions is thus more efficient in representing complexity in the simplest way without loss of generality.

In practice, while all coevolutionary redundancies are cleared out between processes, self-interactions are preserved within each process (e.g. entropy production stemming from internal positive feedbacks). This irreducible coevolutionary behaviour is natural, as a process is redundant with itself. The coevolutionary nature of the problem is then concentrated along the ”diagonal spine” of the dynamical system, i.e. the diagonal of its dynamic interaction tensor, corresponding to self-interactions within each basis component. By taking the spectral view on such components, such self-interactions account for internal scale interactions or spectral coevolution, ranging from separable slow-fast dynamics to non-separable interactions among multiple indistinguishable scales (broadband coevolution).

In a spatiotemporal setting wherein the spatial structure embodies the legacy of long-term dynamics, spectral coevolution then comes down to a fundamental dynamic dependence between space and time. In this work, a coevolution manifold was introduced representing this codependence across the overall dynamical system, enabling space-for-time substitution to be performed at once for non-ergodic complex systems exhibiting diverse coevolution rates across different directions in space and time. This opens new perspectives in multivariate space-for-time trading, e.g. in regional hydrology.

In the hydro-climatic context, DSA was implemented as a means to retrieve a non-redundant set of nonlinearly interacting processes exerting control over precipitation in space and time. For that purpose, DSA has been implemented on the geopotential height fields at 500 hPa, retrieving a nonlinear spatiotemporal functional basis comprising two independent dynamic sources: the Meridional Dynamic Source (MDS) associated to the North Atlantic Oscillation (NAO), and the Zonal Dynamic Source (ZDS) embodying a pressure dipole between the Baltic and Greenland and providing a dynamic basis to the Baltic-Greenland Oscillation (BGO) unveiled by [Perdigão 2004]. The latter generalises and dynamically links the blocking systems found over the respective centres of action, namely those over Greenland and over the Baltic-Scandinavian area.

Whilst the NAO is responsible for the North-South shifting of the storm tracks over the Atlantic (with consequences downstream of the atmospheric flow well into Europe), the BGO does so in the West-East shift, counterbalancing or enhancing the zonal flow and thus being responsible for increasing or reducing the blocking risk in synoptic atmospheric processes. In fact, the BGO is associated to the zonal maritime-continental pressure gradients and associated atmospheric circulation, being affected by differential changes in temperature over land vs. over sea (which affects the zonal land-sea baroclinicity and thus the zonal circulation and progression of the synoptic systems).

The newly obtained dynamical sources were decomposed in space and time, by taking into account the spatiotemporal codependence as expressed by relative celerities in the atmospheric dynamics. For that purpose, a new spatiotemporal decomposition methodology of nonlinearly connected manifolds has been introduced and implemented, resulting in the spatial and temporal structures of the dynamical sources.

These have then been used to evaluate linear and nonlinear dynamical interactions and predictability between precipitation and the geopotential heights. The added value of the nonlinear interactions has been put into evidence, showing that by accounting for these the predictability of precipitation from the dynamical sources can be strongly improved throughout the domain.

On the temporal domain, the dynamic simulation of evolving Precipitation distributions has aptly captured the observational values for the test regions, as the latter were mostly captured by the former. Naturally though, higher[lower] observed values mostly lie on higher[lower] quantiles of the simulated distributions for each monthly record, with upper and lower observational extremes being at the corresponding tails of the simulated distributions.

All in all, Dynamic Source Analysis elicits fundamental processes and interactions at play in complex systems, eliminating coevolutionary redundancies and expressing complexity in the simplest dynamic architecture without loss of generality: a synergistic dynamical system. Moreover, by providing solutions constrained on first principles in physics, DSA ensures physical consistency in the resulting model structure and intervening processes.

With these advances in mind, a diversity of applications can thus be considered for the proposed theory of Dynamic Source Analysis, such as signal processing, feature extraction and inverse dynamic modelling from large datasets (big data), complex network design, dynamic model building and architecture optimisation of dynamical systems.

Further applications to hydro-climate dynamics and analytics are currently in progress, with particular emphasis on improving the dynamical understanding, model design and predictability of coevolutionary regimes, transitions and extremes in hydro-climatic systems.

## A Coevolution vs. Synergies in Information-Theoretical Terms

The fundamental difference between synergy and coevolution is that while the former involves cooperation between independent source components in the production of mutually dependent child processes, the latter entails dependence among the intervening components. In information terms, synergy produces extra information exceeding the sum of the information content of the parts (constructive interference pattern), whereas coevolution entails redundancy, resulting in the total information content being less than the sum of the informations of each individual components.

In order to illustrate these concepts, consider a pair of independent spatiotemporal dynamic sources and a set of observable child processes dynamically forced by the sources, corresponding to a complex dynamical system. The triadic interaction information among them is derived from [Pires and Perdigão 2015]’s equation (9d), yielding:

(0) |

where is Mutual Information (always non-negative) and is Interaction Information [Tsujishita 1995].

While and are independent from each other, they exhibit joint and individual codependence with the child processes . Equation (A) quantifies the difference between their joint contribution and the sum of their individual contributions to the coevolutionary system represented by . If the former exceeds the latter, and we are in the presence of a net synergy. The independence between and ensures that will at least be non-negative. This can be further seen by expressing the triadic interaction information as the difference between conditional and non-conditional information terms:

(0) |

In fact, when the sources are independent, and thus .

Naturally though, if the sources were codependent and their mutual information exceeded the conditional term in the first rhs (right-hand-side) term of (A), then the converse could happen, i.e. expressing net redundancy.

When looking at the overall statistic-dynamic properties of a complex system, coevolutionary and synergistic interactions can be captured in information-theoretical terms and linked to underlying dynamical systems, as done in [Perdigão and Blöschl 2014] and [Pires and Perdigão 2015]. In our proposed theory of Dynamic Source Analysis, a dynamic treatment is sought wherein coevolution and synergies are expressed in dynamic interaction terms. These aim to capture not only the overall nonlinear information as in the illustrative example (A)-(A), but also the underlying dynamic relationships that ultimately explain the diagnosed statistics.

## B Interacting Subspace Retrieval and Decomposition

### b.1 Retrieval Product

Consider the generic -dimensional functional or rank- tensor and the -dimensional functional or rank- tensor , along with a generic codependence -dimensional functional representing their interdependence.

We hereby define the Retrieval Product of -structure from as:

(0) |

where refers to the orthogonal complement.

As example applications, can be a manifold in a functional space or matrix in a vector space, and a sub manifold or sub matrix of .

This retrieval product is clearly not commutative, i.e. the -structure in is not necessarily the same as the -structure in . The exception occurs only when .

### b.2 Decomposition Operator

Let be a generic -dimensional functional or rank- tensor spanning an -dimensional space , and let be generic -dimensional subspace constituents of , with every , and every .

The subspace constituents are not necessarily independent from each other. Instead, they can interact through binding functionals of dimension , where is the set of values quantifying the dimension of the subspace constituents taken into consideration.

We hereby define an Interacting Subspace Decomposition Operator (.) as:

(0) |

where is the functional basis of the subspace spanned by , and is the retrieval product defined in equation (B.1).

The application of to a full partition of denoted as is then given as the set of all within . The partitions are generic and not necessarily disjoint.

## C Interacting Subspace Composition

### c.1 Composition Product

Consider processes represented by generic -dimensional functionals , , with . These processes interact in a -dimensional functional , with . The Composition Product among the functionals is defined as:

(0) |

with

(0) |

Let denote the dimensionality or rank of the lhs (left-hand-side) of (C.1). Then, , i.e. the composition product is rank additive, albeit discounting the dimensionality of the interactions.

### c.2 Composition Operator

Consider now a partition into interacting sub manifolds as defined in the previous Appendix (section B). Our aim is to generate the full manifold from its subspace constituents. This can be done by introducing a Interacting Subspace Composition Operator (or Composition Operator for short) that essentially implements the Composition Product among the partition members (sub manifolds) of :

(0) |

with denoting the partition member and denoting the cardinality of the partition (i.e. the number of its members).

The interacting subspace composition of a spatiotemporal structure generalises the traditional invariant space-time manifolds in Mathematical Physics by taking into account the spatiotemporal codependence quantified by the dynamic codependence functional corresponding to the spatiotemporal coevolution manifold. That is, we operate in a generalised setting where space and time are not necessarily independent dimensions in a structurally invariant manifold, but actually coevolving via nonlinear dynamical interactions [Perdigão 2016].

##### Acknowledgements:

This research was supported by the ERC Advanced Grant ’Flood Change’ project no. 291152. C.P. also acknowledges support from the Portuguese Foundation for Science and Technology (FCT) through the project RECI/GEO-MET/0380/2012 - SHARE - Seamless High-resolution Atmosphere-ocean Research, and through the FCT funding contract UID/GEO /50019/2013 for Instituto Dom Luiz. Our families are gratefully acknowledged for their omnipresent support.

## References

- Aires et al. 1999 Aires F., Chédin A., Nadal J.-P. (1999): Analyse de séries temporelles géophysiques et théorie de lÕinformation: lÕanalyse en composantes indépendantes. C. R. Acad. Sci. Paris, 328, 569-575.
- Aires et al. 2002 Aires W., Rossow B., Chédin A. (2002): Rotation of EOFs by the independent component analysis: Towards a solution of the mixing problem in the decomposition of geophysical time series. J. Atmos. Sci., 59, 111-123.
- Bárdossy and Pegram 2011 Bárdossy A., Pegram G. (2011): Downscaling precipitation using regional climate models and circulation patterns toward hydrology. Water Resources Research, 47(4), W04505.
- Barnston and Livezey 1987 Barnston A.G., Livezey R.E. (1987): Classification, seasonality and persistence of low-frequency atmospheric circulation patterns. Mon. Weather Rev., 115, 1083-1126.
- Basak et al. 2004 Basak J., Sudarshan A., Trivedi D., Santhanam M.S. (2004): Weather Data Mining Using Independent Component Analysis. Journal of Machine Learning Research, 5, 239-253.
- Boers et al. 2013 Boers N., Bookhagen B., Marwan N., Kurths J., Marengo J. (2013): Complex networks identify spatial patterns of extreme rainfall events of the South American Monsoon System. Geophys. Res. Lett., 40, 4386-4392, doi:10.1002/grl.50681, 2013.
- Comon 1994 Comon, P. (1994): Independent component analysis, a new concept? Signal Process., 36, 287-314.
- Courant and Hilbert 1953 Courant R., Hilbert D. (1953): Methods of Mathematical Physics, Vol. I, Interscience.
- Cover and Thomas 1991 Cover T.M., Thomas J.A. (1991): Elements of Information Theory. Wiley, 576 pp.
- Da Silva 2001 Da Silva A.C. (2001): Lectures on symplectic geometry (Vol. 1764). Springer Science & Business Media.
- Descombes 2001 Descombes St. (2001): Convergence of a splitting method of high order for reaction-diffusion systems. Mathematics of Computations, 70, 1481-1501, 2001.
- Donges et al. 2009 Donges J.F., Zou Y., Marwan N., Kurths J. (2009): Complex networks in climate dynamics: Comparing linear and nonlinear network construction methods. The European Physical Journal Special Topics, Vol. 174, Issue 1, pp 157-179.
- Ehrlich and Raven 1964 Ehrlich P.R., Raven P.H. (1964): Butterflies and plants: A study in coevolution. Evolution, 18, 586-608.
- Fodor and Kamath 2003 Fodor I.K., C. Kamath (2003): Using independent component analysis to separate signals in climate data, Proceedings of the Independent Component Analyses, Wavelets, and Neural Networks, SPIE Proceedings, 5102, 25-36. SPIE Press, Bellingham, WA.
- Geiser 2007 Geiser J. (2007): Higher-order difference and higher-order splitting methods for 2D parabolic problems with mixed derivatives. International Mathematical Forum, Vol. 2. No. 67. 2007.
- Grines and Pochinka 2010 Grines V., Pochinka O. (2010): Energy functions for dynamical systems. Regular and Chaotic Dynamics, 15.2, 185.
- Haken 1983 Haken H. (1983): Synergetics, an Introduction: Nonequilibrium Phase Transitions and Self-Organization in Physics, Chemistry, and Biology, 3rd rev. enl. ed. New York: Springer-Verlag.
- Hasselmann 1988 Hasselmann K. (1988): PIPs and POPs: The reduction of complex dynamical systems using principal interaction and oscillation patterns. J. Geophys. Res., 93(D9), 11015-11021, doi:10.1029/JD093iD09p11015.
- Horel 1981 Horel J.D. (1981): A rotated principal component analysis of the inter annual variability of the Northern Hemisphere 500 mb height field. Mon. Wea. Rev., 109, 2080-1092.
- Hyvärinen et al. 2001 Hyvärinen A., Karhunen J., Oja E. (2001): Independent Component Analysis. John Wiley and Sons, 478 pp.
- Hyvärinen and Pajunen 1999 Hyvärinen A., Pajunen P. (1999): Nonlinear Independent Component Analysis: Existence and Uniqueness results. Neural Networks 12(3): 429-439.
- Hocke and Kämpfer 2008 Hocke K., Kämpfer N. (2008): Bispectral analysis of the long-term recording of surface pressure at Jakarta, J. Geophys. Res., 113, D10113, doi:10.1029/2007JD009356.
- Karaa 2006 Karaa S. (2006): High-Order Compact ADI Methods for Parabolic Equations. Journal of Computers and Mathematics with Applications, Vol. 52, Iss. 8-9, 1343-1356, 2006.
- Kalnay et al. 1996 Kalnay E., M. Kanamitsu, R. Kistler, W. Collins, D. Deaven, L. Gandin, M. Iredell, S. Saha, G. White, J. Woollen, Y. Zhu, A. Leetmaa, R. Reynolds, M. Chelliah, W. Ebisuzaki, W. Higgins, J. Janowiak, K. C. Mo, C. Ropelewski, J. Wang, R. Jenne, and D. Joseph (1996): The NCEP/NCAR 40-Year Reanalysis Project. Bull. Amer. Meteor. Soc., 77, 437-471.
- Karl et al. 1982 Karl T.R., Koscielny A.J., Diaz H.F. (1982): Potential errors ion the application of principal component (eigenvector) analysis to geophysical data. J. Appl. Meteor., 21, 1183-1186.
- Kidson and Thompson 1998 Kidson J.W., Thompson C. (1998): A comparison of statistical and model-based downscaling techniques for estimating local climate variations. Journal of Climate, 11, 735-753.
- Komen et al. 1996 Komen G.J., Cavaleri L., Donelan M., Hasselmann K., Hasselmann S., Janssen P.A.E.M. (1996): Dynamics and Modelling of Ocean Waves. Cambridge University Press, 339 pp.
- Lage 1995 Lage E.J.S. (1995): Física Estatística. Fundação Calouste Gulbenkian, 650 pp., ISBN 9723106469.
- Mandelbrot 1974 Mandelbrot B. (1974): Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier. Journal of Fluid Mechanics, 62, 331-358.
- Maskey et al. 2015 Maskey M.L., Puente C.E., Sivakumar B., Cortis, A. (2015): Encoding daily rainfall records via adaptations of the fractal multifractal method. Stochastic Environmental Research and Risk Assessment, 1-15.
- Holton 2004 Holton J.R. (2004): An Introduction to Dynamic Meteorology. 4th ed., Academic Press, 535 pp.
- Murphy and Washington 2001 Murphy S.J., Washington R. (2001): United Kingdom and Ireland precipitation variability and the North Atlantic sea-level pressure field. International Journal of Climatology, 21.8 (2001): 939-959.
- Nicolis and Nicolis 2007 Nicolis G., Nicolis C. (2007): Foundations of Complex Systems. World Scientific.
- Nicolis et al. 2009 Nicolis C., Perdigão R.A.P., Vannitsem S. (2009): Dynamics of prediction errors under the combined effect of initial condition and model errors. Journal of the Atmospheric Sciences, 66, 766-778.
- Obregón et al. 2002 Obregón N., Puente C.E., Sivakumar B. (2002): Modeling high resolution rain rates via a deterministic fractal-multifractal approach. Fractals, 10(3):387-394.
- Oja 1997 Oja E. (1997): The nonlinear PCA learning rule in independent component analysis. Neurocomputing, 17(1):25-46, 1997.
- Ott 2002 Ott E. (2002): Chaos in Dynamical Systems. 2nd ed., Cambridge University Press.
- Pajunen et al. 1996 Pajunen P., Hyvärinen A., Karhunen J. (1996): Non-Linear Blind Source Separation by Self-Organizing Maps. Proc. Int. Conf. on Neural Information Processing, Hong Kong, pp. 1207-1210, 1996.
- Paschalis et al. 2013 Paschalis A., Molnar P., Fatichi S., Burlando P. (2013): A stochastic model for high-resolution space-time precipitation simulation. Water Resour. Res., 49, 8400-8417, doi:10.1002/2013WR014437.
- Peixoto and Oort 1992 Peixoto J.P., Oort A.H. (1992): The Physics of Climate. Springer-Verlag, New York, 1992.
- Perdigão 2016 Perdigão R.A.P. (2016): Lecture Notes in Fluid Dynamical Systems, TU Wien course 222.562, 2016.
- Perdigão 2010 Perdigão R.A.P. (2010): Nonlinear Statistics and Dynamics of Atmospheric Predictability and Downscaling. Doctoral Thesis in Physics, University of Lisbon, Portugal. http://hdl.handle.net/10451/2013.
- Perdigão 2004 Perdigão R.A.P. (2004): Independent Component Analysis of the Geopotential Height Field and its Relevance to Precipitation Regimes over the Euro-Atlantic Area. Thesis in Geophysical Sciences, University of Lisbon, Portugal.
- Perdigão and Blöschl 2014 Perdigão R.A.P., Blöschl G. (2014): Spatiotemporal flood sensitivity to annual precipitation: Evidence for landscape-climate coevolution. Water Resour. Res., 50, 5492-5509, doi:10.1002/2014WR015365.
- Pires and Perdigão 2007 Pires C.A.L., Perdigão R.A.P. (2007): Non-Gaussianity and asymmetry of the winter monthly precipitation estimation from the NAO. Monthly Weather Review, 135, 430-448.
- Pires and Perdigão 2012 Pires C.A.L., Perdigão R.A.P. (2012): Minimum Mutual Information and Non-Gaussianity Through the Maximum Entropy Method: Estimation from finite samples. Entropy 2012, 14(6), 1103-1126, doi: 10.3390/e14061103
- Pires and Perdigão 2013 Pires C.A.L., Perdigão R.A.P. (2013): Minimum Mutual Information and Non-Gaussianity Through the Maximum Entropy Method: Estimation from finite samples. Entropy 2013, 15(3), 721-752, doi:10.3390/ e15030721.
- Pires and Perdigão 2015 Pires C.A.L., Perdigão R.A.P. (2015): Non-Gaussian interaction information: estimation, optimization and diagnostic application of triadic wave resonance. Nonlinear Processes in Geophysics, 22, 87-108, doi:10.5194/npg-22-87-2015.
- Pires and Ribeiro 2016 Pires C.A. Ribeiro A. (2016): Separation of the atmospheric variability into non-Gaussian multidimensional sources by projection pursuit techniques. Climate Dynamics, DOI: 10.1007/s00382-016-3112-9.
- Preisendorfer et al. 1988 Preisendorfer R.W., Mobley C.D., Barnett T.P. (1988): The principal discriminant method of prediction: Theory and evaluation. J. Geophys. Res., 93(D9), 10815-10830, doi:10.1029/JD093iD09p10815.
- Puente 1992 Puente C.E. (1992): Multinomial multifractals, fractal interpolators, and the Gaussian distribution. Phys. Lett. A, 161, 441-447.
- Reed and Simon 1980 Reed M., Simon B. (1980): Methods of modern mathematical physics: Functional analysis (Vol. 1). Gulf Professional Publishing.
- Roman 2005 Roman S. (2005): Advanced linear algebra (Vol. 3). Heidelberg: Springer.
- Ruddell and Kumar 2009 Ruddell B.L., Kumar P. (2009): Ecohydrologic process networks: 1. Identification. Water Resour. Res., 45, W03419, doi:10.1029/2008WR007279.
- Rudolf et al. 1994 Rudolf B., Hauschild H., Rueth W., Schneider U. (1994): Terrestrial Precipitation Analysis: Operational Method and Required Density of Point Measurements. NATO ASI I/26, Global Precipitations and Climate Change (Ed. M. Desbois and F. Desalmand), Springer Verlag Berlin, 173 - 186.
- Rudolf and Schneider 2005 Rudolf B., Schneider U. (2005): Calculation of gridded precipitation data for the global land-surface using in-situ gauge observations. Proceedings of the 2nd Workshop of the International Precipitation Working Group IPWG, Monterey, October 2004, 231-247.
- Salby 2012 Salby M.L. (2012): Physics of the Atmosphere and Climate. Cambridge University Press. 666 pp.
- Shannon 1948 Shannon C.E. (1948): The mathematical theory of communication. Bell Systems Technical Journal, 27, 379-423.
- Sobolev 1938 Sobolev S.L. (1938): Sur un théoréme de l’analyse fonctionelle. Mat. Sb., 46 (4) (1938), 471-496.
- Sobolev 1963 Sobolev S.L. (1963): Applications of Functional Analysis in Mathematical Physics. Amer. Math. Soc. Transl. Math. Mono., Providence, RI.
- Song et al. 2014 Song Y., Li Y., Bates B., Wikle C.K. (2014): A Bayesian hierarchical downscaling model for south-west Western Australia rainfall. Journal of the Royal Statistical Society: Series C (Applied Statistics), 63: 715-736. doi: 10.1111/rssc.12055.
- Spivak 2005 Spivak M. (2005): A comprehensive introduction to differential geometry, vol. I, 3 edition. Publish or Perish, Boston, Mass., 491 pp, ISBN 0914098705.
- Stoer and Burlish 2002 Stoer J., Burlisch R.(2002): Introduction to Numerical Analysis. 3rd. ed., Springer-Verlag, New-York, 2002.
- Trigo et al. 2004 Trigo R.M., Pozo-Vázquez D., Osborn T. J., Castro-Díez Y., Gámiz-Fortis S., Esteban-Parra M. J. (2004): North Atlantic oscillation influence on precipitation, river flow and water resources in the Iberian Peninsula. Int. J. Climatol., 24: 925-944. doi:10.1002/joc.1048.
- Troch et al. 2015 Troch P.A., Lahmers T., Meira A., Mukherjee R., Pedersen J.W., Roy T., Valdes-Pineda R. (2015): Catchment coevolution: A useful framework for improving predictions of hydrological change? Water Resour. Res., 51, doi:10.1002/2015WR017032.
- Tsujishita 1995 Tsujishita T.: On triple mutual information, Adv. Appl. Math., 16, 269-274, 1995.
- Tu et al. 2014 Tu J.H., Rowley C.W., Luchtenburg D.M., Brunton S.L., Kutz J.N. (2014): On dynamic mode decomposition: Theory and applications. Journal of Computational Dynamics, 1, 2, 391.
- Wallace and Gutzler 1981 Wallace J.M., Gutzler D.S. (1981): Teleconnections in the geopotential height field during the Northern Hemisphere winter. Mon. Wea. Rev., 109, 784-812.
- Whipple et al. 1999 Whipple K.X., Kirby E., Brocklehurst S.H. (1999): Geomorphic limits to climate-induced increases in topographic relief. Nature, 401(6748), 39-43.
- Wilks and Wilby 1999 Wilks D.S., Wilby R.L. (1999): The weather generation game: a review of stochastic weather models. Progress in Physical Geography, 23, 3, 329-357, doi: 10.1177/030913339902300302.
- Yu et al. 2014 Yu X., Hu D., and Xu J. (2014): Blind source separation: theory and applications. Wiley, New York, 416 pp. ISBN: 978-1-118-67984-5.