Parametrization of stochastic multiscale triads

We discuss applications of a recently developed method for model reduction based on linear response theory of weakly coupled dynamical systems. We apply the weak coupling method to simple stochastic differential equations with slow and fast degrees of freedom. The weak coupling model reduction method results in general in a non-Markovian system, we therefore discuss the Markovianization of the system to allow for straightforward numerical integration. We compare the applied method to the equations obtained through homogenization in the limit of large time scale separation between slow and fast degrees of freedom. We look at ensemble spread from a fixed initial condition, correlation functions and exit times from domain. The weak coupling method gives better results in all test cases, albeit with a higher numerical cost.


Introduction
Many models of physical systems are too complex to be solved analytically, or even numerically if a large range of temporal and spatial scales is involved. For some high-dimensional dynamical systems it is however possible to derive lower-dimensional reduced models (Givon et al., 2004;Huisinga et al., 2003). The reduced model is easier to solve analytically and faster to integrate numerically, while still preserving some of the essential characteristics of the full system. This line of research lies at the heart of many applications, for example in molecular dynamics (Hijón et al., 2009;Lu and Vanden-Eijnden, 2014) and climate modeling (Lucarini et al., 2014;Imkeller and Von Storch, 2001;Palmer and Williams, 2009).
The derivation of a reduced model is possible, for example, in the presence of a time scale separation between slow resolved and fast unresolved variables, as is assumed in the homogenization method (Pavliotis and Stuart, 2008). This method applies to slow-fast systems of the forṁ in the limit of infinite time scale separation ε → 0, where ξ denotes a standard Brownian motion (i.e. the equations should be considered equivalent to a stochastic integral in the Itô interpretation) (Khas'minskii, 1963;Papanicolaou, 1976). It is evident from the dynamical equation that the y variables evolve on a faster time scale than the x variables. For finite values of ε there is an intricate feedback between the evolution of the x and y variables. The situation simplifies in the limit of ε → 0 where the slow variables do not evolve on the time scales on which y strongly fluctuates. As a result, the slow dynamics converges to a stochastic evolution, where the effect of y is completely replaced by statistical quantities related to the motion of y for a fixed value of x. On a more technical note, the precise expression for the quantities entering in the reduced dynamics can be easily obtained through an expansion in ε of the backward Kolmogorov equation (the adjoint of the Fokker-Planck equation) ∂ t v(x, t) = (L 0 + L 1 /ε + L 2 /ε 2 )v(x, t) of corresponding to the slow-fast dynamics (where L 0 = f 0 ∂ x , L 1 = f 1 ∂ x and L 2 = g 1 ∂ y + (β/2)∂ 2 y ) (Pavliotis and Stuart, 2008). The method of homogenization has found a great number of applications in different fields of physics and mathematics (Pavliotis and Stuart, 2008). Many physical systems, however, do not feature a time scale separation. As an example, the climate system has variability on many different temporal (and spatial) scales, but no clear spectral gaps can be identified. This creates fundamental difficulties in the theoretical investigation of climate dynamics and in the construction of climate models.
As a result, approximate equations are used for dealing with scales of motions belonging to a range of scales of interest, and numerical models are able to resolve explicitly only a fractions of the full range of scales. The dynamics taking place on scales that are too small and/or fast to be resolved need to be parametrized. Consider the case of convective motion in the Earth's atmosphere. Convective clouds are significant for the climate, yet can only be resolved at a spatial resolution of 10-100 m (Sakradzija et al., 2015), whereas climate models only resolve scales of the order of 100 km (Intergovernmental Panel on Climate Change, 2013). Unresolved convective motion however features a so-called "gray zone", a range of time scales overlapping with the dynamical time scales of the resolved large scale flow (Sakradzija et al., 2015), therefore homogenization can not be applied. It is a formidable challenge to derive dimension reduction methods that do not require a time scale separation.
One should underline that when facing a lack of time scale separation, we would like to be able to construct self-adaptive parametrizations as opposed to empirical ones, so that when the resolution of a numerical model is changed we do not need to redo the exercise of fitting a reduced model.
Going beyond the familiar setting of infinite time scale separation requires a novel approach to the derivation of closed equation for the reduced system. Recently, we have developed a model reduction technique that does not rely on the presence of such a separation (Lucarini et al., 2014;Wouters and Lucarini, 2012;Wouters and Lucarini, 2013). The alternative method for model reduction makes use of a weak coupling approach, in which response theory (Ruelle, 2009(Ruelle, , 1997 is used to derive a closure. The systems of interest follow a dynamics determined bẏ where x is the variable of interest. Exploiting the weak coupling form of this equation, response theory can be employed to expand expectation values of x-dependent observable under the invariant measure in orders of ε. This expansion yields a series in terms of ε, reminiscent of the Dyson series in scattering theory, each representing a sequence of interactions between the x and y subsystems, corresponding to a certain Feynman diagram.
The truncation of this series up to a given order yields an approximation of the response of the x subsystem to the coupling to the y subsystem. More importantly, it allows to determine the statistical quantities of the y system that dictate this response.
The first order correction to the dynamics of the x system can be written as the expectation value ε dyψ x (x, y)ρ y (y), where ρ y is the invariant density of the uncoupledẏ = g y (y) dynamics. At second order two correction terms appear, one due to double ψ x interactions from y to x, determined by a correlation function of the uncoupled y dynamics, and a feedback term, determined by a response function of the uncoupled y dynamics. This knowledge can then be exploited to derive a surrogate dynamics for x that reproduces the effect of the coupling of x to y up to second order in ε. While this theory has been originally developed assuming that the uncoupled systems are Axiom A dynamical systems, it can be equally applied in the case where the uncoupled dynamics is stochastic, the only needed requirement being to have a physical measure. Interestingly, the results obtained using response theory match what one can derive by constructing a perturbative expansion of the dynamics of the system using the Mori-Zwanzig projection method (Wouters and Lucarini, 2013).
Previously, we have proposed a surrogate dynamical equation for the x variable that introduces an ε-dependent perturbing term to the dynamics f x to match the response of the statistics of the full system. The perturbing term contains a non-Markovian memory term and a correlated noise, with the memory kernel and correlation functions depending on the statistics of the uncoupled dynamicsẏ = g y . In a recent study of the applicability of the weak coupling approach to a simple ocean-atmosphere system, the method has been shown to give a good result for sufficiently weak coupling between the ocean and the atmosphere (Demaeyer and Vannitsem, 2016), even if it is clear that a systematic investigation of the performance of the weak coupling approach is indeed still needed.
We remark that Chekroun et al. (Chekroun et al., 2015b, a) have recently proved that, indeed, constructing reduced order models entails introducing deterministic, stochastic and memory correction to the dynamics of the variables of interest.
Here we will apply and extend the weak coupling approach of (Wouters and Lucarini, 2012;Wouters and Lucarini, 2013) for the development of parameterizations for various stochastic triad models. Triad interactions arise from quadratic nonlinearities with energy conserving properties (see e.g., (Gluhovsky and Tong, 1999)). The triad models considered here appear in applications of the homogenization technique to construction of parameterizations in climate modeling (see e.g., (Majda et al., 2001(Majda et al., , 2002Franzke et al., 2005;Franzke and Majda, 2006;Achatz et al., 2013;Dolaptchiev et al., 2012)). The non-Markovian memory kernel in the weak coupling approach will be calculated for these simple stochastic multiscale models and approximated by a Markovian stochastic process, in order to allow for easier numerical implementation. The systems we investigate can be written in both the weak coupling form of Eq. 2 and the slow-fast form of Eq. 1, therefore direct comparison is possible and will be performed on a number of metrics, namely initial ensemble spread, correlation functions and exit times from an interval.

The additive triad
The first model we look at is the stochastically forced additive triad. This system is a low-dimensional model that has non-linear interactions reminiscent of those occurring between the Fourier modes of a fluid flow. It is stochastically forced to mimic the interaction with further unresolved modes. The system has three variables, one slow variable x and two fast variables y 1 and y 2 . The fast dynamics is dominated by two independent Ornstein-Uhlenbeck processes. The dynamical equations for this triad The processes ξ i are independent Brownian motions in the Itô sense. Here and below a differential equation featuring a Brownian motion will be interpreted as the equivalent stochastic integral. In addition, we require i B (i) = 0, which guarantees energy conservation in the case γ i = σ i = 0.

Homogenization
On the time scale t, when increasing the time scale separation 1/ε to infinity, we have trivial dynamics of the averaged equationṡ By expanding the backward Kolmogorov equation for the slow-fast system in orders of ε, a Kolmogorov equation for only the slow variables can be derived (see (Pavliotis and Stuart, 2008)). The dynamical equation corresponding to this Kolmogorov equation is in this case a one-dimensional Ornstein-Uhlenbeck process (Majda et al., 2002) ∂x where See Fig. 1 for an illustration of the homogenization principle for the additive triad. The mean and variance of the triad converge to those of the Ornstein-Uhlenbeck process (4) for small ε.

Weak coupling limit
We will now discuss the weak coupling method as described in (Wouters and Lucarini, 2012;Wouters and Lucarini, 2013).
By rescaling the time as τ = ε −1 t we can write the stochastically forced additive triad equation (3) as a two-dimensional Ornstein-Uhlenbeck system weakly coupled non-linearly to a trivial zero-gradient x system: with ψ x (y 1 , y 2 ) = B (0) y 1 y 2 and ψ y (x, y) = (B (1) xy 2 , B (2) xy 1 ) T . The stochastic parametrization derived in (Wouters and Lucarini, 2012;Wouters and Lucarini, 2013) is given by a non-Markovian equation where the the memory kernel R(s,x) and first two moments of the stochastic process σ(τ ) are derived using the weak coupling method to the following statistics of the uncoupled y Ornstein-Uhlenbeck dynamics: where the evolution of y 1 and y 2 into y 1 (τ ) and y 2 (τ ) are taken to be the uncoupled Ornstein-Uhlenbeck dynamics dy i /dτ = −γ i y i + σ i ξ i . We have for the case of the additive triad (3) and

Markovian parametrization
Due to the identical time-scale γ 1 + γ 2 in both memory and noise correlation, the memory equation (6) can be transformed to a Markovian parametrization. We want to find a parametrizing two level Markovian dynamical system of the form such that the second order response of this system to changes in ε is the same as the response of (6). In other words, we want to determine the parameters C 1 , C 2 , γ and σ z in (11) such that the correlation and memory functions of the fast equation in (11) are equal to (9) and (10) respectively. The correlation function C(τ ) and memory function R(τ ) of the fast equation of (11) are where the evolution of z 2 to z 2 (τ ) is now given by dz 2 /dτ = −γz 2 + σ z ξ(τ ). By equating these functions to their counterparts in (9) and (10) we see that by choosing the reduced z 1 dynamics of the parametrized dynamical system in the weak coupling method are of the same form as those of the stochastic triad (3).
This Markovian reduced equation (11) is in fact a reformulation of the non-Markovian equation (6). To see this, we write an explicit solution for z 2 in function of the history of z 1 and ξ as This solution can then be inserted into (11), to obtain which agrees with (6), the first two terms being an Ornstein-Uhlenbeck process with the required correlation plus a memory term with the required memory kernel.
This Markovian formulation allows for a straightforward numerical implementation of the parametrization, compared to the non-Markovian equation (6) which requires one to store the history of the process in memory.
A comparison of the performance of the two model reductions is show in Figure 2. Shown are the spread of an ensemble initiated at a fixed value for the slow variables x = z 1 = −5 and the autocorrelation function of the slow variables. The weak coupling method clearly gives better results.
By correctly rescaling time and taking the limit of ε → 0 in the Markovian parametrization (11) one can furthermore verify that in this limit it converges to the homogenization of the original triad equation (Eq. (4)).

The slowly oscillating additive triad
The additive triad as specified in Eq.
(3) can be generalized to allow for an additional interaction between the y variables on the slow time scale that is independent of x. The dynamical equations for this slowly oscillating triad are dx dt = B (0) y 1 y 2

Homogenization
The homogenized equation is similar to the one for the additive triad, with an added constant forcing C r in the reduced SDE

Weak coupling limit
The coupling functions ψ x and ψ y are now The correlation function (7) of the coupling to x, determining the correlations of the parametrization noise σ is The response function (8) of ψ x to ψ y , determining the memory kernel of the parametrization, is similar to the one for the additive triad, with an added exponential function, the integral of which gives the same constant C r of the homogenized equations R(τ, x) = ψ y (x, y)∂ y ψ x (y(τ )) = exp(−γτ )(D 1 x + D 0 ) Combined, this then results in the following non-Markovian parametrized equations

Markovian parametrization
The non-Markovian equation (16) can again be Markovianized by a two-level Ornstein-Uhlenbeck process of the form The corresponding correlation and memory terms are We can therefore take In the limit ε → 0 in the Markovian parametrization (17) we again recover the homogenized equations.

Exit times
When comparing initial ensemble spread and autocorrelation functions for the slow variable of this system with the weak coupling parametrization and the homogenized system, the results are similar to those presented for the additive triad above.
Additionally, here we perform a comparison of a rare event statistic, the first exit time of the slow variable from an interval [−1, 1] when the slow variable is initialized at 0. The results in Tables 1 and 2 show that the statistics of exit times are significantly better approximated in the weak coupling parametrization.

The rapidly oscillating additive triad
A further generalization of the additive triad (3) is to introduce an interaction between the y variables on the fast time scale (Dolaptchiev et al., 2012). The dynamical equations for the rapidly oscillating triad are Note the difference in scaling on the oscillatory terms ωy i compared to Eq. (15). The invariant measure of the fast system is a correlated Gaussian measure ρ(y) = exp(−y T (2R) −1 y)/Z determined by Homogenization leads to a solvability condition on the system 20 that is fulfilled if either ω = 0 or σ 2 1 /γ 1 = σ 2 2 /γ 2 . The homogenized equation is now given bẏ
The correlation function ψ x (y 1 , y 2 )ψ x (y 1 (t), y 2 (t)) appearing in the weak coupling expansion can again be calculated explicitly. Solutions of the fast Ornstein-Uhlenbeck systemẏ = −Γy + Σξ can be written as

Markovian parametrization
Guided by the Markovian form of the previous triad systems, we again want to derive a Markovian parametrization with a reduced one-level Ornstein-Uhlenbeck system as the fast component: In this case, there is no exact match between the auto-correlation and response functions of this Markovian system and the non-Markovian weak coupling parametrization. The choice of the parametrization parameters is therefore not exactly determined and one needs to choose a parametrization such that the auto-correlation and response functions of the coupling function in the fast component of the full system are approximated in some sense. A further restriction comes from the fact that in the limit ε → 0 the limiting path distribution of the full system is determined by the homogenized equation and we therefore want to retain this limiting behavior in the parametrized system. To have this limiting property, we have the following constraints on the parameters in Eq. (23) where A ω and C ω are the forcing and friction parameters obtained through homogenization. For formal equivalence between the reduced and full equations, we furthermore set C 1 = B (0) . With the remaining free parameters we can match the response and correlation functions in a more precise manner, for example by matching the values of these functions at time t = 0. In this way, we get where h y = h ω /x.
A simulation of the ensemble spread from a fixed initial condition is shown in Figure 3. It demonstrates that the weak coupling parametrization (23) outperforms the homogenized reduced system.

Exit times
The same experiment on exits from an interval has been performed as described in Section 3.3. The results are displayed in Table 3. As before, the weak coupling reduced system gives a much better result when compared to the homogenized system.  Table 4. The relative error on the standard deviation of the exit times |σ1(τ ) − σ0(τ )|/σ0(τ ) where σ0(τ ) is the standard deviation of exit times from [−1, 1] of the full triad system and σ1(τ ) is the standard deviation of exit times of the parametrized systems. The parameters are the same as those used for Fig. 3.

The multiplicative triad
A final type of interactions is given by the multiplicative triad equations (Majda et al., 2002) which describes the interplay between two x modes and a stochastically forced single y mode. In the absence of forcing and dissipation energy conservation is satisfied if i B (i) = 0. In the system (24) the y mode can be eliminated directly by integrating the last equation of (24) Inserting this result in the equations for the x variables, one obtains Note that the first two term on the righthand side. result from a Ornstein-Uhlenbeck process with zero mean and stationary time autocorrelation function given by σ 2 m 2γm e − γm ε t .
The resulting parametrization in the weak coupling approach (Wouters and Lucarini, 2012;Wouters and Lucarini, 2013) reads with a noise term σ i with zero mean and correlation given by The memory kernel has the form Thus (27) can be written as which is exactly the same result as in (25), if we rescale time and assume as initial condition x 1 (t) = x 2 (t) = 0 for t < 0. In this case the weak coupling approach recovers exactly the full model. The original three component system was reduced to a two component non-Markovian system but there is no efficiency gain using the parametrization since the corresponding Markovian system is again a three component one.

Homogenization
Introducing a longer time scale θ = ε 2 τ in (28) and taking the limit ε → 0 one recovers the homogenization result in Stratonivich formulation The latter corresponds to an Itô stochastic differential equation of the form d dθ For a comparison of the statistics of the multiplicative triad and of the homogenization model we refer to (Majda et al., 2002).
In this work we have worked out a first application of the weak coupling response method of (Wouters and Lucarini, 2012;Wouters and Lucarini, 2013) to a multiscale stochastic system. By the choice of system we were able to perform both homogenization and the weak coupling reduction on this system, thereby allowing for direct comparison between the two reductions.
The response method yields a non-Markovian equation, making it cumbersome to integrate numerically. We have demonstrated here that for the systems studied the non-Markovian equation can be further reduced to a Markovian equation. Even with this further reduction the system gives a better match to the original system than the homogenized equations.
In the case of the additive triads, the system that is obtained through the Markovianization procedure is of intermediate complexity, between the full system and the homogenized limit. In the systems studied here, the retention of a fast time scale in the reduced system means that the reduction in simulation complexity is modest (one variable instead of two and a linear coupling instead of a nonlinear one). In the case of the multiplicative triad the weak coupling parametrization recovers exactly the full model and there is no efficiency gain. In many applications of practical relevance, however, one considers situations where the number of degrees of freedom of the unresolved variables is considerably larger than those of the slow variables of interest. A reduction to a system of one or a few variables will constitute a significant reduction in complexity in this case. This approach can be compared to the superparametrization approach to convection, where convection is parametrized by a model that is still dynamical in nature, yet significantly simpler than the full convective motion (Randall et al., 2003;Majda, 2013, 2014).