Climate predictions: the chaos and complexity in climate models

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 Godel's Theorem and Rosen's definition of complexity and predictability is discussed. It is pointed out to occurrence of chaos in computing the environmental interface temperature from the energy balance equation given in a difference form. A coupled system of equations, often used in climate models is analyzed. It is shown that the Lyapunov exponent mostly has positive values allowing presence of chaos in this systems. The horizontal energy exchange between environmental interfaces, which is described by the dynamics of driven coupled oscillators, is analyzed. Their behavior and synchronization, when a perturbation is introduced in the system, as a function of the coupling parameters, the logistic parameter and the parameter of exchange, was studied calculating the Lyapunov exponent under simulations with the closed contour of N=100 environmental interfaces. Finally, we have explored possible differences in complexities of two global and two regional climate models using their output time series by applying the algorithm for calculating the Kolmogorov complexity.


Introduction
Among the most interesting and fascinating phenomena that are predicted/predictable are the chaotic ocean/atmosphere/land system called weather and its long time averageclimate. 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 early sixties of the 20 th century, which was elaborated in pioneering works by Edward N. Lorenz [2][3][4][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 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 a 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 it is numerically considered the predictability of a forced nonlinear system, proposed by Lorenz, as a compelling heuristic model of the mid-latitude global circulation.
The above insight of the predictability underlined in the context of the "environmental predictability" (primarily linked to the climate change issues), we finish with the question: Can we significantly "improve" the weather/climate predictions comparing to the level they currently reached? The answer can not be strictly elaborated with either yes or no. An optimistic and acceptable attitude, that prefers option yes, is concisely written down by Hunt (1999) as the 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 option no 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 the 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 word 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 both systems of entailment. It is the relation 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 the modeling relation. In a precise sense, this theorem asserts that a formalization which all entailment is syntactic entailment is too impoverished in entailment to be congruent to Number Theory, no matter how we try to establish such congruence. This kind of situation is termed complexity by Rosen (1977) [35]. Namely, in this light, Gödel's Theorem says that Number Theory is more complex than any of its formalization, or equivalently, that formalizations, governed by syntactic inference alone, are simpler 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 selforganizing systems; from hierarchical steering to participation,  from conditions of equilibrium to dynamic balances of non equilibrium; 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 emergetism and from being to becoming.
In this paper we address three issues that, to our mind, are important for further improvements in designing the climate models. (1) The phenomenon of chaos in computing the environmental interface temperature from the energy balance equation which will be considered through the question how to replace given differential equations by appropriate difference equations in climate simulations? (Section 2). (2) The synchronization of energy exchange between environmental interfaces in dependence on perturbation of environmental parameters (Section 3) and (3) complexity analysis of the climate model output time series which is elaborated in Section 4. In Section 5 we give concluding remarks.

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 [36]. Traditionally, this approach has worked well in physics and in engeneering 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, i.e. environmental sciences in their broadest context [37], 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 modelling and thus in climate simulations [38].
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 [38]. 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 above-mentioned approximation comes in. Under this kind of circumstances, it seems more natural to 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 [36].
In this section 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 as interface 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 [39]. 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 ( Figure  1) [40]. 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 [41]. 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, i.e. the deeper soil layer temperature was constant, but it can also vary in time making with the energy balance equation a coupled system of equations. 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 section 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 occurence 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 and then we analyze behavior of the solutions of the coupled system and we have examined domains of stability using the Lyapunov exponent.

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 [40]. 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 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 [42], i.e. ( ) 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 where C L is the water vapour transfer coefficient, o -1 b = 0.06337 C , d is parameter which occurs in expanding the series [43]. Further, the conduction of the heat into the soil can be written in the form 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 where C H is the sensible heat transfer coefficient. The prognostic equation for temperature of the deeper soil layer T d is where τ =86400 s. After collecting all terms (2)-(6), the coupled system takes the form More details about the nature and the range of physical parameters C R , C L , C D and C H can be found in [44]. Now, using the time scheme forward in time (n indicates the time step) and dividing both sides of Eqs. (7) and (8) with the constant temperature 0 T (for example, value of mean Earth temperature, i.e. 0 T = 288K ) we get Finally, introducing replacements 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) [0, 4] A∈ 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 with r . All other groups of parameters in the system (13)- (14) 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 (13)- (14) can be written in the form of coupled maps, i.e., Now we analyse the stability of physical solutions of coupled maps (15)-(16), 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 is negative, although for chaotic attractor is positive. Calculating Lyapunov exponent for the coupled system (15) , because we thought that 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 [45]. 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, i.e. domains of stability.

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 in designing the model, the complexity, and thus its reliability, will not be retained in the highest degree. In Section 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, relying on Section 2 and using paper by Mihailović et al. (2012) we analyze the horizontal energy exchange between environmental interfaces which is described by the dynamics of driven coupled oscillators [46]. 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 [46].
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 ( ) the dimensionless temperature of environmental interface and r is a logistic parameter [45,46]. 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 [36,38] .
The dynamics of this map (Eq. (17)) are governed by two parameters, p and r , which 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 exchange p 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 represent the 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 [47]. 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 [47].

