Improving solution accuracy and convergence for stochastic physics parameterizations with colored noise

TL;DR

This paper extends the Itô correction for stochastic differential equations to include colored noise, demonstrating improved accuracy and convergence in numerical weather prediction models with faster, more efficient simulations.

Contribution

A generalized Itô correction for colored noise is derived and applied, enhancing solution accuracy and convergence in stochastic PDEs used in climate modeling.

Findings

Reduces time integration error with colored noise.

Enables longer time steps for same accuracy.

Improves convergence rate in stochastic simulations.

Abstract

Stochastic parameterizations are used in numerical weather prediction and climate modeling to help capture the uncertainty in the simulations and improve their statistical properties. Convergence issues can arise when time integration methods originally developed for deterministic differential equations are applied naively to stochastic problems. (Hodyss et al 2013, 2014) demonstrated that a correction term to various deterministic numerical schemes, known in stochastic analysis as the It\^o correction, can help improve solution accuracy and ensure convergence to the physically relevant solution without substantial computational overhead. The usual formulation of the It\^o correction is valid only when the stochasticity is represented by {\it white} noise. In this study, a generalized formulation of the It\^o correction is derived for noises of any color. The formulation is applied to a test problem described by an advection-diffusion equation forced with a spectrum of fast processes. We present numerical results for cases with both constant and spatially varying advection velocities to show that, for the same time step sizes, the introduction of the generalized It\^o correction helps to substantially reduce time integration error and significantly improve the convergence rate of the numerical solutions when the forcing term in the governing equation is rough (fast varying); alternatively, for the same target accuracy, the generalized It\^o correction allows for the use of significantly longer time steps and hence helps to reduce the computational cost of the numerical simulation.

Equations

$$\frac{\partial u}{\partial t}=D(u)+P(u),$$

$$P(u)=g(u)H(t).$$

$$P_s(u)=g(u)\dot{R}(t),$$

$$\frac{\partial u}{\partial t}=D(u)+P_s(u).$$

$$u(t_2)-u(t_1)$$

$$\displaystyle\mbox{It\^{o} integral: }\int g(u)dB=\lim_{\Delta t\rightarrow 0}\sum_{j}g\big{[}u(t_{j}^{*})\big{]}\Delta B_{j}$$

$$\text{where}{t}_j^{\ast }=t_j\text{ (left endpoint),}$$

$$\displaystyle\mbox{Stratonovich integral: }\int g(u)\circ dB=\lim_{\Delta t\rightarrow 0}\sum_{j}g\big{[}u(t_{j}^{*})\big{]}\Delta B_{j}$$

$$\text{where}{t}_j^{\ast }=\frac{t_j+t_{j+1}}{2}=t_j+\frac{\Delta t}{2}\text{ (midpoint).}$$

$$\displaystyle g\big{[}u(t_{j}^{*})\big{]}$$

$$\dot{R}(t_j)\Delta t\approx \Delta R_j.$$

$$\displaystyle g\big{[}u(t_{j}^{*})\big{]}\approx g\big{[}u(t_{j})\big{]}+\left(\frac{1}{2}\frac{dg(u)}{du}{D}(u)\right)\!\!\bigg{|}_{t_{j}}\!\!\Delta t+\left(\frac{1}{2}\frac{dg(u)}{du}g[u]\right)\!\!\bigg{|}_{t_{j}}\!\!\Delta{R}_{j}$$

$$\displaystyle+\frac{\Delta t^{2}}{8}\left(\frac{d^{2}g[u(t)]}{dt^{2}}\right)\!\!\bigg{|}_{\xi}$$

$$u(t_2)-u(t_1)=\underset{\text{mbox}traditionalintegral}{\underbrace{{\int }_{t_1}^{t_2}D(u)dt}}+\underset{\text{mbox}It\hat{o}integral}{\underbrace{{\int }_{t_1}^{t_2}g(u)dR}}$$

$$\displaystyle+\lim_{\Delta t\rightarrow 0}\sum_{j}\left(\frac{1}{2}\frac{dg(u)}{du}g[u]\right)\!\!\bigg{|}_{t_{j}}\!\!\left(\Delta{R}_{j}\right)^{2}$$

$$\displaystyle\lim_{\Delta t\rightarrow 0}\sum_{j}\left(\frac{1}{2}\frac{dg(u)}{du}g[u]\right)\!\!\bigg{|}_{t_{j}}\!\!\mathbb{E}\big{[}\left(\Delta{R}_{j}\right)^{2}\big{]}\,,$$

$$\mathbb{E}\big{[}\left(\Delta{R}_{j}\right)^{2}\big{]}=\Delta t,$$

$${\int }_{t_1}^{t_2}\left(\frac{1}{2}\frac{dg(u)}{du}g(u)\right)dt,$$

$$\frac{\partial u_i}{\partial t}=D_i(u)+\sum _{l=1}^p{g}_{il}(u){\dot{R}}_l(t),\text{for}i=1,\dots ,n$$

