Correlation and Spectral Properties of a Coupled Nonlinear Dynamical System in the Context of Numerical Weather Prediction and Climate Modeling

Complex dynamical processes occurring in the earth’s climate system are strongly nonlinear and exhibit wave-like oscillations within broad time-space spectrum. One way to imitate essential features of such processes is using a coupled nonlinear dynamical system, obtained by coupling two versions of the well-known Lorenz (1963) model with distinct time scales that differ by a certain time-scale factor. This dynamical system is frequently applied for studying various aspects of atmospheric and climate dynamics, as well as for estimating the effectiveness of numerical algorithms and techniques used in numerical weather prediction, data assimilation, and climate simulation. This paper examines basic dynamic, correlative, and spectral properties of this system and quantifies the influence of the coupling strength on power spectrum densities, spectrograms, and autocorrelation functions.


Introduction
Numerical weather prediction (NWP) and climate simulation are among the most interesting and important problems facing the world today.Dynamical processes occurring in the earth's climate system have turbulent nature and exhibit wave-like oscillations within broad time-space spectrum and, therefore, are characterized by strong nonlinearity [1,2].The earth's climate system is a complex system that consists of the atmosphere, ocean, lithosphere, cryosphere, and biosphere subsystems.Each of these components has unique dynamics and physics and is characterised by specific time-space spectrum of motions, characteristic time, internal variability, and other properties [3][4][5][6].
NWP and climate simulation focus on processes with very different time scales and, more importantly, pursue very different objectives [7][8][9].NWP is a typical initial value problem aiming to predict, as precisely as possible, the future state of the atmosphere taking into account current (initial) weather conditions.To date, due to fundamental limits to the predictability of the atmospheric processes, the practical importance of NWP is limited to a time horizon of 7-10 days.Errors in the NWP are primarily generated by inaccuracies in the initial conditions rather than by the imperfections in mathematical models (e.g., [10]).By contrast, climate change and variability simulations focus on much longer time scales (several months or even years) and involve the study of stability of climate model attractors with respect to external forcing.Nevertheless, deterministic mathematical models, based on the same physical principles and fundamental laws of physics, such as conservation of momentum, energy, and mass, are commonly applied to both NWP and climate simulations.Mathematically, these physical laws and principles are generally expressed in terms of a set of nonlinear partial differential equations (PDEs) and, in their discrete form, contain a large number of state variables and parameters.Performing a detailed simulation of the dynamics of the earth's climate system and its components using such models requires significant computational resources for simulations and considerable time for analysis and interpretation of obtained results.Consequently, loworder models have drawn the attention of scientists as very attractive and convenient tools for studying various aspects of dynamics of natural physical processes and phenomena and also for estimating the effectiveness of numerical algorithms and techniques used in computational simulations.
Deterministic nonlinear dynamical models are highly sensitive to initial conditions, a property which precludes exact forecasting of their future states and can lead to chaotic numerical solutions.This means that, over time, under certain conditions, the behaviour of a simulated system begins to resemble a random process, despite the fact that the system is defined by deterministic laws and described by deterministic equations.This phenomenon of deterministic chaos was first uncovered by Lorenz as he observed the sensitive dependence of atmospheric convection model output on initial conditions [13].Since the publication of his seminal theoretical work, a vast number of studies were implemented in order to analyse various aspects of nonlinear dynamics and deterministic chaos, as well as properties of the original Lorenz model (hereinafter referred to as L63 model) and the broader family of Lorenz-like systems (e.g., [14][15][16][17][18][19][20][21][22][23]).The Lorenz system is derived by strong spectral truncation of Saltzman's equations, which describe the Rayleigh-Bénard convection [24], and consists of three autonomous ordinary differential equations (ODEs) for time-dependent variables , , and : with  corresponding to the intensity of the convective motion in terms of the stream function,  to the temperature differences between rising and descending currents, and  to the departure of the vertical temperature gradient from its equilibrium magnitude.The Lorenz system also contains three positive parameters , , and , with  being the Prandtl number,  the normalized Rayleigh number, and  a geometric parameter characterizing length scale of the convective cell.
Despite its simplicity, the L63 model is capable of imitating some essential properties of the general circulation of the atmosphere and ocean [25][26][27][28].For example, the heat flux from equator to the poles can be represented by variable , which is proportional to meridional temperature gradient that can be represented by parameter .Various academic papers have previously presented successful applications of the L63 model to climate studies (e.g., [25][26][27][28]), NWP and data assimilation (e.g., [29][30][31]), sensitivity analysis, parameter estimations, and predictability studies (e.g., [32][33][34]).The L63 model serves as a prototype for developing coupled nonlinear dynamical systems combining different time scales [11,12,35,36].Coupling of two systems, one with "fast" and another with "slow" time scales, allows imitating the interaction between a fast-oscillating atmosphere and slow-fluctuating ocean.This paper aims to explore essential dynamical, correlative, and spectral properties of this system for time scales applicable to the NWP and climate variability simulations., . . .,   ) is a smooth vector-valued function defined in the domain  ⊆   ,  :  →   , then an autonomous deterministic dynamical system can be described by the following set of ODEs:

