The Modular Arbitrary-Order Ocean-Atmosphere Model: maooam v1.0

Abstract

This paper describes a reduced-order quasi-geostrophic coupled ocean–atmosphere model that allows for an arbitrary number of atmospheric and oceanic modes to be retained in the spectral decomposition. The modularity of this new model allows one to easily modify the model physics. Using this new model, coined the “Modular Arbitrary-Order Ocean-Atmosphere Model” (maooam), we analyse the dependence of the model dynamics on the truncation level of the spectral expansion, and unveil spurious behaviour that may exist at low resolution by a comparison with the higher-resolution configurations. In particular, we assess the robustness of the coupled low-frequency variability when the number of modes is increased. An “optimal” configuration is proposed for which the ocean resolution is sufficiently high, while the total number of modes is small enough to allow for a tractable and extensive analysis of the dynamics.

Introduction

The atmosphere at mid-latitudes displays a variability on a wide range of space scales and timescales, and in particular a low-frequency variability at interannual and decadal timescales as suggested by the analyses of different time series developed in the past years . In contrast to the phenomenon of El Niño–Southern Oscillation (ENSO), of which the driving mechanisms are intensively studied and quite well understood e.g., the origin of mid-latitude low-frequency variability (LFV) remains highly debated, mainly due to the poor ability of state-of-the-art coupled ocean–atmosphere models to simulate it correctly e.g.. The most plausible candidates of this LFV are either the coupling with the local ocean , or teleconnections with the tropical Pacific ocean–atmosphere variability , or both.

Recently the impact of the coupling between the ocean and the atmosphere at mid-latitudes on the atmospheric predictability and the development of the LFV has been explored in a series of low-order coupled ocean–atmosphere systems. However, the limited flexibility of the possible geometries of these previous models led the present authors to develop a series of new model versions. The first of these, OA-QG-WS v1 , for Ocean-Atmosphere–Quasi-Geostrophic–Wind Stress, features only mechanical coupling between the ocean and the atmosphere, and uses 12 atmospheric variables following and four oceanic modes following . In a successor of this model, OA-QG-WS v2, the set of atmospheric variables is extended from 12 to 20 as in . This increase in resolution in the atmosphere was shown to be key to the development of a realistic double gyre in the ocean . A third version of this model, hereafter referred to as vddg in reference to the authors of the model, includes passively advected temperature in the ocean and an energy balance scheme, combined with an extended set of modes for the ocean .

In the vddg model, an LFV associated with the coupling between the ocean and the atmosphere is successfully identified, allowing for extended-range coupled ocean–atmosphere predictions. Moreover, the development of this coupled ocean–atmosphere mode is robust when stochastic forcings are added , or when a seasonal radiative forcing is incorporated into the low-order model . Remarkably the presence of the seasonal radiative input favours the development of the coupled mode due to the amplification of the impact of the wind stress forcing in summer, associated with a drastic reduction of the mixed layer thickness at that period of the year. While these are encouraging results, which suggest the generic character of the coupled ocean–atmosphere mode, they need to be confirmed through the analysis of more sophisticated models, and in particular in higher-resolution coupled systems.

In this article, we present a model that generalizes the vddg model by allowing for an arbitrary number of modes, or basis functions in which the dynamical fields are expanded. The modes can be selected independently for the ocean and the atmosphere, and for the zonal and meridional directions. The modular approach allows one to straightforwardly modify the model physics, such as changing the drag coefficient, introducing new dissipative schemes or adding a seasonal insolation. This model was coined maooam: the Modular Arbitrary-Order Ocean-Atmosphere Model. The model equations and its technical implementation are detailed in Sect. 2. In Sect. 3, maooam is used to investigate the dependence of the model dynamics, i.e. its climatology and the qualitative structure of its attractor, on the number of modes included. Furthermore, the development of the LFV as a function of the spectral truncation is discussed. Key results are summarized in Sect. 4.

Model formulation

The model is composed of a two-layer quasi-geostrophic (QG) atmosphere, coupled both thermally and mechanically to a QG shallow-water ocean layer, in the $\mathit{\beta }$-plane approximation. The atmospheric component is an extension of the QG model, first developed by and further refined by . The equations of motion for the atmospheric streamfunction fields ${\mathit{\psi }}_{\text{a}}^{\mathrm{1}}$ at 250 hPa and ${\mathit{\psi }}_{\text{a}}^{\mathrm{3}}$ at 750 hPa, and the vertical velocity $\mathit{\omega }=\text{d}p/\text{d}t$, read $$\begin{array}{cc}& & {\displaystyle \frac{\partial }{\partial t}}\left({\mathrm{\nabla }}^{\mathrm{2}},{\mathit{\psi }}_{\text{a}}^{\mathrm{1}}\right)+J({\mathit{\psi }}_{\text{a}}^{\mathrm{1}},{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\psi }}_{\text{a}}^{\mathrm{1}})+\mathit{\beta }{\displaystyle \frac{\partial {\mathit{\psi }}_{\text{a}}^{\mathrm{1}}}{\partial x}}& =-k_d^{\prime }{\mathrm{\nabla }}^{\mathrm{2}}({\mathit{\psi }}_{\text{a}}^{\mathrm{1}}-{\mathit{\psi }}_{\text{a}}^{\mathrm{3}})+{\displaystyle \frac{f_{\mathrm{0}}}{\mathrm{\Delta }p}}\mathit{\omega },\\ & & {\displaystyle \frac{\partial }{\partial t}}\left({\mathrm{\nabla }}^{\mathrm{2}},{\mathit{\psi }}_{\text{a}}^{\mathrm{3}}\right)+J({\mathit{\psi }}_{\text{a}}^{\mathrm{3}},{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\psi }}_{\text{a}}^{\mathrm{3}})+\mathit{\beta }{\displaystyle \frac{\partial {\mathit{\psi }}_{\text{a}}^{\mathrm{3}}}{\partial x}}& =+k_d^{\prime }{\mathrm{\nabla }}^{\mathrm{2}}({\mathit{\psi }}_{\text{a}}^{\mathrm{1}}-{\mathit{\psi }}_{\text{a}}^{\mathrm{3}})-{\displaystyle \frac{f_{\mathrm{0}}}{\mathrm{\Delta }p}}\mathit{\omega }-k_d{\mathrm{\nabla }}^{\mathrm{2}}({\mathit{\psi }}_{\text{a}}^{\mathrm{3}}-{\mathit{\psi }}_{\text{o}}).\end{array}$$

The Coriolis parameter $f$ is linearized around a value $f_{\mathrm{0}}$ estimated at latitude ${\mathit{\varphi }}_{\mathrm{0}}=45$${}^{\circ }$ N, $f=f_{\mathrm{0}}+\mathit{\beta }y$, with $\mathit{\beta }=\text{d}f/\text{d}y$. The parameters $k_d^{\prime }$ and $k_d$ quantify the friction between the two atmospheric layers and between the ocean and the atmosphere, respectively, and $\mathrm{\Delta }p=500$ hPa is the pressure difference between the atmospheric layers.

The equation of motion for the streamfunction ${\mathit{\psi }}_{\text{o}}$ of the ocean layer reads cf. $$\begin{array}{cc}& & {\displaystyle \frac{\partial }{\partial t}}\left({\mathrm{\nabla }}^{\mathrm{2}},{\mathit{\psi }}_{\text{o}},-,{\displaystyle \frac{{\mathit{\psi }}_{\text{o}}}{L_{\text{R}}^{\mathrm{2}}}}\right)+J({\mathit{\psi }}_{\text{o}},{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\psi }}_{\text{o}})+\mathit{\beta }{\displaystyle \frac{\partial {\mathit{\psi }}_{\text{o}}}{\partial x}}& =-r{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\psi }}_{\text{o}}+{\displaystyle \frac{C}{\mathit{\rho }h}}{\mathrm{\nabla }}^{\mathrm{2}}({\mathit{\psi }}_{\text{a}}^{\mathrm{3}}-{\mathit{\psi }}_{\text{o}}).\end{array}$$ $L_{\text{R}}$ is the reduced Rossby deformation radius, $\mathit{\rho }$ the density, $h$ the depth, and $r$ the friction at the bottom of the active ocean layer. The rightmost term represents the impact of the wind stress, and is modulated by the drag coefficient of the mechanical ocean–atmosphere coupling, $d=C/\left(\mathit{\rho }h\right)$.