$$\displaystyle\lim_{\Delta t\rightarrow 0}\sum_{j}\left(\frac{1}{2}\sum_{l=1}^{p}\sum_{k=1}^{n}\frac{\partial g_{il}(u)}{\partial u_{k}}g_{kl}[u]\right)\!\!\bigg{|}_{t_{j}}\!\!\mathbb{E}\big{[}\left(\Delta R_{lj}\right)^{2}\big{]}\,.$$

$$u(t_2)-u(t_1)=\underset{\text{mbox}traditionalintegral}{\underbrace{{\int }_{t_1}^{t_2}D(u)dt}}+\underset{\text{mbox}It\hat{o}integral}{\underbrace{{\int }_{t_1}^{t_2}g(u)dR}}$$

$$\displaystyle+\lim_{\Delta t\rightarrow 0}\sum_{j}\left(\lambda\frac{dg(u)}{du}g[u]\right)\!\!\bigg{|}_{t_{j}}\!\!\left(\Delta{R}_{j}\right)^{2}\,.$$

$$\displaystyle\lim_{\Delta t\rightarrow 0}\sum_{j}\left((\frac{1}{2}-\lambda)\frac{dg(u)}{du}g[u]\right)\!\!\bigg{|}_{t_{j}}\!\!\mathbb{E}\big{[}\left(\Delta{R}_{j}\right)^{2}\big{]},$$

$$\frac{\partial u}{\partial t}=-\left[c+\frac{ϵ}{2}\cos (x)\right]\frac{\partial u}{\partial x}+\mu \frac{{\partial }^{2}u}{\partial x^2}+g(u)n(t),$$

$$c=1,\mu =0.1,$$

$$g(u)=\rho \frac{\partial u}{\partial x}\text{with}\rho =0.2.$$

$$I=\frac{1}{2}{\rho }^2\frac{{\partial }^{2}u}{\partial x^2}\underset{N_f\to \infty }{\lim }\frac{1}{N_f}\left[\frac{C({\omega }_0{)}^2}{2}+\sum _{m=1}^{N_f}C({\omega }_m{)}^2\right],$$

$$\underset{N_f\to \infty }{\lim }\frac{1}{N_f}\left[\frac{C({\omega }_0{)}^2}{2}+\sum _{m=1}^{N_f}C({\omega }_m{)}^2\right]=1,$$

$$u(x,t)=\sum _{k\in \mathbb{Z}}{F}_k(t)\exp (ikx)$$

$$u(x,0)=\cos (k_{0}x)\text{with}{k}_0=1.$$

$$\frac{d{F}_k}{dt}=-ikc{F}_k-\mu k^2{F}_k+ik\rho n(t)F_k.$$

$$F_k(t)=A\exp \left[-\left(ick+\mu k^2\right)t+i\rho k{\int }_0^{t}n(t^{\prime })d{t}^{\prime }\right]$$

$$\begin{array}{rl}\sum _{k=-N_x}^{N_x}\frac{d{F}_k(t)}{dt}\exp (ikx)=& \sum _{k=-N_x}^{N_x}\{-ik{F}_k(t)\exp (ikx)\left[c+\frac{ϵ}{2}\cos (x)\right]\\ & -\mu k^2{F}_k(t)\exp (ikx)+ik\rho F_k(t)\exp (ikx)n(t)\frac{}{}\}.\end{array}$$

Full text

1 Introduction

Physical and chemical processes happening in the Earth’s atmosphere span many orders of magnitude in terms of their spatial and temporal scales, which presents great challenges to numerical modeling. For example, in general circulation models, motions or phenomena that can not be resolved in space or time but have significant impact on the large-scale flow motions need to be accounted for using parameterizations (see e.g. [12].)

In recent years, stochastic parameterizations have become an active area of research (see review by [1, 11]). The fundamental principle behind the stochastic formulation is that the state of the unresolved processes at any instant is not entirely determined by the state of the resolved processes. Thus, an element of randomness needs to be introduced to account for this indeterminacy. This randomness can act as a source of roughness in the temporal evolution of the governing equations’ right-hand-side terms as well as in the evolution of the solution. Deterministic time integration schemes used in numerical weather prediction and climate projection models, however, typically assume temporal smoothness of the underlying solutions. When such schemes are applied naively to stochastic parameterizations, the conditions for solution convergence might no longer be satisfied.

There is a large body of work dedicated to the development of stochastic numerical schemes and stochastic versions of deterministic numerical schemes especially targeted for weather and climate applications (see e.g. [15, 18, 16, 3, 4, 6]). The purpose of our work is more modest. We want to study the use of a correction term that can help improve the solution accuracy of deterministic schemes when part of the variables of the model is replaced by noise and in particular colored noise.

