# Climate Predictions: The Chaos and Complexity in Climate Models

Authors: D. T. Mihailović; G. Mimić; I. Arsenić
Journal: Advances in Meteorology
 (2014)

---

## Abstract

Some issues which are relevant for the recent state in climate modeling have been considered. A detailed overview of literature related to this subject is given. The concept in modeling of climate, as a complex system, seen through Gödel’s theorem and Rosen’s definition of complexity and predictability is discussed. Occurrence of chaos in computing the environmental interface temperature from the energy balance equation given in a difference form is pointed out. A coupled system of equations, often used in climate models, was analyzed. It is shown that the Lyapunov exponent mostly has positive values allowing presence of chaos in this system. The horizontal energy exchange between environmental interfaces, which is described by the dynamics of driven coupled oscillators, was analyzed. Their behavior and synchronization, when a perturbation is introduced in the system, as a function of the coupling parameter, the logistic parameter, and the parameter of exchange, were studied calculating the Lyapunov exponent under simulations with the closed contour ofN
=
100 environmental interfaces. Finally, we have explored possible differences in complexities of two global and two regional climate models using their air temperature and precipitation output time series. The complexities were obtained with the algorithm for calculating the Kolmogorov complexity.


---

## Body

## 1. Introduction

Among the most interesting and fascinating phenomena that are predicted/predictable is the chaotic ocean/atmosphere/land system called weather and its longtime average, climate. While weather is not predictable beyond a few days, aspects of the climate may be predictable for years, decades, and perhaps longer [1]. These two phrases clearly summarise the current opinion and state in climate modeling community that deals with the aforementioned subjects. However, the question of the weather and climate modeling and predictability has been initiated in the early sixties of the 20th century, which was elaborated in pioneering works by Lorenz [2–5]. He was the first person in the scientific world who explicitly pointed out the following points related to the nonlinear dynamics in atmospheric motion: (i) question of prediction and predictability, (ii) importance of understanding the nonlinearity in modeling procedure, (iii) demand for discovery of chaos, and (iv) careful consideration of sensitivity of differential equations in modeling system on initial conditions. Subsequent three decades after appearance of these papers have been characterised by strong interest for predictability of weather and climate on theoretical and practical level. The following topics have been set in the focus: (1) dynamics of error growth; (2) linear and nonlinear systems (normal modes, optimal modes, nonlinear geophysical systems, and scale selection in error growth); (3) predictability of systems with many scales; (4) limit of predictability; (5) weather predictability (growth of errors in General Circulation Models (GCMs) based on Lorenz’s analysis); (6) predictability from analogs (targeted observations); (7) climate predictability (predictability of time-mean quantities, predictability of the second kind) and potential predictability; (8) seasonal mean predictability; and (9) El Niño-Southern Oscillation (ENSO) chaos, predictability of coupled models, and decadal modulation of predictability [6–18]. Because the focus of our paper is complexity and predictability in climate modeling, we finish this overview with the comment by Orell (2003) as follows: “Prediction problems have been described by Lorenz as falling into two categories. Problems that depend on the initial condition, such as short- to medium-range weather forecasting, are described as predictions of the first kind, while problems that depend on boundary rather than initial conditions, such as, in many cases, the longer-term climatology, are referred to as predictions of the second kind. Both kinds of prediction will be affected by error in the model equations used to approximate the true system” [19–21].Earth’s atmosphere has evolved into a complex system in which life and climate are intricately interwoven. The interface between Earth and atmosphere as a “pulsating biophysical organism” is a complex system itself. We use the term complex system in Rosen’s sense (Rosen, 1991) as it was explicated in the comment by Collier (2003) as follows: “In Rosen’s sense a complex system cannot be decomposed non-trivially into a set of part [sic] for which it is the logical sum. Rosen’s modeling relation requires this. Other notions of modeling would allow complete models of Rosen style complex systems, but the models would have to be what Rosen callsanalytic, that is, they would have to be a logical product. Autonomous systems must be complex. Other types of systems may be complex, and some may go in and out of complex phases” [22, 23]. Also, we will explain in which sense the termcomplexity will be used in further text. Usually, that is an ambiguous term, sometimes used [22] to refer to systems that cannot be modeled precisely in all respects. However, following Arshinov and Fuchs (2003) the term “complexity” has three levels of meaning [24]. (1) There is self-organization and emergence in complex systems [25]. (2) Complex systems are not organized centrally but in a distributed manner; there are many connections between the system’s parts [25, 26]. (3) It is difficult to model complex systems and to predict their behaviour even if one knows to a large extent the parts of such systems and the connections between the parts [25, 27]. The complexity of a system depends on the number of its elements and connections between the elements (the system’s structure). According to this assumption, Kauffman (1993) defines complexity as the “number of conflicting constraints” in a system [26]; Heylighen (1996) says that complexity can be characterized by a lack of symmetry (symmetry breaking) which means that “no part or aspect of a complex entity can provide sufficient information to actually or statistically predict the properties of the others parts;” [28] Edmonds (1996) defines complexity as “that property of a language expression which makes it difficult to formulate its overall behavior, even when given almost complete information about its atomic components and their inter-relations” [29]. Aspects of complexity are things, people, number of elements, number of relations, nonlinearity, broken symmetry, nonholonic constraints, hierarchy, and emergence [30].Generally, predictability refers to the degree that a correct forecast of a system’s state can be made either qualitatively or quantitatively. For example, while the second law of thermodynamics can tell us about the equilibrium that a system will evolve to and steady states in dissipative systems can sometimes be predicted, there exists no general rule to predict the time evolution of systems far from equilibrium, that is, chaotic systems, if they do not approach some kind of equilibrium. Their predictability usually deteriorates with time. To quantify predictability, the rate of divergence of system trajectories in phase space can be measured (Kolmogorov-Sinai entropy, Lyapunov exponents).Lorenz (1984) discussed several issues in the predictability of weather systems [31]. According to him predictability is defined as the degree of accuracy with which it is possible to predict the state of weather system in the near and also the distant future (predictability in Lorenz’s sense). In this paper it is assumed that weather predictions are made on the basis of imperfect knowledge of weather system’s present and past states. This rather general statement is comprehensively elaborated by Hunt (1999) [20]. He described the fundamental assumptions and current methodologies of the two main kinds of environmental forecast (i.e., weather forecast); the first is valid for a limited period of time into the future and over a limited space-time “target” and is largely determined by the initial and preceding state of the environment, such as the weather or pollution levels, up to the time when the forecast is issued and by its state at the edges of the region being considered; the second kind provides statistical information over long periods of time and/or over large space-time targets so that they only depend on the statistical averages of the initial and “edge” conditions. Environmental forecasts depend on the various ways that models are constructed. These range from those based on the “reductionist” methodology (i.e., the combination of separate, scientifically based models for the relevant processes) to those based on statistical methodologies, using a mixture of data and scientifically based empirical modeling. For example, limitations of the predictability in the world of atmospheric motions are concisely discussed in paper by James (2002) [32]. In this paper the predictability of a forced nonlinear system is numerically considered, proposed by Lorenz, as a compelling heuristic model of the midlatitude global circulation.The above insight of the predictability is underlined in the context of the “environmental predictability” (primarily linked to the climate change issues); we finish with the following question:can we significantly “improve” the weather/climate predictions compared to the level they currently reached? The answer cannot be strictly elaborated with eitheryes orno. An optimistic and acceptable attitude, that prefers optionyes, is concisely written down by Hunt (1999) as the following phrase: “We concluded that philosophical studies of how scientific models develop and of the concept of determinism in science are helpful in considering these complex issues” [20]. If we give advantage to the optionno then we do not close the door for the first option. It only means that there exists limitation of the modeling attempts on an epistemological level. To show that, we will use Gödel’s incompleteness theorem about Number Theory [33]. Basically it says that no matter how one tries to formalize a particular part of mathematics, syntactic truth in the formalization does not coincide with the set of truths about numbers. In other words Gödel’s theorem shows that formalizations are part of mathematics but not all of mathematics. There are many ways to look and “read” Gödel’s theorem. One exclusive way is offered by Rosen (1985) [34]. According to him the first thing to bear in mind is that both Number Theory and any formalization of it are systems of entailment. It is therelation between them, or more specifically the extent to which these schemes of entailment can be brought into congruence, that is of primary interest. The establishment of such congruencies, through the positing of referents in one of them for elements of the other, is the essence of themodeling relation. In a precise sense, this theorem asserts that a formalization in which all entailment is syntactic entailment is tooimpoverished in entailmentto be congruent to Number Theory, no matter how we try to establish such congruence. This kind of situation is termedcomplexity by Rosen (1977) [35]. Namely, in this light, Gödel’s theorem says that Number Theory is morecomplex than any of its formalizations or, equivalently, that formalizations, governed by syntactic inference alone, aresimpler than Number Theory. To reach Number Theory from its formalizations or, more generally, to reach a complex system from simpler one requires some kind of limiting processes.Let us return to the question we were asking ourselves after we had shortly considered climate modeling (i.e., predictability) beyond the complexity. To our mind there is a significant space for “improvement” of models and their capabilities to provide good forecasts. It can be done only if the modeling attempts are directed towards the following steps: from structures and states to processes and functions; from self-correcting to self-organizing systems; from hierarchical steering to participation; from conditions of equilibrium to dynamic balances of nonequilibrium; from single trajectories to bundles of trajectories; from linear causality to circular causality; from predictability to relative chance; from order and stability to instability, chaos, and dynamics; from certainty and determination to a larger degree of risk, ambiguity, and uncertainty; from reductionism to emergentism; from being to becoming.In this paper we address two issues that, to our mind, are important for further improvements in designing the climate models. (1) The phenomenon of chaos: (i) in behavior of the environmental interface temperature computed from the energy balance equation, (ii) in coupling of processes of vertical and horizontal energy transfers in climate models which can result in something that is much more complex than the deterministic chaos of those models (Section2), and (2) complexity analysis of the climate model output time series which is elaborated in Section 3.


## 2. Chaos in Modeling the Global Climate System

### 2.1. The Current Issues in Modeling the Global Climate System

The target of global climate models is the Earth’s climate system, consisting of the physical and chemical components of the atmosphere, ocean, land surface, and cryosphere. In climate simulations, the objective is to correctly simulate the spatial variation of climate conditions in some average sense. There exists a hierarchy of different climate models, ranging from simple energy balance models to the very complex global circulation models. These models attempt to account for as many processes as possible to simulate the detailed evolution of the atmosphere, ocean, cryosphere, and land system, at a horizontal resolution that is typically 100 s of km. Climate model complexity is result of the nonlinearity of the equations, high dimensionality, and the linking of multiple subsystems. However, the secret of understanding of climate model complexity lies in the nonlinear dynamics of the atmosphere and oceans, which is described by the Navier-Stokes equations whose solution is one of the most vexing problems in all of mathematics.We shortly enhance the main current issues related to the modeling of the global climate system [36]. (i)Chaos. Weather can be considered as being in state of deterministic chaos, owing to its sensitivity to initial conditions. The source of the chaos is nonlinearities in the Navier-Stokes equations. Therefore, a consequence of sensitivity to initial conditions is that beyond a certain time the system will no longer be predictable more than seven days. Climate models are also sensitive to initial conditions. However, in addition, in these models, coupling of a nonlinear, chaotic atmospheric model to a nonlinear, chaotic ocean model gives rise to something much more complex than the deterministic chaos of the weather model. Those coupled modes give rise to bifurcation, instability, and chaos. The situation is further complicated because the coupled atmosphere/ocean system cannot be classified by current theories of nonlinear dynamical systems, where definitions of chaos and attractor cannot be invoked in situations involving transient changes of parameter values [36–38]. (ii)Confidence in Climate Models.The relevant issue is how well the climate model reproduces reality, that is, whether the model “works” and is fit for its intended purpose. In the absence of model verification or falsification, Stainforth et al. (2007) describe the challenges of building confidence in predictions using current models and consider the implications for experimental design and the balance of resources in climate modeling research. (iii) We are aware that our understanding of, and ability to simulate, the Earth’s climate is rather limited. That fact causes theclimate model imperfection[39], which is divided into two types: uncertainty and inadequacy. The term model uncertainty means that we cannot reliably choose parameter values (or ensembles of parameter values), which will provide the most informative results. In addition, further complications arise from the choice of parameterization. Finally, model uncertainty is associated with uncertainty in model parameters, subgrid parameterizations, and also initial conditions. It is a well-known problem that numerical models of natural systems cannot be identical to the structure of those systems; that is, they cannot be isomorphic to the real system [39]. In other words, they are inadequate; that is, before we run any simulation of the future, we know in advance that models are unrealistic representations of many relevant aspects of the real-world system [39–41]. And, finally, atmospheric science has played a leading role in the development and use of computer simulation in scientific endeavors. Climate simulations of future states of weather and climate have important societal applications. Thus, we should have in mind this following statement by Heymann (2010): “Computer simulation in the atmospheric sciences has caused a host of epistemic problems, which scientists acknowledge and with which philosophers and historians are grappling with [sic]. But historically practice overruled the problems of epistemology. Atmospheric scientists found and created their proper audiences, which furnished them with legitimacy and authority. Whatever these scientists do, it does not only tell us something about science, it tells us something about the politics and culture within which they thrive…” [42]. In the subsections that follow we have considered the issue related to the chaos.In Section2.2, first we illustrate a possibility of the occurrence of chaos in computing the environmental interface temperature from the energy balance equation. This temperature is used in calculating global annual mean temperature. This is a variable which has been most comprehensively studied, susceptible to relatively simple analysis and in which we might arguably achieve the greatest confidence in prediction [43].In Section2.3, we have considered the dynamics of the horizontal and vertical transfer of energy, for example, advection, convection, and diffusion between environmental interfaces, which is described by the dynamics of driven coupled oscillators and formalism of the category theory. Coupling of processes of vertical and horizontal energy transfers in climate models gives rise to something much more complex than the deterministic chaos of those models.


