{'fraction_non_alphanumeric': 0.03432149344876899, 'fraction_numerical': 0.008309821048969969, 'mean_word_length': 4.769007803790412, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improving solution accuracy and convergence for stochastic physics parameterizations with colored noise

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. It 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.


Introduction
State-of-the-art weather, climate, and Earth system models simulate the time evolution of atmospheric motions by solving a set of differential equations. The atmospheric processes (i.e., motions or phenomena) that have significant impact on the large-scale flow motions but cannot be resolved by the computational mesh are described by the physics parameterizations. These parameterizations provide source and sink terms on the right-hand side of the resolved dynamics equations [8].
In recent years, stochastic parameterizations have become an active area of research (see review by [1]). 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.
As shown in [4], the convergence issue can be investigated through the use of tools from stochastic analysis (see e.g. [9,7]). 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), except for special cases (e.g., the second-order Runge-Kutta scheme analyzed by [4]), will no longer converge to the physically relevant solution. Multiple examples relevant for atmospheric modeling can be found in the paper of [4]. Here, "physically relevant solution" refers to the one corresponding to ordinary calculus (see discussion below).
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. The two most popular choices for discretizing the integral 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). For the physical systems the weather and climate researchers are attempting to model, the physically relevant solution alluded to in the previous paragraph is the one corresponding to the Stratonovich integral. This is because when a sufficiently small time period is simulated using sufficiently fine temporal resolution, the evolution of the atmospheric motion is expected to be deterministic and to follow the basic principles of physics. In other words, the future state of the system depends on its current state; the solution of the atmospheric equations is expected to be the Stratonovich solution. This is unlike the evolution of e.g., economy or population, where there exist aspects of the future evolution that cannot be fully determined by the current state of the system. The Itô interpretation of the time integral and the Itô solution of the stochastic equation are more suitable for the latter case.
It is important to note that most of the popular time integration schemes designed for deterministic problems will converge to the Itô solution when applied to stochastic problems driven by white noise [7]. 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. It is therefore plausible that colored noise, which has non-zero autocorrelation length, can provide a better description of such processes.
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.
Another point worth mentioning is that, as [4] and [5] have pointed out, certain deterministic numerical schemes (e.g., the second-order Runge-Kutta scheme) have a 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 an Euler forward scheme complemented by a correction constructed only for the stochastic term, and hence can be cost-effective.
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.
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.

The generalized Itô correction
We consider the following deterministic differential equation 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 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 whereṘ(t) represents a general noise term. The deterministic equation (1) becomes stochastic: Without loss of generality, we assume We focus on how to solve numerically 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.

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 ∆t and denote the increment of R in each bin as ∆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 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 [9]. The two most popular choices are defined as where t * j = t j (left endpoint), and Stratonovich integral: where t * j = (t j + t j+1 ) /2 = t j + ∆t/2 (midpoint) and ∆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 [9].
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. [9]). Below we show that the same is true for colored noise, although the Itô correction needs to be generalized.
For t * j = (t j + t j+1 ) /2, performing a Taylor expansion of g u(t * j ) about t j and expressing ∂u/∂t using Eq. (4) gives where ξ ∈ [t j , (t j + t j+1 ) /2]. For small ∆t, we writė Hence, Eq. (9) can be approximated as Assuming g(u) is sufficiently smooth and ∆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 Eq. (5) becomes The mathematical expectation of the last term in Eq. (12) is The exact form of the expectation depends on the formulation of R. Expression (13) is the generalized Itô correction in its integral form. For example, for the Wiener process we have hence (13) becomes which is the integral form of the traditional Itô correction (see e.g. Section 3.3 in [9]).

Remarks
We want to make three 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 [10]), 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Ṙ(t) independent of the solution u(t). The Itô correction serves as a way to account for the interaction ofṘ(t) and u(t) during the interval ∆t, similar to the role played by memory terms in model reduction which account for the interaction between resolved and unresolved variables.
Third, if we consider as a discretization point for the integral the combination (1 − λ)t j + λt j+1 , where 0 ≤ λ ≤ 1, we obtain a generalization of (13) given by This is the generalization of the correction formula for white noise that is found in [4] (see also Section 3.5 in [7]). Note that the case λ = 0 corresponds to the Itô interpretation and we recover (13). The case λ = 1/2 corresponds to the Stratonovich interpretation and the correction vanishes.

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 with initial condition u(x, 0) = u 0 (x) and periodic boundary conditions on [0, 2π]. 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 is set to 0, we recover the advection-diffusion equation with constant advection velocity discussed in [4]. The inclusion of 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 = 0 and = 10 −3 . Following [4], we let The stochastic noise process n(t) is the same as described in Appendix A of [4] (also described in Appendix A of this paper). For this choice of g(u) and n(t), the generalized Itô correction is given by where N f , C, ω 0 and ω m are parameters of the noise process n(t) (cf. Appendix A).
To derive analytical solutions for the test problem, we express u(x, t) in the form of a superposition of Fourier modes and transform Eq. (17) into a system of stochastic differential equations.
Here i is the imaginary unit. Like in [4], we assume the initial condition contains only one mode, i.e., u(x, 0) = cos(k 0 x) with k 0 = 1.