As shown in [7], the convergence issue of deterministic numerical schemes when applied to stochastic parameterizations can be investigated through the use of tools from stochastic analysis (see e.g Section 3.3 in [13] and Section 4.9 in [10]). In particular, for the cases when an unresolved process is replaced by a rough random process (e.g. white noise), it is not difficult to construct examples for which popular deterministic numerical schemes (e.g, Euler forward and backward, Adams-Bashforth) will no longer converge to the physically relevant solution except for special cases (e.g., the second-order Runge-Kutta scheme analyzed by [7]). Multiple examples relevant for atmospheric modeling can be found in Figure 3 in the paper of [7]. Here, “physically relevant solution” refers to the one corresponding to ordinary calculus (see discussion below). In the study by [8], ensemble simulations of Hurricane Isaac in the year 2012 were conducted using the Navy Operational Global Atmospheric Prediction System; it was shown that the choice of numerical scheme for the stochastic term can lead to failure in predicting the correct ensemble mean of hurricane intensity.

The mathematical reason for the lack of convergence is that when we replace an unresolved process with white noise, the equations describing the phenomena under investigation make sense only in integral form (not in the usual differential form). The integral form of the equations contains a temporal integral of an expression involving the white noise process. If we try to estimate such an integral through a limiting process involving progressively refined subintervals, different answers will be obtained depending on the manner we choose to discretize the interval of integration (see e.g. Section 4.9 in [10]). The two most-studied discretization methods in stochastic analysis are: i) using the left endpoint of each subinterval (which leads to the Itô integral or Itô interpretation) and ii) using the middle point of each subinterval (which leads to the Stratonovich integral or Stratonovich interpretation). The Stratonovich interpretation leads to ordinary calculus while the Itô interpretation does not (see e.g. Sections 3.9 and 4 in [13]). Recall from the previous paragraph, that the physically relevant solution for the systems the weather and climate researchers are attempting to model is the one corresponding to ordinary calculus. Thus, the solution resulting from the Stratonovich interpretation is the physically relevant one [7].

It is important to note that many popular time integration schemes designed for deterministic problems will converge to the Itô solution when applied to stochastic problems driven by white noise [10]. In other words, naively describing an unresolved process by white noise and solving the stochastic equation with a deterministic numerical scheme can lead to erroneous results even in the limit of infinite temporal resolution. Fortunately, the Itô and Stratonovich interpretations are related, and this relationship can help recover, at least to some extent, the convergence of deterministic numerical schemes to the physically relevant Stratonovich solution. The connection between the two interpretations comes in the form of a correction term called the Itô correction. When the Itô correction is added to the equation, the numerical solution under the Itô interpretation converges to the Stratonovich solution.

While the above-mentioned Itô-Stratonovich correspondence is a basic concept in stochastic analysis, the widely known form of the Itô correction applies only to the case of white noise. A key feature of white noise is that it has zero auto-correlation (and hence no memory). Given the typical time step size of seconds to an hour in weather and climate models, some parameterized processes (e.g., turbulence and cumulus convection) can have characteristic time scales equivalent to multiple time steps. Therefore, colored noise, which has non-zero autocorrelation length, can provide a better description of such processes. In fact, the state of the art in accounting for model uncertainties of Earth systems points to the need of stochastic processes with spatio-temporal correlations (see e.g. Sections 3 and 5 in [11]), which makes our construction more relevant for applications.

A fundamental difference between colored noise and white noise is that colored noise is in principle resolvable while white noise is not. In other words, if one could use small enough step sizes, there would be no distinction between the Itô and Stratonovich interpretations for the case of colored noise. All deterministic numerical schemes will eventually converge to the Stratonovich solution. But even in simple examples, let alone the very complex and expensive systems encountered in weather and climate prediction, the critical timestep that recovers convergence to the Stratonovich solution can be prohibitively small. As a result, for realistically affordable time step sizes, the dichotomy between Itô and Stratonovich interpretations practically exists and needs to be addressed also for the cases of colored noise. In other words, the long time steps in practical applications motivate us to find time-stepping methods with higher accuracy.

Another point worth mentioning is that, as [7] and [8] have pointed out, certain deterministic numerical schemes (e.g., the second-order Runge-Kutta scheme) have (an unaveraged version of) the Itô correction already “built-in” and hence perform better for stochastic problems. Such schemes are typically multi-stage schemes which require multiple evaluations of the right-hand side of the governing equations, making them very expensive for weather and climate models. The Itô correction, in contrast, allows for the use of single-stage schemes (e.g. the Euler forward scheme) to be complemented by a correction constructed only for the stochastic term, and hence can be cost-effective. It should be noted that although the Itô integral can be viewed as applying the Euler forward scheme to the stochastic term, other discretization methods can also fail to converge to the Stratonovich integral, see Section 3 of [7]. This means that the addition of an Itô correction could aid in the restoration and/or acceleration of convergence to the Stratonovich solution also in the case of other numerical schemes.