Observational heterarchy and horizontal energy exchange between environmental interfaces
Observational heterarchy consists of two sets of intra-layer maps, called Intent and Extent perspectives, and inter-layer operations satisfying the following conditions: (1) the inter-layer operations inherit the mixture of intra-and inter-layer operations and (2) there is a procedure by which the inter-layer operation can be regarded as an adjoint functor. If the inter-layer operation satisfies the conditions (1) -(2), it is called a pre-functor [47]. Preserving the above mentioned composition occurs as follows: A pre-functor, : F 〈 〉 Int → Ext is mapping a set, X , to a set, F X 〈 〉 , and map Φ to a map, f f In this sense we call f * pseudo-inverse 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, 1 f f * − = holds, then a pre-functor can become a functor. Given two maps, , : It implies that F preserves the composition of maps, Φ and Ψ . The time development of the environmental surface dynamics , i n x , for two interfaces, is expressed as where: n is the time iteration, , i n x ∈ , c the coupling parameter, f the map representing the horizontal energy exchange between environmental interfaces, Φ is one of maps in the pair ( , ) Ψ Φ whose composition is preserved by a pre-functor 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 then it can be reduced to Eq. (19). We perform our analysis following the procedure described in [47]. First, in this section we address the synchronization of the passive coupling for two environmental interfaces given by Eqs. (19) and (17), and then, in the next section, we will show that perturbation can modify the dynamics and enhance robust behavior in a multi-environmental interface system of active coupling. Synchronization is well-known collective phenomenon in various multi-component physical as well as the climate systems [48][49][50]. 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 [51].We studied the stability of the fixed point by linearzing 2 n ≥ component coupled system, and obtain

3c Simulations of active coupling in a multi-environmental interface system
Here, we address the behavior of active coupling [47], and estimate whether a coupled map system described above can achieve synchronization under influence of perturbations. The dynamics of two-environmental interface system called active coupling [47], used for simulations, is expressed as .
We note, that the dynamical system defined by Eqs. (23a) and (23d) 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 adjunction or the equivalence between Intent and Extent. Compared with passive coupling, the behavior of active coupling is much more complex [47]. In Eqs. (23a) -(23d), because of a pseudo-inverse map, f * , all calculations are defined to be approximations. In simulations, the Intent map was a discontinuous map, expressed by In order to see how perturbation enhances robust behavior in the framework of observational heterarchy in a multi-environmental interface system represented by closed contour of coupled environmental interfaces exchanging the energy horizontally. Then the system of coupled difference equations for N environmental interfaces exchanging the energy, can be written in the form of matrix equation The elements in matrices in Eq. (24a) are Simulations with the active coupling, defined by Eqs. (23a)-(23d), were performed with and without perturbation given as in [47]. 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 p and the logistic parameter r . Simulations were performed with the closed contour of 100 N = interfaces. The Lyapunov exponent was calculated using Eqs. (20) - (21) and the Jacobian of the system given by Eqs. (24b)-(24c) is representing this contour.
In calculating λ , for each c from 0.0 to 1.0 with step 0.005, 4 10 iterations were applied for an initial state, and then first 3 10 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 4 0 r .
= and the fixed value of the parameter of exchange p (Figures 5a -5c); 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 5d -5f). Figures 5a -5c depict that in the chaotic regime ( 4 0 r . = ), regardless of the value p, the Lyapunov exponent is always positive ( 0 λ > ) and therefore the process of the horizontal energy exchange in a multi-environmental interface system is always unsynchronized. However, the stormy perturbation disturbs this state (Figures 5a -5c).
Although the logistic parameter is settled at 4.0 r = for chaotic behavior, the coupling parameter c tunes interaction and leads to synchronization in some intervals, particularly for . This behavior is more pronounced in Figures 5d -5f where p is randomly chosen; here the process of horizontal energy exchange in a multi-environmental 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).

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. 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, 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 Boshetti (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" [52].
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, which an advantage will be given to a tool which gives answers on the 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 Boshetti (2008). For our consideration we will rely on Boshetti (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 [52]. Among the many different measures of complexity available in the literature, for that purpose, he adopted the statistical complexity defined in [53].

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 [54] 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 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 [55]. Depending on the application, other encoding schemes are also used. Step The ) (N C k is a parameter to represent the information quantity contained in a time series, and it is to be a 0 for a periodic or regular time series and to be a 1 for a random time series, if N is large enough. For a non-linear time series, ) (N C k 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 [56,57]. 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 [58] and (b) Eta-POM regional model [56].The Novi Sad data set, for the period 2020-2050, was derived from: (a) the ECHAM5 which is the 5th generation of the ECHAM general circulation model [59] and (b) RegCM regional model [60].
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. From Figure 6a 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 6b 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 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 Figure 6c 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 6d 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 these models to simulate the participation. To our knowledge this complexity analysis has not been used for analyzing the complexity of climate models.

Concluding remarks
We have considered some issues which are relevant for climate modeling. We gave 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 more complex than any of its formalization; further we clearly underlined the Rosen's definition of complexity and predictability. We have emphasized three issues.
Firstly, we have pointed out on 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 and where the solutions of the coupled system are stable.
Secondly, we have analyzed the horizontal energy exchange between environmental interfaces which is described by the dynamics of driven coupled oscillators. To study 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, we have considered dynamics of two maps serving the diffusive coupling. Then, we have performed simulations, calculating the Lyapunov exponent, with the closed contour of N = 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.