Model Equations
where  is known as the phase space of system (1).By introducing the solution operator   having the semigroup property  + =   ⋅   , the system (1) can be written as Thus, the trajectory of the system (1) in its phase space  ∈  is uniquely defined by the initial values of state variables  0 (initial conditions).

Model Equations.
A multiscale nonlinear dynamical system examined in this paper is derived by coupling the fast and slow versions of the original L63 model [13] and can be written as follows [11,12]: where lower case letters , , and  represent the dynamic variables of the fast model, capital letters , , and  denote the state variables of the slow model,  > 0,  > 0, and  > 0 are the parameters of the original L63 model,  is a time-scale factor (if, e.g.,  = 0.1, then the slow subsystem is ten times slower than the fast subsystem),  is a coupling strength for , , , and  variables,   is a coupling strength parameter for  and  variables,  is a "decentring" parameter [11], and  is a parameter representing the amplitude scale factor ( = 1 indicates that slow and fast subsystems have the same amplitude scale).The coupling strength parameters  and   control the interconnection between fast and slow subsystems: the smaller the parameters  and   , the weaker the interdependence between two subsystems.Without loss of generality, one can assume that  = 1,  = 0, and  =   ; then the system of model equations ( 5) can be represented in an operator form as follows: where  = (, , , , , )  is a coupled model state vector,  = (, , , , )  is a model parameter vector, and (; ) is the matrix operator such that System of autonomous ODEs ( 6) has five control parameters (, , , , and ) and together with given initial conditions, represents an initial value problem.

Model Equations in Variational
Form.The model equations in variational form are used to study the system's dynamical properties.These equations can be obtained by linearization of ( 6) around a certain trajectory, which is a particular solution of (6) and represents a vector-function (),  ∈ [0,∞).Assume that the system (6) operates along the trajectory  * () and let () be an infinitesimal perturbation of the state vector such that () = () −  * ().Approximating the right-hand side of ( 6) by a Taylor expansion in the vicinity of  * () and neglecting 2nd order and higher order terms, one can obtain the following set of linear ODEs, the equations in variations: where  = (/)| = * is a Jacobian matrix: Variational equations ( 9) can be rewritten as where   is a linear solution operator.[13,16].These parameter values are used in this study since the atmospheric motions are inherently chaotic.
It is important to note that, for  = 10 and  = 8/3, there is a critical value for parameter , equal to 24.74, and any  larger than 24.74 induces chaotic behaviour of the L63 system [16].The time-scale factor  is taken to be 0.1.Thus, in our study, the main control parameter is the coupling strength , which essentially determines the strength of interactions between fast and slow models and, therefore, the behaviour of the entire coupled system.In numerical experiments values of this parameter have been chosen in accordance with [11,12] and are given in Table 1.Therefore,  ∈ [0.15, 1.0].