For these reasons, we present in this paper a generalization of the Itô correction that is valid for noises of any color. We use an advection-diffusion equation with constant or spatially varying advection velocity to demonstrate that, for both white and colored noises, the generalized Itô correction can accelerate convergence to the Stratonovich solution when added to the Euler forward scheme. We demonstrate that improved convergence means higher accuracy for the same step size or, alternatively, larger step size (and hence lower computational cost) for the same target accuracy. These results from the simple but relevant test problem provide a proof of concept, which motivates further exploration of the generalized Itô correction for the purpose of helping improve the solution accuracy and efficiency in atmospheric models of various complexity, including general circulation models using stochastic parameterizations.

The remainder of the paper is organized as follows: Section 2 presents the derivation of the generalized Itô correction. Section 3 contains a presentation of the test problem, the advection-diffusion equation with constant and spatially varying advection velocity, along with analytical results (supplemental details can be found in the Appendix). Section 4 contains numerical results. Finally, Section 5 contains a discussion of our results as well as suggestions for future work.

2 The generalized Itô correction

We consider the following deterministic differential equation

[TABLE]

where $D(u)$ and $P(u)$ are the resolved dynamics and parameterized physics, respectively. Here we focus on the special case where $P(u)$ takes the form

[TABLE]

This form results from the attempt to eliminate a fast-evoling physical quantity from the original equations and replace it by a time-dependent process. We can consider a more general form where $H$ depends also on the spatial variable, but that generalization will not alter the derivation of the generalized Itô correction below, hence we restrict our attention to the case where $H$ depends only on $t$.

If the time scales associated with ${H}(t)$ are substantially shorter than the time scales of ${D}(u)$, we can approximate ${P}(u)$ by its stochastic counterpart ${P}_{s}(u)$ defined as

[TABLE]

where $\dot{{R}}$(t) represents a general noise term. We note here that replacing $H(t)$ by the noise term, $\dot{{R}}(t),$ may include a limiting process where the function $g(u)$ may be also modified (see e.g. [14]).This does not alter the main line of our derivation and we keep the notation $g(u)$ for the multiplicative factor.

Using Eq. (3), we get the following stochastic counterpart of the deterministic equation originally given by Eq. (1):

[TABLE]

Without loss of generality, we assume $\mathbb{E}[\dot{{R}}(t)]\equiv 0$, where $\mathbb{E}[\cdot]$ denotes the mean over different realizations of the noise process. If $\mathbb{E}[\dot{{R}}(t)\dot{{R}}(t^{\prime})]=\delta(t-t^{\prime})$ where $\delta(\cdot)$ is Dirac’s delta function, then $\dot{{R}}(t)$ is white noise and ${R}(t)$ is a Wiener process; when $\mathbb{E}[\dot{{R}}(t)\dot{{R}}(t^{\prime})]\neq\delta(t-t^{\prime})$, $\dot{{R}}(t)$ is a colored noise.

We focus on how to numerically solve Eq. (4) after its form has been derived; how to construct a good $P_{s}(u)$ to approximate the original $P(u)$ is a separate topic which is outside the scope of the current work.

2.1 Derivation

Let us take the integral over an arbitrary time window $(t_{1},t_{2})$ on both sides of Eq. (4). For ${P}_{s}(u)$, we discretize the time interval into $J$ bins of equal length $\Delta t$ and denote the increment of $R$ in each bin as $\Delta R$. We use $t_{j}^{*}$ to denote the discretization point inside the $j$-th bin, i.e., the instant where the value of ${P}_{s}(u)$ is evaluated for numerical integration. With this notation, the integral of Eq. (4) can be written as

[TABLE]

In the white noise case (i.e., ${R}(t)=B(t)$ where $B(t)$ is the Wiener process), the choice of discretization point for the integral can lead to different results [13]. The two most popular choices are defined as

[TABLE]

and

[TABLE]

Here $\Delta B_{j}=B_{j+1}-B_{j}.$ Because the physical processes represented by the deterministic equation (1) are assumed continuous, the Stratonovich integral should be used in our case (see Section 3.3 in [13]).

Before proceeding further with the derivation, we note that the Itô interpretation for the stochastic integral as shown in Eq. (6) coincides with how the Euler forward approach would treat the stochastic term, but the correspondence between Itô and Stratonovich interpretations that we derive below is not tied to the forward Euler scheme. We will come back to this point in Section 22.3.

It is well known in stochastic analysis that in the white noise case, the Stratonovich integral can be written as the sum of an Itô integral and a correction term called the Itô correction (see, e.g., [13]). Below we show that the same is true for colored noise, although the Itô correction needs to be generalized.

For $t_{j}^{*}=\left(t_{j}+t_{j+1}\right)/2,$ performing a Taylor expansion of g\big{[}u(t_{j}^{*})\big{]} about $t_{j}$ and expressing $\partial u/\partial t$ using Eq. (4) gives

[TABLE]

where $\xi\in\left[t_{j},\left(t_{j}+t_{j+1}\right)/2\right]$. For small $\Delta t$, we write

[TABLE]

Hence, Eq. (9) can be approximated as

[TABLE]

