Comparison of stochastic parameterizations in the framework of a coupled ocean-atmosphere model

A new framework is proposed for the evaluation of stochastic subgrid-scale parameterizations in the context of MAOOAM, a coupled ocean-atmosphere model of intermediate complexity. Two physically-based parameterizations are investigated, the first one based on the singular perturbation of Markov operator, also known as homogenization. The second one is a recently proposed parameterization based on the Ruelle's response theory. The two parameterization are implemented in a rigorous way, assuming however that the unresolved scale relevant statistics are Gaussian. They are extensively tested for a low-order version known to exhibit low-frequency variability, and some preliminary results are obtained for an intermediate-order version. Several different configurations of the resolved-unresolved scale separations are then considered. Both parameterizations show remarkable performances in correcting the impact of model errors, being even able to change the modality of the probability distributions. Their respective limitations are also discussed.

These methods have all in common a rigorous mathematical framework. They provide promising alternatives to other methods such as the ones based on the reinjection of energy from the unresolved scale through backscatter schemes (Frederiksen and Davies (1997); Frederiksen (1999), see also Frederiksen et al. (2017) for a recent review) or on empirical stochastic modeling methods based on autoregressive processes (Arnold et al., 2013).
The usual way to test the effectiveness of a parameterization method is to consider a well-known climate low-resolution model over which other methods have already been tested. For instance, several methods cited above have been tested on the Lorenz'96 model (Lorenz, 1996), see e.g. Crommelin and Vanden-Eijnden (2008); Arnold et al. (2013); Abramov (2015) and Vissio and Lucarini (2016). These approaches have also been tested in more realistic models of intermediate complexity that possess a wide range of scales and possibly a lack of timescale separation 1 , like for instance the evaluation of the MTV parameterization on barotropic and baroclinic models (Franzke et al., 2005;Franzke and Majda, 2006). Due to the blooming of parameterization methods developed with different statistical or dynamical hypothesis 2 , new comparisons are called for.
In this work, we investigate two parameterizations in the context of the MAOOAM ocean-atmosphere coupled model (De Cruz et al., 2016), used as a paradigm for multiscale systems. It is a two-layer baroclinic atmospheric model coupled to a shallow water ocean. It has been shown to exhibit multiple timescales including a low-frequency variability which is realistic for the midlatitude O-A system (Vannitsem et al., 2015). It possesses also a wide range of behaviors depending on the chosen operating resolution (De Cruz et al., 2016). As such, it forms a nice framework for testing parameterization methods on ocean-atmosphere related problems.
The particular problem of the atmospheric impact on the ocean could be addressed in this context as in Arnold et al. (2003) and Vannitsem (2014). It is an elegant problem, with a nice timescale separation, and which dates back to the work of Hasselmann (Hasselmann, 1976). However, in the present framework, other arbitrary decompositions of the model are possible, to address other questions. For instance, one may ask the question of the influence of the fast atmospheric processes on the slower atmospheric mode as well as on the very slow ocean. This problem was already addressed by the authors in Demaeyer and Vannitsem (2017) for a particular decomposition of the atmospheric modes based on the existence of an underlying invariant manifold. The parameterization considered was the one proposed by Wouters and Lucarini (Wouters and Lucarini, 2012), and for this specific invariant-manifold configuration, the stochastic parameterization is greatly simplified. Here, we extend this study by considering also the MTV parameterization (Franzke et al., 2005) with different arbitrary configurations. This in particular allows us to study different cases with or without multiplicative noise. This paper is organized as follow: in Section 2, we introduce briefly the MAOOAM model and its time-evolution equations.
In Section 3, we review the parameterization methods we have applied to MAOOAM , and detail the model decompositions into resolved and unresolved components. The results obtained with these parameterizations, with different configurations, are presented in Section 4. Finally, the conclusions are provided in Section 5 and new work avenues are proposed. 1 Timescale separation, or the existence of a spectral gap, is a crucial ingredient over which numerous parameterization methods rely. 2 i.e. weak coupling hypothesis, large scale separation, . . .

The MAOOAM model
The Modular Arbitrary-Order Ocean-Atmosphere Model is a coupled ocean-atmosphere model for midlatitudes. It is composed of a two-layer atmosphere over a shallow-water ocean layer on a β-plane (De Cruz et al., 2016). The ocean is considered as a closed basin with no-flux boundary conditions, while the atmosphere is defined in a channel, periodic in the zonal direction and with free-slip boundary conditions along the meridional boundaries. The model incorporates both a frictional coupling (momentum transfer) and an energy balance scheme which accounts for radiative and heat fluxes transfers between the ocean and the atmosphere. This model is developed as a basic representation of the processes at play between the ocean and the atmosphere, and to emphasize their impact on the midlatitude climate and weather. In particular, it has been shown to display a prominent low-frequency variability (LFV) in a 36-dimensions model version (Vannitsem et al., 2015).
The dynamical fields of the model include the atmospheric barotropic streamfunction ψ a and temperature anomaly T a = 2 f0 R θ a (with f 0 the Coriolis parameter at midlatitude and R the Earth radius) at 500hPa, as well as the oceanic streamfunction ψ o and temperature anomaly T o . In order to compute the time evolution of these fields, they are expanded in Fourier series with a set of functions satisfying the aforementioned boundary conditions: F K M,P (x , y ) = 2 cos(M nx ) sin(P y ) F L H,P (x , y ) = 2 sin(Hnx ) sin(P y ) for the atmosphere and for the ocean, with integer values of M , H, P , H o , and P o and where x = x/L and y = y/L are the non-dimensional coordinates on the β-plane. These functions are then reordered as specified in De Cruz et al. (2016) and the fields are expanded as follow: with a term that allows for mass conservation in the ocean, but does not play any role in the dynamics.
It results into a set of ODEs for the coefficients Z = ({ψ a,i }, {θ a,i }, {ψ o,j }, {θ o,j }) i∈{1,...,na}, j∈{1,...,no} of the expansion: where N = 2n a +2n o is the total number of coefficients. This model includes thus forcings, linear interactions and dissipations, as well as quadratic interactions representing the advection terms. Note that the full model equations of MAOOAM include quartic terms for the temperature fields, but those terms have been linearized around equilibrium temperatures (Vannitsem et al., 2015).

