Derivation of Delay Equation Climate Models Using the Mori-Zwanzig Formalism

TL;DR

This paper uses the Mori-Zwanzig formalism to systematically derive delay models for climate phenomena like ENSO, providing a more justified approach and resulting in models with improved accuracy in replicating observed periods.

Contribution

It introduces the Mori-Zwanzig formalism to derive delay climate models from PDEs, including methods to approximate unresolved dynamics, and demonstrates its application to ENSO modeling.

Findings

Derived a nonlinear delay model for ENSO with Mori-Zwanzig formalism

The new delay term increases the ENSO period, aligning better with observations

Provides a systematic framework for delay model derivation in climate science

Abstract

Models incorporating delay have been frequently used to understand climate variability phenomena, but often the delay is introduced through an ad-hoc physical reasoning, such as the propagation time of waves. In this paper, the Mori-Zwanzig formalism is introduced as a way to systematically derive delay models from systems of partial differential equations and hence provides a better justification for using these delay-type models. The Mori-Zwanzig technique gives a formal rewriting of the system using a projection onto a set of resolved variables, where the rewritten system contains a memory term. The computation of this memory term requires solving the orthogonal dynamics equation, which represents the unresolved dynamics. For nonlinear systems, it is often not possible to obtain an analytical solution to the orthogonal dynamics and an approximate solution needs to be found. Here, we demonstrate the Mori-Zwanzig technique for a two-strip model of the El Nino Southern Oscillation (ENSO) and explore methods to solve the orthogonal dynamics. The resulting nonlinear delay model contains an additional term compared to previously proposed ad-hoc conceptual models. This new term leads to a larger ENSO period, which is closer to that seen in observations.

Tables

Table 1. Table 1 : The ranges in which the nondimensional parameters are varied. Values of y n = 2 subscript 𝑦 𝑛 2 y_{n}=2 , A 0 = 0.2 subscript 𝐴 0 0.2 A_{0}=0.2 , and θ = 3 𝜃 3 \theta=3 were used in Subsection 5.1 .

Parameter		Dimensional	Dimensionless	Step
Wind forcing factor at $y_n$	$\theta$	-	1.0 - 4.0	0.2
Wind forcing strength	$A_0$	0.5 - 3.0 $\cdot {10}^{-2}$ Pa	0.1 - 0.6	0.05
Latitude Rossby waves	$y_n$	${5.0}^{\circ }$ - ${12.1}^{\circ }$	1.4 - 3.4	0.2

Table 2. Table 2 : Parameter values used in determining c T subscript 𝑐 𝑇 c_{T} and c h ∗ superscript subscript 𝑐 ℎ c_{h}^{*} .

Damping scale Newtonian cooling	${ϵ}_T$	$9.25\cdot {10}^{-8}$ s-1
Basin length	$L$	$1.5\cdot {10}^7$ m
Velocity first baroclinic Kelvin mode	$c_0$	2 m/s
Background wind forcing strength	${\tau }_0$	$2.667\cdot {10}^{-7}$ m/s2
Parametrization constant	$b_w$	$1.026\cdot {10}^2$ s
Depth surface layer	$H_1$	50 m
Depth top layer	$H$	200 m
Depth for temperature gradient	$\tilde{H}$	50 m
Steepness transition subsurface temperature	$H^{\ast }$	30 m
Temperature without dynamics	$T_0$	30 ∘C
Background subsurface temperature	$T_{s0}$	22 ∘C
Rayleigh friction coefficient	$a_M$	$1.3\cdot {10}^{-8}$ s-1
Scaling parameter	$ϵ$	${10}^{-4}$
Reference point in east of basin	$x_E$	0.9

Equations

$$\frac{\mathrm{d}T}{\mathrm{d}t}=T(t)-T(t{)}^3-\alpha T(t-\delta ).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{array}{c}\hat{\varphi }\\ \tilde{\varphi }\end{array}\right)=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)\left(\begin{array}{c}\hat{\varphi }\\ \tilde{\varphi }\end{array}\right)$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\hat{\varphi }(t)=A_{11}\hat{\varphi }(t)+A_{12}{\mathrm{e}}^{A_{22}t}\tilde{x}+{\int }_0^t{A}_{12}{\mathrm{e}}^{A_{22}(t-s)}{A}_{21}\hat{\varphi }(s)\mathrm{d}s.$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\varphi (t)=R(\varphi (t)),\varphi (0)=x,$$

$$\frac{\partial }{\partial t}u(x,t)=\mathcal{L}u(x,t),u(x,0)=g(x),$$

$$\frac{\partial }{\partial t}[{\mathrm{e}}^{t\mathcal{L}}g](x)=[{\mathrm{e}}^{t\mathcal{L}}\mathcal{L}g](x)=[{\mathrm{e}}^{t\mathcal{L}}P\mathcal{L}g](x)+[{\mathrm{e}}^{t\mathcal{L}}Q\mathcal{L}g](x).$$

$${\mathrm{e}}^{t(A+B)}={\mathrm{e}}^{tA}+{\int }_0^t{\mathrm{e}}^{(t-s)(A+B)}B{\mathrm{e}}^{sA}\mathrm{d}s,$$

$$[{\mathrm{e}}^{t\mathcal{L}}Q\mathcal{L}g](x)=[{\mathrm{e}}^{tQ\mathcal{L}}Q\mathcal{L}g](x)+{\int }_0^t[{\mathrm{e}}^{(t-s)\mathcal{L}}P\mathcal{L}{\mathrm{e}}^{sQ\mathcal{L}}Q\mathcal{L}g](x)\mathrm{d}s.$$

$$\frac{\partial }{\partial t}[{\mathrm{e}}^{t\mathcal{L}}g](x)=[{\mathrm{e}}^{t\mathcal{L}}P\mathcal{L}g](x)+[{\mathrm{e}}^{tQ\mathcal{L}}Q\mathcal{L}g](x)+{\int }_0^t[{\mathrm{e}}^{(t-s)\mathcal{L}}P\mathcal{L}{\mathrm{e}}^{sQ\mathcal{L}}Q\mathcal{L}g](x)\mathrm{d}s.$$

$$\frac{\partial }{\partial t}{\varphi }_i(x,t)=R_i(\hat{\varphi }(x,t))+F_i(x,t)+{\int }_0^t{K}_i(\hat{\varphi }(x,t-s),s)\mathrm{d}s,$$

$$F_i(x,t)=[{\mathrm{e}}^{tQ\mathcal{L}}Q\mathcal{L}g](x),K_i(\hat{x},t)=[P\mathcal{L}{F}_i](x,t),$$

$$\frac{\partial }{\partial t}{F}_i(x,t)=Q\mathcal{L}{F}_i(x,t),F_i(x,0)=Q\mathcal{L}{x}_i.$$

$$({\partial }_t+{ϵ}_0)(h_e-h_n)+{\partial }_x{h}_e$$

$$({\partial }_t+{ϵ}_0)h_n-\frac{1}{y_n^2}{\partial }_x{h}_n$$

$${\partial }_t{T}_e+c_T{T}_e-c_h{h}_e$$

$$h_e(0,t)=r_W{h}_n(0,t),h_n(1,t)=r_E{h}_e(1,t).$$

$${\partial }_t{h}_c+{ϵ}_0{h}_c+{\partial }_x{h}_c$$

$${\partial }_t{h}_n+{ϵ}_0{h}_n-\frac{1}{y_n^2}{\partial }_x{h}_n$$

$$\displaystyle\partial_{t}T_{e}+c_{T}T_{e}-c_{h}\Big{(}h_{c}+\frac{1}{1+y_{n}^{2}}h_{n}\Big{)}$$

$$h_{c}(0,t)=\Big{(}r_{W}-\frac{1}{1+y_{n}^{2}}\Big{)}h_{n}(0,t),\qquad r_{E}h_{c}(1,t)=\Big{(}1-\frac{r_{E}}{1+y_{n}^{2}}\Big{)}h_{n}(1,t).$$

$$h_c^0(x,t)=H_c{\mathrm{e}}^{{\sigma }_{k}t}{\mathrm{e}}^{-({\sigma }_k+{ϵ}_0)x},h_n^0(x,t)=H_n{\mathrm{e}}^{{\sigma }_{k}t}{\mathrm{e}}^{({\sigma }_k+{ϵ}_0)y_n^{2}x},$$

$$\sigma_{k}=-\epsilon_{0}+\frac{1}{1+y_{n}^{2}}\Big{(}\ln\Big{(}\frac{r_{E}r_{W}(1+y_{n}^{2})-r_{E}}{(1+y_{n}^{2})-r_{E}}\Big{)}+2\pi ik\Big{)},\qquad k\in\mathbb{N}.$$