Assuming $g(u)$ is sufficiently smooth and $\Delta t$ is small, one can neglect the second and fourth terms on the right-hand side of Eq. (11) but, in general, not the third term. Therefore, with the Stratonovich interpretation of the stochastic integral in Eq. (5), we have

[TABLE]

The mathematical expectation of the last term in Eq. (12) is

[TABLE]

which is the generalized Itô correction in its integral form. The exact form of the expectation in expression (13) depends on the formulation of ${R}.$ For example, when $R$ is the Wiener process, the increment $\Delta R_{j}$ is a Gaussian random variable with mean 0 and variance $\Delta t$, i.e.,

[TABLE]

hence (13) becomes

[TABLE]

which is the integral form of the traditional Itô correction (see, e.g. Section 3.3 in [13]).

The generalized Itô correction (13) can be extended to the case of multiple partial differential equations (PDEs) each containing multiple noise processes. Let us assume a system of $n$ PDEs for the functions ${\bf u}=(u_{1},u_{2},\ldots,u_{n})$:

[TABLE]

where $\dot{R_{l}}(t)=(\dot{R_{1}}(t),\dot{R_{2}}(t),\ldots,\dot{R_{p}}(t))$ is a $p$-dimensional vector noise process with independent components. Then, the expression for the generalized Itô correction for the equation for $u_{i}$ is given by

[TABLE]

where $\Delta R_{lj}=R_{l,j+1}-R_{lj}$ is the increment of the $l$th noise process $R_{l}.$

2.2 Remarks

We want to make two remarks concerning the derivation of the generalized Itô correction (13). First, there is an alternative way to derive the generalized Itô correction. In particular, under the assumption that the correlation time of the noise is short, one can employ the expansion devised by Stratonovich (see Section 4.8 in [17]), through which a stochastic equation driven by colored noise can be rewritten as an effective stochastic equation driven by white noise. Then, one can compute the traditional Itô correction for the resulting white noise driven equation.

Second, the Itô correction, in its traditional or generalized form, can be interpreted as a memory term encountered in model reduction formalisms (see e.g. [2]). By using as discretization point the left endpoint of each interval, the Itô interpretation of the stochastic integral makes the evolution of the stochastic process $\dot{{R}}(t)$ independent of the solution $u(t).$ The Itô correction serves as a way to account for the interaction of $\dot{{R}}(t)$ and $u(t)$ during the interval $\Delta t,$ similar to the role played by memory terms in model reduction which account for the interaction between resolved and unresolved variables.

2.3 Applicability

It has been stated earlier in Section 22.1 that our generalized Itô correction (13) is not tied to the specific discretization method (e.g., Euler forward) that is chosen for the time integral of the stochastic term in Eq. (4). The reason is that for any discretization, as long as an analysis similar to Eqs. (9)–(11) reveals that the discretized integral converges to the Itô integral, expression (13) can be used to obtain numerical results that converge to the Stratonovich solution.

The method of analysis demonstrated by Eqs. (9)–(11) can also be applied to the “decentered” time-stepping methods commonly used in atmospheric models. Since these methods approximate time integrals (or derivatives) using the discretization point $t^{*}_{j}=(1-\lambda)t_{j}+\lambda t_{j+1}=t_{j}+\lambda\Delta t$ where $0\leq\lambda\leq 1$, Eq. (5) becomes

[TABLE]

A comparison of Eq. (18) with Eq. (12) suggests that the correction term linking the “decentered” integral and the Stratonovich integral is

[TABLE]

which is a further generalization of expression (13) for non-zero $\lambda$. In the special case of $\lambda=0$, the “decentered” scheme gives the Itô interpretation of the stochastic integral and we recover (13); the case of $\lambda=1/2$ corresponds to the Stratonovich interpretation of the integral and the correction vanishes. For the case when $R$ is the Wiener process, expression (19) corresponds to the correction formula for white noise that is found in [7] (see equations (2.7) and (3.21) therein) as well as in Section 3.5 of [10].

The example of decentered schemes discussed above can be generalized even further: for a generic discretization method, an analysis similar to Eqs. (9)–(12), followed by a comparison with the desired interpretation of the equation (i.e., Eq. (12) for Stratonovich or Eq. (12) without the last right-hand-side term for Itô), can lead to the correction term needed to obtain numerical results converging to the desired type of solution (Stratonovich or Itô).

2.4 Integral versus differential form

The expressions for the generalized Itô correction presented so far have been obtained from the integral form of the stochastic equation, while in the literature on stochastic analysis, the Itô correction conventionally denoted by the symbol $I$ is typically the term that is added to the differential form of the stochastic equation. The differential form of our generalized Itô correction in the test problem discussed below is given in Section 3 and in Appendix Appendix A. An example showing how the differential form is obtained from the integral form can be found in Appendix Appendix A (Eqs. 49–53).

3 Test problem

In the remainder of the paper, we use an example to demonstrate the impact of the generalized Itô correction. We consider the following stochastic differential equation

[TABLE]

