Stochastic parameterization of subgrid-scale processes: A review of recent physically-based approaches

TL;DR

This paper reviews recent physically-based stochastic parameterization methods for subgrid-scale processes in climate models, emphasizing their impact on accurately representing atmospheric variability and the importance of stability properties.

Contribution

It provides a comprehensive review of recent approaches and demonstrates their effectiveness using a stochastic triad system relevant to climate dynamics.

Findings

Stability of subgrid processes significantly affects parameterization performance.

Physically-based stochastic methods improve climate variability representation.

Stochastic triad system illustrates key dynamics of subgrid processes.

Abstract

We review some recent methods of subgrid-scale parameterization used in the context of climate modeling. These methods are developed to take into account (subgrid) processes playing an important role in the correct representation of the atmospheric and climate variability. We illustrate these methods on a simple stochastic triad system relevant for the atmospheric and climate dynamics, and we show in particular that the stability properties of the underlying dynamics of the subgrid processes has a considerable impact on their performances.

Figures

Equations

$$\dot{Z}=T(Z)$$

$$\left\{\begin{array}[]{lcl}\dot{X}&=&F(X,Y)\\ \dot{Y}&=&H(X,Y)\end{array}\right.$$

$$Y=\Xi (X,t)$$

$$\left\{\begin{array}[]{lcl}\dot{X}&=&F_{X}(X)+\Psi_{X}(X,Y)\\ \dot{Y}&=&F_{Y}(Y)+\Psi_{Y}(X,Y)\end{array}\right.$$

$$\dot{X}=F_X(X)+G(X,t)+\sigma (X)\cdot \tilde{\xi }(t)$$

$$\tilde{\rho }={\rho }_0+{\delta }_{\Psi }{\rho }^{(1)}+{\delta }_{\Psi ,\Psi }{\rho }^{(2)}+O({\Psi }^3)$$

$$\dot{X}=F_X(X)+M_1(X)+M_2(X,t)+M_3(X,t)$$

$$M_{1}(X)=\Big{\langle}\Psi_{X}(X,Y)\Big{\rangle}_{\rho_{0,Y}}$$

$$\Big{\langle}\sigma_{R}(X,t)\otimes\sigma_{R}(X,t^{\prime})\Big{\rangle}=\boldsymbol{\rm g}(X,s)=\Big{\langle}\Psi_{X}^{\prime}(X,Y)\otimes\Psi_{X}^{\prime}\big{(}\phi^{s}_{X}(X),\phi^{s}_{Y}(Y))\big{)}\Big{\rangle}_{\rho_{0,Y}}$$

$$M_3(X,t)={\int }_0^{\infty }\mathrm{d}sh(X(t-s),s).$$

$$h(X,s)=\Big{\langle}\Psi_{Y}(X,Y)\cdot\boldsymbol{\rm\nabla}_{Y}\Psi_{X}\big{(}\phi^{s}_{X}(X),\phi^{s}_{Y}(Y)\big{)}\Big{\rangle}_{\rho_{0,Y}}$$

$$-\frac{\partial {\rho }^{\delta }}{\partial s}=\left[\frac{1}{{\delta }^2}{\mathcal{L}}_1+\frac{1}{\delta }{\mathcal{L}}_2+{\mathcal{L}}_3\right]{\rho }^{\delta }$$

$$-\frac{\partial {\rho }^0}{\partial s}=\bar{\mathcal{L}}{\rho }^0$$

$$\dot{X}=F_X(X)+G(X)+\sqrt{2}{\sigma }_{\mathrm{MTV}}(X)\cdot \tilde{\xi }(t)$$

$$G(X)$$

$${\sigma }_{\mathrm{MTV}}(X)$$

$$\dot{X}=\bar{F}(X)=\big{\langle}F(X,Y)\big{\rangle}_{\rho_{Y|X}}$$

$$\dot{Y}=H(X,Y)$$

$$\big{\langle}F(X,Y)\big{\rangle}_{\rho_{Y|X}}=F(X,\bar{Y}(X))+\frac{1}{2}\frac{\partial^{2}F}{\partial Y^{2}}(X,\bar{Y}(X)):\boldsymbol{\rm\Sigma}(X)$$

$$\bar{Y}(X)$$

$$\mathrm{\Sigma }(X)$$

$$\dot{X}=\bar{F}(X)+{\sigma }_{\mathrm{A}}(X)\cdot \xi (t)$$

$$\sigma_{\rm A}(X)=\left(2\int_{0}^{\infty}\,\mathrm{d}s\,\left\langle\big{(}F(X,\phi_{Y|X}^{s}Y)-\bar{F}(X)\big{)}\,\big{(}F(X,Y)-\bar{F}(X)\big{)}\right\rangle_{\rho_{Y|X}}\right)^{1/2}$$

$$\dot{X}=F_X(X)+\mathcal{U}(X)$$

$$\mathcal{U}(X)={\mathcal{U}}_{\det }(X)+e(X(t),t).$$

$$e(X(t),t)=\varphi \frac{{\sigma }_e(X(t))}{{\sigma }_e(X(t-\Delta t))}e(X(t-\Delta t),t-\Delta t)+{\sigma }_e(X(t))(1-{\varphi }^2{)}^{1/2}z(t)$$

$$\left\{\begin{array}[]{lcl}\dot{X}&=&-D\,X+q\,\xi(t)+\frac{\varepsilon}{\delta}Y^{\mathsf{T}}\cdot\boldsymbol{\rm C}\cdot Y\\ \dot{Y}&=&\frac{1}{\delta^{2}}\Big{(}\boldsymbol{\rm A}\cdot Y+\delta\,\boldsymbol{\rm B}_{Y}\cdot\xi_{Y}(t)\Big{)}+\frac{\varepsilon}{\delta}X\,\boldsymbol{\rm V}\cdot Y\end{array}\right.$$

$$Y=\left[\begin{array}[]{c}y_{1}\\ y_{2}\end{array}\right]\quad.$$

$$\boldsymbol{\rm C}=\left[\begin{array}[]{cc}0&B\\ B&0\end{array}\right]\quad,\quad\boldsymbol{\rm A}=\left[\begin{array}[]{cc}-a&\beta\\ -\beta&-a\end{array}\right]\quad,\quad\boldsymbol{\rm V}=\left[\begin{array}[]{cc}0&B_{1}\\ B_{2}&0\end{array}\right]\quad\mathrm{and}\quad\boldsymbol{\rm B}_{Y}=\left[\begin{array}[]{cc}q_{Y}&0\\ 0&q_{Y}\end{array}\right]$$

$$2B+B_1+B_2=0.$$

$$\dot{Y}=\mathrm{T}\cdot Y+\mathrm{B}\cdot {\xi }_Y(t)$$

Full text

I Introduction

From a global point of view, the Earth system is composed of a myriad of different interacting components. These components can be regrouped in compartments like the atmosphere, the hydrosphere, the lithosphere, the biosphere and the cryosphere Olbers (2001)111More recently, a new compartment have appeared, whose effect is not negligible at all and which is not predictable nor descriptive by evolution equations, namely the impact of the human activities.. Those compartments play a role on different timescales from seconds to ice ages. In this perspective, the resulting Earth’s climate is a “concert” at which each compartment seems to play its own partition with its own tempo. Their respective contribution to the total variability of an observable, say e.g. the global temperature, is, however, the outcome of complex interactions between the different components, leading to an emergent dynamics far from the one that could be generated by a linear additive superposition principle Nicolis and Nicolis (1981, 2012).

A paradigmatic example is provided in the work of Hasselmann, detailed in his seminal paper of 1976, which states precisely that the slowly evolving components of the climate system, besides their own dynamics due their own physical processes, also integrate the impact of the faster components Hasselmann (1976). Hasselmann describes this process using the analogy of the Brownian motion where a macro-particle in a liquid integrates the effect of the collisions with the fluid’s micro-particles, leading to the erratic trajectory of the former. The interest of this framework is that it provides a natural description of the “red noise” spectral density observed in most climatic records and observations Ghil et al. (2002); Lovejoy and Schertzer (2013). Subsequently, during the following decade, stochastic modeling for meteorology and climatology became an important research topic Lemke (1977); Frankignoul and Hasselmann (1977); Frankignoul (1979); Frankignoul and Müller (1979); Lemke et al. (1980); Nicolis and Nicolis (1981); Nicolis (1981, 1982); Penland (1989) before falling into disuse in what has been described as a “lull” of work in this field Arnold et al. (2003). However, during that period, the idea that correct parameterizations of sub-grid processes are important to improve climate and weather models gained popularity Penland and Matrosova (1994); Penland (1996); Newman et al. (1997). Stochastic parameterizations for the “turbulent” closure in 2-D large-eddy simulations on the sphere have also been considered Frederiksen and Davies (1997); Frederiksen (1999). It led recently to the implementation of stochastic schemes to correct the model errors Nicolis (2003, 2004) made in large numerical weather prediction (NWP) models Buizza et al. (1999); Shutts (2005), improving the reliability of probabilistic forecasts and correcting partially their variability Nicolis (2005); Doblas-Reyes et al. (2009). The relation between multiplicative noise and the non-Gaussian character observed in some geophysical variables has also been considered Sura et al. (2005); Sardeshmukh and Penland (2015), as well as stochastic models for the climate extremes Sura (2013).

Since the beginning of the 21st century, a revival of the interest in stochastic parameterization methods have occurred, due to the availability of new mathematical methods to perform the stochastic reduction of ODEs systems. Almost simultaneously, a rigorous mathematical framework for the Hasselmann “program” was devised Kifer (2001); Arnold (2001); Kifer (2003) and a new method based on the singular perturbation theory of Markov processes Majda et al. (2001) was proposed. The latter approach is currently known as the Majda-Timofeyev-Vanden-Eijnden (MTV) method. Both methods have been tested and implemented successfully in geophysical models Franzke et al. (2005); Culina et al. (2011); Arnold et al. (2003); Vannitsem (2014). The revival of the Hasselmann program has also stressed the need to consider the occurrence of very rare events triggered by the noise that allow for the solutions of the system to jump from one local attractor to another one Arnold (2001). Such events display recurrence timescales that are few orders greater than the timescale of the climate variables considered, and thus induce transitions between different climates. The statistics of these transitions is then given by the so-called Large Deviations theory Freidlin and Wentzell (1984) (for recent developments on this matter, see Ref. Bouchet et al. (2016)). In addition to these two methods, the modeling of the effects of subgrid scale through conditional Markov chain has been considered Crommelin and Vanden-Eijnden (2008) and recently, new stochastic parameterization techniques have been proposed, based on an expansion of the backward Kolmogorov equation Abramov (2015) and on the Ruelle response theory Wouters and Lucarini (2012). The latter has been tested on a simple coupled ocean-atmosphere model Demaeyer and Vannitsem (2016), on stochastic triads Wouters et al. (2016) and on an adapted version of the Lorenz’96 model Vissio and Lucarini (2016).

This renewal of interest for stochastic modeling and reduction methods illustrates how fruitful was the original idea of Hasselmann. However, in view of the availability of several possible approaches, one might wonder about their efficiency in different situations. Indeed, depending on the specific purpose that it needs to fulfill, some parameterizations might perform better than others. The present review aims to shed some light on these questions and to illustrate some of the aforementioned parameterization methods on a simple model for which most of the calculations can be made analytically.

In Section II, we will present the general framework in which the problem of stochastic parameterizations is posed. In Section III, we present the different parameterizations that we shall consider for the analysis model. The stochastic triad model used here and the comparison are presented in Section IV. Finally, the conclusions are given in Section. V.

II The parameterization problem

Consider the following ordinary differential system of equations (ODEs):

[TABLE]

where $Z\in\mathbb{R}^{d}$ is a set of variables relevant for the problem under interest for which the tendencies $T(Z)$ are known. And suppose that one wants to separate this set of variables into two different subset $Z=(X,Y)$, with $X\in\mathbb{R}^{m}$ and $Y\in\mathbb{R}^{n}$. In general, such a decomposition is made such that the subset $X$ and $Y$ have strongly differing response times $\tau_{Y}\ll\tau_{X}$ Arnold et al. (2003), but we will assume here that this constraint is not necessarily met. System (1) can then be expressed as:

[TABLE]

The timescale of the $X$ sub-system is typically (but not always) longer than the one of the $Y$ sub-system, and it is often materialized by a parameter $\delta=\tau_{Y}/\tau_{X}\ll 1$ in front of the time derivative $\dot{Y}$. The $X$ and the $Y$ variables represent respectively the resolved and the unresolved sub-systems. The resolved sub-system is the part of the full system that we would like to simulate, i.e. generate explicitly and numerically its time-evolution. The general problem of model reduction consists thus to approximate the resolved component $X$ as accurately as possible by obtaining a closed equation for the system $X$ alone 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 $X$ component of system (2). We can also ask that the closed system trajectories remains as close as possible to the trajectories of the full system for long times.

In general a parameterization of the sub-system $Y$ is a relation $\Xi$ between the two sub-systems:

[TABLE]

which allows to effectively close the equations for the sub-system $X$ while retaining the effect of the coupling to the $Y$ sub-system.

The problem of the model reduction is not new, and was considered first in celestial mechanics. Famous mathematician have considered it and contributed to what is known nowadays as the theory of averaging Sanders and Verhulst (1985) and which forms the first step of the Hasselmann program Arnold (2001). The mathematical framework was set in the 1960s by the influential contribution of Bogoliubov and Mitropolski Bogoliubov and Mitropolski (1961). However, this averaging technique is a deterministic method which does not take into account the deviations from the average. The proposition of Hasselmann was thus to take into account these deviations by considering stochastic parameterization where the relation (3) can be considered in a statistical sense. In that framework, the $Y$ sub-system and its effect on the sub-system $X$ can be considered as a stochastic process, which possibly depends upon the state of the $X$ sub-system. Different methods to achieve this program are now discussed in Section III.

III The parameterization methods

Let us now write system (2) as:

[TABLE]

where the coupling and the intrinsic dynamics are explicitly specified. In the present work, following the Hasselmann program, we shall focus on parameterizations that are defined in terms of stochastic processes. We will consider methods based

• •

on the Ruelle response theory Wouters and Lucarini (2012); Demaeyer and Vannitsem (2016); Wouters et al. (2016).

• •

on the singular perturbation theory of Markov processes Majda et al. (2001); Franzke et al. (2005).

• •

on the Hasselmann averaging methods Kifer (2003); Arnold et al. (2003); Culina et al. (2011); Vannitsem (2014).

• •

on empirical methods Arnold et al. (2013).

All these parameterizations can be written in the following form:

[TABLE]

where the matrix $\boldsymbol{\rm\sigma}$, the deterministic function $G$ and the random processes $\tilde{\xi}(t)$ have to be determined. The mathematical definition of these quantities obtained through averaging procedure and the measure being used to perform the averaging are usually both differing between the methods. These different choices are rooted in their different underlying hypothesis, as it will be discussed below. Specifically, the response theory method uses the measure of the uncoupled unresolved sub-system $\dot{Y}=F_{Y}(Y)$, the singular perturbation method uses the measure of the perturbation, and the averaging methods use the measure of the full unresolved sub-system $\dot{Y}=H(X,Y)$ with $X$ considered as “frozen” (constant). Finally, the empirical methods use in general the output of the full unresolved $Y$ sub-system, conditional or unconditional on the state $X$.

In the rest of the section, we shall describe more precisely each of the above methods.

III.1 The method based on response theory

This method is based on the Ruelle response theory Ruelle (1997, 2009) and was proposed by Wouters and Lucarini Wouters and Lucarini (2012, 2013). In this context, system (4) must be considered as two intrinsic sub-dynamics for $X$ and $Y$ that are weakly coupled. The response theory quantifies the contribution of the “perturbation” $\Psi_{X}$, $\Psi_{Y}$ to the invariant measure222The theory assumes that for the system under consideration, a SRB measure Young (2002) exists (e.g. an Axiom-A system). $\tilde{\rho}$ of the fully coupled system (4) as:

[TABLE]

where $\rho_{0}$ is the invariant measure of the uncoupled system which is also supposed to be an existing, well defined SRB measure. As shown in Wouters and Lucarini (2012), this theory gives the framework to parameterize the effect of the coupling on the component $X$. The parameterization is based on three different terms having a response similar, up to order two, to the couplings $\Psi_{X}$ and $\Psi_{Y}$:

[TABLE]

where

[TABLE]

is an averaging term. $\rho_{0,Y}$ is the measure of the uncoupled system $\dot{Y}=F_{Y}(Y)$. The term $M_{2}(X,t)=\sigma_{R}(X,t)$ is a correlation term:

[TABLE]

where $\otimes$ is the outer product, $\Psi_{X}^{\prime}(X,Y)=\Psi_{X}(X,Y)-M_{1}(X)$ is the centered perturbation and $\phi^{s}_{X}$, $\phi^{s}_{Y}$ are the flow of the uncoupled system $\dot{X}=F_{X}(X)$ and $\dot{Y}=F_{Y}(Y)$. The process $\sigma_{R}$ is thus obtained by taking the square root of the matrix $\boldsymbol{\rm g}$, which is here accomplished by a decomposition of $\Psi_{X}^{\prime}$ on a proper basis Wouters and Lucarini (2012). The $M_{3}$ term is a memory term:

[TABLE]

involving the memory kernel

[TABLE]

All the averages are thus taken with $\rho_{0,Y}$, the invariant measure of the unperturbed system $\dot{Y}=F_{Y}(Y)$. This particular choice of the measure is due to the perturbative nature of the method and simplify the averaging procedure in many case. 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 perturbation $\Psi_{X}$ and $\Psi_{Y}$. Consequently, this ensures that for a weak coupling, the response of the parameterization (7) on the observables will be approximately the same as the coupling.

The advantages of this simplified averaging procedure (by using $\rho_{0,Y}$) should be tempered by the additional cost induced by the computation of the memory term, the latter implying that this parameterization is a non-Markovian one Chekroun et al. (2015). However, the integral (10) in this memory term must only be evaluated from [math] up to the timescale $\tau_{Y}$ of the fast variable, due to the exponential decrease of the integrand. Moreover, in some cases, this non-Markovian parameterization can be effectively replaced by a Markovian one Wouters et al. (2016).

III.2 Singular perturbation theory method

Singular perturbation methods were developed in the 1970s for the analysis of the linear Boltzmann equation in some asymptotic limit Grad (1969); Ellis and Pinsky (1975); Papanicolaou (1976); Majda et al. (2001). Here, these methods are applicable if the problem can be cast into a Fokker-Planck equation. The procedure described in Ref. Majda et al. (2001) requires assumptions on the timescales of the different terms of system (4). In term of the small parameter $\delta=\tau_{Y}/\tau_{X}$ defined in Section II, the fast variability of the unresolved component $Y$ are considered of order $O(\delta^{-2})$ and modeled as an Ornstein-Uhlenbeck process. The Markovian nature of the process defined by Eq. (4) and its singular behavior in the limit of an infinite timescale separation ($\delta\to 0$) allow then to apply the method.

More specifically, the parameter $\delta$ serves to distinguish terms with different timescales and is then used as a small perturbation parameter. In this setting, the backward Fokker-Planck equation reads Majda et al. (2001):

[TABLE]

where the function $\rho^{\delta}(s,X,Y|t)$ is defined with the final value problem $f(X)$: $\rho^{\delta}(t,X,Y|t)=f(X)$. The function $\rho^{\delta}$ can be expanded in term of $\delta$ and inserted in Eq. (12). The zeroth order of this equation $\rho^{0}$ can be shown to be independent of $Y$ and its evolution given by a closed, averaged backward Fokker-Planck equation Kurtz (1973):

[TABLE]

This equation is obtained in the limit $\delta\to 0$ and gives the sought limiting, averaged process $X(t)$. Note that this procedure does not necessarily require the presence of the explicit small parameter $\delta$ in the original equation (4). Since $\delta$ disappears from Eq. (13), one can simply use the parameter to identify the fast terms to be considered, and eventually consider $\delta=1$ Franzke et al. (2005).

The parameterization obtained by this procedure is given by Franzke et al. (2005):

[TABLE]

with

[TABLE]

with the same notation as in the previous subsection. The measure $\tilde{\rho}$ is the measure of the $O(\delta^{-2})$ perturbation, i.e. the source of the fast variability of the unresolved $Y$ component. This measure thus depends on which terms of the unresolved component are considered as “fast”, and some assumptions should here be made. For instance, it is customary to consider as the fast terms the quadratic terms in $Y$ and replaced them by Ornstein-Uhlenbeck processes whose measures is used to compute the averages Majda et al. (2001); Franzke et al. (2005).

Finally, if one assume that the source of the fast variability in the sub-system is given by the “intrinsic” term $F_{Y}(Y)$ (such that $\tilde{\rho}=\rho_{0,Y}$) and if the perturbation $\Psi_{X}$ only depends on $Y$, this parameterization is simply given by the integration of the function $g(s)$ and $h(X,s)$ of the response theory parameterization given by Eqs. (9) and (11). This can be interpreted as an averaging of the latter parameterization when the timescale separation is infinite and $X$ can thus be considered as constant over the timescale of the integrand. Therefore, $M_{2}$ can be modeled as a white noise and the memory term is Markovian.

III.3 Hasselmann averaging method

Since the initial work of Hasselmann in the 1970s Hasselmann (1976), various approaches have been considered to average directly the effects of the “fast” evolving variables on the “slow” ones. These methods assume in general a sufficient timescale separation between the resolved and unresolved components of the systems, and a direct average can be performed as,

[TABLE]

where $\rho_{Y|X}$ is the measure of the system

[TABLE]

conditional on the value of $X$. In this approach, $X$ is thus viewed as a constant parameter for the unresolved dynamics. In other words, this particular framework assumes that since $X$ is slowly evolving with respect to the typical timescale of $Y$, it can be considered as “frozen” while $Y$ evolves. With some rigorous assumptions, this approach has been mathematically justified Kifer (2003) and applied successfully to idealized geophysical models Arnold et al. (2003) with non-trivial invariant measures. In the same vein, an approximation has been proposed in Ref. Abramov (2013) for the average (17), assuming that $F$ is at most quadratic,

[TABLE]

where “:” means the element-wise matrix product with summation and where

[TABLE]

The approximation to the second order is particularly well suited for the application to atmospheric and climate flows for which the quadratic terms are usually the main non-linearities associated with the advection in the system.

In Ref. Abramov (2013), an approach based on the fluctuation-dissipation theorem was proposed to estimate the mean state $\bar{Y}(X)$ and the covariance matrix $\boldsymbol{\rm\Sigma}(X)$.

The deterministic parameterization (17) can be recast in a stochastic parameterization following the same principle. Such a parameterization has been derived in Refs. Arnold et al. (2003); Abramov (2015) and reads

[TABLE]

with

[TABLE]

where $\phi_{Y|X}^{s}$ is the flow of the system (18) for $X$ constant (“frozen”). A drawback of such an approach is that it requires that the measure $\rho_{Y|X}$ exists and is well-defined (ideally a SRB measure Arnold et al. (2003)). Such a requirement may not be always fulfilled, like for instance if the fast system conditional on the state $X$ is unstable and does not possess any attractor (see Section IV for an example).

III.4 Empirical methods

The empirical methods are generally based on the statistical analysis of the timeseries $Y$ of the full system (4). Many procedures exist as discussed in Section I but we will consider here a method based on state dependent AR(1) processes proposed in Ref. Arnold et al. (2013). In this case, a timeserie $r(t)$ of the coupling part $\Psi_{X}$ of the $X$ tendency must first be computed with (4). The parameterization is then given by

[TABLE]

with

[TABLE]

The function $\mathcal{U}_{\det}(X)$ represents the deterministic part of the parameterization and is obtained by a least-squares fit of the timeserie $r(t)$ versus the timeserie $X(t)$ with the cubic function $\mathcal{U}_{\det}(X)=p_{0}+p_{1}\,X+p_{2}\,X^{2}+p_{3}\,X^{3}$. The “stochastic” part $e(X(t),t)$ is then given by the following state-dependent AR(1) process

[TABLE]

where $z(t)$ is a standard Gaussian white noise process. The parameters of the process $e$ are determined by considering the residual timeserie $r(t)-\mathcal{U}_{\det}(X(t))$ to compute the lag-1 autocorrelation $\phi$ and the state dependent standard deviation $\sigma_{e}(X)$ which is modeled as $\sigma_{e}(X)=\sigma_{0}+\sigma_{1}\,|X|$ with the parameters $\sigma_{0}$ and $\sigma_{1}$ given by a binning procedure. The parameter $\Delta t$ is the time step of integration of Eq. (24). Other empirical parameterizations have been proposed by the authors of Ref. Arnold et al. (2013), notably one with the function $\mathcal{U}(X)=(1+e(t))\,\mathcal{U}_{\det}(X)$ which resembles the SPPT333Acronym for Stochastically Perturbed Parameterization Tendencies Scheme. parameterization used in the ECMWF444Acronym for European Center for Medium-Range Weather Forecasts. Numerical Weather Prediction model Buizza et al. (1999). However, the study shows no substantial differences with the parameterization (25).

IV Applications and results

In this section, we will illustrate the various parameterizations described in Section III to the following example

[TABLE]

where $D>0$, $q>0$ and

[TABLE]

with $X$, $y_{1}$, $y_{2}\in\mathbb{R}$. The matrices involved are defined as

[TABLE]

with $a,\beta,q_{Y}>0$. The process $\xi(t)$ and $\xi_{Y}(t)$ are uncorrelated standard Gaussian white noise processes.

The $X$ and $Y$ variables represent respectively the resolved and the unresolved sub-systems. The parameter $\delta>0$ quantify the time-scale separation of the terms of the tendencies of the two components, with the three timescales $O(1)$,$O(\delta^{-1})$ and $O(\delta^{-2})$ as in Ref. Majda et al. (2001) (see also Subsection III.2). Additionally, the parameter $\varepsilon>0$ controls the coupling strength between the two sub-systems. In this setup the coupling are thus proportional to the ratio $\varepsilon/\delta$, and therefore the characterization of the coupling as “weak” depends directly on the time-scale separation.

The deterministic part of Eq. (27) is a well known model called a triad and encountered in fluid dynamics Waleffe (1992); Ohkitani and Kida (1992); Smith and Waleffe (1999), and in simplified geophysical flows e.g. Majda et al. (2001); Wouters et al. (2016). Due to the presence of invariant manifolds, its mathematical structure can be found in higher-order model. See Ref. Demaeyer and Vannitsem (2016) for an example of such structure in the framework of a coupled ocean-atmosphere model. In the present study, the interest of the stochastic triad model (27) is that, $H(X,Y)$ being linear in $Y$, the measure $\rho_{0,Y}$ and $\rho_{Y|X}$ can be analytically computed since both $\dot{Y}=F_{Y}(Y)$ and $\dot{Y}=H(z,Y)|_{z=X}$ are two-dimensional Ornstein-Uhlenbeck processes. Therefore, for this simple case, the set of methods proposed in the previous section can be applied exactly without resorting to a binning procedure of the output of the $Y$ sub-system555Except for the empirical methods which by definition use this kind of procedures..

As energy conservation is a rule in physical systems in the absence of dissipation and fluctuations, we will adopt this rule for the current system. With the dissipation and stochastic terms set to zero, the system (27) conserves the “energy” $(X^{2}+y_{1}^{2}+y_{2}^{2})/2$ if the coefficient $B$, $B_{1}$ and $B_{2}$ are chosen such that Smith and Waleffe (1999); Majda et al. (2001)

[TABLE]

It allows for the following configurations of their signs : $(+,-,-)$, $(+,+,-)$, $(+,-,+)$, $(-,+,+)$, $(-,+,-)$, $(-,-,+)$. These different configurations are associated with different kinds of energy exchange scenario and different stability properties Waleffe (1992).

We will focus on the two configurations $(-,-,+)$ and $(-,+,+)$, with parameters

$B=-0.0375$, $B_{1}=-0.025$, $B_{2}=0.1$ 2. 2.

$B=-0.0375$, $B_{1}=0.025$, $B_{2}=0.05$

and consider various values of the parameters $\delta$ and $\varepsilon$. The other parameters are fixed to $a=0.01$, $D=0.01$ and $\beta=0.01/12$. Once the parameterizations have been developed, the different model version have been integrated over $4.5\times 10^{5}$ timeunits with a timestep $\Delta t=0.01$ after a transient period of $5.0\times 10^{4}$ timeunits to let the system relax to its stationary state. The state $X$ has been recorded every $0.1$ timeunit, giving a dataset of $4.5\times 10^{6}$ points for the analysis. The parameterizations given by Eqs. (7), (14) and (22) have been integrated with a second order Runge-Kutta (RK2) stochastic scheme which converges to the Stratonovich calculus Hansen and Penland (2006); Rüemelin (1982). The equation (24) has been integrated with a deterministic RK2 scheme where the stochastic forcing $e(X,t)$ is considered constant during the timestep. The memory term $M_{3}$ appearing in the parameterization (7) and given by the integral (10) over the past of $X$ has been computed numerically at each timestep. Although it increases considerably the integration time, this method is adopted in order to clarify the memory effect in Eq. (7). A Markovianization of this parameterization is possible Wouters et al. (2016) but in the present case it would have required some assumptions that would blur the comparison of the methods.

The relative performances of the parameterizations can be tested in multiple ways, by comparing the climatology (the average state) or the variability (variance) of the systems Nicolis (2005). Another method is to look at the predictive skill score of the models, that is the ability of the parameterizations to provide skillful forecast compared to original system, as in Arnold et al. (2013); Wouters et al. (2016). On longer term, the good representation of the “climate” of a model by the parameterizations can be assessed by looking at the stationary probability densities and comparing them using some score Franzke et al. (2005); Crommelin and Vanden-Eijnden (2008); Abramov (2012, 2013, 2015). The decorrelation properties of the models and the parameterizations can also be tested, to provide information about the correct representation of the timescales of the models. All those different aspects can be significant, depending on the purpose of the parameterization scheme. However, for the brevity of the present work, we shall focus on the probability densities and whether or not they are correctly reproduced by the parameterizations.

We present now the results obtained the proposed methods and consider first the different measures used for averaging in system (27).

IV.1 Stability and measures

All the ingredients needed to compute the parameterizations presented in Section III can be derived with the help of the covariance and the correlation of the $Y$ variables in the framework of two different systems related to the unresolved dynamics, namely the unperturbed dynamics $\dot{Y}=F_{Y}(Y)$ and the unresolved dynamics $\dot{Y}=H(X,Y)$ with $X$ frozen. The measure of the former is necessary to derive the response theory and the singular perturbation-based parameterizations, while the latter is needed for the Hasselmann averaging method. These two systems are both two-dimensional Ornstein-Uhlenbeck processes of the form

[TABLE]

for which respectively $\boldsymbol{\rm T}=\boldsymbol{\rm A}/\delta^{2}$ and $\boldsymbol{\rm T}=\boldsymbol{\rm A}/\delta^{2}+\left(\varepsilon\,X/\delta\right)\,\boldsymbol{\rm V}\cdot Y$. In both cases, we have $\boldsymbol{\rm B}=\boldsymbol{\rm B}_{Y}/\delta$. Their measure is then given by Wouters et al. (2016)

[TABLE]

where $\mathcal{Z}$ is a normalization factor and where $\boldsymbol{\rm\Sigma}$ is the covariance matrix solution of

[TABLE]

In order for theses processes to be stable, the real part of the eigenvalues of the matrix $\boldsymbol{\rm T}$ must be negative Gardiner (2009) for every state $X$ that the full coupled system (27) can possibly achieve. The eigenvalues of the system $\dot{Y}=F_{Y}(Y)$ are $\lambda_{\pm}=(-a\pm{i\mkern 1.0mu}\beta)/\delta^{2}$ and it is thus always stable (since $a>0$ and $\beta\in\mathbb{R}$). On the other hand, the system $\dot{Y}=H(X,Y)$ has the eigenvalues $\lambda_{\pm}=(-a\pm\sqrt{\Delta(X))}/\delta^{2}$ with

[TABLE]

Therefore, if

[TABLE]

for some $X$, the Ornstein-Uhlenbeck process is unstable, and it is then called an explosive process. For any initial condition, the process diverges, and thus the only possible stationary measure is the trivial null one. Consequently, equation (33) gives nonphysical solutions, the stationary covariance matrix does not exist and the parameterizations depending upon cannot be derived.

For the system (27), if $\mathrm{sgn}(B_{1}\,B_{2})=-1$, as in case 1, then the process is stable for every $X$ if $a^{2}>-\frac{\left(B_{1}+B_{2}\right){}^{2}\beta^{2}}{4B_{1}B_{2}}$. For case 1, this inequality is satisfied, and thus the process (31) is stable for every $X$. Moreover, depending on the sign of $\Delta(X)$, the process for $X$ fixed is a stochastic focus (if $\Delta(X)>0$) or a stochastic damped oscillator (if $\Delta(X)<0$). Here, it is a focus if

[TABLE]

and an oscillator otherwise. That is, for the considered $\varepsilon$ and $\delta$ parameters value, the system (27) is an oscillator for most of the $X$ values.

If $\mathrm{sgn}(B_{1}\,B_{2})=1$, as in case 2, then the condition (35) must be satisfied for every state $X$. For case 2, this inequality was not satisfied for every state $X$ for most of the values of the $\varepsilon$ and $\delta$ parameters considered (see Subsection IV.3 below). The stability is therefore reversed as the system is non-oscillating for most of the $X$ values.

To summarize, if $B_{1}$ and $B_{2}$ are of opposite sign, the dynamics of $\dot{Y}=H(X,Y)$ is stable and generally oscillatory. If $B_{1}$ and $B_{2}$ have the same sign, then the dynamics is unstable in most cases and generally hyperbolic. This is a consequence of the well known difference of stability of the triads depending on their energy exchange properties Waleffe (1992).

For the interested reader, the exact calculation of the parameterization of Section III using the covariance and correlation matrices are detailed in the appendix (see Section A).

IV.2 The $(-,-,+)$ stochastic triad (case 1)

Let us now consider case 1 corresponding to the $(-,-,+)$ stochastic triad for two different values of the timescale separation $\delta=0.1$ and $0.4$. For each of these timescale separation, we considered three values of the coupling strength $\varepsilon$: $0.05$, $0.125$ and $0.4$. The probability densities associated with these different systems are represented on Figs. 1 and 3. For a timescale separation $\delta=0.1$, the fully coupled dynamics given by Eq. (27) is quite well represented by all the proposed parameterizations. Since it is hard to distinguish the different density curves, a score such as the Hellinger distance Arnold et al. (2013)

[TABLE]

between the distribution $P$ of the full coupled system and the distribution $Q$ of the parameterizations is worth computing to quantify the differences. This distance quantifies the similarity between the densities. It is depicted for $\delta=0.1$ on Fig. (2), and it shows that for a very small coupling parameter $\varepsilon=0.05$, the best parameterization is the response theory given by Eq. (7). For larger values of $\varepsilon$, it is the Hasselmann averaging method which performs best. The empiric method gives a good correction of the uncoupled dynamics for $\varepsilon=0.125$ but diverges for $\varepsilon=0.4$. This may be due to instabilities introduced by the cubic deterministic parameterization $\mathcal{U}_{\det}(X)$ or to the inadequacy of the fitting function $\sigma_{0}+\sigma_{1}\,|X|$ for the standard deviation $\sigma_{e}(X)$ in the AR(1) process (26). Indeed, in general, this model fits quite well the statistics in the neighborhood of $X=0$, but the standard deviation reaches a plateau for higher values of $X$. A more complicated fitting function would thus be necessary to get a stable dynamics. For an decreased timescale separation $\delta=0.4$, the same conclusions are reached as for $\delta=0.1$, but the singular perturbation method performs not very well in all cases, as illustrated in Fig. 4 that for $\varepsilon=0.125$ and $0.4$. The response based and singular perturbation methods are even less effective than the uncoupled dynamics. It is not surprising for the latter since it is supposed to be valid in the limit $\delta\to 0$.

IV.3 The $(-,+,+)$ stochastic triad (case 2)

We now consider the parameters of case 2, for which the system (27) is a $(-,+,+)$ stochastic triad. The probability densities are depicted on Figs. 5 and 7, and the Hellinger distances are shown on Figs. 6 and 8. First, we must remark that the parameterization based on the Hasselmann’s averaging method is not defined for most of the $\delta$ and $\varepsilon$ parameters values. It is due to the fact that the dynamics of the unresolved component $Y$ with $X$ considered as a parameter is unstable, as shown in Subsection IV.1. Indeed, this linear system undergoes a bifurcation at some value $X^{\star}$ which destabilizes the dynamics $\dot{Y}=H(X,Y)$ with $X$ frozen. Therefore, the measure $\rho_{Y|X}$ is not defined for some ranges of the full $X$ dynamics and the method fails. The only case where this destabilization does not occur is for $\delta=0.1$ and $\varepsilon=0.05$, but the parameterization does not perform well. For these parameter values, the only parameterization that performs very well is the one based on response theory. For the other values of the parameters $\delta$ and $\varepsilon$, all the parameterizations have good performances. A particularly unexpected result is the very good correction provided by the response theory and singular perturbation based methods for the extreme case $\delta=0.4$ and $\varepsilon=0.4$ (see the bottom panel of Fig. 8). This have to be contrasted with their bad performances in the case of the other triad (see the bottom panel of Fig. 4). Note that for this extreme case, the direct averaging method fails and the empirical method is unstable and diverges.

IV.4 Discussion

The results obtained so far with these two types of triads highlight the utility of the parameterization schemes discussed here. First, the empiric parameterization gives usually good results when it does not destabilize the dynamics. However, this method requires a case by case time-consuming statistical analysis whose complexity increases with the dimensionality of the problem considered. Physically based parameterizations do not require such an analysis, and the best approach in the present system is the Hasselmann averaging one, but it requires that the dynamics of the unresolved system be stable. It was thus very effective to correct the dynamics of the $(-,-,+)$ triad, but not the other triad $(-,+,+)$. In this latter case, the perturbative methods like the singular perturbation method or the response theory method gives very good results. This difference is quite intriguing and interesting. It indicates that different physically based parameterization should be considered depending on the kind of problems encountered. In particular, the stability properties of the system considered seems to play an important role. This conclusion holds whatever the timescale separation and for the most realistic values of the coupling strength between the components ($\varepsilon=0.125$ and $0.4$). However, for very small values of the coupling strength, the response based method seems to be the best approach in all cases.

A question that is left open in the present work is to determine precisely which stability property is giving the contrasting observed result. In extenso, is it the hyperbolic instability of the $(-,+,+)$ triad which makes the perturbative approach and the response based parameterization perform so well, or is it simply the fact that it is unstable? On the other hand, is it the damped oscillatory behavior of the $(-,-,+)$ triad which makes the Hasselmann’s method works well, or is it simply the fact that it is stable? Such questions should be addressed in the case of a more complex, globally stable system, which allows to have locally stable and unstable fast dynamics.

V Conclusions

The parameterization of subgrid-scale processes is an important tool available for meteorologist and climatologist, for stochastic modeling, model reduction, or to improve the statistical properties of the forecasting systems. The variety of different approaches available bear witness of the richness of the field but at the same time can also lead to questions on the best choice for the problem at hand. The purpose of the present review was to describe briefly some of the most recent methods and to illustrate them on a simple stochastic triad example. The methods covered included perturbative methods like the Ruelle response theory Wouters and Lucarini (2012) or the singular perturbation theory Majda et al. (2001), averaging methods like the Hasselmann method Hasselmann (1976); Arnold et al. (2003) or empirical methods Arnold et al. (2013). As expected, these parameterizations provided good results depending on the timescale separation and on the coupling between the resolved variables and the subgrid one. But more importantly, our results in the context of this simple triad stress the importance of the underlying stability properties of the unresolved system. It thus confirms a known result that the structure of the Jacobian and of the Hessian of a given system control the behavior and performance of model error parameterizations Nicolis (2005).

Further comparisons of the different methods are needed in the context of more sophisticated systems in order to analyze the role of the stability properties of the sub-grid scale processes on their performances. This type of analysis is currently under way in the context of a coupled ocean-atmosphere system De Cruz et al. (2016).

Appendix A Appendix : Practical computation of the parameterizations

In the following section, for illustrative purposes, we detail the computation that we have made to obtain the result of the present review. We start with the method based on response theory.

A.1 Response theory method

We consider the system (27) with the form (4) in mind. In this case, the influence of the $Y$ sub-system on the $X$ sub-system is parameterized as:

[TABLE]

where then terms $M_{1}$,$M_{2}$ and $M_{3}$ are respectively given by Eqs. (8), (9 and (10). The average in these formula are performed with the measure $\rho_{0,Y}$ of the unperturbed $Y$ dynamics $\dot{Y}=F_{Y}(Y)$. Since this latter is an Ornstein-Uhlenbeck process, its measure is the Wiener measure

[TABLE]

where $\boldsymbol{\rm\Sigma}$ is the covariance matrix solution of

[TABLE]

and $\mathcal{Z}$ is a normalization factor.

The covariance and correlation of the stationary process $\dot{Y}=F_{Y}(Y)$ are thus straightforward to compute Gardiner (2009):

[TABLE]

where $\boldsymbol{\rm I}$ is the identity matrix, $\phi_{Y}^{t}$ is the flow of $\dot{Y}=F_{Y}(Y)$ and the matrix $\boldsymbol{\rm E}(t)$ is the exponential

[TABLE]

The various terms $M_{i}$ are then computed as follow.

A.1.1 The term $M_{1}$

It is the average term:

[TABLE]

We have thus

[TABLE]

by using Eq. (41).

A.1.2 The term $M_{2}$

It is the noise/correlation term which is defined here as:

[TABLE]

with

[TABLE]

and the correlation function

[TABLE]

where $\Psi_{X}^{\prime}(Y)=\Psi_{X}(Y)-M_{1}$. The result in the present case is given by the formula (see Ref. Demaeyer and Vannitsem (2016)) :

[TABLE]

The term $M_{2}$ must thus be devised as a process with the same correlation.

A.1.3 The term $M_{3}$

This is the memory term, defined by

[TABLE]

with the memory kernel

[TABLE]

which in the present case is given by the formula (see Ref. Demaeyer and Vannitsem (2016))

[TABLE]

The fact that the memory kernel (54) and the correlation function (50) present the same form imply that a Markovian parameterization is available Wouters et al. (2016) even if by definition, Eq. (38) is a non-Markovian parameterization.

A.2 The singular perturbation method

With this parameterization, the parameter $\delta$ serves to distinguish terms with different timescale and is then used as a small perturbation parameter Majda et al. (2001); Franzke et al. (2005). The parameterization is given by:

[TABLE]

with notably $\langle\xi(t)\tilde{\xi}(t^{\prime})\rangle=0$ and

[TABLE]

We see that the quantities appearing in this parameterization can easily be obtained from the functions $h$ and $g$ of the section A.1. Indeed we have,

[TABLE]

where we notice that the parameter $\delta$ has disappeared, since this parameterization is valid in the limit $\delta\to 0$.

A.3 Averaging method

In this approach, we consider the system (2) and the parameterization from Ref. Abramov (2013):

[TABLE]

with

[TABLE]

and

[TABLE]

where $\rho_{Y|X}$ is the measure of the system $\dot{Y}=H(X,Y)$ with $X$ “frozen”. It is the measure of an Ornstein-Uhlenbeck process

[TABLE]

where $\mathcal{Z}$ is a normalization factor and $\boldsymbol{\rm\Sigma}(X)$ is the stationary covariance matrix solution of

[TABLE]

with

[TABLE]

With the help of $\bar{Y}(X)=0$ and $\boldsymbol{\rm\Sigma}(X)$, we can now rewrite Eq. (61) as

[TABLE]

This forms a deterministic averaging parameterization. It can be extended into a stochastic parameterization Abramov (2015) as follow:

[TABLE]

with

[TABLE]

and

[TABLE]

We thus have

[TABLE]

where we have extended the result of Eq. (50) to the stationary Ornstein-Uhlenbeck process $\dot{Y}=H(X,Y)$ for X “frozen”. The function $S(X)$ can be computed analytically using mathematical softwares but is a very complicate function that is not worth displaying in this review. This can however be provided upon query to Dr. J. Demaeyer.