Numerical
Integration Procedure.The system of ( 6) is numerically integrated by applying a fourth order Runge-Kutta algorithm with a time step Δ = 0.01.To begin with, (6) are transformed into a discrete-time form (3) and then integrated.This integration produces time series for each of the dynamic variables at equally spaced time-points starting from  0 and denoted by {  = (  ),   = Δ,  = 0, . . .,  − 1}, where Δ is the integration time step, also known as the sampling interval and   = (Δ) −1 is the sampling frequency, also known as the sampling rate.To discard the initial transient period the numerical integration starts at time  − = −2 14 Δ with the initial conditions  − = (0.01; 0.01; 0.01; 0.02; 0.02; 0.02)  (12) and finishes at time  0 .In addition, this ensures that the calculated vector of dynamic variables  0 = ( 0 ) is on the system's attractor.The state vector  0 is then used as the initial conditions for further numerical experiments.Note that, for Δ = 0.01, the numerical integration with length of 100 time steps corresponds to one nondimensional unit of time.
The structure of the resulting attractor depends on the coupling strength parameter .Figures 1 and 2 illustrate phase portraits in -, -, and - phase planes of both fast and slow subsystems for weak ( = 0.15) and strong ( = 0.8) coupling.The evolutions of fast and slow dynamic variables for these values of  are shown in Figures 3 and 4. It is well known that the L63 model produces chaotic oscillations of a switching type: the structure of its attractor contains two regions divided by the stabile manifold of a saddle point in the origin.For relatively small coupling strength parameter ( < 0.5), the attractor for both fast and slow subsystems maintains a chaotic structure, which is inherent in the original L63 attractor.As the parameter  increases, the attractor for both fast and slow subsystems undergoes structural changes breaking the patterns of the original L63 attractor.Fast and slow subsystems affect each other through coupling terms, and at some value of the coupling strength parameter ( > 0.5) a chaotic behavior is destroyed and dynamic variables begin to exhibit some sophisticated motions which are not obviously periodic.Moreover, qualitative examination shows that the evolution through time of both subsystems becomes, to large degree, synchronous (however, phase synchronization requires specific analysis which is not within the scope of this paper).For example, for  = 0.8, the plane phase portraits -, -, and - of the slow subsystem represent closed curves which are mostly smooth and have no visible kinks.These portraits indicate that the motion possesses periodic properties.At the same time, the attractor of the slow subsystem for  = 0.8 has a more complex structure.

Equilibrium
Points.Equilibrium (fixed) points of the system of ODEs (5) occur whenever u = 0.Because this system is homogeneous, it has at least a trivial solution, which is also one of the equilibrium points.In addition, for given parameter values, system (6) has eight nontrivial equilibrium points, all of which can be found numerically.Coordinates of fixed points depend on the parameter  and are listed in Table 2 for  = (0.8; 0.15).The stability of the system in the local vicinity of equilibrium points can be studied by evaluating the Jacobian matrix (10) at each of the fixed points of the system (6) and then finding the resulting eigenvalues.If at least one eigenvalue has a positive real part, the equilibrium  is an unstable node.For instance, for  = 1.0 at the origin, the Jacobian has the following six eigenvalues: Since two of the eigenvalues are positive, the origin is an unstable node and since the remaining four of the eigenvalues are negative this point is saddle.Stability of the system in the vicinity of the remaining fixed points can be examined analogously.
(a) For  <  0 (underdumped oscillator) the general solution of (21) has the following form In this case the solution oscillates with the amplitude  0  − exponentially decreasing to zero and the angular frequency  is given by  =  0 √ 1 − (/ 0 ) 2 .The initial amplitude  0 and the initial phase  0 are defined by the initial conditions.(b) For  >  0 (overdumped oscillator) the general solution of ( 21) is represented by where  1,2 =  ±  0 √ (/ 0 ) 2 − 1 and ,  are unknown constants of integration.In this case, the solution asymptotically tends to the equilibrium  = 0 without oscillations.The values of  and  are fixed by imposing specific initial condition.(c) For  =  0 (critical damping) the general solution of ( 21) is as follows: where ,  are arbitrary constants of integration.Similarly to the overdamped oscillator case, there are no oscillations; only monotonous damping is observed.In our example a damping constant  ≈ 1.47 and the natural frequency  0 depend on the coupling strength parameter  (see Table 3).Critical frequency   0 , at which  =  0 , occurs when  ≈ 1.2.Since for our experiments  ∈ [0.15; 1.0], oscillator ( 21) is overdumped ( >  0 ) and its general solution is given by (23). Figure 5 shows computed solutions to (18) and (19) with initial conditions ( 0 ) = ( 0 ) = 1 for  = 0.8.As  → ∞, the dynamic variables  and  tend to zero.

Dissipativity.
Let  be the volume of some region of phase space.The rate of volume contraction is given by Lie derivative: Since the rate of volume contraction is always negative under the chosen model parameters, the coupled system ( 5) is dissipative and, therefore, volumes in phase space shrink exponentially with time.This is represented mathematically by () = ( 0 )  .