with initial condition $u(x,0)=u_{0}(x)$ and periodic boundary conditions on $[0,2\pi].$ In the context of atmosphere modeling, the first two terms on the right-hand side represent the resolved dynamics and the last term represents fast varying physics parameterizations. When the parameter $\epsilon$ is set to 0, we recover the advection-diffusion equation with constant advection velocity discussed in [7]. The inclusion of $\epsilon\cos(x)/2$ in the first right-hand-side term makes the advection velocity spatially varying. In Section 4, numerical results are shown for both $\epsilon=0$ and $\epsilon=10^{-3}$. Following [7], we let

[TABLE]

The stochastic noise process $n(t)$ is the same as described in Appendix A of [7] (also described in Appendix A of this paper). For this choice of $g(u)$ and $n(t)$, the generalized Itô correction is given by

[TABLE]

where $N_{f}$, $C$, $\omega_{0}$ and $\omega_{m}$ are parameters of the noise process $n(t)$ (cf. Appendix A). The expression for $C$ contains a parameter $\alpha$ that controls the color of the Fourier spectrum of $n(t)$, with $\alpha=0$ corresponding to white noise, and larger $\alpha$ values corresponding to noise spectra that are more red.

Before we proceed, we want to make an important remark about the formula for $I.$ For the case of white noise ($\alpha=0$), we obtain

[TABLE]

and we recover the usual Itô correction expression. However, for the case of colored noise ($\alpha\neq 0$) with exponentially decaying spectrum, we have $\lim_{N_{f}\rightarrow\infty}\frac{1}{N_{f}}\!\left[\frac{C(\omega_{0})^{2}}{2}+\sum_{m=1}^{N_{f}}C(\omega_{m})^{2}\right]=0.$ This result is not surprising. As we have explained also in the introduction, in the case of colored noise, the distinction between the Itô and Stratonovich interpretations disappears in the limit of infinite temporal resolution. The reason is that in the limit of infinite temporal resolution a non-white colored noise is resolved and thus all discretizations of the integral of the stochastic term give the same answer. However, any numerical experiment that one conducts always has finite temporal resolution. In this case, the generalized Itô correction for the case of colored noise is no longer zero. Moreover, as we show with our numerical results, it can play a significant role in restoring or accelerating convergence to the Stratonovich solution. This remark is particularly pertinent for weather and climate applications where due to computational limitations we are always forced to use larger timesteps than the shortest timescales present in the solution.

To derive analytical solutions for the test problem, we express $u(x,t)$ in the form of a superposition of Fourier modes

[TABLE]

and transform Eq. (20) into a system of stochastic differential equations. Here $i$ is the imaginary unit. Like in [7], we assume the initial condition contains only one mode, i.e.,

[TABLE]

3.1 Case with constant advection velocity ($\epsilon=0$)

As pointed out by [7], when the advection velocity is constant, the Fourier modes are uncoupled. The ordinary differential equation (ODE) for $F_{k}(t)$ reads

[TABLE]

The analytical solution of Eq. (26) takes the form

[TABLE]

with $A$ being any complex constant. Initial condition (25) implies that $A=1$ in Eq. (27); only one Fourier mode (the one corresponding to $k=1$) is sufficient to represent the solution (the Fourier mode for $k=-1$ is also needed but due to the solution of (20) being real, it is the complex conjugate of the solution for the Fourier mode with $k=1.$)

3.2 Case with spatially varying advection velocity ($\epsilon\neq 0$)

For cases with nonzero $\epsilon$, even if the initial condition has a single Fourier mode, the spatially dependent component of the advection velocity causes the representation of the solution to require more than one mode. This is an elementary way to introduce coupling between different Fourier modes but still keep Eq. (20) linear.

We truncate the Fourier series in Eq. (24) to retain only modes with appreciable magnitudes and denote the largest remaining wavenumber as $N_{x}$. Substituting Eq. (24) for the unknown $u$ in (20) gives

[TABLE]

By multiplying Eq. (28) with $\exp(-ikx)$ and integrating over $[0,2\pi]$, we get the following coupled equations:

• •

for $k=-N_{x}+1,\dots,N_{x}-1$,

[TABLE]

• •

for $k=N_{x}$,

[TABLE]

• •

for $k=-N_{x}$,

[TABLE]

Using the notation defined in Appendix B, we can write the above stochastic ODE system for the Fourier mode coefficients $F_{k}$ in matrix form as

[TABLE]

The analytical solution reads

[TABLE]

4 Numerical Results

In this section we use numerical results to show how noise $n(t)$ of different color (roughness in time) can affect the solution convergence. We also demonstrate how the inclusion of the generalized Itô correction can help in restoring and/or accelerating convergence. As explained in Section 22.3, the validity of our generalized Itô correction is not tied to any specific time-stepping method. For simplicity, we use in our numerical experiments the forward Euler scheme as an illustrating example.

4.1 Definition of solution error

The error of a numerical solution is evaluated after two time units of integration using the $L_{2}$ norm ([7] and personal communication):