The time evolution of the atmosphere and ocean temperatures $T_{\text{a}}$ and $T_{\text{o}}$ obeys the following equations: $$\begin{array}{cc}& & {\mathit{\gamma }}_{\text{a}}\left({\displaystyle \frac{\partial T_{\text{a}}}{\partial t}},+,J,(,{\mathit{\psi }}_{\text{a}},,,T_{\text{a}},),-,\mathit{\sigma },\mathit{\omega },{\displaystyle \frac{p}{R}}\right)& =-\mathit{\lambda }(T_{\text{a}}-T_{\text{o}})+{\mathit{ϵ}}_{\text{a}}{\mathit{\sigma }}_{\text{B}}{T}_{\text{o}}^{\mathrm{4}}-\mathrm{2}{\mathit{ϵ}}_{\text{a}}{\mathit{\sigma }}_{\text{B}}{T}_{\text{a}}^{\mathrm{4}}+R_{\text{a}},\\ & & {\mathit{\gamma }}_{\text{o}}\left({\displaystyle \frac{\partial T_{\text{o}}}{\partial t}},+,J,(,{\mathit{\psi }}_{\text{o}},,,T_{\text{o}},)\right)& =-\mathit{\lambda }(T_{\text{o}}-T_{\text{a}})-{\mathit{\sigma }}_{\text{B}}{T}_{\text{o}}^{\mathrm{4}}+{\mathit{ϵ}}_{\text{a}}{\mathit{\sigma }}_{\text{B}}{T}_{\text{a}}^{\mathrm{4}}+R_{\text{o}}.\end{array}$$ Here, ${\mathit{\gamma }}_{\text{a}}$ and ${\mathit{\gamma }}_{\text{o}}$ are the heat capacities of the atmosphere and the active ocean layer. ${\mathit{\psi }}_{\text{a}}=({\mathit{\psi }}_{\text{a}}^{\mathrm{1}}+{\mathit{\psi }}_{\text{a}}^{\mathrm{3}})/\mathrm{2}$ is the atmospheric barotropic streamfunction. $\mathit{\lambda }$ is the heat transfer coefficient at the ocean–atmosphere interface, and $\mathit{\sigma }$ is the static stability of the atmosphere, taken to be constant. The quartic terms represent the long-wave radiation fluxes between the ocean, the atmosphere, and outer space, with ${\mathit{ϵ}}_{\text{a}}$ the emissivity of the grey-body atmosphere and ${\mathit{\sigma }}_{\text{B}}$ the Stefan–Boltzmann constant. By decomposing the temperatures as $T_{\text{a}}=T_{\text{a}}^{\mathrm{0}}+\mathit{\delta }{T}_{\text{a}}$ and $T_{\text{o}}=T_{\text{o}}^{\mathrm{0}}+\mathit{\delta }{T}_{\text{o}}$, the quartic terms are linearized around spatially uniform temperatures $T_{\text{a}}^{\mathrm{0}}$ and $T_{\text{o}}^{\mathrm{0}}$, as detailed in Appendix B of . $R_{\text{a}}$ and $R_{\text{o}}$ are the short-wave radiation fluxes entering the atmosphere and the ocean that are also decomposed as $R_{\text{a}}=R_{\text{a}}^{\mathrm{0}}+\mathit{\delta }{R}_{\text{a}}$ and $R_{\text{o}}=R_{\text{o}}^{\mathrm{0}}+\mathit{\delta }{R}_{\text{o}}$.

The hydrostatic relation in pressure coordinates $(\partial \mathrm{\Phi }/\partial p)=-\mathrm{1}/{\mathit{\rho }}_{\text{a}}$ where the geopotential height ${\mathrm{\Phi }}^i=f_{\mathrm{0}}{\mathit{\psi }}_{\text{a}}^i$ and the ideal gas relation $p={\mathit{\rho }}_{\text{a}}R{T}_{\text{a}}$ allow one to write the spatially dependent atmospheric temperature anomaly $\mathit{\delta }{T}_{\text{a}}=\mathrm{2}{f}_{\mathrm{0}}{\mathit{\theta }}_{\text{a}}/R$, with ${\mathit{\theta }}_{\text{a}}\equiv ({\mathit{\psi }}_{\text{a}}^{\mathrm{1}}-{\mathit{\psi }}_{\text{a}}^{\mathrm{3}})/\mathrm{2}$ often referred to as the baroclinic streamfunction. $R$ is the ideal gas constant. This can be used to eliminate the vertical velocity $\mathit{\omega }$ from Eqs. ()–() and (). This reduces the independent dynamical fields to the streamfunction fields ${\mathit{\psi }}_{\text{a}}$ and ${\mathit{\psi }}_{\text{o}}$, and the spatially dependent temperatures $\mathit{\delta }{T}_{\text{a}}$ and $\mathit{\delta }{T}_{\text{o}}$.

The prognostic equations for these four fields are then non-dimensionalized by dividing time by $f_{\mathrm{0}}^{-\mathrm{1}}$, distance by a characteristic length scale $L$, pressure by the difference $\mathrm{\Delta }p$, temperature by $f_{\mathrm{0}}^{\mathrm{2}}{L}^{\mathrm{2}}/R$, and streamfunction by $L^{\mathrm{2}}{f}_{\mathrm{0}}$. A more detailed discussion of the model equations and their non-dimensionalization can be found in and .

All the parameters of the model equations used in the present work are listed in Table .

Values of the parameters of the model that are used in the analyses of Sect. .

Parameter (unit)	Value	Parameter (unit)	Value
$n=\mathrm{2}{L}_y/L_x$	$1.5$	$L_{\text{R}}$ (km)	$19.93$
$L_y=\mathit{\pi }L$ (km)	$5.0\times {10}^{\mathrm{3}}$	$\mathit{\rho }$ (kg m${}^{-\mathrm{3}}$)	$1000$
$f_{\mathrm{0}}$ (s${}^{-\mathrm{1}}$)	$1.032\times {10}^{-\mathrm{3}}$	${\mathit{\sigma }}_{\text{B}}$ (W m${}^{\mathrm{2}}$ K${}^{-\mathrm{4}}$)	$5.6\times {5.610}^{-\mathrm{8}}$
$\mathit{\lambda }$ (W m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$)	$15.06$	$\mathit{\sigma }$ (m${}^{\mathrm{2}}$ s${}^{-\mathrm{2}}$ Pa${}^{-\mathrm{2}}$)	$2.16\times {10}^{-\mathrm{6}}$
$r$ (s${}^{-\mathrm{1}}$)	$1.0\times {10}^{-\mathrm{7}}$	$\mathit{\beta }$ (m${}^{-\mathrm{1}}$ s${}^{-\mathrm{1}}$)	$1.62\times {10}^{-11}$
$d$ (s${}^{-\mathrm{1}}$)	$1.1\times {10}^{-\mathrm{7}}$	$R$ (J kg${}^{-\mathrm{1}}$ K${}^{-\mathrm{1}}$)	$287$
$C_{\text{o}}$ (W m${}^{-\mathrm{2}}$)	$310$	${\mathit{\gamma }}_{\text{o}}$ (J m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$)	$5.46\times {10}^{\mathrm{8}}$
$C_{\text{a}}$ (W m${}^{-\mathrm{2}}$)	$C_{\text{o}}/\mathrm{3}$	${\mathit{\gamma }}_{\text{a}}$ (J m${}^{-\mathrm{2}}$ K${}^{-\mathrm{1}}$)	$1.0\times {10}^{\mathrm{7}}$
$k_d$ (s${}^{-\mathrm{1}}$)	$3.0\times {10}^{-\mathrm{6}}$	$T_{\text{a}}^{\mathrm{0}}$ (K)	$289.30$
$k_d^{\prime }$ (s${}^{-\mathrm{1}}$)	$3.0\times {10}^{-\mathrm{6}}$	$T_{\text{o}}^{\mathrm{0}}$ (K)	$301.46$
$h$ (m)	$136.5$	${\mathit{ϵ}}_{\text{a}}$	$0.7$

Expansion of the dynamical fields

In non-dimensionalized coordinates $x^{\prime }=x/L$ and $y^{\prime }=y/L$, the domain is defined by $(\mathrm{0}\le x^{\prime }\le \frac{\mathrm{2}\mathit{\pi }}{n},\mathrm{0}\le y^{\prime }\le \mathit{\pi })$, with $n=\mathrm{2}{L}_y/L_x$ the aspect ratio between its meridional and zonal extents (see Table  for the value used here). The atmospheric flow is defined in a zonally periodic channel with no-flux boundary conditions in the meridional direction ($\partial {\cdot }_{\text{a}}/\partial x^{\prime }\equiv \mathrm{0}$ at $y^{\prime }=\mathrm{0},\mathit{\pi }$), whereas the oceanic flow is confined within an ocean basin by imposing no-flux boundaries in both the meridional ($\partial {\cdot }_{\text{o}}/\partial x^{\prime }\equiv \mathrm{0}$ at $y^{\prime }=\mathrm{0},\mathit{\pi }$) and zonal ($\partial {\cdot }_{\text{o}}/\partial y^{\prime }\equiv \mathrm{0}$ at $x^{\prime }=\mathrm{0},\mathrm{2}\mathit{\pi }/n$) directions. These boundary conditions limit the functions used in the Fourier expansion of the dynamical fields. With the proper normalization, the basis functions for the atmosphere must be of the following form, following the nomenclature of : $$\begin{array}{}& & F_P^A(x^{\prime },y^{\prime })=\sqrt{\mathrm{2}}\cos \left(P{y}^{\prime }\right)& & F_{M,P}^K(x^{\prime },y^{\prime })=\mathrm{2}\cos \left(Mn{x}^{\prime }\right)\sin \left(P{y}^{\prime }\right)& & F_{H,P}^L(x^{\prime },y^{\prime })=\mathrm{2}\sin \left(Hn{x}^{\prime }\right)\sin \left(P{y}^{\prime }\right).\end{array}$$

Analogously, the oceanic basis functions must be of the form $${\mathit{\varphi }}_{H_{\text{o}},P_{\text{o}}}(x^{\prime },y^{\prime })=\mathrm{2}\sin \left({\displaystyle \frac{H_{\text{o}}n}{\mathrm{2}}}{x}^{\prime }\right)\sin \left(P_{\text{o}}{y}^{\prime }\right),$$ with integer values of $M$, $H$, $P$, $H_{\text{o}}$, and $P_{\text{o}}$.

For example, the spectral truncation used by can be specified as Eqs. ()–() with $M=H=\mathrm{1}$; $P\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$. extend this set by two blocks of two functions each, and the resulting set can be specified as $M,H\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$; $P\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$. The vddg model has $M,H\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$; $P\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$ and $H_{\text{o}}\in \mathit{\{}\mathrm{1},\mathrm{2}\mathit{\}}$; $P_{\text{o}}\in \mathit{\{}\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}\mathit{\}}$. Note that, for consistency, the ranges for $M$ and $H$ should be the same. The distinction between $M$ and $H$ is, however, required to avoid ambiguities in the formulae of the inner products, as specified in Appendix .

For the given ranges of $\mathrm{1}\le P\le P^{max}$, $\mathrm{1}\le (M,H)\le H^{max}$ and $\mathrm{1}\le P_{\text{o}}\le P_{\text{o}}^{max}$, and $\mathrm{1}\le H_{\text{o}}\le H_{\text{o}}^{max}$, the number of basis functions can be calculated as $$n_{\text{a}}=P^{max}(\mathrm{2}{H}^{max}+\mathrm{1});n_{\text{o}}=P_{\text{o}}^{max}{H}_{\text{o}}^{max}.$$