$$\begin{split}\mathcal{L}&=\Big{(}-\big{(}\epsilon_{0}+\partial_{x}\big{)}h_{c}(x,0)+\mu\Big{(}1-\frac{\theta}{1+y_{n}^{2}}\Big{)}g(x)T_{e}(x_{E},0)\Big{)}\partial_{h_{c}}\\ &\quad+\Big{(}-\big{(}\epsilon_{0}-\frac{1}{y_{n}^{2}}\partial_{x}\big{)}h_{n}(x,0)-\mu\frac{\theta}{y_{n}^{2}}g(x)T_{e}(x_{E},0)\Big{)}\partial_{h_{n}}\\ &\quad+\Big{(}-c_{T}(x)T_{e}(x,0)+c_{h}(x,\epsilon t)\big{(}h_{c}(x,0)+\frac{1}{1+y_{n}^{2}}h_{n}(x,0)\big{)}\Big{)}\partial_{T_{e}}.\end{split}$$

$$P(f(T_e,h_c,h_n))=f(T_e,0,0)=:\hat{f}(T_e),$$

$$[{\mathrm{e}}^{t\mathcal{L}}P\mathcal{L}{T}_e](x,0)=-c_T(x)T_e(x,t).$$

$$\begin{array}{rl}{\partial }_t{h}_c^Q(x,t)& =-({ϵ}_0+{\partial }_x)h_c^Q(x,t),\\ {\partial }_t{h}_n^Q(x,t)& =-({ϵ}_0-\frac{1}{y_n^2}{\partial }_x)h_n^Q(x,t),\end{array}$$

$$\begin{split}F_{T_{e}}(x,t)&=c_{h}(x,\epsilon t)\big{(}h_{c}^{Q}(x,t)+\frac{1}{1+y_{n}^{2}}h_{n}^{Q}(x,t)\big{)}\\ &=c_{h}(x,\epsilon t)\big{(}e^{-(\epsilon_{0}+\partial_{x})t}h_{c}(x,0)+\frac{1}{1+y_{n}^{2}}e^{-(\epsilon_{0}-\frac{1}{y_{n}^{2}}\partial_{x})t}h_{n}(x,0)\big{)}.\end{split}$$

$$\begin{split}K_{T_{e}}(T_{e}(x,0),t)&=c_{h}(x,\epsilon t)\Big{(}\mu\Big{(}1-\frac{\theta}{1+y_{n}^{2}}\Big{)}\mathrm{e}^{-(\epsilon_{0}+\partial_{x})t}-\mu\frac{\theta}{y_{n}^{2}}\frac{1}{1+y_{n}^{2}}\mathrm{e}^{-(\epsilon_{0}-\frac{1}{y_{n}^{2}}\partial_{x})t}\Big{)}\\ &\qquad\cdot g(x)T_{e}(x_{E},0).\end{split}$$

$$\frac{\mathrm{d}{T}_e}{\mathrm{d}t}(x,t)$$

$$\displaystyle\quad+c_{h}(x,\epsilon t)\Big{(}\mathrm{e}^{-(\epsilon_{0}+\partial_{x})t}h_{c}(x,0)+\frac{1}{1+y_{n}^{2}}\mathrm{e}^{-(\epsilon_{0}-\frac{1}{y_{n}^{2}}\partial_{x})t}h_{n}(x,0)\Big{)}$$

$$\displaystyle\quad+\int_{0}^{t}c_{h}(x,\epsilon(t-s))\Big{(}\mu\big{(}1-\frac{\theta}{1+y_{n}^{2}}\big{)}\mathrm{e}^{-(\epsilon_{0}+\partial_{x})(t-s)}$$

$$\displaystyle\quad\qquad-\mu\frac{\theta}{y_{n}^{2}}\frac{1}{1+y_{n}^{2}}\mathrm{e}^{-(\epsilon_{0}-\frac{1}{y_{n}^{2}}\partial_{x})(t-s)}\Big{)}g(x)T_{e}(x_{E},s)\mathrm{d}s.$$

$${\partial }_{t}f+c{\partial }_{x}f=-{ϵ}_{0}f,$$

$${\mathrm{e}}^{-({ϵ}_0+c{\partial }_x)(t-s)}f(x,s)={\mathrm{e}}^{-{ϵ}_0(t-s)}f(c(t-s)+x_s,s),$$

Full text

1Institute for Marine and Atmospheric Research Utrecht, Department of Physics, Utrecht University, Utrecht, The Netherlands

2Mathematics of Planet Earth Program, University of Reading, Reading, UK

3Department of Mathematics, University of Exeter, Exeter, UK

4CSIRO Oceans and Atmosphere, Hobart, TAS, AU

5Mathematical Institute, Utrecht University, Utrecht, The Netherlands

6Centre for Complex Systems Studies, Faculty of Science, Utrecht University, Utrecht, The Netherlands

Corresponding author: s.k.j.falkena@student.reading.ac.uk

Subject: Applied mathematics, Oceanography, Differential equations

Keywords: Delay models, Mori-Zwanzig, Reduction methods, El Niño Southern Oscillation, Conceptual models, Feedback effects

1 Introduction

To study climate variability and climate change, a hierarchy of models is currently used [5]. At the low end of this hierarchy are conceptual climate models, which contain only the necessary features for specific phenomena to occur and thus represent the dominant physical processes. These models are often used to study physical mechanisms in their purest forms, such as the causal chain behind a specific oscillation and how the period of this oscillation depends on the different processes involved. At the high end of the hierarchy, there are modern multi-process, multi-scale global climate models (GCMs) that aim to represent the total climate system in substantial detail. These models are used, for example, to make projections of future climate change as documented by assessments of the Intergovernmental Panel on Climate Change (IPCC).

One special class of conceptual models consists of differential delay models, which are particularly useful for compactly representing physical phenomena. Compared to ordinary differential equation models, delay models can potentially convey more information – they are infinite-dimensional dynamical systems – but can still be formulated in terms of functions of a single variable. This allows for an easier mathematical treatment than would a partial differential equation model, while such a delay model can still represent complex physical processes. A useful way in which these models can be analysed is through bifurcation analysis. Such analysis allows for the distinction of different dynamical regimes and the dependence on a few parameters can be investigated relatively easy compared to GCMs.

Keane et al. [18] provided an overview of delay models used to describe climate processes. The two main areas in the climate system upon which delay models so far have focused are Energy Balance Models and models for the El Niño Southern Oscillation (ENSO). There exist many positive and negative feedbacks in the climate system. Some of these feedbacks are delayed by a non-negligible amount of time, for example by transport through an ocean basin. The presence of such a delayed feedback can sometimes be determined from data [22]. In models incorporating delay it is not necessary to resolve all of the processes involved in the feedback. A parametrization by the resulting delay time is sufficient. Usually delay models for climate phenomena are derived from more complex models by making strong assumptions about the system, see e.g. [25, 1]. In many studies, the delay model is introduced in an ad-hoc way, usually through physical reasoning or semi-empirical indications. For instance, in many processes in climate, propagating waves play an important role and a wave-basin crossing time is used as a delay.

A prominent case of delayed feedback through wave propagation is the ENSO variability in the Tropical Pacific. During an El Niño event, the sea surface temperature in the eastern part of the basin is warmer than usual, resulting in severe weather disturbances in countries on both sides of the Pacific Ocean. Its counterpart is La Niña, when the sea surface temperature is colder than usual. These two events alternate with intermediate phases in between, resulting in an irregular oscillation with a period of four to seven years. One of the most successful models of ENSO is that by Zebiak and Cane [29]. The view of the behaviour of ENSO in terms of normal modes resulted in the so-called delayed-oscillator mechanism of ENSO [15, 16]. The delay mechanism here is related to the propagation of equatorial Kelvin and off-equatorial Rossby waves. These waves take time to travel through the basin, resulting in a delayed arrival of a temperature anomaly. Already in 1988, Suarez and Schopf proposed a delay model for the sea surface temperature $T$ in the eastern Pacific Ocean of the form [25]:

[TABLE]

Here $\delta$ is the delay time and $\alpha$ a parameter indicating the strength of the delayed feedback. Note that this model is scaled to contain as few parameters as possible. Other delay models based on the same mechanism have been proposed and studied, see e.g. [27, 19].

The delay model of Equation (1) gives oscillations with a period of two to three years for realistic values of $\alpha$ and $\delta$. This is on the short side with respect to the observed period of ENSO, indicating that some aspect is missing in the model. The physics behind the delay is well understood, justifying the linear terms in Equation (1). However, the nonlinear term in the model by Suarez and Schopf is proposed ad-hoc and no physical justification is given for it in their article. Battisti and Hirst provided some arguments for the form of the nonlinearity [1], but no thorough mathematical derivation has yet been provided.