[TABLE]

Here $\widehat{u}$ and $u$ are the discrete and analytical solutions, respectively. To ensure the accuracy of the analytical solution computed for our error evaluation, the time integral of the noise process in Eqs. (27) and (34) is calculated analytically.

4.2 Case with constant advection velocity ($\epsilon=0$)

For the case with $\epsilon$ = 0, the discretization of Eq. (26) using forward Euler with the Itô correction included is given by

[TABLE]

where $n(t_{j})$ is the colored noise at $t=t_{j}$ and $I_{k}(t_{j})$ is the Itô correction for $F_{k}$ at $t_{j}$,

[TABLE]

We note that although the expression for the Itô correction in Eq. (23) involves a limiting process, the limit is not present in the expression in Eq. (37) because we have discretized the equation and thus have picked a finite timestep.

Panel (a) of Figure 1 (which appears also in [7]) shows the effect of different noises on the convergence of the Euler scheme without the Itô correction. The thick dots are the $l_{2}$ error of the numerical solution averaged over 100 realizations of the noise process; the error bars denote the standard deviation around the average. Using the terminology of stochastic analysis, this plot (and the rest of them in the paper) shows the strong convergence of the numerical solution. 111As a reminder, we note that strong convergence is measured by the mean of the solution error of individual realizations of the stochastic equation while weak convergence is measured by the error of the mean solution.

We make two observations. First, for the case of white noise ($\alpha=0$, purple line), the Euler scheme without the Itô correction fails to converge to the analytical solution no matter how small the step size is (it converges to the Itô solution, cf. Section 10.2 in [10]). Second, for the case of colored noise ($\alpha\neq 0$, blue, green, orange and red lines), the Euler scheme without Itô correction will start converging to the analytical solution with order 1 (as predicted by deterministic numerical analysis, see e.g. Chapter I.7 in [5]) when the step size becomes smaller than some critical step size which depends on the color of the noise (value of $\alpha$). The more red the noise is (larger $\alpha$), the larger is the critical step size (see also [7] for a discussion and estimation of the critical stepsize).

The right panel in Figure 1 shows the effect of including the Itô correction. We want to make again two observations. First, for the case of white noise ($\alpha=0$, purple line), the Euler scheme with the Itô correction does converge to the analytical solution with order $1/2$ (see [9] for an explanation of this convergence rate). We note that this numerical result was mentioned in [7] although not illustrated by any graphic there. Second, for the case of colored noise ($\alpha\neq 0$, blue, green, orange and red lines), the Euler scheme with the generalized Itô correction starts converging to the analytical solution with order 1 for larger step sizes than the Euler scheme without the Itô correction. Thus, the addition of the Itô correction can help restore and/or accelerate convergence of the forward Euler scheme.

4.3 Case with spatially varying advection velocity ($\epsilon\neq 0$)

We continue with the case of a spatially-dependent advection velocity with $\epsilon=10^{-3}.$ A small value was chosen for $\epsilon$ because the forward Euler scheme is explicit and only first-order. As such, it needs a very large number of steps in order to reach the asymptotic convergence regime for larger values of $\epsilon$ due to the need to resolve steepening gradients associated with the oscillatory nature of the spatial perturbation of the advection velocity. Moreover, the cost of evaluation of the noise $n(t),$ which depends quadratically on the number of timesteps, becomes very large when $\epsilon$ is large. For practical purposes (computational cost), we chose a small $\epsilon$ for the demonstration here.

We discretized Eq. (33) using the Euler scheme with the generalized Itô correction, i.e.,

[TABLE]

where the Itô correction reads

[TABLE]

with the matrix $\boldsymbol{\mathsf{G}}$ being

[TABLE]

The truncation wavenumber $N_{x}$ was chosen empirically: a test simulation was conducted using Eq. (38) with a large $N_{x}$; an inspection of the magnitude of the resulting $\widehat{F}_{k}$ revealed $N_{x}=5$ was sufficient to retain all modes with $|\widehat{F}_{k}|>10^{-4}$. Hence, $N_{x}=5$ was then used to obtain the results shown in Figure 2.

Figure 2 shows that for the case of white noise ($\alpha=0$, purple line), the forward Euler scheme without the Itô correction fails to converge to the analytical solution as expected. Figure 2 shows how the inclusion of the Itô correction can restore convergence with order $1/2.$ We note that the standard deviation bars around the mean appear larger than in the case with $\epsilon=0$ because of the logarithmic scale of the plot. This figure demonstrates that for the case of colored noise ($\alpha\neq 0$), the use of the generalized Itô correction again accelerates the establishment of the order 1 convergence regime predicted by deterministic numerical analysis.

5 Conclusions

Stochastic parameterizations are increasing in popularity in numerical weather prediction and climate modeling as a way to improve the statistical representation of the studied phenomena. Naive implementation of such parameterizations with deterministic numerical time integration schemes can cause serious convergence issues. Such issues can be alleviated by the addition of certain correction terms (known as the Itô correction in stochastic analysis) to the deterministic numerical schemes. However, the Itô correction was originally derived only for the special case when the stochastic process is represented by white noise. For numerical weather prediction and climate modeling it will be useful to have the option to properly handle colored noise.