Model decompositions and parameterizations
Now let us consider a more general system of ordinary differential equations (ODEs): where Z ∈ R N represents the set of variables of a given model. A parameterization of the model supposes first to define a decomposition of this set of variables into two different subsets Z = (X, Y ), with X ∈ R N X and Y ∈ R N Y . In general, this decomposition is made such that the subsets X and Y have strongly differing response times: τ Y τ X (Arnold et al., 2003).
However, we will assume here that this constraint is not necessarily met, allowing for a more arbitrary split of the system variables. System (10) can then be expressed as: The timescale of the X sub-system is typically (but not always) much longer than the one of the Y sub-system, and it is sometimes materialized by a parameter δ = τ Y /τ X 1 in front of the time derivativeẎ . The X and the Y variables represent respectively the resolved and the unresolved components of the system. A parameterization is a reduction of the system (11) into a closed evolution equation for X alone such that this reduced system approximates the resolved component (Arnold et al., 2003). The term "accurately" here can have several meanings, depending on the kind of problem to solve. For instance we can ask that the closed system for X has statistics that are very close to the ones of the resolved component of system (11). We can also ask that the trajectories of the closed system remains as close as possible to the trajectories of the full system for long times, but this latter question will not be addressed in the present work.
More precisely a parameterization of the sub-system Y is thus a relation Ξ between the two sub-systems: which allows to effectively close the equations for the sub-system X while retaining the effect of the coupling to the Y sub-system. Since the work of Hasselmann (1976), various stochastic parameterizations have been introduced (Demaeyer and Vannitsem (2018) and Frederiksen et al. (2017) for reviews). They allow for a better representation of the variability of the processes and consider the relation (12) in a statistical sense. In this case, the aforementioned closed system for X becomes a stochastic differential equation (SDE). If the system (11) is rewritten: these SDEs are usually written: where the matrix σ, the deterministic function G and the vector of random processes ξ(t) have to be determined. We now detail in the rest of this section the two parameterizations that we will consider, namely the WL and the MTV approach.

The Wouters-Lucarini (WL) parameterization
The first one is based on the Ruelle's response theory (Ruelle, 1997(Ruelle, , 2009 and is valid when the resolved and the unresolved components are weakly coupled. This method proposed by Wouters and Lucarini (2012) is connected to the Mori-Zwanzig formalism (Wouters and Lucarini, 2013). It has already been applied to stochastic triads (Wouters et al., 2016;Demaeyer and Vannitsem, 2018), to the Lorenz'96 model (Vissio and Lucarini, 2016) and to the MAOOAM model in the 36-variables configuration displaying LFV (Demaeyer and Vannitsem, 2017). It considers the following decomposition: where ε is a small parameter accounting for the weak coupling and the functions F X and F Y account for all the terms involving only X and Y , respectively.
As said above, the parameterization is based on the Ruelle's response theory, which quantifies the contribution of the "perturbations" Ψ X and Ψ Y to the invariant measure ρ of the fully coupled system (13) as: where ρ 0 is the invariant measure of the uncoupled system. The two measures ρ and ρ 0 are supposed to be existing, well-defined Sinai-Ruelle-Bower (SRB) measures (Young, 2002). As a result of Eq. (17), the parameterization is based on three different terms having a response similar, up to order two, to the couplings Ψ X and Ψ Y : where is an averaging term. ρ 0,Y is the measure of the uncoupled systemẎ = F Y (Y ). The term M 2 (X, t) is a correlation term: where ⊗ is the outer product, is the centered perturbation and φ s X , φ s Y are the flows of the uncoupled systemẊ = F X (X) andẎ = F Y (Y ). The M 3 term is a memory term: involving the memory kernel All the averages are thus taken with ρ 0,Y , and the terms M 1 , M 2 and M 3 are derived (Wouters and Lucarini, 2012) such that their responses up to order two match the response of the perturbations Ψ X and Ψ Y . Consequently, this ensures that for a weak coupling, the response of the parameterization (18) on the observables will be approximately the same as the response of the full coupled system.

