An eddifying Stommel model: Fast eddy effects in a two-box ocean

TL;DR

This paper develops a stochastic two-box ocean model incorporating fast eddy effects, comparing non-Gaussian, Gaussian, and averaged approximations to understand their impact on variability and equilibrium transitions.

Contribution

It introduces a novel stochastic model with non-Gaussian eddy forcing and compares it to Gaussian and averaged models using fast averaging and homogenization techniques.

Findings

Averaged model underestimates variability in single equilibrium regime.

Gaussian model captures second-order statistics but not skew or rare events.

All models show similar stationary distributions in multi-equilibrium regime, but jump rates differ.

Abstract

A system of stochastic differential equations is formulated describing the heat and salt content of a two-box ocean. Variability in the heat and salt content and in the thermohaline circulation between the boxes is driven by fast Gaussian atmospheric forcing and by ocean-intrinsic, eddy-driven variability. The eddy forcing of the slow dynamics takes the form of a colored, non-Gaussian noise. The qualitative effects of this non-Gaussianity are investigated by comparing to two approximate models: one that includes only the mean eddy effects (the `averaged model'), and one that includes an additional Gaussian white-noise approximation of the eddy effects (the `Gaussian model'). Both of these approximate models are derived using the methods of fast averaging and homogenization. In the parameter regime where the dynamics has a single stable equilibrium the averaged model has too little variability. The Gaussian model has accurate second-order statistics, but incorrect skew and rare-event probabilities. In the parameter regime where the dynamics has two stable equilibria the eddy noise is much smaller than the atmospheric noise. The averaged, Gaussian, and non-Gaussian models all have similar stationary distributions, but the jump rates between equilibria are too small for the averaged and Gaussian models.

Equations

$$\rho ={\rho }_0[1+{\alpha }_S(S-S_0)-{\alpha }_T(T-T_0)]$$

$$\frac{\text{d}{T}_1}{\text{d}t}=-\frac{1}{{\tau }_T}(T_1-T_1^{\ast })-\frac{F_q}{{\rho }_0{c}_p{V}_1},\frac{\text{d}{T}_2}{\text{d}t}=-\frac{1}{{\tau }_T}(T_2-T_2^{\ast })+\frac{F_q}{{\rho }_0{c}_p{V}_2}$$

$$\frac{\text{d}{S}_1}{\text{d}t}=\frac{1}{2}F(t)-F_S,\frac{\text{d}{S}_2}{\text{d}t}=-\frac{1}{2}F(t)+F_S$$

$$\frac{\text{d}\Delta T}{\text{d}t}=-\frac{1}{{\tau }_T}(\Delta T-\Delta T^{\ast })-2\left[\frac{1}{{\rho }_0{c}_p{V}_1}+\frac{1}{{\rho }_0{c}_p{V}_2}\right]F_T$$

$$\frac{\text{d}\Delta S}{\text{d}t}=F(t)-2F_S.$$

$$\text{d}\Delta T=\left[-\frac{1}{{\tau }_T}(\Delta T-\overline{\Delta T^{\ast }})-2\left[\frac{1}{{\rho }_0{c}_p{V}_1}+\frac{1}{{\rho }_0{c}_p{V}_2}\right]F_T\right]\text{d}t+\frac{{\sigma }_{\Delta T}}{\sqrt{{\tau }_T}}\text{d}{W}_{\Delta T}$$

$$\text{d}\Delta S=\left[\overline{F}-2F_S\right]\text{d}t+\frac{{\sigma }_S}{\sqrt{{\tau }_d}}\text{d}{W}_S.$$

$$2\left[\frac{1}{{\rho }_0{c}_p{V}_1}+\frac{1}{{\rho }_0{c}_p{V}_2}\right]F_T=\left(\frac{1}{{\tau }_d}+\frac{1}{{\tau }_a{\rho }_0{\alpha }_T\overline{\Delta T^{\ast }}}|\Delta \rho |\right)\Delta T$$

$$2F_S=\left(\frac{1}{{\tau }_d}+\frac{1}{{\tau }_a{\rho }_0{\alpha }_T\overline{\Delta T^{\ast }}}|\Delta \rho |\right)\Delta S$$

$$\Delta \rho ={\rho }_0[{\alpha }_S\Delta S-{\alpha }_T\Delta T]$$

$$2\left[\frac{1}{{\rho }_0{c}_p{V}_1}+\frac{1}{{\rho }_0{c}_p{V}_2}\right]F_T=\left(\frac{1}{{\tau }_d}+\frac{1}{{\tau }_a({\rho }_0{\alpha }_T\overline{\Delta T^{\ast }}{)}^2}\Delta {\rho }^2\right)\Delta T$$

$$2F_S=\left(\frac{1}{{\tau }_d}+\frac{1}{{\tau }_a({\rho }_0{\alpha }_T\overline{\Delta T^{\ast }}{)}^2}\Delta {\rho }^2\right)\Delta S.$$

$$2\left[\frac{1}{{\rho }_0{c}_p{V}_1}+\frac{1}{{\rho }_0{c}_p{V}_2}\right]F_T=\left(\frac{1}{{\tau }_d}+\frac{1}{{\tau }_a({\rho }_0{\alpha }_T\Delta T^{\ast }{)}^2}\Delta {\rho }^2\right)\Delta T+\left[\frac{1}{\ell }+\frac{1}{L_y-\ell }\right]vT$$

$$2F_S=\left(\frac{1}{{\tau }_d}+\frac{1}{{\tau }_a({\rho }_0{\alpha }_T\Delta T^{\ast }{)}^2}\Delta {\rho }^2\right)\Delta S+\left[\frac{1}{\ell }+\frac{1}{L_y-\ell }\right]vS.$$

$$\text{d}v=-\frac{1}{{\tau }_e}v\text{d}t+\sqrt{\frac{2}{{\tau }_e}}{\sigma }_v\text{d}{W}_v$$

$$\text{d}v$$

$$\frac{\text{d}T}{\text{d}t}=-\frac{T}{{\tau }_e}-v\frac{2\Delta T}{L_y},$$

$$\frac{\text{d}S}{\text{d}t}=-\frac{S}{{\tau }_e}-v\frac{2\Delta S}{L_y}.$$

$$\text{d}x$$

$$\text{d}y$$

$$\text{d}v$$

$$\text{d}T$$

$$\text{d}S$$

$${ϵ}_T=\frac{{\tau }_T}{{\tau }_d},ϵ=\frac{{\tau }_e}{{\tau }_d},P_a=\frac{{\tau }_d}{{\tau }_a},P_e=\frac{{\sigma }_v{\tau }_d}{L_y},P=\sqrt{ϵ}{P}_e.$$

$${ϵ}_T=\frac{1}{400},ϵ=\frac{1}{5000},P_a=6,P_e=80,{\sigma }_x=0.005,{\sigma }_y=\mathrm{0.15.}$$

$$\text{d}v$$

$$\text{d}T$$

$$\text{d}S$$

$$\left[\begin{array}[]{c c c}1&-P^{2}x&-P^{2}y\\ -P^{2}x&2P^{4}x^{2}+\sigma_{\epsilon}^{2}&2P^{4}xy\\ -P^{2}y&2P^{4}xy&2P^{4}y^{2}+\sigma_{\epsilon}^{2}\end{array}\right].$$

$$\text{d}x$$

$$\text{d}y$$

$$\text{\bf M}=-\frac{1}{\epsilon}\left[\begin{array}[]{c c c}1&0&0\\ 2P^{2}x&1&0\\ 2P^{2}y&0&1\end{array}\right],\;\;\text{\bf G}=\sqrt{\frac{2}{\epsilon}}\left[\begin{array}[]{c c c}1&0&0\\ 0&\sigma_{\epsilon}&0\\ 0&0&\sigma_{\epsilon}\end{array}\right]$$

$$\mathbf{Y}(\tau )=e^{\text{M}\tau }{\mathbf{Y}}_0+{\int }_0^{\tau }{e}^{\text{M}(\tau -s)}\text{G}\text{d}\mathbf{W}.$$

$$\bm{f}(x,y,\bm{Y})=\left(\begin{array}[]{c}vT+4P^{2}x\\ vS+4P^{2}y\end{array}\right).$$

Full text

In the parameter regime where the dynamics has a single stable equilibrium the averaged model has too little variability. The Gaussian model has accurate second-order statistics, but incorrect skew and rare-event probabilities. In the parameter regime where the dynamics has two stable equilibria the eddy noise is much smaller than the atmospheric noise. The averaged, Gaussian, and non-Gaussian models all have similar stationary distributions, but the jump rates between equilibria are too small for the averaged and Gaussian models.

1 Introduction

H. Stommel (1961) developed a conceptual model of the global ocean thermohaline circulation that consists of a system of ordinary differential equations modeling the heat and salt content of two containers (‘boxes’). One box models the equatorial ocean, and the other models the extra-tropical ocean. The boxes exchange heat and salt with each other and with the atmosphere. The rate of flow between the boxes is proportional to the density difference between the boxes, and a major result of Stommel’s investigation was that in some parameter regimes the system exhibits two equilibria: one analogous to the current climate, with dense cold water sinking at high latitudes, and one corresponding to a very different regime with dense salty water sinking in the equatorial ocean. In general, the goal of studies using extremely simplified models like Stommel’s is to observe and understand qualitative features that might inform and guide subsequent studies using more complete and more complex models. The qualitative predictions of Stommel’s model have since been verified using more complete ocean models, e.g. Rahmstorf (1995) and Deshayes et al. (2013).

The present investigation develops a model closely related to Stommel’s where the slow, density-driven exchange of heat and salt between the boxes is augmented by fast, non-Gaussian stochastic processes representing eddy-driven heat and salt transport. Eddies smaller than the grid scale of comprehensive numerical ocean (and atmosphere) models can have significant impacts on the global circulation, and modeling the impacts of these unresolved eddies is a topic of continuing research; Berner et al. (2016) contains a review of stochastic models of eddy effects from an operational modeling perspective. The second author recently proposed a non-Gaussian model of the heat and salt transport associated with unresolved ocean eddies (Grooms, 2016). In this model, the eddy velocity and density fields (the latter linearly related to temperature and salinity) are modeled as centered Gaussian random fields, and the transports are modeled as the product of eddy velocity and density. The product of centered, jointly Gaussian random variables has a distinctive, non-Gaussian probability density, with a logarithmic singularity at the origin and skewed, algebraically-modulated exponential decay in the tails. This non-Gaussian model is significantly different from recent Gaussian stochastic models of eddy transport, e.g. Andrejczuk et al. (2016); Williams et al. (2016) and Juricke et al. (2017). The present investigation is motivated by the desire to observe the qualitative effects of this kind of non-Gaussian transport in an extremely simple model, in particular by comparison to Gaussian stochastic models, with the expectation of informing future investigations using more complex models.

Several authors have developed stochastic versions of Stommel’s model to investigate the slow response of the ocean thermohaline circulation to fast atmospheric forcing, e.g. Cessi (1994); Vélez-Belchı et al. (2001); Monahan (2002); Monahan et al. (2002) and Monahan and Culina (2011). In these stochastic Stommel models the atmospheric heat and salt fluxes in Stommel’s model are replaced by Gaussian stochastic noise terms, resulting in a system of stochastic differential equations (SDEs). The model developed here attempts to understand a qualitatively different physical process: fast eddy transport. Since the eddies are typically faster than the global thermohaline circulation, the new model has the form of a slow-fast system, where eddy variables evolve on a fast time scale and converge towards a jointly Gaussian distribution conditioned on the slow variables. The slow variables (the heat and salt difference between the boxes) are impacted by quadratic products of fast variables modeling the fast eddy transport. The formal theory of fast averaging (Papanicolaou and Kohler, 1974; Pavliotis and Stuart, 2008; Freidlin and Wentzell, 2012), is used to generate approximate slow systems for comparison: one with a drift correction and one with both drift and diffusion corrections derived from the eddy dynamics. These approximate systems qualitatively represent more complete ocean models with, respectively, deterministic and Gaussian stochastic models of the eddy transport.

A new stochastic Stommel model is developed in §2. The two approximate models of the slow system are derived in §3. The numerical methods and experimental configuration are described in §4 and the results of these simulations are described in §5. A slightly different model with two stable equilibria is formulated and simulated in §6. The results and their implications are discussed in §7.

2 Formulating a Slow-Fast Two-Box Stochastic Ocean Model

Consider a domain $[0,L_{x}]\times[0,L_{y}]\times[0,H]$ representing an ocean basin, and let this domain be partitioned into two subdomains $[0,L_{x}]\times[0,\ell]\times[0,H]$ and $[0,L_{x}]\times[(\ell,L_{y}]\times[0,H]$ with volumes $V_{1}=L_{x}\ell H$ and $V_{2}=L_{x}(L_{y}-\ell)H$. The first box (index 1) will represent the equatorial side of the ocean basin, and the second will represent the poleward side. The domain is filled with a fluid whose density is related to its temperature and salinity via

[TABLE]

where $\rho_{0}=1029$ kg/m3 is a constant reference density, $T_{0}=5$ C and $S_{0}=35$ psu are a constant reference temperature and salinity (psu are practical salinity units; for the present purposes it is reasonable to use the simplification 1 psu = 1 g/kg), and $\alpha_{S}=7.5\times 10^{-4}$ psu*-1* and $\alpha_{T}=1.7\times 10^{-4}$ C*-1* are coefficients of haline and thermal expansion. The conservation equations for heat are in the form of a system of two differential equations

[TABLE]

where $T_{1}$ and $T_{2}$ are the mean temperature in each box, $\tau_{T}$ is the timescale of relaxation towards an externally-specified atmospheric temperature $T_{i}^{*}$, $c_{p}=4000$ J/kg is the heat capacity of seawater (e.g. $\rho_{0}c_{p}V_{1}T_{1}$ is the heat content of the equatorial box), and $F_{q}$ is the heat flux from the equatorial box to the poleward box. The total heat content $\rho_{0}c_{p}(V_{1}T_{1}+V_{2}T_{2})$ thus depends only on the external forcing.

Similarly, the conservation equations for salt are

[TABLE]

where $2F(t)$ is the external freshwater forcing in the equatorial box (e.g. rain, runoff, ice melt) and $F_{S}$ is the salt flux from the equatorial box to the poleward box. The external freshwater forcing is assumed not to change the net salt content, so that $S_{1}+S_{2}$ remains constant in time.

Following Stommel (1961), the heat and salt fluxes between the boxes are assumed to depend only on the temperature and salinity differences between the boxes. As a result, the temperature and salinity differences between the boxes decouple from the net heat and salt content. Defining $\Delta T=T_{1}-T_{2}$ and $\Delta S=S_{1}-S_{2}$,

[TABLE]

[TABLE]

Similar to Cessi (1994) and Vélez-Belchı et al. (2001), the atmospheric temperature difference $\Delta T^{*}$ and external freshwater forcing are here modeled as constant mean terms plus Gaussian white noise, leading to

[TABLE]

[TABLE]

The amplitude of the noise forcing is here related to $\tau_{d}$, a diffusive time scale, defined below.

In Stommel’s original model the fluxes between the boxes consists of diffusive fluxes proportional to the temperature and salinity differences, and advective fluxes associated with the large-scale ocean circulation whose rate is proportional to the magnitude of the density difference between the boxes

[TABLE]

[TABLE]

where $\tau_{d}$ is the time scale of diffusive transport, $\tau_{a}$ is the time scale of advective transport, and

[TABLE]

is the density difference between the boxes. Cessi (1994) used a smoother formulation, which does not qualitatively change the results

[TABLE]

[TABLE]

The novel contribution to the model made here consists of the addition of fast variables crudely representing eddy velocity $v$, temperature $T$, and salinity $S$ anomalies at the interface between the boxes. The eddy-induced fluxes between the boxes will be modeled as an addition to the slow diffusive and advective fluxes

[TABLE]

[TABLE]

The prefactors of $\ell^{-1}+(L_{y}-\ell)^{-1}$ account for the fact that the boxes need not have equal volume, and that total heat and salt need to be conserved. For simplicity, only $\ell=L_{y}/2$ is considered from here on.

In general the flux between the boxes should be described by $\int_{0}^{L_{x}}\int_{0}^{H}vT$d$z$d$x$ where $v$ and $T$ are evaluated at $y=\ell$. Our formulation amounts to a severe simplification that ignores the spatial structure of the eddy velocity and temperature perturbations between the boxes, and considers them only as zero-mean jointly-Gaussian variables. This level of simplification is consistent with the simplification of the ocean to two well-mixed boxes in the original Stommel model, and is guided by the desire to investigate the qualitative effects of Gaussian-product noise, since eddy noise with this structure was recently observed by Grooms (2016).

The fast eddy velocity will be modeled as an Ornstein-Uhlenbeck process

[TABLE]

where $\tau_{e}$ is the eddy time scale and $\sigma_{v}$ is the eddy velocity scale, chosen to be $15$ days and $10$ cm/s, respectively Stammer (1997). The eddy velocity can be thought of as being set by wind-driven processes independent of the density difference between the boxes. This is a simplification of the more complex reality where eddy kinetic energy depends also on the large-scale density gradient. The following model of the eddy dynamics is perhaps more qualitatively appropriate

[TABLE]

where $\mu>0$ is a parameter representing the sensitivity of the eddy variance to the large-scale density gradient. This model is not pursued further here, in part because of the difficulties in guaranteeing its ergodicity and in finding a robust numerical method for its solution.

The eddy temperature and salinity anomalies will be modeled as resulting from eddy transport across the large-scale gradients

[TABLE]

[TABLE]

The relaxation towards zero on a time scale of $\tau_{e}$ qualitatively represents the full range of dissipative processes acting on temperature and salinity anomalies: cascade towards small scales and eventual diffusion, and atmospheric damping of thermal anomalies, etc. The time scale $\tau_{e}$ should not be associated with any particular physical process, but instead guarantees decorrelation of eddy anomalies on the time scale $\tau_{e}$.

The governing equations are nondimensionalized using the diffusive time scale, the external constant atmospheric temperature difference $\overline{\Delta T^{*}}\approx 20$ C for large-scale temperature, and the convenient salinity scale $\alpha_{T}\overline{\Delta T^{*}}/\alpha_{S}\approx$ 4.5 psu for large-scale salinity. It will be convenient to scale the eddy temperature and salinity variables differently; specifically, $T$ will have dimensions $\overline{\Delta T^{*}}L_{y}/(\sigma_{v}\tau_{d})$ and $S$ will have dimensions $\alpha_{T}\overline{\Delta T^{*}}L_{y}/(\alpha_{S}\sigma_{v}\tau_{d})$. The reason for this unexpected scaling will be commented on shortly.

Following traditional notation, the nondimensional temperature difference will be denoted $x$ and the nondimensional salt difference will be denoted $y$; risking confusion, the nondimensional eddy variables will have the same notation as their dimensional counterparts. The complete nondimensional system is therefore

[TABLE]

where

[TABLE]

$P_{a}$ and $P_{e}$ are Péclet numbers comparing the time scales of large-scale advective transport and fast eddy transport to the time scale of diffusion, respectively.

The following parameter estimates are drawn from Cessi (1994) and Vélez-Belchı et al. (2001), and are consistent with the more recent observational analysis of Schmitt (2008). The diffusive time scale $\tau_{d}$ is approximately 219 years, and the time scale of large scale advection $\tau_{a}$ is approximately 35 years. Cessi estimates $\tau_{T}$ to be 25 days Cessi (1994), but Vélez-Belchi et al. argue convincingly that large-scale temperature anomalies are damped on a slower time scale of approximately $220$ days. Vélez-Belchi et al. used salinity noise whose nondimensional amplitude is here $\sigma_{y}=0.15$, and assuming that fast atmospheric temperature fluctuations lead to perturbations on the order of 0.07 C implies nondimensional thermal noise has amplitude $\sigma_{x}=0.005$. Finally, using a length scale appropriate to the global oceans $L_{y}\approx 8,250$ km leads to the following set of parameters which are adopted for the remainder of the investigation

[TABLE]

The reason for scaling $S$ and $T$ differently from $\Delta S$ and $\Delta T$ should now be clear: $2P_{e}^{2}$ is the same order of magnitude as $\epsilon^{-1}$, implying that both terms in the evolution equations for $S$ and $T$ are of comparable magnitude.

For the parameters (2) the system (1) has three equilibria, two of which are stable. The equilibria all have $v,T,S=0$, and the stable equilibria occur at $(x,y)\approx(0.989,0.22)$ and $(x,y)\approx(0.998,1.00)$. In the absence of eddy dynamics, one would expect small atmospheric noise to lead to jumping between the two stable equilibria of the system; this was the focus of Cessi (1994); Monahan (2002); Monahan et al. (2002) and Monahan and Culina (2011). The existence of stable equilibria is intrinsically tied to the nonlinear terms that model slow advective exchange between the boxes. As the exchange between the boxes becomes dominated by diffusion instead of advection ($P_{a}\to 0$) one of the stable equilibria disappears in a reverse saddle-node bifurcation leaving a single stable equilibrium.

Equations (1d) and (1e) lack noise terms, implying that the classical conditions for ergodicity (Khasminskii, 2012) do not apply. Conditions for ergodicity of this type of system of SDEs can be found in Mattingly et al. (2002). The first condition is that there is an inner-product norm $\|\bm{\cdot}\|$ such that $\langle\bm{u},\bm{F}(\bm{u})\rangle\leq\alpha-\beta\|\bm{u}\|^{2}$ for some $\alpha,\beta>0$ where $\bm{u}$ is a vector containing the dependent variables and $\bm{F}(\bm{u})$ is the drift. It is straightforward to verify that $\|\bm{u}\|^{2}=x^{2}+y^{2}+\epsilon v^{2}+(2\epsilon/P^{2})(T^{2}+S^{2})$ satisfies this condition. The second condition is that the vectors $\left\{\bm{\rho}_{i},[[\bm{F},\bm{\rho}_{j}],\bm{\rho}_{k}]\right\}$ span $\mathbb{R}^{5}$ where $\bm{\rho}_{i}$, $i=1,\,2,\,3$ are the columns of the diffusion matrix, which are here proportional to the first three standard basis vectors, and $[\mathbf{\cdot},\mathbf{\cdot}]$ is a Lie bracket. Since $[[\bm{F},\bm{\rho}_{1}],\bm{\rho}_{3}]$ and $[[\bm{F},\bm{\rho}_{2}],\bm{\rho}_{3}]$ are proportional to the fourth and fifth standard basis vectors, respectively, the system satisfies the conditions of Mattingly et al. (2002) for ergodicity.

3 Two Approximate Slow Models

In this section two systems of SDEs are derived approximating the evolution of the slow variables $x$ and $y$ in (1). The system of SDEs (1) with parameters (2) has three time scales since $\epsilon<\epsilon_{T}\ll 1$: $x$ evolves significantly more quickly than $y$, yet slower than the eddy variables $v$, $T$, and $S$. Many previous investigations (which lacked the eddy variables) accounted for the scale separation somewhat crudely by setting $x=1$, and focused on the dynamics of the slowest variable $y$, e.g. Cessi (1994); Monahan (2002); Monahan et al. (2002), and Monahan et al. (2008). The analysis of Monahan and Culina (2011) is more careful, employing the same methods used here but for the system without eddy variables and in the limit $\epsilon_{T}\to 0$. This section considers the limit $\epsilon\to 0$ while holding $\epsilon_{T}$ fixed.

The two approximate models are derived using standard approximations for slow-fast systems (Papanicolaou and Kohler, 1974; Pavliotis and Stuart, 2008; Freidlin and Wentzell, 2012). The presentation here follows the convenient review found in Bouchet et al. (2016); the formulas are derived in a straightforward manner using formal asymptotic methods applied to the backwards Kolmogorov equation for the system (for details, see the appendices of Bouchet et al. (2016)).

The first approximation is derived via simple averaging. In the limit $\epsilon\to 0$ the eddy variables are well approximated as solutions to (1c)–(1e) with $x$ and $y$ considered constant. Curiously, although the full system (1) has a smooth invariant measure the system (1c)–(1e) does not: the long-time limiting distribution of $v$, $T$, and $S$ is jointly Gaussian with a singular covariance matrix. In light of this, the following noise-augmented system is considered instead

[TABLE]

and the limit $\sigma_{\epsilon}\to 0$ is taken after the fact.

The invariant measure of (3) is Gaussian with zero mean and covariance

[TABLE]

The averages of the terms $vT$ and $vS$ in the slow equations with respect to the invariant measure of the fast system are simply $-P^{2}x$ and $-P^{2}y$, respectively. It is worth noting that these values are independent of the auxiliary noise amplitude $\sigma_{\epsilon}$. Inserting these into the slow equations leads to the following approximate model

Deterministic Approximation

[TABLE]

The model (5) is referred to as the ‘deterministic’ or ‘averaged’ approximation since it models the eddy terms $vT$ and $vS$ as deterministic functions of $x$ and $y$. It is straightforward to verify that this model is ergodic under the classical conditions of Khasminskii (2012).

As described in Bouchet et al. (2016), one can derive equations that approximate the variations of the true solution to (1) around the solution of the approximate model (5). Combining the equations for the variations with the deterministic approximation leads to further corrections in both the drift and diffusion, of order $\epsilon$ and $\sqrt{\epsilon}$, respectively. The drift correction is significantly smaller than the leading-order drift. But the leading-order diffusion terms in the $x$ and $y$ equations are of order $\approx 0.1$, and corrections of order $\sqrt{\epsilon}$ may be of comparable magnitude.

In order to compute the diffusion corrections, it is convenient to define some notation. Let $\bm{Y}=(v,T,S)^{T}$ denote the solution to the noise-augmented system (3). Define constant matrices

[TABLE]

such that the fast system (3) may be written d$\bm{Y}=$ M$\bm{Y}+$ Gd$\bm{W}$. The solution is thus

[TABLE]

The deviations of the eddy terms $vT$ and $vS$ from their conditional means are denoted

[TABLE]

According to Bouchet et al. (2016), the diffusion-corrected model for the slow variables has the form

[TABLE]

where the matrix

[TABLE]

is any square root of the following symmetric positive definite matrix

[TABLE]

The matrix C is the integral of the time-lagged auto-covariance of $\bm{f}$ with $x$ and $y$ considered constant. In the above expression, $\mathbb{E}^{\bm{Y}(\tau)}$ denotes the expectation on $\bm{Y}(\tau)$ conditioned on the initial condition $Y_{0}$; the distribution is Gaussian with mean and covariance implied by (6). $\mathbb{E}^{\bm{Y}_{0}}$ denotes expectation on $\bm{Y}_{0}$ whose distribution is the stationary distribution of the fast process, in this case a zero-mean Gaussian with covariance (4). The calculation for the system under consideration here is particularly straightforward since it requires only higher moments of jointly-Gaussian variables. The matrix C is found to have the form

[TABLE]

In this case (unlike the leading-order drift term) the limit $\sigma_{\epsilon}\to 0$ is singular in the sense that the matrix C becomes positive semi-definite. Nevertheless, a square root matrix A exists; in the limit $\sigma_{\epsilon}\to 0$ it has the form

[TABLE]

The model for the slow variables with leading-order drift and diffusion corrections (but ignoring the order-$\epsilon$ drift correction) is thus

Gaussian Stochastic Approximation

[TABLE]

For $x\approx 1$ the noise amplitude associated with the eddies is $\approx 0.16$, which is slightly larger than the ‘atmospheric’ noise $\sigma_{\epsilon}/\sqrt{\epsilon_{T}}=0.1$. The order-$\epsilon$ drift corrections have also been calculated, but they are small in comparison with the leading-order terms, and have been left out of the model for simplicity. This system of SDEs is interpreted in the Ito sense; while the drift corrections in slow-fast systems with one slow degree of freedom can be interpreted as a correction from Stratonovich to Ito, this is no longer generally true in systems with multiple slow degrees of freedom (Pavliotis and Stuart, 2008; Freidlin and Wentzell, 2012). It is straightforward to verify that this model is ergodic under the classical conditions of Khasminskii (2012).

It is interesting to note that the Gaussian stochastic model replaces the eddy terms $vT$ and $vS$ by $-4P^{2}x($d$t+\sqrt{5\epsilon}$d$\hat{W})$ and $-4P^{2}y($d$t+\sqrt{5\epsilon}$d$\hat{W})$. This form of subgrid-scale parameterization is qualitatively the same as that proposed in Buizza et al. (1999), where it was proposed to multiply a deterministic parameterization (here $-4P^{2}x$) by a stochastic process (here $1+\sqrt{5\epsilon}\dot{W}$). This style of stochastic parameterization has been widely used in atmospheric models (Berner et al. (2016) provides a review), and much has been made of the role of multiplicative noise by, e.g. Sura et al. (2005). The above derivation gives an example where this style of ad hoc parameterization is rigorously justified, though multiplicative noise with a linear coefficient is certainly not the universal form of eddy-induced noise (see e.g. Monahan and Culina, 2011, for a counterexample).

Recall that for the parameters (2) the system (1) has only three equilibria, two of which are stable. The equilibria all have $v,T,S=0$, and the stable equilibria occur at $(x,y)\approx(0.989,0.22)$ and $(x,y)\approx(0.998,1.00)$. The deterministic and Gaussian stochastic models have the same drift, which has only one equilibrium at $(x,y)\approx(0.974,0.093)$. As will be verified by the results in §5, the inclusion of nonlinear eddy effects completely changes the regime of the ocean model from a regime of multiple equilibria to a regime with a single stable equilibrium.

The averaged drift has a single stable equilibrium for all $P$ greater than approximately 0.117; below this value the drift undergoes a saddle-node bifurcation that creates a pair of equilibria near $x=1$ and $y=1$. To achieve such small values of $P$ would require reducing the eddy velocity scale from 10 cm/s to 1 cm/s, which is unrealistically small. The approximate models derived in this section show that the mean effect of eddies is linear and diffusive. Since a linear diffusive effect is already present in the equations (the terms $-x$ and $-y$ in (1a) and (1b)), the mean eddy effect could be viewed as a double-counting of eddy-induced diffusive exchange between the boxes. This can be rectified by eliminating the mean diffusion terms, and such a model is formulated and studied in §6. By avoiding a double-counting of diffusive exchange, the model in §6 allows multiple equilibria with small, yet realistic eddy amplitudes.

4 Numerical Methods

Numerical methods are needed to compare the qualitative behavior of the three models (1), (5), and (7) . Many methods are derived based on the assumption that the drift is globally-Lipschitz (Kloeden and Platen, 1992), which is not the case here. Several more recent investigations have analyzed numerical methods for SDEs whose drift satisfies a one-sided Lipschitz condition (e.g. Higham et al. (2002) and Mao and Szpruch (2013)), but none of the models in consideration here satisfy such a condition. A method appropriate to polynomial drifts is derived by Lamba et al. (2007), but their analysis requires an invertible diffusion matrix, which the model (1) does not have. The Euler-Maruyama method may be appropriate, but is known to behave poorly in problems with polynomial drift (Mattingly et al., 2002; Hutzenthaler et al., 2011). In light of this, the ‘backward Euler’ (BE) method is used here for all three models. For a general system of SDEs of the form

[TABLE]

the BE method takes the following form

[TABLE]

where $\Delta t$ is the time step. In every simulation presented here $\Delta t=2\times 10^{-6}$, which is significantly smaller than the smallest time scale of the system $\epsilon=2\times 10^{-4}$. Mattingly et al. (2002) prove that the method is ergodic (for sufficiently small $\Delta t$) and that the invariant measure of the numerical method converges to that of the SDE as $\Delta t\to 0$. Though the analysis of Mattingly et al. (2002) focuses on models with additive noise, the BE method is nevertheless applied here to the model (7) with multiplicative noise.

For the model (1), a two-step process is used to generate solutions of the nonlinear system of equations (8). First, an asymptotic approximation in the limit $\Delta t\to 0$ is computed that has the form $\bm{X}_{*}=\bm{X}_{n}+\mathbf{\Sigma}(\bm{X}_{n})\Delta\bm{W}_{n}+\mathcal{O}(\Delta t)$; this approximation is followed by a single Newton step. For the systems (5) and (7), approximate solutions to the the nonlinear systems were generated using 10 fixed-point iterations started at $\bm{X}_{n}$. Given the small step size, the resulting approximations solve their respective nonlinear systems with high accuracy; the residuals are typically on the order of $10^{-11}$.

5 Results

5.1 Climatology

A suite of 10,000 independent simulations was run starting from $x,\,y,\,v,\,T,\,S=0$ at $t=0$. Data were saved for the time interval $t\in[4,10]$, saving every 100${}^{\text{th}}$ time step for a spacing of $2\times 10^{-4}$. The mean and covariance appeared to have stabilized by $t=4$, suggesting that the data in $t\in[4,10]$ represents the stationary climatological distribution of the system. Recalling that the dimensional time unit is 219 years, this amounts to 1,314 years of data saved approximately twice per month. The three models all have the same mean of $(x,y)\approx(0.974,0.094)$, which is very close to the equilibrium of the deterministic and Gaussian stochastic models at $(0.974,0.093)$. All three models have the same marginal standard deviation of $y$ approximately equal to $0.034$. This can be explained by the fact that the amplitude of the eddy noise in the $y$ equation is estimated in the Gaussian stochastic model to be $4\sqrt{5\epsilon}P^{2}y\approx 0.015$ for $y=0.093$, which is much less than the atmospheric noise with amplitude $\sigma_{y}=0.15$.

The climatological distributions of the models differ in other respects. For example, the marginal standard deviation of $x$ is 0.0063, 0.0035, and 0.0065 in the full, deterministic, and Gaussian stochastic models, respectively. The eddy noise in the $x$ equation is of comparable size to the atmospheric noise, and has a significant impact on the variability; the deterministic model lacks this eddy noise, and has too little variability. The lack of eddy noise in the $x$ equation of the deterministic model also leads to an overestimate of the correlation between $x$ and $y$: the full and Gaussian stochastic models have correlations $0.15$ and $0.14$, respectively, while the deterministic model has correlation $0.23$. The most-probable values of the distributions are $(x,y)\approx(0.976,0.092)$ for the full model, $(0.974,0.091)$ for the deterministic model, and $(0.973,0.093)$ for the Gaussian stochastic model; the differences in the $y$ value are negligible, but the differences in the $x$ value are up to half of a standard deviation.

Time-lagged correlation functions were computed, for example Corr$[x(t),x(t+\tau)]=C(\tau)$ (stationarity is assumed). The correlation functions are all very similar across the models (not shown). The correlation functions all decay monotonically to zero, so it is natural to define a decorrelation time by $\int_{0}^{\infty}C(\tau)\text{d}\tau$. The correlation functions for $y$ in all three models are very similar, with decorrelation time approximately 22 years. The correlation function for $x$ exhibits similar rapid initial decay in all three models. The correlation function for $x$ in the deterministic model has a long tail, with larger long-lag correlations than the other two models, leading to a decorrelation time of 1.6 years, which is longer than the decorrelation times of the full model and Gaussian stochastic model, both of which are approximately 1 year.

A simple binning procedure was used to generate approximations to the climatological probability density function (pdf) for each model; results are shown in Fig. 1, with panels (a)–(c) presenting the full model, deterministic model, and Gaussian stochastic model, respectively. It has already been noted that the three models have the same marginal variance for $y$, and indeed the range of $y$ in the three models is quite similar. The deterministic model is clearly under-dispersed with respect to $x$. The climatological distribution of the Gaussian stochastic model has a more-accurate core, but is not skewed in the same way as the full model.

It is possible that minor deficiencies near the core of the distribution could be corrected by adding order-$\epsilon$ corrections to the drift of the Gaussian stochastic model, but the results of Bouchet et al. (2016) indicate that such corrections will not generate correct rare-event probabilities even in the limit $\epsilon\to 0$. To emphasize differences in the rare event probabilities, the probabilities of $x\leq 0.96$ and $x\geq 0.985$ were calculated for the three models (these $x$ values are indicated by vertical lines in Fig. 1). The small-event probabilities are 0.039 for the full model, less than $10^{-4}$ for the deterministic model, and $0.022$ for the Gaussian stochastic model. The large-event probabilities are 0.016 for the full model, less than $10^{-3}$ for the deterministic model, and $0.048$ for the Gaussian stochastic model. Not surprisingly, the deterministic approximation has too-small rare event probabilities. The Gaussian-stochastic model is more accurate, but is still incorrect by more than a factor of 2 in each case.

In order to verify that the system does not remain near the stable equilibrium of (1) at $(x,y)\approx(1,1)$, a set of 1,000 simulations of (1) was run with initial condition $(x,y,v,T,S)=(1,1,0,0,0)$. These simulations were again run for the interval $t\in[0,10]$, saving the output from $t\in[4,10]$. The stationary distribution did not display a secondary peak near $(1,1)$, indicating that the two stable equilibria of the full model are largely irrelevant to the dynamics of the system.

5.2 Rare event forecasting

The previous section examined only the stationary climatological distributions of the three models. Within a climate prediction scenario, short-term behavior is also important. Given that the climatological distributions differ mainly in their rare event probabilities, a separate set of experiments was used to investigate the ability of the models to predict rare events over a shorter time interval. The goal was to test how accurately the approximate models forecast the probability of the unusually large and small $x$ values over a range of forecast lead times. Two trajectories of the system (1) were selected out of the 10,000 discussed above: one with a value of $x\leq 0.96$ and one with a value of $x\geq 0.985$. These trajectories are shown in Fig. 2 panels (a) and (b). Note that the large-$x$ trajectory passes the threshold of $0.985$ approximately half a year before the final time, whereas the small-$x$ trajectory crosses the $0.96$ threshold only at the last time step. For a single lead time, e.g. 2 years before the rare event, ensembles of 10,000 independent forecasts for all three models were initialized from the true trajectory. At the time of the rare event these ensembles were used to estimate the probability of $x\leq 0.96$ (or $x\geq 0.985$). The probabilities shown in Fig. 2c correspond to the large-$x$ trajectory, and those in Fig. 2c correspond to the small-$x$ trajectory. Since the large-$x$ trajectory crosses the threshold nearly half a year before the final time, all 10,000 of the forecasts initialized at any lead time less than half a year in advance are already above threshold; nevertheless, the probability at the final time is less than one because many of the trajectories cross the threshold back towards smaller values of $x$.

In both cases the forecast by the deterministic model is severely untrustworthy, assigning extremely low rare-event probability even when initialized close to or even beyond the threshold. The rare-event probability forecast by the Gaussian stochastic model, in contrast, begins to increase from its climatological value at approximately the same time that the true forecast probability begins to increase, between 0.8 and 0.6 years in advance for the large-$x$ event and around 0.3 years in advance for the small-$x$ event. Although the actual probability assigned by the Gaussian stochastic model at relatively long lead times is incorrect, the fact that it begins to increase at the right time could still be used qualitatively to predict whether the model is getting close to a rare event. Once the probability of a rare event increases past about 20%, the Gaussian stochastic model uniformly under-predicts the correct probability, despite having over-predicted the climatological probability for $x>0.985$. For example, with a lead time of about 2.5 months the Gaussian stochastic model predicts the large-$x$ event with probability only 53% while the true probability is in fact 67%; with a lead time of half a month the Gaussian stochastic model predicts the small-$x$ event with probability only 21% while the true probability is 61%.

In summary, the deterministic model is essentially useless for rare-event forecasting, while the Gaussian stochastic model is only qualitatively useful, predicting whether a rare event is more likely but not with a robust uncertainty estimate.

6 A model without mean diffusion

As noted at the end of §3, the averaged effect of the eddies is linear and diffusive. Linear diffusive terms are already included in the budgets of heat and salt, with the result that the averaged models have only one stable equilibrium unless the eddies are assumed to be extremely weak, with velocities on the order of 1 cm/s. If one assumes that linear diffusive exchange between the boxes is entirely eddy-driven then one can drop the mean diffusion terms from the governing equations of the full model, i.e. equations (1a) and (1b) are changed to

[TABLE]

and

[TABLE]

respectively. The eddy reductions proceed as before, so that the $-x$ and $-y$ terms are similarly dropped from the deterministic (5) and Gaussian (7) models. The resulting model is much more amenable to multiple equilibria. For $P$ greater than about 0.514 there is a single stable equilibrium with $x\approx 1$ and $y\approx 0.25$. Below this value of $P$ the system undergoes a saddle-node bifurcation that creates a pair of equilibria near $(x,y)=(1,1)$; the saddle then moves down towards the original equilibrium, which it joins in a reverse saddle-node bifurcation at $P$ approximately 0.301, below which there remains only a single equilibrium. We investigate the system at a value of $P_{e}=32$, i.e. $P\approx 0.45$, where there are three equilibria: a stable one at $(.99,.24)$, a saddle at $(1.00,.65)$, and another stable one at $(1.00,1.11)$.

6.1 Ergodicity

Recall that there are two conditions for ergodicity of hypoelliptic SDEs in Mattingly et al. (2002). The first condition is that there is an inner-product norm $\|\bm{\cdot}\|$ such that $\langle\bm{u},\bm{F}(\bm{u})\rangle\leq\alpha-\beta\|\bm{u}\|^{2}$ for some $\alpha,\beta>0$ where $\bm{u}$ is a vector containing the dependent variables and $\bm{F}(\bm{u})$ is the drift. The second condition is that the vectors $\left\{\bm{\rho}_{i},[[\bm{F},\bm{\rho}_{j}],\bm{\rho}_{k}]\right\}$ span $\mathbb{R}^{5}$ where $\bm{\rho}_{i}$, $i=1,\,2,\,3$ are the columns of the diffusion matrix, and $[\mathbf{\cdot},\mathbf{\cdot}]$ is a Lie bracket. It is straightforward to verify that the second condition is met in this model in the same way that it is met in the original model (1).

The first condition is more difficult. We will use the inner product $\langle\bm{u},\bm{v}\rangle=u_{1}v_{1}+u_{2}v_{2}+\epsilon u_{3}v_{3}+(2\epsilon/P^{2})(u_{4}v_{4}+u_{5}v_{5})$, so we must show that there are $\alpha,\beta>0$ such that

[TABLE]

i.e.

[TABLE]

The terms involving the eddy variables ($v$, $T$, and $S$) will clearly pose no problem provided that $\beta<\epsilon^{-1}$. It therefore remains to see whether one can choose $\alpha,\beta$ such that

[TABLE]

Consider the behavior along a line through the origin in the $(x,y)$ plane: along any line except $y=x$ the function is a quartic polynomial that can be made negative by choosing $\alpha$ sufficiently large. Along the line $y=x$ the condition reduces to

[TABLE]

As long as $\beta<\epsilon_{T}/2$ it will be possible to choose $\alpha$ sufficiently large that this condition is met. The model without mean diffusion terms is therefore still ergodic. Ergodicity is important because it implies that there is a single climatological distribution independent of the initial condition; the conditions of Mattingly et al. (2002) further guarantee that the distribution collapses exponentially quickly towards the climatological distribution.

6.2 Numerical experiments

Ensemble simulations for the three models without mean diffusion were run with 1000 ensemble members each. The deterministic and Gaussian approximate models were initialized with $x=1$, $y=0.6$, while the full model was initialized with $x=1$, $y=0.65$, and $v,T,S=0$. After a burn-in of 4 nondimensional time units, the simulations were run for 500 more time units, i.e. about 110,000 years. Although the models are geometrically ergodic, with distributions collapsing exponentially quickly towards the invariant distribution, this was not enough time for the approximate models to reach the invariant distribution. These models were then extended for a further 500 time units, during which time their distributions converged. The full model was initialized closer to the saddle point, so its distribution converged within the first 504 time units.

The climatological distributions of the slow variables are shown in Fig. 3. Panel (a) shows the marginal $y$ distributions of the three models, while panels (b)–(d) show the joint $(x,y)$ distributions. The three models are remarkably similar. Though $\epsilon$ is the same as in the previous case, $P$ is smaller. The diffusion correction in the $x$ equation of the Gaussian-stochastic approximation has amplitude $4\sqrt{5\epsilon}P^{2}x\approx 0.026$ which is smaller than the atmospheric noise amplitude $\sigma_{x}/\sqrt{\epsilon_{T}}=0.1$; the atmospheric noise similarly dominates the $y$ equation. As a result, the effects of eddy noise are not seen in the equilibrium distributions of the three models.

The noise levels are low enough that the system trajectories make rare transitions between the neighborhoods of the two stable equilibria; Fig. 4 panel (a) shows a system trajectory $y$ from the full model that jumps between regimes. The rates and paths of these transitions are the subject of large deviation theory (Freidlin and Wentzell, 2012). The methods of Bouchet et al. (2016) to analyze the transitions do not seem to apply directly here because of the inclusion of noise forcing in the slow dynamics. In any case, it is not difficult to estimate the transition probabilities from simulations. For practical purposes it was convenient to estimate the following probabilities $p_{01}(\tau)=P(y(t+\tau)>0.8\,|\,y(t)<0.5)$ and $p_{10}(\tau)=P(y(t+\tau)<0.5\,|\,y(t)>0.8)$. These transition probabilities are plotted for the three models in panels (b) and (c), respectively. The effects of differences in the eddy noise are clear: the deterministic model has the lowest transition probabilities; the Gaussian stochastic model has higher transition probabilities; the full model has the highest transition probabilities.

7 Conclusions

This paper formulates a stochastic two-box ocean model modeled after Stommel’s (1961); the model consists of a system of 5 SDEs (1). Previous stochastic Stommel models (e.g. Cessi, 1994; Vélez-Belchı et al., 2001; Monahan et al., 2002; Monahan, 2002; Monahan and Culina, 2011), modeled the atmospheric heat and freshwater forcing as Gaussian stochastic processes, and the exchange of heat and salt between the boxes as a nonlinear drift term corresponding to the large-scale overturning thermohaline circulation. The novelty of the formulation here is that a fast, eddy-driven component is added to the the exchange between the boxes. The terms modeling the eddy-driven exchange are quadratic products of approximately Gaussian random variables; products of jointly-Gaussian random fields were recently found to be an accurate model of eddy-driven exchanges in Grooms (2016).

In more complete and complex ocean models, fast eddy effects are frequently modeled deterministically. Stochastic parameterizations have recently been developed that multiply these deterministic eddy parameterizations by Gaussian random fields (Andrejczuk et al., 2016; Juricke et al., 2017), which is a popular approach for atmospheric models based on the work of Buizza et al. (1999) and Sura et al. (2005). Using methods of averaging and homogenization for slow-fast systems (Pavliotis and Stuart, 2008; Freidlin and Wentzell, 2012; Bouchet et al., 2016), two models were derived approximating the evolution of the slow components (the difference in heat and salt content of the two boxes). The first model (5) replaces the fast eddy-driven exchange terms by a fixed ‘deterministic’ drift term, analogous to the standard approach of deterministic parameterization in more complex ocean models. The second model (7) adds an additional multiplicative noise term accounting for fast variations in the eddy-driven flux. A suite of simulations of each of the three models was used to compare their qualitative behavior. Two regimes were considered: one with a single stable equilibrium, and one with two stable equilibria.

The main results are as follows. There is little qualitative difference in the core of the stationary distributions of the full, non-Gaussian model and the Gaussian multiplicative approximation. In the regime with a single equilibrium the deterministic model has too little variability, but the Gaussian model gives an accurate climatological mean and covariance. In the regime with two stable equilibria the climatological distribution of the three models is nearly the same. In the regime with two stable equilibria the amplitude of the eddies is smaller than in the regime with a single equilibrium, which could perhaps account for the fact that the deterministic model is more accurate in the former regime. Observational estimates suggest that up to 30% of the variability of the Atlantic Meridional Overturning Circulation (AMOC) is driven by ocean eddies, with the rest driven by atmospheric noise Hirschi et al. (2013); Sonnewald et al. (2013).

Though the Gaussian stochastic model gives a good approximation of the core of the climatological distribution, the rare event probabilities are inaccurate. In the single-equilibrium regime there is no clear trend in the behavior: the Gaussian model overestimates rare event probabilities on one side of the mean, and underestimates on the other side. This inaccuracy manifests for short time, transient behavior too: even with a short lead time, the Gaussian model gives inaccurate predictions of the probability of a rare event.

In the regime with two stable equilibria the rare events of interest are the transitions between the two. Despite the reduced amplitude of the eddies in this regime, clear differences were observed in the rates of transition from the neighborhood of one equilibrium to another: the deterministic model had the rarest transitions, and the Gaussian model still made transitions less frequently than the full model.

The goal of the investigation was to investigate the qualitative impacts of non-Gaussian eddy noise of the type observed by Grooms (2016) in a simple model, and to compare to models with Gaussian noise and without eddy noise. The results suggest that Gaussian stochastic parameterizations may be able to successfully produce the day-to-day variability associated with the core of the climatological distribution, but that more accurate non-Gaussian models may be needed to correctly model rare events. Such rare events include extreme behavior like droughts and heat waves, as well as abrupt transitions between climate regimes. The impact of stochastic parameterizations on rare event distributions in climate models has only recently begun to be investigated Tagle et al. (2016).

Acknowledgments

IG is grateful to F. Bouchet for a discussion of the methods from Bouchet et al. (2016) that are applied here. WB is supported by NSF EXTREEMS grant 1407340.