Using Machine Learning for Climate Modelling: Application of Neural Networks to a Slow-Fast Chaotic Dynamical System as a Case Study

Abstract

This paper explores the capabilities of two types of recurrent neural networks, unidirectional and bidirectional long short-term memory networks, to build a surrogate model for a coupled fast–slow dynamic system and predicting its nonlinear chaotic behaviour. The dynamical system in question, comprising two versions of the classical Lorenz model with a small time-scale separation factor, is treated as an atmosphere–ocean research simulator. In numerical experiments, the number of hidden layers and the number of nodes in each hidden layer varied from 1 to 5 and from 16 to 256, respectively. The basic configuration of the surrogate model, determined experimentally, has three hidden layers, each comprising between 16 and 128 nodes. The findings revealed the advantages of bidirectional neural networks over unidirectional ones in terms of forecasting accuracy. As the forecast horizon increases, the accuracy of forecasts deteriorates, which was quite expected, primarily due to the chaotic behaviour of the fast subsystem. All other things being equal, increasing the number of neurons in hidden layers facilitates the improvement of forecast accuracy. The obtained results indicate that the quality of short-term forecasts with a lead time of up to 0.75 model time units (MTU) improves most significantly. The predictability limit of the fast subsystem (“atmosphere”) is somewhat greater than the Lyapunov time.

1. Introduction

Predicting climate change caused by natural and human-made factors is one of the key problems of contemporary climatology [1,2]. Since the Earth’s climate system (ECS) has a number of rather specific features [3], to predict and project future climate, researchers use fairly complex climate models, which describe, among other things, the general circulation of the atmosphere and oceans [4,5,6]. In general terms, climate models represent systems of partial differential equations (DEs), which are mathematical expressions of the laws of physics, chemistry, thermodynamics, and fluid mechanics. Since the solution to these equations cannot be found analytically, numerical methods are usually used for this purpose. Building and running climate models require enormous human and computational resources [7]. Researchers have always looked for ways to overcome these critical issues. Recent advancements in artificial intelligence (AI), and machine learning (ML) in particular, have allowed scientists to create so-called surrogate models (emulators or meta-models), which are comparable in accuracy to traditional “physically based” climate models but are computationally much more efficient (e.g., [8,9,10]). Surrogate models represent quite simplistic versions of complex models of physical and other systems. Typically, these models map input data to output data, assuming that the inputs and outputs are (nonlinearly) related to each other. Supervised ML is a suitable and effective tool for building emulators.

Over the last few years, there has been a steady increase in the number of publications considering various aspects of the application of ML—and in particular deep learning (DL), which is a subset of ML—to solve different problems in the fields of meteorology, hydrology, and climatology. This demonstrates the growing interest of climate modelling experts in AI. The main areas of AI application in climate research and numerical weather prediction (NWP) have been discussed in a number of review articles (e.g., [11,12,13,14,15,16,17,18,19,20,21,22]). Experts agree that the application of AI technologies, and ML algorithms in particular, have great potential for bias correction of contact and remote observations, data assimilation, developing parameterization schemes of sub-grid physical processes, processing satellite images, estimating NWP and climate model parameters, post-processing of NWP and climate model outputs, building surrogate weather prediction and climate models, etc.