Ordering the basis functions as in Eqs. ()–(), along increasing values of $M=H_{\text{(o)}}$ and then $P_{\text{(o)}}$, allows one to write the set as $\left\{F_i,(,x^{\prime },,,y^{\prime },),,,,{\mathit{\varphi }}_j,(,x^{\prime },,,y^{\prime },)\right\}$ $(\mathrm{1}\le i\le n_{\text{a}},\mathrm{1}\le j\le n_{\text{o}})$. The dynamical fields can then be written as the following truncated series expansions: $$\begin{array}{cc}& {\mathit{\psi }}_{\text{a}}(x^{\prime },y^{\prime },t)& =\sum _{i=\mathrm{1}}^{n_{\text{a}}}{\mathit{\psi }}_{\text{a},i}\left(t\right)F_i(x^{\prime },y^{\prime }),& \mathit{\delta }{T}_{\text{a}}(x^{\prime },y^{\prime },t)& =\sum _{i=\mathrm{1}}^{n_{\text{a}}}\mathit{\delta }{T}_{\text{a},i}\left(t\right)F_i(x^{\prime },y^{\prime }),& =\mathrm{2}{\displaystyle \frac{f_{\mathrm{0}}}{R}}\sum _{i=\mathrm{1}}^{n_{\text{a}}}{\mathit{\theta }}_{\text{a},i}\left(t\right)F_i(x^{\prime },y^{\prime }),\\ & {\mathit{\psi }}_{\text{o}}(x^{\prime },y^{\prime },t)& =\sum _{j=\mathrm{1}}^{n_{\text{o}}}{\mathit{\psi }}_{\text{o},j}\left(t\right)\left({\mathit{\varphi }}_j\right(x^{\prime },y^{\prime })-\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_j}),& \mathit{\delta }{T}_{\text{o}}(x^{\prime },y^{\prime },t)& =\sum _{j=\mathrm{1}}^{n_{\text{o}}}\mathit{\delta }{T}_{\text{o},j}\left(t\right){\mathit{\varphi }}_j(x^{\prime },y^{\prime }).\end{array}$$

Furthermore, the short-wave radiation or insolation is determined by $\mathit{\delta }{R}_{\text{a}}=C_{\text{a}}{F}_{\mathrm{1}}$; $\mathit{\delta }{R}_{\text{o}}=C_{\text{o}}{F}_{\mathrm{1}}$. In Eq. (), a term $\stackrel{\mathrm{‾}}{{\mathit{\varphi }}_j}$ is added to the oceanic basis function ${\mathit{\varphi }}_j(x^{\prime },y^{\prime })$ in order to give it a vanishing spatial average. This is required to guarantee mass conservation in the ocean , but otherwise does not affect the dynamics. Indeed, it can be added a posteriori when plotting the field ${\mathit{\psi }}_{\text{o}}(x^{\prime },y^{\prime },t)$. This term is non-zero for odd $P_{\text{o}}$ and $H_{\text{o}}$, $$\begin{array}{cc}& \stackrel{\mathrm{‾}}{{\mathit{\varphi }}_j}& ={\displaystyle \frac{n}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\frac{\mathrm{2}\mathit{\pi }}{n}}{\int }}{\mathit{\varphi }}_j(x^{\prime },y^{\prime })\text{d}{x}^{\prime }\text{d}{y}^{\prime }& =\mathrm{2}{\displaystyle \frac{\left(\right(-\mathrm{1}{)}^{H_{\text{o}}}-\mathrm{1}\left)\right((-\mathrm{1}{)}^{P_{\text{o}}}-\mathrm{1})}{H_{\text{o}}{P}_{\text{o}}{\mathit{\pi }}^{\mathrm{2}}}}.\end{array}$$ The mass conservation is automatically satisfied for ${\mathit{\psi }}_{\text{a}}(x^{\prime },y^{\prime },t)$, as the spatial averages of the atmospheric basis functions $F_i(x^{\prime },y^{\prime })$ are zero.

Substituting the fields in Eqs. ()–() and projecting on the different basis functions yield $\mathrm{2}(n_{\text{a}}+n_{\text{o}})$ ordinary differential equations (ODEs) for as many variables. Due to the linearization of the quartic temperature fields in Eqs. () and (), these equations are at most bilinear (due to the advection term) in the variables ${\mathit{\psi }}_{\text{a},i},{\mathit{\theta }}_{\text{a},i},{\mathit{\psi }}_{\text{o},j}$, and $\mathit{\delta }{T}_{\text{o},j}$, which will henceforth jointly be referred to as ${\mathit{\eta }}_i$, the components of the state vector $\mathit{\eta }$.

To construct the dynamical equations of these variables, one has to compute the various projections or inner products with the basis functions, for which the following shorthand notation will be used: $$⟨S,G⟩\equiv {\displaystyle \frac{n}{\mathrm{2}{\mathit{\pi }}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }/n}{\int }}S(x^{\prime },y^{\prime })G(x^{\prime },y^{\prime })\text{d}{x}^{\prime }\text{d}{y}^{\prime }.$$

As described by , the inner products for the atmosphere can be computed as purely algebraic formulae of the wave numbers $P$, $M$, and $H$. We reiterate these algebraic formulae in Sect.  of Appendix  and extend them with the formulae for both the ocean–atmosphere coupling terms and the ocean inner products in Sect. . The inner products can be represented as either two-dimensional or three-dimensional tensors, which are sparse but generally not diagonal.

Technical implementation

Substituting the fields by Eqs. ()–() and calculating the coefficients using the expressions for the inner products as in Appendix  yields a set of $N\equiv \mathrm{2}(n_{\text{a}}+n_{\text{o}})$ prognostic ordinary differential equations. These equations are at most bilinear in the variables ${\mathit{\eta }}_i$ $(\mathrm{1}\le i\le N$) due to the linearization of the radiative terms around a reference temperature present in Eqs. ()–(). This system of ODEs can therefore be most generically expressed as the sum of a constant, a matrix multiplication, and a tensor contraction: $$\begin{array}{cc}& {\displaystyle \frac{\text{d}{\mathit{\eta }}_i}{\text{d}t}}& =c_i+\sum _{j=\mathrm{1}}^N{m}_{i,j}{\mathit{\eta }}_j& +\sum _{j,k=\mathrm{1}}^N{t}_{i,j,k}{\mathit{\eta }}_j{\mathit{\eta }}_k(\mathrm{1}\le i\le N).\end{array}$$ This expression can be further simplified by adding a dummy variable that is identically equal to one: ${\mathit{\eta }}_{\mathrm{0}}\equiv \mathrm{1}$. This extra variable allows one to merge $c_i$, $m_{i,j}$, and $t_{i,j,k}$ into the tensor ${\mathcal{T}}_{i,j,k}$, in which the linear terms are represented by ${\mathcal{T}}_{i,j,\mathrm{0}}$ and the constant term by ${\mathcal{T}}_{i,\mathrm{0},\mathrm{0}}$: $${\displaystyle \frac{\text{d}{\mathit{\eta }}_i}{\text{d}t}}=\sum _{j=\mathrm{0}}^N\sum _{k=\mathrm{0}}^N{\mathcal{T}}_{i,j,k}{\mathit{\eta }}_j{\mathit{\eta }}_k(\mathrm{1}\le i\le N).$$ The elements of the tensor ${\mathcal{T}}_{i,j,k}$ are specified in Appendix . Recasting the system of ordinary differential equations for ${\mathit{\eta }}_i$ in the form of a tensor contraction has certain advantages, as we will clarify below. The symmetry of Eq. () allows for a unique representation of ${\mathcal{T}}_{i,j,k}$, if it is taken to be upper triangular in the last two indices (${\mathcal{T}}_{i,j,k}\equiv \mathrm{0}$ if $j>k$). Since ${\mathcal{T}}_{i,j,k}$ is known to be sparse, it is stored using the coordinate list representation, i.e. a list of tuples $(i,j,k,{\mathcal{T}}_{i,j,k})$. This representation renders the computation of the tendencies $\text{d}{\mathit{\eta }}_i/\text{d}t$ computationally very efficient as well as conveniently parallelizable.

Two implementations of maooam are provided as a Supplement: one in Lua and one in Fortran. The Lua code is optimized for LuaJIT, a just-in-time compiler for Lua , and runs about 20 % slower than the Fortran version. By default, the model equations are numerically integrated using the Heun method. We have tested higher-accuracy methods, but these did not significantly change the results. The integration method can easily be changed; as an example, a fourth-order Runge–Kutta integrator is also included in the Lua implementation.

Derivation of Jacobian, tangent linear, and adjoint models