### 2.2. Energy Balance Equation and Chaotic Behavior of Environmental Interface Temperature

#### 2.2.1. Background

Traditional mathematical analysis of physical systems tacitly assumes that integers and all real numbers, no matter how large or how small, are physically possible and all mathematically possible trajectories are physically possible [44]. Traditionally, this approach has worked well in physics and in engineering, but it does not lead to a very good understanding of chaotic systems, which, as is now known, are extremely important in the study of real-world phenomena ranging from weather to biological systems. In this paper we deal with one issue in modeling pathways in meteorology as well as in physics, biology, and chemistry, that is, environmental sciences in their broadest context [45], in particular in autonomous dynamical systems, which are common subject under consideration in climate modeling. Namely, we consider how to replace given differential equations by appropriate difference equations in environmental modeling and thus in climate simulations [46].According to van der Vaart many models for environmental problems have been and will be built in the form of differential equations or systems of such equations [46]. With the advent of computers one has been able to find (approximate) solutions for equations that used to be intractable. Many of the mathematical techniques have been used in this area to replace given differential equations by appropriate difference equations. So a huge effort has been invested into choice of appropriate difference equations whose solutions are “good” approximations to the solutions of the given differential equations. This question includes a requirement for better understanding of the fundamental problem: interrelations between classical continuum mathematics and reality in different sciences. For many atmospheric phenomena the “continuum” type of thinking, that is, at the basis of any differential equation, is not natural to the phenomenon but rather constitutes an approximation to a basically discrete situation: in much work of this type the “infinitesimal step lengths” handled in the reasoning which lead us to the differential equation are not really thought of as infinitesimally small but as finite; yet, in the last stage of such reasoning, where the differential equation rises from the differentials, these “infinitesimal” step lengths go to zero, that is, where the above-mentioned approximation comes in. Under this kind of circumstances, it seems more naturalto build the model as a discrete difference equation from the start, without going through the painful, doubly approximative process of first, during the modeling stage, finding a differential equation to approximate a basically discrete situation and then, for numerical computing purposes, approximating that differential equation by a difference scheme [44].In this subsection we analyze the energy balance equation in procedure of computing the environmental interface temperature and the deeper soil layer temperature commonly used in climate models. The environmental interface is defined asinterface between two biotic or abiotic environments that are in relative motion and exchange energy, matter, and information through physical, biological, and chemical processes, fluctuating temporally and spatially regardless of space and time scale[47].There are a lot of examples of environmental interfaces in the nature, but here we deal with the ground surface, where there exist all three mechanisms of energy transfer: incoming and outgoing radiation, convection of heat and moisture into the atmosphere, and conduction of heat into deeper soil layers of ground (Figure1) [48]. Parameterization of these processes is of great importance for environmental models of different spatial and temporal scales and thus climate ones. In the paper by Mihailović and Mimić (2012) it is shown that ground surface is treated as a complex system in which chaotic fluctuations occur while we compute its temperature [49]. This system, as an actual dynamic system, is very sensitive to initial conditions and arbitrarily small perturbation of the current trajectory that may lead to its unpredictable behavior. In the aforementioned paper the lower boundary condition, that is, the deeper soil layer temperature, was constant, but it can also vary in time. That system, often used in environmental models, is of interest to be analyzed by the methods of nonlinear dynamics. Having in mind those facts, in this, we (i) perform a nonlinear dynamical analysis of coupled system for computing the environmental interface temperature and the deeper soil layer temperature and (ii) examine behavior of the coupled system in dependence on the main system parameters, in order to show the possible occurrence of the chaos in computing the environmental interface temperature. Firstly, we consider difference form of the energy balance equation and deeper soil layer temperature equation transforming them into the coupled system with the corresponding parameters, then we analyze behavior of the solutions of the coupled system, and we have examined domains of stability using the Lyapunov exponent.Figure 1
Energy balance equation terms.


#### 2.2.2. Physical Background and Derivation of the Coupled System

One of the most important conditions for functioning of any complex system is a proper supply of the system with energy. Dynamics of energy flow is based on the energy balance equation [48]. As we mentioned before, environmental interface is a complex system. General difference form of energy balance equation for the ground surface as an environmental interface is
(1)
C
g
Δ
T
g
Δ
t
=
R
net
-
H
-
λ
E
-
G
,
where T
g is the ground surface temperature, Δ
t is the time step, C
g is the soil heat capacity, R
net is the net radiation, H is the sensible heat flux, λ
E is the latent heat flux, and G is the heat flux into the ground. First, we assume that the net radiation is given as in [50]; that is,
(2)
R
net
=
C
R
(
T
g
-
T
a
)
,
where T
a is the air temperature at some reference level and C
R is the coefficient for the net radiation term. Second, we make expansion of the exponential term in the expression for latent heat flux as follows:
(3)
λ
E
=
C
L
d
[
b
(
T
g
-
T
a
)
+
b
2
2
(
T
g
-
T
a
)
2
]
,
where C
L is the water vapour transfer coefficient, b
=
0.06337
°
C
-
1, and d is parameter which occurs in expanding the series [51]. Further, the conduction of the heat into the soil can be written in the following form:
(4)
G
=
C
D
(
T
g
-
T
d
)
,
where C
D is the coefficient of the heat conduction while T
d is the temperature of the deeper soil layer. The sensible heat flux H can be parameterized as
(5)
H
=
C
H
(
T
g
-
T
a
)
,
where C
H is the sensible heat transfer coefficient. The prognostic equation for temperature of the deeper soil layer T
d is
(6)
Δ
T
d
Δ
t
=
1
τ
(
T
g
-
T
d
)
,
where τ
=
86400 s. After collecting all terms (2)–(6), the coupled system takes the following form:
(7)
C
g
Δ
T
g
Δ
t
=
C
R
(
T
g
-
T
a
)
-
C
H
(
T
g
-
T
a
)
-
C
L
d
[
b
(
T
g
-
T
a
)
+
b
2
2
(
T
g
-
T
a
)
2
]
-
C
D
(
T
g
-
T
d
)
,
Δ
T
d
Δ
t
=
1
τ
(
T
g
-
T
d
)
.
More details about the nature and the range of physical parameters C
R, C
L, C
D, and C
H can be found in [52]. Now, using the time scheme forward in time (n indicates the time step) and dividing both sides of (7) with the constant temperature T
0 (e.g., value of mean Earth temperature, i.e., T
0
=
288 K) we get
(8)
T
g
n
+
1
-
T
a
T
0
=
T
g
n
-
T
a
T
0
+
Δ
t
C
g
C
R
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
H
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
L
b
d
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
L
d
T
0
b
2
2
(
T
g
n
-
T
a
)
2
T
0
2
-
Δ
t
C
g
C
D
T
g
n
-
T
a
T
0
+
Δ
t
C
g
C
D
T
d
n
-
T
a
T
0
,
T
d
n
+
1
-
T
a
T
0
=
T
d
n
-
T
a
T
0
+
Δ
t
τ
T
g
n
-
T
a
T
0
-
Δ
t
τ
T
d
n
-
T
a
T
0
.
Finally, introducing replacements x
=
(
T
g
-
T
a
)
/
T
0 and y
=
(
T
d
-
T
a
)
/
T
0, where x is the dimensionless environmental interface temperature and y is the dimensionless deeper soil layer temperature, we reach the coupled system as follows:
(9)
x
n
+
1
=
A
x
n
-
B
x
n
2
+
C
y
n
,
y
n
+
1
=
D
x
n
+
(
1
-
D
)
y
n
,
where A
=
1
+
(
Δ
t
/
C
g
)
(
C
R
-
C
H
-
C
L
b
d
-
C
D
), B
=
C
L
d
T
0
(
b
2
Δ
t
/
2
C
g
), C
=
Δ
t
(
C
D
/
C
g
), and D
=
Δ
t
/
τ. Introducing the replacement x
1
,
n
=
A
x
n
/
B, where x
1 is the modified dimensionless environmental interface temperature and x
2
,
n
=
y
n, we can write the following:
(10)
x
1
,
n
+
1
=
A
x
1
,
n
(
1
-
x
1
,
n
)
+
C
B
A
x
2
,
n
,
x
2
,
n
+
1
=
D
A
B
x
1
,
n
+
(
1
-
D
)
x
2
,
n
.
Analysis of values of parameters A, B, C, and D, based on a large number of energy flux outputs from the land surface scheme runs, indicates that their values are ranged in the following intervals: (i) A
∈
[
0,4
] and (ii) B, C, and D are ranged in the interval [
0,1
]. Thus, A is the logistic parameter, which from now will be denoted by r. All other groups of parameters in the system (10) have the values in the same interval [
0,1
]. Let us underline that under some circumstances those parameters can be equal. Correspondingly, we replaced all of them by introducing the coupling parameter c.Finally, system (10) can be written in the form of coupled maps; that is,
(11)
x
1
,
n
+
1
=
r
x
1
,
n
(
1
-
x
1
,
n
)
+
c
x
2
,
n
,
x
2
,
n
+
1
=
c
(
x
1
,
n
+
x
2
,
n
)
.Now we analyse the stability of physical solutions of coupled maps (11), using the Lyapunov exponent, which is a measure of convergence or divergence of near trajectories in phase space. Sign of Lyapunov exponent is characteristic of attractor type and for stable fixed point it is negative, although for chaotic attractor it is positive. Calculating Lyapunov exponent for the coupled system (11) with values of parameters c
∈
(
0.05,0.1
) and r
∈
(
3.6,3.8
), because we thought that it will be interesting to investigate behavior of the system for small values of coupling parameter and high values of logistic parameter, we got results depicted in Figure 2 [53]. It is shown that the Lyapunov exponent mostly has positive values approving presence of chaos in this system, but there are still some strait regions where the Lyapunov exponent is negative and where the solutions of the coupled system are stable, that is, domains of stability.Figure 2
Lyapunov exponent of the coupled system (11), which shows presence of strait regions of stability in highly developed chaos.Irregularities in solution of the system of difference equation (11) can come from two reasons. They are (i)numerical, that is, because we try to choose appropriate difference equation whose solution is “good” approximation to the solution of the given partial differential equation and (ii)physical, that is, occurrence of chaotic fluctuations in the considered system because the environmental interface cannot oppose an enormous radiative forcing, suddenly reaching the interface. Therefore, it raises the question whether we can find either domain or domains where physically meaningful solutions exist [44]. We do that by considering the stability of physical solution of (11). Stability, in mathematics, is condition in which a slight disturbance in a system does not produce too disrupting effect on that system. In terms of the solution of a differential equation, a function f
(
ζ
) is said to be stable if any other solution of the equation that starts out sufficiently close to it when ζ
=
0 remains close to it for succeeding values of ζ. If the difference between the solutions approaches zero as ζ increases, the solution is called asymptotically stable. If a solution does not have either of these properties, it is called unstable. We consider the stability of physical solution of (11) in sense of the Lyapunov exponent (λ). Thus, the system is stable if λ
<
0 and unstable if λ
>
0. This analysis, set into context of the climate modeling, points out the fact that there exists set of domains where the environmental interface temperature cannot be calculated by the physics of currently designed climate models, because of occurrence of chaotic phenomena on the environmental interfaces.


### 2.3. Horizontal Energy Exchange between Environmental Interfaces

#### 2.3.1. Background

There are three major sets of processes that must be considered when constructing a climate model: (i) radiative (the transfer of radiation through the climate system, e.g., absorption and reflection); (ii) dynamic (the horizontal and vertical transfer of energy, e.g., advection, convection, and diffusion); and (iii) surface process (inclusion of processes involving land/ocean/ice and the effects of albedo, emissivity, and surface-atmosphere energy exchanges). If the nonlinearities in these processes are treated improperly, then while designing the model, the complexity and thus its reliability will not be retained in the highest degree. In Section2.2 we have considered surface-atmosphere energy exchanges with cadence on the phenomenon of a possible occurrence of the chaos in solving the energy balance equation for calculating the environmental interface temperature in climate models. Here, following Mihailović et al. (2012) we analyze the horizontal energy exchange between environmental interfaces which is described by the dynamics of driven coupled oscillators [54]. In order to study their behavior, when a perturbation is introduced in the system, as a function of the coupling parameter, the logistic parameter, and the horizontal energy exchange intensity (parameter of exchange, in further text), we considered dynamics of two maps serving the diffusive coupling [54].As noted above, the horizontal exchange of energy between environmental interfaces is considered as diffusion-like process. The dynamics of energy exchange behavior on environmental interface are typically expressed as a logistic mapΦ
(
x
)
=
r
x
(
1
-
x
), where x is the dimensionless temperature of environmental interface and r is a logistic parameter [53, 54]. However, we use an alternative form of this map, which includes a parameter p that represents the horizontal energy exchange intensity (Figure 3). By introducing this parameter we formalize an intrinsic property of the environmental interfaces, which depends on the nature of the interface. The environmental interface dynamics are expressed here as a difference equation, so we avoid the double approximation of (i) finding a differential equation to approximate an essentially discrete process (during the modeling stage) and then (ii) approximating that differential equation by a difference scheme for numerical computing purposes [44, 46]; that is,
(12)
Φ
(
x
i
,
n
)
=
r
x
i
,
n
p
(
1
-
x
i
,
n
p
)
.
The dynamics of this map (12) are governed by two parameters, p and r, which express intrinsic property of the environmental interfaces and the influence of the environment, respectively.Figure 3
Schematic diagram of horizontal energy exchange between two environmental interfaces. Parametersp and r express intrinsic property of the environmental interfaces and the influence of the environment, respectively.Since these and many other processes on environmental interface are defined as diffusion-like, we will explore (i) how these processes can be better represented in climate models by introducing parameter of exchangep in the diffusive coupling associated with the horizontal energy exchange and (ii) how the horizontal energy exchange intensity dynamics are affected by the perturbation of parameters that represents influence of the environment, environmental interface coupling, and horizontal energy exchange intensity.In considering these problems we have to include observational heterarchy, a challenging topic when dealing with complex systems. Essentially, observational heterarchy reveals that it is impossible to unambiguously determine to which subsystems an element belongs [55]. Therefore, the dynamics of the complex system are articulated in terms of two kinds of dynamics, Intent and Extent dynamics and the interaction between them, where Intent corresponds to an attribute of a given phenomenon and Extent corresponds to a collection of objects satisfying that phenomenon [55].