The Majda-Timofeyev-Vanden Eijnde (MTV) parameterization
This approach is based on the singular perturbation methods that were developed for the analysis of the linear Boltzmann equation in an asymptotic limit (Grad, 1969;Ellis and Pinsky, 1975;Papanicolaou, 1976;Majda et al., 2001). These methods are applicable for parameterization purposes if the problem can be cast into a backward Kolmogorov equation. Additionally, the procedure, described in Majda et al. (2001), requires some assumptions on the timescales of the different terms of the system (13). In term of the timescale separation parameter δ = τ Y /τ X , the fast variability of the unresolved component Y is considered of order O(δ −2 ), and the coupling terms are acting on an intermediary timescale of order O(δ −1 ): The intrinsic dynamics F X is thus the climate O(1) timescale which one is trying to improve with the parameterization. The term F Y is defined based on which terms are considered as part of the fast variability. It is then assumed that this fast variability can be modeled as an Ornstein-Uhlenbeck process. The Markovian nature of the process defined by Eqs. (23)-(24) and its singular behavior in the limit of an infinite timescale separation (δ → 0) allow then to apply the method. More specifically, the parameter δ serves to distinguish terms with different timescales and is then used as a small perturbation parameter. In this setting, the backward Kolmogorov equation reads (Majda et al., 2001): where the function ρ δ (s, X, Y |t) is defined with the final value problem f (X): ρ δ (t, X, Y |t) = f (X). The function ρ δ can be expanded in term of δ and inserted in Eq. (25). The zeroth order of this equation ρ 0 can be shown to be independent of Y and its evolution given by a closed, averaged backward Kolmogorov equation (Kurtz, 1973): This equation is obtained in the limit δ → 0 and gives the sought limiting, averaged process X(t). The parameterization obtained by this procedure is detailed in Franzke et al. (2005) and takes the form (14). As stated above, the parameterization itself depends on which terms of the unresolved component are considered as "fast", and some assumptions should here be made. It is the subject of the next section.

Model decompositions
As MAOOAM is a model whose nonlinearities consist solely of quadratic terms, the decomposition of Eq. (9) into a resolved and an unresolved component can be done based on the constant, linear and quadratic terms of its tendencies, as in Majda et al. (2001) and Franzke et al. (2005) : The vectors H are the constant terms of the tendencies. The matrices L give the linear dependencies through matrix-vector products: The symbol ⊗ is the outer product and is used to define matrices and tensors, e.g.: and ":" means an element-wise product with summation over the last and first two indices of its first and second arguments respectively 3 . For a rank-3 tensor and a matrix, it gives for example: As we have seen in Section 3.1, the two parameterization methods rely on a weak coupling between the components for the WL method and on a clear three-stages timescale separation for the MTV method. Moreover, as we will see below, their decomposition of the equations (27)-(28) is not necessarily the same, depending on modeling choices.

Decomposition of the resolved component
For both parameterizations, the decomposition of the X component is the same, and we have: and

Decompositions of the unresolved component
The definition of F Y and Ψ Y is not necessarily the same for both, but it is of particular importance since it is the measure of the system whose tendencies are given by F Y (Y ) over which the averages are performed (Demaeyer and Vannitsem, 2018).
Indeed, in the framework of the WL method, it is always the measure of the intrinsic Y -dynamics: that is used to perform the averaging.
In the framework of the MTV method, the measure of the singular systemẎ = 1 δ 2 F Y (Y ) is used for the averaging and it is usually assumed that the quadratic Y -terms of the unresolved component tendencies represent the fastest timescale of the system (see Majda et al. (2001); Franzke et al. (2005)): and are the ones over which the averaging has to be done. The others terms 4 belong then to Ψ Y . It is interesting to note that there is no a priori justification for this assumption. For instance, the decomposition of the unresolved dynamics could be based on the full intrinsic dynamics (as in the WL method) and F Y would then be given by the expression (34). We consider these two different assumptions for the MTV in the following.
The MTV and WL parameterizations described above are presented in more details in the appendices A and B, respectively.

Noisy model
To take into account model errors or the impact of smaller scales, the present implementation of MAOOAM allows for the addition of Gaussian white noise in each components: resolved and unresolved, for both the ocean and the atmosphere. It also includes the timescale separation parameter δ of the MTV framework (see Eqs. (23) and (24)) and the coupling parameter ε of the WL framework (see Eqs. (15) and (16)). The full equation (11) then reads: and the noise amplitude can hence be different for each components. The dW 's vectors are standard Gaussian white noise vectors. In the framework of stochastic parameterization, the presence of noise is sometimes required to smooth the averaging measure (Colangeli and Lucarini, 2014) or to increase the small-scale variability to address the problem of "dead" scales. The ) function can be specified by either Eq. (35) or (34) (only (34) for the WL parameterization).
This flexible setup allows for testing both parameterizations in the same framework, with or without additional noise. We now present the results obtained by applying them to the MAOOAM model.