The form of Eq. () allows one to easily compute the Jacobian matrix of this system of ODEs. Indeed, denoting the right-hand side of Eq. () as $\text{d}{\mathit{\eta }}_i/\text{d}t=f_i$, the expression reduces to $$\begin{array}{cc}& J_{i,j}={\displaystyle \frac{\text{d}{f}_i}{\text{d}{\mathit{\eta }}_j}}& =\text{d}\left(\sum _{k,l=\mathrm{0}}^N{\mathcal{T}}_{i,k,l}{\mathit{\eta }}_k{\mathit{\eta }}_l\right)/\text{d}{\mathit{\eta }}_j& =\sum _{k=\mathrm{0}}^N\left({\mathcal{T}}_{i,k,j},+,{\mathcal{T}}_{i,j,k}\right){\mathit{\eta }}_k(\mathrm{1}\le i,j\le N).\end{array}$$ The differential form of the tangent linear (TL) model for a small perturbation $\mathit{\delta }{\mathit{\eta }}^{\text{TL}}$ of a trajectory ${\mathit{\eta }}^{\ast }$ is then simply $$\begin{array}{cc}& {\displaystyle \frac{\text{d}\mathit{\delta }{\mathit{\eta }}_i^{\text{TL}}}{\text{d}t}}& =\sum _{j=\mathrm{1}}^N{J}_{i,j}^{\ast }\mathit{\delta }{\mathit{\eta }}_j^{\text{TL}}& =\sum _{j=\mathrm{1}}^N\sum _{k=\mathrm{0}}^N\left({\mathcal{T}}_{i,k,j},+,{\mathcal{T}}_{i,j,k}\right){\mathit{\eta }}_k^{\ast }\mathit{\delta }{\mathit{\eta }}_j^{\text{TL}}(\mathrm{1}\le i\le N).\end{array}$$ To obtain the differential form of the adjoint model along the trajectory ${\mathit{\eta }}^{\ast }$, the Jacobian is transposed to yield the following equations for the adjoint variable $\mathit{\delta }{\mathit{\eta }}^{\text{AD}}$: $$\begin{array}{cc}& -{\displaystyle \frac{\text{d}\mathit{\delta }{\mathit{\eta }}_i^{\text{AD}}}{\text{d}t}}& =\sum _{j=\mathrm{1}}^N{J}_{j,i}^{\ast }\mathit{\delta }{\mathit{\eta }}_j^{\text{AD}}& =\sum _{j=\mathrm{1}}^N\sum _{k=\mathrm{0}}^N\left({\mathcal{T}}_{j,k,i},+,{\mathcal{T}}_{j,i,k}\right){\mathit{\eta }}_k^{\ast }\mathit{\delta }{\mathit{\eta }}_j^{\text{AD}}(\mathrm{1}\le i\le N).\end{array}$$

Model dynamics

This section details some key results obtained with the model for various levels of spectral truncation, with the set of parameter values given in Table . The parameter values for $L$, $L_{\text{R}}$, $\mathit{\lambda }$, $r$, $d$, $C_{\text{o}}$, $C_{\text{a}}$, $k_d$, and $k_d^{\prime }$ were selected as detailed in . The same value was chosen for $k_d$ and $k_d^{\prime }$, as was done in ; see also . Unless otherwise stated, all the following results are obtained after first integrating the model for a transient period of 30 726.5 years. The model is subsequently integrated for another 92 179.6 years to obtain a sufficiently long trajectory from which good statistics can be extracted.

For the atmospheric part of the model, a previous study , referred to as CT in the following, has shown that spurious chaos and a too large variability in the modes near the spectral cut-off could take place if the resolution is not high enough. These manifestations of spurious behaviour can lead to solutions that differ significantly from the solutions of the full partial differential equations (PDEs, here Eqs. –). These findings lead us to the important question of convergence: to what degree has the solution of the truncated equations converged towards the solution of the PDEs? Although we do not have access to the latter, one can infer how the solutions are altered when the resolution is increased. Therefore, it cannot be asserted that convergence has been reached, and this point was also clearly stated in CT. However, we can reasonably suppose that when the solutions stabilize, they give an insight into the full dynamics.

Number of variables, transient time, and effective runtime of the runs (in years).

Resolution	No. of	Transient	Eff. runtime
	variables	($y$)	($y$)
atm. $\mathrm{2}x$–$\mathrm{2}y$ oc. $\mathrm{2}x$–$\mathrm{4}x$	$36$	30 726.5	92 179.6
atm. $\mathrm{2}x$–$\mathrm{4}y$ oc. $\mathrm{2}x$–$\mathrm{4}x$	$56$	30 726.5	92 179.6
atm. $\mathrm{3}x$–$\mathrm{3}y$ oc. $\mathrm{3}x$–$\mathrm{3}x$	$60$	30 726.5	92 179.6
atm. $\mathrm{4}x$–$\mathrm{4}y$ oc. $\mathrm{4}x$–$\mathrm{4}x$	$104$	30 726.5	92 179.6
atm. $\mathrm{5}x$–$\mathrm{5}y$ oc. $\mathrm{5}x$–$\mathrm{5}x$	$160$	30 726.5	92 179.6
atm. $\mathrm{6}x$–$\mathrm{6}y$ oc. $\mathrm{6}x$–$\mathrm{6}x$	$228$	30 726.5	92 179.6
atm. $\mathrm{7}x$–$\mathrm{7}y$ oc. $\mathrm{7}x$–$\mathrm{7}x$	$308$	30 726.5	92 179.6
atm. $\mathrm{8}x$–$\mathrm{8}y$ oc. $\mathrm{8}x$–$\mathrm{8}x$	$400$	30 726.5	92 179.6
atm. $\mathrm{9}x$–$\mathrm{9}y$ oc. $\mathrm{9}x$–$\mathrm{9}y$	$504$	30 726.5	92 179.6
atm. $10x$–$10$ oc. $10x$–$10y$	$620$	30 726.5	92 179.6
atm. $\mathrm{5}x$–$\mathrm{5}y$ oc. $12x$–$12y$	$398$	15 363.3	92 179.6
atm. $12x$–$12y$ oc. $12x$–$12y$	$888$	15 363.3	74 972.7

Cross section of the attractors for various model resolutions. The atmospheric and oceanic resolutions are both indicated above each panel. The parameters are given in Table .

[Figure omitted. See PDF]

Cross section of the attractors for various model resolutions (continued from Fig. ).

[Figure omitted. See PDF]

This question is now addressed for the maooam coupled atmosphere–ocean model. Figures  and display cross sections of the attractors of the model for different resolutions. The three variables selected in this projection are ${\mathit{\psi }}_{\text{a},\mathrm{1}}$, ${\mathit{\psi }}_{\text{o},\mathrm{2}}$, and ${\mathit{\theta }}_{\text{o},\mathrm{2}}$, which have already been used to represent the large-scale variability of the model . We use the same notation as in CT to specify the resolution of each component: $\left(H^{max}\right)x$–$\left(P^{max}\right)y$ for the atmosphere and $\left(H_{\text{o}}^{max}\right)x$–$\left(P_{\text{o}}^{max}\right)y$ for the ocean, with $M^{max}=H^{max}$. All the model configurations used are listed in Table . To alleviate the notation in the following, a model configuration denoted simply by $H^{max}x$–$P^{max}y$ indicates that the resolution is the same in both components: $H_{\text{o}}^{max}=H^{max}$ and $P_{\text{o}}^{max}=P^{max}$.

Variance distributions of the ${\mathit{\psi }}_{\text{a},i}$ variables in percents for various model resolutions. For the variables associated with the $A$-type basis functions, the wave numbers $M$ and $H$ are not defined (nd).

[Figure omitted. See PDF]

The first panel of Fig. , with the atm. $\mathrm{2}x$–$\mathrm{2}y$ oc. $\mathrm{2}x$–$\mathrm{4}y$ resolution, shows the typical attractor geometry found in  and  with a noisy, seemingly periodic orbit associated with the development of a large low-frequency signal. However, as the resolution is increased in both the ocean and atmosphere components, this structure destabilizes and we obtain more compact, noisy attractors in Figs.  and . The cause of this structural change is an interesting question in itself, which is worth exploring further in the future, as it is associated with the problem of structural stability of models, but is beyond the scope of the present work.

Regarding the question of convergence, the variability of the atmospheric variables becomes quite stable as the resolution increases beyond $\mathrm{6}x$–$\mathrm{6}y$. Indeed, the bounds of the attractors on the vertical axis (${\mathit{\psi }}_{\text{a},\mathrm{1}}$) stabilize at this resolution. This result is in agreement with the findings of CT. On the other hand, the convergence is not yet reached for the oceanic variables whose variability is strongly affected by adding further modes as in the $\mathrm{7}x$–$\mathrm{7}y$ and $\mathrm{8}x$–$\mathrm{8}y$ resolutions.

The impact of the resolution on the solutions can also be examined by computing the variance of each variable of the barotropic and baroclinic streamfunctions, since these are associated with the kinetic and potential energy of the system . The presence of spurious behaviour can then be detected through substantial changes in this variability. The distributions of the total variance of the variables ${\mathit{\psi }}_{\text{a},i}$ and ${\mathit{\psi }}_{\text{o},i}$ are depicted in Figs. –. The results show that the variance distribution does not change much beyond the $\mathrm{4}x$–$\mathrm{4}y$ resolution for the atmospheric component. However, for the oceanic component, the variance distribution is strongly modified when the resolution increases, and therefore one cannot conclude from Fig.  that some sort of convergence is reached at the $\mathrm{8}x$–$\mathrm{8}y$ resolution. To interpret this specific property, one must recall an important feature of two-dimensional quasi-geostrophic turbulence, namely the presence of a specific space scale, the Rhines scale, which delimits the two regimes associated with a wave-dominated dynamics and a turbulent dynamics. This space scale is given by

Variance distributions of the ${\mathit{\psi }}_{\text{a},i}$ variables in percents for various model resolutions (continued from Fig. ). For the variables associated with the $A$-type basis functions, the wave numbers $M$ and $H$ are not defined (nd).

[Figure omitted. See PDF]

$$L_{\mathit{\text{Rh}}}=\sqrt{{\displaystyle \frac{U}{\mathit{\beta }}}},$$ where $U$ represents the root-mean-square velocity of the energy-containing scales  and $\mathit{\beta }=\text{d}f/\text{d}y$ is the meridional derivative of the Coriolis parameter $f$. If one takes the typical velocity of the order of a few metres per second and a few centimetres per second within the atmosphere and the ocean at large scales, the typical length scales will be of the order of 1000 and 100 km, respectively. Therefore the highest wave numbers necessary to resolve the wave-dominated part within the atmosphere and the ocean differ by a factor of 10. Coming back to our analysis, if this limit is reached for the atmosphere in our model at $H$, $P=\mathrm{4}$–5, we should suspect that a value of $H_{\text{o}}/\mathrm{2}\approx P_{\text{o}}\approx 40$–50 should be used for the ocean. This of course imposes strong constraints on our reduced-order model and would considerably limit its flexibility.