#### 2.3.2. Observational Heterarchy and Horizontal Energy Exchange between Environmental Interfaces

Observational heterarchy consists of two sets of intralayer maps, called Intent and Extent perspectives, and interlayer operations satisfying the following conditions. (1) The interlayer operations inherit the mixture of intra- and interlayer operations and (2) there is a procedure by which the interlayer operation can be regarded as an adjoint functor. If the interlayer operation satisfies the conditions (1) and (2), it is called a prefunctor [55]. Preserving the above-mentioned composition occurs as follows: a prefunctor, 〈
F
〉
: Intent → Extent is mapping a set, X, to a set, 〈
F
〉
X, and map Φ to a map, f
*
Φ
f, where f
*
f
(
x
)
=
x for all x in f
(
X
) with f
(
X
)
:
〈
F
〉
X
→
X. In this sense we call f
* pseudoinverse of  f. Because applying a prefunctor to a map is expressed as composition of maps, it satisfies the conditions (1) and (2). The approximation is defined by the assumption that f is a one-to-one onto map. If one accepts the approximation, f
*
=
f
-
1 holds, then a prefunctor can become a functor. Given two maps,Φ
,
Ψ
:
X
→
X,
(13)
〈
F
〉
(
Φ
)
〈
F
〉
(
Ψ
)
=
(
f
*
Φ
f
)
(
f
*
Ψ
f
)
=
f
*
Φ
(
f
f
*
)
Ψ
f
=
f
*
Φ
(
f
f
-
1
)
Ψ
f
=
f
*
Φ
Ψ
f
=
〈
F
〉
(
Φ
f
)
.
It implies that F preserves the composition of maps, Φ and Ψ.The time development of the environmental surface dynamicsx
i
,
n, for two interfaces, is expressed as
(14)
x
i
,
n
+
1
=
(
1
-
c
)
Φ
(
x
i
,
n
)
+
f
(
Φ
(
x
j
,
n
)
)
,
where n is the time iteration, i
,
j
=
1,2, x
i
,
n
∈
[
0,1
], c is the coupling parameter, f is the map representing the horizontal energy exchange between environmental interfaces, and Φ is one of maps in the pair (
Ψ
,
Φ
) whose composition is preserved by a prefunctor 〈
F
〉. Here, we apply the framework of an observational heterarchy to the two environmental interface systems. If Intent and Extent are denoted by Φ and Ψ, respectively, the time development of the concentration is expressed as x
i
,
n
+
1
=
(
1
-
c
)
Φ
(
x
i
,
n
)
+
Ψ
(
x
j
,
n
). In this expression, if Ψ
(
X
)
=
f
(
Φ
(
x
)
), then it can be reduced to (14).We perform our analysis following the procedure described in [55]. First, in this subsection we address the synchronization of the passive coupling for two environmental interfaces given by (14) and (12), and then, in the next subsection, we will show that perturbation can modify the dynamics and enhance robust behavior in a multienvironmental interface system of active coupling. Synchronization is well-known collective phenomenon in various multicomponent physical as well as the climate systems [56–58]. The exchange of information (coupling) among the components can be either global or local. Here, we consider that chaotic systems are synchronized only when the largest Lyapunov exponent of the driven system is negative. It was calculated by approach proposed in [59]. We studied the stability of the fixed point by linearizing n
≥
2 component coupled system and obtain Z
n
+
1
=
ζ
n
Z
n where ζ
n is the Jacobian of this system evaluated in (
0,0
,
…
,
0
) and ζ
n
=
(
x
1
,
n
,
x
2
,
n
,
…
,
x
N
,
n
). By iterating we obtain
(15)
Z
n
+
1
=
(
∏
s
=
0
n
ζ
s
)
Z
0
,
and thus we get Lyapunov exponent as follows:
(16)
λ
=
lim
⁡
n
→
∞
⁡
(
ln
⁡
∥
∏
s
=
0
n
ζ
s
∥
)
n
.
Figure 4 depicts the diagram of normalized frequency of synchronization F
p (λ
<
0) for system of two environmental interfaces passively coupled ((12) and (14)), as a function of parameter of exchange p, averaged over all values of the coupling parameter c and logistic parameter r. The value of the normalized frequency of synchronization F
p is calculated as
(17)
F
p
=
∑
N
n
(
λ
<
0
)
∑
N
n
(
λ
<
0
)
+
∑
N
p
(
λ
>
0
)
,
where N
n
(
λ
<
0
) and N
p
(
λ
>
0
) are numbers of negative and positive values of the Lyapunov exponent, respectively. These numbers were calculated for the fixed value of p, while c and r were changing in intervals (
0,1
) and (
1,4
), respectively, with the step of 0.05. From this figure it is seen that after p
>
0.2, F
p becomes lower, indicating a decrease of number of states, which are synchronized.Figure 4
Normalized frequency of synchronization,F
p
(
λ
<
0
), for system of environmental interfaces passively coupled (12)–(14) as a function of parameter of exchange p. An averaging was done over all values of coupling parameter c and logistic parameter r.


#### 2.3.3. Simulations of Active Coupling in a Multienvironmental Interface System