Results
The relative performance and the interesting features of the parameterizations described in the previous section require to consider multiple versions and resolutions of the model. We thus shall consider in the following two different resolutions. The first one is the 36-variables model version considered in Vannitsem et al. (2015) for which the maximum values for M , H and P in Eqs.
(1)-(3) is 2. It corresponds to a "2x-2y" resolution for the atmosphere as referred to in De Cruz et al. (2016). For the ocean, the maximum values for H o and P o in Eq. (4) are respectively 2 and 4, and the resolution is therefore noted "2x-4y" in the same notation system. The model version for this first case is thus noted "atm-2x-2y oc-2x-4y" and we shall abbreviate it the acronym "VDDG" from the name of the authors in Vannitsem et al. (2015). The other resolution that we shall consider is "atm-5x-5y oc-5x-5y" that can be abbreviated "5 × 5" since we include Fourier modes up to the wavenumber 5 in both the ocean and the atmosphere. In this latter model version, the model possesses 160 variables.
We shall also consider different sets of parameters for the model configuration. To control the results obtained with the code implementation provided as supplementary material, we will compare our results with those obtained in Demaeyer and  Cruz et al. (2016) and Demaeyer and Vannitsem (2017).  with the first set of parameter defined therein. We will refer to this set of parameter DV2017. A second set of parameters from De Cruz et al. (2016) is also considered and will be denoted DDV2016. For both "DV2017" and "DDV2016", MAOOAM has been shown to depict coupled ocean-atmosphere low-frequency (LFV) variability. In consequence, we will also consider a third parameter set where no LFV is present. The LFV has been removed by reducing by an order of magnitude the wind friction parameter d between the ocean and the atmosphere. The coupled-mode dynamics then disappears. This parameter set is referred as noLFV. In addition, the parameters δ and ε appearing in Eqs. (36)-(39) will be set to 1, meaning that we consider the natural timescale separations and coupling strengths of the model. Nevertheless, the study of the impact of these parameters is important (Demaeyer and Vannitsem, 2018) and should be carried out in forthcoming works.
Finally, for a given resolution and a given parameter set, multiple different parameterization experiments can be designed, by using for example the unresolved dynamics (35) or (34) for the MTV parameterization. However, to simplify the study and to be able to compare both the MTV and the WL parameterizations, we will here consider only the dynamics of (34) to perform the averages. Another way of defining different parameterization experiments is by changing the resolved and unresolved components. We shall detail these different experiments and the results obtained in the following subsections.
Unless otherwise specified, the subsequent results were obtained by considering long trajectories lasting 2.8 × 10 6 days 7680 years, generated with a timestep of 96.9 s, and after an equivalent transient time. Such long trajectories were needed to be able to sample sufficiently the long timescales present in the model ( 20-30 years).
We now propose different partitions of these variables into different resolved and unresolved sets. We will detail, for each of these parameterization experiment, the results obtained for the different aforementioned parameter sets.

Parameterization based on the invariant manifold
We first consider a parameterization based on the presence in the VDDG model of an genuine invariant manifold. As stated in Demaeyer and Vannitsem (2017), this manifold is due to the existence of a subset of the Fourier modes that is let invariant by the Jacobian term of the partial differential equation of the system: In the same spirit, we consider the atmospheric variables outside of this invariant manifold to be unresolved, all the other variables being resolved. The modes F 2 , F 3 , F 4 , F 7 and F 8 are outside of this invariant set, and therefore the variables ψ a,2 , ψ a,3 , ψ a,4 , ψ a,7 , ψ a,8 , θ a,2 , θ a,3 , θ a,4 , θ a,7 , θ a,8 are considered as unresolved. Using the same parameterizations, parameters and methods as in Demaeyer and Vannitsem (2017), we should recover the same results with our current implementation 5 . This would thus provide a first mandatory check for the correctness of the current implementation. We show the results obtained on Fig. 2 for the "DV2017" parameter set, where the probability density functions (PDFs) of three important variables of the dynamics are displayed. These three variables are ψ a,1 , ψ o,2 and θ o,2 , and were shown in Vannitsem and Ghil (2017) to account for respectively 18%, 42% and 51% of the variability in some key reanalysis 2-dimensional fields over the North-Atlantic ocean. In Vannitsem et al. (2015), it was also shown that these three variables are dominant in the LFV pattern found in the VDDG model version of MAOOAM . In addition to the PDFs of the full coupled system (13) (34) is a multidimensional Ornstein-Uhlenbeck, for which the measure is well known and thus analytical expressions for the correlations can be provided to the code. This analytical solution is also used in Demaeyer and Vannitsem (2017) and is one of the setting that we used to compute the parameterization. This setting is also compared with a setting using a numerical least-square regression method to obtain the correlation of variables of the system (34), assuming that these are of the simplified form where t is the time lag and τ is the decorrelation time. The results obtained for these two different settings of the correlation (analytical vs. numerical) are shown in Fig. 2 with "check" indicating the results obtained with the analytical expressions for the correlations. The latter expressions improve the correction of the dynamics for both the WL and the MTV methods and we recover the same results as in Demaeyer and Vannitsem (2017). 6 On the other hand, in the case where the correlations are specified by the results of the numerical analysis, the parameterizations perform less well, indicating that these methods can be very sensitive to a correct representation of the correlation structure of the unresolved dynamics.
Both the MTV and the WL methods correct better the oceanic variables whereas the atmospheric variables seem to display too different dynamical behavior between the coupled and uncoupled systems which are difficult to correct. As stated in Demaeyer and Vannitsem (2017), it may be due to the huge dimension of the unresolved system, with half of the atmospheric mode for the X and Y components. The densities of the full coupled system (13) and of the uncoupled system are depicted for the "DV2017" parameter set, as well as the ones of the MTV and WL parameterizations and their verifications. The "check" label indicates that the correlations used for the parameterizations were analytical, while for the other they were obtained by numerical analysis.
of the advection operator (42). It leads to simplified coupling and is thus computationally advantageous. Within this framework, an adaptation based on the next order of both parameterization methods or the consideration of other parameterizations could lead to a very efficient correction for both the ocean and atmosphere.
As the present implementation allows for an arbitrary selection of the resolved-unresolved components, we shall now consider cases of the VDDG model version with different unresolved components.