The formal model reduction approach proposed by Mori [20] and Zwanzig [32] potentially provides a way to place the derivation of delay models upon a stronger mathematical foundation. Using the Mori-Zwanzig approach, one formally reduces the dimension of a system of ordinary differential equations by projecting the dynamics onto a select subset of resolved dependent variables [3, 11]. Closure is attained by replacing dependence on the unresolved variables by a memory integral and a term referred to as the noise term. This noise term only represents noise when the unresolved variables are chaotic on their own. Its statistics then are determined by the initial values of the unresolved variables.

A common application of the formalism is stochastic modelling [12], for example for systems with large scale differences, which in certain cases can be reduced to Markovian systems [28]. For Markov chains there are several approaches to model reduction [2]. Also for non-Markovian systems some results exist [4]. Hamiltonian systems are another class of systems to which the formalism has been applied [3, 30]. The component of the rewritten system that is focused on in this paper is the memory term. When the memory term is expressed as the convolution of a kernel function with the history of the system, it is in the form of a distributed delay. Under some approximations, this memory term can be simplified to a term with a discrete delay.

In this paper the Mori-Zwanzig formalism is applied to a spatially extended model of ENSO, which serves as a test case for the derivation of delay models using this formalism. In Section 2 the theory behind the Mori-Zwanzig formalism is shortly recapitulated. In Sections 3 and 4 the formalism is applied to linear and nonlinear two-strip models of ENSO and new delay models are derived. The new nonlinear models are analyzed in Section 5 and the results are summarized and discussed in Section 6.

2 Mori-Zwanzig Formalism

The Mori-Zwanzig formalism reformulates a set of ordinary differential equations (ODEs) into a reduced system for the resolved variables that still retains all of the dynamics of the original system. In this section, we outline the formalism. We adopt the notation and follow the approach of Chorin et al. [3] in the formulation of the reduction, which is based on the work by Mori (1965) [20] and Zwanzig (1973) [32]. The theory discussed in this section holds for ODEs, whereas the system we will consider is a partial differential equation (PDE). The formalism is not straightforward to apply to general systems of PDEs, which are infinite-dimensional. Problems arising when treating PDEs are discussed in Section 3.2.

We start by presenting a simple linear example with constant coefficients, to illustrate the idea behind the formalism. Consider the following linear system of ODEs for $\phi=(\hat{\phi},\tilde{\phi}):\mathbb{R}\rightarrow\mathbb{R}^{n}$ continuously differentiable:

[TABLE]

Here, we call $\hat{\phi}\in\mathbb{R}^{m}$ the resolved variables and $\tilde{\phi}\in\mathbb{R}^{n-m}$ the unresolved variables. We have $A_{11}\in\mathbb{R}^{m\times m}$, $A_{12}\in\mathbb{R}^{m\times(n-m)}$, $A_{21}\in\mathbb{R}^{(n-m)\times m}$ and $A_{22}\in\mathbb{R}^{(n-m)\times(n-m)}$ and initial conditions $\phi(0)=(\hat{x},\tilde{x})$. The goal is to derive an equation for the resolved variables $\hat{\phi}$ only. In this example, this can be done by solving the equation for $\tilde{\phi}$ using variation of constants and substituting the result into the equation for $\hat{\phi}$:

[TABLE]

The system in Equation (2) has been reduced to one equation for the resolved variables $\hat{\phi}$, with the only dependence on $\tilde{\phi}$ being the initial condition $\tilde{x}$. This reduced system is equivalent to the full system and exhibits the same behaviour. The first term on the right-hand side is referred to as the Markovian term, the second as the noise term (due to possible uncertainty in the initial condition) and the last as the memory term.

For the generalization of this idea to nonlinear systems, consider the following $n$-dimensional system of ODEs:

[TABLE]

where $\phi(t)\in\mathbb{R}^{n}$ is a continuously differentiable function of $t\in\mathbb{R}_{+}$, $x\in\mathbb{R}^{n}$ denotes the initial condition, and $R:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ has components $R_{i}$. To every initial condition $x$ there corresponds a trajectory $\phi(t)=\phi(x,t)$, $\phi:\mathbb{R}^{n}\times\mathbb{R}_{+}\rightarrow\mathbb{R}^{n}$, whose existence is assumed for all $t>0$. We consider the evolution of an observable $u(x,t):=g(\phi(x,t))$ along a solution of Equation (4), where $g$ is defined on $\mathbb{R}^{n}$. The quantity $u(x,t)$, $u:\mathbb{R}^{n}\times\mathbb{R}_{+}\rightarrow\mathbb{R}^{n}$, satisfies the PDE

[TABLE]

where $[\mathcal{L}u](x)=\sum_{i=1}^{n}R_{i}(x)\partial_{x_{i}}u(x)$ is the generator associated with vector field $R$ of Equation (4). This generator is called the Liouville operator [21].

The goal, as for the linear system, is to construct a system of equations for a select subset of $m$ resolved variables $\hat{\phi}\in\mathbb{R}^{m}$. As before, the unresolved variables are denoted by $\tilde{\phi}\in\mathbb{R}^{n-m}$, such that $\phi=(\hat{\phi},\tilde{\phi})$. To reduce the system from $n$ components to the desired $m$ components a projection operator $P:C(\mathbb{R}^{n},\mathbb{R}^{k})\rightarrow C(\mathbb{R}^{m},\mathbb{R}^{k})$ is needed. Here, $k$ is the dimension of an arbitrary function $f$ to which the projection is applied. Examples of projection operators include the conditional expectation (infinite-rank) [7] and the linear projection, defined by $[Pf](\hat{x})=f(\hat{x},0)=:\hat{f}(\hat{x})$, which sets all unresolved variables to zero and retains only the resolved components. The complement of $P$ is denoted $Q=I-P$, where $I$ is the identity operator.

We denote the solution of the linear PDE (5) as $u(x,t)=[\mathrm{e}^{t\mathcal{L}}g](x)$, where $\mathrm{e}^{t\mathcal{L}}$ is also referred to as the evolution operator [4]. In particular for $g(x)=x_{i}$ we find $\phi_{i}(x,t)=\mathrm{e}^{t\mathcal{L}}x_{i}$. Using this notation, Equation (5) can be written as

[TABLE]

Note that $\mathcal{L}$ and $\mathrm{e}^{t\mathcal{L}}$ commute. We now consider the second term in the right-hand side of this equation, $[\mathrm{e}^{t\mathcal{L}}Q\mathcal{L}g](x)$. This component gives the evolution of the unresolved variables. The Dyson formula [21],

[TABLE]

applied for $A=Q\mathcal{L}$ and $B=P\mathcal{L}$ gives

[TABLE]

Substitution into Equation (6) yields

[TABLE]

In particular for $g(x)=x_{i}$ we find $[\mathrm{e}^{t\mathcal{L}}g](x)=\phi_{i}(x,t)$ and Equation 9 becomes the generalized Langevin equation:

[TABLE]

where we use the shorthand notation $R_{i}(\hat{\phi}(x,t))=R_{i}([\hat{\phi}(x,t),0])=[PR_{i}](\phi(x,t))$ and

[TABLE]

with notation $K_{i}(\hat{x},t)=K_{i}([\hat{x},0],t)$. Note that $F_{i}(x,t)$ is the solution to the orthogonal dynamics equation:

[TABLE]

In general it is not known whether this initial value system is well-posed, but in specific cases (approximate) solutions can be obtained. The three terms on the right-hand side of the Langevin equation (10) are called the Markovian term $R_{i}(\hat{\phi}(x,t))$, the noise term $F_{i}(x,t)$ and the memory term, defined as the integral over the memory integrand $K_{i}(\hat{\phi}(x,t-s),s)$. This memory integrand consists of a memory kernel applied to the resolved variables. Note that the integral over $K_{i}(\hat{\phi}(x,s^{\prime}),t-s^{\prime})$ is equal to the memory term in Equation (10) after a change of variables $s^{\prime}=t-s$. From this point onward we will use the integral in $s^{\prime}$, dropping the prime for conciseness. For linear systems solving the orthogonal dynamics system simplifies to the case of Equation (2) and the Langevin equation is equivalent to the result obtained by using variation of constants.

The generalized Langevin system (10) is still exact, but it is not necessarily simpler. If solving the orthogonal dynamics equation (12) is as difficult as solving the full system, there is no use in applying the formalism. The applicability thus depends on the particular system and whether a suitable projection exists. Such a projection would yield an orthogonal dynamics system that is relatively straightforward to solve or approximate in a good way. Applications to slow-fast [28], Markovian [2], non-Markovian [4], and Hamiltonian systems [3], as well as systems with an orthogonal basis of eigenfunctions [26], have been considered in the literature [11]. For non-Hamiltonian systems of PDEs, which are considered here, less is known. The standard approach when considering PDE systems is to expand in different modes of the system [26]. More recently an expansion based on the Faber series has been proposed [31].