Variance distributions of the ${\mathit{\psi }}_{\text{o},i}$ variables in percents for various model resolutions.

[Figure omitted. See PDF]

Let us now focus on the development of the LFV in these different model configurations, and let us define the geopotential height difference $\mathit{\delta }z$ between the locations ($\mathit{\pi }/n,\mathit{\pi }/\mathrm{4}$) and ($\mathit{\pi }/n,\mathrm{3}\mathit{\pi }/\mathrm{4}$) of the model's non-dimensional domain: $$\begin{array}{cc}& \mathit{\delta }z\left(t\right)=z(\mathit{\pi }/n,\mathit{\pi }/\mathrm{4},t)-z(\mathit{\pi }/n,\mathrm{3}\mathit{\pi }/\mathrm{4},t),\\ & z(x^{\prime },y^{\prime },t)={\displaystyle \frac{f_{\mathrm{0}}}{g}}{\mathit{\psi }}_{\text{a}}(x^{\prime },y^{\prime },t),\end{array}$$ where $z$ is the geopotential height field, as in . The results shown in Figs.  and indicate that the LFV, present for atm. $\mathrm{2}x$–$\mathrm{2}y$ oc. $\mathrm{2}x$–$\mathrm{4}y$ as in , is a very weak signal at intermediate resolutions, but develops again when the number of modes is increased, as shown by the 1-year and 5-year running means. It suggests that the LFV previously found in low-resolution versions (see Fig. , panel atm. $\mathrm{2}x$–$\mathrm{2}y$ oc. $\mathrm{2}x$–$\mathrm{4}y$) is a robust feature of the model. Moreover, at high resolutions this LFV is weaker than for the VDDG model version, but it seems closer to the actual dynamics found for the North Atlantic Oscillation (NAO) as discussed in and .

Variance distributions of the ${\mathit{\psi }}_{\text{o},i}$ variables in percents for various model resolutions (continued from Fig. ).

[Figure omitted. See PDF]

Time series of the geopotential height difference (m) between locations ($\mathit{\pi }/n,\mathit{\pi }/\mathrm{4})$ and ($\mathit{\pi }/n,\mathrm{3}\mathit{\pi }/\mathrm{4})$ of the model's non-dimensional domain for different resolutions. Running averages for $\mathit{\tau }=\mathrm{1}y$ (black) and $\mathit{\tau }=\mathrm{5}y$ (red) are also provided, highlighting the LFV signal present in the series.

[Figure omitted. See PDF]

Time series of the geopotential height difference (continued from Fig. ).

[Figure omitted. See PDF]

The climatologies of the atmospheric barotropic streamfunction expressed in geopotential height further highlight the changes in the statistical properties of the model as a function of resolution. As shown in Figs.  and , the convergence is pretty fast toward an averaged zonal atmospheric circulation as the model resolution is increased. By contrast, the convergence for the oceanic streamfunction ${\mathit{\psi }}_{\text{o}}$ is less clear (Figs.  and ), although a recurrent “global” double gyre is present for each resolution. As for the LFV, the topology of the gyres at high resolutions and their small-scale structures also seem to depend on whether $H^{max}$, $M^{max}$, $H_{\text{o}}^{max}$ and $P_{\text{o}}^{max}$, $P^{max}$ are even or odd numbers.

The previous results point toward the important question of the optimal resolution of the oceanic component needed to get a sufficiently low-resolution model while keeping a dynamics with strong similarities to a very high-resolution model. To answer this question, we have performed some higher-resolution integrations, but on shorter time spans. The time span for each integration is given in Table .

The variance distributions of the oceanic streamfunction variables (see Fig. ) have decreased at the spectral cut-off's edges compared to the distributions of the lower-resolution model configurations shown in Fig. . However, this decrease is not sufficient, and apparently spurious effects are still present. For instance, the decay is not identical in both directions, with a slower decay rate as the zonal wave number $H_{\text{o}}$ increases. We can even notice a peak in the distribution around $H_{\mathrm{0}}=H_{\mathrm{0}}^{max}$ and $P_{\mathrm{0}}=\mathrm{2}$ for all these higher model resolutions. This indicates that in fact we are still far from a quantitatively representative solution in the ocean. It confirms that, as stated previously, a resolution of the order of the Rhines scale is needed to achieve a good convergence. For the ocean, it corresponds to a 100 km resolution which would then require roughly 2000 modes. Such a model will of course be very computationally expensive and cannot be considered a “reduced”-order model anymore.

Climatologies for the geopotential height field $z=\frac{f_{\mathrm{0}}}{g}{\mathit{\psi }}_{\text{a}}$ (m) presented on the non-dimensional model domain, as obtained using 92 179.6 years of model integration.

[Figure omitted. See PDF]

Climatologies for the geopotential height field $z=\frac{f_{\mathrm{0}}}{g}{\mathit{\psi }}_{\text{a}}$ (m) (continued from Fig. ).

[Figure omitted. See PDF]

However, the comparison between the atm. $\mathrm{5}x$–$\mathrm{5}y$ oc. $12x$–$12y$ model configuration and the $10x$–$10y$ or $12x$–$12y$ model configurations shows that the former displays a large-scale behaviour close to the latter two, but with a reduced complexity and computational cost. This similarity can be assessed by considering the climatologies of these higher-resolution runs displayed in Fig.  and by watching the corresponding videos (see below). We therefore believe that the atm. $\mathrm{5}x$–$\mathrm{5}y$ oc. $12x$–$12y$ model configuration is a good candidate when investigating more realistic dynamics than the one presented in vddg. It must however be stressed that the vddg model is still an important tool in this hierarchy of models since it already contains the basic mechanisms leading to low-frequency variability. In addition, the climatologies shown in Fig.  confirm the dependence of the dynamics on whether $H^{max}$, $M^{max}$, $H_{\text{o}}^{max}$ and $P_{\text{o}}^{max}$, $P^{max}$ are even or odd, and also the presence of a global double-gyre in the ocean.

Climatologies for the oceanic streamfunction field ${\mathit{\psi }}_{\text{o}}$ (m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$) presented on the non-dimensional model domain, as obtained using 92 179.6 years of model integration.

[Figure omitted. See PDF]

Climatologies for the oceanic streamfunction field ${\mathit{\psi }}_{\text{o}}$ (m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$) (continued from Fig. ).

[Figure omitted. See PDF]

Finally, the dynamics of the model for the various resolutions are also illustrated in the videos provided as supplementary material. These videos depict the time evolution of the streamfunction and temperature fields, as well as the geopotential height difference and the three-dimensional phase-space projection shown in Figs.  and . They give an insight into the high-frequency atmospheric and low-frequency oceanic variability, and also show the interesting time evolution of the oceanic gyres. In these videos, a striking feature is the presence of a westward wave propagation within the ocean while the LFV is developing in the coupled system. This feature has been associated with the propagation of Rossby-like waves .

Conclusions

A new reduced-order coupled ocean–atmosphere model is presented, extending the low-resolution versions previously published . It is referred to as maooam, for Modular Arbitrary-Order Ocean-Atmosphere Model. This new model retains the main features of the previous versions but allows for the selection of an arbitrary resolution within the ocean and the atmosphere. Besides the potential utility of this new functionality for evaluating the impact of the number of modes on the dynamics (as has been done here), it opens the possibility of addressing several new questions in a very flexible way, such as the development of a consistent stochastic parameterization scheme through scales, the understanding of the predictability problem at multiple scales and the role of model error, or the implementation of a data-assimilation scheme for the coupled ocean–atmosphere system.

Variance distributions of the ${\mathit{\psi }}_{\text{o},i}$ variables in percents for the high-resolution runs. The information on the different runtimes is gathered in Table .

[Figure omitted. See PDF]

In the present work, we have studied the impact of the resolution on the model solution's dynamics, by investigating the properties of the attractors and the variance distributions in both the oceanic and atmospheric components. The conclusion that can be drawn is that the convergence of the atmospheric component of the system is quite fast (as noted in ), with variance distributions decreasing rapidly as a function of scale. However, the convergence of the oceanic component is much slower. Consequently, none of the solutions presented so far have satisfactorily converged toward a dynamics that correctly reflects the wave-dominated regime of the coupled ocean–atmosphere system. This regime corresponds to a resolution associated with the Rhines scale (which for the ocean is equal to 100 km or, equivalently, to wave numbers of the order of $H_{\text{o}}^{max}/\mathrm{2}\approx P_{\text{o}}^{max}\approx 50$). This stresses the need for high-resolution oceanic models to correctly represent the full coupled dynamics. One coupled model configuration which could, however, be recommended so far is the atm. $\mathrm{5}x$–$\mathrm{5}y$ oc. $12x$–$12y$ configuration, which seems to display some robustness in the ocean climatology as compared to the full $10x$–$10y$ and $12x$–$12y$ configurations. This conclusion requires further investigation with even higher resolutions, together with the use of more advanced tools of analysis like the computation of the Lyapunov exponents as in . These can be computed using the tangent linear model version for which an implementation is also provided. This will be the subject of a future investigation.

Climatologies for the oceanic streamfunction ${\mathit{\psi }}_{\text{o}}$ field (m${}^{\mathrm{2}}$ s${}^{-\mathrm{1}}$) presented on the non-dimensional model domain for the high-resolution runs displayed in Fig. .

[Figure omitted. See PDF]

The robustness of the LFV pattern, one of the most interesting features of the model, has also been explored. As it turns out, a LFV is still present in a large portion of the model configurations explored (not in $\mathrm{2}x$–$\mathrm{4}y$, $\mathrm{3}x$–$\mathrm{3}y$, and $\mathrm{4}x$–$\mathrm{4}y$), but a weaker LFV signal is found when high-resolution configurations are used. A dominant signal is found with a wide variety of periods ranging from 1 to 100 years, depending on the model configuration. A more detailed analysis of the underlying structure of the system's attractor is needed to clarify the origin of this diversity, for instance through a bifurcation analysis as in . Note that the VDDG model is still an important tool in this hierarchy of models, since it already contains the basic mechanisms leading to the LFV.