Parameterization of the wavenumber 2 atmospheric variables
We consider now a smaller set of unresolved variables Y , composed of the two higher resolution modes of the model: The unresolved variables Y are thus the following ones: ψ a,9 , ψ a,10 , θ a,9 and θ a,10 .
A first comment on the results obtained with that particular configuration is that the WL method seems to be unstable for all the parameter sets investigated. As stated in the Appendix section B4, the equation (18) is integrated with a Heun algorithm where the term M 3 is considered constant during the corrector-predictor steps. As this could lead to instabilities, a backward differentiation formulae (BDF) designed for stiff integro-differential equations has been tested to integrate (18) (Lambert, 1991;Wolkenfelt, 1982), but the instabilities were still present, leading us to suspect a genuine instability in the parameterized system. Indeed, we found that it is the cubic term M(s) in the memory term M 3 (see Eq. (B30)) which causes the divergence, in particular the interactions with the F 4 mode. On the contrary, the MTV parameterization is stable and performs well, despite having a similar cubic term. It indicates that the correlations induced by the memory term are a possible cause of the instability.
More research effort to understand this stability issue is needed.
The quality of the solutions obtained with the MTV method alone, applied to the system with the three parameter sets of Table 1, is evaluated using the Kullback-Leibler divergence between the distribution of p(x) the full coupled system and the distribution q(x) of the parameterized system that is tested. It is related to the information lost when this parameterization is used instead of the full system. The bigger the divergence, the greater is the amount of information lost. For clarity, we show only the divergence of the marginal distributions averaged by component. In every case and for every component, the MTV parameterization reduces dramatically the divergence and thus corrects well the models, making it clearly a better choice than the "Uncoupled" dynamics. These divergences are compiled in particular, one can notice on Fig. 3 a change of the distributions modality due to the parameterization. It raises the question about the mechanism of this rather drastic modification: Is it the noise that changes the dynamics or does the noise simply trigger and shift a bifurcation of the unresolved system, which then induces the change of modality? To give some insight about this latter possibility, we considered the "noLFV" parameter set and increased the ocean-atmosphere wind stress coupling parameter d to 5 × 10 −9 . With this change of parameter, the system now lies nearby the Hopf bifurcation generating the long periodic orbit which forms the backbone of the LFV. In other words, the system does not exhibit a LFV but a small parameter perturbation could induce it. As seen on Fig. 7, the MTV parameterization then induces a LFV which is not present in both the resolved and the full system. In that case, the parameterization wrongfully "triggers" a bifurcation and thus leads to a false reduced dynamics. This interesting issue, also considered by recent works on the effect of the noise on models dynamics, will be further discussed in the conclusion of this paper (see Section 5).
The impact of the MTV parameterization on the correlation is also particularly important, as seen on Figs. 4 and 6, correcting the covariance (the value at the time lag t = 0) and inducing or suppressing oscillations. It shows that the parameterization not only affects the attractors and the climatologies of the models, but also the temporal dynamics.
Additional marginal PDF and autocorrelation function figures are available for each variable of the system in the supplementary material. The Kullback-Leibler divergences (45) are available as well.
Since the WL parameterization does not work in the current case, we cannot properly compare both methods. To do so, we shall consider two other parameterization configurations in the following sections.

Parameterization of the wavenumber 2 atmospheric baroclinic variables
Let us now consider that only the two baroclinic variables θ a,9 and θ a,10 are unresolved. This particular case allows for the comparison of the WL and the MTV methods. Indeed, in this configuration, the destabilizing cubic interactions with the mode F 4 are suppressed and the WL parameterization is now stable. The results are summarized in Table 3 where the Kullback-Leibler divergence (45) of the parameterizations marginal distributions with respect to the full coupled system ones are indicated.
The divergence of the uncoupled system with respect to the full system is less pronounced than in the previous section. The    parameterization. It is shown for the parameter set "noLFV" but with the wind-stress parameter set to d = 0.5 × 10 −9 and for which the MTV parameterization depicts a false LFV.
Barotropic  Table 3. Component averaged Kuhlback-Leibler divergence with respect to the distributions of the full coupled system, in the case of the wavenumber 2 atmospheric baroclinic variables parameterization.
-The parameterizations correct quite well the ocean components, except for the MTV parameterization of the baroclinic component in the noLFV case. The MTV parameterization is better in the DDV2016 case, and the WL one is better in the two other cases.
-The barotropic component of the atmosphere is only well corrected in the DV2017 case. Additionally, the MTV fails to correct also the baroclinic PDFs in the two other cases. In fact, the MTV method seems only to performs well in the DV2017 case. Looking to the divergence for every variables (see the supplementary material), we note that those under-performance are due to the incorrect representation of the small-scale wavenumber 2 barotropic variables, namely the ψ a,6 -ψ a,10 (ψ a,5 -ψ a,10 ) variables for the WL (MTV) method.
Interestingly, the PDFs of the dominant variables ψ a,1 , ψ o,2 and θ o,2 are well corrected as shown in Fig. 8 for the DV2017 case and slightly better corrected in Fig. 10 for the DDV2016 case. We remark that in general the WL parameterization is better at correcting the LFV than the MTV one, but the situation can be more complicated, like in Fig. 10(b), where the WL parameterization captures well one mode of the distribution, and the MTV parameterization captures well another mode.
Regarding the correlation functions, a first general comment is that in this configuration, the decorrelation time of the large scales (mode F 1 ) does not appear to be significantly affected by the absence of the unresolved variables. On the other hand, the impact on the other modes is noticeable, and both parameterizations improve in general the correlations of the resolved variables. Finally, a general observation for all the parameter sets is the bad correction of the variables ψ a,9 and ψ a,10 (see Fig. 9(d)). This is not surprising since these variables are strongly coupled to the unresolved θ a,9 and θ a,10 variables, and have roughly the same decorrelation timescale. It also may explain the poor scores of the atmospheric zonal wavenumber 2 modes noted above. It implies that by parameterizing baroclinic variables at certain scales, one should not expect these methods to perform well for the barotropic variables at the same scale.