As mentioned above, all existing physically based (classical) NWP and climate models are built on a unified methodological basis, since DEs used in these models to describe the evolution of the ECS and its components (e.g., the atmosphere and oceans) are essentially the same. Meanwhile, in recent times, in the study of complex systems, an interdisciplinary approach called surrogate modelling has become increasingly widespread (e.g., [23,24,25,26,27,28]. The essence of this approach is to build an emulator that could replace a physically based model of a system, with a comparable degree of adequacy. The process of building an emulator is based on ML, including the use of multilayer artificial neural networks (ANN). Meta-models vastly (by orders of magnitude) reduce the computational time. In other words, they have significantly greater computational efficiency compared to physically based models.

To date, the most notable research advances have been achieved in developing ML-based sub-grid-scale parameterization of physics, data pre-processing, quality control and assimilation, and post-processing of modelling and forecasting results (e.g., [8,9,10,19,20,21,22,29] and references therein). In recent years, a number of studies have demonstrated the capabilities of DL to model global temperature trends (e.g., [30]) and phenomena, such as the El Niño–Southern Oscillation using coupled atmosphere–ocean models (e.g., [31,32,33,34]). However, science is never settled; therefore, the potential for using new ML technologies to build surrogate NWP and climate models, comparable in quality to classical models, is actively explored on an ongoing basis (e.g., [8,9,35,36,37,38,39,40,41,42,43,44]).

It is known that the ECS is completely nonlinear by nature. Its components, such as the atmosphere and oceans, as well as various internal dynamical, physical, and chemical processes and cycles, are also nonlinear. For this reason, building a neural network-based ECS model is a challenging problem that requires, among other things, high computing power and tuning, carried out via a large number of numerical experiments [7].

A fairly large, continuously expanding variety of ANNs have been developed to date. All of them can be conditionally divided into classes concerning types of architecture and topology, methods of training, and the determination of the weights of the ANN, as well as the problems solved using the ANN, etc. (see [45,46,47] and references therein). Meanwhile, ANNs of each class may differ from each other in terms of the number of hidden layers and neurons in them, the activation function, the optimization algorithm, and other details that establish the nature of processing input data in relation to output data. In general, the choice of ANN structure depends on the problem being solved, available data, accuracy requirements, experience, and other factors. In this context, reduced-order climate models can be viewed as testing tools that allow you to determine and evaluate the capabilities of ANNs with different architectures to build reliable data-driven surrogate models of the climate system (e.g., [48]).

It is important to note that some reduced-order models describing the essential properties of the ECS can also reproduce the complex and chaotic behavior of the atmosphere, which is the most unstable and fluctuating component of the climate system [49,50].

This paper aims to explore the capabilities of two types of recurrent neural networks (RNN), namely, long short-term memory (LSTM) and bidirectional LSTM (BiLSTM), to mimic the nonlinear chaotic behaviour of a coupled fast–slow dynamic system, which can be viewed as an atmosphere–ocean research simulator [51,52,53,54,55,56]. In doing so, we also demonstrate a methodological approach to constructing an emulator that replaces the original, physically based model of the dynamical system, reproducing processes with two different time scales.

2. Materials and Methods

2.1. Coupled Nonlinear Dynamical System with Two Time Scales

In this paper, a low-order coupled fast–slow nonlinear dynamical system is considered as a physically based representation of the atmosphere–ocean system, for which a corresponding surrogate is constructed, in order to estimate the ability of neural networks to predict nonlinear and chaotic dynamics of complex systems such as the ocean and atmosphere. The dynamical system used combines two versions of the well-known Lorenz model [49] (Lorenz’63) with the small time-scale separation factor τ. Using lowercase letters x, y, and z for the fast subsystem variables and uppercase letters X, Y, and Z for the slow subsystem variables, we write the basic set of model equations as follows [51,52,53,54,55,56]:

$$dx/dt=\sigma \left(y-x\right)-c\left(aX+k\right),$$

(1)

$$dy/dt=rx-y-xz+c\left(Y+k\right),$$

(2)

$$dz/dt=xy-bz-c_{z}Z,$$

(3)

$$dX/dt=\tau \sigma \left(Y-X\right)-c\left(ax+k\right),$$

(4)

$$dY/dt=\tau \left(rX-Y-aXZ\right)+c\left(y+k\right),$$

(5)

$$dZ/dt=\tau \left(aXY-bZ\right)-c_{z}z,$$

(6)

where $\sigma >0$, $b>0$, and $r>0$ are the Lorenz model parameters (recalling that their standard values, providing chaotic dynamics, are $\sigma =10$, $b=8/3$, and $r=28$); $c$ is the coupling strength parameter of the x, y, X, and Y variables while $c_z$ is the coupling strength parameter of the z and Z variables; $k$ is an offset (“uncentering”) parameter; and a is an amplitude scale factor (if $a=1$ then the fast and slow systems have the same amplitude scale). Parameters $c$ and $c_z$ regulate the interconnection strength between the fast and slow subsystems; the larger the parameters $c$ and $c_z$, the stronger the interdependence between two subsystems.

The dynamic behaviour of the coupled system is controlled by the values of parameters $\sigma , b, r, c, c_z, k,a$, and τ. In calculations, the default values of 10, 8/3, and 28 are taken respectively for parameters $\sigma , b$, and $r$ [49]. With no loss of generality, we can assume that $k=1$, $a=1$, and $c=c_z$. The time-scale separation factor τ is taken to be 0.1, meaning that the fast subsystem is 10 times “faster” than the slow subsystem.

Particular attention should be paid to the coupling strength parameter $c$ due to its decisive role in shaping the dynamic behaviour of the system (1)–(6) [54,56]. It is known that the stability of dynamical systems in their phase space is characterized by Lyapunov exponents, the number of which is equal to the dimensionality of the system’s phase space. The existence of a positive exponent in the Lyapunov spectrum indicates chaotic behaviour of the dynamical system in question. Let us note in passing that phase-space trajectories of climatic systems are unstable in the sense of Lyapunov [57]. Since the dynamical system under consideration is six dimensional, the general patterns of its behaviour are determined by the six-component Lyapunov spectrum $\left\{{\lambda }_1, \dots , {\lambda }_6\right\}$, depending on parameter $c$ [54,56]. Calculations show that for $c\in \left[0.1, 0.9\right]$ the maximal (largest) Lyapunov exponent ${\lambda }_1$ is positive, and hence the dynamical system (1)–(6) exhibits the behaviour known as deterministic chaos. To build the emulator, we chose the value of 0.15 for parameter $c$ [54], for which ${\lambda }_1=8.60\times {10}^{-1}$.

2.2. Building a Surrogate Model of a Nonlinear Dynamical System with Two Time Scales: Problem Statement

Various ML techniques are used to build surrogate models, ranging from relatively simple (e.g., regression) to complex (e.g., neural networks). We are building a surrogate model based on neural networks, applying a supervised ML approach to “tune” the model that we trying to build. Supervised ML requires training and validation datasets. Sometimes a testing dataset is also required. In this paper, the training and validation datasets are generated by numerically integrating the set of Equations (1)–(6), performed using the fourth order Runge–Kutta method with a time step of $\Delta t={10}^{-2}$. At first, the Equations (1)–(6) are integrated from $t_{-T}={-2}^{14}\Delta t$ to $t_0=0$ with a random initial condition. During this time period, the coupled system reaches its attractor. This procedure allows the initial transition period to be excluded from consideration. Then, state vector $x_0=\left(x_0, y_0, z_0, X_0, Y_0, Z_0\right)$, calculated at time $t_0$, is taken as the initial condition when solving Equations (1)–(6)over the given time interval $\left[t_0, t_{end}\right]$, where $t_{end}=0.85\times {10}^4\Delta t=85$ model time units (MTUs). The obtained numerical solution represents dataset $\mathcal{D}$, which will be split into training and validation samples. It should be noted that $\mathcal{D}$ is a multivariate time series consisting of six interrelated variables x, y, z, X, Y, and Z, each of which is an equally spaced (regular) one-dimensional time series. We use 20% of the data for validation and the remaining 80% of the data for training.

Supervised learning is based on the assumption that there is some unknown relationship between the inputs and outputs of the model under test. Only a finite set of precedents (model input–output pairs), called the training set, is known. It is the set of precedents that is used to determine the relationship between the input and output of the model, which, if successfully trained, can then serve predictive purposes [58].

The emulator building problem is as follows: let a training sample $S_n=\left\{\left(x_i,y_i\right), i=1,\dots ,n\right\}$ be given, obtained using model $y=f\left(x\right)$, where $x\in X\subset {\mathbb{R}}^{\mathit{k}}$ is the input vector, $y\in Y\subset {\mathbb{R}}^{\mathit{l}}$ is the output vector, and $\left(x_i,y_i\right)\in X\times Y$ is a precedent (an “input–output” pair). It is required to determine an unknown empirical relationship $\widehat{f}\left(x\right)=\widehat{f}\left(x|S_{\mathit{n}}\right)$ based on the training sample $S_n$, such that $\widehat{f}\left(x\right)=f\left(x\right)+ϵ$ for all $x\in X$, where $ϵ$ is the “random” error term, which can be considered to be small.

The empirical relationship obtained in this way is considered suitable for practical purposes if generalization error $E_{gen}$, between predicted values and actual values of validation sample $S_m=\left\{\left(x_j,y_j\right), j=1,\dots ,m\right\}$

$$E_{gen}={\displaystyle \frac{1}{m}}{{\displaystyle \sum }}_{j=1}^m{\left[y_j-\widehat{f}\left(x_j\right)\right]}^2$$

(7)

is fairly small.

To build the surrogate model, we use LSTM and BiLSTM neural networks, which are types of RNN [59,60]. Some previous studies have examined LSTMs in varying detail, including their advantages and disadvantages (e.g., [61,62,63,64]). Without going into detail, we note that one of the main advantages of LSTMs is their ability to detect nonlinear patterns and long-term dependencies in data. The accuracy of time series forecasts produced using neural networks is influenced by the volume and quality of training data, the number of hidden layers and the number of neurons in them, the number of epochs (iterations), the type of activation function, the forecasting time horizon, etc. The complexity of a neural network and its performance are largely determined by the number of hidden layers $N_h$ and the number of nodes in each layer $N_{h,n}$. Unfortunately, there are no clear theoretically grounded approaches to determining the size and number of hidden layers in a neural network. For this reason, parameters $N_h$ and $N_{h,n}$ are usually determined heuristically, using elements of the “trial and error” method, i.e., experimentally. Some recommendations on this matter can be found, for example, in [65]. Thus, to build the emulator we should specify the number of hidden layers and their sizes, the type of activation function, the optimization algorithm, the forecast lead time, and the number of training cycles (epochs). The surrogate model is built using the TensorFlow platform with its high-level Keras API.

3. Results

As mentioned earlier, the surrogate model is built based on RNN architectures, such as unidirectional LSTM and its extension known as bidirectional LSTM. It should be noted that it is impossible to theoretically prove the advantages of BiLSTM over LSTM since we have no appropriate proof methods. When determining the basic hyperparameters of neural networks used to simulate the nonlinear dynamics of the system (1)–(6), we assume that there is no rigorous theoretical proof for their choice to date. This applies, in particular, to determining the number of hidden layers of neural networks $N_h$ and their sizes $N_{h,n}$. Therefore, the number of hidden layers in numerical experiments varied from 1 to 5, with the number of neurons in them ranging from 16 to 256. Along with the hidden layers, the neural networks contain two attributive layers: the input layer and the output layer. The input layer has six neurons (according to the number of dependent variables x, y, z, X, Y, and Z of the model under consideration), while the output layer has three neurons, because we are interested in the “atmospheric” variables x, y, and z, only, i.e., in a “weather prediction”. The learning accuracy and rate of neural networks depend on the optimizer used to update the weights of the neural network. The Keras library has several built-in optimizers. The calculations were made based on the “Adam” (Adaptive Moment Estimation) optimizer, which is one of the most efficient optimization algorithms for recurrent neural network training. The advantages of this algorithm are the effective handling of steep and low gradients arising during neural network training and the ability to partially overcome the over-training problem.

The sigmoid function is used as an activation function, which is defined by the formula $f\left(\phi \right)=1/\left(1+e^{-\phi }\right)$, where $\phi$ is the input variable. Sigmoid function compresses inputs into a bounded range between 0 and 1. It should be noted that that the activation function type affects the neural network learning rate, since its weights change during the learning process, and the extent of these changes depends on the value of the activation function $f\left(\phi \right)$ and its derivative $f\prime \left(\phi \right)$ at the given point $\phi ={\phi }_{*}$. The neural network training rate also depends on the selected initial values of the weights. Since no universal way of initializing the neural network weights exists, random numbers from the interval $\left[-0.5;0.5\right]$ are taken as the initial values of the weights in the problem under consideration. The number of epochs was determined experimentally based on the analysis of the behaviour of the loss function, represented by the mean squared error (MSE):

$$\mathrm{MSE}={\displaystyle \frac{1}{k}}{{\displaystyle \sum }}_{i=1}^k{\left({\phi }_{A,i}-{\widehat{\phi }}_{A,i}\right)}^2$$

(8)

where ${\phi }_A={\left(x,y,z\right)}^T$ and ${\widehat{\phi }}_A={\left(\widehat{x},\widehat{y},\widehat{z}\right)}^T$ are the vectors of “observed” and predicted values of “atmospheric” variables, respectively; and $k$ is the number of forecasts being estimated.

The forecast time horizon (forecast lead time), i.e., the time period over which predictions of “atmospheric” variables were made, varied from 0.05 to 1.5 MTU. It has been experimentally established that a further increase in the forecast lead time makes no sense, since the quality of forecasts deteriorated sharply. This was quite expected, since we considered the chaotic regime of the dynamic system (1)–(6). It should be noted that the predictability limit of chaotic systems is characterized by the Lyapunov time $T_L$, which is the reciprocal of the largest Lyapunov exponent. In our case, this is $T_L=1/{\lambda }_1=1.15$ MTU.

It should be recalled that to obtain the training and validation samples, Equations (1)–(6) were integrated over 8500 forward steps, which corresponds to 85 MTU. Thus, for each variable $x, y, z, X, Y, \mathrm{and}Z$, the training and validation sample sizes were 6800 and 1700, respectively. Before model training or predictions, data scaling was applied to bring all values into the range [0, 1].

It was found that with the selected hyperparameters, the neural network over-training effect is not observed, and the applied dropout procedure worsens results, so it was decided to abandon its use.

Let us consider some of the results obtained in the present study. Figure 1 shows graphs illustrating the time variation of variables x, y, z, X, Y, and Z obtained as a result of numerical integration of Equations (1)–(6). We will categorize these data as “observational data” used to form the training and verification samples. The corresponding attractors of the “fast” and “slow” subsystems are shown in Figure 2. Model testing (in other words, evaluating forecast accuracy) was performed only for the “atmospheric” variables x, y, and z. Forecast accuracy was estimated using the standard deviation $\mathrm{RMSE}=\sqrt{\mathrm{MSE}}$ and the mean absolute error (MAE):

$$\mathrm{MAE}={\displaystyle \frac{1}{k}}{{\displaystyle \sum }}_{i=1}^k\left|{\phi }_{A,i}-{\widehat{\phi }}_{A,i}\right|.$$

(9)

Mean absolute percentage error (MAPE) is used quite often to estimate prediction quality. However, for values close to zero, MAPE is not suitable, since when the denominator in scaling formula (not shown here) is close to zero, the division result tends to infinity. Therefore, only metrics $\mathrm{RMSE}$ and MAE were used to estimate forecast quality.

It is well known that absolutely accurate weather forecasts are fundamentally impossible, although their quality is continuously improving. All other things being equal, forecast accuracy depends largely on lead time; the shorter the lead time, the higher the accuracy of weather forecasts. Therefore, our attention was focused on how forecast accuracy changes with increasing forecast time horizon.

To determine the basic configuration of neural networks used for forecasting, a large number of experiments were conducted with different numbers of hidden layers and nodes per hidden layers. Neural networks with three hidden layers were chosen as the basic networks following the analysed results of these experiments. In particular, Figure 3 illustrates that the three-layer LSTM and BiLSTM networks appear preferable to the others. This figure summarizes the results of testing LSTMs and BiLSTMs with different numbers of hidden layers.

The accuracy metrics of forecasts made over different time horizons using LSTM and BiLSTM in the basic configuration are given in

Table 1. Results of testing the model built on LSTM with 3 hidden layers, each comprising 16 and 128 nodes.

Forecast Lead Time in MTU	16 Nodes		128 Nodes
	RMSE	MAE	RMSE	MAE
5	0.010	0.011	0.010	0.005
10	0.025	0.014	0.009	0.025
50	0.136	0.070	0.136	0.057
75	0.141	0.075	0.144	0.070
100	0.202	0.137	0.151	0.092
150	0.253	0.192	0.221	0.157

and

Table 2. Results of testing the model built on BiLSTM with 3 hidden layers, each comprising 16 and 128 nodes.

Forecast Lead Time	16 Nodes		128 Nodes
	RMSE	MAE	RMSE	MAE
5	0.006	0.004	0.002	0.002
10	0.005	0.004	0.004	0.002
50	0.051	0.026	0.037	0.012
75	0.110	0.058	0.100	0.041
100	0.118	0.058	0.136	0.059
150	0.237	0.165	0.223	0.126

. The findings from the analysis of these tables revealed the advantages of bidirectional neural networks over unidirectional ones in terms of forecasting accuracy. This result seems to be very important since BiLSTM are mainly used for natural language processing rather than for simulating the chaotic behaviour of dynamic systems. The data in Table 1 and Table 2 also show that forecast accuracy decreases with increasing lead time, which is quite expected, mainly due to the chaotic behaviour of the fast subsystem. Thus, the accuracy of “weather” predictions tends to decrease as the forecast time horizon extends, and forecasts beyond 1.5 MTU lead times have no value due to their low accuracy. Overall, it can be concluded that the predictability limit of the fast system (“atmosphere”) is somewhat greater than the Lyapunov time and is about 1.5 MTU.

All other things being equal, increasing the number of neurons in hidden layers facilitates the improvement of forecast accuracy. The obtained results indicate that the quality of short-term forecasts with a lead time of up to 0.75 MTU improves most significantly.

Figure 4, Figure 5 and Figure 6 provide a clear indication of the correspondence between the “observed” values of dependent variables and the predicted values produced by the LSTM and BiLSTM models for different forecast horizons. Although the forecasted series generally approximate the original ones quite well, the quality of forecasts obtained using BiLSTM models is better than forecasts obtained using LSTM models, which is especially noticeable as the forecast horizon increases. To determine the weights of neural networks, 50–100 epochs are sufficient, since a further increase in their number leads to the over-training of neural networks.

4. Discussion and Conclusions

Developing and using climate models to improve our understanding of the climate system and predict future climate change require enormous computational and intellectual resources. However, a powerful technique, called surrogate modelling, significantly reduces these resource requirements, since these detailed and physically based climate models are approximated using much simpler algorithms trained on simulation data. As we know, neural network-based deep learning technology is an important tool for building surrogate models of various physical systems and processes described by differential equations. In this context, it is of the utmost importance to study the neural networks’ capabilities to predict solutions to sets of differential equations of the form $dx/dt=f\left(x\right)$ based on a sample of the state vector x, i.e., without knowing the function $f\left(x\right)$. It should be highlighted that DL can lead to positive results if the development of surrogate models is based on an understanding of the physics of the processes that actually shape the climate. This fact emphasizes the need to continue and intensify research into all aspects of surrogate modelling of climate processes and phenomena. In this regard, low-order models of dynamic systems that describe the essential features of the object under study, in particular its chaotic behaviour, can serve as a useful research tool. Simplified climate models can be used on their own (at least as a qualitative research tool) to study physical climate processes, and can also be used as an instrument for testing new approaches (for example, machine learning technologies and algorithms, and in particular deep learning) to modelling complex systems and processes.

We have shown the procedure, and in essence the methodology, of building a surrogate model (meta-model) that “replaces” the original, physically based, simplified model of the coupled “ocean–atmosphere” system whose behaviour becomes chaotic under certain conditions. The dynamical system, comprising two versions of the classical Lorenz’63 model with a small time-scale separation factor, is treated as a physically based model for which a surrogate model is built. A brief and intelligible description of the process of building a meta-model allows one to form a holistic understanding of the main and related problems that may arise.

The capabilities of two types of recurrent neural networks, unidirectional and bidirectional LSTMs, to build a surrogate model and predicting its nonlinear chaotic behaviour were explored. The results of numerous numerical experiments have shown that LSTMs can be successfully used to build surrogate models of nonlinear systems and predict solutions to differential equations without knowing their right-hand sides. Based on this, it can be concluded that the ability of neural networks to predict the future behaviour of nonlinear dynamic systems over certain time periods stimulates their use for the development of emulators of the climate system and its individual components, as well as parameterization schemes for use in complex models employed to simulate the global climate system. However, the accuracy of surrogate models is affected not only by the architecture and structure (hyperparameters) of neural networks but also by the availability of sufficient training data and its quality. Climate models project climate using observational data and some other factors as inputs, which are inherently uncertain. It is therefore important to assess how the uncertainty of training data affects the quality of climate modelling and projection. In this regard, exploring the capabilities of deep learning methods to represent climate model outputs in probabilistic form appears to be a very significant problem. However, this will be the subject of further research.