Correlation and Spectral Properties
4.1.General Notes.Numerical integration of (5) produces the time series of dynamic variables hereinafter referred to as signals or oscillations.We can conduct the diagnosis of the coupled dynamical system by analysing the signals using standard tools such as autocorrelation functions and the distribution of power density in the frequency domain from signals obtained in the time domain.
Autocorrelation functions (ACFs) enable one to distinguish between regular and chaotic processes and to detect transition from order to chaos.In particular, for chaotic motions, ACF decreases in time, in many cases exponentially, while, for regular motions, ACF is unchanged or oscillating.In general, however, the behaviour of ACFs of chaotic oscillations is frequently very complicated and depends on many factors (e.g., [37][38][39]).Autocorrelation functions can also be used to define the so-called typical time memory (typical timescale) of a process [40].If it is positive, ACF is considered to have some degree of persistence: a tendency for a system to remain in the same state from one moment in time to the next.
For a given discrete-time signal {  } −1 =0 the power spectrum density (PSD) characterizes the signal intensity (power) per unit of bandwidth.For a wide-sense stationary process, the Wiener-Khinchin theorem relates the ACF to the PSD by means of a Fourier transform (i.e., PSD is a Fourier transform of ACF) and provides information about correlation structure of the time series generated by the system.The term "power spectral density function" is frequently shortened to spectrum.The units of PSD are  2 /Hz, irrespective of what the units of  are.Oscillations of different types have specific spectral properties and, therefore, can be characterized by their PSD.For instance, a periodic motion consisting of the sum of finite number of sine curves has a set of lines in its spectrum, whereas a chaotic motion has a continuous spectral density function.The ACF and spectrum represent different characterization of the same time series information.However, the ACF analyzes information in the time domain and the spectrum in the frequency domain.