Case with constant advection velocity
As pointed out by [4], when the advection velocity is constant, the Fourier modes are uncoupled. The ordinary differential equation The analytical solution of Eq. (23) takes the form with A being any complex constant. Initial condition (22) implies that A = 1 in Eq. (24); 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 (17) being real, it is the complex conjugate of the solution for the Fourier mode with k = 1.)

Case with spatially varying advection velocity
For cases with nonzero , 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. (17) linear. We truncate the Fourier series in Eq. (21) to retain only modes with appreciable magnitudes and denote the largest remaining wavenumber as N x . Substituting Eq. (21) for the unknown u in (17) gives By multiplying Eq. (25) with exp(−ikx) and integrating over [0, 2π], we get the following coupled equations: (26) (28) 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 The analytical solution reads

Numerical Results
In this section we use numerical results to show how noise n(t) of different color (roughness in time) can affect the convergence of the simplest explicit scheme, namely the forward Euler scheme. In addition, we demonstrate how the inclusion of the generalized Itô correction can help in restoring and/or accelerating convergence.

Definition of solution error
The error of a numerical solution is evaluated after two time units of integration using the L 2 norm ( [4] and personal communication): Here 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. (24) and (30) is calculated analytically.

Case with constant advection velocity
For the case with = 0, the discretization of Eq. (23) using forward Euler with the Itô correction included is given by 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 , Panel (a) of Figure 1 (which appears also in [4]) 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. 1 We make two observations. First, for the case of white noise (α = 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. [7]). Second, for the case of colored noise (α = 0, Figure 1: Error in the numerical solution of the 1D advection-diffusion equation with constant advection velocity ( = 0 in Eq. 17) and the dependency on time step size (x-axis) and characteristics of the noise term (α = 0, 10 −6 , 10 −5 , 10 −4 , or 1, shown in different colors). The left and right panels show results obtained using the forward Euler scheme without and with the generalized Itô correction, respectively. Simulations were performed for 100 realizations of the noise process and the l 2 solution error was calculated separately for each realization using Eq. (31). The thick dots are the mean error of the 100 realizations; the vertical bars denote the standard deviation around the mean. 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 [3]) when the step size becomes smaller than some critical step size which depends on the color of the noise (value of α). The more red the noise is (larger α), the larger is the critical step size (see also [4] 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 (α = 0, purple line), the Euler scheme with the Itô correction does converge to the analytical solution with order 1/2 (see [6] for an explanation of this convergence rate). We note that this numerical result was mentioned in [4] although not illustrated by any graphic there. Second, for the case of colored noise (α = 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.

Case with spatially varying advection velocity
We continue with the case of a spatially-dependent advection velocity with = 10 −3 . A small value was chosen for 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 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 time steps, becomes very large when is large. For practical purposes (computational cost), we chose a small for the demonstration here.
We discretized Eq. (29) using the Euler scheme with the generalized Itô correction, i.e., where the Itô correction reads with the matrix G G G being The truncation wavenumber N x was chosen empirically: a test simulation was conducted using Eq. (34) with a large N x ; an inspection of the magnitude of the resulting F k 's revealed N x = 5 was sufficient to retain all modes with | F k | > 10 −4 . N x = 5 was then used to obtain the results shown in Figure 2. Figure 2 shows that for the case of white noise (α = 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 = 0 because of the logarithmic scale of the plot. This figure demonstrates that for the case of colored noise (α = 0), the use of the generalized Itô correction again accelerates the establishment of the order 1 convergence regime predicted by deterministic numerical analysis.

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. In addition, as was mentioned in the end of Section 2, we plan to exploit the connection between the Itô correction and memory terms appearing in model reduction formalisms.
where N is the number of discrete time levels per unit time, including the starting and ending time levels. α is a parameter controlling the color of the Fourier spectrum of n(t). To construct different realizations of the noise process, we take values of the coefficients a m and b m from the normal distribution N (0, 1). It should be noted that n(t) is an approximate random noise. The difference between n(t) and the noise termṘ(t) in Section 2 is thatṘ(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Ṙ(t). Let us define This notation allows us to write Eqs. (26)-(28) as Eq. (25) and the analytical solution as Eq. (30).