Another interesting finding is the change of structure of the climatologies of the ocean gyres when choosing even or odd wave numbers ($H^{max}$, $M^{max}$, $H_{\text{o}}^{max}$ and $P_{\text{o}}^{max}$, $P^{max}$). Is this feature purely associated with the convergence toward a spatially continuous field, or does it reflect specific properties of the dynamical equations, such as symmetries or invariance? These questions are still open and will be the subject of a future investigation that should allow one to clarify the best set of modes needed for the ocean description.

Finally, the aim of the model is to study the effects of specific physical interaction mechanisms between the ocean and the atmosphere on the mid-latitude climate, both at large and intermediate scales. The modular design of the code of the model is adapted to such purposes, with the possibility of implementing new components, such as oceanic active transport, time-dependent forcings, or salinity fields.

Code availability

Maooam v1.0 is freely available for research purposes in the Supplement and is also available at http://github.com/Climdyn/MAOOAM. In addition, the code is archived at http://dx.doi.org/10.5281/zenodo.47507. A version of the Lua implementation which is parallelized using MPI is also available at http://github.com/Climdyn/MAOOAM/tree/mpi. The parallelized version is archived at http://dx.doi.org/10.5281/zenodo.47510.

Formulae to compute the inner products

In the formulae of the inner products of the atmospheric modes, use the following helper functions: $$\begin{array}{}& B_{\mathrm{1}}(u,v,w)=& {\displaystyle \frac{w+v}{u}},& B_{\mathrm{2}}(u,v,w)=& {\displaystyle \frac{w-v}{u}},& & \mathit{\lambda }\left(r\right)=\mathrm{0}\left(r\text{even}\right)\text{or}\mathrm{1}\left(r\text{odd}\right),& S_{\mathrm{1}}(u,v,w,z)=& -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(zu+wv),& S_{\mathrm{2}}(u,v,w,z)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(wv-zu),& S_{\mathrm{3}}(u,v,w,z)=& -S_{\mathrm{1}}(u,v,w,z),& S_{\mathrm{4}}(u,v,w,z)=& S_{\mathrm{2}}(u,v,w,z).\end{array}$$ The same notation will be used in this appendix. In what follows, ${\mathit{\delta }}_{ij}$ is the Kronecker delta, so that ${\mathit{\delta }}_{ij}=\mathrm{1}$ if $i=j$, and 0 otherwise. Likewise, the function $\mathit{\delta }\left(x\right)$ used in this appendix is defined as $$\mathit{\delta }\left(x\right)=\left\{\begin{array}{ll}\mathrm{1},& \text{if}x=\mathrm{0}\\ \mathrm{0},& \text{otherwise.}\end{array}\right.$$ Using these functions, the various coefficients of the model are calculated, starting with the internal atmosphere coefficients.

Atmospheric coefficients

In the following, we consider the ordering of the basis function used in Eqs. ()–(). For the sake of clarity, we add an extra informative upper index specifying the type of the atmospheric function in the definitions below. However, the inner products are completely defined by the lower indices alone. The atmospheric functions are thus noted: $$F_i^{\mathit{\alpha }}(x^{\prime },y^{\prime })=\left\{\begin{array}{lc}\sqrt{\mathrm{2}}\cos \left(P_i{y}^{\prime }\right)& \mathrm{if}\mathit{\alpha }=A,\\ \mathrm{2}\cos \left(M_{i}n{x}^{\prime }\right)\sin \left(P_i{y}^{\prime }\right)& \mathrm{if}\mathit{\alpha }=K,\\ \mathrm{2}\sin \left(H_{i}n{x}^{\prime }\right)\sin \left(P_i{y}^{\prime }\right)& \mathrm{if}\mathit{\alpha }=L,\end{array}\right.$$ and the oceanic functions $${\mathit{\varphi }}_i(x^{\prime },y^{\prime })=\mathrm{2}\sin \left({\displaystyle \frac{H_{\text{o},i}n}{\mathrm{2}}}{x}^{\prime }\right)\sin \left(P_{\text{o},i}{y}^{\prime }\right).$$

The $a_{i,j}$ coefficients

These coefficients correspond to the eigenvalues of the Laplacian operator acting on the spectral expansion basis functions: $$a_{i,j}^{\mathit{\alpha },\mathit{\beta }}=⟨F_i^{\mathit{\alpha }},{\mathrm{\nabla }}^{\mathrm{2}}{F}_j^{\mathit{\beta }}⟩,\mathit{\alpha },\mathit{\beta }\in \mathit{\{}A,K,L\mathit{\}},$$ which are given for each case by $$\begin{array}{}& & a_{i,j}^{A,A}=-{\mathit{\delta }}_{ij}{P}_i^{\mathrm{2}},& & a_{i,j}^{K,K}=-{\mathit{\delta }}_{ij}(n^{\mathrm{2}}{M}_i^{\mathrm{2}}+P_i^{\mathrm{2}}),& & a_{i,j}^{L,L}=-{\mathit{\delta }}_{ij}(n^{\mathrm{2}}{H}_i^{\mathrm{2}}+P_i^{\mathrm{2}}).\end{array}$$

The $c_{i,j}$ coefficients

These coefficients are needed to evaluate the contribution of the $\mathit{\beta }$-terms, and only involve the $K$- and $L$-type base functions. $$c_{i,j}^{\mathit{\alpha },\mathit{\beta }}=⟨F_i^{\mathit{\alpha }},{\partial }_{x^{\prime }}{F}_j^{\mathit{\beta }}⟩,\mathit{\alpha },\mathit{\beta }\in \mathit{\{}K,L\mathit{\}}.$$ We have that $$\begin{array}{}& & c_{i,j}^{K,K}=c_{i,j}^{L,L}=\mathrm{0},& & c_{i,j}^{K,L}=M_i\mathit{\delta }(M_i-H_j)\mathit{\delta }(P_i-P_j)=-c_{j,i}^{L,K}.\end{array}$$

The $g_{i,j,k}$ coefficients

These coefficients are given by $$g_{i,j,k}^{\mathit{\alpha },\mathit{\beta },\mathit{\gamma }}=⟨F_i^{\mathit{\alpha }},J(F_j^{\mathit{\beta }},F_k^{\mathit{\gamma }})⟩,\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\in \mathit{\{}A,K,L\mathit{\}},$$ and the non-zero ones are given by $$\begin{array}{}& & \begin{array}{l}{g}_{i,j,k}^{A,K,L}=-{\displaystyle \frac{\mathrm{2}\sqrt{\mathrm{2}}}{\mathit{\pi }}}{M}_j\left({\displaystyle \frac{B_{\mathrm{1}}(P_i,P_j,P_k{)}^{\mathrm{2}}}{B_{\mathrm{1}}(P_i,P_j,P_k{)}^{\mathrm{2}}-\mathrm{1}}},-,{\displaystyle \frac{B_{\mathrm{2}}(P_i,P_j,P_k{)}^{\mathrm{2}}}{B_{\mathrm{2}}(P_i,P_j,P_k{)}^{\mathrm{2}}-\mathrm{1}}}\right),\\ \times \mathit{\delta }(M_j-H_k)\mathit{\lambda }(P_i+P_j+P_k),\end{array}& & \begin{array}{l}{g}_{i,j,k}^{K,K,L}=S_{\mathrm{1}}(P_j,P_k,M_j,H_k)\mathit{\left\{}\mathit{\delta }\right(M_i-H_k-M_j\left)\mathit{\delta }\right(P_i-P_k+P_j)\\ -\mathit{\delta }(M_i-H_k-M_j)\mathit{\delta }(P_i+P_k-P_j)\\ +\left[\mathit{\delta }\right(H_k-M_j+M_i)+\mathit{\delta }(H_k-M_j-M_i\left)\right]\mathit{\delta }(P_k+P_j-P_i)\mathit{\}}\\ +S_{\mathrm{2}}(P_j,P_k,M_j,H_k)\mathit{\left\{}\mathit{\delta }\right(M_i-H_k-M_j\left)\mathit{\delta }\right(P_i-P_k-P_j)\\ +\left[\mathit{\delta }\right(H_k-M_j-M_i)+\mathit{\delta }(M_i+H_k-M_j\left)\right]\\ \times \left[\mathit{\delta }\right(P_i-P_k+P_j)-\mathit{\delta }(P_k-P_j+P_i\left)\right]\mathit{\}},\end{array}& & \begin{array}{l}{g}_{i,j,k}^{L,L,L}=S_{\mathrm{3}}(P_j,P_k,H_j,H_k)\mathit{\left\{}\mathit{\delta }\right(H_k+H_j-H_i\left)\mathit{\delta }\right(P_k-P_j+P_i)\\ +\left[\mathit{\delta }\right(H_k-H_j-H_i)-\mathit{\delta }(H_k-H_j+H_i\left)\right]\mathit{\delta }(P_k+P_j-P_i)\\ -\mathit{\delta }(H_k+H_j-H_i)\mathit{\delta }(P_k-P_j-P_i)\mathit{\}}\\ +S_{\mathrm{4}}(P_j,P_k,H_j,H_k)\mathit{\left\{}\mathit{\delta }\right(H_k+H_j-H_i\left)\mathit{\delta }\right(P_k-P_j-P_i)\\ +\left[\mathit{\delta }\right(H_k-H_j+H_i)-\mathit{\delta }(H_k-H_j-H_i\left)\right]\\ \times \left[\mathit{\delta }\right(P_k-P_j-P_i)-\mathit{\delta }(P_k-P_j+P_i\left)\right]\mathit{\}}\mathrm{for}k\ge j\ge i,\end{array}\end{array}$$ where we have used the functions defined at the beginning of this appendix. All the other permutations can be obtained thanks to $$g_{i,j,k}^{\mathit{\alpha },\mathit{\beta },\mathit{\gamma }}=-g_{j,i,k}^{\mathit{\beta },\mathit{\alpha },\mathit{\gamma }}=g_{k,i,j}^{\mathit{\gamma },\mathit{\alpha },\mathit{\beta }}=g_{j,k,i}^{\mathit{\beta },\mathit{\gamma },\mathit{\alpha }}.$$

The $b_{i,j,k}$ coefficients

These coefficients are given by $$b_{i,j,k}^{\mathit{\alpha },\mathit{\beta },\mathit{\gamma }}=⟨F_i^{\mathit{\alpha }},J(F_j^{\mathit{\beta }},{\mathrm{\nabla }}^{\mathrm{2}}{F}_k^{\mathit{\gamma }})⟩,\mathit{\alpha },\mathit{\beta },\mathit{\gamma }\in \mathit{\{}K,L\mathit{\}}.$$ Therefore we obtain $$b_{i,j,k}^{\mathit{\alpha },\mathit{\beta },\mathit{\gamma }}=a_{k,k}^{\mathit{\gamma },\mathit{\gamma }}⟨F_i^{\mathit{\alpha }},J(F_j^{\mathit{\beta }},F_k^{\mathit{\gamma }})⟩=a_{k,k}^{\mathit{\gamma },\mathit{\gamma }}{g}_{i,j,k}^{\mathit{\alpha },\mathit{\beta },\mathit{\gamma }}.$$

The $s_{i,j}$ coefficients

These coefficients encode the inner products between the atmospheric and oceanic basis functions: $$s_{i,j}^{\mathit{\alpha }}=⟨F_i^{\mathit{\alpha }},{\mathit{\varphi }}_j⟩,\mathit{\alpha }\in \mathit{\{}A,K,L\mathit{\}},$$ which gives $$\begin{array}{}& & s_{i,j}^A=\mathrm{8}\sqrt{\mathrm{2}}{P}_{\text{o},j}{\displaystyle \frac{\mathit{\lambda }\left(H_{\text{o},j}\right)\mathit{\lambda }(P_{\text{o},j}+P_i)}{{\mathit{\pi }}^{\mathrm{2}}{H}_{\text{o},j}(P_{\text{o},j}^{\mathrm{2}}-P_i^{\mathrm{2}})}},& & s_{i,j}^K=\mathrm{4}{H}_{\text{o},j}{\displaystyle \frac{\mathit{\lambda }(\mathrm{2}{M}_i+H_{\text{o},j})\mathit{\delta }(P_{\text{o},j}-P_i)}{\mathit{\pi }(H_{\text{o},j}^{\mathrm{2}}-\mathrm{4}{M}_i^{\mathrm{2}})}},& & s_{i,j}^L=\mathit{\delta }(P_{\text{o},j}-P_i)\mathit{\delta }(\mathrm{2}{H}_i-H_{\text{o},j}).\end{array}$$

The $d_{i,j}$ coefficients

These coefficients are related to the forcing of the ocean on the atmosphere. They are given by the formula $$d_{i,j}^{\mathit{\alpha }}=⟨F_i^{\mathit{\alpha }},{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\varphi }}_j⟩=M_{j,j}{s}_{i,j}^{\mathit{\alpha }},\mathit{\alpha }\in \mathit{\{}A,K,L\mathit{\}},$$ where the $M_{j,j}$ are given by the eigenvalues of the Laplacian operator acting on the oceanic basis functions (see next section).

Oceanic coefficients

The $K_{i,j}$ coefficients

These coefficients are related to the forcing of the atmosphere on the ocean. They are given by $$K_{i,j}^{\mathit{\alpha }}=⟨{\mathit{\varphi }}_i,{\mathrm{\nabla }}^{\mathrm{2}}{F}_j^{\mathit{\alpha }}⟩=a_{j,j}^{\mathit{\alpha },\mathit{\alpha }}{s}_{j,i}^{\mathit{\alpha }},\mathit{\alpha }\in \mathit{\{}A,K,L\mathit{\}}.$$

The $M_{i,j}$ coefficients

These coefficients identify with the eigenvalues of the Laplacian acting on the oceanic basis functions: $$M_{i,j}=⟨{\mathit{\varphi }}_i,{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\varphi }}_j⟩=-{\mathit{\delta }}_{ij}(n^{\mathrm{2}}{H}_{\text{o},i}^{\mathrm{2}}/\mathrm{4}+P_{\text{o},i}^{\mathrm{2}}).$$