Autocorrelation Function. The ACF for a given discrete dynamic variable {𝑢
where the angular brackets denote ensemble averaging.Assuming time series originates from a stationary and ergodic process, ensemble averaging can be replaced by time averaging over a single normal realization: Signal analysis commonly uses the normalized ACF, defined as () = ()/(0).ACF plots for realizations of dynamic variables  and  and  and  calculated for different values of the coupling strength parameter  are presented in Figures 6 and 7, respectively.For relatively small parameter  ( < 0.4), the ACFs for both  and  variables decrease fairly rapidly to zero, consistently with the chaotic behaviour of the coupled system.However, as expected, the rate of decay of the ACF of the slow variable  is less than that of the fast variable .The ACFs for variables  and  (really, their envelopes) also decay almost exponentially from the maximum to zero.For coupling strength parameter on the interval 0.4 <  < 0.6 the ACF of the fast variable  becomes smooth and converges to zero.At the same time, the envelopes of the ACFs of variables , , and  demonstrate a fairly rapid fall, indicating the chaotic behaviour.As the parameter  increases, the ACFs become periodic and their envelopes decay slowly with time, indicating transition to regularity.For  > 0.8 calculated ACFs show periodic signal components.

Power Density Spectra.
There are several methods, both parametric and nonparametric, for spectrum estimation [41,42].This paper uses periodogram, which is the most common where  is a discrete normalized frequency.Then the spectrum can be represented as follows: The spectrum   can be plotted on a dB scale, relative to the reference amplitude  ref = 1; therefore  dB  = 10 ⋅ log 10 (  ) ,  = 0, . . .,  − 1.
The frequency   corresponding to point  of the DFT is The spectrum estimates for fast ( and ) and slow ( and ) variables as well as for weak ( = 0.15) and strong ( = 0.8) coupling are shown in Figures 8 and 9, respectively.For weak coupling, the signal power for all dynamic variables decreases exponentially from low frequencies toward higher frequencies and distinctive energy peaks are not present for almost all variables.The only exception is the fast variable , for which a local peak is observed at frequency ∼1.2 Hz.For weak coupling, the spectrum of the fast subsystem is similar to the spectrum of the L63 model: the fast subsystem generates a broadband spectrum reminiscent of random noise corresponding to irregular aperiodic oscillations.At the same time, the low-frequency component strongly dominates in the spectrum of slow subsystem.As the coupling strength increases, the power spectrum of both fast and slow subsystems shifts toward the low frequencies, which predominate in the spectra.The periodogram is an asymptotically unbiased estimator of the true spectrum; however in statistical terms it is inconsistent since the periodogram variance does not tend to zero as the time series length increases and tends to infinity.The improved estimator of the PSD, according to Welch's method [43], reduces the periodogram variance by averaging.This approach first splits the time series data into overlapping segments, then calculates a modified periodogram for each section, and finally averages periodograms of all segments.The spectrum estimates obtained using Welch's method for fast ( and ) and slow ( and ) variables for weak ( = 0.15) and strong ( = 0.8) coupling are shown in Figure 10.In this figure, the PSDs are represented by curves having no background noise oscillations.Meanwhile, the spectrum patterns remain the same: the signal energy is still concentrated in the low-frequency domain of the spectrum.

Spectrogram.
Spectrogram is another powerful technique used in many applications for estimating the spectrum of the time series data.Spectrogram provides information about power as a function of frequency and time [42] and is generally presented as plot with the frequency of the signal shown on the vertical axis, time on the horizontal axis, and signal power on a colour-scale.Thus, for a given time frame the spectrogram provides the information about frequency content of a signal.Normalized spectrograms for fast (, , and ) and slow (, , and ) variables for weak ( = 0.15) and strong ( = 0.8) coupling are shown in Figures 11 and  12, respectively, with red color representing the highest signal power and blue the lowest.Calculated spectrograms are fully consistent with the ACFs and PSPs, providing additional information about dominant and minor frequencies in the spectrum for a given time.

Concluding Remarks
Low-order system discussed in this paper represents a powerful tool to study various physical and computational aspects of numerical weather prediction, data assimilation, and climate simulation.However, as mentioned previously, the NWP and climate modeling pursue very different objectives and are focused on dynamical processes of significantly different spatial and time scales.The integration time  of the system equations can be classified based on its duration as short, intermediate, long, and very long [32], with the corresponding values of  set to  = 0.1,  = 0.44,  = 2.26, and  = 131.36,respectively.The short integration times traverse some portion of a trajectory along the attractor, the intermediate integrations correspond to complete circle around the attractor, the long integrations complete several circles, and the very long integrations correspond to movement along the attractor of about 100 times.The time step Δ = 0.01 used in the numerical integrations is equivalent to 1.2 hours of a real time [13].Therefore, intermediate and longtime intervals defined above correspond to 2.2 and 11.3 days, respectively, which are consistent with the NWP and data assimilation time of integrations.In turn, the very long integration intervals correspond to climate modeling time scales.
This paper analyzed the basic dynamic, spectral, and correlative properties of the discrete-time nonlinear coupled dynamical system consisting of two versions of the L63 model.The autocorrelation functions, power spectrum densities, and spectrograms for dynamic variables of the fast and slow subsystems were computed by numerical integration of the system equations.The influence of the coupling strength parameter on the ACFs, PSDs, and spectrograms of system's dynamic variables was estimated.By changing the coupling strength parameter , one can obtain the system behaviour that reflects the major dynamical patterns of weather and climate for given natural conditions.

Figure 1 :
Figure 1: Phase portraits for fast and slow subsystems for  = 0.15.

Figure 2 :
Figure 2: Phase portraits for fast and slow subsystems for  = 0.8.

Figure 3 :
Figure 3: Time evolution of fast and slow dynamic variables for  = 0.15.

Figure 4 :
Figure 4: Time evolution of fast and slow dynamic variables for  = 0.8.

Figure 5 :
Figure 5: Time evolution of dynamic variables  and  for  = 0.8.

Figure 6 :Figure 7 :Figure 8 :
Figure 6: Autocorrelation functions for dynamic variables  and  for different parameter .

Figure 9 :
Figure 9: Power spectral density estimate of fast ( and ) and slow ( and ) dynamic variables for  = 0.8.

Figure 10 : 1 ∑
Figure 10: Power spectral density estimate obtained by Welch method for fast ( and ) and slow ( and ) dynamic variables for  = 0.15 (a) and  = 0.8 (b).

Spectrogram of variable xFigure 11 :
Figure 11: Spectrogram for fast and slow dynamic variables for  = 0.15.
Continuous-time deterministic dynamical systems are commonly specified by a set of either autonomous or nonautonomous ODEs; thus, one can consider the problem of the evolution of state variables in time as an initial value problem.If  = ( 1 ,  2 , . . .,   ) is a set of dynamic variables,  is time, and  = ( 1 , 2.1.Preliminary Notes. 2

Table 1 :
[11,12]of the coupling strength parameter  used in the numerical experiments[11,12].The time evolution of coupled nonlinear system considered in this paper is conditioned by a set of ODEs (5) and control parameters , , , , and .By setting parameter  equal to zero, one can restore the original L63 model.Standard values of the L63 parameters corresponding to chaotic behaviour are  = 10,  = 8/3, and  = 28 3.1.Model Parameters.

Table 3 :
Natural frequency  0 for different values of the coupling strength parameter.