The parameterization of the atmospheric wavenumber 2 and F 4 modes
Another possibility to test the WL parameterization is to remove the aforementioned cubic interactions by considering that the wavenumber 2 modes and the meridional mode   are unresolved. In this case the variables, ψ a,4 , θ a,4 ,ψ a,9 , ψ a,10 , θ a,9 and θ a,10 are considered as unresolved. The WL method is then stable and this configuration allows to compare both parameterizations.
The global averaged Kuhlback-Leibler divergence are given in Table 4 for two parameter sets. These results show a great disparity. For the DV2017 parameter set, the WL method does a better job at correcting the atmosphere while the MTV method corrects better the oceanic modes. Considering the other parameter set without a LFV, we notice that both methods have troubles at improving the oceanic component. Note that with this particular wind stress reduced configuration, this component dynamics is less important since it can basically be modeled as an atmospheric fluctuations integrator (Hasselmann, 1976). The MTV method has also troubles at modeling the atmospheric dynamics correctly.
The PDFs of the three main variables (see Fig. 11) show that both parameterizations induce a change of modality for ψ a,1 and a good correction of the LFV. We note that the MTV method corrects better the LFV signal in the dominant oceanic modes ψ o,2 and θ o,2 (as also shown in Table 4).

The 5 × 5 model version
A higher-resolution test with the MTV method has also been performed by considering the 5 × 5 resolution system discussed at section 4. In this version, the resolution goes up to wavenumber 5 in every spatial direction and in every component. We do a quite drastic reduction of the dimensionality of the system by considering every mode above the wavenumber 2 in the atmosphere as unresolved. Consequently, the uncoupled system is an "atm-2x-2y oc-5x-5y" model version. The result of the parameterization on the dominant ψ o,2 and θ o,2 is shown on Fig. 12. The MTV method significantly corrects the 2D marginal distribution compared to the uncoupled model, as seen on the anomaly plots. Therefore, the oceanic section of the attractor of the MTV system is closer to the full coupled system than the uncoupled one. The atmosphere is not very well corrected by the parameterization. This could have multiple cause. First, the trajectory computed at this resolution is shorter (roughly one million days) and thus some long term equilibration of the dynamics might still take place. In this case, much longer computation might be undertaken. Secondly, the 5 × 5 resolution represents a particular case, where the atmosphere's Rhines scale is attained. At this scale, the dynamics is quite different from the higher or lower resolution model version De Cruz et al. (2016). If this has an influence on the parameterization performances, then it calls for even higher-resolution tests for confirmation. Finally, the parameterization may genuinely well correct the ocean while having trouble at improving the atmosphere representation, like for instance with the invariant manifold parameterization described in section 4.1.1 (see Fig. 2(c)).

