Lyapunov Exponents of Two Stochastic Lorenz 63 Systems

TL;DR

This paper compares the effects of two types of stochastic perturbations, SALT and FD noise, on the Lorenz 63 system, developing a stochastic Lyapunov exponent algorithm and revealing distinct impacts on system volume contraction.

Contribution

It introduces a stochastic Lyapunov exponent algorithm and analyzes how SALT and FD noise differently affect the Lorenz 63 system's dynamics.

Findings

FD noise alters the sum of Lyapunov exponents significantly.

SALT noise preserves the original Lorenz 63 Lyapunov sum.

LA SALT yields a deterministic subsystem isomorphic to Lorenz 63.

Abstract

Two different types of perturbations of the Lorenz 63 dynamical system for Rayleigh-Benard convection by multiplicative noise -- called stochastic advection by Lie transport (SALT) noise and fluctuation-dissipation (FD) noise -- are found to produce qualitatively different effects, possibly because the total phase-space volume contraction rates are different. In the process of making this comparison between effects of SALT and FD noise on the Lorenz 63 system, a stochastic version of a robust deterministic numerical algorithm for obtaining the individual numerical Lyapunov exponents was developed. With this stochastic version of the algorithm, the value of the sum of the Lyapunov exponents for the FD noise was found to differ significantly from the value of the deterministic Lorenz 63 system, whereas the SALT noise preserves the Lorenz 63 value with high accuracy. The Lagrangian averaged version of the SALT equations (LA SALT) is found to yield a closed deterministic subsystem for the expected solutions which is found to be isomorphic to the original Lorenz 63 dynamical system. The solutions of the closed chaotic subsystem, in turn, drive a linear stochastic system for the fluctuations of the LA SALT solutions around their expected values.

Tables

Table 1. Table 1: The individual NLEs and sum for σ = 10 𝜎 10 \sigma=10 , r = 28 𝑟 28 r=28 and b = 8 / 3 𝑏 8 3 b=8/3 as computed with the Cayley method and those found in literature.

	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$	${\sum }_{i=1}^3{\lambda }_1$
Cayley method with forward Euler	0.8739	-0.0798	-14.4604	-13.6665
Values according to [Spr03, Spa12]	0.9056	0	-14.5721	-13.6665

Table 2. Table 2: The individual NLEs and sum for σ = 16 𝜎 16 \sigma=16 , r = 45.92 𝑟 45.92 r=45.92 and b = 4 𝑏 4 b=4 as computed with the Cayley method and those found in literature.

	${\lambda }_1$	${\lambda }_2$	${\lambda }_3$	${\sum }_{i=1}^3{\lambda }_1$
Cayley method with forward Euler	1.4858	-0.0721	-22.4135	-20.9998
Values according to [WSSV85]	1.50	0	-22.46	-22.96

Equations

$$\frac{d}{d\tau }X$$

$$\frac{d}{d\tau }Y$$

$$\frac{d}{d\tau }Z$$

$$v(x,t,dW):=u(x,t)dt+\sum _{i=1}^N{\xi }_i(x)\circ d{W}^i.$$

$$\frac{\partial }{\partial t}\text{u}+\text{u}\cdot \nabla \text{u}$$

$$\frac{\partial }{\partial t}T+\text{u}\cdot \nabla T$$

$$\nabla \cdot \text{u}$$

$${\eta }_{t}X:=\eta (X,t)\in {\mathbb{R}}^3,\text{for initial position}x(X,0)={\eta }_{0}X=X.$$

$$\frac{d{\eta }^j(X,t)}{dt}=u^j(\eta (X,t),t)={\eta }_t^{\ast }{u}^j(X,t)=u^j({\eta }_t^{\ast }X,t).$$

$$\frac{d}{dt}{\eta }_t^{\ast }T(X,t)=\frac{d}{dt}T({\eta }_{t}X,t)=\frac{\partial }{\partial t}T(\eta (X,t),t)+\frac{\partial T}{\partial {\eta }^j}\frac{\partial {\eta }^j(X,t)}{\partial t}={\eta }_t^{\ast }\left(\frac{\partial }{\partial t}T+\frac{\partial T}{\partial x^j}{u}^j\right)(X,t).$$

$$\frac{d{\eta }^j(X,t)}{dt}={\tilde{u}}^j(\eta (X,t),t):={\eta }_t^{\ast }{\tilde{u}}^j(X,t)={\tilde{u}}^j({\eta }_t^{\ast }X,t),$$

$$\frac{\partial }{\partial t}\text{u}+\tilde{\text{u}}\cdot \nabla \text{u}$$

$$\frac{\partial }{\partial t}T+\tilde{\text{u}}\cdot \nabla T$$

$$\nabla \cdot \tilde{\text{u}}$$

$$I(t)={\oint }_{c(t)}\text{u}\cdot d\text{x},$$

$$\frac{d}{dt}I(t)={\oint }_{c(t)}\alpha gTdz$$

$$\frac{d}{dt}I(t)$$

$$={\oint }_{c(0)}{\eta }_t^{\ast }\left(\left(\frac{\partial }{\partial t}\text{u}+\tilde{\text{u}}\cdot \nabla \text{u}+u_j\nabla {\tilde{u}}^j\right)\cdot d\text{x}\right)$$

$$={\oint }_{c(t)}\left(\frac{\partial }{\partial t}\text{u}+\tilde{\text{u}}\cdot \nabla \text{u}+u_j\nabla {\tilde{u}}^j\right)\cdot d\text{x}$$

$$={\oint }_{c(t)}\left(-\nabla p-u_j\nabla {\tilde{u}}^j+\alpha gT\nabla z\right)\cdot d\text{x}$$

$$\frac{d}{dt}I(t)$$

$$={\oint }_{c(t)}-d\left(p+\frac{1}{2}|\text{u}{|}^2\right)+\left(\alpha gT\right)dz$$

$$={\oint }_{c(t)}\alpha gTdz.$$

$$\frac{d}{dt}I(t)={\oint }_{c(t)}(\alpha gT\nabla z+\nu \Delta \text{u})\cdot d\text{x},$$

$$I(t)={\oint }_{c(t)}\text{u}\cdot d\text{x},$$