The $N_{i,j}$ coefficients

These coefficients are needed to evaluate the contribution of the $\mathit{\beta }$-terms and are given by $$\begin{array}{cc}& N_{i,j}& =⟨{\mathit{\varphi }}_i,{\partial }_{x^{\prime }}{\mathit{\varphi }}_j⟩& =-\mathrm{2}n{H}_{\text{o},i}{H}_{\text{o},j}{\displaystyle \frac{\mathit{\delta }(P_{\text{o},i}-P_{\text{o},j})\mathit{\lambda }(H_{\text{o},i}+H_{\text{o},j})}{\mathit{\pi }(H_{\text{o},j}^{\mathrm{2}}-H_{\text{o},i}^{\mathrm{2}})}}.\end{array}$$

The $O_{i,j,k}$ coefficients

These coefficients are given by $$O_{i,j,k}=⟨{\mathit{\varphi }}_i,J({\mathit{\varphi }}_j,{\mathit{\varphi }}_k)⟩$$ with $$\begin{array}{l}{O}_{i,j,k}={\displaystyle \frac{n}{\mathrm{2}}}\left[S_{\mathrm{3}}\right(P_{\text{o},j},P_{\text{o},k},H_{\text{o},j},H_{\text{o},k}\left)\mathit{\right\{}\left[\mathit{\delta }\right(H_{\text{o},k}-H_{\text{o},j}-H_{\text{o},i})\\ -\mathit{\delta }(H_{\text{o},k}-H_{\text{o},j}+H_{\text{o},i})\left]\mathit{\delta }\right(P_{\text{o},k}+P_{\text{o},j}-P_{\text{o},i})\\ +\mathit{\delta }(H_{\text{o},k}+H_{\text{o},j}-H_{\text{o},i})\left[\mathit{\delta }\right(P_{\text{o},k}-P_{\text{o},j}+P_{\text{o},i})\\ -\mathit{\delta }(P_{\text{o},k}-P_{\text{o},j}-P_{\text{o},i})\left]\mathit{\right\}}+S_{\mathrm{4}}(P_{\text{o},j},P_{\text{o},k},H_{\text{o},j},H_{\text{o},k})\\ \times \mathit{\left\{}\right[\mathit{\delta }(H_{\text{o},k}+H_{\text{o},j}-H_{\text{o},i})\mathit{\delta }(P_{\text{o},k}-P_{\text{o},j}-P_{\text{o},i})]\\ +\left[\mathit{\delta }\right(H_{\text{o},k}-H_{\text{o},j}+H_{\text{o},i})-\mathit{\delta }(H_{\text{o},k}-H_{\text{o},j}-H_{\text{o},i}\left)\right]\\ \times \left[\mathit{\delta }\right(P_{\text{o},k}-P_{\text{o},j}-P_{\text{o},i})-\mathit{\delta }(P_{\text{o},k}-P_{\text{o},j}+P_{\text{o},i}\left)\right]\mathit{\left\}}\right]\\ \mathrm{for}k\ge j\ge i.\end{array}$$

The $C_{i,j,k}$ coefficients

These coefficients are given by $$C_{i,j,k}=⟨{\mathit{\varphi }}_i,J({\mathit{\varphi }}_j,{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\varphi }}_k)⟩=M_{k,k}{O}_{i,j,k}.$$

The $W_{i,j}$ coefficients

These coefficients are related to the short-wave radiative forcing of the ocean and are given by $$W_{i,j}^{\mathit{\alpha }}=⟨{\mathit{\varphi }}_i,F_j^{\mathit{\alpha }}⟩=s_{j,i}^{\mathit{\alpha }},\mathit{\alpha }\in \mathit{\{}A,K,L\mathit{\}}.$$

Definition of the tensor ${\mathcal{T}}_{i,j,k}$

The system of non-dimensionalized ODEs for the model variables is encoded in the model tensor ${\mathcal{T}}_{i,j,k}$, of which the complete definition is given in this appendix. Tensor elements that are not listed below are equal to zero. To alleviate the notations, we use a shorthand notation for the indices of the different variables, $$\begin{array}{ccc}& & {\mathit{\psi }}_i=i& (\mathrm{1}\le i\le n_{\text{a}}),& {\mathit{\theta }}_i=i+n_{\text{a}}& (\mathrm{1}\le i\le n_{\text{a}}),\\ & {\mathrm{\Psi }}_i=i+\mathrm{2}{n}_{\text{a}}& (\mathrm{1}\le i\le n_{\text{o}}),\\ & {\mathrm{\Theta }}_i=i+\mathrm{2}{n}_{\text{a}}+n_{\text{o}}& (\mathrm{1}\le i\le n_{\text{o}}).\end{array}$$ Furthermore, we suppress the upper indices which indicate the atmospheric function types but are otherwise not needed to unambiguously specify the inner products.

Atmosphere equations