We have derived a generalized Itô correction for the case of colored noise and applied it to a test problem of an advection-diffusion equation driven by noise of different colors. Our results indicate that the generalized Ito correction can substantially reduce the time discretization error, significantly improve the convergence rate of the numerical solutions and allow for the use of significantly larger step sizes.

While our derivation started from a stochastic differential equation, the fact that colored noise is in principle resolvable by sufficiently small step sizes implies that the generalized Itô correction can also be useful for deterministic problems for the purpose of improving solution convergence, accuracy, and efficiency. Compared to higher-order schemes like the Runge-Kutta family, the Itô correction is less costly in terms of computing time; compared to implicit methods that may provide better stability, the Itô correction is less intrusive in terms of the code modification it requires.

In the future, we plan to apply the current framework to more realistic atmospheric modeling problems, e.g. simplified versions of the atmospheric general circulation models or their parameterizations. We acknowledge the fact that the parameterizations will likely not be given directly in the multiplicative form required by our formulation. However, there is hope that new approaches e.g. training a neural network to represent the function $g(u)$ in Eq. (2) will allow us to still use our construction for more involved parameterizations.

The generalized Itô correction for the case of colored noise still assumes that the subgrid phenomena evolve on a significantly shorter timescale than the phenomena we resolve explicitly. However, as was mentioned in Section 2.2, the Itô correction can be interpreted as a memory term in model reduction formalisms. This opens the possibility of exploiting more sophisticated types of memory terms which correspond to more nuanced and realistic noise processes.

Acknowledgements

The authors thank Drs. Christopher J. Vogl (LLNL), Carol S. Woodward (LLNL), and Shixuan Zhang (PNNL) for helpful discussions on the numerical examples shown in this paper, and Dr. Daniel Hodyss (NRL) for clarifications regarding his earlier work that inspired our study. This work was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Advanced Scientific Computing Research (ASCR) and Office of Biological and Environmental Research (BER), Scientific Discovery through Advanced Computing (SciDAC) program. Pacific Northwest National Laboratory is operated by Battelle Memorial Institute for DOE under Contract DE-AC05-76RL01830.

Appendix A

Following [7], we use the following specification of the noise process $n(t)$:

[TABLE]

where $N$ is the number of discrete time levels per unit time, including the starting and ending time levels. $\alpha$ is a parameter controlling the color of the Fourier spectrum of $n(t)$ ($\alpha=0$ corresponds to white noise while $\alpha\neq 0$ to colored noise). To construct different realizations of the noise process, we sample, for $m=0,\ldots,N_{f},$ the coefficients $a_{m}$ and $b_{m}$ independently from the normal distribution $\mathcal{N}(0,1)$. It should be noted that $n(t)$ is an approximate random noise. The difference between $n(t)$ and the noise term $\dot{{R}}(t)$ in Section 2 is that $\dot{{R}}(t)$ would contain an infinite number of Fourier modes while $n(t)$ only has a finite number of modes. Nevertheless, in numerical modeling, we can use $n(t)$ to approximate $\dot{{R}}(t)$.

Let us define

[TABLE]

and

[TABLE]

and consider $\Delta\beta_{j}$ as an approximation to $\Delta{R}_{j}$ (recall that $\Delta{R}_{j}=R_{t_{j}+\Delta t}-R_{t_{j}}$ is the increment of the stochastic process $R$). For the above-defined $n(t)$, we find

[TABLE]

For small $\Delta t,$ we can approximate \mathbb{E}\big{[}\left(\Delta\beta_{j}\right)^{2}\big{]} as

[TABLE]

Taking into account the independence among the coefficients $a_{m}$ and $b_{m}$, we have

[TABLE]

Also note that per construction, we have

[TABLE]

for any $m=0,1,\dots,N_{f}$. Therefore

[TABLE]

The expression for \mathbb{E}\big{[}\left(\Delta\beta_{j}\right)^{2}\big{]} can be used to obtain the integral form of the generalized Itô correction, namely

[TABLE]

Thus, the differential form of the generalized Itô correction reads

[TABLE]

The reason we keep in the final formula the limit $N_{f}\rightarrow\infty$ is that $N_{f}$ depends on the length of the subinterval $\Delta t.$ So, when $\Delta\rightarrow 0,$ the number of frequencies $N_{f}\rightarrow\infty.$ We have provided in the main text a short discussion on the behavior of this limit.

We have computed the autocorrelation and the e-folding time of the noise process $n(t)$ for different values of the parameter $\alpha$ (see Fig. 3).

Appendix B

Let us define

[TABLE]

and

[TABLE]

[TABLE]

This notation allows us to write Eqs. (29)–(32) as Eq. (28) and the analytical solution as Eq. (34).