Here, we address the behavior of active coupling [55] and estimate whether a coupled map system described above can achieve synchronization under influence of perturbations. The dynamics of two environmental interface systems called active coupling [55], used for simulations, is expressed as(18a)
x
i
,
n
+
1
=
(
1
-
c
)
Φ
n
(
x
i
,
n
)
+
Ψ
n
(
x
j
,
n
)
,
(18b)
Ψ
n
+
1
=
f
Φ
n
f
*
f
,
(18c)
Φ
n
+
1
=
f
*
Ψ
n
+
1
,
(18d)
Φ
n
(
x
i
,
n
)
=
r
x
i
,
n
p
(
1
-
x
i
,
n
p
)
.We note that the dynamical system defined by (18a) and (18d) is called the passive coupling, and that is a usual coupled map system. The active coupling can be approximated to passive coupling, where the approximation is defined by the adjunction or the equivalence between Intent and Extent. Compared with passive coupling, the behavior of active coupling is much more complex [55]. In (18a)–(18d), because of a pseudoinverse map, f
*, all calculations are defined to be approximations. In simulations, the Intent map was a discontinuous map, expressed by Φ
n
+
1
=
f
*
Ψ
n
+
1.In order to see how perturbation enhances robust behavior in the framework of observational heterarchy in a multienvironmental interface system represented by closed contour of coupled environmental interfaces exchanging the energy horizontally, then the system of coupled difference equations forN environmental interfaces exchanging the energy can be written in the form of the following matrix equation:(19a)
X
N
1
=
(
A
+
B
)
·
X
N
.
The elements in matrices in (19a) are
(19b)
X
N
1
i
,
n
+
1
=
x
i
,
n
+
1
,
X
N
i
,
n
=
x
i
,
n
,
A
i
,
k
=
(
1
-
c
)
Φ
n
(
x
i
,
n
)
δ
i
,
k
,
(19c)
B
i
,
k
=
{
Ψ
n
(
x
k
,
n
)
,
k
=
i
+
1
,
i
<
N
0
,
k
≠
i
+
1
,
i
<
N
Ψ
n
(
x
k
,
n
)
,
k
=
1
,
i
=
N
0
,
k
≠
1
,
i
=
N
,where i
=
1,2
,
…
,
N and δ
i
,
k is the Kronecker symbol.Simulations with the active coupling, defined by (18a)–(18d), were performed with and without perturbation given as in [56]. The results of simulations are shown in Figure 5. In this figure Lyapunov exponent λ is plotted against coupling parameter c for active coupling with perturbation (black line) compared to the passive coupling (red line), for different values of the parameter of exchange pand the logistic parameter r. Simulations were performed with the closed contour of N
=
100 interfaces. The Lyapunov exponent was calculated using (15)-(16) and the Jacobian of the system given by (19b)-(19c) is representing this contour.Figure 5
Diagram of Lyapunov exponent,λ, against coupling parameter c for the fluctuated active coupling defined by (18a)–(18d) with P (black line) compared to passive coupling and N (red line) for different values of parameter of exchange p and logistic parameter r. In (a)–(c) p takes the fixed values (1.0, 0.5, 0.2), while r
=
4. In (d)–(f) p is randomly chosen, while r takes values 4.0, 3.82, and 3.6, respectively. Simulations were performed with the closed contour of N
=
100 environmental interfaces.In calculatingλ, for each c from 0.0 to 1.0 with step 0.005, 1
0
4 iterations were applied for an initial state, and then the first 1
0
3 steps were abandoned. In order to see how the active coupling modifies the synchronization of horizontal energy exchange between environmental interfaces, we performed two kinds of simulations. Firstly, we used r
=
4.0 and the fixed value of the parameter of exchange p (Figures 5(a)–5(c)); secondly, we used a randomly chosen p and a logistic r parameter with the values of 4.0, 3.82, and 3.6, respectively (Figures 5(d)–5(f)). Figures 5(a)–5(c) depict that, in the chaotic regime (r
=
4.0), regardless of the value p, the Lyapunov exponent is always positive (λ
>
0) and therefore the process of the horizontal energy exchange in a multienvironmental interface system is always unsynchronized. However, the stormy perturbation disturbs this state (Figures 5(a)–5(c)). Although the logistic parameter is settled at r
=
4.0 for chaotic behavior, the coupling parameter c tunes interaction and leads to synchronization in some intervals, particularly for p
=
1.0 and p
=
0.5. This behavior is more pronounced in Figures 5(d)–5(f) where p is randomly chosen; here the process of horizontal energy exchange in multienvironmental interfaces exhibits a strong tendency towards the synchronization, even though the logistic parameter r is in chaotic region (r
=
4.0,3.82 and 3.6).


## 2.1. The Current Issues in Modeling the Global Climate System

The target of global climate models is the Earth’s climate system, consisting of the physical and chemical components of the atmosphere, ocean, land surface, and cryosphere. In climate simulations, the objective is to correctly simulate the spatial variation of climate conditions in some average sense. There exists a hierarchy of different climate models, ranging from simple energy balance models to the very complex global circulation models. These models attempt to account for as many processes as possible to simulate the detailed evolution of the atmosphere, ocean, cryosphere, and land system, at a horizontal resolution that is typically 100 s of km. Climate model complexity is result of the nonlinearity of the equations, high dimensionality, and the linking of multiple subsystems. However, the secret of understanding of climate model complexity lies in the nonlinear dynamics of the atmosphere and oceans, which is described by the Navier-Stokes equations whose solution is one of the most vexing problems in all of mathematics.We shortly enhance the main current issues related to the modeling of the global climate system [36]. (i)Chaos. Weather can be considered as being in state of deterministic chaos, owing to its sensitivity to initial conditions. The source of the chaos is nonlinearities in the Navier-Stokes equations. Therefore, a consequence of sensitivity to initial conditions is that beyond a certain time the system will no longer be predictable more than seven days. Climate models are also sensitive to initial conditions. However, in addition, in these models, coupling of a nonlinear, chaotic atmospheric model to a nonlinear, chaotic ocean model gives rise to something much more complex than the deterministic chaos of the weather model. Those coupled modes give rise to bifurcation, instability, and chaos. The situation is further complicated because the coupled atmosphere/ocean system cannot be classified by current theories of nonlinear dynamical systems, where definitions of chaos and attractor cannot be invoked in situations involving transient changes of parameter values [36–38]. (ii)Confidence in Climate Models.The relevant issue is how well the climate model reproduces reality, that is, whether the model “works” and is fit for its intended purpose. In the absence of model verification or falsification, Stainforth et al. (2007) describe the challenges of building confidence in predictions using current models and consider the implications for experimental design and the balance of resources in climate modeling research. (iii) We are aware that our understanding of, and ability to simulate, the Earth’s climate is rather limited. That fact causes theclimate model imperfection[39], which is divided into two types: uncertainty and inadequacy. The term model uncertainty means that we cannot reliably choose parameter values (or ensembles of parameter values), which will provide the most informative results. In addition, further complications arise from the choice of parameterization. Finally, model uncertainty is associated with uncertainty in model parameters, subgrid parameterizations, and also initial conditions. It is a well-known problem that numerical models of natural systems cannot be identical to the structure of those systems; that is, they cannot be isomorphic to the real system [39]. In other words, they are inadequate; that is, before we run any simulation of the future, we know in advance that models are unrealistic representations of many relevant aspects of the real-world system [39–41]. And, finally, atmospheric science has played a leading role in the development and use of computer simulation in scientific endeavors. Climate simulations of future states of weather and climate have important societal applications. Thus, we should have in mind this following statement by Heymann (2010): “Computer simulation in the atmospheric sciences has caused a host of epistemic problems, which scientists acknowledge and with which philosophers and historians are grappling with [sic]. But historically practice overruled the problems of epistemology. Atmospheric scientists found and created their proper audiences, which furnished them with legitimacy and authority. Whatever these scientists do, it does not only tell us something about science, it tells us something about the politics and culture within which they thrive…” [42]. In the subsections that follow we have considered the issue related to the chaos.In Section2.2, first we illustrate a possibility of the occurrence of chaos in computing the environmental interface temperature from the energy balance equation. This temperature is used in calculating global annual mean temperature. This is a variable which has been most comprehensively studied, susceptible to relatively simple analysis and in which we might arguably achieve the greatest confidence in prediction [43].In Section2.3, we have considered the dynamics of the horizontal and vertical transfer of energy, for example, advection, convection, and diffusion between environmental interfaces, which is described by the dynamics of driven coupled oscillators and formalism of the category theory. Coupling of processes of vertical and horizontal energy transfers in climate models gives rise to something much more complex than the deterministic chaos of those models.


## 2.2. Energy Balance Equation and Chaotic Behavior of Environmental Interface Temperature

### 2.2.1. Background

Traditional mathematical analysis of physical systems tacitly assumes that integers and all real numbers, no matter how large or how small, are physically possible and all mathematically possible trajectories are physically possible [44]. Traditionally, this approach has worked well in physics and in engineering, but it does not lead to a very good understanding of chaotic systems, which, as is now known, are extremely important in the study of real-world phenomena ranging from weather to biological systems. In this paper we deal with one issue in modeling pathways in meteorology as well as in physics, biology, and chemistry, that is, environmental sciences in their broadest context [45], in particular in autonomous dynamical systems, which are common subject under consideration in climate modeling. Namely, we consider how to replace given differential equations by appropriate difference equations in environmental modeling and thus in climate simulations [46].According to van der Vaart many models for environmental problems have been and will be built in the form of differential equations or systems of such equations [46]. With the advent of computers one has been able to find (approximate) solutions for equations that used to be intractable. Many of the mathematical techniques have been used in this area to replace given differential equations by appropriate difference equations. So a huge effort has been invested into choice of appropriate difference equations whose solutions are “good” approximations to the solutions of the given differential equations. This question includes a requirement for better understanding of the fundamental problem: interrelations between classical continuum mathematics and reality in different sciences. For many atmospheric phenomena the “continuum” type of thinking, that is, at the basis of any differential equation, is not natural to the phenomenon but rather constitutes an approximation to a basically discrete situation: in much work of this type the “infinitesimal step lengths” handled in the reasoning which lead us to the differential equation are not really thought of as infinitesimally small but as finite; yet, in the last stage of such reasoning, where the differential equation rises from the differentials, these “infinitesimal” step lengths go to zero, that is, where the above-mentioned approximation comes in. Under this kind of circumstances, it seems more naturalto build the model as a discrete difference equation from the start, without going through the painful, doubly approximative process of first, during the modeling stage, finding a differential equation to approximate a basically discrete situation and then, for numerical computing purposes, approximating that differential equation by a difference scheme [44].In this subsection we analyze the energy balance equation in procedure of computing the environmental interface temperature and the deeper soil layer temperature commonly used in climate models. The environmental interface is defined asinterface between two biotic or abiotic environments that are in relative motion and exchange energy, matter, and information through physical, biological, and chemical processes, fluctuating temporally and spatially regardless of space and time scale[47].There are a lot of examples of environmental interfaces in the nature, but here we deal with the ground surface, where there exist all three mechanisms of energy transfer: incoming and outgoing radiation, convection of heat and moisture into the atmosphere, and conduction of heat into deeper soil layers of ground (Figure1) [48]. Parameterization of these processes is of great importance for environmental models of different spatial and temporal scales and thus climate ones. In the paper by Mihailović and Mimić (2012) it is shown that ground surface is treated as a complex system in which chaotic fluctuations occur while we compute its temperature [49]. This system, as an actual dynamic system, is very sensitive to initial conditions and arbitrarily small perturbation of the current trajectory that may lead to its unpredictable behavior. In the aforementioned paper the lower boundary condition, that is, the deeper soil layer temperature, was constant, but it can also vary in time. That system, often used in environmental models, is of interest to be analyzed by the methods of nonlinear dynamics. Having in mind those facts, in this, we (i) perform a nonlinear dynamical analysis of coupled system for computing the environmental interface temperature and the deeper soil layer temperature and (ii) examine behavior of the coupled system in dependence on the main system parameters, in order to show the possible occurrence of the chaos in computing the environmental interface temperature. Firstly, we consider difference form of the energy balance equation and deeper soil layer temperature equation transforming them into the coupled system with the corresponding parameters, then we analyze behavior of the solutions of the coupled system, and we have examined domains of stability using the Lyapunov exponent.Figure 1
Energy balance equation terms.


### 2.2.2. Physical Background and Derivation of the Coupled System

One of the most important conditions for functioning of any complex system is a proper supply of the system with energy. Dynamics of energy flow is based on the energy balance equation [48]. As we mentioned before, environmental interface is a complex system. General difference form of energy balance equation for the ground surface as an environmental interface is
(1)
C
g
Δ
T
g
Δ
t
=
R
net
-
H
-
λ
E
-
G
,
where T
g is the ground surface temperature, Δ
t is the time step, C
g is the soil heat capacity, R
net is the net radiation, H is the sensible heat flux, λ
E is the latent heat flux, and G is the heat flux into the ground. First, we assume that the net radiation is given as in [50]; that is,
(2)
R
net
=
C
R
(
T
g
-
T
a
)
,
where T
a is the air temperature at some reference level and C
R is the coefficient for the net radiation term. Second, we make expansion of the exponential term in the expression for latent heat flux as follows:
(3)
λ
E
=
C
L
d
[
b
(
T
g
-
T
a
)
+
b
2
2
(
T
g
-
T
a
)
2
]
,
where C
L is the water vapour transfer coefficient, b
=
0.06337
°
C
-
1, and d is parameter which occurs in expanding the series [51]. Further, the conduction of the heat into the soil can be written in the following form:
(4)
G
=
C
D
(
T
g
-
T
d
)
,
where C
D is the coefficient of the heat conduction while T
d is the temperature of the deeper soil layer. The sensible heat flux H can be parameterized as
(5)
H
=
C
H
(
T
g
-
T
a
)
,
where C
H is the sensible heat transfer coefficient. The prognostic equation for temperature of the deeper soil layer T
d is
(6)
Δ
T
d
Δ
t
=
1
τ
(
T
g
-
T
d
)
,
where τ
=
86400 s. After collecting all terms (2)–(6), the coupled system takes the following form:
(7)
C
g
Δ
T
g
Δ
t
=
C
R
(
T
g
-
T
a
)
-
C
H
(
T
g
-
T
a
)
-
C
L
d
[
b
(
T
g
-
T
a
)
+
b
2
2
(
T
g
-
T
a
)
2
]
-
C
D
(
T
g
-
T
d
)
,
Δ
T
d
Δ
t
=
1
τ
(
T
g
-
T
d
)
.
More details about the nature and the range of physical parameters C
R, C
L, C
D, and C
H can be found in [52]. Now, using the time scheme forward in time (n indicates the time step) and dividing both sides of (7) with the constant temperature T
0 (e.g., value of mean Earth temperature, i.e., T
0
=
288 K) we get
(8)
T
g
n
+
1
-
T
a
T
0
=
T
g
n
-
T
a
T
0
+
Δ
t
C
g
C
R
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
H
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
L
b
d
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
L
d
T
0
b
2
2
(
T
g
n
-
T
a
)
2
T
0
2
-
Δ
t
C
g
C
D
T
g
n
-
T
a
T
0
+
Δ
t
C
g
C
D
T
d
n
-
T
a
T
0
,
T
d
n
+
1
-
T
a
T
0
=
T
d
n
-
T
a
T
0
+
Δ
t
τ
T
g
n
-
T
a
T
0
-
Δ
t
τ
T
d
n
-
T
a
T
0
.
Finally, introducing replacements x
=
(
T
g
-
T
a
)
/
T
0 and y
=
(
T
d
-
T
a
)
/
T
0, where x is the dimensionless environmental interface temperature and y is the dimensionless deeper soil layer temperature, we reach the coupled system as follows:
(9)
x
n
+
1
=
A
x
n
-
B
x
n
2
+
C
y
n
,
y
n
+
1
=
D
x
n
+
(
1
-
D
)
y
n
,
where A
=
1
+
(
Δ
t
/
C
g
)
(
C
R
-
C
H
-
C
L
b
d
-
C
D
), B
=
C
L
d
T
0
(
b
2
Δ
t
/
2
C
g
), C
=
Δ
t
(
C
D
/
C
g
), and D
=
Δ
t
/
τ. Introducing the replacement x
1
,
n
=
A
x
n
/
B, where x
1 is the modified dimensionless environmental interface temperature and x
2
,
n
=
y
n, we can write the following:
(10)
x
1
,
n
+
1
=
A
x
1
,
n
(
1
-
x
1
,
n
)
+
C
B
A
x
2
,
n
,
x
2
,
n
+
1
=
D
A
B
x
1
,
n
+
(
1
-
D
)
x
2
,
n
.
Analysis of values of parameters A, B, C, and D, based on a large number of energy flux outputs from the land surface scheme runs, indicates that their values are ranged in the following intervals: (i) A
∈
[
0,4
] and (ii) B, C, and D are ranged in the interval [
0,1
]. Thus, A is the logistic parameter, which from now will be denoted by r. All other groups of parameters in the system (10) have the values in the same interval [
0,1
]. Let us underline that under some circumstances those parameters can be equal. Correspondingly, we replaced all of them by introducing the coupling parameter c.Finally, system (10) can be written in the form of coupled maps; that is,
(11)
x
1
,
n
+
1
=
r
x
1
,
n
(
1
-
x
1
,
n
)
+
c
x
2
,
n
,
x
2
,
n
+
1
=
c
(
x
1
,
n
+
x
2
,
n
)
.Now we analyse the stability of physical solutions of coupled maps (11), using the Lyapunov exponent, which is a measure of convergence or divergence of near trajectories in phase space. Sign of Lyapunov exponent is characteristic of attractor type and for stable fixed point it is negative, although for chaotic attractor it is positive. Calculating Lyapunov exponent for the coupled system (11) with values of parameters c
∈
(
0.05,0.1
) and r
∈
(
3.6,3.8
), because we thought that it will be interesting to investigate behavior of the system for small values of coupling parameter and high values of logistic parameter, we got results depicted in Figure 2 [53]. It is shown that the Lyapunov exponent mostly has positive values approving presence of chaos in this system, but there are still some strait regions where the Lyapunov exponent is negative and where the solutions of the coupled system are stable, that is, domains of stability.Figure 2
Lyapunov exponent of the coupled system (11), which shows presence of strait regions of stability in highly developed chaos.Irregularities in solution of the system of difference equation (11) can come from two reasons. They are (i)numerical, that is, because we try to choose appropriate difference equation whose solution is “good” approximation to the solution of the given partial differential equation and (ii)physical, that is, occurrence of chaotic fluctuations in the considered system because the environmental interface cannot oppose an enormous radiative forcing, suddenly reaching the interface. Therefore, it raises the question whether we can find either domain or domains where physically meaningful solutions exist [44]. We do that by considering the stability of physical solution of (11). Stability, in mathematics, is condition in which a slight disturbance in a system does not produce too disrupting effect on that system. In terms of the solution of a differential equation, a function f
(
ζ
) is said to be stable if any other solution of the equation that starts out sufficiently close to it when ζ
=
0 remains close to it for succeeding values of ζ. If the difference between the solutions approaches zero as ζ increases, the solution is called asymptotically stable. If a solution does not have either of these properties, it is called unstable. We consider the stability of physical solution of (11) in sense of the Lyapunov exponent (λ). Thus, the system is stable if λ
<
0 and unstable if λ
>
0. This analysis, set into context of the climate modeling, points out the fact that there exists set of domains where the environmental interface temperature cannot be calculated by the physics of currently designed climate models, because of occurrence of chaotic phenomena on the environmental interfaces.


## 2.2.1. Background

Traditional mathematical analysis of physical systems tacitly assumes that integers and all real numbers, no matter how large or how small, are physically possible and all mathematically possible trajectories are physically possible [44]. Traditionally, this approach has worked well in physics and in engineering, but it does not lead to a very good understanding of chaotic systems, which, as is now known, are extremely important in the study of real-world phenomena ranging from weather to biological systems. In this paper we deal with one issue in modeling pathways in meteorology as well as in physics, biology, and chemistry, that is, environmental sciences in their broadest context [45], in particular in autonomous dynamical systems, which are common subject under consideration in climate modeling. Namely, we consider how to replace given differential equations by appropriate difference equations in environmental modeling and thus in climate simulations [46].According to van der Vaart many models for environmental problems have been and will be built in the form of differential equations or systems of such equations [46]. With the advent of computers one has been able to find (approximate) solutions for equations that used to be intractable. Many of the mathematical techniques have been used in this area to replace given differential equations by appropriate difference equations. So a huge effort has been invested into choice of appropriate difference equations whose solutions are “good” approximations to the solutions of the given differential equations. This question includes a requirement for better understanding of the fundamental problem: interrelations between classical continuum mathematics and reality in different sciences. For many atmospheric phenomena the “continuum” type of thinking, that is, at the basis of any differential equation, is not natural to the phenomenon but rather constitutes an approximation to a basically discrete situation: in much work of this type the “infinitesimal step lengths” handled in the reasoning which lead us to the differential equation are not really thought of as infinitesimally small but as finite; yet, in the last stage of such reasoning, where the differential equation rises from the differentials, these “infinitesimal” step lengths go to zero, that is, where the above-mentioned approximation comes in. Under this kind of circumstances, it seems more naturalto build the model as a discrete difference equation from the start, without going through the painful, doubly approximative process of first, during the modeling stage, finding a differential equation to approximate a basically discrete situation and then, for numerical computing purposes, approximating that differential equation by a difference scheme [44].In this subsection we analyze the energy balance equation in procedure of computing the environmental interface temperature and the deeper soil layer temperature commonly used in climate models. The environmental interface is defined asinterface between two biotic or abiotic environments that are in relative motion and exchange energy, matter, and information through physical, biological, and chemical processes, fluctuating temporally and spatially regardless of space and time scale[47].There are a lot of examples of environmental interfaces in the nature, but here we deal with the ground surface, where there exist all three mechanisms of energy transfer: incoming and outgoing radiation, convection of heat and moisture into the atmosphere, and conduction of heat into deeper soil layers of ground (Figure1) [48]. Parameterization of these processes is of great importance for environmental models of different spatial and temporal scales and thus climate ones. In the paper by Mihailović and Mimić (2012) it is shown that ground surface is treated as a complex system in which chaotic fluctuations occur while we compute its temperature [49]. This system, as an actual dynamic system, is very sensitive to initial conditions and arbitrarily small perturbation of the current trajectory that may lead to its unpredictable behavior. In the aforementioned paper the lower boundary condition, that is, the deeper soil layer temperature, was constant, but it can also vary in time. That system, often used in environmental models, is of interest to be analyzed by the methods of nonlinear dynamics. Having in mind those facts, in this, we (i) perform a nonlinear dynamical analysis of coupled system for computing the environmental interface temperature and the deeper soil layer temperature and (ii) examine behavior of the coupled system in dependence on the main system parameters, in order to show the possible occurrence of the chaos in computing the environmental interface temperature. Firstly, we consider difference form of the energy balance equation and deeper soil layer temperature equation transforming them into the coupled system with the corresponding parameters, then we analyze behavior of the solutions of the coupled system, and we have examined domains of stability using the Lyapunov exponent.Figure 1
Energy balance equation terms.


## 2.2.2. Physical Background and Derivation of the Coupled System

One of the most important conditions for functioning of any complex system is a proper supply of the system with energy. Dynamics of energy flow is based on the energy balance equation [48]. As we mentioned before, environmental interface is a complex system. General difference form of energy balance equation for the ground surface as an environmental interface is
(1)
C
g
Δ
T
g
Δ
t
=
R
net
-
H
-
λ
E
-
G
,
where T
g is the ground surface temperature, Δ
t is the time step, C
g is the soil heat capacity, R
net is the net radiation, H is the sensible heat flux, λ
E is the latent heat flux, and G is the heat flux into the ground. First, we assume that the net radiation is given as in [50]; that is,
(2)
R
net
=
C
R
(
T
g
-
T
a
)
,
where T
a is the air temperature at some reference level and C
R is the coefficient for the net radiation term. Second, we make expansion of the exponential term in the expression for latent heat flux as follows:
(3)
λ
E
=
C
L
d
[
b
(
T
g
-
T
a
)
+
b
2
2
(
T
g
-
T
a
)
2
]
,
where C
L is the water vapour transfer coefficient, b
=
0.06337
°
C
-
1, and d is parameter which occurs in expanding the series [51]. Further, the conduction of the heat into the soil can be written in the following form:
(4)
G
=
C
D
(
T
g
-
T
d
)
,
where C
D is the coefficient of the heat conduction while T
d is the temperature of the deeper soil layer. The sensible heat flux H can be parameterized as
(5)
H
=
C
H
(
T
g
-
T
a
)
,
where C
H is the sensible heat transfer coefficient. The prognostic equation for temperature of the deeper soil layer T
d is
(6)
Δ
T
d
Δ
t
=
1
τ
(
T
g
-
T
d
)
,
where τ
=
86400 s. After collecting all terms (2)–(6), the coupled system takes the following form:
(7)
C
g
Δ
T
g
Δ
t
=
C
R
(
T
g
-
T
a
)
-
C
H
(
T
g
-
T
a
)
-
C
L
d
[
b
(
T
g
-
T
a
)
+
b
2
2
(
T
g
-
T
a
)
2
]
-
C
D
(
T
g
-
T
d
)
,
Δ
T
d
Δ
t
=
1
τ
(
T
g
-
T
d
)
.
More details about the nature and the range of physical parameters C
R, C
L, C
D, and C
H can be found in [52]. Now, using the time scheme forward in time (n indicates the time step) and dividing both sides of (7) with the constant temperature T
0 (e.g., value of mean Earth temperature, i.e., T
0
=
288 K) we get
(8)
T
g
n
+
1
-
T
a
T
0
=
T
g
n
-
T
a
T
0
+
Δ
t
C
g
C
R
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
H
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
L
b
d
T
g
n
-
T
a
T
0
-
Δ
t
C
g
C
L
d
T
0
b
2
2
(
T
g
n
-
T
a
)
2
T
0
2
-
Δ
t
C
g
C
D
T
g
n
-
T
a
T
0
+
Δ
t
C
g
C
D
T
d
n
-
T
a
T
0
,
T
d
n
+
1
-
T
a
T
0
=
T
d
n
-
T
a
T
0
+
Δ
t
τ
T
g
n
-
T
a
T
0
-
Δ
t
τ
T
d
n
-
T
a
T
0
.
Finally, introducing replacements x
=
(
T
g
-
T
a
)
/
T
0 and y
=
(
T
d
-
T
a
)
/
T
0, where x is the dimensionless environmental interface temperature and y is the dimensionless deeper soil layer temperature, we reach the coupled system as follows:
(9)
x
n
+
1
=
A
x
n
-
B
x
n
2
+
C
y
n
,
y
n
+
1
=
D
x
n
+
(
1
-
D
)
y
n
,
where A
=
1
+
(
Δ
t
/
C
g
)
(
C
R
-
C
H
-
C
L
b
d
-
C
D
), B
=
C
L
d
T
0
(
b
2
Δ
t
/
2
C
g
), C
=
Δ
t
(
C
D
/
C
g
), and D
=
Δ
t
/
τ. Introducing the replacement x
1
,
n
=
A
x
n
/
B, where x
1 is the modified dimensionless environmental interface temperature and x
2
,
n
=
y
n, we can write the following:
(10)
x
1
,
n
+
1
=
A
x
1
,
n
(
1
-
x
1
,
n
)
+
C
B
A
x
2
,
n
,
x
2
,
n
+
1
=
D
A
B
x
1
,
n
+
(
1
-
D
)
x
2
,
n
.
Analysis of values of parameters A, B, C, and D, based on a large number of energy flux outputs from the land surface scheme runs, indicates that their values are ranged in the following intervals: (i) A
∈
[
0,4
] and (ii) B, C, and D are ranged in the interval [
0,1
]. Thus, A is the logistic parameter, which from now will be denoted by r. All other groups of parameters in the system (10) have the values in the same interval [
0,1
]. Let us underline that under some circumstances those parameters can be equal. Correspondingly, we replaced all of them by introducing the coupling parameter c.Finally, system (10) can be written in the form of coupled maps; that is,
(11)
x
1
,
n
+
1
=
r
x
1
,
n
(
1
-
x
1
,
n
)
+
c
x
2
,
n
,
x
2
,
n
+
1
=
c
(
x
1
,
n
+
x
2
,
n
)
.Now we analyse the stability of physical solutions of coupled maps (11), using the Lyapunov exponent, which is a measure of convergence or divergence of near trajectories in phase space. Sign of Lyapunov exponent is characteristic of attractor type and for stable fixed point it is negative, although for chaotic attractor it is positive. Calculating Lyapunov exponent for the coupled system (11) with values of parameters c
∈
(
0.05,0.1
) and r
∈
(
3.6,3.8
), because we thought that it will be interesting to investigate behavior of the system for small values of coupling parameter and high values of logistic parameter, we got results depicted in Figure 2 [53]. It is shown that the Lyapunov exponent mostly has positive values approving presence of chaos in this system, but there are still some strait regions where the Lyapunov exponent is negative and where the solutions of the coupled system are stable, that is, domains of stability.Figure 2
Lyapunov exponent of the coupled system (11), which shows presence of strait regions of stability in highly developed chaos.Irregularities in solution of the system of difference equation (11) can come from two reasons. They are (i)numerical, that is, because we try to choose appropriate difference equation whose solution is “good” approximation to the solution of the given partial differential equation and (ii)physical, that is, occurrence of chaotic fluctuations in the considered system because the environmental interface cannot oppose an enormous radiative forcing, suddenly reaching the interface. Therefore, it raises the question whether we can find either domain or domains where physically meaningful solutions exist [44]. We do that by considering the stability of physical solution of (11). Stability, in mathematics, is condition in which a slight disturbance in a system does not produce too disrupting effect on that system. In terms of the solution of a differential equation, a function f
(
ζ
) is said to be stable if any other solution of the equation that starts out sufficiently close to it when ζ
=
0 remains close to it for succeeding values of ζ. If the difference between the solutions approaches zero as ζ increases, the solution is called asymptotically stable. If a solution does not have either of these properties, it is called unstable. We consider the stability of physical solution of (11) in sense of the Lyapunov exponent (λ). Thus, the system is stable if λ
<
0 and unstable if λ
>
0. This analysis, set into context of the climate modeling, points out the fact that there exists set of domains where the environmental interface temperature cannot be calculated by the physics of currently designed climate models, because of occurrence of chaotic phenomena on the environmental interfaces.


## 2.3. Horizontal Energy Exchange between Environmental Interfaces

### 2.3.1. Background

There are three major sets of processes that must be considered when constructing a climate model: (i) radiative (the transfer of radiation through the climate system, e.g., absorption and reflection); (ii) dynamic (the horizontal and vertical transfer of energy, e.g., advection, convection, and diffusion); and (iii) surface process (inclusion of processes involving land/ocean/ice and the effects of albedo, emissivity, and surface-atmosphere energy exchanges). If the nonlinearities in these processes are treated improperly, then while designing the model, the complexity and thus its reliability will not be retained in the highest degree. In Section2.2 we have considered surface-atmosphere energy exchanges with cadence on the phenomenon of a possible occurrence of the chaos in solving the energy balance equation for calculating the environmental interface temperature in climate models. Here, following Mihailović et al. (2012) we analyze the horizontal energy exchange between environmental interfaces which is described by the dynamics of driven coupled oscillators [54]. In order to study their behavior, when a perturbation is introduced in the system, as a function of the coupling parameter, the logistic parameter, and the horizontal energy exchange intensity (parameter of exchange, in further text), we considered dynamics of two maps serving the diffusive coupling [54].As noted above, the horizontal exchange of energy between environmental interfaces is considered as diffusion-like process. The dynamics of energy exchange behavior on environmental interface are typically expressed as a logistic mapΦ
(
x
)
=
r
x
(
1
-
x
), where x is the dimensionless temperature of environmental interface and r is a logistic parameter [53, 54]. However, we use an alternative form of this map, which includes a parameter p that represents the horizontal energy exchange intensity (Figure 3). By introducing this parameter we formalize an intrinsic property of the environmental interfaces, which depends on the nature of the interface. The environmental interface dynamics are expressed here as a difference equation, so we avoid the double approximation of (i) finding a differential equation to approximate an essentially discrete process (during the modeling stage) and then (ii) approximating that differential equation by a difference scheme for numerical computing purposes [44, 46]; that is,
(12)
Φ
(
x
i
,
n
)
=
r
x
i
,
n
p
(
1
-
x
i
,
n
p
)
.
The dynamics of this map (12) are governed by two parameters, p and r, which express intrinsic property of the environmental interfaces and the influence of the environment, respectively.Figure 3
Schematic diagram of horizontal energy exchange between two environmental interfaces. Parametersp and r express intrinsic property of the environmental interfaces and the influence of the environment, respectively.Since these and many other processes on environmental interface are defined as diffusion-like, we will explore (i) how these processes can be better represented in climate models by introducing parameter of exchangep in the diffusive coupling associated with the horizontal energy exchange and (ii) how the horizontal energy exchange intensity dynamics are affected by the perturbation of parameters that represents influence of the environment, environmental interface coupling, and horizontal energy exchange intensity.In considering these problems we have to include observational heterarchy, a challenging topic when dealing with complex systems. Essentially, observational heterarchy reveals that it is impossible to unambiguously determine to which subsystems an element belongs [55]. Therefore, the dynamics of the complex system are articulated in terms of two kinds of dynamics, Intent and Extent dynamics and the interaction between them, where Intent corresponds to an attribute of a given phenomenon and Extent corresponds to a collection of objects satisfying that phenomenon [55].


### 2.3.2. Observational Heterarchy and Horizontal Energy Exchange between Environmental Interfaces

Observational heterarchy consists of two sets of intralayer maps, called Intent and Extent perspectives, and interlayer operations satisfying the following conditions. (1) The interlayer operations inherit the mixture of intra- and interlayer operations and (2) there is a procedure by which the interlayer operation can be regarded as an adjoint functor. If the interlayer operation satisfies the conditions (1) and (2), it is called a prefunctor [55]. Preserving the above-mentioned composition occurs as follows: a prefunctor, 〈
F
〉
: Intent → Extent is mapping a set, X, to a set, 〈
F
〉
X, and map Φ to a map, f
*
Φ
f, where f
*
f
(
x
)
=
x for all x in f
(
X
) with f
(
X
)
:
〈
F
〉
X
→
X. In this sense we call f
* pseudoinverse of  f. Because applying a prefunctor to a map is expressed as composition of maps, it satisfies the conditions (1) and (2). The approximation is defined by the assumption that f is a one-to-one onto map. If one accepts the approximation, f
*
=
f
-
1 holds, then a prefunctor can become a functor. Given two maps,Φ
,
Ψ
:
X
→
X,
(13)
〈
F
〉
(
Φ
)
〈
F
〉
(
Ψ
)
=
(
f
*
Φ
f
)
(
f
*
Ψ
f
)
=
f
*
Φ
(
f
f
*
)
Ψ
f
=
f
*
Φ
(
f
f
-
1
)
Ψ
f
=
f
*
Φ
Ψ
f
=
〈
F
〉
(
Φ
f
)
.
It implies that F preserves the composition of maps, Φ and Ψ.The time development of the environmental surface dynamicsx
i
,
n, for two interfaces, is expressed as
(14)
x
i
,
n
+
1
=
(
1
-
c
)
Φ
(
x
i
,
n
)
+
f
(
Φ
(
x
j
,
n
)
)
,
where n is the time iteration, i
,
j
=
1,2, x
i
,
n
∈
[
0,1
], c is the coupling parameter, f is the map representing the horizontal energy exchange between environmental interfaces, and Φ is one of maps in the pair (
Ψ
,
Φ
) whose composition is preserved by a prefunctor 〈
F
〉. Here, we apply the framework of an observational heterarchy to the two environmental interface systems. If Intent and Extent are denoted by Φ and Ψ, respectively, the time development of the concentration is expressed as x
i
,
n
+
1
=
(
1
-
c
)
Φ
(
x
i
,
n
)
+
Ψ
(
x
j
,
n
). In this expression, if Ψ
(
X
)
=
f
(
Φ
(
x
)
), then it can be reduced to (14).We perform our analysis following the procedure described in [55]. First, in this subsection we address the synchronization of the passive coupling for two environmental interfaces given by (14) and (12), and then, in the next subsection, we will show that perturbation can modify the dynamics and enhance robust behavior in a multienvironmental interface system of active coupling. Synchronization is well-known collective phenomenon in various multicomponent physical as well as the climate systems [56–58]. The exchange of information (coupling) among the components can be either global or local. Here, we consider that chaotic systems are synchronized only when the largest Lyapunov exponent of the driven system is negative. It was calculated by approach proposed in [59]. We studied the stability of the fixed point by linearizing n
≥
2 component coupled system and obtain Z
n
+
1
=
ζ
n
Z
n where ζ
n is the Jacobian of this system evaluated in (
0,0
,
…
,
0
) and ζ
n
=
(
x
1
,
n
,
x
2
,
n
,
…
,
x
N
,
n
). By iterating we obtain
(15)
Z
n
+
1
=
(
∏
s
=
0
n
ζ
s
)
Z
0
,
and thus we get Lyapunov exponent as follows:
(16)
λ
=
lim
⁡
n
→
∞
⁡
(
ln
⁡
∥
∏
s
=
0
n
ζ
s
∥
)
n
.
Figure 4 depicts the diagram of normalized frequency of synchronization F
p (λ
<
0) for system of two environmental interfaces passively coupled ((12) and (14)), as a function of parameter of exchange p, averaged over all values of the coupling parameter c and logistic parameter r. The value of the normalized frequency of synchronization F
p is calculated as
(17)
F
p
=
∑
N
n
(
λ
<
0
)
∑
N
n
(
λ
<
0
)
+
∑
N
p
(
λ
>
0
)
,
where N
n
(
λ
<
0
) and N
p
(
λ
>
0
) are numbers of negative and positive values of the Lyapunov exponent, respectively. These numbers were calculated for the fixed value of p, while c and r were changing in intervals (
0,1
) and (
1,4
), respectively, with the step of 0.05. From this figure it is seen that after p
>
0.2, F
p becomes lower, indicating a decrease of number of states, which are synchronized.Figure 4
Normalized frequency of synchronization,F
p
(
λ
<
0
), for system of environmental interfaces passively coupled (12)–(14) as a function of parameter of exchange p. An averaging was done over all values of coupling parameter c and logistic parameter r.


### 2.3.3. Simulations of Active Coupling in a Multienvironmental Interface System

Here, we address the behavior of active coupling [55] and estimate whether a coupled map system described above can achieve synchronization under influence of perturbations. The dynamics of two environmental interface systems called active coupling [55], used for simulations, is expressed as(18a)
x
i
,
n
+
1
=
(
1
-
c
)
Φ
n
(
x
i
,
n
)
+
Ψ
n
(
x
j
,
n
)
,
(18b)
Ψ
n
+
1
=
f
Φ
n
f
*
f
,
(18c)
Φ
n
+
1
=
f
*
Ψ
n
+
1
,
(18d)
Φ
n
(
x
i
,
n
)
=
r
x
i
,
n
p
(
1
-
x
i
,
n
p
)
.We note that the dynamical system defined by (18a) and (18d) is called the passive coupling, and that is a usual coupled map system. The active coupling can be approximated to passive coupling, where the approximation is defined by the adjunction or the equivalence between Intent and Extent. Compared with passive coupling, the behavior of active coupling is much more complex [55]. In (18a)–(18d), because of a pseudoinverse map, f
*, all calculations are defined to be approximations. In simulations, the Intent map was a discontinuous map, expressed by Φ
n
+
1
=
f
*
Ψ
n
+
1.In order to see how perturbation enhances robust behavior in the framework of observational heterarchy in a multienvironmental interface system represented by closed contour of coupled environmental interfaces exchanging the energy horizontally, then the system of coupled difference equations forN environmental interfaces exchanging the energy can be written in the form of the following matrix equation:(19a)
X
N
1
=
(
A
+
B
)
·
X
N
.
The elements in matrices in (19a) are
(19b)
X
N
1
i
,
n
+
1
=
x
i
,
n
+
1
,
X
N
i
,
n
=
x
i
,
n
,
A
i
,
k
=
(
1
-
c
)
Φ
n
(
x
i
,
n
)
δ
i
,
k
,
(19c)
B
i
,
k
=
{
Ψ
n
(
x
k
,
n
)
,
k
=
i
+
1
,
i
<
N
0
,
k
≠
i
+
1
,
i
<
N
Ψ
n
(
x
k
,
n
)
,
k
=
1
,
i
=
N
0
,
k
≠
1
,
i
=
N
,where i
=
1,2
,
…
,
N and δ
i
,
k is the Kronecker symbol.Simulations with the active coupling, defined by (18a)–(18d), were performed with and without perturbation given as in [56]. The results of simulations are shown in Figure 5. In this figure Lyapunov exponent λ is plotted against coupling parameter c for active coupling with perturbation (black line) compared to the passive coupling (red line), for different values of the parameter of exchange pand the logistic parameter r. Simulations were performed with the closed contour of N
=
100 interfaces. The Lyapunov exponent was calculated using (15)-(16) and the Jacobian of the system given by (19b)-(19c) is representing this contour.Figure 5
Diagram of Lyapunov exponent,λ, against coupling parameter c for the fluctuated active coupling defined by (18a)–(18d) with P (black line) compared to passive coupling and N (red line) for different values of parameter of exchange p and logistic parameter r. In (a)–(c) p takes the fixed values (1.0, 0.5, 0.2), while r
=
4. In (d)–(f) p is randomly chosen, while r takes values 4.0, 3.82, and 3.6, respectively. Simulations were performed with the closed contour of N
=
100 environmental interfaces.In calculatingλ, for each c from 0.0 to 1.0 with step 0.005, 1
0
4 iterations were applied for an initial state, and then the first 1
0
3 steps were abandoned. In order to see how the active coupling modifies the synchronization of horizontal energy exchange between environmental interfaces, we performed two kinds of simulations. Firstly, we used r
=
4.0 and the fixed value of the parameter of exchange p (Figures 5(a)–5(c)); secondly, we used a randomly chosen p and a logistic r parameter with the values of 4.0, 3.82, and 3.6, respectively (Figures 5(d)–5(f)). Figures 5(a)–5(c) depict that, in the chaotic regime (r
=
4.0), regardless of the value p, the Lyapunov exponent is always positive (λ
>
0) and therefore the process of the horizontal energy exchange in a multienvironmental interface system is always unsynchronized. However, the stormy perturbation disturbs this state (Figures 5(a)–5(c)). Although the logistic parameter is settled at r
=
4.0 for chaotic behavior, the coupling parameter c tunes interaction and leads to synchronization in some intervals, particularly for p
=
1.0 and p
=
0.5. This behavior is more pronounced in Figures 5(d)–5(f) where p is randomly chosen; here the process of horizontal energy exchange in multienvironmental interfaces exhibits a strong tendency towards the synchronization, even though the logistic parameter r is in chaotic region (r
=
4.0,3.82 and 3.6).


## 2.3.1. Background

There are three major sets of processes that must be considered when constructing a climate model: (i) radiative (the transfer of radiation through the climate system, e.g., absorption and reflection); (ii) dynamic (the horizontal and vertical transfer of energy, e.g., advection, convection, and diffusion); and (iii) surface process (inclusion of processes involving land/ocean/ice and the effects of albedo, emissivity, and surface-atmosphere energy exchanges). If the nonlinearities in these processes are treated improperly, then while designing the model, the complexity and thus its reliability will not be retained in the highest degree. In Section2.2 we have considered surface-atmosphere energy exchanges with cadence on the phenomenon of a possible occurrence of the chaos in solving the energy balance equation for calculating the environmental interface temperature in climate models. Here, following Mihailović et al. (2012) we analyze the horizontal energy exchange between environmental interfaces which is described by the dynamics of driven coupled oscillators [54]. In order to study their behavior, when a perturbation is introduced in the system, as a function of the coupling parameter, the logistic parameter, and the horizontal energy exchange intensity (parameter of exchange, in further text), we considered dynamics of two maps serving the diffusive coupling [54].As noted above, the horizontal exchange of energy between environmental interfaces is considered as diffusion-like process. The dynamics of energy exchange behavior on environmental interface are typically expressed as a logistic mapΦ
(
x
)
=
r
x
(
1
-
x
), where x is the dimensionless temperature of environmental interface and r is a logistic parameter [53, 54]. However, we use an alternative form of this map, which includes a parameter p that represents the horizontal energy exchange intensity (Figure 3). By introducing this parameter we formalize an intrinsic property of the environmental interfaces, which depends on the nature of the interface. The environmental interface dynamics are expressed here as a difference equation, so we avoid the double approximation of (i) finding a differential equation to approximate an essentially discrete process (during the modeling stage) and then (ii) approximating that differential equation by a difference scheme for numerical computing purposes [44, 46]; that is,
(12)
Φ
(
x
i
,
n
)
=
r
x
i
,
n
p
(
1
-
x
i
,
n
p
)
.
The dynamics of this map (12) are governed by two parameters, p and r, which express intrinsic property of the environmental interfaces and the influence of the environment, respectively.Figure 3
Schematic diagram of horizontal energy exchange between two environmental interfaces. Parametersp and r express intrinsic property of the environmental interfaces and the influence of the environment, respectively.Since these and many other processes on environmental interface are defined as diffusion-like, we will explore (i) how these processes can be better represented in climate models by introducing parameter of exchangep in the diffusive coupling associated with the horizontal energy exchange and (ii) how the horizontal energy exchange intensity dynamics are affected by the perturbation of parameters that represents influence of the environment, environmental interface coupling, and horizontal energy exchange intensity.In considering these problems we have to include observational heterarchy, a challenging topic when dealing with complex systems. Essentially, observational heterarchy reveals that it is impossible to unambiguously determine to which subsystems an element belongs [55]. Therefore, the dynamics of the complex system are articulated in terms of two kinds of dynamics, Intent and Extent dynamics and the interaction between them, where Intent corresponds to an attribute of a given phenomenon and Extent corresponds to a collection of objects satisfying that phenomenon [55].


## 2.3.2. Observational Heterarchy and Horizontal Energy Exchange between Environmental Interfaces

Observational heterarchy consists of two sets of intralayer maps, called Intent and Extent perspectives, and interlayer operations satisfying the following conditions. (1) The interlayer operations inherit the mixture of intra- and interlayer operations and (2) there is a procedure by which the interlayer operation can be regarded as an adjoint functor. If the interlayer operation satisfies the conditions (1) and (2), it is called a prefunctor [55]. Preserving the above-mentioned composition occurs as follows: a prefunctor, 〈
F
〉
: Intent → Extent is mapping a set, X, to a set, 〈
F
〉
X, and map Φ to a map, f
*
Φ
f, where f
*
f
(
x
)
=
x for all x in f
(
X
) with f
(
X
)
:
〈
F
〉
X
→
X. In this sense we call f
* pseudoinverse of  f. Because applying a prefunctor to a map is expressed as composition of maps, it satisfies the conditions (1) and (2). The approximation is defined by the assumption that f is a one-to-one onto map. If one accepts the approximation, f
*
=
f
-
1 holds, then a prefunctor can become a functor. Given two maps,Φ
,
Ψ
:
X
→
X,
(13)
〈
F
〉
(
Φ
)
〈
F
〉
(
Ψ
)
=
(
f
*
Φ
f
)
(
f
*
Ψ
f
)
=
f
*
Φ
(
f
f
*
)
Ψ
f
=
f
*
Φ
(
f
f
-
1
)
Ψ
f
=
f
*
Φ
Ψ
f
=
〈
F
〉
(
Φ
f
)
.
It implies that F preserves the composition of maps, Φ and Ψ.The time development of the environmental surface dynamicsx
i
,
n, for two interfaces, is expressed as
(14)
x
i
,
n
+
1
=
(
1
-
c
)
Φ
(
x
i
,
n
)
+
f
(
Φ
(
x
j
,
n
)
)
,
where n is the time iteration, i
,
j
=
1,2, x
i
,
n
∈
[
0,1
], c is the coupling parameter, f is the map representing the horizontal energy exchange between environmental interfaces, and Φ is one of maps in the pair (
Ψ
,
Φ
) whose composition is preserved by a prefunctor 〈
F
〉. Here, we apply the framework of an observational heterarchy to the two environmental interface systems. If Intent and Extent are denoted by Φ and Ψ, respectively, the time development of the concentration is expressed as x
i
,
n
+
1
=
(
1
-
c
)
Φ
(
x
i
,
n
)
+
Ψ
(
x
j
,
n
). In this expression, if Ψ
(
X
)
=
f
(
Φ
(
x
)
), then it can be reduced to (14).We perform our analysis following the procedure described in [55]. First, in this subsection we address the synchronization of the passive coupling for two environmental interfaces given by (14) and (12), and then, in the next subsection, we will show that perturbation can modify the dynamics and enhance robust behavior in a multienvironmental interface system of active coupling. Synchronization is well-known collective phenomenon in various multicomponent physical as well as the climate systems [56–58]. The exchange of information (coupling) among the components can be either global or local. Here, we consider that chaotic systems are synchronized only when the largest Lyapunov exponent of the driven system is negative. It was calculated by approach proposed in [59]. We studied the stability of the fixed point by linearizing n
≥
2 component coupled system and obtain Z
n
+
1
=
ζ
n
Z
n where ζ
n is the Jacobian of this system evaluated in (
0,0
,
…
,
0
) and ζ
n
=
(
x
1
,
n
,
x
2
,
n
,
…
,
x
N
,
n
). By iterating we obtain
(15)
Z
n
+
1
=
(
∏
s
=
0
n
ζ
s
)
Z
0
,
and thus we get Lyapunov exponent as follows:
(16)
λ
=
lim
⁡
n
→
∞
⁡
(
ln
⁡
∥
∏
s
=
0
n
ζ
s
∥
)
n
.
Figure 4 depicts the diagram of normalized frequency of synchronization F
p (λ
<
0) for system of two environmental interfaces passively coupled ((12) and (14)), as a function of parameter of exchange p, averaged over all values of the coupling parameter c and logistic parameter r. The value of the normalized frequency of synchronization F
p is calculated as
(17)
F
p
=
∑
N
n
(
λ
<
0
)
∑
N
n
(
λ
<
0
)
+
∑
N
p
(
λ
>
0
)
,
where N
n
(
λ
<
0
) and N
p
(
λ
>
0
) are numbers of negative and positive values of the Lyapunov exponent, respectively. These numbers were calculated for the fixed value of p, while c and r were changing in intervals (
0,1
) and (
1,4
), respectively, with the step of 0.05. From this figure it is seen that after p
>
0.2, F
p becomes lower, indicating a decrease of number of states, which are synchronized.Figure 4
Normalized frequency of synchronization,F
p
(
λ
<
0
), for system of environmental interfaces passively coupled (12)–(14) as a function of parameter of exchange p. An averaging was done over all values of coupling parameter c and logistic parameter r.


## 2.3.3. Simulations of Active Coupling in a Multienvironmental Interface System

Here, we address the behavior of active coupling [55] and estimate whether a coupled map system described above can achieve synchronization under influence of perturbations. The dynamics of two environmental interface systems called active coupling [55], used for simulations, is expressed as(18a)
x
i
,
n
+
1
=
(
1
-
c
)
Φ
n
(
x
i
,
n
)
+
Ψ
n
(
x
j
,
n
)
,
(18b)
Ψ
n
+
1
=
f
Φ
n
f
*
f
,
(18c)
Φ
n
+
1
=
f
*
Ψ
n
+
1
,
(18d)
Φ
n
(
x
i
,
n
)
=
r
x
i
,
n
p
(
1
-
x
i
,
n
p
)
.We note that the dynamical system defined by (18a) and (18d) is called the passive coupling, and that is a usual coupled map system. The active coupling can be approximated to passive coupling, where the approximation is defined by the adjunction or the equivalence between Intent and Extent. Compared with passive coupling, the behavior of active coupling is much more complex [55]. In (18a)–(18d), because of a pseudoinverse map, f
*, all calculations are defined to be approximations. In simulations, the Intent map was a discontinuous map, expressed by Φ
n
+
1
=
f
*
Ψ
n
+
1.In order to see how perturbation enhances robust behavior in the framework of observational heterarchy in a multienvironmental interface system represented by closed contour of coupled environmental interfaces exchanging the energy horizontally, then the system of coupled difference equations forN environmental interfaces exchanging the energy can be written in the form of the following matrix equation:(19a)
X
N
1
=
(
A
+
B
)
·
X
N
.
The elements in matrices in (19a) are
(19b)
X
N
1
i
,
n
+
1
=
x
i
,
n
+
1
,
X
N
i
,
n
=
x
i
,
n
,
A
i
,
k
=
(
1
-
c
)
Φ
n
(
x
i
,
n
)
δ
i
,
k
,
(19c)
B
i
,
k
=
{
Ψ
n
(
x
k
,
n
)
,
k
=
i
+
1
,
i
<
N
0
,
k
≠
i
+
1
,
i
<
N
Ψ
n
(
x
k
,
n
)
,
k
=
1
,
i
=
N
0
,
k
≠
1
,
i
=
N
,where i
=
1,2
,
…
,
N and δ
i
,
k is the Kronecker symbol.Simulations with the active coupling, defined by (18a)–(18d), were performed with and without perturbation given as in [56]. The results of simulations are shown in Figure 5. In this figure Lyapunov exponent λ is plotted against coupling parameter c for active coupling with perturbation (black line) compared to the passive coupling (red line), for different values of the parameter of exchange pand the logistic parameter r. Simulations were performed with the closed contour of N
=
100 interfaces. The Lyapunov exponent was calculated using (15)-(16) and the Jacobian of the system given by (19b)-(19c) is representing this contour.Figure 5
Diagram of Lyapunov exponent,λ, against coupling parameter c for the fluctuated active coupling defined by (18a)–(18d) with P (black line) compared to passive coupling and N (red line) for different values of parameter of exchange p and logistic parameter r. In (a)–(c) p takes the fixed values (1.0, 0.5, 0.2), while r
=
4. In (d)–(f) p is randomly chosen, while r takes values 4.0, 3.82, and 3.6, respectively. Simulations were performed with the closed contour of N
=
100 environmental interfaces.In calculatingλ, for each c from 0.0 to 1.0 with step 0.005, 1
0
4 iterations were applied for an initial state, and then the first 1
0
3 steps were abandoned. In order to see how the active coupling modifies the synchronization of horizontal energy exchange between environmental interfaces, we performed two kinds of simulations. Firstly, we used r
=
4.0 and the fixed value of the parameter of exchange p (Figures 5(a)–5(c)); secondly, we used a randomly chosen p and a logistic r parameter with the values of 4.0, 3.82, and 3.6, respectively (Figures 5(d)–5(f)). Figures 5(a)–5(c) depict that, in the chaotic regime (r
=
4.0), regardless of the value p, the Lyapunov exponent is always positive (λ
>
0) and therefore the process of the horizontal energy exchange in a multienvironmental interface system is always unsynchronized. However, the stormy perturbation disturbs this state (Figures 5(a)–5(c)). Although the logistic parameter is settled at r
=
4.0 for chaotic behavior, the coupling parameter c tunes interaction and leads to synchronization in some intervals, particularly for p
=
1.0 and p
=
0.5. This behavior is more pronounced in Figures 5(d)–5(f) where p is randomly chosen; here the process of horizontal energy exchange in multienvironmental interfaces exhibits a strong tendency towards the synchronization, even though the logistic parameter r is in chaotic region (r
=
4.0,3.82 and 3.6).


## 3. How to Face the Complexity of Climate Models

### 3.1. Background

In the introduction we have considered the complex ocean/atmosphere/land dynamical system, called weather and its long time average climate, as a complex one. This system is modeled by climate models having different levels of sophistication. An important concept in climate system modeling is that of a spectrum of models of differing levels of complexity, each being optimum for answering specific questions. It is not meaningful to judge one level as being better or worse than another independently of the context of analysis. What is important is that each model be asked questions appropriate for its level of complexity and quality of its simulation [60]. In this paper they comprehensively considered the following: (i) Earth system models of intermediate complexity, that is, reduced-resolution models that incorporate most of the processes represented by AOGCMs (Atmosphere-Ocean General Circulation Models) and models of reduced complexity. However, in this paper the model complexity will be analyzed on the basis of calculation of the maximum complexity that can be generated by a model. Given a time series and the problem of choosing among a number of climate models to study it, we suggest that models whose maximum complexity is lower than the time series complexity should be disregarded because of being unable to reconstruct some of the structures contained in the data. The increasing complexity of those models is a growing concern in the modeling community. They are used to integrate and process knowledge from different parts of the system and in doing so allow us to test system understanding and create hypotheses about how the system will respond to the virtual numerical experiments. However, if we strive to design our models to be more “realistic,” we have to include more and more parameters and processes. Then, within this approach the model complexity increases, and thus we are less able to manage and understand the model behavior. Obviously, the question about model complexity could be considered from the standpoint of a practitioner who sees it as a compromise between complexity and manageability. His\her question is basically very simple: “How can I check if this model is appropriate to study this problem with this data set?” According to Boschetti (2008): “As a result, the ability of a model to simulate complex dynamics is no more an absolute value in itself, rather a relative one: we need enough complexity to realistically model a process, but not so much that we ourselves cannot handle” [61].Clearly, an answer to the above question requires (i) a definition and a measure of complexity and (ii) that this measure is equally applicable to the model and to the data, because some sort of comparison is necessary. It is a hard task to find that measure even approximately. However, intuitively we can put a cadence on a view of complexity which is more related to a model’s dynamical properties rather than its architecture. Thus, we can say that, in developing tools, an advantage will be given to a tool which gives answers to the following questions: (i) what is the maximal dynamical complexity a given model can generate? and (ii) what kind of different dynamical behaviors can a given model generate? as it is underlined by Boschetti (2008). For our consideration we will rely on Boschetti (2008) who defined the complexity of an ecological model as the statistical complexity of the output it produces that allows a direct comparison between data and model complexity [61]. Among the many different measures of complexity available in the literature, for that purpose, he adopted the statistical complexity defined in [62].


### 3.2. An Example of Comparison between Complexities of a Global and Regional Model

In this subsection we will illustrate an example of comparison between complexities of global and regional model. Here, we do not deal with statistical complexity of the global and regional models. Our intention is just to show possible differences in complexities of time series of precipitation as well as air temperature for both models, applying the algorithm for calculating the Kolmogorov complexity.We have calculated the Kolmogorov complexity following Lempel and Ziv [63] who developed an algorithm for calculating the measure of complexity. It can be considered as a measure of the degree of disorder or irregularity in a time series. This algorithm performs the Kolmogorov complexity analysis of a time series {
x
i
}, i
=
1,2
,
3,4
,
…
,
N, in the following way.Step 1.
Encode the time series by constructing a sequenceS of the characters 0 and 1 written as {
s
(
i
)
}
,
i
=
1,2
,
3,4
,
…
,
N, according to the following rule:
(20)
s
(
i
)
=
{
0
,
x
i
<
x
*
,
1
,
x
i
≥
x
*
.
Here x
* is a chosen threshold. We use the mean value of the time series to be the threshold. The mean value of the time series has often been used as the threshold [64]. Depending on the application, other encoding schemes are also used.Step 2.
Calculate the complexity counterc
(
N
). The c
(
N
) is defined as the minimum number of distinct patterns contained in a given character sequence. The complexity counter c
(
N
) is a function of the length of the sequence N. The value of c
(
N
) is approaching an ultimate value b
(
N
) as N approaches infinity; that is,
(21)
c
(
N
)
=
O
(
b
(
N
)
)
,
b
(
N
)
=
N
lo
g
2
N
.Step 3.
Calculate the normalized complexity measureC
k
(
N
), which is defined as
(22)
C
k
(
N
)
=
c
(
N
)
b
(
N
)
=
c
(
N
)
lo
g
2
N
N
.
The C
k
(
N
) is a parameter to represent the information quantity contained in a time series, and it is to be 0 for a periodic or regular time series and to be 1 for a random time series, if N is large enough. For a nonlinear time series, C
k
(
N
) is to be between 0 and 1.In order to calculate complexities of model time series we have used (i) air temperature and (ii) precipitation time series which are outputs from climate simulations for Belgrade and Novi Sad in Serbia [65, 66]. The Belgrade data set, for the period 2071–2100, was derived from (a) the SINTEX-G which is a coupled Atmosphere-Ocean General Circulation Model [67] and (b) the Eta-POM regional model [65]. The Novi Sad data set, for the period 2020–2050, was derived from (a) the ECHAM5 which is the 5th generation of the ECHAM5 general circulation model [68] and (b) the RegCM regional model [69].We have calculated the Kolmogorov complexity for each time series obtained when each sample, in the original time series, is used as a threshold (N
=
10800 for Belgrade and N
=
11323 for Novi Sad). The results are depicted in Figure 6. We also have calculated Kolmogorov complexity (KL) and its maximal value (KLM) of time series from Figure 6. Results of those calculations are given in Table 1.Table 1
Kolmogorov complexity (KL and its maximum: KLM) values for the precipitation and air temperature time series for Belgrade and Novi Sad, in Serbia, obtained from climate simulations using different models.
Quantity
Measure
Model
Global
Regional
SINTEX-G
ECHAM5
Eta-POM
RegCM
Temperature
KL
0.176
0.207
(Belgrade)
KLM
0.326
0.331
Temperature
KL
0.241
0.251
(Novi Sad)
KLM
0.318
0.354
Precipitation
KL
0.705
0.671
(Belgrade)
KLM
0.834
0.793
Precipitation
KL
0.265
0.871
(Novi Sad)
KLM
0.289
0.935Kolmogorov complexity for the (a) precipitation and (b) air temperature time series for Belgrade and (c) precipitation and (d) temperature for Novi Sad, in Serbia, obtained from climate simulations using different models. Onx-axis are depicted the values of the time series normalized as x
i
=
(
X
i
-
X
min
⁡
)
/
(
X
max
⁡
-
X
min
⁡
), where {
X
i
} is the time series of the precipitation or air temperature, X
max
⁡
=
max
⁡
{
X
i
}, and X
min
⁡
=
min
⁡
{
X
i
}.
(a)
(b)
(c)
(d)From Figure6(a) it is seen that there is no difference between complexities of the precipitation time series for Belgrade obtained by both models (global SINTEX-G and regional Eta-POM) over all amplitudes in time series. Moreover, the SINTEX-G model has slightly higher complexity. In contrast to that, Figure 6(b) depicts that the Eta-POM model mostly has the higher complexity than the SINTEX-G one for the air temperature time series. From Table 1 we can see that for air temperature time series the KL for the Eta-POM model (0.207) is higher than for the SINTEX-G model (0.176), while the KLM values are practically the same (0.331 and 0.326). Note that all of these complexities are pronouncedly low. Further inspection of this table clearly shows that the precipitation time series obtained by the SINTEX-G model has higher complexities (KL: 0.705 and KLM: 0.834) than those obtained by the Eta-POM model (KL: 0.671 and KLM: 0.793). This analysis indicates that the SINTEX-G and Eta-POM models, in particular for precipitation, have approximately the same level of complexity.Now, we analyze the air temperature and precipitation time series for Novi Sad obtained by the global ECHAM5 and regional RegCM models. From Figure6(c) it is seen that there is a large difference between complexities of the precipitation time series over all amplitudes in time series. Moreover, the RegCM model has pronouncedly higher complexity. Figure 6(d) depicts that the RegCM and ECHAM5 models mostly have very similar complexities for the air temperature time series. From Table 1 we can see that for air temperature time series the KL for the RegCM model (0.251) is higher than for the ECHAM5 model (0.241) and also for the KLM values, 0.354 and 0.318, respectively. Similarly, as for the above analyzed models, these values of complexity are still low. Further inspection of this table clearly shows that the precipitation time series obtained by the ECHAM5 model has lower complexities (KL; 0.265 and KLM: 0.289) than those obtained by the RegCM model (KL: 0.871 and KLM: 0.935). This analysis indicates the ECHAM5 and RegCM models have approximately the same level of complexity in simulation of the air temperature. In contrast to that, there is a large difference in capabilities of these models to simulate the participation. Note that a higher value of the KL points out the presence of stochastic influence of different factors on a time series. In this paper we suggest that climate models whose maximum complexity is lower than the time series complexity should be disregarded because of being unable to reconstruct some of the structures contained in the data. To our knowledge this complexity analysis has not been used for analyzing the complexity of climate models. However, for more reliable conclusion that could be given we need to test outputs of many different GCM and RegCM models.


## 3.1. Background

In the introduction we have considered the complex ocean/atmosphere/land dynamical system, called weather and its long time average climate, as a complex one. This system is modeled by climate models having different levels of sophistication. An important concept in climate system modeling is that of a spectrum of models of differing levels of complexity, each being optimum for answering specific questions. It is not meaningful to judge one level as being better or worse than another independently of the context of analysis. What is important is that each model be asked questions appropriate for its level of complexity and quality of its simulation [60]. In this paper they comprehensively considered the following: (i) Earth system models of intermediate complexity, that is, reduced-resolution models that incorporate most of the processes represented by AOGCMs (Atmosphere-Ocean General Circulation Models) and models of reduced complexity. However, in this paper the model complexity will be analyzed on the basis of calculation of the maximum complexity that can be generated by a model. Given a time series and the problem of choosing among a number of climate models to study it, we suggest that models whose maximum complexity is lower than the time series complexity should be disregarded because of being unable to reconstruct some of the structures contained in the data. The increasing complexity of those models is a growing concern in the modeling community. They are used to integrate and process knowledge from different parts of the system and in doing so allow us to test system understanding and create hypotheses about how the system will respond to the virtual numerical experiments. However, if we strive to design our models to be more “realistic,” we have to include more and more parameters and processes. Then, within this approach the model complexity increases, and thus we are less able to manage and understand the model behavior. Obviously, the question about model complexity could be considered from the standpoint of a practitioner who sees it as a compromise between complexity and manageability. His\her question is basically very simple: “How can I check if this model is appropriate to study this problem with this data set?” According to Boschetti (2008): “As a result, the ability of a model to simulate complex dynamics is no more an absolute value in itself, rather a relative one: we need enough complexity to realistically model a process, but not so much that we ourselves cannot handle” [61].Clearly, an answer to the above question requires (i) a definition and a measure of complexity and (ii) that this measure is equally applicable to the model and to the data, because some sort of comparison is necessary. It is a hard task to find that measure even approximately. However, intuitively we can put a cadence on a view of complexity which is more related to a model’s dynamical properties rather than its architecture. Thus, we can say that, in developing tools, an advantage will be given to a tool which gives answers to the following questions: (i) what is the maximal dynamical complexity a given model can generate? and (ii) what kind of different dynamical behaviors can a given model generate? as it is underlined by Boschetti (2008). For our consideration we will rely on Boschetti (2008) who defined the complexity of an ecological model as the statistical complexity of the output it produces that allows a direct comparison between data and model complexity [61]. Among the many different measures of complexity available in the literature, for that purpose, he adopted the statistical complexity defined in [62].


## 3.2. An Example of Comparison between Complexities of a Global and Regional Model

In this subsection we will illustrate an example of comparison between complexities of global and regional model. Here, we do not deal with statistical complexity of the global and regional models. Our intention is just to show possible differences in complexities of time series of precipitation as well as air temperature for both models, applying the algorithm for calculating the Kolmogorov complexity.We have calculated the Kolmogorov complexity following Lempel and Ziv [63] who developed an algorithm for calculating the measure of complexity. It can be considered as a measure of the degree of disorder or irregularity in a time series. This algorithm performs the Kolmogorov complexity analysis of a time series {
x
i
}, i
=
1,2
,
3,4
,
…
,
N, in the following way.Step 1.
Encode the time series by constructing a sequenceS of the characters 0 and 1 written as {
s
(
i
)
}
,
i
=
1,2
,
3,4
,
…
,
N, according to the following rule:
(20)
s
(
i
)
=
{
0
,
x
i
<
x
*
,
1
,
x
i
≥
x
*
.
Here x
* is a chosen threshold. We use the mean value of the time series to be the threshold. The mean value of the time series has often been used as the threshold [64]. Depending on the application, other encoding schemes are also used.Step 2.
Calculate the complexity counterc
(
N
). The c
(
N
) is defined as the minimum number of distinct patterns contained in a given character sequence. The complexity counter c
(
N
) is a function of the length of the sequence N. The value of c
(
N
) is approaching an ultimate value b
(
N
) as N approaches infinity; that is,
(21)
c
(
N
)
=
O
(
b
(
N
)
)
,
b
(
N
)
=
N
lo
g
2
N
.Step 3.
Calculate the normalized complexity measureC
k
(
N
), which is defined as
(22)
C
k
(
N
)
=
c
(
N
)
b
(
N
)
=
c
(
N
)
lo
g
2
N
N
.
The C
k
(
N
) is a parameter to represent the information quantity contained in a time series, and it is to be 0 for a periodic or regular time series and to be 1 for a random time series, if N is large enough. For a nonlinear time series, C
k
(
N
) is to be between 0 and 1.In order to calculate complexities of model time series we have used (i) air temperature and (ii) precipitation time series which are outputs from climate simulations for Belgrade and Novi Sad in Serbia [65, 66]. The Belgrade data set, for the period 2071–2100, was derived from (a) the SINTEX-G which is a coupled Atmosphere-Ocean General Circulation Model [67] and (b) the Eta-POM regional model [65]. The Novi Sad data set, for the period 2020–2050, was derived from (a) the ECHAM5 which is the 5th generation of the ECHAM5 general circulation model [68] and (b) the RegCM regional model [69].We have calculated the Kolmogorov complexity for each time series obtained when each sample, in the original time series, is used as a threshold (N
=
10800 for Belgrade and N
=
11323 for Novi Sad). The results are depicted in Figure 6. We also have calculated Kolmogorov complexity (KL) and its maximal value (KLM) of time series from Figure 6. Results of those calculations are given in Table 1.Table 1
Kolmogorov complexity (KL and its maximum: KLM) values for the precipitation and air temperature time series for Belgrade and Novi Sad, in Serbia, obtained from climate simulations using different models.
Quantity
Measure
Model
Global
Regional
SINTEX-G
ECHAM5
Eta-POM
RegCM
Temperature
KL
0.176
0.207
(Belgrade)
KLM
0.326
0.331
Temperature
KL
0.241
0.251
(Novi Sad)
KLM
0.318
0.354
Precipitation
KL
0.705
0.671
(Belgrade)
KLM
0.834
0.793
Precipitation
KL
0.265
0.871
(Novi Sad)
KLM
0.289
0.935Kolmogorov complexity for the (a) precipitation and (b) air temperature time series for Belgrade and (c) precipitation and (d) temperature for Novi Sad, in Serbia, obtained from climate simulations using different models. Onx-axis are depicted the values of the time series normalized as x
i
=
(
X
i
-
X
min
⁡
)
/
(
X
max
⁡
-
X
min
⁡
), where {
X
i
} is the time series of the precipitation or air temperature, X
max
⁡
=
max
⁡
{
X
i
}, and X
min
⁡
=
min
⁡
{
X
i
}.
(a)
(b)
(c)
(d)From Figure6(a) it is seen that there is no difference between complexities of the precipitation time series for Belgrade obtained by both models (global SINTEX-G and regional Eta-POM) over all amplitudes in time series. Moreover, the SINTEX-G model has slightly higher complexity. In contrast to that, Figure 6(b) depicts that the Eta-POM model mostly has the higher complexity than the SINTEX-G one for the air temperature time series. From Table 1 we can see that for air temperature time series the KL for the Eta-POM model (0.207) is higher than for the SINTEX-G model (0.176), while the KLM values are practically the same (0.331 and 0.326). Note that all of these complexities are pronouncedly low. Further inspection of this table clearly shows that the precipitation time series obtained by the SINTEX-G model has higher complexities (KL: 0.705 and KLM: 0.834) than those obtained by the Eta-POM model (KL: 0.671 and KLM: 0.793). This analysis indicates that the SINTEX-G and Eta-POM models, in particular for precipitation, have approximately the same level of complexity.Now, we analyze the air temperature and precipitation time series for Novi Sad obtained by the global ECHAM5 and regional RegCM models. From Figure6(c) it is seen that there is a large difference between complexities of the precipitation time series over all amplitudes in time series. Moreover, the RegCM model has pronouncedly higher complexity. Figure 6(d) depicts that the RegCM and ECHAM5 models mostly have very similar complexities for the air temperature time series. From Table 1 we can see that for air temperature time series the KL for the RegCM model (0.251) is higher than for the ECHAM5 model (0.241) and also for the KLM values, 0.354 and 0.318, respectively. Similarly, as for the above analyzed models, these values of complexity are still low. Further inspection of this table clearly shows that the precipitation time series obtained by the ECHAM5 model has lower complexities (KL; 0.265 and KLM: 0.289) than those obtained by the RegCM model (KL: 0.871 and KLM: 0.935). This analysis indicates the ECHAM5 and RegCM models have approximately the same level of complexity in simulation of the air temperature. In contrast to that, there is a large difference in capabilities of these models to simulate the participation. Note that a higher value of the KL points out the presence of stochastic influence of different factors on a time series. In this paper we suggest that climate models whose maximum complexity is lower than the time series complexity should be disregarded because of being unable to reconstruct some of the structures contained in the data. To our knowledge this complexity analysis has not been used for analyzing the complexity of climate models. However, for more reliable conclusion that could be given we need to test outputs of many different GCM and RegCM models.


## 4. Concluding Remarks

We have considered climate predictions through two issues: (i) occurrence of chaos and (ii) complexity in climate models. We have given a detailed overview of literature related to this subject. Then, we considered the climate modeling through the light of Gödel’s theorem that says that Number Theory is morecomplex than any of its formalizations; further we have underlined Rosen’s definition of complexity and predictability. In that sense the following points can be enhanced.First, we have pointed out occurrence of chaos in computing the environmental interface temperature from the energy balance equation when the given differential equation is replaced by a difference equation. For that purpose we have analyzed a coupled system of equations, often used in climate models. It is shown that the Lyapunov exponent mostly has positive values approving presence of chaos in this system, but there are still some strait regions where the Lyapunov exponent is negative, that is, where there exist physically meaningful solutions. This analysis, set into context of the climate modeling, points out the fact that there exists set of domains where the environmental interface temperature cannot be calculated by the physics of currently designed climate models.Second, we have analyzed the coupling of processes of vertical and horizontal energy exchange between environmental interfaces which is described by the dynamics of driven coupled oscillators. To study that coupling, when a perturbation is introduced in the system, as a function of the coupling parameter, the logistic parameter, and the parameter of exchange, we have considered dynamics of two maps serving the diffusive coupling. Then, we have performed simulations, calculating the Lyapunov exponent, with the closed contour ofN
=
100 environmental interfaces.Finally, we have explored possible differences in complexities of two global and two regional climate models using their output time series for the precipitation and air temperature. We have applied the algorithm for calculating the Kolmogorov complexity on those time series. We have found differences in the level of complexity among models. However, for more reliable conclusion that could be given we need to test outputs of many different GCM and RegCM models.



---

*Source: 878249-2014-04-03.xml*