The main challenge when applying the Mori-Zwanzig formalism is the solution of the orthogonal dynamics system (12). The choice of projection operator $P$ is an important factor in determining the form and complexity of the orthogonal dynamics equation. The projection should be chosen such that the orthogonal dynamics system is stable, meaning you need to retain stabilizing factors in the unresolved dynamics, and less complex than the original system. Alternatively, one may approximate the orthogonal dynamics equation by a less complex system. An example is the pseudo-orthogonal dynamics approximation derived by Gouasmi et al. [13], which is discussed in Appendix A. This approximation is applied in Sections 3.2 and 4.3. Preferably, the orthogonal dynamics system decays at a faster rate than the full system. For linear constant coefficient systems this means that the largest eigenvalue of the orthogonal dynamics system is smaller than that of the full system. In that case one may justify neglecting the noise term $F_{i}(x,t)$ in the Langevin equation.

3 Linear ENSO Model

In this section we apply the Mori-Zwanzig formalism to a linear model of ENSO variability. The model we study is a system of PDEs in one space variable describing the dynamics on two strips, one at the equator and one at higher latitude, in the Pacific Ocean. The Mori-Zwanzig formalism and the use of characteristics reduce this model to a linear delay equation similar to the model by Suarez and Schopf [25].

3.1 Model Formulation

ENSO is a coupled ocean-atmosphere phenomenon, where the variations in the sea surface temperature (SST) induce wind stress anomalies which drive ocean circulation changes affecting the SST. The interaction between the wind stress and the ocean is described in the two-strip ocean model derived and studied by Jin [15, 16]. This two-strip model is derived from the dimensionless shallow water equations in a normalized equatorial basin [5]. Assuming a parabolic dependence of the thermocline on latitude near the equator simplifies the shallow water equations to a system of equations for the thermocline depth at the equator ($h_{e}$) and at some latitude $y_{n}$ between 5∘N and 15∘N ($h_{n}$). Since the model used is scaled, the basin is of length one with $x=0$ being the western boundary and $x=1$ the eastern boundary of the Pacific ocean at the equator. The details of this scaling can be found in Appendix B.

The thermocline depth at both latitudes responds to a wind forcing whose strength depends on the ocean temperature at the equator. To describe this coupling we use a simplified version of the Gill atmosphere model [5]. An equation for the SST perturbations $T_{e}$ at the equator, see e.g. [6], completes the system. The resulting two-strip model describing the dynamics of ENSO is [5]

[TABLE]

for $x\in[0,1]$. The boundary conditions are

[TABLE]

Here $\epsilon_{0}$ is a linear damping coefficient, $\mu$ a coupling coefficient for the wind forcing and $\theta$ an order one coefficient representing the difference in the effect of wind-stress between the equator and higher latitudes. In the equation for SST the coefficient $c_{T}$ represents local damping and $c_{h}$ represents the effect of thermocline depth on temperature through background upwelling. Both $c_{T}$ and $c_{h}$ can depend on time, space or any of the variables ($h_{c}$, $h_{n}$, $T_{e}$), which can result in nonlinearities. Furthermore, $r_{W}$ and $r_{E}$ are a measure of the mass flux at the western and eastern boundaries respectively. The function $g(x)$ gives the pattern of the wind forcing in the zonal direction. The wind forcing depends on the SST anomaly at $x=x_{E}$, in the east of the basin. Because the SST only feeds back into the thermocline equations for $x=x_{E}$, it is sufficient for solving the system to only consider $T_{e}$ at that location. The $\partial_{x}$-terms represent the advection of anomalies in the thermocline by Kelvin waves ($h_{e}$) and Rossby waves ($h_{n}$).

The left-hand side of Equations (13a) and (13b) can be decoupled by introducing a new variable $h_{c}=h_{e}-\frac{1}{1+y_{n}^{2}}h_{n}$ and keeping $h_{n}$. The new variable $h_{c}$ is dominated by the thermocline depth at the equator, but includes some influence of the higher latitudes as well. In the new variables the equations are

[TABLE]

with boundary conditions

[TABLE]

If $r_{E}=0$, then $h_{n}(1,t)=0$, meaning no reflection occurs at the eastern boundary.

The solution to the homogeneous equations for thermocline depth in the rewritten system of Equations (15), that is without wind forcing ($\mu=0$), can be expanded in eigenmodes for rates $\sigma_{k}$:

[TABLE]

with

[TABLE]

The boundary conditions (16) imply that H_{c}=\Big{(}r_{W}-\frac{1}{1+y_{n}^{2}}\Big{)}H_{n}, where $H_{n}$ is arbitrary. For $r_{E}=0$ the solution of the homogeneous system is trivial after one round trip: $h_{c}^{0}=h_{n}^{0}=0$. The solutions in Equation (17) are the eigensolutions of the two-strip model. Note that these eigensolutions are not orthogonal, meaning they are not convenient to use as a basis on which can be projected.

3.2 Mori-Zwanzig Formalism

Starting from the rewritten version of the two-strip model in Equation (15), the goal is to derive a delay equation describing ENSO using the Mori-Zwanzig formalism. The resulting model is expected to be similar to that by Suarez and Schopf [25] as given in Equation (1). The Mori-Zwanzig formalism works on ODEs, while the two-strip model is a system of PDEs. A route that is often taken in such situations is to expand the general solution in a basis of eigensolutions of the PDE [4, 11, 26]. Truncating this expansion and projecting along the adjoint eigensolutions reduces the PDE to a system of ODEs. However, for the two-strip model this is impractical, since the eigensolutions may be degenerate. Other projection methods, such as the use of Fourier exponentials or orthogonal polynomials as a basis [23], will converge very slowly (thus, requiring high truncation order) due to incompatible boundary conditions, making them unsuitable for analytical computation.

Here we treat $\partial_{x}$ as an operator and consider the system as ODEs on the space of continuous functions (in $x$). The result of the Mori-Zwanzig formalism will contain terms of the form $\mathrm{e}^{t\partial_{x}}$, which are solution operators of scalar wave equations. In this section we consider the linear version of the two-strip model, meaning the coefficients $c_{T}$ and $c_{h}$ are allowed to depend on space ($x$) and possibly time ($t$), but not on any of the variables $h_{c}$, $h_{n}$ or $T_{e}$. Specifically, we permit dependencies $c_{T}(x)$ and $c_{h}(x,\epsilon t)$, where $\epsilon$ represents the option for a slowly varying background effect of thermocline anomalies on SST. Note that by considering only linear terms, we cannot expect to find the cubic term in the model by Suarez and Schopf in Equation (1). Extension to a nonlinear model is discussed in Section 4.

The first step in applying the Mori-Zwanzig formalism is to identify the Liouville operator following its definition in Section 2:

[TABLE]

The second step is to choose a projection operator $P$. The model by Suarez and Schopf is a delay equation for the temperature at the equator in the east of the basin. Therefore, we choose the equatorial temperature $T_{e}$ as the resolved variable and use the linear projection as defined in Section 2,

[TABLE]

with $P:C(\mathbb{R}^{3},\mathbb{R})\rightarrow C(\mathbb{R},\mathbb{R})$, reducing the number of dependent variables from three to one. This projection is infinte-rank as $T_{e}$ is still a function of $x$. Applying the formalism thus results in an equation for only the resolved variable $T_{e}$.

Now we have the information needed to apply the formalism and find the Langevin equation (10) for the linear two-strip model. The different terms are computed and discussed separately. Firstly the Markovian term is computed:

[TABLE]

As expected, this is the right-hand side dependence on the resolved variable $T_{e}$ in Equation (15c), which gives the evolution of the resolved variable in the model.

To compute the noise and memory term we need to solve the orthogonal dynamics equation. One problem that can arise with infinite-rank projections is that it has not been generally proven that $\mathrm{e}^{tQ\mathcal{L}}$ is well-posed. The projection we consider has infinite rank, however since $P\mathcal{L}Q$ is bounded, the operator $\mathrm{e}^{tQ\mathcal{L}}$ can be bounded by $M\mathrm{e}^{\omega t}$ for some $\omega$ [30] (note that $g(x)$ in Equation (19) is empirically derived and will be bounded - see Section 3.4). This implies the contribution of the memory term is bounded and $\mathrm{e}^{tQ\mathcal{L}}$ is well-posed.

Because the operator $Q\mathcal{L}$ is linear in this case, we can use the method of Gouasmi et al. [13] to simplify the orthogonal dynamics system in Equation (12). For an explanation of this simplification the reader is referred to Appendix A. The resulting orthogonal dynamics are advection equations with complete boundary conditions. The equations for the thermocline depth are

[TABLE]

with boundary conditions (16) for $h_{c}^{Q}$ and $h_{n}^{Q}$. The equation for $T_{e}^{Q}$ decouples and is not needed for computation of the noise or memory term. Here $Q$ is used to denote the variables in the orthogonal dynamics equation. These two equations are independent of each other and have exponential solutions. The solutions are the same as those of the homogeneous two-strip model as given in Equation (17). Having solved the orthogonal dynamics equation, the noise term can be computed following Equation (11):