Conclusions
In the present work, we have introduced a new framework to test different stochastic parameterization methods in the context of the ocean-atmosphere coupled model MAOOAM. We have implemented two methods: a homogenization method based on the singular perturbation of Markovian operator and known as MTV (Majda et al., 2001;Franzke et al., 2005), and a method based on the Ruelle's response theory (Wouters and Lucarini, 2012)  Additionally, we have found that these methods are able to change correctly the modality of the distributions in some cases.
However, in some other cases, they can also trigger a LFV that is absent from the full system. This leads us to underline the profound impact that a stochastic parameterization, and noise in general, can have on models. For instance, Kwasniok (2014) has shown that the noise can influence the persistence of dynamical regimes and can thus have a non-trivial impact on the PDFs, whose origin is the modification of the dynamical structures by the noise. In the present study as well, the perturbation of the dynamical structures by the noise is a very plausible explanation for the observed change of modality, and for the good performances of the parameterizations in general. However, if these perturbations can lead to a correct representation of the full dynamics, they can also generate regimes that are not originally present. This may happen near a sensitive bifurcation, as it is the case with a wrong LFV regime, which develops around a long-period periodic orbit (Vannitsem et al., 2015) arising through a Hopf bifurcation. If the main parameter (here for instance the wind stress parameter d) is set close to this bifurcation, the noise may thus trigger it while it is not present in the full system (Horsthemke and Lefever, 1984).
The MTV parameterization has also been tested in a intermediate-order version of the model, showing that this parameterization reduces the anomaly of the PDF of the two dominant oceanic modes. The atmospheric modes are however less well corrected. In this case the number of modes that are removed is large and one can wonder whether reducing this number or increasing the resolution will help. More work needs to be done to assess the impact of the parameterizations on higher-order version of the models in the future.
The MTV method is simpler, less involved than the WL one, with no memory term estimations needed, and thus no integrals are being computed at every step. The memory term could however be Markovianized as in Wouters et al. (2016). The interest of the WL method is that it requires only a weak coupling between the resolved and the unresolved components, but no timescale separation between them, which is a desirable feature for a parameterization. Consequently, an improvement of the present framework would be to make the MTV method less dependent on the timescale separation. One way to do that is to consider the next order in δ, the timescale separation parameter. This can be done effectively using the so-called Edgeworth expansion formalism, as shown by Wouters and Gottwald (2017). A next step would thus be to compute this expansion for the present coupled ocean-atmosphere system.
Finally, it would be interesting to consider that the unresolved dynamics used to perform the averaging may have non-Gaussian statistics. In the present work, as stated in Appendix A, the statistics of the neglected variables are assumed to come from a Gaussian. However, depending on the terms regrouped in this discarded part of the system, the statistics may well be non-Gaussian, and the resulting parameterization developed in this Appendix is then only an approximation. Indeed, unresolved variables with different linear damping terms can lead to such non-Gaussianity (see Sardeshmukh and Penland (2015)). An example of it is the MTV parameterization of the wavenumber 2 in the 36-variable system (see section 4.1.2), where the statistics are weakly non-Gaussian, because the linear damping of the baroclinic unresolved variables is not the same as the barotropic ones. But other configurations are concerned as well. Taking this into account could probably improve further the results obtained in the present study.
Code availability. The source code for MAOOAM v1.2 is available on GitHub at http://github.com/Climdyn/MAOOAM. The stochastic parameterization code is available in fortran in the stoch branch of this repository, in the fortran subdirectory. It is also provided as supplementary material to the present article.