$$\displaystyle\widehat{\textbf{u}}\rightarrow{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{y}_{t}=\textbf{u}\,dt+\sum_{i}{\boldsymbol{\xi}}_{i}\circ dW_{t}^{i}\,,$$

$$\displaystyle{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{u}+{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{y}_{t}\cdot\nabla\textbf{u}+u_{j}\nabla{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}y_{t}^{j}$$

$$\displaystyle{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}T+{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{y}_{t}\cdot\nabla T$$

$$\displaystyle\nabla\cdot{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{y}_{t}$$

$$T(x,z,t)=\overline{T}(z,t)+T^{\prime }(x,z,t),$$

$$\overline{T}(z,t)=\overline{T}(0,t)-\frac{T_{\Delta }}{H}z+{\overline{T}}^{\prime \prime }(z,t),$$

$$T(x,z,t)=\left(\overline{T}(0,t)-\frac{T_{\Delta }}{H}z\right)+T^{\prime }(x,z,t)+{\overline{T}}^{\prime \prime }(z,t).$$

$$\varphi (x,z,t):=T^{\prime }(x,z,t)+{\overline{T}}^{\prime \prime }(z,t).$$

$${\rm{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\phi+{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{y}_{t}\cdot\nabla\phi=\left(\frac{T_{\Delta}}{H}w+\gamma\Delta\phi\right)dt,$$

$${\rm{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\omega+{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{y}_{t}\cdot\nabla\omega=\left(\nu\Delta\omega+\alpha g\phi_{x}\right)dt.$$

$$\displaystyle{\rm{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\omega+\left|\frac{\partial(\widetilde{\psi},\omega)}{\partial(x,z)}\right|$$

$$\displaystyle{\rm{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\phi+\left|\frac{\partial(\widetilde{\psi},\phi)}{\partial(x,z)}\right|$$

Full text

1 Introduction

A fundamental challenge in geophysical fluid dynamics (GFD) is the estimation of combined measurement error and model uncertainty. These errors can arise, for example, from unknown or neglected physical scales, incomplete information in the data and incomplete formulation of the theoretical models. The modern methodology of stochastic data-driven modelling has been developed for managing the loss of predictability associated with error and uncertainty. A fundamental tenet of our approach to this methodology is that stochastic data-driven models should respect the physical principles underlying their deterministic counterparts. For example, if the deterministic physical model follows from Hamilton’s variational principle, one may introduce stochastic data dependence as a constraint on the deterministic variational principle. Recently, [Hol15] followed this approach for fluid dynamics by introducing Stochastic Advection by Lie Transport (SALT) as a constraint on Hamilton’s principle for ideal fluid dynamics. The SALT variational structure preserves both the physical conservation laws and their mathematical implications, such as the Kelvin circulation theorem. Thus, data-driven stochastic fluid flow models based on the SALT approach will obey the same laws of physics as their corresponding deterministic models, although the time dependence for advection by fluid flow will become stochastic.

Now, the solutions of all GFD models of climate and weather tend to be chaotic in nature. That is, some aspects of their solutions tend to show sensitive dependence on initial conditions and on time dependence of forcing. For dynamical systems, this sensitive dynamical dependence is exhibited by positive Lyapunov exponents, meaning that initially nearby solutions diverge away from each other exponentially in time. In addition, GFD models tend to admit bifurcations of solution behaviour depending on their model parameters. Thus, for GFD one must take on the challenge of modelling measurement error and uncertainty which is also nonlinearly concatenated with chaotic deterministic dynamics, is therefore sensitive to initial conditions as well as forcing, and is prone to bifurcations. The most famous example of chaos in GFD is the chaotic dynamics on the strange attractor of the deterministic Lorenz 63 system [Lor63], which may be represented for $(X,Y,Z)^{T}\in\mathbb{R}^{3}$ as the dynamical system,

[TABLE]

where $\sigma$, $r$, $b$, are real constants and $\tau$ is time. The physical meaning of the Lorenz 63 system (1) and some of its main properties will be recalled below. The paper begins with a review of its derivation from the Rayleigh-Bénard model of thermal convection of fluids heated from below which includes an explanation of how its deterministic derivation accommodates the transition to Stochastic Advection by Lie Transport (SALT).

Brief summary. This paper aims to incorporate SALT into the Rayleigh-Bénard model of thermal convection and then follow the standard historical approach to derive the stochastic Lorenz 63 system augmented by SALT. It also compares the dynamical simulation results for the numerical Lyapunov exponents (NLEs) of the SALT Lorenz 63 model with those of the stochastic Lorenz 63 system investigated in [CSG11]. The individual NLEs of the two cases appear to be almost identical for each realisation of the noise. However, the sums are different, so the total phase-space volume contraction rates are different. Namely, the SALT approach preserves the phase-space contraction rate of the underlying deterministic Lorenz 63 system, whereas the stochastic Lorenz 63 system in [CSG11] varies with each realisation. The implications of this difference in temporal behaviour of total contraction rates are elicited in numerical simulations. After a brief Conclusion section which summarises the main resutls, the paper’s Outlook section briefly discusses a toy climate model based on the SALT Lorenz 63 system. In particular, the relations between ‘what you expect’ and the dynamical evolution of the statistics in ‘what you get’ are discussed in the climate context.

Background

In [Hol15], SALT was introduced by using Hamilton’s variational principle for fluids, constrained to enforce stochastic Lagrangian fluid trajectories, arising as characteristic curves of the following stochastic Eulerian velocity vector field

[TABLE]

The velocity vector field (2) is regarded as a decomposition into a deterministic drift velocity $u(x,t)$ and a sum over spatially dependent but temporally constant vector field amplitudes of Stratonovich stochastic terms. The $\xi_{i}(x)$ are taken as stationary eigenvectors of the velocity-velocity correlation tensor. In practice, the $\xi_{i}(x)$ need to be supplied a priori. See [CCH*+*18, CCH*+*19a, CCH*+*19b] for how the $\xi_{i}(x)$ are determined in practical applications. Stratonovich stochastic calculus is preferred in the derivation of the SALT equations, because its standard calculus properties (product rule and chain rule) also admit variational calculus. The decomposition (2) was shown in [CGH17] to also arise from homogenisation of the velocity obtained by writing the deterministic time-dependent Lagrange-to-Euler diffeomorphism as the composition of two diffeomorphisms, a fluctuating map with two time scales $(t,t/\epsilon)$ where $\epsilon\ll 1$, applied to a mean flow map depending only on the slow time scale $t$. In this analysis, the $\xi_{i}(x)$ should be understood as empirical orthogonal functions which correspond to the different modes of the fast flow. The rigorous application of homogenisation theory requires added assumptions of mildly chaotic fast small-scale deterministic dynamics, as well as a centering condition, according to which the mean of the fluctuating deviations is small, when pulled back to the mean flow. The homogenisation analysis performed in [CGH17] gives rise to Itô noise, which is then transformed to the Stratonovich type.

The work of [CGH17] justified regarding the Eulerian vector field in (2) as a genuine decomposition of the fluid velocity into a sum of drift and stochastic parts, rather than simply a perturbation of the dynamics meant to model unknown effects in uncertainty quantification. Consequently, one should expect that the properties of the fluid equations with SALT should closely track the properties of the unapproximated solutions of the fluid equations. For example, if the unapproximated model equations are Hamiltonian, then the model equations with SALT should also be Hamiltonian. This was shown in [Hol15].

In [CCH*+*18] and [CCH*+*19a], respectively, a 2-layer quasi-geostrophic model and the 2D Euler equations were each perturbed using the SALT velocity decomposition in (2). These investigations were aimed at developing a new methodology for uncertainty quantification and data assimilation for ideal fluids. The new methodology tested in these papers was found to be suitable for coarse graining in both cases. Specifically, uncertainty quantification tests of the velocity decomposition of [CGH17] were performed by comparing ensembles of coarse-grid realisations of solutions of the resulting stochastic partial differential equations with the “true solutions” of the deterministic fluid partial differential equations obtained by computing the same problem at higher resolution. The time discretisations used for approximating the solution of the stochastic partial differential equations were shown to be consistent in each case and comprehensive numerical tests confirmed the non-Gaussianity of the flows under SALT dynamics.

In [CFH19] it is shown that the Euler equations for an ideal, inviscid fluid with SALT are locally well-posed in regular spaces and a Beale-Kato-Majda type blow-up criterion is proved. Thus the analytical properties of the 3D Euler fluid equations with SALT closely mimic the corresponding analytical properties of the original deterministic 3D Euler equations.

Inspired by spatiotemporal observations from satellites of the trajectories of objects drifting near the surface of the ocean in the National Oceanic and Atmospheric Administration’s (NOAA) Global Drifter program, in [GBH18] a data-driven stochastic model of geophysical fluid dynamics was developed. Here non-stationary spatial correlations represent the effects of advection by ocean currents. The methods in this paper are similar to those in [Hol15], where the models were derived using reduction by symmetry of stochastic variational principles. The corresponding momentum maps, conservation laws and Lie-Poisson bracket structures were used in developing the new stochastic Hamiltonian models of geophysical fluid dynamics with nonlinearly evolving stochastic properties.

For motivation of the present work, we recall that the complex bifurcation behaviour shown in the Rayleigh-Bénard model of thermally driven convection is quite challenging to simulate numerically, as discussed, e.g., in [KBF*+*18, PJMS16, OSG*+*13, SW07]. One might hope that the SALT approach could play a useful role for investigating this convection model further. In particular, the large amount of available experimental and numerical data for thermally driven convection could, in principle, facilitate determination of the necessary correlation eigenvectors $\xi_{i}(x)$ in the decomposition of the fluid velocity vector field above in equation (2). One would expect progress in this matter to produce stochastic models which are numerically much less expensive than the corresponding deterministic models. In this paper, we shall derive the SALT modification of the Rayleigh-Bénard convection equations and then derive the corresponding Lorenz 63 model. We hope this study will stimulate future research in using SALT to develop efficient stochastic numerical simulations of Rayleigh-Bénard convection. For example, the effects of stochasticity implemented by the SALT approach on the bifurcation structure of Rayleigh-Bénard convection could be an avenue for further investigations of interest in fluid dynamics for both theory and applications. However, such an investigation is beyond the scope of the present paper.

Plan of the paper

In section 2 of this paper, SALT is introduced into the formulation of the Rayleigh-Bénard convection, by using the methods from geometric mechanics [Hol15]. In particular, the Kelvin circulation theorem is used to introduce SALT into the Oberbeck-Boussinesq (OB) equations which govern the Rayleigh-Bénard model. In section 3, the OB equations with SALT are restricted to a vertical slice and by means of a truncated spatial Fourier series expansion are projected onto low wavenumber modes, following the approach of [Lor63]. This process gives rise to a stochastic version of the Lorenz 63 system. In section 4, we compare this stochastic Lorenz 63 model derived using the SALT approach with an alternative stochastic Lorenz 63 model found in the literature [CSG11], in which linear multiplicative noise is added into the equation for each variable in (1). By using the theory of Lyapunov exponents and Oseledet’s multiplicative ergodic theorem, analytical statements can be derived to describe the average rate of contraction of phase-space volume. In section 5, a numerical method to compute the numerical Lyapunov exponents (NLEs) is introduced, based on a QR method involving the Cayley-transform, as shown in [UvB01]. This method is generalised to deal with systems of stochastic differential equations and is shown to give accurate and robust results. Section 6 summarises our conclusions. Section 7 briefly discusses an intriguing outlook for modifications of the present approach which could conceivably apply in climate science.

2 Rayleigh-Bénard convection and Lorenz 63

In this paper, the goal is to look at the effect of SALT on the Lorenz 63 system in (1). The Lorenz 63 system can be systematically derived by expanding the equations of Rayleigh-Bénard convection into Fourier series and then truncating. The Rayleigh-Bénard model may be interpreted as a simple, local, weather system, where only heat and wind are involved. One of the important ingredients that is missing from the model is the effect of rotation, that gives rise to the Coriolis force and hence we can only think of this as a local model for weather effects. The domain is a three dimensional box, where the bottom plate is being heated and the top plated is being cooled at a constant temperature. This corresponds to the earth heating the surface of the earth, which heats the air directly above it. The hot air has a lower density than the cool air which is above it. Gravity acts to restore stability, causing the hot air to rise and the cool air to descend. The top plate is being cooled, corresponding to the cold air above the domain which acts as a heat sink. The fluid in the box is modelled with the Navier-Stokes equations under the Oberbeck-Boussinesq approximation. The fluid is assumed to have a constant heat capacity, so the advection-diffusion equation for the heat can be expressed in terms of temperature. The governing set of equations are

[TABLE]

The buoyancy force $\textbf{F}=\alpha gT\hat{\textbf{z}}$ in the top equation in (3) acts in the vertical direction and depends on the thermal expansion coefficient $\alpha$, gravity $g$ and the local temperature $T$. This model is dissipative in nature due to the linear diffusion of momentum per unit mass u, by viscosity $\nu$, and heat per unit mass $T$, by heat diffusivity $\gamma$. The dissipation makes the dynamics irreversible, which does not fit the mathematical framework that is necessary to introduce SALT. At this stage, let us drop the dissipative terms temporarily, to illustrate the framework. Suppose $T$ is the heat per unit mass of a Lagrangian fluid parcel following the flow given by the smooth, invertible, time-dependent Lagrange-to-Euler map,

[TABLE]

By using the pullback operation, denoted by $\eta_{t}^{*}$, the temperature equation can be rewritten in terms of the Lagrange-to-Euler map (4). In particular, for the parcel occupying spatial position $x\in\mathbb{R}^{3}$ at time $t$, we have $\eta_{t}^{*}T(X,t)=T(\eta(X,t),t)$. The Eulerian velocity vector along a Lagrangian path may be written in terms of the flow $\eta_{t}$ and its pullback $\eta_{t}^{*}$ as

[TABLE]

The time derivative of the pullback of $\eta_{t}$ for the scalar function $T(X,t)$ is given by the chain rule as

[TABLE]

Equation (6) shows the geometric nature of advection of heat per unit mass along a Lagrangian particle path. A similar relation can be shown to exist for the momentum per unit mass u, although this relation involves an additional term, as we shall see below in the proof of Theorem 2.2.

Remark 2.1** (Distinguishing between flow velocity and momentum per unit mass)**

In the momentum equation in (3), two quantities with the dimensions of velocity appear, both noted as u. In Cartesian coordinates on $\mathbb{R}^{n}$ equipped with the Euclidean metric, when the kinetic energy is given by the $L^{2}$ metric, this is valid. In different metrics or different spaces, one has to distinguish between the contravariant vector field $u^{j}\partial_{j}$, which transports fluid properties and the covariant momentum per unit mass $u_{j}dx^{j}$, which is transported. To emphasize this difference, we will henceforth denote the vector field, or transport velocity, by $\widetilde{\textbf{u}}$ and write momentum per unit mass by u. Then equation (5) is written as

[TABLE]

and we can write the ideal Oberbeck-Boussinesq equations for the Rayleigh-Bénard model as

[TABLE]

The distinction between $\widetilde{\textbf{u}}$ and u is particularly important in the Kelvin circulation theorem, where $\widetilde{\textbf{u}}$ transports the fluid loop and u will be the circulation velocity, which is integrated around the loop.

Theorem 2.2** (Kelvin circulation theorem for the ideal Oberbeck-Boussinesq equations)**

For the circulation integral

[TABLE]

where $c(t)$ is a closed loop that is moving with velocity $\widetilde{\textbf{u}}$, the ideal Oberbeck-Boussinesq equations given in (8) imply the following circulation dynamics,

[TABLE]

Proof.

[TABLE]

Since the space that we are working in is $\mathbb{R}^{3}$, with Cartesian coordinates, an Euclidean metric and the kinetic energy is given by the $L^{2}$ metric, $\textbf{u}_{j}$ and $\textbf{u}^{j}$ can be identified as vectors. The identity $u_{j}\nabla u^{j}=\frac{1}{2}\nabla|\textbf{u}|^{2}$ is then valid and can be used to find

[TABLE]

In the last step we have used the fundamental theorem of calculus.

The dissipative Oberbeck-Boussinesq equations (3) possess the following circulation dynamics.

Corollary 2.3** (Kelvin circulation theorem for the dissipative Oberbeck-Boussinesq equations)**

The dissipative Oberbeck-Boussinesq equations (3) possess the following circulation dynamics,

[TABLE]

for the circulation integral given by

[TABLE]

where $c(t)$ is a closed loop moving with transport velocity $\widetilde{\textbf{u}}$.

Remark 2.4

The Kelvin theorem (9) is secretly Newton’s second law $\frac{d\mathbf{P}}{dt}=\mathbf{F}$ for the momentum per unit mass $\mathbf{P}$ and the force per unit mass $\mathbf{F}$ for masses distributed on a space of loops. From this viewpoint, the circulation is momentum per unit mass and the right hand side of (9) is the force per unit mass.

We now replace the velocity $\widetilde{\textbf{u}}$ that transports the closed loop $c(t)$ with a stochastic process for the Lagrangian trajectory given by

[TABLE]

in which the assumption is made that $\nabla\cdot{\boldsymbol{\xi}}_{i}=0$ for all $i=1,\ldots,N$. This replacement introduces stochasticity into the dissipative Oberbeck-Boussinesq equations (3) as follows

[TABLE]

In stochastic fluid dynamics with SALT, the advective velocity field transforms, as we have seen. The framework for this transformation can be found in [Hol15]. The Eulerian velocity field becomes the stochastic vector field {\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{y}_{t}, which also advects the loop in the Kelvin circulation theorem. Here, this change was applied to the Oberbeck-Boussinesq equations to derive their stochastic analogue with SALT. The development of model equations for stochastic fluid dynamics revolves around the choice of forces appearing in Newton’s second law and Kelvin’s circulation theorem. A similar approach can be found in [Mém14, Res17], where the goal is stochastic turbulence modelling. In [CDMR18] a stochastic version of Lorenz 63 can be found with a similar structure, as the noise was introduced carefully in the transport terms.

2.1 Lorenz 63 equations

In order to be able to project the equations onto the Fourier modes used in [Lor63], the OB equations in (11) need to written in terms of scalar vorticity and temperature profile, as in [Sal62]. The temperature $T(x,z,t)$ at constant $y$ can be expanded into a horizontal mean value and a departure from that mean

[TABLE]

where $\overline{T}$ is the horizontal mean and $T^{\prime}$ is the departure from the mean. Also the mean temperature can be decomposed into a linear difference between the lower and upper boundary and a perturbation of this linear difference. This gives

[TABLE]

where $T_{\Delta}$ is the constant temperature difference between the lower and upper boundary, $H$ is the height between the vertical boundaries and $\overline{T}^{\prime\prime}$ is the perturbation of the linear difference. The temperature can therefore be written as

[TABLE]

The temperature profile is defined as

[TABLE]

Substitution of equation (14) into the convection-diffusion equation for the heat yields

[TABLE]

where $w$ is the vertical component of the velocity field. The vorticity equation is obtained by taking the curl of the momentum equation in (11). This eliminates the pressure term and since the domain is two dimensional, the vorticity is a scalar function. This gives

[TABLE]

The velocity terms in (16) and (17) can be rewritten in terms of streamfunctions. The streamfunction associated to the deterministic velocity field u will be denoted by $\psi$ and the streamfunction associated to the stochastic vector field {\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}\textbf{y}_{t} will be denoted by $\widetilde{\psi}$,

[TABLE]

This is the same formulation of the Rayleigh-Bénard convection problem as in [Sal62, Lor63]. The nonlinear terms have been rewritten in terms of a determinant of a Jacobian. Lorenz used the following truncated Fourier series

[TABLE]

In these Fourier expansions, $k$ is the wave number, $R_{a}=\alpha gH^{3}T_{\Delta}\nu^{-1}\gamma^{-1}$ is the Rayleigh number and $R_{c}=\pi^{4}k^{-2}(1+k^{2})^{3}$ is the critical value of the Rayleigh number. These scaling constants are introduced so that the resulting equations take a compact form. We take the truncated Fourier series of the stochastic streamfunction $\widetilde{\psi}$ to be the same as the expansion of the deterministic streamfunction $\psi$, but with a stochastic coefficient. The motivation for this choice is that from a physical point of view, we do not want the stochasticity to give rise to types of motion other than rolls between the two plates. The expansion of the stochastic streamfunction then is

[TABLE]

Carrying out the projection yields the Lorenz 63 system, augmented with SALT,

[TABLE]

where $\sigma=\gamma\nu^{-1}$ is the Prandtl number, $r=R_{a}R_{c}^{-1}$ is a scaled Rayleigh number and $b=(4(1+k^{2}))^{-1}$ is a parameter related to the wavenumber. Note that there is no stochasticity in the term proportional to $r$. The equations in (21) were obtained by projecting (11) onto Fourier modes. In (11), the stochasticity is in the nonlinear terms and not in the term that is proportional to $\alpha g$. The term proportional to $\alpha g$ is the term that corresponds to the term proportional to $r$ in (21). The time $\tau=\pi^{2}(1+k^{2})\gamma tH^{-2}$ is dimensionless. Hereafter, we will not distinguish between $\tau$ and $t$ and just write $t$. In (18), observe that the noise appears only in the nonlinear transport terms. We can formally rewrite (21) to show that also the associated Lorenz equations have this property,

[TABLE]

upon introducing the notation, $\widetilde{X}=X+\beta\circ d\dot{W}_{t}$. The nonlinear transport terms in (21) and (22) represent circulation as rigid rotation of the $Y,Z$ plane around the $X$-axis at angular velocity $\widetilde{X}$. Physically, then, the SALT modification may be interpreted as a stochastic transport around the $X$-axis with angular velocity of circulation given by $\beta\circ dW_{\tau}$. The Itô form of (21) and (22) is

[TABLE]

The first equation in each of the equation sets (21), (22) and (23) contains no noise. This is because the Jacobian projects to zero for this Fourier expansion. A stochastic L63 system similar to that in (23) was also derived using the “location uncertainty” approach in [CDMR18].

In [CSG11], a stochastic version of the Lorenz system was introduced and found to possess a pullback attractor that supports a random Sinai-Ruelle-Bowen (SRB) measure. We will not try to compute the SRB measure for our version of a stochastic Lorenz system. Instead, we will compute the numerical Lyapunov exponents for both these systems. In [CSG11], linear multiplicative Itô noise is added in each variable, and the stochastic Lorenz system is formulated as

[TABLE]

When comparing linear terms in (24), one sees that the stochasticity is paired with the dissipation. Hence we call this noise fluctuation-dissipation noise.

3 Lyapunov Exponents

We now want to compare the two version of stochastic Lorenz 63 in terms of their Lyapunov exponents. Firstly, we will determine the sum of the Lyapunov exponents analytically, since this is possible for the Lorenz 63 system. To be able to compute Lyapunov exponents, a number of conditions are necessary. First of all, the system of equations needs to generate a random dynamical system (RDS), which is a tuple $(\varphi,\vartheta)$, where $\varphi$ is a cocycle, the solution of the dynamical system $\vartheta$. Additionally, an integrability condition,

[TABLE]

where $\log^{+}x:=\max(0,\log x)$, needs to be satisfied to make sure that the linear equation associated to the system that we are considering is wellposed. The integrability condition is sufficient for Oseledet’s multiplicative ergodic theorem (MET), which implies the regularity and existence of Lyapunov exponents. A stochastic differential equation is locally wellposed if its solutions exist locally and are unique up to indistinguishability. Sufficient conditions are local Lipschitz continuity and a linear growth condition. The coefficients of the stochastic differential equations in (21) and (24) are of polynomial type, and therefore smooth. This implies local Lipschitz continuity. Both systems satisfy linear growth. After introducing the following notation, the following theorem implies that both (21) and (24) generate a RDS. In general, a Stratonovich stochastic differential equation can be written as

[TABLE]

with the convention $dW_{t}^{0}=dt$ to allow for this shorthand. From [Arn03], we have the following theorem.

Theorem 3.1** (RDS from Stratonovich SDE)**

Let $f_{0}\in\pazocal{C}^{k,\delta}_{b}$, $f_{1},\ldots,f_{m}\in\pazocal{C}^{k+1,\delta}_{b}$ and $\sum_{j=1}^{m}\sum_{i=1}^{d}f_{j}^{i}\frac{\partial}{\partial x_{i}}f_{j}\in\pazocal{C}^{k,\delta}_{b}$ for some $k\geq 1$ and $\delta>0$. Here $\pazocal{C}^{k,\delta}_{b}$ is the Banach space of $\pazocal{C}^{k}$ vector fields on $\mathbb{R}^{d}$ with linear growth and bounded derivatives up to order $k$ and the $k$-th derivative is $\delta$-Hölder continuous. Then:

1. i)

[TABLE]

generates a unique $\pazocal{C}^{k}$ RDS $\varphi$ over the dynamical system (DS) describing Brownian Motion (the background theory for this can be found in **[Arn03, Elw78]**). For any $\epsilon\in(0,\delta)$, $\varphi$ is a $\pazocal{C}^{k,\epsilon}$-semimartingale cocycle and $(t,x)\to\varphi(t,\omega)x$ belongs to $\pazocal{C}^{0,\beta;k,\epsilon}$ for all $\beta<\frac{1}{2}$ and $\epsilon<\delta$. 2. ii)

The RDS $\varphi$ has stationary independent (multiplicative) increments, i.e. for all $t_{0}\leq t_{1}\leq\ldots\leq t_{n}$, the random variables

[TABLE]

are independent and the law of $\varphi(t+h)\circ\varphi(t)^{-1}$ is independent of $t$. Here $\circ$ means composition. 3. iii)

If $D\varphi(t,\omega,x)$ denotes the Jacobian of $\varphi(t,\omega)$ at $x$, then $(\varphi,D\varphi)$ is a $\pazocal{C}^{k-1}$ RDS uniquely generated by (27) together with

[TABLE]

Hence $D\varphi$ uniquely solves the variational Stratonovich SDE on $\mathbb{R}$

[TABLE]

and is thus a matrix cocycle over $\Theta=(\vartheta,\varphi)$. 4. iv)

The determinant $\det D\varphi(t,\omega,x)$ satisfies Liouville’s equation on $\mathbb{R}$

[TABLE]

and is thus a scalar cocycle over $\Theta$.

The background theory and proof of this theorem can be found in [Arn03]. The variational equation is commonly referred to as the tangent model. The matrix solutions of the tangent model are commonly called the fundamental matrix of solutions. Although both (21) and (24) satisfy the theorem above, an additional observation is required to guarantee that the solutions do not blow up. This is can be shown by a Lyapunov function argument, which will imply that deterministic Lorenz equations have a global attractor set. Together with the local existence and uniqueness of strong solutions to the system of SDEs, this is enough to satisfy the integrability condition. To prove the existence of a globally attracting set, one considers the Lyapunov function [Spa12]

[TABLE]

Dividing the total time derivative of $V(\textbf{X})$ by $2r^{2}\sigma b$ yields the equation for an ellipsoid

[TABLE]

This shows that $\dot{V}$ is negative outside of the ellipsoid and positive inside the ellipsoid given by

[TABLE]

So inside the ellipsoid the dynamics are unstable, as there is no converging behaviour. Outside of the ellipsoid, where $\dot{V}<0$, the dynamics converge towards the ellipsoid. Hence $V(\textbf{X})$ is a Lyapunov function outside of an ellipsoid. This proves that no finite time blow-up can occur for the deterministic case. Since the transport noise and fluctuation-dissipation noise Lorenz systems both satisfy linear growth, the stochastic versions also do not blow up. We are now able to prove that both (21) and (24) satisfy the integrability condition (25). For all finite systems of SDEs (so no stochastic partial differential equations), the Jacobian of the dynamics is a square matrix. Since in $\mathbb{R}^{d\times d}$ all norms are equivalent, the condition is satisfied or dissatisfied for all norms simultaneously. The Jacobian for (21) is

[TABLE]

and for (24) the Jacobian is given by

[TABLE]

For Oseledet’s MET to be valid we require $\|\sum_{j=0}^{m}Df_{j}\|\in L^{1}$. This is true if all elements of the matrices are in $L^{1}$. This conditions is violated if any of the elements of the matrix is unbounded, since then the argument of the logarithm would become unbounded. The Lyapunov function has proven that the dynamics have a global attractor and we have established local existence and uniqueness, so for any initial condition, the dynamics stay bounded. Therefore the integrability condition (25) is satisfied and regularity and existence of Lyapunov exponents is guaranteed. Oseledet’s MET now states that, for $v_{t}$ the solution of (28), $\lim_{t\to\infty}(v_{t}(\omega)^{T}v_{t}(\omega))^{1/2t}=:\Phi(\omega)\geq 0$ exists and the logarithm of the eigenvalues of $\Phi$ are the Lyapunov exponents.

3.1 Sum of Lyapunov exponents

By definition of Lyapunov exponents and by using Liouville’s equation (30) (also called Abel-Jacobi-Liouville formula), the following important fact is derived.

Lemma 3.2

If the trace of the Jacobian $Df_{0}$ is constant and the trace of $Df_{j}$ for $j\geq 1$ is zero, then the sum of the Lyapunov exponents is equal to the trace of $Df_{0}$.

Proof. By taking the determinant of $\Phi$ and using several properties of the determinant for square matrices, we can show

[TABLE]

First, for any square matrix $A$, $\det(A^{T})=\det(A)$, which lets us write $\det(v_{t}^{T}v_{t})=\det(v_{t}^{2})$. Secondly, for any square matrix $A$ and $B$, we have $\det(AB)=\det(A)\det(B)$, which allows us to write $\det(v_{t}^{2})=\det(v_{t})^{2}$. Finally, the determinant of a square matrix is related to the eigenvalues of that square matrix by $\det(A)=\prod_{i}\kappa_{i}$, where $\kappa_{i}$ are the eigenvalues of A. Now let $e^{\gamma}_{i}$ be the eigenvalues of $v_{t}$. Then $\gamma_{i}$ are the unaveraged Lyapunov exponents. Using Liouville’s equation (30) and the right hand side of (36)

[TABLE]

Since the trace of the Jacobian $Df_{0}$ was assumed to be constant and the trace of Jacobian $Df_{j}$ is zero for all $j\geq 1$, by taking the logarithm, we find

[TABLE]

where $\lambda_{i}$ are the Lyapunov exponents by definition.

This lemma applies to the Lorenz 63 system with SALT (21), where it implies that the sum of the Lyapunov exponents is equal to that of the deterministic system

[TABLE]

This conclusion does not hold for the Lorenz 63 system with fluctuation dissipation noise (24), because for the latter the trace of $Df_{1}$ is nonzero. Still one can use the Liouville equation (30) and a similar computation to that in the proof of the lemma, to obtain the sum of the Lyapunov exponents for (24)

[TABLE]

The sum of Lyapunov exponents represents the average rate of expansion or contraction of phase-space volume. Hence this result shows on a theoretical level that the phase-space contraction (or expansion) properties of the two systems are different for any finite time. In the limit, the two systems have the same properties.

4 Method to compute Lyapunov exponents

To compute the finite time approximation to the Lyapunov exponents (which we will refer to as numerical Lyapunov exponents, or NLEs), one needs to simultaneously solve the governing dynamics and the corresponding variational equation. When the dynamics are given by a system of differential equations, the variational equation becomes a matrix differential equation. The appropriate initial condition for the variational equation is usually the identity matrix, as this corresponds to evolving the unit ball along the linearised dynamics. The ball deforms and its average deformation is associated to the NLEs. However, directly solving the variational equation is problematic, as the vectors associated to the NLEs tend to align along the direction of largest increase. Regularly orthonormalising avoids this issue. The well known QR method is therefore the go-to option for solving the variational equation. The QR method dictates that the matrix $v_{t}$ is decomposed into an orthogonal matrix $Q\in O(n):=\{Q\in\mathbb{R}^{n\times n}:\det Q=\pm 1\}$ and an upper triangular matrix $R$ at every time step. Consider the Stratonovich SDE on $\mathbb{R}^{n}$ given by

[TABLE]

where the functions $f_{j}$ are sufficiently regular to guarantee wellposedness of the SDE. Again the convention $dW_{t}^{0}=dt$ is used. The corresponding variational equation is

[TABLE]

where $J_{j}:=Df_{j}(Y_{t})$ is the Jacobian of the dynamical system and $I\in\mathbb{R}^{n\times n}$ is the identity matrix. Applying the QR decomposition, multiplying by $Q^{T}$ from the left and by $R^{-1}$ from the right then yields

[TABLE]

The first matrix on the left hand side is skew-symmetric and the second matrix is upper triangular. This fact will be used in solving for $Q$ and $R$ independently. It is necessary to take a new $QR$-decomposition once every so often, otherwise the matrix $Q$ may lose its orthogonality and cause the algorithm to break down. In addition, for high-dimensional dynamical systems, by means of the Cayley transform, in [UvB01] the $QR$-method is slightly adapted to gain a small speed-up. The Cayley transform will provide an orthogonal matrix with determinant +1 as long as the eigenvalues do not approach -1. This may occur during the time integration, so a restarting procedure can avoid this potential problem. The restarting procedure is possible due to the following lemma.

In (28), for $t>t_{0}$ set $v_{t_{0}}=Q_{0}R_{0}$, where $Q_{0}$ is orthogonal and $R_{0}$ is upper triangular with all diagonal elements positive. As in [UvB01], the real line is divided into subintervals $t_{i}\leq t\leq t_{i+t}$ for $i=1,2,\ldots,$ so that each interval has length $\Delta t_{i}=t_{i+1}-t_{i}$. The solution $v_{t_{i}}$ to the variational equation (28) at time $t_{i}$ can be decomposed as $v_{t_{i}}=Q_{i}R_{i}$ for $i=1,2,\ldots$. We can now introduce the following lemma.

Lemma 4.1

At any time $t=t_{i}+\tau$, $0\leq\tau\leq\Delta t_{i}$, for $i=1,2,\ldots$, the solution of the variational equation (28) can be expressed as

[TABLE]

where $\widetilde{v}_{\tau}$ is the solution to the differential equation given by

[TABLE]

with $Q_{0}=I$, $R_{0}=I$ and $\widetilde{J}_{j}(\tau)=Q_{i}^{T}J_{j}(t_{i}+\tau)Q_{i}$.

The proof of this lemma can be found in [UvB01]. Although in that proof the variational equation is deterministic, the stochastic case is straightforwardly found from the deterministic one, as the only change is the variational equation itself. The Cayley transform is defined as

[TABLE]

where $I\in\mathbb{R}^{n\times n}$ is the identity matrix and $K\in\mathbb{R}^{n\times n}$ is a skew-symmetric matrix. An important feature of $(I-K)$ and $(I+K)^{-1}$ is that they commute. The transformation (46) is valid as long as none of the eigenvalues of $Q$ are equal to -1. By deriving the stochastic differential equation for $K$, the Cayley method takes form. Since the initial condition for $Q(0)=I$, the initial condition for $K(0)=0$. By taking the stochastic evolution differential of $Q$, we obtain

[TABLE]

Using the skew-symmetry of $K$, the distributive property of the transpose, the fact that for any invertible matrix $A$ we have that $(A^{T})^{-1}=(A^{-1})^{T}$ and by writing $(I+K)=-(I-K)+2I$ and setting $H:=(I+K)^{-1}$, one obtains the following expression

[TABLE]

Upon introducing the notation $G:=(I-K)$ and using $H$ as before, the right hand side of (43) is expressed as

[TABLE]

Substitution of (48) and (49) into (43) yields the following differential equation

[TABLE]

Recall that the first matrix on the left hand side is skew-symmetric and the second matrix is upper triangular. Let S=H^{T}{\sf{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}d}}KH so that

[TABLE]

which determines the differential equation for $K$ as

[TABLE]

Observe that since $K$ is skew-symmetric, it is determined by the lower triangular part of $G^{T}SG$. Now that $K$ is known, the Lyapunov exponents are determined as the averages of the solutions of the differential equation for $\rho_{a}:=\log(R_{aa})$,

[TABLE]

where $h_{a}$ are the columns of $H=[h_{1}\quad h_{2}\quad\ldots\quad h_{n}]$. The Lyapunov exponents are then found as

[TABLE]

Recall that this solution method is valid as long as the eigenvalues of $Q$ do not equal -1. As the initial condition is $Q(0)=I$, there is always an interval of time $0\leq t\leq t_{0}$ in which the condition for the Cayley transform (46) is not violated. The following condition for restarting the algorithm is introduced: let $\eta\in(0,1)$ be chosen by the user such that $\|K\|\leq\eta<1$ for some suitable norm. At time $t_{0}$, when the norm of $K$ equals $\eta$, $Q(t_{0})=:Q_{0}$ is computed and stored. The algorithm is then restarted at time $t_{0}$, where due to the lemma we have

[TABLE]

for $t_{0}\leq\tau$. Besides the adapted Jacobian, this is the same equation and hence can be solved using the Cayley method as long as the norm of $K$ is smaller than our chosen value for $\eta$. The initial condition for equation (53) changes to $\rho(0)=\rho(t_{0})$.

5 Numerical Results

The Lorenz system has been studied intensively with the standard parameter values $\sigma=10$, $r=28$ and $b=8/3$, [Lor63, AS01], though in the latter two papers for an adapted version of the Lorenz system. Lorenz showed that for these parameter values, the deterministic Lorenz has a strange attractor. In [WSSV85], the Lorenz system is studied for the nonstandard parameter values $\sigma=16$, $r=45.92$ and $b=4$. Upon introducing stochasticity, one can no longer speak of an attractor in the standard sense, since the noise, due to the unbounded variation nature of the Wiener process, will push the dynamics out of any bounded set almost surely. We set the initial condition to $(X(0),Y(0),Z(0))=(0,1,0)$ and evolve the system for 50000 time steps with time step size $\Delta t=0.001$ with the Euler-Maruyama method. This sets the initial condition for determining the numerical Lyapunov exponents. We solve the SDEs in the Cayley method with the same time step $\Delta t$ for $10^{5}$ iterations in total. The norm tolerance for the matrix $K$ is set to $\eta=0.8$. It is known that the Euler-Maruyama method has poor convergence, so the individual exponents can be computed much more accurately if one were to implement a better numerical method. For the deterministic case, the individual exponents agree reasonably well with the existing literature. The sum of the numerical Lyapunov exponents turns out to be a very robust value, as even the Euler-Maruyama method establishes the correct value to high accuracy.

5.1 Deterministic Case

When there is no noise, the Liouville equation implies that sum of the Lyapunov exponents is equal to the trace of the Jacobian of the dynamics. For the standard parameter values $\sigma=10,r=28$ and $b=8/3$, the sum of the Lyapunov exponents is

[TABLE]

The individual NLEs for the deterministic Lorenz 63 system have been determined by [Spr03] who used a 4th order Runge-Kutta method with a fixed step size of 0.001, performed over $10^{9}$ iterations. The Cayley method allows us to determine the individual NLEs, where we solve the dynamics with a forward Euler scheme. We find the following values

The values are not exactly the same for the individual NLEs. This is due to the poor convergence of the numerical methods used here (Euler-Maruyama has order 1/2 convergence). However, the sum is the same in four decimal places. As an additional comparison, we also study the deterministic system for the nonstandard parameter values used in [WSSV85], $\sigma=16,r=45.92$ and $b=4$, where the sum is

[TABLE]

In this situation, we find the following values

Again the values are not the same for the individual NLEs, but the sum is accurate. For both sets of parameter values, the individual NLEs are computed in good agreement with those found in literature.

5.2 Stochastic Cases

Here the Lorenz system is studied with SALT. The noise amplitude is chosen to be $\beta=0.5$. In figure 1 it can be seen that the stochastic dynamics give rise to a perturbed butterfly shaped object in phase space. The numerical Lyapunov exponents converge. The $x$-axis in the convergence plot of the NLEs shows time, which is simply the number of iterations multiplied by the time step. The sum of the individual NLEs evaluates to -13.6665, which is the value produced by the deterministic algorithm as well. This agrees with the analysis done in section 3. The individual NLEs for both stochastic versions of the Lorenz 63 system behave very similarly when the noise amplitude is increased. In both cases, the top two exponents decrease when the noise amplitude increases, whereas the bottom exponent increases. For the system with SALT, the bottom exponent increases faster than in the fluctuation-dissipation system. Thus the sum of the NLEs is independent of the noise, while the individual NLEs are not. Analytical results for the individual NLEs are difficult to establish and often only estimates are possible. We therefore restrict ourselves here to studying the value of the sum of the exponents, rather than the individual NLEs.

The Lorenz system with fluctuation-dissipation noise, with noise amplitude $\beta=0.5$, gives rise to the plots shown in figure 2. The NLEs converge and their sum for this particular realisation is -13.7636. This is a discrepancy in the first decimal place compared to both the deterministic and the SALT case, which agree.

From figures 1 and 2 it can be seen in the left panels that both systems give rise to a locus of phase space trajectories that vaguely resemble a butterfly. The realisation of the Wiener process is the same for both figures. It can be seen that the two types of noise give rise to different trajectories in phase space. The right panels show that for both systems the NLEs have converged. The constancy of the sum for SALT becomes especially clear in figure 3, which shows a hundred computations for increasing noise amplitude. Each time the noise amplitude is increased, a new realisation of the Wiener process is used. The solid black line is the sum of the NLEs for fluctuation-dissipation noise, which depends on each realisation, whereas the dashed line, which shows the sum of the NLEs for SALT, is constant at -13.6665 and independent of noise amplitude and realisation of the Wiener process. For a fixed realisation of the Wiener process, the theory predicts a linear behaviour. This can be seen in figure 4.

The individual Lyapunov exponents for systems (21) and (24) are almost identical for each realisation, although the sums are different, since the phase-space volume contraction rate is preserved for the system (21).

6 Conclusion

We used the Kelvin circulation theorem for the ideal Oberbeck-Boussinesq equations to modify the equation of motion and the advection equation for the heat so as to include stochastic advection by Lie transport (SALT). We then added viscosity in the motion equation and heat diffusivity in the advection equation to obtain the Oberbeck-Boussinesq equations with SALT. By using a specific truncated Fourier series expansion of these equations, the corresponding Lorenz 63 system with SALT was obtained in equation set (23). This low-dimensional system of stochastic differential equations was then compared with a Lorenz 63 system (24), which was perturbed using fluctuation-dissipation noise (linear multiplicative noise in each variable), which has been previously studied in [CSG11].

By applying methods from the theory of random dynamical systems, we were able to show that the two types of systems (23) and (24) have rather different qualitative properties. In particular the linear multiplicative noise in each variable of the fluctuation-dissipation system changes the average rate of contraction or expansion of phase space volume, relative to the deterministic system; whereas the system with SALT conserves this rate. This means that the type of noise introduced into low-dimensional dynamical systems can affect properties of the underlying deterministic system, so one should examine the effects of stochasticity on a qualitative level. For example, when a system of deterministic equations is Hamiltonian, introducing arbitrary noise would destroy the Hamiltonian structure by altering the average rate of phase space volume contraction, which for Hamiltonian systems is exactly zero. SALT would preserve this property. See [Hol15, ADCH15, HT16, Hol17, CGH17, ADCH18, AGH18, CHR18, HT18, CFH19].

For the numerical verification of our analytical results for the stochastic Lorenz systems, we introduced a stochastic generalisation of the Cayley method for the computation of numerical Lyapunov exponents (NLEs). This method is a QR-based algorithm in which the orthogonal matrix is determined via the Cayley transform. The method turned out to be robust and stable for the stochastic case. Improvements to the numerical calculations are possible by using better numerical methods. However, the numerical results we obtained here agree with our analytical results and the individual NLEs for the deterministic Lorenz system agree reasonably well with those found in literature. In the stochastically perturbed Lorenz systems, the numerical method is able to distinguish clearly between the two types of noise.

7 Outlook for a new climate science paradigm

In this paper we have investigated the effect of Stochastic Advection by Lie Tranport (SALT) on Rayleigh-Bénard convection and how the underlying Lorenz 63 system changes when the SALT approach is applied. In a celebrated unpublished paper [Lor95] Ed Lorenz defined the statistical approach to climate science by the adage, “Climate is what you expect, weather is what you get.” That is, climate science is fundamentally probabilistic. This realisation motivates the use of so-called Lagrangian Averaged SALT (LA-SALT) approach, first introduced in [DH18] and developed further in [DHL19] and [AOdLHT19]. Briefly put, the LA-SALT approach alters the stochastic material velocity $\widehat{\textbf{u}}$ of the circulation loop in Kelvin’s theorem for SALT in equation (10) by replacing the drift velocity u by its expected value $\mathbb{E}[\textbf{u}]$. Namely,

[TABLE]

where $\mathbb{E}[\,\cdot\,]$ denotes expectation. The altered velocity of the circulation loop in (58) allows one to take the expectation of the entire Kelvin circulation theorem. Consequently, the equations governing the expected solution comprise a deterministic subsystem of the stochastic Kelvin circulation theorem for the stochastic velocity vector field in (58). By following the approach illustrated in the present paper from equation (11) onward in the present paper, one may derive the LA-SALT version of the Lorenz 63 system. Perhaps not surprisingly, after replacing the transport drift velocity $X$ in the nonlinear terms in equations (22) by its expectation $\mathbb{E}[X]$ these equations become

[TABLE]

Taking the expectation of each of the equations in (59) then yields a deterministic dynamical system for the expected solutions,

[TABLE]

The system (60) is isomorphic to the deterministic Lorenz 63 system with a slightly renormalised dissipation. Hence this system has a strange attractor and one can make the following analogy to Lorenz’s quote. The Lorenz 63 system with SALT in (21) is like the weather. (It’s what you get.) and the Lorenz 63 system with LA-SALT in (60) is like the climate. (It’s what you expect.) Since the LA-SALT approach recovers the deterministic Lorenz 63 equations as a subsystem, all of the dynamical systems analysis which has been done for the deterministic Lorenz 63 system previously now carries over to the LA-SALT expectation equations in (60).

Finally, if one subtracts the deterministic expected LA-SALT equations in (60) from the stochastic LA-SALT equations in (59), one finds a system of linear stochastic equations driven by the expected solutions of the system (60) for the differences,

[TABLE]

regarded as stochastic fluctuations away from the Lorenz 63 system. Namely,

[TABLE]

We would like to dedicate future work to the further investigation of the systems in equations (59)–(62) and their implications for the differences between ‘what you expect’ in (60) and the dynamical evolution of the statistics of ‘what you get’ in (62). In particular, the accuracy of the prediction of ‘what you expect’ in the chaotic solutions on the strange attractor of the equation (60) depends on the amplitude of the noise representing ‘what you don’t know’ in (58). Moreover, the time-dependent variances $\mathbb{E}[(X^{\prime})^{2}]$, $\mathbb{E}[(Y^{\prime})^{2}]$ and $\mathbb{E}[(Z^{\prime})^{2}]$ of ‘what you get’ derived from the equation set (62) depend on the amplitude of the noise in (58) and they are driven by the chaotic evolution of ‘what you expect’ in (60). Even with the simplifications that (i) the expectations of the solutions decouple as a closed nonlinear chaotic system isomorphic to the Lorenz 63 system in (1), and (ii) the stochastic evolution of the fluctuations in (61) is linear, the dynamics of the statistics of the fluctuation will be intricate and challenging to characterise fully.

Acknowledgements

We are enormously grateful for many suggestions for improvements offered in encouraging discussions with M. Chekroun, J. de Cloet, D.O. Crisan, M. Engel, V. Lucarini, J.C. McWilliams, E. Mémin, M. Rasmussen, V. Resseguier and S. Takao. The work of DDH was partially supported by EPSRC Standard Grant [grant number EP/N023781/1] and EL was supported by [grant number EP/L016613/1]. EL is grateful for the warm hospitality at the Imperial College London EPSRC Centre for Doctoral Training in Mathematics of Planet Earth (mpecdt.org) during the course of this work.