[TABLE]

Here the solutions of the orthogonal dynamics system (22) have been substituted to find the final expression.

The last component of the Langevin equation that needs to be computed is the memory term. To do this, first the memory integrand is computed following Equation (11):

[TABLE]

Here we exploit that $P$ and $\mathcal{L}$ commute with $\mathrm{e}^{-(\epsilon_{0}+\partial_{x})t}$ and $\mathrm{e}^{-(\epsilon_{0}-\frac{1}{y_{n}^{2}}\partial_{x})t}$. This can be verified by comparing with the result of applying variation of constants (Equation (3)). As expected $K_{T_{e}}$ is linear in $T_{e}$. The memory term in the Langevin equation is found by substituting the computed memory integrand into the integral.

The resulting Langevin equation is obtained by substituting the computed Markovian (21), noise (23) and memory term (24) into Equation (10). We find the following equation for the temperature at the equator:

[TABLE]

Note the change of variables discussed in Section 2 has been applied to arrive at the memory term (25c). The partial derivative to $x$ is still present in the exponential terms. In the noise term (25b) this is not an issue, since the terms are exactly the solutions to the homogeneous system as given in Equation (17). In the memory term (25c) it is less clear how to evaluate these exponential operators. In the next section these terms are simplified using the method of characteristics.

3.3 Evaluation along Characteristics

In this section we focus on the exponential $\partial_{x}$-terms in the memory integral (25c) using the method of characteristics [10]. Both components in the memory kernel are of the form $e^{-(\epsilon_{0}+c\partial_{x})(t-s)}f(x,s)$, for either $c=1$ or $c=-1/y_{n}^{2}$. This expression is the solution to the PDE

[TABLE]

with initial conditions given at $t=s$. The characteristic curves of this equation are $x-x_{0}=c(t-t_{0})$, along which $f(x,t)$ is constant apart from damping caused by the $\epsilon_{0}$-term. This gives

[TABLE]

where $x_{s}$ is the location at time $s$. This expression is valid as long as the argument of $f$ lies in the domain. Otherwise boundary effects need to be considered.

In the memory integral the two exponential terms have different corresponding characteristics. They represent the eastward traveling equatorial Kelvin waves for $c=1$ and the westward traveling Rossby waves for $c=-1/y_{n}^{2}$. Note that the Rossby waves take longer to cross the basin than the Kelvin waves as $y_{n}>1$. In Figure 1 the characteristics of the memory term are shown. The red line shows what happens to a signal emitted from $x=0.7$ at time zero until it arrives at the eastern boundary.

Since the domain of the two-strip model is bounded ($x\in[0,1]$), the effect of these boundaries on the temperature signal needs to be discussed. Following the boundary conditions in Equation (16), we can express the equation for $T_{e}$ at $x_{b}=0,1$ in terms of either $h_{n}(x_{b},t)$ for the westward Rossby wave, or $h_{c}(x_{b},t)$ for the eastward Kelvin wave. If we consider $r_{E}$ and $r_{W}$ in Equation (16) to be nonzero (i.e. allowing for energy transfer between strips at the boundaries) then reflection of the characteristics must be considered. The fraction between the incoming and outgoing wave, corrected for their respective effects on $T_{e}$ (Equation (15c)), determine the reflection coefficients for the western and eastern boundary which are given by, respectively,

[TABLE]

The characteristics of the system are used to get an expression for the memory term. The goal is to find a result for the temperature in the east of equatorial basin ($x=x_{E}$), where the model by Suarez and Schopf is defined [25]. This also is the location of the temperature on which the wind forcing depends. Looking at the signal at one location allows for the following of characteristics from a source given by $g(x)$ to that one location. The memory integral at $x=x_{E}$ is of the form

[TABLE]

The first term in the memory integral (25c) gives the waves that are traveling westward at $t=s$. Setting $x_{E}=1$ for the characteristics, these waves need a time $t=1-x$ to arrive at the eastern boundary starting from an initial location $x$. Assuming reflection takes place at the eastern and western boundary, the signal arrives a second time after $t=1-x+(y_{n}^{2}+1)$. The signal keeps reflecting through the basin, arriving at the eastern boundary after times $t_{k}=1-x+k(y_{n}^{2}+1)$ for $k=0,1,2,...$. At each reflection the wave loses energy by a factor $A_{rE}$ at the eastern boundary and a factor $A_{rW}$ at the western boundary.

The result of these reflections through the basin is found by combining Equation (27) and the above discussion in Equation (25c). The resulting expression for the first part of the memory term at the eastern boundary is

[TABLE]

where $K_{max}(t)=\lfloor\frac{t-1}{y_{n}^{2}+1}\rfloor$ for $t\geq 1$ is the number of reflections that have occurred by time $t$. Note that there are time intervals for which this term has no effect at the eastern boundary, since it only represents one of the two characteristics. To get this part of the memory integral in a form which shows more of the delay behaviour, a change of coordinates can be applied. Let $x=1+k(y_{n}^{2}+1)-(t-s)$, for which $\frac{\mathrm{d}x}{\mathrm{d}s}=1$. Changing coordinates from $s$ to $x$, yields a memory integral for $T_{e}^{E}(t-(1+k(y_{n}^{2}+1)-x))$. This shows that the memory term contains a component that depends on past states of the resolved variable $T_{e}^{E}$. How strong the effect is at a certain time depends on the function $g(x)$, which gives the spatial distribution of the wind forcing.

The result for the waves that first travel towards the western boundary is achieved in a similar way. These waves need times $t_{k}=y_{n}^{2}x+1+k(1+y_{n}^{2})$ for $k=0,1,2,...$ to arrive at the eastern boundary. Now $x=-\frac{1}{y_{n}^{2}}(1+k(y_{n}^{2}+1)-(t-s))$ is used for the change of coordinates. Going through the same steps as before, the total memory integral becomes

[TABLE]

This expression shows that there are multiple delays present in the two-strip model. The exact form of the delay (distributed or discrete) and delay times are determined by the spatial pattern of the wind forcing $g(x)$.

3.4 Delay Model

Using the results from the two previous sections, a delay equation for the evolution of the SST in the east of the basin can be obtained. We consider Equation (25) with the expression for the memory term in Equation (31) at $x=x_{E}$. The effect of the components of the memory term decreases for higher $k$ by energy loss at reflection. To simplify this expression we make two assumptions.

First, we assume that there is no reflection at the eastern boundary, meaning $r_{E}=0$ and thus $A_{rE}=0$. The only two components in the sum of the memory term (Equation (31)) that remain with this assumption are the components for $k=0$. As noted in Section 3.1, the homogeneous solution in that case is identically zero after finite time, such that the noise term (25b) vanishes.

The function $g(x)$ in the memory integral is unspecified until now. It determines the form of the memory kernel by setting the pattern of the wind forcing. More specifically, $g(x)$ indicates where the effect of the wind is strong and weak. Following Jin [15] the wind dominantly has an effect near the centre of the basin. Away from this location the effect is small, meaning that the wind forcing acts quite locally. We approximate this local effect of the wind forcing by $g(x)=A_{0}\delta_{x_{w}}(x)$, a delta function of height $A_{0}$ at $x=x_{w}$. This delta function leaves only the effect of that one location on the integral, such that Equation (25) considered at $x=x_{E}$ simplifies to

[TABLE]

This is a linear equation with discrete delay for the temperature. Because we assumed $c_{h}$ and $c_{T}$ to be independent of $h_{c}$, $h_{n}$ and $T_{e}$ (see 3.2) no nonlinearity is found.

Equation (32) does not yet resemble the model by Suarez and Schopf from Equation (1) in the linear terms. Instead of one, there are two delay times present. However, $1-x_{w}\ll 1+y_{n}^{2}x_{w}$ for realistic parameters (see Appendix B) indicating that this effect can be considered to be immediate. Thus it is assumed that $T_{e}^{E}(t-(1-x_{w}))\approx T_{e}^{E}(t)$. This approximation yields the final linear delay model:

[TABLE]

where

[TABLE]

This model (after rescaling) gives the linear part of the model by Suarez and Schopf in Equation (1) [25]. The delay is due to the propagation of Rossby waves caused by a wind forcing which depends on the temperature near the eastern boundary. These waves travel to the western boundary, where they reflect in the form of Kelvin waves. This delay model does not yet account for the nonlinearity in the model by Suarez and Schopf, without which no stable oscillation will occur in the model. To get a more realistic result a nonlinear version of the two-strip model is considered as the starting point for applying the Mori-Zwanzig formalism in the next section.

4 Nonlinear ENSO Model