Appendix A: The MTV method
We now consider Eqs. (23)-(24) and assume that the dynamics induced by the order 1/δ 2 term of the Y -dynamics can be replaced by (or behaves like) a multidimensional Ornstein-Uhlenbeck process: where W Y is a vector of independent Wiener processes. This assumption is described in Majda et al. (2001) as a working assumption for stochastic modeling. This assumption is related to other underlying assumptions, i.e. that the dynamics of this singular terṁ is ergodic and mixing with an integrable decay of correlation (Franzke et al., 2005). We will also assume that this dynamics generates Gaussian distributions, such that the odd-order moments vanish and that the even-order moments are related to the second-order one. In this case, the singular perturbation theory mentioned in section 3.1.2 gives the following result (see Papanicolaou (1976)) in the limit δ 1 for the parameterization of the resolved component: and σ 2 (X) = P(X) with and where W is a vector of independent Wiener processes. Note that the terms G(X) and P( The measureρ Y is the one of the system (A2) or the measure of the Ornstein-Uhlenbeck process (A1) replacing it. Similarly, is the result of the evolution of the flow φ s Y of the system (A2) or the non-stationary solution of the Ornstein-Uhlenbeck process: In section A1, we sketch the derivation of Eq. (A3), assuming that the Y -dynamics is an Ornstein-Uhlenbeck process. Furthermore, if the Y -dynamics is an Ornstein-Uhlenbeck process like (A1), the results (A3)-(A7) can also be expanded to give formula in terms of the covariance and correlation matrices of this process. Therefore, assuming that the dynamics of system (A2) has Gaussian statistics, its measure can be used asρ Y and these formula for the process (A1) can then be applied directly using the covariance and correlation matrices of system (A2). It gives a way to practically implement the parameterization as we detail it in section A2.

A1 Brief sketch of the parameterization derivation
As stated in Section 3.1, the MTV parameterization is based on the singular perturbation theory of Markovian operator. We follow the derivation proposed in Majda et al. (2001) and assume that the singular O(1/δ 2 ) term is an Ornstein-Uhlenbeck process like (A1). The backward Kolmogorov equation for Eqs. (23)-(24) for the density ρ δ (s, X, Y |t), where t is the final time, reads: with the final condition given by some function of X: The operator appearing in Eq. (A10) are given by: where T denotes the matrix transpose operation. Since δ 1, we can now perform an expansion procedure and recast it in Eq. (A10), equating term by term at each order, to get: The solvability condition for theses equations is that the right-hand side belong to the range of L 1 and thus that their average with respect to the invariant measure of (A1). The first solvability conditions is obviously satisfied but the the second is not necessarily satisfied 7 since where P is the expectation with respect to the measure of the process (A1) and σ Y is its covariance matrix. It indicates that 1/δ effects have to be taken into account in the parameterization. This can be done by following the method to treat "fast-waves" effects described in Majda et al. (2001). From here we thus depart from the standard derivation, to include those 1/δ effects in the parameterization. We need to assume that two timescales are present in the parameterization and consider them in the Kolmogorov equation (A10). Hence, we modify it explicitly by setting: and ρ δ (s, τ, X, Y |t) = ρ 0 (s, τ, X, Y |t) + δ ρ 1 (s, τ, X, Y |t) + δ 2 ρ 2 (s, τ, X, Y |t) .

A2 Practical implementation
We will now derive explicitly the MTV parameterization for the system (27)-(28). For the time of the derivation, we will again assume that the Y -dynamics is given by the process (A1), but we suppose that the final results apply as well with the measure, covariance and correlation matrices of system (A2).
We consider first the case defined by Eq. (35) where the intrinsic dynamics is considered to be given by the quadratic terms of the tendencies alone. We define the matrix and note that = 0 and where σ Y is the covariance matrix of the process (A1). These results can be explicitly obtained using the non-stationary solution (A9). We also define which is thus a rank-4 tensor. Note that both Σ and Σ 2 are O(δ 2 ) objects. With those definitions, it can be shown that the parameterization (A3) becomeṡ The product " : · " is similar to the definition (31) but with a rank-4 tensor and three variables : All the products involved concern the last and the first indices of respectively their first and second arguments. M and V are rank-4 tensors, U is a rank-3 tensor and the Q (i) 's are matrices. The H (i) 's are vectors, the L (i) 's are matrices and the B (i) 's are rank-3 tensors. The formula of these quantities are given below: where · is the product and summation over the last and the first indices of respectively its first and second arguments 8 . The product : is defined as in Eq. (31). The tensors whose indices are indicated define a lower rank tensor, for instance the rank-3 tensor B XXY i hence noted defines a matrix whose elements are B XXY i j k . The symbol then defines the following permuted product and summation of two given rank-4 tensors C and D: With the notable exception of H (0) , H (3) and L (0) , one can easily check that the formulation (A39)-(A52) gives back the formulas of the appendix in Franzke et al. (2005) when σ Y is diagonal. It is these particular tensors that are implemented and computed by the code provided with the present article as supplementary material. They are computed using the covariance matrix σ Y and the integrated correlation matrices Σ and Σ 2 as input. The formulas presented here are valid for an Ornstein-Uhlenbeck process, but as stated above, the covariance and correlations of the dynamics (A2) can be used directly as well,

A3 Technical details
The equation (A32) is integrated with a stochastic Heun algorithm described in Greiner et al. (1988) and which converges toward the Stratonovitch statistic (Hansen and Penland, 2006). In particular, the R(X) and G(X) terms and the P 1 (X) and P 2 matrix are computed by performing sparse-matrix products. The square root of the matrix P 2 is computed once at initialization with a singular-value decomposition (SVD) method to obtain the matrix σ (2) . The square root of the state-dependent matrix P 1 (X) is computed each mnuti time 9 to obtain the matrix σ (1) (X), again with a SVD decomposition. In general, in the present work, we have set mnuti equal to the Heun algorithm timestep.

Appendix B: The WL method
The Wouters-Lucarini method is considered with the decomposition (34) and consists of three terms (Wouters and Lucarini, 2012) : X = F X (X) + ε M 1 (X) + ε 2 M 2 (X, t) + ε 2 M 3 (X, t) which we detail now in the following. 8 Thus for vectors and matrices it is the standard product. 9 The notation mnuti indicates a parameter of the code implementation, please read the code documentation for more information.
where div(a, b) is the integer division of a by b and n is the second dimension of A. The vector vec (Y ⊗ Y − σ Y ) is thus a N 2 Y length vector. The object a X ij in Eq. (B7) can then be rewritten as a matrix: where we define Ψ X,1,1 (X) · a X 1j = 1 · a X 1j = a X 1j = L XY mat B XY Y and The symbol mat denotes an operation transforming rank-3 tensors into matrices, e.g.: where n = N Y is the second dimension of B XY Y . With this notation, we have thus for instance: Here mat B XY Y is thus a N X × N 2 Y matrix. The notation vec and mat introduced here are defined such that the expression (B7) can now be written in the compact form with matricial products: Ψ X (X, Y ) = Ψ X,1 (X) · a X ij · Ψ X,2 (Y ).

B3 The M 3 term
It is the memory term which is defined as where X s = φ s X (X) is the result of the time-evolution of the flow of the uncoupled X-dynamics given by the systemẊ = F X (X). Detailing each term in this expression, M 3 can be rewritten as (1) because the dynamics of (34) is assumed to be Gaussian and where and the expression of this last definition can be inferred from Eq. (A30).

B4 Technical details
The tendencies appearing in Eq. (B1) are computed according to the same Heun algorithm described in Sec. A3. A small difference with the latter here is that the term M 3 , due to its computational cost, is computed every muti time 10 and remains constant in between. The term M 1 is simply obtained by computing sparse-matrix products. The process σ(t) of the term M 2 is emulated by a Multidimentional n-th order Auto-Regressive process (MAR(n)) described in Penny and Harrison (2007). The parameters and the order n of this process can be assessed for instance by using the variational Bayesian approach described in this article. The integrand of the term M 3 is obtained by computing sparse-matrix products and integrating for each timestep the uncoupled X-dynamics forward, to obtainX s . The integral is then performed with a simple trapezoidal rule with a timestep muti over a time interval which is set by the parameter meml. This interval should be set so that the absolute value of the integrand has sufficiently decreased at the end of it. In practice, it is sufficient to set meml proportional to the longest decorrelation time of the uncoupled Y -dynamics. 10 The notation muti indicates a parameter of the code implementation, please read the code documentation for more information.