The components of the tensor for the atmosphere streamfunction are given by $$\begin{array}{cc}& & {\mathcal{T}}_{{\mathit{\psi }}_i,{\mathit{\psi }}_j,\mathrm{0}}=-{\displaystyle \frac{c_{i,j}{\mathit{\beta }}^{\prime }}{a_{i,i}}}-{\displaystyle \frac{k_d}{\mathrm{2}}}{\mathit{\delta }}_{i,j},& {\mathcal{T}}_{{\mathit{\psi }}_i,{\mathit{\theta }}_j,\mathrm{0}}={\displaystyle \frac{k_d}{\mathrm{2}}}{\mathit{\delta }}_{i,j},\\ & {\mathcal{T}}_{{\mathit{\psi }}_i,{\mathit{\psi }}_j,{\mathit{\psi }}_k}={\mathcal{T}}_{{\mathit{\psi }}_i,{\mathit{\theta }}_j,{\mathit{\theta }}_k}=-{\displaystyle \frac{b_{i,j,k}}{a_{i,i}}},\\ & {\mathcal{T}}_{{\mathit{\psi }}_i,{\mathrm{\Psi }}_j,\mathrm{0}}={\displaystyle \frac{k_d{d}_{i,j}}{\mathrm{2}{a}_{i,i}}}.\end{array}$$

The atmospheric temperature equations are determined by the tensor elements $$\begin{array}{cc}& & {\mathcal{T}}_{{\mathit{\theta }}_{\mathrm{1}},\mathrm{0},\mathrm{0}}={\displaystyle \frac{C_{\text{a}}^{\prime }}{\mathrm{1}-a_{\mathrm{1},\mathrm{1}}{\mathit{\sigma }}_{\mathrm{0}}}},& {\mathcal{T}}_{{\mathit{\theta }}_i,{\mathit{\psi }}_j,\mathrm{0}}={\displaystyle \frac{a_{i,j}{k}_d{\mathit{\sigma }}_{\mathrm{0}}}{\mathrm{2}{a}_{i,i}{\mathit{\sigma }}_{\mathrm{0}}-\mathrm{2}}},\\ & {\mathcal{T}}_{{\mathit{\theta }}_i,{\mathit{\theta }}_j,\mathrm{0}}=-{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{0}}\left(\mathrm{2},,c_{i,j},,{\mathit{\beta }}^{\prime },+,a_{i,j},(,k_d,+,\mathrm{4},,k_d^{\prime },)\right)-\mathrm{2}(S_{\text{B,a}}^{\prime }+{\mathit{\lambda }}_{\text{a}}^{\prime }){\mathit{\delta }}_{i,j}}{\mathrm{2}{a}_{i,i}{\mathit{\sigma }}_{\mathrm{0}}-\mathrm{2}}},\\ & {\mathcal{T}}_{{\mathit{\theta }}_i,{\mathit{\psi }}_j,{\mathit{\theta }}_k}={\displaystyle \frac{g_{i,j,k}-b_{i,j,k}{\mathit{\sigma }}_{\mathrm{0}}}{a_{i,i}{\mathit{\sigma }}_{\mathrm{0}}-\mathrm{1}}},\\ & {\mathcal{T}}_{{\mathit{\theta }}_i,{\mathit{\theta }}_j,{\mathit{\psi }}_k}={\displaystyle \frac{b_{i,j,k}{\mathit{\sigma }}_{\mathrm{0}}}{\mathrm{1}-a_{i,i}{\mathit{\sigma }}_{\mathrm{0}}}},\\ & {\mathcal{T}}_{{\mathit{\theta }}_i,{\mathrm{\Psi }}_j,\mathrm{0}}={\displaystyle \frac{k_d{d}_{i,j}{\mathit{\sigma }}_{\mathrm{0}}}{\mathrm{2}-\mathrm{2}{a}_{i,i}{\mathit{\sigma }}_{\mathrm{0}}}},\\ & {\mathcal{T}}_{{\mathit{\theta }}_i,{\mathrm{\Theta }}_j,\mathrm{0}}=s_{i,j}{\displaystyle \frac{\mathrm{2}{S}_{\text{B,o}}^{\prime }+{\mathit{\lambda }}_{\text{a}}^{\prime }}{\mathrm{2}-\mathrm{2}{a}_{i,j}{\mathit{\sigma }}_{\mathrm{0}}}},\end{array}$$ where we used the non-dimenionalized quantities $$\begin{array}{cc}& {\mathit{\beta }}^{\prime }=\mathit{\beta }L/f_{\mathrm{0}},\\ & C_{\text{a}}^{\prime }=C_{\text{a}}R/\left(\mathrm{2}{\mathit{\gamma }}_{\text{a}}{f}_{\mathrm{0}}^{\mathrm{3}}{L}^{\mathrm{2}}\right),\\ & {\mathit{\sigma }}_{\mathrm{0}}=\mathit{\sigma }\mathrm{\Delta }{p}^{\mathrm{2}}/\left(\mathrm{2}{L}^{\mathrm{2}}{f}_{\mathrm{0}}^{\mathrm{2}}\right),\\ & {\mathit{\lambda }}_{\text{a}}^{\prime }=\mathit{\lambda }/\left({\mathit{\gamma }}_{\text{a}}{f}_{\mathrm{0}}\right),\\ & S_{\text{B,o}}^{\prime }=\mathrm{2}{\mathit{ϵ}}_{\text{a}}{\mathit{\sigma }}_{\text{B}}{\left(T_{\text{o}}^{\mathrm{0}}\right)}^{\mathrm{3}}/\left({\mathit{\gamma }}_{\text{a}}{f}_{\mathrm{0}}\right),\\ & S_{\text{B,a}}^{\prime }=\mathrm{8}{\mathit{ϵ}}_{\text{a}}{\mathit{\sigma }}_{\text{B}}{\left(T_{\text{a}}^{\mathrm{0}}\right)}^{\mathrm{3}}/\left({\mathit{\gamma }}_{\text{a}}{f}_{\mathrm{0}}\right).\end{array}$$

Ocean equations

The components of the tensor for the ocean streamfunction are $$\begin{array}{cc}& & {\mathcal{T}}_{{\mathrm{\Psi }}_i,{\mathit{\psi }}_j,\mathrm{0}}=-{\mathcal{T}}_{{\mathrm{\Psi }}_i,{\mathit{\theta }}_j,\mathrm{0}}={\displaystyle \frac{K_{i,j}{d}^{\prime }}{M_{i,i}+\mathit{\gamma }}},& {\mathcal{T}}_{{\mathrm{\Psi }}_i,{\mathrm{\Psi }}_j,\mathrm{0}}=-{\displaystyle \frac{N_{i,j}{\mathit{\beta }}^{\prime }+{\mathit{\delta }}_{i,j}{M}_{i,i}(r^{\prime }+d^{\prime })}{M_{i,i}+\mathit{\gamma }}},\\ & {\mathcal{T}}_{{\mathrm{\Psi }}_i,{\mathrm{\Psi }}_j,{\mathrm{\Psi }}_k}=-{\displaystyle \frac{C_{i,j,k}}{M_{i,i}+\mathit{\gamma }}},\end{array}$$ with $\mathit{\gamma }=-L/L_{\text{R}}$, $d^{\prime }=d/f_{\mathrm{0}}$, and $r^{\prime }=r/f_{\mathrm{0}}$.

Finally, the equations for the ocean temperature are determined by $$\begin{array}{cc}& & {\mathcal{T}}_{{\mathrm{\Theta }}_i,\mathrm{0},\mathrm{0}}=C_{\text{o}}^{\prime }{W}_{i,\mathrm{1}},& {\mathcal{T}}_{{\mathrm{\Theta }}_i,{\mathit{\theta }}_j,\mathrm{0}}=W_{i,j}(\mathrm{2}{\mathit{\lambda }}_{\text{o}}^{\prime }+{\mathit{\sigma }}_{\text{B,a}}^{\prime }),\\ & {\mathcal{T}}_{{\mathrm{\Theta }}_i,{\mathrm{\Theta }}_j,\mathrm{0}}=-{\mathit{\delta }}_{i,j}({\mathit{\lambda }}_{\text{o}}^{\prime }+{\mathit{\sigma }}_{\text{B,o}}^{\prime }),\\ & {\mathcal{T}}_{{\mathrm{\Theta }}_i,{\mathrm{\Psi }}_j,{\mathrm{\Theta }}_k}=-O_{i,j,k},\end{array}$$ where the following non-dimensionalized quantities are used: $$\begin{array}{cc}& C_{\text{o}}^{\prime }=C_{\text{o}}R/\left({\mathit{\gamma }}_{\text{a}}{f}_{\mathrm{0}}^{\mathrm{3}}{L}^{\mathrm{2}}\right),\\ & {\mathit{\lambda }}_{\text{o}}^{\prime }=\mathit{\lambda }/\left({\mathit{\gamma }}_{\text{o}}{f}_{\mathrm{0}}\right),\\ & {\mathit{\sigma }}_{\text{B,o}}^{\prime }=\mathrm{4}{\mathit{\sigma }}_{\text{B}}{\left(T_{\text{o}}^{\mathrm{0}}\right)}^{\mathrm{3}}/\left({\mathit{\gamma }}_{\text{o}}{f}_{\mathrm{0}}\right),\\ & {\mathit{\sigma }}_{\text{B,a}}^{\prime }=\mathrm{8}{\mathit{ϵ}}_{\text{a}}{\mathit{\sigma }}_{\text{B}}{\left(T_{\text{a}}^{\mathrm{0}}\right)}^{\mathrm{3}}/\left({\mathit{\gamma }}_{\text{o}}{f}_{\mathrm{0}}\right).\end{array}$$

Information about the Supplement

The time evolution of the model dynamics is illustrated by a set of videos, which are available online at 10.5446/18552, 10.5446/18553, 10.5446/18554, 10.5446/18555, 10.5446/18556, 10.5446/18557, 10.5446/18558, 10.5446/18559, 10.5446/18560, 10.5446/18561, 10.5446/18562, and 10.5446/18563.

The Supplement related to this article is available online at doi:10.5194/gmd-9-2793-2016-supplement.

Acknowledgements

This work is partly supported by the Belgian Federal Science Policy Office under contract BR/121/A2/STOCHCLIM. The figures and videos have been prepared with the Matplotlib software . Edited by: O. Marti Reviewed by: T. Sengul and one anonymous referee