In this section we start by deriving a nonlinear variation of the two-strip model describing the ENSO dynamics. This nonlinearity allows for a more realistic modelling of ENSO. The resulting nonlinear two-strip model is studied using the Mori-Zwanzig formalism and variation of constants.

4.1 Model Formulation

In the two-strip model (15), as considered in Section 3, the thermocline feedback coefficient $c_{h}$ in the temperature equation did not depend on any of the variables $T_{e}$, $h_{e}$ or $h_{n}$. We now introduce a realistic state dependence into $c_{h}$, resulting in a nonlinear PDE system. The expression for $c_{h}$ is given by [6]:

[TABLE]

where $f_{h}(x)$ is the background wind forcing and $T_{s}$ the subsurface temperature at the equator as a function of thermocline depth $h$. For the parametrization of the subsurface temperature the result by Hao et al. is used [14]:

[TABLE]

Here $T_{0}$ is the ocean equilibrium temperature in absence of dynamics, $h_{0}$ an offset value for the thermocline and $T_{s0}$ the temperature at $h=-h_{0}$. The parameter $H^{*}$ determines the steepness of the transition when $h$ passes through $-h_{0}$. The derivative of $T_{s}$ to $h$, which is needed to get an expression for $c_{h}$ (35), is proportional to 1-\Big{(}\frac{T_{s}-T_{s0}}{T_{0}-T_{s0}}\Big{)}^{2}.

We assume that $T_{e}$ is proportional to $T_{s}-T_{s0}$, meaning perturbations in the equatorial sea surface temperature are proportional to perturbations in the subsurface temperature. Using this assumption, Equation (35) has the form

[TABLE]

where $c_{se}$ is the proportionality constant. This will introduce a cubic nonlinearity in the temperature equation of the two-strip model, which was not present in the linear version considered in Section 3. To check the validity of the assumption, buoy data of the equatorial Pacific Ocean is considered. Since the model only contains feedback of the temperature in the east of the basin, it is sufficient to consider only buoys in the eastern part of the Pacific Ocean. There are ten locations in the equatorial Pacific where buoy data is available. To avoid coastal boundary layer effects the second most eastern buoy, which is located at 110 degrees west, is chosen. In Figure 2 anomalies of the SST versus those of the subsurface temperature (40 m) are shown for this buoy. The correlation between the two datasets is 0.83, indicating there is a strong relation. The slope between the two temperature anomalies is $c_{se}\approx 1$. We note that also at other depths correlation is strong. However, the slope between the temperature anomalies decreases slightly with depth.

Using Equation (37) the nonlinear two-strip model (in the form of Equation (15)) becomes

[TABLE]

where $c_{h}^{*}(x)=f_{h}(x)\frac{T_{0}-T_{s0}}{H^{*}}$ and \beta=\Big{(}\frac{1}{T_{0}-T_{s0}}\Big{)}^{2}. We emphasize that this result is in principle only valid in the eastern part of the basin, which is sufficient for the model presented here. At other locations the correlation between anomalies in SST and subsurface temperatures is quite strong as well, but the proportionality constant is different. In the following we discuss two ways to apply the Mori-Zwanzig formalism to the nonlinear system of Equation (38). The first method is based on variation of constants, the second on an approximation to the Mori-Zwanzig formalism.

4.2 Variation of Constants

Since the equations for $h_{c}$ and $h_{n}$ are still linear, we may apply the Mori-Zwanzig formalism in its elementary form (3) based on the variation of constants. This yields the following equation for $T_{e}$:

[TABLE]

The terms within the integrals are the same as those in the memory integral of Equation (25c). Therefore, we can use the evaluation along characteristics in Section 3.3 to rewrite them in the form of Equation (31).

Then, under the same assumptions as made in Section 3.3, that is no reflection at the eastern boundary and a localized wind forcing, we may apply equivalent simplifications to obtain a delay equation for $T_{e}^{E}$. Assuming, as before, that the short delay is instantaneous, yields a nonlinear delay model for the temperature in the east of the basin, which is an exact reduction of the nonlinear two-strip PDE (38):

[TABLE]

where

[TABLE]

This is a nonlinear equation with discrete delay for the temperature at the equator, including two cubic terms. The difference with the model by Suarez and Schopf in Equation (1) is a fourth term, which is proportional to $T_{e}^{E}(t)^{2}T_{e}^{E}(t-d)$. The effect of this additional term on the dynamics of the delay model is studied in Section 5. Before discussing the behaviour of this extended delay model, first the application of the Mori-Zwanzig formalism in the form of the Langevin equation (10) to the nonlinear model in Equation (38) is considered.

4.3 Mori-Zwanzig Formalism

The Mori-Zwanzig formalism is valid for both linear and nonlinear equations. The challenge when considering nonlinear equations arises in solving the orthogonal dynamics equation. We use the same linear projection as for the linear model (20). Here, it could be tempting to use a Taylor expansion to find a formal expression for the noise term. However, because the nonlinear operator $\mathcal{L}$ is unbounded, this does not apply [17]. An alternative is the Faber expansion [31], which yields a numerically computable solution where the expansion is truncated after a certain number of terms. However, we aim at finding an analytic expression for the resulting scalar equation for $T_{e}^{E}$ as found in the linear case (33). For this purpose the truncation approximation is not suitable. Another option is to approximate the orthogonal dynamics equation by the pseudo-orthogonal dynamics (POD) equation as derived by Gouasmi et al. [13]. The conditions for this approximation are not met for the nonlinear two-strip model, but it can be used as a first estimate.

Here we will apply this POD approximation to derive an approximate solution. This derivation is given in Appendix A. The resulting equation is simplified by using the method of characteristics, assuming no reflection takes place at the eastern boundary and considering a localized wind forcing, just as was done in Sections 3.3 and 3.4. Considering the short delay as being instantaneous, the resulting nonlinear delay equation for the temperature in the east of the basin is

[TABLE]

where

[TABLE]

This equation is almost the same as Equation (40), which was derived by applying variation of constants. The difference is the fourth term. Since Equation (40) is exact, this is due to errors introduced by the POD approximation.

5 Analysis of the Delay Models

Section 3 shows that the presence of a delay follows already from the analysis of the two-strip model at the linear level, resulting in Equation (32). The nonlinearity due to temperature dependence of the thermocline feedback enters the resulting delay model in different terms depending on the reduction method (Section 4). Equation (40) is derived by applying variation of constants to the thermocline equations and is therefore exact. Equation (42) is derived using the POD approximation of Gouasmi et al. [13] (Appendix A) and contains some approximation error. Both models contain an extra term compared to the model proposed by Suarez and Schopf as given in Equation (1). In the rest of this section we will refer to Equation (40) as the VoC model, Equation (42) as the MZ model and Equation (1) as the S&S model.

Before studying the qualitative behaviour of the delay models (40) and (42) in more detail, we scale temperature and time to reduce the number of parameters. In this section the sub- and superscripts of temperature are omitted, such that $T$ is written for $T_{e}^{E}$. Time is scaled by $\tilde{t}=(c_{S}^{*}-c_{T}(x_{E}))t$ and temperature by $\tilde{T}=\sqrt{\frac{\beta c_{S}^{*}}{c_{S}^{*}-c_{T}(x_{E})}}T$, where $\tilde{t}$ and $\tilde{T}$ are the scaled quantities. The scaled equation for the VoC model of Equation (40) (omitting the tildes for simplicity) is

[TABLE]

where

[TABLE]

For $\gamma=0$ Equation (44) reduces to the S&S model (1). Note that the scaled parameters do not depend on the strength of the nonlinearity in the thermocline feedback $\beta$. Only the temperature scale depends on $\beta$. Since $c_{S}^{*}>c_{T}(x_{E})$ all scaled parameters are positive and $\gamma<1$. The scaling for the MZ model (42), obtained by the approximate Mori-Zwanzig formalism, is the same, resulting in:

[TABLE]

Figure 3 shows time profiles for the three different nonlinear models for realistic parameters. The profile of the MZ model has a very different shape compared to the other two models. Both newly derived models have a longer period than the S&S model, as well as a smaller amplitude. These two aspects are closer to measurements of ENSO. The main focus in this section is on the similarities and differences between the S&S model by Suarez and Schopf and the VoC model derived by applying variation of constants.

5.1 Bifurcation Analysis

Following the analysis by Suarez and Schopf [25], we first determine the steady states of Equations (44) and (46). Depending on the value of $\alpha$, the models have either one or three equilibria:

[TABLE]

The trivial equilibrium, $T_{00}=0$, undergoes a pitchfork bifurcation at $\alpha=1$, with two equilibria emerging for $\alpha<1$. For $\alpha<1$, $T_{00}$ is always unstable, while for $\alpha>1$ the stability depends on the delay.

