A Delay Equation Model for the Atlantic Multidecadal Oscillation
Swinda K.J. Falkena, Courtney Quinn, Jan Sieber, Henk A. Dijkstra
AA Delay Equation Model for the Atlantic MultidecadalOscillation
Swinda K.J. Falkena , , Courtney Quinn , , Jan Sieber , and Henk A. Dijkstra , Institute for Marine and Atmospheric Research Utrecht, Department of Physics, Utrecht University,Utrecht, The Netherlands CSIRO Oceans and Atmosphere, Hobart, TAS, AU Department of Mathematics, University of Exeter, Exeter, UK Department of Mathematics and Statistics, University of Reading, Reading, UK Centre for Complex Systems Studies, Faculty of Science, Utrecht University, Utrecht, The Nether-lands
Corresponding author : [email protected]
Abstract
A new technique to derive delay models from systems of partial differential equations, basedon the Mori-Zwanzig formalism, is used to derive a delay difference equation model for theAtlantic Multidecadal Oscillation. The Mori-Zwanzig formalism gives a rewriting of the originalsystem of equations which contains a memory term. This memory term can be related to a delayterm in a resulting delay equation. Here the technique is applied to an idealized, but spatiallyextended, model of the Atlantic Multidecadal Oscillation. The resulting delay difference modelis of a different type than the delay differential model which has been used to describe the ElNi˜no- Southern Oscillation. In addition to this model, which can also be obtained by integrationalong characteristics, error terms for a smoothing approximation of the model have been derivedfrom the Mori-Zwanzig formalism. Our new method of deriving delay models from spatiallyextended models has a large potential to use delay models to study a range of climate variabilityphenomena. a r X i v : . [ m a t h . D S ] A ug Introduction
To better understand climate variability and climate change often conceptual climate models are used.These models capture the dominant physical processes behind climate phenomena, allowing for animproved understanding. Delay equation models form one class of conceptual climate models. Thesetype of models are infinite dimensional, but often depend only on a few variables and parameters. Thismeans that they can potentially describe more complex behaviour compared to ordinary differentialequation (ODE) models, while still being easier to study than multi-dimensional partial differentialequation (PDE) models.Delay models have already been used to describe certain climate phenomena, particularly for theEl Ni˜no-Southern Oscillation (ENSO) and Earth’s Energy Balance [23]. Recently a new method ofderiving delay equation models has been proposed [18], allowing for a potential extension of the useof delay models to study other climate phenomena. This method of deriving delay models is basedon the Mori-Zwanzig (MZ) formalism, which allows for the reduction of high-dimensional systems toreduced-order models [4]. These reduced-order models are simpler to study, while still describing thephysical processes present in the original high-dimensional model. So far the method in [18] has onlybeen applied to a PDE model of ENSO, for which ad-hoc delay models were already proposed [31].Here we apply the MZ formalism to a PDE model of the Atlantic Multidecadal Oscillation (AMO),to investigate whether this phenomena can be described by a delay model as well.
Year -0.6-0.4-0.200.20.40.60.8 A M O i nde x ( ° C ) Figure 1: The AMO index for the last 160 years. The index is computed as the devia-tions of the area-weighted average SST over the North Atlantic. In black the 12-monthly run-ning mean is shown. Index computed by Enfield et al.
The Mori-Zwanzig (MZ) formalism provides a way of reducing a high-dimensional model to a reduced-order, more tractable system. The formalism is based on the work by Mori [26] and Zwanzig [36] instatistical mechanics. It has been reformulated to be suitable for constructing reduced-order modelsfor systems of ODEs [5, 4, 8]. When studying applications in climate one often has to deal withPDEs, making the application of the formalism challenging. In this section we start with a generaloverview of the MZ formalism based on [4], followed by a more detailed discussion of the particularsof applying the formalism to PDEs.For a vector of state variables φ ( x, t ) ∈ R n which are continuously differentiable in t ∈ R + andinitial conditions x ∈ R n , we consider the system of ODEs defining the dynamics: dd t φ ( x, t ) = R ( φ ( x, t )) , φ ( x,
0) = x, (1)where R : R n → R n is the vector-valued function of the specific system with components R i . Nowconsider the evolution of an observable u ( x, t ) := g ( φ ( x, t )) along a trajectory φ : R n × R + → R n .This observable satisfies the PDE ∂∂t u ( x, t ) = L u ( x, t ) , u ( x,
0) = g ( x ) , (2)where L is the Liouville operator (or generator) [27] given by [ L u ]( x ) = n (cid:88) i =1 R i ( x ) ∂ x i u ( x ) . (3)Note that for a linear system where R is defined by a matrix A having elements A ij , the Liouvilleoperator reads [ L u ]( x ) = (cid:80) ni =1 (cid:80) nj =1 A ij x j ∂ x i u ( x ) .To arrive at a reduced-order model for the dynamics governing φ ( x, t ) , one needs to decide onthe resolved variables ˆ φ ∈ R m . In our illustration we take ˆ φ as a subset of components φ i for someindices i . We also make a choice for an appropriate projection operator P : C ( R n , R k ) → C ( R m , R k ) onto these variables. Examples of projection operators are the linear projection, setting all unresolvedvariables to zero, and the conditional expectation [4]. Let Q = I − P denote the complement of P (with I the identity operator). Furthermore, we use the notation [ P R i ]( φ ( x, t )) = R i ([ ˆ φ ( x, t ) , R i ( ˆ φ ( x, t )) . We consider the choice of setting unresolved variables to zero for our projection P ,such that, for an arbitrary observable g ∈ C ( R n , R k ) (with arbitrary k ≥ ), the projection P isdefined as [ P g ]( φ ( x, t )) := g ([ ˆ φ ( x, t ) , . For example, for a component R i of the right-hand side, [ P R i ]( φ ( x, t )) := R i ([ ˆ φ ( x, t ) , .Having chosen a set of resolved variables and a projection operator P , the reduced-order modelcorresponding to the full system (1) is given by the generalized Langevin equation (see Chorin et al. for its derivation [4]): ∂∂t φ i ( x, t ) = R i ([ ˆ φ ( x, t ) , F i ( x, t ) + (cid:90) t K i ([ ˆ φ ( x, t − s ) , , s )d s, (4)3here φ i ( x, t ) is one of the resolved variables. The functions F i and K i are defined as F i ( x, t ) = [ e tQ L Q L ]( x ) , K i (ˆ x, t ) = [ P L F i ](ˆ x, t ) , (5)where ˆ x denotes the resolved part of the initial conditions x . The terms on the right-hand side of theLangevin equation are often referred to as the Markovian term R i ( ˆ φ ( x, t )) , the noise term F i ( x, t ) and the memory term, being the integral over K i ( ˆ φ ( x, t − s ) , s ) . For a linear system this memoryintegrand is obtained by applying a memory kernel to the resolved variables, i.e. K i ( ˆ φ ( x, t − s ) , s ) =ˆ K i ( s )[ ˆ φ ( x, t − s )] .The main difficulty in the application of the MZ formalism is calculating the terms F i , in which thenoise and the memory term enter, and which are the solutions of the orthogonal dynamics equations: ∂∂t F i ( x, t ) = Q L F i ( x, t ) , F i ( x,
0) = Q L x i . (6)In general it is not known if the system (6) is well-posed. However, for specific cases it is possible tofind approximate solutions. The possibility and difficulty of finding these solutions strongly dependson the choice of the resolved variables and projection operator. A suitable choice would yield anorthogonal dynamics system which can be solved in a more straightforward manner than the fullsystem. In some cases, the choice for the resolved variables and corresponding projection can bemotivated by physical arguments for the specific system. For other models the choice might not beas straightforward and one cannot be certain that a suitable reduced-order model exists.When the ODE system studied is linear, i.e. R ( φ ( x, t )) = Aφ ( x, t ) where A is a constant matrix,finding a suitable set of resolved variables can be done by looking at the eigenvalues of the orthogonaldynamics system. When the system (1) is linear, so is the orthogonal dynamics system (6), meaningits behaviour can be obtained by studying its eigenvalues, as (6) is linear in F i . A set of resolvedvariables is suitable if the eigenvalues of the orthogonal dynamics system show significantly morestability, i.e. have more negative real parts, compared to the full system. If this is not the casethe problem of solving the full system is transferred to the equally difficult problem of solving theorthogonal dynamics system.Up until now we have not discussed the difficulties arising when the MZ formalism is applied to aPDE system instead of ODEs. When the system is Hamiltonian some results exist (e.g. [5]), however,when this is not the case often the system is expanded in a basis of typically orthonormal functions(e.g. [8, 34]) to numerically find a solution. If the aim of applying the MZ formalism is to obtain aset of reduced-order model equations and not only a numerical result, this method is not suitable.Another approach, relying on integration along characteristics, has been explored by Falkena et al. (2019) and yielded an exact reduced- order (delay) model for the system studied [18]. Here we buildon this work to see whether delay type models can be derived for other systems of wave equations.In particular, we focus on a PDE model that describes thermal Rossby wave propagation related tothe AMO. This model is introduced in the next section. The thermal Rossby wave mechanism, suggested to be responsible for the AMO [9, 32], is summarizedin Figure 2. When there is a positive temperature anomaly ( T (cid:48) ) in the northern-central part of thebasin, the meridional temperature gradient increases with respect to the background state. Thisresults in a zonal overturning anomaly with westward surface flow through thermal wind balance(Figure 2a). The negative zonal flow transports the positive temperature anomaly towards thewestern boundary, creating a zonal temperature gradient. Again through thermal wind balance, thisnow leads to anomalies in the meridional overturning circulation (Figure 2b). This flow transportscold water from near the poles southward, reducing the meridional temperature gradient. Thissmaller north-south temperature gradient causes a positive (eastward) zonal flow, after which thesame pattern as described above is followed with a sign change. Hence the variability associated withthe AMO relies on the transport of heat and the flow response through the thermal wind balance[12].Low-order ad-hoc ODE models of the AMO have been studied by, among others, Broer et al.[2]. More recently, S´evellec and Huck [29] developed an idealized PDE model of the AMO model4 a) (b) Figure 2: Schematic diagram of the physical mechanism responsible for the AMO with two phasesa quarter period apart in panels (a) and (b). Figure taken from [12].Figure 3: A schematic diagram of the three-layer ocean basin considered in the AMO-model. Thedashed arrows show the background AMOC circulation, which is taken into account in the backgroundstate of the model discussed in Section 3.2.to which we apply the MZ formalism. This PDE model roughly captures the thermal Rossby wavemechanism described before, but with a simplification of the associated wave dynamics. In thissection, we present this model and briefly discuss its derivation. In addition we investigate the effectof using a more realistic strictly positive meridional overturning circulation as the background stateon the behaviour of the [29] model.
The AMO model by S´evellec and Huck (2015) is a three-layer model describing the evolution oftemperature perturbations in the North Atlantic Ocean [29]. The model describes the temperatures( T i , i = 1 , , ) as a function of latitude ( x ) and time ( t ). For convenience we consider the non-dimensional version of the model with latitude-scale W (basin width) and time-scale Y (a year).The scaled model is ∂ t T = a ∂ x T + b ∂ x T + c ∂ x T + κ s ∂ xx T ,∂ t T = a ∂ x T + b ∂ x T + c ∂ x T + κ s ∂ xx T ,∂ t T = κ s ∂ xx T , (7)with boundary conditions T i ( x | W est = 0) = − T i ( x | East = 1) , i = 1 , , . (8)5able 1: The values of the parameters in the AMO model by S´evellec and Huck [29].Thickness layer one h
600 mThickness layer two h
600 m Vertical temperature gradientThickness layer three h H ∂ z ¯ T = − Ch + h (cid:0) ∆ T − α S α T ∆ S (cid:1) Zonal basin size W L C Time scale (year) Y . · s Standard C = 1 Horizontal diffusivity κ · m /sAcceleration of gravity g Coriolis parameter f − s − β effect β . · − (ms) − Meridional temperature gradientThermal expansion coefficient α T · − K − Haline contraction coefficient α S · − psu − ∂ y ¯ T = L (cid:0) ∆ T − α S α T ∆ S (cid:1) Meridional temperature diff. ∆ T -20 KMeridional salinity diff. ∆ S -1.5 psuZonal velocity ¯ u − m/sThe constants in the model are all positive for physically realistic values and defined by a = YW (cid:16) α T g Hf (cid:16) − h ( h + h ) ∂ y ¯ T + β f h ( h + h ) ∂ z ¯ T (cid:17) − ¯ u (cid:17) ,b = YW α T g Hf (cid:16) − h ( h + 2 h ) ∂ y ¯ T + β f h h ( h + 2 h ) ∂ z ¯ T (cid:17) ,c = YW α T g Hf (cid:16) − h ∂ y ¯ T + β f h h ∂ z ¯ T (cid:17) ,a = YW α T g Hf (cid:16) h ∂ y ¯ T + β f h ( h + 2 h ) ∂ z ¯ T (cid:17) ,b = YW (cid:16) α T g Hf (cid:16) − h ( h − h ) ∂ y ¯ T + β f (4 h h h + h ( h + h )) ∂ z ¯ T (cid:17) − ¯ u (cid:17) ,c = YW α T g Hf (cid:16) − h ∂ y ¯ T + β f h (2 h + h ) ∂ z ¯ T (cid:17) .κ s = κ YW (9)The values of the parameters needed to compute these constants are given in Table 1.The derivation of these equations can be found in [29]. Here we briefly discuss that derivationand assumptions made to get to the above system of equations (7). The derivation starts from anadvection-diffusion equation for temperature, geostrophic balance (a balance between the horizontalpressure gradients and the Coriolis force), hydrostatic balance (a balance between the vertical pres-sure gradient and gravity) and a linear dependence of density on temperature. The equations fortemperature are linearized around a fixed background state comprised of zonal flow ¯ u and temper-ature gradients in the meridional ∂ y ¯ T and vertical ∂ z ¯ T . Note that this means that the overturningcirculation ( ¯ v , ¯ w ) is neglected because of its weakness with respect to the zonal flow. The linearizedtemperature equation is then discretized over three layers assuming no flow through the surface andbottom and no background flow or temperature gradients in the bottom layer, which results in system(7). The boundary conditions (8) are derived by considering a boundary layer at both sides of thebasin with free-slip and zero heat flux assumptions at the inner edge of the boundary layer. The fullderivation can be found in the appendix of [29].Here two additional simplifications to the model (7) are made. Firstly, we note that the only termacting in the third layer is diffusion and the two top layers do not couple into it. As a result anyperturbation in that layer eventually damps out. For this reason, and to simplify the mathematicaltreatment of the system, perturbations in the bottom layer are neglected (i.e. T = 0 ). Secondly, weapproximate the diffusion terms by linear damping with damping coefficient α . The system (7) then6able 2: The numerical values of the parameters. a a b b κ s . · − simplifies to a two-layer system: ∂ t T = a ∂ x T + b ∂ x T − αT ,∂ t T = a ∂ x T + b ∂ x T − αT . (10)This is the AMO model to which we apply the MZ formalism. Also note that this temperature modelexplains changes in the overturning circulation as well, via thermal wind balance, the continuityequation and Sverdrup balance, which is discussed in the Supplementary Information. The parametervalues used for the numerical results in the remainder of the section and coming sections are givenin Table 2.Before looking into the application of the MZ formalism to this AMO model, we illustrate itsbehaviour by simulating it for α = 0 . We use an upwind discretization scheme for the x -derivatives anda forward Euler scheme in time. Note that this discretization includes numerical diffusion leading toartificial damping effects. The result is shown in Figure 4. Note the opposite sign of the temperaturein the two layers, which is due to the baroclinic nature of the waves [7]. The model shows acombination of two oscillations with different periods. First, there is a long period of approximatelysixty years, which corresponds to a thermal Rossby wave responsible for driving the AMO. Secondly,there is a higher frequency oscillation with a period of around five years.The occurrence of the shorter period is at first sight surprising as it is not found in more detailedPDE models. This oscillation does not correspond to a planetary Rossby wave, as one might expect,since decreasing β does not result in a disappearance of these oscillations. It is a thermal Rossbywave, just as the one responsible for the AMO oscillation. The dominant appearance of this thermalRossby wave in the model is undesired to study the AMO. A possible improvement of the model,resulting in the damping of this high-frequency mode, is discussed in the next section. x T e m pe r a t u r e ( ° C ) T (0,x)T (0,x) (a) The initial conditions. Time (year) -1-0.500.51 T e m pe r a t u r e ( ° C ) T T (b) Evolution of model (10) at x = 0 . Figure 4: Model simulation of the temperature anomaly in the AMO model (Equation (10)) foran initial positive Gaussian temperature perturbation in the centre of the basin in the first layer( ∆ t = ∆ x = 0 . ). The AMO model by S´evellec and Huck [29] described in the previous section does not contain anoverturning circulation in the background state, as the background meridional ( ¯ v ) and vertical ( ¯ w )velocities are neglected. This means that in the model the overturning circulation, which can beinferred from the temperature evolution (details are given in the Supplementary Information), canbecome negative. To prevent this from happening in the model we consider an extended background7tate which retains meridional ¯ v and vertical ¯ w flow. With this different background state an extendedtwo-layer temperature model for AMO can be derived following the same steps as in [29]. Thedetails of this derivation can be found in the Supplementary Information. The resulting system fortemperature in the two upper layers is ∂ t T = a ∂ x T + b ∂ x T − ( β + α ) T − β T ,∂ t T = a ∂ x T + b ∂ x T − ( β + α ) T , (11)where β = Y · (cid:16) βf ¯ v + 2 h ¯ w (cid:17) , β = − Y h ¯ w, β = Y · (cid:16) βf ¯ v + 2 h ¯ w (cid:17) . (12)The difference to (10) is that there are additional linear terms in both equations. Note that not allthe additional terms have a dampening effect, as some of the β i -terms can be negative.A model simulation for ¯ v = 0 . · − m/s and ¯ w = − . · − m/s is shown in Figure5, where the values are chosen for plotting purposes being within a realistic range. Note thatif − w/h , (cid:29) β ¯ v/f , we have that β , become strongly negative leading to possible unstablesolutions or at least amplifying effects within the solution. The result of adding the backgroundoverturning circulation is a damping of the high frequency oscillation, as can be seen in Figure5. This short period oscillation is absorbed by the background overturning circulation while the
400 450 500 550 600 650 700
Time (year) -1-0.500.51 T e m pe r a t u r e ( ° C ) T T T T Figure 5: Simulations of the temperature in the two layers without (red, blue) and with (black, cyan)a background overturning circulation ( β = β = 1 . · − , β = 7 . · − ).long period oscillation persists. The amplitude of the oscillation corresponding to the AMO is notnoticeably affected by the damping. This can be due to the presence of an amplifying effect of thebackground overturning in some parts of the equations as mentioned previously. For simplicity weapply the MZ formalism to the AMO model as given in Equation (10) instead of the extended modeldiscussed here. However, one has to keep in mind that the high frequency oscillation present in themodel does not sustain in the presence of a background overturning state. The aim of applying the MZ formalism to the AMO model described in Section 3 is to arrive at aprojected model describing the same dynamics as the full model and analyze the effect of memory inthat system, with the potential of deriving a delay model for the phenomenon. A similar procedurehas been applied to a model of the El Ni˜no Southern Oscillation (ENSO) by Falkena et al. (2019)[18]. A difference is that for the AMO no previously proposed delay model is known. Therefore itis not immediately clear how to choose a projection, nor how to deal with solving the subsequentorthogonal dynamics equation (6). Preferably we arrive at an equation for the temperature at onelocation in space, to remove the explicit dependency on x in the system, but it is not clear from theonset whether or not this is feasible. 8ur first approach, as taken by Falkena et al. (2019) in the application of the MZ formalismto the ENSO model [18], is to use variation of constants to reduce the system to one dependentvariable. Since the model is linear it is possible to formally write down a solution for T (or T ) usingvariation of constants, which than can be substituted into the equation for T (or T ). This yields ∂ t T = a ∂ x T − αT + b ∂ x (cid:16) e ( b ∂ x − α ) t T (0) + (cid:90) t e ( b ∂ x − α )( t − s ) a ∂ x T ( s ) ds (cid:17) . (13)In Falkena et al. (2019) the next step taken was to use integration along characteristics to arriveat an equation that was localized in space. This method, which removes the exponentials in (13),is not suitable here because the exponential terms act on ∂ x T instead of T itself (which meansthat they refer to solutions of second order PDEs). The way in which we proceed is to first convertthe system of PDEs (10) into a set of ODEs by discretization. To this high-dimensional system ofODEs the MZ formalism is then applied. In the following sections this procedure is described. Afterdeciding on the discretization to use, possible (sets of) resolved variables are explored, followed by adiscussion of the different terms in the Langevin equation (4). The first step is to find a stable discretization of the AMO model in Equation (10). A grid of ( N +1) -points in space with distance dx = N is used. Because all parameters in the model are positive weknow that all waves travel westward. Therefore we use an upwind scheme to discretize the model.The discretized equations are ∂ t T n = a dx ( T n +11 − T n ) + b dx ( T n +12 − T n ) − αT n ,∂ t T n = a dx ( T n +11 − T n ) + b dx ( T n +12 − T n ) − αT n , (14)for n = 0 , ..., N (such that T ki ≈ T i ( k/N ) and dx = 1 /N ), with boundary conditions T N = − T , T N = − T . (15)By the circular nature of the boundary conditions this is a N -dimensional system (there is no need tosolve the dynamical equations for discretization points N ). Letting N → ∞ recovers the PDE modelexactly. This system (14) can be written as a matrix equation for (cid:126)T = ( T , T , ..., T N − , T N − ) : ∂ t (cid:126)T = M (cid:126)T . (16)The construction of M is straightforward from system (14).The stability of the solution of this discretized system of ODEs is verified by computing theeigenvalues of the matrix M . These are shown in Figure 6a for two different values of N . Foreach N two curves are visible with approximate equal spacing in the imaginary part between thesubsequent eigenvalues. For increasing N both curves of eigenvalues approach a line with real part − α , being the eigenvalue of the continuous system. Since all eigenvalues are negative for every N the discretization is stable. In the following sections we go into the application of the MZ formalismto this system of ODEs. When applying the MZ formalism the first step is to choose the resolved variables and correspondingprojection. As discussed in Section 2 this is an essential choice determining the final expressions inthe Langevin equation (4). From the modelling perspective the aim is to find a system of equationsfor the temperature at one location in order to remove the x -dependency of the system. With thisin mind there are three possible choices for the resolved variables; T , T , and both. Note thatbecause all waves travel in the same direction without loss of energy it does not matter on whichlocation in space the temperatures are projected. For convenience we choose to project onto thewestern boundary ( n = 0 ), but note that the result for any other location is the same. The moststraightforward way to project onto one of those sets of resolved variables is to use the linear projection9 - - - - - - - - ( λ k ) I m ( λ k ) (a) Full model for N = 200 (blue, circle) and N =400 (yellow, square). - - - - - ( λ k ) I m ( λ k ) (b) Full model for N = 400 (blue, circle) and usinga projection onto either T (yellow, square) or both T and T (green, diamond). Figure 6: The eigenvalues of the discretized AMO model ( α = 0 . ) (a) for different N and (b)different sets of resolved variables. Note that not all eigenvalues are shown (e.g. only one of the twoeigenvalues for the projection onto T and T , the second one is more negative). P , setting all unresolved variables to zero. Its complement Q thus sets the resolved variables to zerosince the system considered is linear. The corresponding orthogonal dynamics equation (6) can bewritten as a matrix equation ∂ t (cid:126)T Q = M Q (cid:126)T Q , (17)where (cid:126)T Q = ( T Q , T Q , ..., T N − Q , T N − Q ) represents the unresolved variables. The matrix M Q isthe same as the matrix for the full system M , but with the rows and columns corresponding to theresolved variable(s) removed.To assess the quality of the three possible projections we study the eigenvalues of the correspond-ing orthogonal dynamics systems. A sufficient decay of the eigenvalues of the orthogonal dynamicssystem indicates the corresponding resolved variables and projections are suitable as discussed inSection 2. In Figure 6b the eigenvalues of the full system and projected systems are shown for twodifferent projections on boundary variables. When one projects onto only T the eigenvalues arequite similar to those of the original system. Similar results are found for projection onto only T .If both T and T are taken as the resolved variables the result is noticeably better. The orthogonaldynamics system has only two eigenvalues which become increasingly negative with increasing N .Therefore this projection onto both T and T at the boundary x = 0 is chosen, where we noteonce again that the choice of the specific location is arbitrary. In the following section we focus onthe derivation of the noise and memory term in the MZ formalism (5). Here we briefly discuss theMarkovian terms. They are given by the projection of the right-hand side of the equation for T and T (Equation (14) for n = 0 ): P (cid:104) a dx ( T − T ) + b dx ( T − T ) − αT (cid:105) = − a dx T − b dx T − αT ,P (cid:104) a dx ( T − T ) + b dx ( T − T ) − αT (cid:105) = − a dx T − b dx T − αT . (18)This simply is the right-hand side dependence of Equation (14) on the resolved variables, as a linearprojection is used. 10 .3 Noise and Memory Terms To compute the noise and memory term we focus on the orthogonal dynamics system (17). For thechosen resolved variables T and T with the linear projection the matrix reads M Q = − a dx − α − b dx a dx b dx − a dx − b dx − α a dx b dx − a dx − α − b dx . . . − a dx − b dx − α . . .. . . a dx b dx . . . a dx b dx − a dx − α − b dx − a dx − b dx − α . (19)Note this matrix is block upper diagonal with all blocks on the diagonal being the same. To solvethe orthogonal dynamics system we have to find the eigenvalues and (generalized) eigenvectors ofthis matrix M Q . As discussed in Section 4.2 there are only two eigenvalues λ ± = − α − l ± dx , (20)with l ± = 12 (cid:16) a + b ± (cid:113) a + b − a b + 4 a b (cid:17) , (21)each of multiplicity N − . Note that l ± yield the characteristics of the original PDE system (10).The corresponding generalized eigenvectors for i = 1 , ..., N − are (cid:126)v i ± = (cid:16) dxl ± (cid:17) i − · (0 , ..., , w ± , , , ..., , (22)where the non-zero values are located on the coordinates corresponding to location i . Here w ± = 12 a (cid:16) a − b ± (cid:113) a + b − a b + 4 a b (cid:17) . (23)Having computed the eigenvalues and eigenvectors we can write down the solutions (cid:126)T Q of theorthogonal dynamics equation (e.g. [35]). Here we only note that (cid:126)T Q is a linear combination of theeigenvectors, meaning it is relatively straightforward to identify the solution at one location. The fullexpressions are given in the Supplementary Information, together with the use of initial conditions todetermine the constants involved.Now that we have the solution to the orthogonal dynamics equation we can write down the noiseterms and subsequently compute the memory terms of the discretized AMO system (14). The noiseterms are defined by F T ( t ) = a dx T Q ( t ) + b dx T Q ( t ) ,F T ( t ) = a dx T Q ( t ) + b dx T Q ( t ) . (24)Note that only the terms of the solution (cid:126)T Q which contain the eigenvectors (cid:126)v ± contribute to the noiseterm, as all other eigenvectors have zeros in the direction of T Q and T Q . The resulting expressions,following the solution of Equation (17), are F T ( t ) = N N − (cid:88) i =1 (cid:16) ( a w + + b ) e λ + t c i + + ( a w − + b ) e λ − t c i − (cid:17) t i − ( i − ,F T ( t ) = N N − (cid:88) i =1 (cid:16) ( a w + + b ) e λ + t c i + + ( a w − + b ) e λ − t c i − (cid:17) t i − ( i − , (25)11here c i + = (cid:16) l + dx (cid:17) i − · T i (0) − w − T i (0) w + − w − ,c i − = − (cid:16) l − dx (cid:17) i − · T i (0) − w + T i (0) w + − w − , (26)depend on the initial conditions of the unresolved variables.To compute the memory terms (as defined in (5)) we first look at the effect of applying theoperator P L to each of the initial conditions. This is sufficient for computation of the memory termsbecause the noise terms (25) are linear in the initial conditions. We find P L (cid:0) T (0) , T (0) , ..., T i (0) , T i (0) , ..., T N − (0) , T N − (0) (cid:1) = P (cid:0) ..., a dx ( T i +11 (0) − T i (0)) + b dx ( T i +12 (0) − T i (0)) − αT i (0) a dx ( T i +11 (0) − T i (0)) + b dx ( T i +12 (0) − T i (0)) − αT i (0) , ... (cid:1) , = (cid:0) , ..., , − a dx T (0) − b dx T (0) , − a dx T (0) − b dx T (0) (cid:1) . (27)We see that only terms that initially depend on T N − (0) and T N − (0) are non-zero after applicationof P L . Combining this result with the noise term (25) and replacing dx by N , the memory integrand(5) becomes K T (( T (0) , T (0) , t ) = N t N − ( N − e − αt (cid:16) ( l + N ) N − e − l + Nt (cid:0) A T (0) + B T (0) (cid:1) + ( l − N ) N − e − l − Nt (cid:0) A − T (0) + B − T (0) (cid:1)(cid:17) ,K T (( T (0) , T (0) , t ) = N t N − ( N − e − αt (cid:16) ( l + N ) N − e − l + Nt (cid:0) A T (0) + B T (0) (cid:1) + ( l − N ) N − e − l − Nt (cid:0) A − T (0) + B − T (0) (cid:1)(cid:17) , (28)with A = ( a w + + b )( − a + w − a ) w + − w − , A − = − ( a w − + b )( − a + w + a ) w + − w − ,B = ( a w + + b )( − b + w − b ) w + − w − , B − = − ( a w − + b )( − b + w + b ) w + − w − ,A = ( a w + + b )( − a + w − a ) w + − w − , A − = − ( a w − + b )( − a + w + a ) w + − w − ,B = ( a w + + b )( − b + w − b ) w + − w − , B − = − ( a w − + b )( − b + w + b ) w + − w − . (29)Now all components of the Langevin equation (4), being the Markovian terms (18), the noiseterms (25) and the memory integrands (28), are known. Thus we can write down the result of12pplying the MZ formalism to the discretized AMO system (14): ∂ t T = − a N T − b N T − αT + N e − αt N − (cid:88) i =1 (cid:16) ( a w + + b ) e − l + Nt c i + + ( a w − + b ) e − l − Nt c i − (cid:17) t i − ( i − (cid:90) t N ( t − s ) N − ( N − e − α ( t − s ) (cid:16) ( l + N ) N − e − l + N ( t − s ) (cid:0) A T ( s ) + B T ( s ) (cid:1) + ( l − N ) N − e − l − N ( t − s ) (cid:0) A − T ( s ) + B − T ( s ) (cid:1)(cid:17) d s,∂ t T = − a N T − b N T − αT + N e − αt N − (cid:88) i =1 (cid:16) ( a w + + b ) e − l + Nt c i + + ( a w − + b ) e − l − Nt c i − (cid:17) t i − ( i − (cid:90) t N ( t − s ) N − ( N − e − α ( t − s ) (cid:16) ( l + N ) N − e − l + N ( t − s ) (cid:0) A T ( s ) + B T ( s ) (cid:1) + ( l − N ) N − e − l − N ( t − s ) (cid:0) A − T ( s ) + B − T ( s ) (cid:1)(cid:17) d s. (30)This system depends on the discretization, or more precisely, on the number of points N . Ideallywe would like to find the equations for the continuous model. Therefore the limiting behaviour as N → ∞ , or /N = (cid:15) → , of the different terms is studied in the Supplementary Material. Theequation obtained after taking this limit can be written as (cid:15) d T d t = − a T ( t ) − b T ( t ) + A τ + e − ατ + T (cid:0) t − τ + (cid:1) + B τ + e − ατ + T (cid:0) t − τ + (cid:1) + A − τ − e − ατ − T (cid:0) t − τ − (cid:1) + B − τ − e − ατ − T (cid:0) t − τ − (cid:1) + (cid:15)f (cid:15) ( t ) + O ( (cid:15) ) ,(cid:15) d T d t = − a T ( t ) − b T ( t ) + A τ + e − ατ + T (cid:0) t − τ + (cid:1) + B τ + e − ατ + T (cid:0) t − τ + (cid:1) + A − τ − e − ατ − T (cid:0) t − τ − (cid:1) + B − τ − e − ατ − T (cid:0) t − τ − (cid:1) + (cid:15)f (cid:15) ( t ) + O ( (cid:15) ) , (31)where f (cid:15) ( t ) = − αT ( t ) + A τ + e − ατ + g (cid:15) + (cid:0) T (cid:1) + B τ + e − ατ + g (cid:15) + (cid:0) T (cid:1) + A − τ − e − ατ − g (cid:15) − (cid:0) T (cid:1) + B − τ − e − ατ − g (cid:15) − (cid:0) T (cid:1) ,f (cid:15) ( t ) = − αT ( t ) + A τ + e − ατ + g (cid:15) + (cid:0) T (cid:1) + B τ + e − ατ + g (cid:15) + (cid:0) T (cid:1) + A − τ − e − ατ − g (cid:15) − (cid:0) T (cid:1) + B − τ − e − ατ − g (cid:15) − (cid:0) T (cid:1) , (32)for g (cid:15) ± (cid:0) T (cid:1) = τ ± (cid:16)(cid:0) ( l ± + α ) − l ± (cid:1) T (cid:0) t − τ ± (cid:1) + 2( l ± + α ) T (cid:48) (cid:0) t − τ ± (cid:1) + T (cid:48)(cid:48) (cid:0) t − τ ± (cid:1)(cid:17) , (33)with τ ± = 1 /l ± . Note that we have dropped the superscript 0 in the notation, as the found systemis valid at every location throughout the basin through a simple coordinate transformation. This isthe final result of the application of the MZ formalism to the AMO model (10) as an expansion interms of order (cid:15) . Letting (cid:15) → a set of delay difference equations is found, giving the exact reducedmodel of the AMO. The Mori-Zwanzig formalism is in principle applicable also when the coefficients a j and b j are spacedependent. However, in that case it may not be possible to derive explicit expressions for delaysand coefficients in (31). For spatially constant coefficients a j and b j and equal damping α in allcomponents we may also derive the leading orders of (31) by integration along wave characteristics.This approach is similar to that taken in [18] and we refer to it as the method of characteristics(MoC). The damped free-wave solutions of the two-layer system, ∂ t T = a ∂ x T + b ∂ x T − αT , ∂ t T = a ∂ x T + b ∂ x T − αT , (34)13 + l − τ + τ − ∂ t ˜ T ± − l ± ∂ x ˜ T ± + α ˜ T ± = 0 , (35)along the wave characteristics, where the characteristic speeds and delays l ± = 12 (cid:104) a + b ± (cid:112) ( a + b ) − a b + 4 a b (cid:105) , τ ± = 1 /l ± , (36)are the same as for the Mori-Zwanzig formalism in (21). The new variables ˜ T ± are related back tothe original variables through the transformation P , (cid:126)T = P (cid:20) ˜ T + ˜ T − (cid:21) , where P = (cid:20) a − l + a a − l − a (cid:21) , P − = (cid:34) − a l + − l − a − l − l + − l − a l + − l − l + − a l + − l − (cid:35) . (37)The damped wave equations (35) have the general solutions ˜ T ± ( t, x ) = ˜ T ± ( t + xτ ± )e αxτ ± , (38)where the arbitrary profiles ˜ T ± are constrained by the boundary conditions ˜ T ± ( t,
0) = − ˜ T ± ( t, .These boundary conditions enforce the delay-difference equations for ˜ T ± , ˜ T ± ( t ) = − e − ατ ± ˜ T ± ( t − τ ± ) , (39)after shifting t by τ ± and multiplying both sides by e − ατ ± in the boundary conditions. Transforming ˜ T ( t ) back using the transformation P − gives the coupled delay equations (cid:126)T ( t ) = e − ατ + C (cid:126)T ( t − τ + ) + e − ατ − C (cid:126)T ( t − τ − ) , (40)(dropping the superscript of ˜ T for ease of notation) with C = (cid:34) l − − a l + − l − ( l + − a )( l − − a ) a ( l + − l − ) − a l + − l − − l + − a l + − l − (cid:35) , C = (cid:34) − l + − a l + − l − − ( l + − a )( l − − a ) a ( l + − l − ) a l + − l − l + − a l + − l − (cid:35) . The above is valid at every location in the basin according to the general solution form (38), asdiscussed in Section 4.3. The MZ derived system in Equation (31) for (cid:15) = 0 can be rewritten to theabove system.System (40) is a delay-difference system. In order to explore solutions of delay-difference system,we convert (40) to a system of DDEs y regularising it with a small time derivative (cid:15) d (cid:126)T / d t : (cid:15) d (cid:126)T ( t ) dt = − (cid:126)T ( t ) + e − ατ + C (cid:126)T ( t − τ + ) + e − ατ − C (cid:126)T ( t − τ − ) , (41)with (cid:15) (cid:28) . The choice of (cid:15) is related to the discretization of the original PDE system (34) through (cid:15) = 1 /N , where N is the number of discretization steps using an ‘upwind’ scheme (discretizing inthe direction of the wave). Table 3 shows the approximate resulting wave speeds and delays whenusing the parameter values in Table 2 for our numerical solutions and spectral analysis, which arediscussed in the next section. 14 Analysis of Delay Models
In this section we analyze the solutions of the delay difference models derived via the MZ formalism(31) and the MoC (41). We start with a discussion of the asymptotic spectrum and a spectral analysisof the MoC model (41). Next, we discuss the error terms as computed using the MZ formalism,followed by a comparison of the MZ delay model, the MoC delay model and the numerical PDEsolution from which the MZ model is derived.
We observe that the delays occurring in the MoC system are of different magnitude: τ + (cid:28) τ − , where τ + is of order O (1) in the time scale of (41). For hierarchical large delays [25] and [28] provide asimple approximation of the spectrum for (41) which captures the range of possible curvatures of thecurves along which the eigenvalues shown in Fig. 6 align. Any eigenvalue λ of (41) satisfies det (cid:104) − (cid:15)λI − I + C e − ( α + λ ) τ + + C e − ( α + λ ) τ − (cid:105) = 0 (42)with (cid:15) (cid:28) , τ + = O (1) and τ + (cid:28) τ − . Hence the term e − ( α + λ ) τ − is negligibly small unless ατ − and Re λτ − are of order or less. If we assume that this is the case, we may introduce α − = α/τ − (which is of order or less) and look for eigenvalues λ of the form λ = γ/τ + + i ω/(cid:15) (called thepseudo-continuous spectrum in [25]). Then ( γ + α − ) τ + /τ − and (cid:15)γ/τ − are small. Dropping theseterms and introducing the phases φ ± = ωτ ± /(cid:15) and z = e − ( γ + α − )+i φ − simplifies (42) to det (cid:2) − i ωI − I + C e i φ + + C z (cid:3) = 0 (43)Equation (43) is in our case a quadratic equation in the complex number z , giving two roots, eachdepending on ω and φ + , which one may express as z ± ( ω, φ + ) . From this root pair one may derivethe damping depending on the frequency, γ ± ( ω, φ + ) = − α − + log z ± ( ω, φ + ) , and, in the originalscaling for eigenvalue λ Re λ = − α + τ + log z ± ( (cid:15) Im λ, φ + ) .This relation determines the curves along which the eigenvalues align for positive small (cid:15) and τ + (cid:28) τ − .The phase φ + is treated here as an independent parameter. It is very sensitive with respect to smallchanges of τ + (since φ + = ωτ + /(cid:15) ) such that the location of the eigenvalue curves will vary stronglydepending on τ + or (cid:15) within the range given by φ + ∈ [0 , π ] . Ruschel and Yanchuk’s [28] analysisshows in general that for hierarchically large delays the spectrum “fills an area” of the complex planeunder small parameter variations. For the spectral analysis of the MoC delay model (41) (with α = 0 ) we compute the trajectories of T and T . To compute the history needed for the difference equation, the PDE system (34) is solvednumerically for an initial profile of the basin using an upwind discretization scheme for τ − years. Wetake a Gaussian initial distribution profile (same as Figure 4a). The DDE system (41) is then evolvedfor a further 200 time steps. Figure 7 shows the results.A spectral analysis is performed on the resulting trajectories to identify the most prominentoscillation periods. A dominant signal of a τ − year cycle is obtained, along with a smaller signal fora τ + year cycle. These two most prominent signals correspond to period doubling of the two delayvalues which arises naturally from the boundary conditions. There is a smaller peak correspondingto a cycle of approximately τ − years. The signals corresponding to τ − (53.07) and τ − (17.77)year cycles align with the literature regarding possible cycle lengths of the AMO [10, 6]. The signalcorresponding to the τ + year cycle is much less pronounced in the surface temperature comparedto the subsurface temperature. We start this section with an evaluation of the theoretical error terms f (cid:15),i for i = 1 , (Equation (32))as computed via the MZ formalism. An example of the evolution of these terms has been plotted in15
50 100 150 200
Time (year) -1-0.500.51 T e m pe r a t u r e ( ° C ) (a) frequency (years -1 ) s pe c t r a l den s i t y T T (b) Figure 7: Numerical results for (41) with parameters from Table 3. (a) Trajectory for 200 years. (b)Power spectral density.Figure 8a. Also shown is the decrease of their maximum amplitude with increasing N (decreasing (cid:15) ),which corresponds to the effect expected from increasing the number of steps in a discretized PDE.Thus the theoretical error of the delay system derived using the MZ formalism (31) indicates thatthe delay model exhibits an error similar to that of the discretized PDE. In Figure 8b the spectraldensity of f (cid:15)i for i = 1 , is shown. We find a peak at a frequency of 0.53 years − which emergesdue to second derivatives in the computation of f (cid:15)i . The small low-frequency peaks correspond tothose of the exact delay model, as shown in Figure 7b, as there is a delayed contribution from thetemperature itself(and its first derivative) to the error term as well.
50 60 70 80 90 100
Time (year) -202 e rr o r ( ° C ) -3 N=2000 N -3 -2 f (t) f (t) (a) Frequency (years -1 ) S pe c t r a l D en s i t y N=2000 f f (b) Figure 8: The error terms f (cid:15) , f (cid:15) as computed using the MZ formalism. (a) The change ofamplitude of the terms with N is shown together with an example of the error terms for N = 2000 .(b) The spectrum of f (cid:15) and f (cid:15) .Next we compare the performance of the two DDE models derived via the MZ projection (Sections4.1, 4.2, 4.3) and the MoC (Section 4.4) with the exact and discretized PDE models. We do thisthrough the calculation of the eigenfunctions of the respective models. Figure 9a shows the realpart of the eigenfunctions for each component of the exact PDE calculated using the Chebfun open-source software [15]. We also compute the eigenfunctions for each approximation to the exact PDE:discretized PDE, delay via method of characteristics (MoC), and delay via Mori-Zwanzig (MZ). Wescale all the eigenfunctions such that V ( x = 0) = 1 . In order to compare the relative approximations,we calculate the error with respect to the eigenfunctions satisfying the PDE boundary conditions (8):16 x -1-0.500.51 r ea l ( V ) V (x)V (x) (a) -3 -2 -4 -3 -2 -1 e rr o r i MZMoCdisc PDE (b)
Figure 9: (a) The real part of the eigenfunctions of the PDE and (b) the errors of the discretizedPDE (plotted for (cid:15) = N ), MoC and MZ model as calculated using the eigenfunctions. error i := | V i ( x = 0) + V i ( x = 1) | for i ∈ { disc PDE , MoC , MZ } (44)As the approximate systems approach the exact PDE, the eigenfunctions of the respective modelsare expected to converge and therefore satisfy the PDE boundary conditions. We then would expect(44) to approach zero as N ( (cid:15) ) is increased (decreased) if the models are good approximations ofthe exact PDE. We calculate (44) for a range of N and (cid:15) values and plot the respective errors inFigure 9b. It can be seen that decreasing the (cid:15) term in the delay equations has the same effecton the approximation to the exact PDE as increasing N in the discretization. In other words, theerror introduced through the (cid:15) term in the ‘smoothing’ approximation of delay difference equationsis proportional to the error introduced by discretization methods of wave equation PDEs ( (cid:15) ∝ /N ). A delay model for the AMO has been derived from a three-layer model by S´evellec and Huck [29] usingthe MZ formalism. This formalism gives a rewriting of a system of ordinary differential equations [4]which contains a Markovian, a noise and a memory term. The advantage of this delay model is thatit precisely shows the propagating wave nature of the AMO, through the (inverse) travel times l ± ,hence providing more support for the thermal Rossby wave mechanism proposed in [32]. In a similarway as the ENSO delayed oscillator model [31, 1], this model can also be used to study effects ofbackground state, non-stationary forcing, noise, and possibly state-dependent delay versions, on thebehavior of the AMO [33, 24].The derived model for the AMO is a first order delay difference model, in contrast to the delaydifferential model for ENSO [18]. This means that the current state is fully determined by paststates. This type of model can exhibit an increasing switching frequency between states [20], makingit physically unrealistic. Hence, an (cid:15) dd t -term was added to prevent this behaviour and allow for betternumerical treatment. We were able to derive an error term for this approximation using the MZformalism, and relate this error to the upwind discretization scheme used in solving the original PDEmodel. For the non-damped version of the AMO model ( α = 0 ) the MZ formalism is not strictlynecessary for deriving a delay model. The method of characteristics yields the same delay differenceequation, as shown in Section 4.4. We showed that the way in which the smoothing approximation (cid:15) dd t -term is added affects the size of the error of the delay-model. Furthermore, the error introducedthrough the smoothing (cid:15) dd t -term is proportional to that introduced by discretization methods, asdiscussed in Section 5.3.The PDE model of the AMO by S´evellec and Huck [29], the starting point for deriving the delaymodel, does not contain a background overturning circulation. This results in a high frequencymodel oscillation which is undesired to study the AMO. As discussed in Section 3.2, adding themeridional overturning circulation to the background state of the model results in a damping of17his high frequency oscillation. It would be interesting to see what the effect of this backgroundoverturning is on the resulting delay equations. The addition of the overturning results in a changeof the eigenvalues of the discretized system. These eigenvalues then show a stronger decay for thehigh frequency mode and a weaker decay for the low frequency mode. It is expected that this willbe expressed in the coefficients of the delay terms in the resulting delay model. We do not expectan effect on the delay times, but this remains to be verified. The application of the Mori-Zwanzigformalism to this extended AMO model is beyond the scope of the present paper.The method of deriving delay equations as applied to a PDE model of the AMO (and alsoto a PDE model of ENSO [18]) can be generalized. It is expected that any diagonalizable linearsystem of wave equations can be rewritten in the form of a delay difference system. Integrationalong characteristics would yield the dominant terms, but the MZ formalism additionally gives errorterms to a smoothening approximation for solving the delay difference system. The necessity ofthe diagonalizability remains to be investigated. For non-diagonalizable systems it is not yet clearwhether a similar result would hold as the computations become more involved.For nonlinear models the method applied here to the AMO model has to be generalized. Thedifficulty lies in solving the orthogonal dynamics equation as noted extensively in the literature [4, 8,34]. For linear models the pseudo-orthogonal dynamics approximation can be used to arrive at thefinal result [21]. For nonlinear models this approximation cannot be used without first showing itsaccuracy. In [18] it was shown to result in a significant error for the nonlinear ENSO model studied.A better approximation for nonlinear models first needs to be proposed before the method describedhere can be accurately generalized. This step is necessary for a reliable and accurate application ofthe MZ formalism to nonlinear models of climate phenomena.Many PDE models used to describe climate phenomena contain some type of wave dynamics.We have shown in this study that projecting a system of wave equations onto one location yieldsa delay model. This would imply that more climate variability phenomena could be described by adelay equation. Once better methods of approximation of the orthogonal dynamics are available,it may be possible to derive accurate nonlinear delay models of climate phenomena, thus clarifyingdynamical mechanisms and allowing for further analysis. Author Contributions
S.K.J.F., C.Q. and H.A.D. designed the study. The work was carried out mainly by S.K.J.F. and C.Q..J.S. designed the analysis and comparison of the delay models in Section 5. All authors contributedto the work, discussed the results and read and approved the manuscript.
Funding
This work was supported by funding from the European Union Horizon 2020 research and innovationprogramme for the ITN CRITICS under Grant Agreement no. 643073 (C.Q., J.S. and H.A.D.).S.K.J.F. was supported by the Centre for Doctoral Training in Mathematics of Planet Earth, UKEPSRC funded (grant EP/L016613/1) and J.S. by EPSRC grants no. EP/N023544/1 and no.EP/N014391/1.
References [1] D S Battisti and A C Hirst. Interannual Variability in a Tropical Atmosphere-Ocean Model:Influence of the Basic State, Ocean Geometry and Nonlinearity.
Journal of the AtmosphericSciences , 46(12):1687–1712, 1989.[2] H W Broer, H A Dijkstra, C Sim´o, A E Sterk, and R Vitolo. The Dynamics of a Low-OrderModel for the Atlantic Multidecadal Oscillation.
Discrete and Continuous Dynamical Systems ,16:73–107, 2011.[3] Anson H Cheung, Michael E Mann, Byron A Steinman, Leela M Frankcombe, Matthew HEngland, and Sonya K Miller. Comparison of low-frequency internal climate variability in cmip5models and observations.
Journal of Climate , 30(12):4763–4776, 2017.184] A J Chorin, O H Hald, and R Kupferman. Optimal prediction with memory.
Physica D , (166),2002.[5] Alexandre J. Chorin, Raz Kupferman, and Doron Levy. Optimal Prediction for HamiltonianPartial Differential Equations.
Journal of Computational Physics , 162(1):267–297, 2000.[6] Petr Chylek, Chris K. Folland, Henk A. Dijkstra, Glen Lesins, and Manvendra K. Dubey. Ice-coredata evidence for a prominent near 20 year time-scale of the Atlantic Multidecadal Oscillation.
Geophysical Research Letters , 38:L13704, 2011.[7] B Cushman-Roisin and J Beckers.
Introduction to Geophysical Fluid Dynamics , volume 101.Academic Press, 2011.[8] E Darve, J Solomon, and A Kia. Computing generalized Langevin equations and generalizedFokker-Planck equations.
Proceedings of the National Academy of Sciences , 106(27):10884–10889, 2009.[9] C de Verdi`ere and T Huck. Baroclinic Instability: An Oceanic Wavemaker for InterdecadalVariability.
Journal of Physical Oceanography , 29:893–910.[10] T. L. Delworth and M. E. Mann. Observed and simulated multidecadal variability in the NorthernHemisphere.
Climate Dynamics , 16(9):661–676, 2000.[11] Clara Deser, Michael a. Alexander, Shang-Ping Xie, and Adam S. Phillips. Sea Surface Temper-ature Variability: Patterns and Mechanisms.
Annual Review of Marine Science , 2(1):115–143,2010.[12] H A Dijkstra. Interaction of SST Modes in the North Atlantic Ocean.
Journal of PhysicalOceanography , 36:286–299, 2006.[13] H A Dijkstra.
Dynamical Oceanography . Springer, 2008.[14] Henk A Dijkstra.
Nonlinear Physical Oceanography , volume 28. Springer, 2nd revise edition,2005.[15] T. A Driscoll, N. Hale, and L. N. Trefethen.
Chebfun Guide . Pafnuty Publications, 2014.[16] D.B. Enfield, A.M. Mestas-Nunez, and P.J. Trimble. The Atlantic Multidecadal Oscillation andits relationship to rainfall and river flows in the continental U.S.
Geophysical Research Letters ,28:2077–2080, 2001.[17] A Erd´elyi.
Asymptotic Expansions . Dover Publications, 1956.[18] Swinda K. J. Falkena, Courtney Quinn, Jan Sieber, Jason Frank, and Henk A. Dijkstra. Deriva-tion of delay equation climate models using the Mori-Zwanzig formalism.
Proceedings of theRoyal Society A: Mathematical, Physical and Engineering Sciences , 475(2227):20190075, 2019.[19] Leela M. Frankcombe, Anna von der Heydt, and Henk A. Dijkstra. North Atlantic MultidecadalClimate Variability: An Investigation of Dominant Time Scales and Processes.
Journal OfClimate , 23(13):3626–3638, 2010.[20] M Ghil, I Zaliapin, and B Coluzzi. Boolean delay equations: A simple way of looking at complexsystems.
Physica D , 237:2967–2986, 2008.[21] A Gouasmi, E J Parish, and K Duraisamy. A priori estimation of memory effects in reduced-order models of nonlinear systems using the Mori-Zwanzig formalism.
Proceedings of the RoyalSociety A , (20170385), 2017.[22] Zhe Han, Feifei Luo, Shuanglin Li, Yongqi Gao, Tore Furevik, and Lea Svendsen. Simulationby cmip5 models of the atlantic multidecadal oscillation and its climate impacts.
Advances inAtmospheric Sciences , 33(12):1329–1342, 2016.[23] A Keane, B Krauskopf, and C M Postlethwaite. Climate models with delay differential equations.
Chaos , 27(114309), 2017. 1924] Andrew Keane, Bernd Krauskopf, and Henk A. Dijkstra. The effect of state dependence in adelay differential equation model for the El Ni˜no Southern Oscillation.
Philosophical Transactionsof the Royal Society A: Mathematical, Physical and Engineering Sciences , 377(2153), 2019.[25] Mark Lichtner, Matthias Wolfrum, and Serhiy Yanchuk. The spectrum of delay differentialequations with large delay.
SIAM journal on mathematical analysis , 43(2):788–802, 2011.[26] H Mori. Transport, Collective Motion and Brownian Motion.
Progress of Theoretical Physics ,33(3):423–455, 1965.[27] Gary P Morriss and Denis J Evans.
Statistical Mechanics of Nonequilbrium Liquids . ANU Press,2013.[28] Stefan Ruschel and Serhiy Yanchuk. The spectrum of delay differential equations with multiplehierarchical large delays. arXiv preprint arXiv:1902.00404 , 2019.[29] F S´evellec and T Huck. Theoretical Investigation of the Atlantic Multidecadal Oscillation.
Journal of Physical Oceanography , 45, 2015.[30] M. Srokosz, M. Baringer, H. Bryden, S. Cunningham, T. Delworth, S. Lozier, J. Marotzke, andR. Sutton. Past, present, and future changes in the atlantic meridional overturning circulation.
Bulletin of the American Meteorological Society , 93(11):1663–1676, 2012.[31] M J Suarez and P S Schopf. A Delayed Action Oscillator for ENSO.
Journal of the AtmosphericSciences , 45(21):3283–3287, 1988.[32] L A Te Raa and H A Dijkstra. Instability of the Thermohaline Ocean Circulation on InterdecadalTimescales.
Journal of Physical Oceanography , 32:138–160, 2002.[33] Eli Tziperman, Lewi Stone, Mark A. Cane, and Hans Jarosh. El Ni˜no chaos: Overlapping ofresonances between the seasonal cycle and the pacific ocean-atmosphere oscillator.
Science ,264(5155):72–74, 1994.[34] Yuanran Zhu and Daniele Venturi. Faber approximation of the Mori–Zwanzig equation.
Journalof Computational Physics , 372:694–718, 2018.[35] Dennis G. Zill and Michael R. Cullen.
Differential Equations with Boundary Value Problems .Brooks/Cole Cengage Learning, 7 edition, 1997.[36] R Zwanzig. Nonlinear Generalized Langevin Equations.
Journal of Statistical Physics , 9(3):215–220, 1973. 20 upplementary Material
A AMO Model
A.1 Relation to Overturning Circulation
The AMO temperature model (10) also provides information about oscillations in the zonal andmeridional overturning circulation. Using thermal wind balance f ∂ z v = α T g∂ x T, (45)the vertical shear in the meridional flow can be related to the zonal temperature gradient. Thisindicated there is a quarter phase difference between the two (as seen in [32]). The vertical shear inthe meridional flow indicates whether there is a positive or negative perturbation in the meridionaloverturning. A positive perturbation in ∂ z v means more northward flow at the surface comparedto the bottom, and thus a positive perturbation in the overturning circulation. Similarly a negativeperturbation in ∂ z v corresponds to a negative overturning perturbation. Note that in the modelconsidered these flow perturbations do not result in temperature perturbations by the assumption ofdominant zonal dynamics.To get an idea of the behaviour of perturbations in the zonal overturning the y -averaged continuityequation is considered: ∂ x u + ∂ z w = 0 , (46)where it is assumed that there is no flow through the boundaries of the basin. Taking the z -derivativeof this equation and using Sverdrup balance βv = f ∂ z w, (47)yields the following equation ∂ x (cid:0) ∂ z u (cid:1) = − βf ∂ z v = − βf α T gf ∂ x T. (48)This indicates there is a difference of half a phase between zonal overturning perturbations andtemperature perturbations. It also implies that the phase difference between the zonal and meridionaloverturning perturbations is a quarter phase, as is expected from literature [32].In Figure S1 the evolution of the perturbations in the vertical shear in the zonal and meridionaldirections is shown as computed from the temperature simulation. For both the short and longperiod oscillations, the quarter and half phase difference between the temperature oscillations andthe meridional and zonal overturning perturbations respectively can be seen. A positive peak in tem-perature coincides with a negative zonal overturning perturbation, which transports this perturbationwestward. Physically a short delay between these two would be expected, but due to the assumptionof immediate thermal wind balance this is not present in the model. The resulting zonal temperaturegradient causes a negative meridional overturning perturbation with equatorward surface flow laggingby a quarter phase. This is followed by a negative temperature perturbation and a positive zonaloverturning, inducing a positive meridional overturning perturbation after which the oscillation startsagain. This represents the physical process of the AMO as described in Section 1. A.2 Derivation of Extended Background AMO Model
The temperature equation, including background overturning is: ∂ t T = − ¯ u∂ x T − ¯ v∂ y T − ¯ w∂ z T − v∂ y ¯ T − w∂ z ¯ T + κ∂ xx T, (49)where κ is the horizontal eddy diffusivity coefficient.21
20 140 160 180 200
Time (year) -1-0.500.51 z v (0.1 s -1 ) z u (0.1 s -1 ) T ( ° C) Figure S1: The evolution of the vertical shear in the meridional and zonal direction (blue) togetherwith the temperature oscillations (red) in the top layer.We start by looking at the y -derivative of temperature. ∂ y T . The thermal wind balance equationsfor the meridional and zonal vertical shear are f ∂ z v = α T g∂ x T, f ∂ z u = − α T g∂ y T. (50)The second equation means ∂ y T can be replaced by ∂ z u in the temperature equation. This gives ∂ t T = − ¯ u∂ x T + ¯ v fα T g ∂ z u − ¯ w∂ z T − v∂ y ¯ T − w∂ z ¯ T + κ∂ xx T. (51)Now we move or attention to the z -derivative of zonal wind ∂ z u . We use the Sverdrup balance,given in Equation (47) and plug this into the continuity equation. This yields ∂ x u + ∂ y v + βf v = 0 . (52)Integrating this equation over y and assuming that there is no flow through the northern and southernboundary, so v | S = v | N = 0 , gives ∂ x u = − βf v, (53)or ∂ x ∂ z u = − βf ∂ z v = − βf α T gf ∂ x T. (54)This implies that ∂ z u + βf α T gf T constant is in x . Assuming it is zero initially, the result for thetemperature equation is ∂ t T = − ¯ u∂ x T − ¯ v βf T − ¯ w∂ z T − v∂ y ¯ T − w∂ z ¯ T + κ∂ xx T. (55)The next step is to discretize this system over three layers, as in [29]. For this discretization wefollow the same steps as [29], so we have for vv − v = 12 ( h ∂ z v + h ∂ z v ) ,v − v = 12 ( h ∂ z v + h ∂ z v ) . (56)22ne more condition is needed to be able to express v in terms of ∂ z u (and thus in terms of T ) and ∂ z T in terms of T . For v we use the baroclinic condition, which in the discretized case reads h v + h v + h v = 0 . (57)Solving the system of the above three equations for all v i yields v = 12 H ( h ( h + h ) ∂ z v + h ( h + 2 h ) ∂ z v + h φ ) ,v = 12 H ( − h ∂ z v + h ( h − h ) ∂ z v + h ∂ z v ) ,v = 12 H ( − h ∂ z v − h (2 h + h ) ∂ z v − h h ∂ z v ) . (58)Here thermal wind balance can be used to express v i in terms of ∂ x T i .To express the vertical velocity w in terms of temperature we use the Sverdrup balance (47) toexpress it in terms of v . With the assumption of zero flow at the bottom, so w | − H = 0 , this yieldsfor the first two layers w = − β f h v ,w = − β f (2 h v + h v ) . (59)We assume the third layer is at rest, meaning the value in the third layer is not needed.For the discretization of temperature over the three layers the equations are: T − T = 12 ( h ∂ z T + h ∂ z T ) ,T − T = 12 ( h ∂ z T + h ∂ z T ) , (60)where it is assumed that there are no temperature perturbations at the bottom of the basin. We findthe vertical gradients of temperature to be ∂ z T = 2 h ( T − T + 2 T ) ,∂ z T = 2 h ( T − T ) ,∂ z T = 2 h T . (61)Combining the results above and assuming the third layer is at rest the system we find is ∂ t T = a ∂ x T − (¯ v βf + ¯ w h ) T + b ∂ x T + ¯ w h T + c ∂ x T − ¯ w h T + κ∂ xx T ,∂ t T = a ∂ x T + b ∂ x T − (¯ v βf + ¯ w h ) T + c ∂ x T + ¯ w h T + κ∂ xx T ,∂ t T = κ∂ xx T , (62)with boundary conditions T i | W est = − T i | East , i = 1 , , . (63)23he constants here are all positive for physically realistic values and defined by a = α T g Hf (cid:16) − h ( h + h ) ∂ y ¯ T + β f h ( h + h ) ∂ z ¯ T (cid:17) − ¯ u,b = α T g Hf (cid:16) − h ( h + 2 h ) ∂ y ¯ T + β f h h ( h + 2 h ) ∂ z ¯ T (cid:17) ,c = α T g Hf (cid:16) − h ∂ y ¯ T + β f h h ∂ z ¯ T (cid:17) ,a = α T g Hf (cid:16) h ∂ y ¯ T + β f h ( h + 2 h ) ∂ z ¯ T (cid:17) ,b = α T g Hf (cid:16) − h ( h − h ) ∂ y ¯ T + β f (4 h h h + h ( h + h )) ∂ z ¯ T (cid:17) − ¯ u,c = α T g Hf (cid:16) − h ∂ y ¯ T + β f h (2 h + h ) ∂ z ¯ T (cid:17) . (64)If, as done for the original system, we neglect the third layer the system becomes ∂ t T = a ∂ x T − (¯ v βf + ¯ w h ) T + b ∂ x T + ¯ w h T + κ∂ xx T ,∂ t T = a ∂ x T + b ∂ x T − (¯ v βf + ¯ w h ) T + κ∂ xx T . (65) B Application of the Mori-Zwanzig Formalism
B.1 Solution to the Orthogonal Dynamics System
The solution to the orthogonal dynamics system (17) for the AMO model can be written using theeigenvalues and generalized eigenvectors of M Q (19) as given in Equations (20) and (22). It is (cid:126)T Q ( t ) = e λ + t (cid:16) c (cid:126)v + c ( t(cid:126)v + (cid:126)v ) + ... + c i + (cid:16) t i − ( i − (cid:126)v + t i − ( i − (cid:126)v + ... + (cid:126)v i + (cid:17) + ... + c N − (cid:16) t N − ( N − (cid:126)v + ... + (cid:126)v N − (cid:17)(cid:17) + e λ − t (cid:16) c − (cid:126)v n + c − ( t(cid:126)v − + (cid:126)v − ) + ... + c i − (cid:16) t i − ( i − (cid:126)v − + t i − ( i − (cid:126)v − + ... + (cid:126)v i − (cid:17) + ... + c N − − (cid:16) t N − ( N − (cid:126)v − + ... + (cid:126)v N − − (cid:17)(cid:17) . (66)Here the constants c i ± are determined by the initial conditions. Each generalized eigenvector (22)has only components in the directions of T i and T i . This means that, to find expressions for theconstants, the following system has to be solved for each i : c i + w + (cid:16) dxl + (cid:17) i − + c i − w − (cid:16) dxl − (cid:17) i − = T i (0) ,c i + (cid:16) dxl + (cid:17) i − + c i − (cid:16) dxl − (cid:17) i − = T i (0) . (67)The solution is c i + = (cid:16) l + dx (cid:17) i − · T i (0) − w − T i (0) w + − w − ,c i − = − (cid:16) l − dx (cid:17) i − · T i (0) − w + T i (0) w + − w − . (68)This way an analytical solution to the orthogonal dynamics equation for general initial conditions hasbeen found. 24 .2 Limiting Behaviour In studying the limiting behaviour of Equation (30) we focus on the memory integral first, as this isthe term that is expected to result in a delay [18]. The dependence on N of all terms in the memorykernel can be described by the function f N ( t ) = N t N − ( N − e − αt ( µN ) N − e − µNt , (69)for µ = l ± . In Figure S2 this function is plotted for several N and fixed µ showing a peak at /µ which increases with height as N gets larger. This behaviour is caused by the projection ontoone location, which prevents waves from travelling through the basin in the orthogonal dynamicssystem. The result is an accumulation of energy at the location of the resolved variables. Becauseone integrates over the memory kernel to obtain the memory term the blow up in peak height doesnot necessarily imply a blow up of the memory term itself. ( year ) f N Figure S2: The function f N ( t ) with µ = 0 . for N = 25 (blue, solid), N = 50 (yellow, dashed) and N = 100 (green, dotted).To find an approximate expression for the memory integral we use Laplace’s approximation [17],which states that for a large number N , a smooth function h ( x ) , and a twice differentiable function g ( x ) , the following integral can be approximated by (cid:90) ba h ( x ) e Ng ( x ) d x = (cid:115) πN | g (cid:48)(cid:48) ( x ) | e Ng ( x ) · (cid:16) h ( x ) + 1 N (cid:16) − h (cid:48)(cid:48) ( x )2 g (cid:48)(cid:48) ( x ) + h ( x ) g (cid:48)(cid:48)(cid:48)(cid:48) ( x )8( g (cid:48)(cid:48) ( x )) + h (cid:48) ( x ) g (cid:48)(cid:48)(cid:48) ( x )2( g (cid:48)(cid:48) ( x )) − h ( x )( g (cid:48)(cid:48)(cid:48) ( x )) g (cid:48)(cid:48) ( x )) (cid:17) + O (cid:16) N (cid:17)(cid:17) , (70)provided that there is an x ∈ ( a, b ) such that g ( x ) is only close to g ( x ) if x is close to x . At thatpoint g ( x ) is required to have a maximum, i.e. g (cid:48)(cid:48) ( x ) < .All components in the memory integrals (30) can be written in the form h ij ± ( s ) e Ng ± ( s ) , for i = 1 , , j = 1 , , with g ± ( s ) = log( l ± ( t − s )) − l ± ( t − s ) ,h ij ± ( s ) = N N ( N − l ± ( t − s )) − e − α ( t − s ) C j ± T i ( s ) , (71)where C j ± is either A j ± or B j ± as given in Equation (29). The function h ij ± is smooth meaningthe first of the conditions for Laplace’s approximation is satisfied. To verify the other two conditionsthe first two derivatives of g ± ( s ) are required. We find the first derivative exists and is only zero for s = t − l ± indicating there is only one extreme. This means that g ± ( s ) is only close to g ± ( s ) for s s , satisfying the second condition. The second derivative is g (cid:48)(cid:48)± ( s ) = − l ± < , indicating thisextreme is a maximum meaning the third condition is satisfied. Therefore we are justified in applyingLaplace’s approximation, which yields: (cid:90) t h ij ± ( s ) e Ng ± ( s ) d s = (cid:115) πN l ± e − N N N ( N − e − αl ± C j ± (cid:16) T i (cid:0) t − l ± (cid:1) + 1 N l ± (cid:16) (2 l ± ( l ± + α ) + α ) T i (cid:0) t − l ± (cid:1) + 2( l ± + α ) T (cid:48) i (cid:0) t − l ± (cid:1) + T (cid:48)(cid:48) i (cid:0) t − l ± (cid:1)(cid:17) + O (cid:16) N (cid:17)(cid:17) . (72)In this equation delay terms have emerged due to the location of the maximum of g ± ( s ) at s = t − l ± .This means that the AMO can be modelled by some type of delay equation, similar to the resultsfound for ENSO [18].To simplify this expression we expand the first part of Equation (72) using a Taylor expansion.Expanding in /N around zero we find (cid:115) πN l ± e − N N N ( N − Nl ± (cid:16) − N + O (cid:16) N (cid:17)(cid:17) , (73)allowing for the expansion of the memory integrals in terms of order (cid:15) = N . Equation (72) thusbecomes (cid:90) t h ij ± ( s ) e Ng ± ( s ) d s = Nl ± e − αl ± C j ± (cid:16) T i (cid:0) t − l ± (cid:1) + (cid:15) l ± (cid:16)(cid:0) ( l ± + α ) − l ± (cid:1) T i (cid:0) t − l ± (cid:1) + 2( l ± + α ) T (cid:48) i (cid:0) t − l ± (cid:1) + T (cid:48)(cid:48) i (cid:0) t − l ± (cid:1)(cid:17) + O ( (cid:15) ) (cid:17) . (74)Here we keep one factor of N = (cid:15) to align with the Markovian term in Equation (30). We now haveconcise expressions for the memory terms in Equation (30).The behaviour of the noise terms for increasing N also needs to be investigated. Each componentof the sum in the noise terms is proportional to N e − l ± Nt ( l ± N t ) N − k ( N − k )! , (75)for some k = 2 , ..., N . Dividing by N the resulting function peaks for t = l ± , just as the memorykernels. At this t the expression in Equation (75) can be expanded in terms of order (cid:15) , similarly towhat has been done for the memory terms: N e − N ( N ) N − k ( N − k )! = 1 √ π N √ (cid:15) (1 − (cid:15) + O ( (cid:15) )) . (76)We find that the noise terms decay an order √ (cid:15) faster than the memory terms (at first order). Notethat this is an effect occurring at one time only as a remnant of the initial conditions. At all othertimes the noise term does not have any effect since Equation (75) approaches zero faster than (cid:15) for all other times. Therefore we neglect the noise term in the final equations of applying the MZformalism to the AMO model.To arrive at the final reduced model for the AMO we divide Equation (30) by N and use Equation(72) to write down expressions for the memory terms. We arrive at (cid:15) d T d t = − a T ( t ) − b T ( t ) + A l + e − αl + T (cid:0) t − l + (cid:1) + B l + e − αl + T (cid:0) t − l + (cid:1) + A − l − e − αl − T (cid:0) t − l − (cid:1) + B − l − e − αl − T (cid:0) t − l − (cid:1) + (cid:15)f (cid:15) ( t ) + O ( (cid:15) ) ,(cid:15) d T d t = − a T ( t ) − b T ( t ) + A l + e − αl + T (cid:0) t − l + (cid:1) + B l + e − αl + T (cid:0) t − l + (cid:1) + A − l − e − αl − T (cid:0) t − l − (cid:1) + B − l − e − αl − T (cid:0) t − l − (cid:1) + (cid:15)f (cid:15) ( t ) + O ( (cid:15) ) , (77)26here f (cid:15) ( t ) = − αT ( t ) + A l + e − αl + g (cid:15) + (cid:0) T (cid:1) + B l + e − αl + g (cid:15) + (cid:0) T (cid:1) + A − l − e − αl − g (cid:15) − (cid:0) T (cid:1) + B − l − e − αl − g (cid:15) − (cid:0) T (cid:1) ,f (cid:15) ( t ) = − αT ( t ) + A l + e − αl + g (cid:15) + (cid:0) T (cid:1) + B l + e − αl + g (cid:15) + (cid:0) T (cid:1) + A − l − e − αl − g (cid:15) − (cid:0) T (cid:1) + B − l − e − αl − g (cid:15) − (cid:0) T (cid:1) , (78)for g (cid:15) ± (cid:0) T (cid:1) = 12 l ± (cid:16)(cid:0) ( l ± + α ) − l ± (cid:1) T (cid:0) t − l ± (cid:1) + 2( l ± + α ) T (cid:48) (cid:0) t − l ± (cid:1) + T (cid:48)(cid:48) (cid:0) t − l ± (cid:1)(cid:17) ..