Using the MATLAB toolbox DDE-BIFTOOL [8, 9, 24] we compute two-parameter bifurcation diagrams in $\delta$ and $\alpha$, for each of the three nonlinear models. Figure 4 shows the result for the original S&S model (note that $\gamma=0$ makes the VoC and MZ models equivalent to the S&S model). For $\alpha>1$ and increasing $\delta$, $T_{00}$ undergoes a supercritical Hopf bifurcation with a family of stable periodic orbits emerging. This large-amplitude symmetric periodic orbit corresponds to the ENSO behaviour. This family of ENSO periodic orbits extends to the region $\alpha<1$, where it emerges from a connecting orbit. The area between the connecting orbit and the Hopf curve for $\alpha<1$ is a multistable region, with the attractors being the ENSO periodic orbit and two stable equilibria. We also show level curves for the period of the ENSO orbit throughout the $\alpha$-$\delta$-plane. The values of the delay $\delta$ and the period length have been converted to the dimensional quantities for better comparison to observations (see Appendix B for scaling values).

As shown in Section 4, the extra terms involving $\gamma$ in the VoC and MZ model appear directly through the derivation from the ENSO PDE model. Using the scaling in Appendix B, this would correspond to $\gamma=0.49$. In Figure 5 we show the effect of incorporating non-zero $\gamma$ in the S&S model through the same two-parameter bifurcation diagrams. We observe that in the VoC method (which is an exact reduction of the PDE), the region of parameter space for $\alpha<1$ where stable oscillations exist is slightly reduced. In the region between the Hopf curve and the connecting orbit curve in Figures 4 and 5 the ENSO oscillation coexists with two non-zero stable equilibria. However, the length of the period has increased (shown by the shift in the level curves), as already suggested by Figure 3.

We can also get an idea of the dynamical consequences of the error in the POD approximation of Gouasmi et al. [13] on the nonlinear ENSO model. As mentioned in Section 2, the Mori-Zwanzig method is highly dependent on the ability to solve, or approximate to a good degree, the orthogonal dynamics equation. The POD approximation of Gouasmi et al. [13] (outlined in Appendix A) is applied in Section 4.3 to derive the MZ delay model. We see through Figure 5 that this approximation introduces an error. The error causes an underestimation of the parameter region of stable oscillations for $\alpha<1$, and an overestimation of the period of oscillation. In this particular case where one can find the exact projected equation through variation of constants, it is clear that the POD approximation is not sufficient and a different approximation should be proposed. This is outside the scope of this study and should be considered in future work.

5.2 Dependence of Period on Physical Parameters

We now turn to the scaling and nondimensionalization process to analyze the dependence on variable physical parameters in some detail. Only the VoC model is considered when comparing with the S&S model, since the MZ model does not show periodic behaviour for most realistic parameter values. The main uncertainties are in the latitude at which the Rossby waves travel ($y_{n}$), the wind forcing factor at $y_{n}$ ($\theta$), and the overall strength of the wind forcing ($A_{0}$). In Section 5.1 we used the values $y_{n}=2$, $A_{0}=0.2$, and $\theta=3$ for the reference scaling. In Figures 6 and 7 these three parameters are varied according to Table 1 to study the dependence of the period. Note that these values differ from the those in Jin [15] due to a different scaling. All other parameter values are kept constant and listed in Appendix B.

Varying the parameters in the ranges given in Table 1, the period of oscillation is computed for both the VoC model and the S&S model. Note that the considered points can lie in the multistable region between the connecting orbit and Hopf curve. Figure 6 compares the period of oscillation of both models. The period of the VoC model is larger whenever an oscillation is present. There are locations in parameter space where no oscillation occurs in the VoC model, where it does occur in the S&S model. This can also be seen by comparing Figures 4 and 5(a), where both the Hopf curve and connecting orbit curve have shifted upward for the VoC model. This upward shift reduces the region in parameter space in which periodic orbits occur.

The dependence of the period on the different parameters is shown in Figure 7. The different points for each parameter value are due to the variation of the other two parameters. For increasing $\theta$, that is, when the effect of the wind forcing at higher latitudes increases, the period slowly increases. There turns out to be a minimum value for $\theta$ around 1.7 below which no oscillations occur. In that case the signal at higher latitudes is too weak to have a significant effect at the eastern boundary. Considering the strength of the wind forcing $A_{0}$, a stronger wind results in a shorter period. This could be due to the larger absolute difference between the effect at the equator and at higher latitudes, leading to a weaker effect of the latter. This decrease with increasing wind strength appears to be approximately exponential. For a realistic $A_{0}$ in the centre of the range, the period of the oscillation is approximately 2.5 to 3.5 years. This is still smaller than that of ENSO.

Looking at the latitude $y_{n}$ at which the Rossby waves travel, instead of the latitude itself, $1/y_{n}^{2}$ is plotted, since this gives the velocity of the Rossby wave traveling at that latitude. For higher velocities, so lower latitudes, the oscillations have a smaller period. The faster the wave travels, the shorter the delay is, resulting in a shorter period. Similarly, slow waves result in longer periods.

6 Summary, Discussion and Conclusion

Delay models are useful as conceptual climate models due to their infinite dimensional nature and reduced number of parameters. They are suited for mathematical analysis and therefore can add to the physical understanding of the processes involved. In this paper the Mori-Zwanzig formalism has been investigated as a method to derive delay equations. The two-strip model of the El Niño Southern Oscillation (ENSO) was used as a test case for the application of the technique. The reason being that models incorporating delay have already been proposed for this phenomenon.

The Mori-Zwanzig formalism gives an exact rewriting of a system of ordinary differential equations [3]. The rewritten equation contains a Markovian, noise, and memory term. Here the focus was on the memory term, since this is an integral over the history of the system, just as a delay term represents this history. The memory integral can be rewritten in the form of a distributed delay term by making several approximations. Under some further approximations this can be reduced to a discrete delay. The memory kernel determines the times at which the integral is large and this way gives the delay time. Because of the integral, a peak in the memory kernel does not necessarily coincide with the dominant delay time.

For the model of ENSO, in addition to the use of the Mori-Zwanzig formalism, also the method of variation of constants has been employed. Starting from the linear two-strip model [15], the linear part of the delay model by Suarez and Schopf has been derived [25]. In this case the Mori-Zwanzig formalism is equivalent to using variation of constants (discussed in Section 2). Furthermore, a nonlinear version of the two-strip model was derived by assuming that the sea surface temperature is proportional to the subsurface temperature. From this nonlinear two-strip model (nonlinear) delay models have been derived. Here the two methods do not yield the same result. This is due to the approximations needed to obtain a closed-form equation from the Mori-Zwanzig formalism; whereas the method of variation of constants is exact.

Both derived nonlinear delay models contain an extra cubic delay term compared to the model by Suarez and Schopf [25]. In both cases this additional term results in an increased period of the model oscillation. Another consequence of the additional term is the decrease of the area in parameter space where stable oscillations occur. For the model derived using the Mori-Zwanzig formalism, this decrease is so large that no stable periodic behaviour occurs for most realistic values of the parameters. The exact derived model does show oscillations for these parameter values. The period of this model derived using variation of constants is closer to the real period of ENSO than the model proposed by Suarez and Schopf. However, its period is still smaller than what is seen in observations. One approach to improve the match between model period and data, could be to no longer assume a delta-function for the spatial pattern of the wind forcing, but rather take a more realistic pattern. As a consequence, the resulting delay model will no longer contain a discrete delay, but a distributed delay. Another option for possible improvement is to incorporate additional nonlinearities, for example in the thermocline equations. For the ENSO model it is not necessary to use the Mori-Zwanzig formalism to arrive at a delay equation. It does not give additional understanding compared to the method of variation of constants. This holds true for all models that are linear in the unresolved variables. This method yields a delay equation for the resolved variable (temperature) in both the linear and nonlinear case, because the thermocline depth equations are linear in both models. The Mori-Zwanzig formalism gives the same result for the linear system, but not for nonlinear systems. There the formalism does not give accurate results, since approximations are needed.

When the equations of the model considered are also nonlinear in the unresolved variables, the Mori-Zwanzig formalism is the only method that will give a result. In such a nonlinear case the orthogonal dynamics system has to be approximated. This approximation needs to be an improvement on the pseudo-orthogonal dynamics approximation [13], since this approximation was shown to be inaccurate for the ENSO model. It does not yield accurate results because the time scale of the unresolved variables is of the same order as that of the resolved variables. Only if its accuracy can be shown for a specific model, it can be justified to apply this approximation method. Additionally, a comment should be made on the computational cost of the pseudo-orthogonal dynamics approximation. If the corresponding ODE system derived for the pseudo-orthogonal dynamics is of high dimension, then solving such a system may not prove to be an improvement over solving the original full ODE model. However, with regard to the justification of the use of DDE models for parameter studies, one is only interested in whether or not the resulting memory kernels have pronounced peaks. Hence, the pseudo-orthogonal dynamics approximation may only need to be solved once. We expect this approach to be fruitful in climate models whenever the unresolved subsystem describes approximately linear wave transport. The derivation of improved approximations is a first step that needs to be taken to apply the Mori-Zwanzig formalism accurately to nonlinear models. This is a necessary step to be able to reliably derive nonlinear delay models for climate models.

Since many climate models are wave equations in one form or another, it is expected that they can also be represented by a delay system. This would imply that there is an abundance of phenomena in the climate system that can be described by a delay equation, possibly enabling the understanding of complex behaviour using relatively simple conceptual models.

Data access: The codes supporting this article have been uploaded as part of the supplementary material. They can also be found on the online repository figshare: https://doi.org/10.6084/m9.figshare.8085683.v1

Authors contributions: SF, CQ and HD designed the study. The work was carried out mainly by SF, with the exception of the bifurcation analysis in Section 5.1 which was done by CQ. All authors contributed to the work, discussed the results and read and approved the manuscript.

Competing Interests: The author(s) declare that they have no competing interests .

Funding: This work was supported by funding from the European Union Horizon 2020 research and innovation programme for the ITN CRITICS under Grant Agreement Number 643073 (CQ, JS and HD), the Mathematics of Planet Earth program (project No. 657.014.006) of the Dutch Science Foundation (NWO, JF) and EPSRC Grants No. EP/N023544/1 and No. EP/N014391/1 (JS).

Acknowledgments: We thank the two anonymous reviewers for their constructive comments. SF thanks the University of Exeter for hosting her for 5 months in 2018.

Appendix A Pseudo-Orthogonal Dynamics Approximation

A.1 Formulation

To simplify the issue of solving the orthogonal dynamics system, Gouasmi et al. [13] derived the pseudo-orthogonal dynamics (POD) equation. Under certain assumptions this is an exact rewriting of the orthogonal dynamics equation. In this rewritten form the orthogonal dynamics system can be more easily solved. The main assumption in the approach is the commutativity of $\mathrm{e}^{tQ\mathcal{L}}$ and $R$:

[TABLE]

For linear systems commutativity holds and this relation is exact, as shown in [13]. For nonlinear systems the above relation may be used as an approximation. The accuracy of this approximation is not a priori clear and requires verification. With the assumption of commutativity, the orthogonal dynamics equation can be reformulated into the POD equation:

[TABLE]

Note that this equation can be implemented more straightforward in numerical codes compared to the original orthogonal dynamics equation. When this equation is solved, the noise term corresponding to the resolved component $\phi_{i}$ is given by

[TABLE]

This is the part of the POD system corresponding to the respective resolved variable. The noise term thus can be retrieved directly when solving the POD equation. Defining $R_{Q}(x)=R(x)-[PR](x)$, the error made in this approximation is:

[TABLE]

The POD approximation simplifies solving the orthogonal dynamics system. However having a solution to this system is not necessarily sufficient to also formulate an expression for the memory term. When the noise term is a complicated function of the solution to the orthogonal dynamics equation, it can be quite difficult to get an expression for the memory kernel. Therefore, it is useful to look into approximations of the kernel. The following approximation has been derived by Gouasmi et al. [13]. The first step is to consider the $n$ components of $\mathcal{L}F_{i}(x,t)$ as the partial derivative of $F_{i}(x,t)$ in the direction of $\bar{R}(x)=R(x)/||R(x)||$, instead of $n$ separate derivatives in directions $R_{j}(x)$ for each coordinate $x_{j}$ with $j=1,\dots,n$. Here $\bar{R}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ and $||\cdot||$ is the $l^{2}$-norm, i.e. the standard Euclidean norm on $\mathbb{R}^{n}$. This yields

[TABLE]

Here the right-hand side is written as the limit corresponding to the derivative of $F_{i}$ at $x$ in the direction of $\bar{R}(x)$. The memory kernel $P\mathcal{L}F_{i}$ is the projection of this equation onto the resolved variables. Since $F_{i}(x,t)$ is the solution to the orthogonal dynamics equation, it only depends on the unresolved variables. Therefore, the second term in the numerator of (52) disappears after projection onto the resolved variables, i.e. $PF_{i}(x,t)=F_{i}(\hat{x},t)=0$. Applying $P$ to Equation (52) and using finite differences, results in the following approximation for the memory kernel:

[TABLE]

The exact result can be recovered if the limit $\epsilon\rightarrow 0$ exists. The memory integral can be approximated by the rectangle rule or another approximation method for integrals.

The derived approximation is most useful when working numerically. To obtain analytical expressions using this method can be cumbersome. However, if the result of the approximation remains tractable and the limit can be computed, it can result in an exact expression for the memory kernel.

Gouasmi et al. [13] applied this approximation to numerically compute the memory for the POD system to the Burger’s equation and the Kuramoto-Sivashinsky equation. They found that the approximation in the Burger’s equation was accurate for high wave numbers, but not for low wave numbers. This means it gives a good approximation when the memory is determined mainly by fast dynamics and thus is relatively short. However, for systems where slow dynamics affects the memory term, resulting in a longer memory, it appears to be less accurate.

A.2 Application to Nonlinear ENSO Model

Since it is not feasible to analytically solve the orthogonal dynamics equation for the nonlinear ENSO model as discussed in Section 4, we apply the POD approximation to find an approximate solution. The conditions for this approximation are not met for the nonlinear two-strip model, but it can be used as a first estimate.

The POD equations for the nonlinear two-strip model of equation (13) are

[TABLE]

The first two equations have exponential functions as solutions. To find a solution for $T_{e}^{Q}$ we substitute the solutions for $h_{c}^{Q}$ and $h_{n}^{Q}$ into the equation for $T_{e}^{Q}$. The solution for $T_{e}^{Q}$ with the condition that $(T_{e}^{Q})^{2}<\frac{1}{\beta}$ is

[TABLE]

If $(T_{e}^{Q})^{2}>\frac{1}{\beta}$ the $\tanh$ has to be replaced by a $\coth$, and when $(T_{e}^{Q})^{2}=\frac{1}{\beta}$ the result is a constant $T_{e}^{Q}$, since then $\partial_{t}T_{e}^{Q}(x,t)=0$. The initial conditions determine which of the solutions should be used. Most likely is that $(T_{e}^{Q})^{2}<\frac{1}{\beta}$, as $\beta$ is small and $T_{e}$ is of order one. Therefore, in the following Equation (55) is used.

The noise term is given by the right-hand side of the equation for $T_{e}^{Q}$ in Equation (54), for which now a closed expression is known. The next step is to compute the memory kernel. By the presence of the hyperbolic tangent and several nonlinearities, Equation (53) is used to approximate the memory kernel. Taking the limit $\epsilon\rightarrow 0$, yields the memory kernel. The equation for $T_{e}$ in the POD approximation becomes

[TABLE]

Similar to the procedure followed for the linear model in Sections 3.3 and 3.4, this equation can be simplified. The desired result is an equation for the temperature in the east of the basin. Applying the method of characteristics one gets rid of the exponential $\partial_{x}$-terms. Assuming no reflection takes place at the eastern boundary, the noise term vanishes. The memory term further simplifies to two delay terms when a localized wind forcing is assumed. Considering the short delay as being instantaneous, the resulting nonlinear delay equation is

[TABLE]

where

[TABLE]

Appendix B Scaling ENSO Model

In this appendix we give an overview of the parameters and scaling used in the final delay model of Equation (44), which is repeated here:

[TABLE]

where

[TABLE]

To arrive at this equation time and temperature are scaled by:

[TABLE]

Note that in Section 5 we drop the tildes after scaling.

The parameters defining $\alpha$, $\gamma$ and $\delta$ are:

[TABLE]

where

[TABLE]

with background wind forcing

[TABLE]

for $x_{0}=0.57$ and parameters

[TABLE]

The values of the involved parameters above are given in Table 2. Here the eastern boundary is not considered at $x=1$, but at $x=x_{E}$ to avoid boundary effects of the model. The other dimensionless variables in the definitions of $c_{S}^{*}$ and $c_{L}^{*}$ are $\mu$, $\epsilon_{0}$, $A_{rW}$, $\theta$, $y_{n}$, $A_{0}$ and $x_{w}$. Based on the book Nonlinear Physical Oceanography by Dijkstra [5] , we set

[TABLE]

with $r_{W}=3/5$. The other parameters are set as $\theta=3$, $y_{n}=2$, and $A_{0}=0.2$ in Section 5.1 and varied as discussed in Section 5.2 to study the dependence of the period.

We also note the nondimensionalisation of the shallow water equations used in [5]:

[TABLE]

where $\rho$ the density, $\beta$ the beta-plane parameter and the tildes denote the dimensional quantities. This scaling is also needed for the dimensionalisation of the parameter values in Table 1.