<figure><img src="image_1.png" /></figure>

HAL
open science

# Numerical discretization causing error variance loss and the need for inflation

Richard Ménard, Sergey Skachko, Olivier Pannekoucke

► To cite this version:

Richard Ménard, Sergey Skachko, Olivier Pannekoucke. Numerical discretization causing error variance loss and the need for inflation. Quarterly Journal of the Royal Meteorological Society, 2021, 147 (740), pp.3498-3520. 10.1002/qj.4139 . hal-04739157

HAL Id: hal-04739157

https://cnrs.hal.science/hal-04739157v1

Submitted on 16 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical discretization causing error variance loss and the need for inflation

Richard Ménard<fnref n="1" /> | Sergey Skachko<fnref n="2" /> | Olivier Pannekoucke<fnref n="3" />,<fnref n="4" />,<fnref n="5" />

**Correspondence**

R. Ménard, Air Quality Research Division,
Environment and Climate Change
Canada, Dorval, Quebec H9P 1J3, Canada
Email: Richard.Ménard@canada.ca

**Funding information**

LEFE INSU, Grant/Award Number:
KAPA grant

**Abstract**

The effects of model discretization errors on the propagation of error covariance have a more complex nature than the effect on the state variable. The analysis is carried out for the advection transport equation, where the continuous (in space and time) propagation of the related error covariance function can be written, solved and compared with the discrete model applied to the covariance matrix. The numerical analysis of the problem is carried out with a 1D-problem, but is also illustrated with a 3D chemical transport model (CTM) used for chemical data assimilation. It is shown that variance loss (compared to the continuous propagation solution) depends on the covariance function itself as well as the numerical discretization scheme. The variance loss is particularly sensitive to the correlation length. In a simple first-order discretization, an analytical expression is obtained and is used to derive an analytical expression for inflation. Experiments show, for example, that following an analysis over a dense network of observations (with spatially uncorrelated errors) a significant variance loss occurs in the propagation step. With the (variance) inflation scheme, we are able to restore the variance lost at each grid point, and at each timestep, during the entire integration. The variance inflation scheme applied to an EnKF can be formulated to change the variance spread of the ensemble or to act directly on the state.

**KEYWORDS**

data assimilation, inflation, loss of variance, advection transport

## 1 | INTRODUCTION

A key characteristic of ensemble Kalman filtering (EnKF) is to produce realistic and accurate forecast-error covariance (Evensen, 2003). However, practical implementations of EnKF result in underestimated ensemble spread when using real observations. To deal with this problem, the forecast-error variance is augmented either by inflating

the ensemble spread by scaling or adding perturbations to ensemble members or by adding stochastic realizations of a forecast-error covariance (Houtekamer and Zhang 2016 give a review). Sources of error variance loss in ensemble data assimilation (DA) may be related to both analysis and forecast stages. In the analysis stage, the limited ensemble sampling (Whitaker and Hamill, 2002; Sacher and Bartello, 2008; Bocquet, 2011) and ‘inbreeding’ or lack of

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

© 2021 Her Majesty the Queen in Right of Canada. *Quarterly Journal of the Royal Meteorological Society* published by John Wiley & Sons Ltd on behalf of Royal Meteorological Society. Reproduced with the permission of the Minister of Environment and Climate Change Canada.

<footnote n="1">¹Air Quality Research Division,
Environment and Climate Change
Canada, Dorval, Quebec, Canada</footnote>

<footnote n="2">²Meteorological Research Division,
Environment and Climate Change
Canada, Dorval, Quebec, Canada</footnote>

<footnote n="3">³INPT-ENM, Toulouse, France</footnote>

<footnote n="4">⁴CNRM, Université de Toulouse,
Météo-France, CNRS, Toulouse, France</footnote>

<footnote n="5">⁵CERFACS, Toulouse, France</footnote>
<figure><img src="image_3.png" /></figure>

cross-validation (Houtekamer and Mitchell, 1998; Mitchell and Houtekamer, 2009) are responsible for the error variance loss. The forecast stage has three issues:

(a) unaccounted random model error due to missing physical processes (e.g., Meng and Zhang 2007; Buizza et al., 1999),

(b) subgrid-scale upscale propagation (SKEB; Berner et al., 2009), and

(c) numerical model discretization error effect in the propagation of the forecast-error covariance in space and time.

Although discretization error on the model state is usually investigated, little information about its impact on error covariance is known. The aim of this article is to reveal this relationship, and to explain how the contribution of the discretization differs between the model state and the forecast-error covariance.

In the absence of an analytical solution for the propagation of the error covariance, the effect of discretizing the model on the error covariance can only be perceived indirectly. One way is to perform DA using different spatial resolutions. For example Hamill and Whitaker (2005) conducted a simulation using a spectral T31, two-layer model with T127 resolution for the truth, and showed that an additive inflation based on differences between T63 and T31 reduces the model error variance by about 40%, as well as reducing the correlation length. In another experiment, using the ECMWF hybrid DA system, Bonavita et al. (2016) showed that increasing the model resolution of the ensemble DA component from T399 (≈50 km) to T639 (≈30 km) increased the 4D-Var error variance, and thus required significantly less inflation. It was argued that the increase in variance was due to the increase of dynamical variability allowed by higher-resolution models. However, the increase can also be associated with smaller model discretization error on the error covariance which, as we shall see in this study, leads to a smaller variance loss. Some studies have also examined the effect of reduced arithmetic precision (Hatfield et al., 2018) on DA with a general circulation model or the effect of numerical discretization error on ordinary differential equations in solving the Lorenz-96 model (Grudzien et al., 2020). It is worth noting that Houtekamer et al. (2014) also concluded that increasing the horizontal resolution of the ensemble members in the operational ensemble Kalman filter (EnKF) at the Canadian Meteorological Center significantly improved the performance of the DA system. However, since a number of modifications were conducted simultaneously in their experiment, we cannot attribute the effect solely to the increased resolution.

It is also known that realistic DA systems require larger inflation over densely observed regions (Zhang et al., 2004; Miyoshi, 2011; Whitaker and Hamill, 2012).

Accounting for the fact that, over these regions, the correlation length is usually smallest, it could be an indication that variance loss may also depend (inversely) on the correlation length-scale. In particular, Miyoshi (2011) developed an adaptive inflation scheme used in simulating the radiosondes network, where they showed that an inflation of 1.05–1.25 was needed over continents (highest values over densely observed areas) while over the ocean the optimal inflation was neutral, with some regions even requiring deflation.

However, a direct and detailed analysis of the effect of spatial discretization errors on the predictability error covariance is possible for the (pure) advection transport model (Lin and Rood, 1996; Lyster et al., 2004), where we can establish continuous dynamics of forecast-error covariance both spatially and temporally and obtain solutions. The formalism was introduced in the seminal paper of Cohn (1993) after Roger Daley (personal communication, 1993) noted an interesting error covariance property in the context of a simple 1D advection problem. In the continuum (space and time) representation, covariance matrices are continuous covariance functions of a pair of spatial coordinates and time. The time integration of a covariance function can be split into separate terms accounting for the dynamics with respect to each spatial coordinate. In this case there are two operators that need to be integrated; this is accomplished by using an operator splitting approach (Strang, 1968; Cohn, 1993; Beljaars et al., 2018). In the case of the advection problem, a solution of the original equation without splitting can also be obtained by using the method of characteristics (Cohn, 1993). The original equation with no model error has the important property that the error variance simply obeys the advection equation.

One of the first implementations of a KF scheme using satellite observations was achieved for the assimilation of long-lived chemical species in the stratosphere, where the flow is nearly two-dimensional and generally follows isentropic surfaces (Ménard et al., 2000). Two configurations of the assimilation system were possible. One configuration used the standard KF. In the other, the forecast-error variance of the KF (i.e., between analysis times), is forced to obey the advection dynamics while simultaneously imposing conservation of total covariance (i.e., $d[P]/dt = 0$ where $[P]$ is the Frobenius norm of the error covariance) as the continuous theory of error covariance evolution would suggest. In the case of pure transport, that is, with no assimilation of observations, differences were particularly striking in regions of shear flow. We should emphasize that these experiments were not performed with an ensemble formulation and there are thus no issues with sampling, localisation, or non-Gaussian effects since the model equation is linear. In the case of assimilation of
observations, the assimilation system was equipped with an estimation of a residual climatological model error growth, by assuring there is no time tendency in the $\chi^2$ metric of the innovation-covariance fit. The results showed that, when observations are assimilated with the standard KF, the residual model error growth is two times larger than the model error growth obtained with the continuous variance evolution Ménard and Chang (2000). The origin of the lack of efficiency in the case of the standard KF was not clearly understood at that time, but the use of the continuous covariance evolution equation leads to the development of the so-called sequential filter (Khattatov et al., 2000) and to the Lagrangian filter (Lyster et al., 2004).

Recently Pannekoucke et al. (2016) further extended this approach by developing evolution equations of error covariance parameters (a method called the parametric KF), specifically the error variance and the metric tensor (related to the correlation length). It was applied to a 1D advection-diffusion problem and also to the nonlinear Burger's advection equation (Pannekoucke et al., 2018). Using this parametric approach, a method to compute the model error covariance (growth) due to numerical discretization was also developed (Pannekoucke et al., 2021).

The organization of the article is as follows. First, we present the theoretical aspects. We recall the continuous and discrete representations for the linear advection transport model and discuss how spatial splitting can yield solutions for both formulations. Then, using a first-order discretization, we obtain an analytical expression for the loss of variance. Following Cohn (1993), a local covariance coordinate is introduced such that the transport is described along the main diagonal and main anti-diagonal of the covariance matrix. The numerical results of variance loss is examined for different covariance propagation formulations (Section 3), for several numerical schemes (Section 4), and for a state-of-the-art 3D chemical transport model (Section 5). Finally in Section 6, we develop and test a variance inflation scheme, both deterministic (based on spread and correlation length) and stochastic. In particular, we examine how to restore the error variance after an analysis generated from a dense observation network over a limited region. Conclusions then follow.

## 2 | THEORY

The continuous and discrete formulation of the error covariance evolution for the linear advection equation is outlined in this section. Some theoretical derivations and numerical details are given in the Appendices. We begin this section by first recalling important properties in continuous space derived by Cohn (1993). We then discuss how spatial splitting leads to the standard KF form of the

predictability error covariance. We recall the analytical solution of Cohn (1993) based on the method of characteristics which we should point out does not involve spatial splitting. By comparing the continuous and discrete solutions, we derive an analytic expression of the variance loss for a simple numerical scheme. We end this section by presenting an alternate spatial covariance space introduced by Cohn (1993) and we show how it can be implemented numerically.

### 2.1 | Background on the continuous and discrete error covariance evolution

Let us consider as a dynamical model the linear advection equation. In 1D it has the form

$$ \frac{\partial c}{\partial t} + u(x, t) \frac{\partial c}{\partial x} = 0. \qquad (1) $$

Since this equation also represents the transport of a tracer, c can also be the tracer concentration expressed as a mixing ratio and $u(x, t)$ the wind.

Assuming the errors are due only to initial tracer concentrations (i.e., no wind error), the state error $\epsilon$ then obeys the linear partial differential equation (PDE) of the form

$$ \frac{\partial \epsilon}{\partial t} + u(x, t) \frac{\partial \epsilon}{\partial x} = 0. \qquad (2) $$

Here the error is defined with respect to the truth, which is assumed to obey the same linear advection equation. The error $\epsilon(x, t)$ is a scalar-valued function of the spatial and temporal continuum. Because the error evolution Equation (2) is linear, the continuous evolution equation of the covariance function has a closed form (i.e., without involving higher moments).

The error covariance function $P(x_1, x_2, t)$ at time $t$ is the covariance of the tracer concentration error at a pair of points, one with coordinate $x_1$ and the other at the coordinate $x_2$,

$$ P(x_1, x_2, t) = E[\epsilon(x_1, t) \epsilon(x_2, t)], \qquad (3) $$

where $E$ is the expectation operator. The covariance function is defined over a tensor product of the spatial space (that is a $2n$-dimensional space if $n$ is the dimension of the spatial space). In most examples considered here, the spatial space is one-dimensional, so the covariance evolution equation will be in a two-dimensional spatial space. As in Cohn (1993), we first introduce an extended (i.e., tensor product) spatial space $(x_1, x_2)$ where the error $\epsilon_1 = \epsilon_1(x_1, t)$ obeys the advection Equation (2) with respect to the spatial coordinate $x_1$ and the error $\epsilon_2 = \epsilon_2(x_2, t)$ obeys the advection equation with respect to the spatial coordinate $x_2$.
<figure><img src="image_4.png" /></figure>

To obtain the evolution equation for the covariance function, we multiply the evolution error equation with respect to the coordinate $x_1$ with the error $\epsilon_2 = \epsilon_2(x_2, t)$ and similarly, for the error evolution equation with respect to the coordinate $x_2$ multiplied by $\epsilon_1(x_1, t)$. Adding these resulting equations and taking the expectation we get

$$ \frac{\partial P}{\partial t} + L_1(P) + L_2(P) = 0, \quad (4) $$

where $P = P(x_1, x_2, t)$ is the (error) covariance function, and $L_1, L_2$ are linear operators defined as

$$ L_i = u(x_i, t) \frac{\partial}{\partial x_i}, \quad (i = 1, 2), \quad (5) $$

which are applied on the covariance function $P$. Cohn (1993) gives further details and generalizations for higher-dimensional space and for several variables.

### 2.1.1 | Solution by operator splitting

The typical solution of Equation (4) is obtained by operator splitting (Strang, 1968). Since the operators $L_1, L_2$ commute, that is, for any function $f(x_1, x_2, t)$ we have $L_1 L_2 f = L_2 L_1 f$, the solution of Equation (4) can be obtained (without loss of accuracy) by solving in sequence

$$ \frac{\partial P}{\partial t} + L_1(P) = 0, \quad (6) $$

$$ \frac{\partial P}{\partial t} + L_2(P) = 0. \quad (7) $$

A proof and an example are given in Appendix A. In a discrete formulation, the initial error is a vector $\epsilon(0)$ that is propagated in time as $\epsilon(t) = M_{t \leftarrow 0} \epsilon(0)$, where $M_{t \leftarrow 0}$ is the model matrix from time 0 to time $t$. When we represent the initial error covariance as a matrix, that is, $P(0) = E[\epsilon(0) \epsilon^T(0)]$, the elements of the columns of the matrix correspond to a position in $x_1$ and the rows elements to a position in $x_2$. Thus solving Equation (6) in a discretized space takes the form

$$ P^*(t) = M P(0). \quad (8) $$

Similarly, solving Equation (7) using as initial condition the result of Equation (6), that is, Equation (8), actually takes the form

$$ P(t) = M [P^*(t)]^T, \quad (9) $$

since we are now operating on the $x_2$ spatial coordinate. Combining Equation (8) with Equation (9), we then get

$$ P(t) = M[M P(0)]^T = M P(0) M^T, \quad (10) $$

the standard KF form for a discrete space and discrete model. A step-by-step derivation of this result is given in Appendix A.

### 2.1.2 | Solution with the method of characteristics

Let $X$ represent the position of a fluid particle at time $t = 0$, that can also be used to label the particle. As the particle moves with the flow, its position $x$ at a later time is given as

$$ x(t; X) = \Phi_t(X). \quad (11) $$

The flow $\Phi_t$ or trajectory is obtained by solving the characteristics equation,

$$ \frac{dx}{dt} = u(x, t), \quad (12) $$

whose solution is given as

$$ x(t; X) = \int_{0}^{t} u[x(\tau), \tau] d\tau, \quad (13) $$

where $u$ is the wind. We recall that the initial particle position given as $x(0; X) = X$.

Equation (4) is a hyperbolic equation for a pair of fluid particles ($X_1, X_2$) with trajectories

$$ x_1(t; X_1) = \Phi_t(X_1), \quad x_2(t; X_2) = \Phi_t(X_2), \quad (14) $$

which are solutions of the characteristics

$$ \frac{dx_1}{dt} = u(x_1, t), \quad \frac{dx_2}{dt} = u(x_2, t). \quad (15) $$

Note that the trajectories are driven by the same wind field but with different initial positions.

Equation (4) then says that the covariance function between a pair of fluid particles, that is,

$$ P\{x_1(t; X_1), x_2(t; X_2), t\} = \tilde{P}(t; X_1, X_2), \quad (16) $$

satisfies the conservation equation

$$ \frac{d\tilde{P}}{dt} = 0. \quad (17) $$

In other words, the covariance function between a pair of fluid particles is conserved along the flow.

An important property arises for the error variance. Consider the case where the position of the two fluid particles coincides initially, that is, $x_1(0) = x_2(0)$. Then according to Equation (15), the two positions will remain
coincident at any later time, $x_1(t) = x_2(t)$ (the characteristic line $x_1 = x_2$ is an invariant manifold) since in this case the winds in Equation (15) are identical along both trajectories. Because an error covariance evaluated at the same coordinate $x_1 = x_2 = x$ is actually an error variance, which we denote by

$$ \tilde{P}(x, x) = V(x), \quad (18) $$

we conclude that error variance is conserved and obeys the advection equation,

$$ \frac{dV}{dt} = \frac{\partial V}{\partial t} + u(x, t) \frac{\partial V}{\partial x} = 0. \quad (19) $$

These properties were first derived by Cohn (1993).

To obtain a numerical solution of the error covariance propagation Equation (17), we simply use a Lagrangian formulation. The characteristics (Equations (15)) can be solved by using a Runge-Kutta method or, as we did in this study, using a simple first-order trajectory (as long as the total time integration remains short enough to have an accurate result, typically less than a few days). We should note, in passing, that a pure Lagrangian DA has also been developed on these principles (Lyster et al., 2004), where it is found that it takes several days of integration before clustering of the material points becomes an issue; this then requires a remapping of the error covariance.

## 2.2 | An analytical expression for variance loss

An analytical expression of variance loss for the (first-order) semi-Lagrangian discretization scheme with linear interpolation is derived in this section. To simplify we assume that the Courant number value is between 0 and 1. The (first-order) upstream displacement is simply computed as $\delta x = u_i^n \Delta t$, where $u$ is the wind at grid point $i$ and timestep $n$. The discretized tracer concentration $c$ then evolves as

$$ \begin{aligned} c_{i}^{n+1} &= c_{i}^{n} + \frac{u_{i}^{n} \Delta t}{\Delta x} [c_{i-1}^{n} - c_{i}^{n}] \\ &= (1 - \alpha_{i}^{n}) c_{i}^{n} + \alpha_{i}^{n} c_{i-1}^{n}, \end{aligned} \quad (20) $$

where $\alpha_i^n = u_i^n \Delta t / \Delta x$ is the local Courant number. In the standard KF, we apply the advection first on the $x_1$ spatial coordinate and then, using operator splitting, on the $x_2$ coordinate. Let us denote by $P^n(i, j)$ the error covariance at time $n$ for the pair of grid points $(x_1(i), x_2(j))$. Applying the transport on the column of the matrix $P$ is equivalent to applying the transport operator on the coordinate $x_1$ only, that is for the first part of the operator splitting scheme

$$ P^*(i, j) = (1 - \alpha_i^n)P^n(i, j) + \alpha_j^n P^n(i - 1, j). \quad (21) $$

Then, applying the transport of $P^*$ along $x_2$, consistent with the second step of operator splitting, we get

$$ P^{n+1}(i, j) = (1 - \alpha_j^n)P^*(i, j) + \alpha_j^n P^*(i, j - 1). \quad (22) $$

Combining Equation (21) with Equation (22) gives

$$ P^{n+1}(i, j) = (1 - \alpha_j^n)[(1 - \alpha_i^n)P^n(i, j) + \alpha_i^n P^n(i - 1, j)] \\ + \alpha_j^n[(1 - \alpha_i^n)P^n(i, j - 1) + \alpha_i^n P^n(i - 1, j - 1)]. $$

For the error variance $P^n(i, i) = V^n(i)$ we get, in particular,

$$ \begin{aligned} P^{n+1}(i, i) = P^n(i, i) &- \alpha_i^n P^n(i, i) - \alpha_i^n P^n(i, i) \\ &+ \alpha_i^n \alpha_i^n P^n(i, i) \\ &+ \alpha_i^n (1 - \alpha_i^n) P^n(i - 1, i) \\ &+ \alpha_i^n (1 - \alpha_i^n) P^n(i, i - 1) \\ &+ \alpha_i^n \alpha_i^n P^n(i - 1, j - 1) \end{aligned} $$

which, if rewritten as follows,

$$ \begin{aligned} P^{n+1}(i, i) = P^n(i, i) &- \alpha_i^n P^n(i, i) - \alpha_i^n P^n(i - 1, i - 1) \\ &- \alpha_i^n P^n(i, i) + \alpha_i^n \alpha_i^n P^n(i, i) \\ &+ \alpha_i^n (1 - \alpha_i^n) P^n(i - 1, i) \\ &+ \alpha_i^n (1 - \alpha_i^n) P^n(i, i - 1) \\ &- \alpha_i^n P^n(i - 1, i - 1) + \alpha_i^n \alpha_i^n P^n(i - 1, j - 1), \end{aligned} $$

yields the following expression

$$ V^{n+1}(i) = (1 - \alpha_i^n)V^n(i) + \alpha_i^n V^n(i - 1) + \text{HOT}, \quad (23) $$

where the High-Order Terms (HOT) are given by

$$ \begin{aligned} \text{HOT} &= L_n^D(i) \\ &= -\alpha_i^n(1 - \alpha_i^n) \\ & \quad \times \{V^n(i) - P^n(i, i - 1) - P^n(i - 1, i) + V^n(i - 1)\}. \end{aligned} \quad (24) $$

For the purpose of the discussion let us ignore the HOT just for now. Equation (23) expresses the discrete form of the advection of error variance, as suggested by the error variance Equation (19). The difference between a discrete solution of the continuous propagation Equation (4) and the solution obtained by operator splitting (as we just did above), is the HOT.

The HOT together produce a variance loss which we will denote by $L_n^D$. Indeed, the quantity in braces is generally positive since usually $P(i, i - 1) < V(i)$ and
<figure><img src="image_5.png" /></figure>

$P(i-1, i) < V(i-1)$ (for error variances that do not vary too strongly in space). Also, the term $\alpha(1-\alpha)$ varies between 0 to a maximum value of 1/4 for $\alpha = 1/2$. Thus, the HOT or $L_n^D(i)$ is negative, so that the predicted error variance is reduced compared to the transport of variance, that is, Equation (23) without this HOT term.

We should also note that the terms in braces depend on the correlation length. Indeed, since $P(i, i-1) = \sigma(i)\sigma(i-1)C(i, i-1)$ where $\sigma$ is the error standard deviation and $C$ is the correlation, evaluating $C(i, i-1)$ using a Taylor series expansion evaluated at $(i, i)$ gives

$$V(i) - P(i, i-1) = \sigma^2(i) - \sigma(i)\sigma(i-1)C(i, i-1) \\ \approx \sigma^2(i) \left( \frac{\Delta x^2}{2L_c^2} \right), \qquad (25)$$

where we used that fact that the second derivative of the correlation is equal to $-1/(L_c)^2$, where $L_c$ is the correlation length. Note that there is another term in Equation (24) that is similar to Equation (25). Thus, in the case when $\sigma(i) = \sigma(i-1)$ the term in braces in Equation (24) becomes equal to $\sigma^2(i)(\Delta x^2/L_c^2)$, that is, twice the quantity in Equation (25). In summary, the variance loss is larger when the ratio of the grid size to the correlation length is large. That happens, for example, after assimilating observations with spatially uncorrelated errors, where $L_c$ reaches very small values. The variance loss is also a multiplicative error that depends on the error variance. In the extreme case where the correlation length is equal to the grid size, that is, $\Delta x = L_c$, the loss of variance is maximized and the term in braces in Equation (24) is then equal to $\sigma^2(i)$.

## 2.3 | Forecast covariance in local covariance coordinates

There is another approach to solve numerically the continuous covariance evolution equation for the advection equation which involves a change of the coordinate system of the covariance propagation (Equation (4)). This change of coordinate was introduced by Cohn (1993) and could be called the local covariance space. What is interesting about this new formulation is that, with the same numerical discretization scheme used in the standard operator splitting approach of the KF, the propagation of error covariance in the local covariance coordinates produces significantly different numerical results. Thus, the accuracy of the transport of error variance is not just dependent on the numerical scheme used for the state, but depends also on the covariance coordinates (e.g., Lagrangian, local, and standard splitting) and the corresponding scheme to evolve the covariance function.

The local covariance coordinates $(\chi, \xi)$ are obtained from a transformation of the physical spatial coordinates $(x_1, x_2)$ as

$$\chi = \frac{1}{2}(x_1 + x_2), \quad \xi = \frac{1}{2}(x_1 - x_2). \qquad (26)$$

$\chi$ can be regarded as the centred position of a pair of points and $\xi$ as a separation distance measured between a pair of points. In a covariance matrix (1D univariate) the $\chi$ coordinate corresponds to positions parallel to the diagonal, and the $\xi$ coordinate to positions perpendicular to the diagonal. This coordinate system is implicitly invoked when we describe the error covariance matrix in terms of variance (along $\xi = 0$) and correlation length (related to the second derivative with respect to $\xi$ of the error correlation).

Following Cohn (1993), the covariance function can be rewritten in the $(\chi, \xi)$ local coordinates,

$$\tilde{P}(\chi, \xi, t) = P[x_1(\chi, \xi), x_2(\chi, \xi), t] \\ = P(\chi + \xi, \chi - \xi, t). \qquad (27)$$

It evolves in the local coordinates $(\chi, \xi)$ as

$$\frac{\partial \tilde{P}}{\partial t} + L_{\chi}(\tilde{P}) + L_{\xi}(\tilde{P}) = 0, \qquad (28)$$

where

$$\tilde{L}_{\chi} = \frac{1}{2}[u(\chi + \xi, t) + u(\chi - \xi, t)] \frac{\partial}{\partial \chi}, \qquad (29)$$

$$\tilde{L}_{\xi} = \frac{1}{2}[u(\chi + \xi, t) - u(\chi - \xi, t)] \frac{\partial}{\partial \xi}. \qquad (30)$$

As the reader can verify, the operators $\tilde{L}_{\chi}$ and $\tilde{L}_{\xi}$ commute, and thus the solution using operator splitting

$$\frac{\partial \tilde{P}}{\partial t} + \tilde{L}_{\chi}(\tilde{P}) = 0, \qquad (31)$$

$$\frac{\partial \tilde{P}}{\partial t} + \tilde{L}_{\xi}(\tilde{P}) = 0, \qquad (32)$$

is exact. This set of equations leads naturally to the error variance conservation (Equation (19)). Indeed, letting $\xi = 0$, we have $\tilde{L}_{\chi} = u(\chi, t)\partial_{\chi}$ and $\tilde{L}_{\xi} = 0$. Thus the evolution of $\tilde{P}(\chi, \xi = 0, t) = V(\chi, t)$ obey the advection equation.

The spatial discretization in this local covariance space is somewhat non-trivial and requires attention to details. The interested reader can consult Appendix B for details. However, it is not necessary to know these details in order to follow the experiments presented in the next section.
<figure><img src="image_6.png" /></figure>

Before we present the results of the experiments, we should note that for small $\xi$ and for a smooth wind field, $[u(\chi + \xi, t) - u(\chi - \xi, t)]$ becomes small compared to $[u(\chi + \xi, t) + u(\chi - \xi, t)]$, so that the 'advection' term along the $\xi$ coordinate or parallel to the main diagonal, that is Equation (32), produces a small increment. In fact, in the case of no shear flow, that is, $[u(\chi + \xi, t) - u(\chi - \xi, t)] = 0$ there is conservation of $\tilde{P}$ in the $\xi$ direction, that is, parallel to the main diagonal. Illustration of this will be presented in the next section.

## 3 | VARIANCE LOSS FOR DIFFERENT COVARIANCE PROPAGATION FORMULATIONS

As we have developed in the previous section, there are more than one numerical approaches that can be used to propagate the error covariance for the advection dynamics. Numerical experiments using the standard KF form (Equation (10)) and the local covariance coordinates form (Equation (31)), both using the same Eulerian model discretization, are compared with the Lagrangian algorithm, Section 2.1.2. The Lagrangian algorithm has no spatial splitting and small discretization errors, and can be considered to be close to the true propagation of the error covariance.

The Eulerian model used here is a first-order semi-Lagrangian scheme with linear interpolation (other schemes are considered in Section 4). The domain is periodic of length $L = 10,000$ km and spatially discretized with $N = 101$ grid points, which gives $\Delta x = 99.01$ km or with a Fourier representation of a total of 50 wavenumbers. The transport is driven by a non-uniform but stationary wind field having the shape $u(x) = -U_{\max} \sin(2\pi x/L)$, so that the maximum convergence is at $x = 0$. The amplitude of the wind sinusoid, $U_{\max}$, was chosen to be $16.0433 \text{ m} \cdot \text{s}^{-1}$ so we can see the effect of wind shear over a time period of 24 hr. The timestep was chosen such that the maximum Courant number is 0.2917 giving $\Delta t = 0.5$ hr based on the $U_{\max}$. The initial error covariance is a homogeneous SOAR (second-order autoregressive) covariance model with uniform variance of 10 with a correlation length-scale parameter of $p_{\text{SOAR}} = 296.4$ km, which corresponds to a correlation length of $L_{\text{rmc}} \approx 400$ km, that is four times the grid size $\approx 4\Delta x$ (Ménard et al., 2016 or Section 4 for discussions on the relationship between the length-scale parameter and the correlation length).

The results after a 48 timesteps ($\approx 24$ hr of integration) are illustrated in Figure 1 where we compare (a) the physical space propagation (KF), (b) the local covariance space propagation, and (c) the propagation along characteristics, that is, Lagrangian formulation. The covariances

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 1** Comparison of error covariance propagation using (a) physical space KF, and (b) local covariance space. (c) shows the covariance propagation using the characteristics or Lagrangian model</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIGURE 2 (a) The error variance and (b) the correlation near $x = 0$ along a constant $\xi$ for KF formulation (green), the local covariance propagation (blue), and the Lagrangian formulation (red)</figcaption></figure>

displayed in Figure 1 have been normalized by the initial (uniform) variance in all three panels. In theory, with the continuous formulation, a uniform initial variance should have a forecast variance that remains uniform at later times. However, with the standard KF experiment, the normalized error variance varies from 0.6 to 1. Note that the normalized forecast-error variance remains equal to one at $x = 0$, which is not a numerical artefact but a result from having $u(0) = 0$, resulting in $c^{n+1} = c^n$. In (b), the variance remains constant and uniform with the local covariance coordinate formulation while using the same discretization scheme as in the KF. The variance is also nearly uniform using the Lagrangian formulation.

We should remark that loss of error variance has already been reported in past studies. Loss of variance was observed in global 2D standard KF using a state-of-the-art transport model (Ménard et al., 2000) for the assimilation of long-lived species in the stratosphere (on isentropic coordinates).

Figure 2 shows the error variance of the three formulations corresponding to Figure 1, and the correlation near $x/L = -0.05$ where the variance loss is maximum. We note that close to $x/L = -0.05$ the correlation length in the local covariance space scheme is nearly identical to the length-scale given by the Lagrangian scheme, while the physical space (KF) has increased the correlation length. We also note that the variance loss is always associated with an increase of correlation length.

We have also compared the error covariance eigenvalues in the physical spatial space (KF) and with the local covariance space. The result after only one timestep integration shows that the eigenvalues for the large scales

remain the same, but for the smallest resolved scales is decreased by nearly a factor of 10 in the physical space (KF) formulation (results not shown).

These experiments indicate that in a KF, shorter correlation length gives rise to larger variance loss. However, the local covariance space coordinate formulation produces (in these experiments) no variance loss and we observe only slight increase of correlation length compared to the Lagrangian solution. The standard KF, which integrates a field (slice of a covariance along columns and rows) that is peaked around the diagonal and thus integrates a field that varies considerably over a short distance. However, the local covariance coordinate scheme integrates slices parallel to the diagonal which are lot smoother than across an $x_1$ or $x_2$ slice, followed by a wind shear advection along anti-diagonal; Equations (30) and (32). The integration parallel to the anti-diagonal does not create large numerical errors, since over a short $\xi$ distance the difference of winds $u(\chi + \xi, t) - u(\chi - \xi, t)$ is small, and so the discretization error along the $\xi$ directions remain small. Overall, there is less numerical diffusion effectively used in the local covariance coordinate formulation.

## 4 | VARIANCE LOSS FOR SEVERAL NUMERICAL SCHEMES

Let us now turn our investigation to different numerical model discretization schemes and their associated loss of variance. Again we consider a 1D advection equation on a periodic domain. We examine four basic discretizations of two-time-level schemes. In flux form, the donor cell (first order) and the van Leer scheme (second order as in Daley (1995)), and semi-Lagrangian schemes, of first-order with linear interpolation (as used previously; Bates and McDonald (1982)) and with the cubic spline interpolation (Riijsbøgaard et al., 1998). We will also explore a dispersive scheme; the spectral advection scheme with a fourth-order dissipation, i.e., $\nabla^4$.

First we should realize that, although we are solving a linear advection problem, some of these discretized schemes are actually nonlinear (i.e., van Leer, semi-Lagrangian cubic spline). To evolve the error covariance, we have adopted an ensemble KF approach with a large number of ensembles. In fact, and as we will now show, the continuous propagation Equation (4) can actually be derived by using a finite number of ensemble members, and not by just using the expectation operator. Indeed, Equation (2) can be rewritten for each member of an ensemble of errors $\{\epsilon^{[n]}(x, t) \ n = 1, \dots, N\}$ and sum up over all ensemble members, giving
<figure><img src="image_11.png" /></figure>

$$ \frac{\partial}{\partial t} \sum_n (\varepsilon^{[n]}(x_1, t) \varepsilon^{[n]}(x_2, t)) \\ + u(x_1, t) \frac{\partial}{\partial x_1} (\varepsilon^{[n]}(x_1, t) \varepsilon^{[n]}(x_2, t)) \\ + u(x_2, t) \frac{\partial}{\partial x_2} (\varepsilon^{[n]}(x_1, t) \varepsilon^{[n]}(x_2, t)) = 0, \quad (33) $$

that is

$$ \frac{\partial P^e}{\partial t} + u(x_1, t) \frac{\partial P^e}{\partial x_1} + u(x_2, t) \frac{\partial P^e}{\partial x_2} = 0, \quad (34) $$

where $P^e$ is the ensemble covariance.

The parameters of these 1D experiments were chosen to be similar to the zonal conditions of the 3D-transport experiment which will be presented in the next section. Considering the length of the periodic domain $L = 40,000$ km equal to the circumference of the Earth, and a correlation length $L_c = 800$ km as used in the 3D transport model, we have $L_c/L = 0.02$. With a grid size of $3.75^\circ$ in the zonal direction, the whole domain thus contains 96 grid boxes. This corresponds to a wave number, in the zonal direction, of $N = 48$ such that the number of grid points $= 2N + 1 = 97$. We choose a timestep of 1 hr (someWHAT comparable to the timesteps used in Section 5), and we make sure that the winds amplitude is such that the (average) Courant number is around 0.1.

To simulate more realistic winds, we use a realization of a stochastic process. Specifically, we have chosen a spatially correlated model based on a First-Order Auto-Regressive (FOAR) model with a correlation length equal to the initial concentration correlation length used in the simulations of the previous section. For error correlation models, it is important to realize that the correlation length and the length-scale parameter in a correlation model do not match. As proposed by Daley (1991) and Pannekoucke et al. (2008), the correlation length $L_c$ can be defined based on the curvature of the correlation function at the origin, that is, $L_c^2 = -1/C''(0)$. For a Gaussian correlation model, the length-scale parameter $p_G$ is identical to the correlation length, $L_c = p_G$. But for a Second-Order Auto-Regressive (SOAR) correlation model we have $L_c = 1.3494 p_{SOAR}$, and for a compact support Gaspari and Cohn correlation function we have $L_c = 0.5756 p_{GC}$. For non-differentiable functions at the origin, we can define the correlation length as the distance for which the correlation reaches a value of $e^{-1/2} \approx 0.606$ based on the Gaussian correlation model. Thus for a FOAR correlation model we then have $L_c = 0.5005 p_{FOAR}$ (Ménard et al., 2016). With these definitions, $p_{FOAR} = 2.72 p_{SOAR}$.

The results of a 48 hr integration are shown in Figure 3. We note that the first-order schemes, that is, the donor cell and the semi-Lagrangian with linear interpolation, are the most diffusive and lose variance the most. In fact the numerics of both schemes in 1D are identical. About

<figure><img src="image_12.png" /><figcaption>FIGURE 3 EnKF error variance as a function of time for several numerical schemes: flux-form schemes (first-order red solid; second-order van Leer red dashed), semi-Lagrangian schemes (linear interpolation pale blue circles; cubic spline interpolation blue solid), spectral schemes with $\nabla^4$ diffusion (same diffusion coefficient as in other schemes black solid; with augmented diffusion black dotted)</figcaption></figure>

30% of the initial variance is lost in 2 days. The van Leer (second-order) scheme and semi-Lagrangian with cubic spline interpolation perform better, having less dissipation and conserve more closely the error variance. Finally, the spectral model with non-uniform winds is usually dispersive and we used a diffusion term of order 4 to selectively dissipate the smallest scales. In solid black we use the same diffusion coefficient as in the semi-Lagrangian with cubic spline interpolation (required when the winds are not uniform), and in dotted black we augment the diffusion coefficient by several orders of magnitude in order to reach a nearly neutral loss of error variance.

These experiments show that the error variance is generally not conserved, but varies between schemes. In the next section we will examine the loss of variance in a 3D transport model that uses a third-order scheme, that is, the ppm (piece-wise parabolic method), known to yield realistic transport simulations.

## 5 | VARIANCE LOSS IN A 3D CHEMICAL TRANSPORT MODEL

Ensemble model experiments were conducted using a state-of-the art chemical transport model (CTM) used for stratospheric chemical assimilation (Errera et al., 2008; Skachko et al., 2014; Skachko et al., 2016). The numerical scheme is a flux-form semi-Lagrangian model using a third-order PPM (Piecewise Parabolic Method) for
<figure><img src="image_13.png" /></figure>

FIGURE 4 Comparison of error correlation for different ensemble sizes: (a) 20 members, (b) 100 members, and (c) 500 members. The reference point has been chosen in the surf-zone during a Rossby wave breaking event displaying very strong anisotropy. No localization has been applied

<figure><img src="image_14.png" /></figure>

subgrid-scale representation (Lin and Rood, 1996). For the purposes of this study, the chemistry is turned off and only the advection of passive tracer is considered. The CTM is driven by winds and temperatures obtained from the ECMWF ERA-Interim reanalysis (Dee *et al.*, 2011). The horizontal resolution of the model grid is 3.75° longitude by 2.5° latitude. Vertically, the model uses a subset of 37 levels of the (60-level) ERA-Interim reanalysis from the surface up to 0.1 hPa. The model timestep was set to 30 min.

Here we consider the propagation of the error covariance with no assimilation of observations. The

initialization of the ensemble follows Skachko *et al.* (2014), where at initial time the realizations are spatially correlated homogeneous isotropic covariance using spherical harmonic basis functions. In the present study, the spatial correlation are Gaussian correlations in both horizontal and vertical directions with horizontal correlation length-scales fixed to 400, 800 and 1,600 km and correlated over one level vertically. The standard deviation has a value of 20% from the initial state. The ensemble size used in the experiments was 20, 100 and 500 members.

Figure 4 depicts the spatial correlation after one day of model integration at 50 hPa. All three panels display
<figure><img src="image_15.png" /><figcaption>FIGURE 5 Total error variance evolution using ERA-Interim winds for different initial correlation length using 20 and 100 members</figcaption></figure>

the correlation between one point in the Southern Hemisphere surf zone and the rest of the domain. Each panel presents the result using a different ensemble size -20, 100 and 500 members. Note that no localization is applied here. The reference point was chosen at a time where a Rossby wave breaking event took place, thus displaying strong anisotropic correlation structures, due to the strong wind shear. It is remarkable that with only 500 members it appears that spurious correlations are less than 0.2 in absolute value.

The evolution of the total error variance (properly weighted by the area of each grid box) and normalized with the initial total error variance as a function of ensemble size and initial correlation length is shown in Figure 5. We note that the total variance loss is not sensitive to ensemble size, but depends primarily on the initial correlation length. For 400 km, it is being reduced by 50% in 5 days whereas it takes 26 days (nearly a month) to obtain a similar reduction of the total error variance for an initial correlation length of 1,600 km.

In another experiment we prescribed the (zonal) winds to simulate a solid-body rotation about the polar axis. The total variance loss presented in Figure 6 is for an ensemble of 20 members initialized at 800 km, for the same time period. The result using solid-body winds is depicted with the red curve, and we have added the result of the previous experiment using ERA-Interim reanalysis as the blue curve (corresponding to the green curve in Figure 5). The presence of wind shear, which is absent from the

solid-body rotation winds, does contribute significantly to the variance loss.

These experiments shows that variance loss depends primarily on the correlation length and occurs even without wind shear, although wind shear greatly amplify the loss rate. We also note that the loss of variance with the 3D transport model after 2 days with a correlation length of ≈ 4Δx is somewhat comparable to (perhaps a little larger than) the loss of variance for 1D semi-Lagrangian with cubic spline interpolation (with a uniform grid).

There is thus an indication that, between analyses of dense observations with spatially uncorrelated errors, an important loss of variance is likely to occur during the propagation of error covariance. More will be presented in next section.

## 6 | INFLATION FOR NUMERICAL DISCRETIZATION ERROR

So far, we have examined how variance loss is produced during a forecast. We will now turn to an important consideration on how to fix the variance loss using an inflation. The reader should note that inflation has been used for several purposes as discussed in the Introduction. However, here we will consider inflation related only to numerical discretization.
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 6 Total error variance evolution using solid-body winds and using 20 members (red), and using initial error correlation length of 800 km (blue, same as green circles in Figure 5)</figcaption></figure>

To conduct this study we will use (as an illustration) the first-order semi-Lagrangian scheme with linear interpolation. A full exploration of inflation (and variance loss) for other numerical schemes could be a significant development which would require a study on its own, which we will not perform at this point.

To verify how the inflation works, we have conducted a KF (and EnKF) experiment with a first-order semi-Lagrangian transport scheme with linear interpolation:

(a) In the context of a convergent/divergent wind flow using a deterministic inflation formulation,

(b) As (a), but using an EnKF approach with a stochastic inflation formulation (derived in Appendix C), and

(c) In the context of observing a large region using observations with spatially uncorrelated errors.

We should remark that inflation is often found to be significant in relation to the assimilation of observations, and thus case (c) could be important. Inflation for the parametric KF formulation is also be discussed in Appendix D.

In this study, the inflation is computed from the variance loss $I_n^D(i)$ as it is needed to compensate for this loss, and thus is simply given by

$$I_n^D(i) = -L_n^D(i), \quad (35)$$

where $I_n^D(i)$ is the inflation of variance at grid point $i$ at timestep $n$.

## 6.1 | Inflation in convergent/divergent wind flow

We conducted an experiment as in Section 3, with the same parameters. We considered a 1D periodic domain of $N = 50$ wavenumber (i.e., 101 grid points) with the domain length of $1.0 \times 10^7$ m and thus a grid size of 99.01 km. We chose a timestep of a half-hour; we integrate for one day using 48 timesteps. The initial correlation is a SOAR model with a correlation length $L_c = 400$ km $\approx 4\Delta x$ (i.e., $p_{SOAR} = 296.4$ km). The wind was chosen as a sinusoid with a wind amplitude of $16.0433 \text{ m} \cdot \text{s}^{-1}$ thus giving a maximum Courant number of 0.2917.

The inflation was applied as follows. The inflation/variance loss $I_n^D(i)$ (Equations (24),(35)) for the semi-Lagrangian scheme with linear interpolation was added (at each model timestep), on the computed error variance (or diagonal of the covariance matrix – not the entire covariance matrix),

$$\mathbf{P}(t_{n+1}) = \mathbf{M} \mathbf{P}(t_n) \mathbf{M}^T + \text{diag}(I_n^D). \quad (36)$$

First we first note that with this (inflated) error variance the KF error variance remains equal to 1 at each each grid point. The error covariance after 48 timesteps (one day) is displayed in Figure 7.

It is remarkable that, although only the error variance has been inflated at each timestep, the spatial correlation is not far from to the truth covariance, that
<figure><img src="image_18.png" /><figcaption>FIGURE 7 Kalman filter error covariance with inflation, after 24 timesteps</figcaption></figure>

is the covariance propagated along the characteristics; Figure 1c. A closer examination reveals that in fact it lies somewhere between the covariance propagated in local covariance space and the covariance propagated along the characteristics. Although the inflation was only applied along the diagonal of the matrix, it appears that the effect of model diffusion at each timestep smooths out the spatially uncorrelated structure created by the variance inflation (along the diagonal of the covariance matrix), with the net result that it changes the spatial correlation.

It may also be interesting to see how this inflation varies in time and space. We displayed $I_n^D$ in Figure 8 each 6 hr from 0 to 24 hr. We note that the inflation is different at each grid point and varies in time differently. The inflation values depend on the covariance values themselves.

## 6.2 | Inflation for an ensemble Kalman filter

An inflation for an EnKF should be expressed in terms of the state vector or statistics of it. In Appendix C, we show how this is done. First, we note that the analytical expression (24) can be rewritten using an ensemble of model state as

$$I_n^D(i) = \begin{cases} A_i \langle (c_i^n - c_{i-1}^n)^2 \rangle & \text{for } u_i^n \ge 0, \\ A_i \langle (c_i^n - c_{i+1}^n)^2 \rangle & \text{for } u_i^n < 0, \end{cases} \quad (37)$$

where $\langle \rangle$ is the ensemble mean, and

$$A_i^n = |\alpha_i^n|(1 - |\alpha_i^n|). \quad (38)$$

<figure><img src="image_19.png" /><figcaption>FIGURE 8 Inflation $I_n^D$ at $n = 0, 6, 12, 18, 24$ hr</figcaption></figure>

To have an algorithm that does not use the statistics of the ensemble but rather uses the model state itself can be obtained by taking the square root of the above expression (without the ensemble mean) and multiplying the corrections with a spatially uncorrelated random variable (and independent of the state). The statistics of the stochastic inflation then becomes spatially uncorrelated with a variance equal to (37). The derivation is made in Appendix C and takes the form

$$c_i^{n+1} = c_i^n + \alpha_i^n (c_{i-1}^n - c_i^n) + \sqrt{A(i)} (c_i^n - c_{i-1}^n) U_i \quad (39)$$

for positive winds, and where $U$ is a spatially uncorrelated Gaussian random variable of unit variance,

$$U \sim N(0, I). \quad (40)$$

Using the same parameter values as in Section 3 (i.e., for Figure 7) we have conducted an ensemble KF with a very large ensemble size $N_e = 10,000$. The results after 48 timesteps (i.e., 24 hr) are illustrated in Figure 9. We observe that the ensemble error variance remains close to 1 at all times, and the correlation structure is close to the KF covariance with inflation (Figure 7).

This inflation for an EnKF is in fact a doubly stochastic process, which is introduced by multiplying a spatially uncorrelated random process with the state variable.

## 6.3 | Inflation after analysis over a densely observed region

It is observed that the need for inflation is most significant when observations are assimilated (Zhang et al., 2004; Whitaker and Hamill, 2012; Houtekamer and Zhang,
<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>FIGURE 9 Ensemble Kalman filter error covariance with inflation, after 48 timesteps</figcaption></figure>

2016). It is usually argued that errors in the analysis step itself, such as undersampling and the lack of cross-validation, contribute to variance loss and the need for inflation. Here, we will see that inflation is also needed in the propagation part, that is, immediately after the analysis step and up to the next analysis step.

As in previous subsections, we conduct the forecast of the error covariance with the same model and grid parameters, but this time using a non-homogeneous initial error covariance, obtained after an analysis where half of the domain (from $-L/4$ to $L/4$) is observed at each grid point. The observation-error variance is set to be $1/4$ of the background-error variance and the observation error is assumed to be spatially uncorrelated. The wind is also uniform and constant. We will examine the forecast-error covariance immediately after an analysis and let it be propagated for 12 hr for the next analysis. During that time period, the error covariance should (in theory) simply be a horizontal translation of the initial error covariance (obtained immediately after the first analysis). After 12 hr (or 24 half-hour timesteps) and given the flow parameters as in the previous subsections, this should result in a total displacement of exactly seven grid points.

The result of propagating the error variance with a semi-Lagrangian scheme with linear interpolation, for a period of 12 hr immediately after an analysis and up to the next analysis, is illustrated in Figure 10. Figure 10a shows the results with no inflation. The initial error variance is depicted with a dashed blue line. In principle, the initial error variance should be translated over seven grid points and 12 hr as a result of a uniform constant wind, depicted with a solid blue line, but this is not what the propagation of error covariance with a KF gives. The error

<figure><img src="image_22.png" /><figcaption>FIGURE 10 Error variance 12 hr after an analysis over a densely observed region. The winds are uniform and constant. The dashed blue line shows the error variance right after an analysis. The solid blue line is the translation of the initial error variance equivalent to a transport for 12 hr. The solid red line is the error variance produced by a Kalman filter. (a) shows results with no inflation and (b) with inflation as given by Equations (24) and (35). The green dashed curve is the advection of variance using the discretized model. The green dashed curve is superimpose on the red curve</figcaption></figure>

variance of the KF is displayed by a red line. There is a reduction of error variance from 1.0 to about 0.8 in the unobserved region, and a stronger reduction of variance from 0.07 to 0.03 in the observed region (about a half). Since the effect of using observations with spatially uncorrelated errors is to reduce the correlation length, then according to Equation (25) the reduction of error variance is larger in the observed region than in the unobserved region. We emphasize that this non-uniform reduction of error variance is done *between* the analyses, during the propagation step. We should add that there is no sampling error or cross-validation error in these experiments, since they were conducted using a KF and not with an ensemble DA.
<figure><img src="image_23.png" /><figcaption>FIGURE 11 As Figure 10b, but after 3 days. Note that the winds have pushed the error variance valley to nearly half of the periodic domain after 3 days</figcaption></figure>

We then added the deterministic inflation (Equations (24) and (35)) at each timestep during that time period. The result is illustrated in Figure 10b. The KF error variance plus inflation (red curve) is nearly identical to the translated solution. Both maximum and minimum error variance are the same, but there appears to be some diffuse effects near the sharp gradient of error variance. This is as expected, and we should emphasize that adding back the lost variance as derived in Section 2.2 is to get the advection of variance using the discrete model, in occurrence the semi-Lagrangian model with linear interpolation. In Figure 10b we superposed the solution of advection of initial variance using the discrete model with a dashed green line. We then observe that the error variance obtained with the KF plus inflation is identical to the discrete advection of variance (green dash curve). An identical experiment over a period of 3 days, which corresponds to the revisit time of some atmospheric composition satellites (e.g., Greenhouse Gases Observing Satellite, GOSAT), has also been conducted (Figure 11). The match of error variance of the KF plus inflation and advection of variance is still maintained, but the accumulated effect of diffusion on error variance gradient is more important.

We also briefly examined an inflation scheme which would recover the continuous (space/time) advection (conservation) of error variance (instead of recovering the solution of the discretized advection of error variance). The development is based on the work of Pannekoucke et al. (2021) on obtaining the continuous (space/time) evolution of error variance of a discrete numerical scheme. The inflation method and some preliminary results are presented in Appendix D. It appears that some more work

may be required to obtain a non-dispersive solution, so further development has been postponed to another study.

The results presented in this section may help to understand the results of other EnKF experiments. For example Miyoshi (2011) (figure 11 therein) using a local adaptive inflation scheme found that an inflation factor much greater than 1 was needed over a region densely observed by radiosondes, while over the ocean the inflation factor was found to be around 1 and a little lower in certain areas.

## 7 | CONCLUSIONS

The main conclusion from this study is that the standard formulation for computing the predictability error covariance within an ensemble approach ($XX^T$) or deterministic KF ($MPM^T$) introduces a significant variance loss related to the model discretization error. The effects of numerical discretization on the state and on the error covariance are different because the error covariance field and the state field have different structures. This can be understood from the computation carried out with the KF form. In a KF, the model matrix operates on columns of the covariance matrix $P$ and then on the rows. Yet, the structure of each of these columns (and rows) is for the most part flat, close to zero (if not equal to zero) and then abruptly changes near the diagonal of the $P$ matrix (i.e., the error variance). This type of structure in the error covariance (columns or rows) may be strongly diffused even with a slightly diffusive numerical scheme. The maximum value (i.e., error variance) is then lowered, thus giving rise to a variance loss.

We are able to draw this conclusion and interpretation using the advection equation as a model where the evolution of the error covariance function is known (Cohn, 1993) and two other solution methods could be used. The first method derives from a solution by the method of characteristics which Cohn (1993) eloquently applied to the covariance evolution equation and showed that error variance should also obey the advection dynamics. The second method derives from the change of coordinates of the covariance function and corresponding evolution equation written in terms of $\chi = (x_1 + x_2)/2$ and $\xi = (x_1 - x_2)/2$ (Cohn, 1993). In this new coordinate system, called local covariance coordinates, the evolution of error covariance function (which takes the form of an advection equation) is made parallel to the main diagonal and is then followed by an evolution parallel to the main anti-diagonal, in the form of an advection equation. In this second step, the advecting winds are actually equal to zero at the diagonal. As a result, both solution methods are nearly diffusion-free and conserve the error variance in the sense that they obey the advection
dynamics. Furthermore, the same numerical discretiza-
tion can be applied to local covariance propagation and
the standard form for the predictability. Yet, the numer-
ical solutions can be very different. Thus, it is not just
only a matter of numerical discretization error, but also
how the covariance evolution equation is solved, that
dictates the numerical discretization error of the error
covariance.

Using the first-order semi-Lagrangian scheme with
linear interpolation in 1D as a numerical model, we
derived an analytical expression for the variance loss in
the standard predictability computation. We show that
the variance loss depends on diagonal and the nearby
off-diagonal elements of the error covariance matrix. For
homogeneous isotropic error covariances, the analytical
expression depends simply on the ratio of the grid spac-
ing to the correlation length. Smaller correlation length
gives rise to larger variance loss. To illustrate the gener-
ality of this issue, we conducted an experiment using an
ensemble of a state-of-the-art 3D transport model using a
flux-form semi-Lagrangian scheme (i.e., tpcore of Lin and
Rood 1996) with a piecewise parabolic method (PPM) hor-
izontal discretization and a Prather scheme in the vertical
 Errera et al., 2008; Skachko et al., 2014). ERA-Interim
winds were used to conduct the experiments. We simu-
lated nearly a month of error covariance propagation in
September and where large wind shear is produced in
the Southern Hemisphere surfzone giving rise to highly
anisotropic associated correlation structures known to
give rise to large loss of error variance (Ménard et al.,
2000). Depending on the initial correlation length, the
variance loss in 5 days can vary from 50%, when the ini-
tial correlation length is 400 km (just slightly larger than
the grid size) to 10% for an initial correlation length of
1,600 km. Changing the ensemble size (from 20 to 100)
did not have much influence on the results. We also noted
that variance loss also occurs in the absence of shear
flow, for example, when using a solid-body rotation wind
field. What is striking from these realistic experiments
is that the loss of variance is much larger than what the
1D model simulation seems to suggest. We should note
that the parameters used in nearly all 1D experiments
conducted in this study were chosen to mimic the (zonal
direction) condition with the 3D transport model with an
initial correlation length-scale of 800 km.

We also derived, for first-order numerical schemes in
1D, a variance inflation scheme for the (deterministic)
KF, and a stochastic inflation scheme on the state that
exactly compensate for the variance loss. We found that
inflating only the error variance still provide a reason-
ably accurate error covariance field, as the model diffusion
spreads the variance inflation to spatial correlations. We
show that in a 1D experiment with convergent/divergent

winds, the variance loss can be recovered at each timestep
and each grid point, using our inflation formulation. We
also observed an improvement of the error correlation. We
also conducted experiments with a constant and uniform
wind, where half of the domain had observations at each
grid point. Because we assumed spatially uncorrelated
observation errors, the analysis-error correlation is signifi-
cantly reduced in the observed region. We then propagated
the error covariance for 12 hr (as in a 12-hr assimilation
cycle). After 12 hr, the propagated error covariance could
be though as the background-error covariance to process
the next set of observations. We showed that, after 12 hrs
over the observed region, the computed error variance of
the background-error covariance is about half of its true
value. By applying our inflation scheme, we recover the
variance loss almost entirely. In fact, adding the infla-
tion results in a variance field which obeys the numerical
solution advection of the initial non-uniform error vari-
ance, leading to some diffusion near the boundaries of
the observed/non-observed region. To recuperate the true
error variance by inflation requires the use the continuous
equation of variance loss developed by Pannekoucke et al.
(2021). We report preliminary experiments with this case
continuous inflation scheme in Appendix D.

We wish to point out that in meteorological DA with
ensemble KF, the need for inflation is particularly impor-
tant in the context of the analysis where sampling errors
and “inbreeding” effects contribute to underestimating
the error variance. For example in the work of Miyoshi
(2011) on adaptative inflation, it was found that an infla-
tion coefficient could reach a value of 1.25 over a region
densely observed by radiosondes, but elsewhere would be
nearly neutral, that is, 1, where there are no radiosondes.
Considering the fact that, over a dense radiosonde net-
work, the (analysis-)error correlation length is very small,
we may argue that the propagation step (i.e., between
analyses) could also contribute to an underestimation of
error variance.

Overall, it thus appears that the standard formula-
tion for computing the predictability error covariance,
although general, is one that can introduce a significant
variance loss. We should emphasize that variance loss
due to numerical discretization error arises whether a KF
approach or an ensemble KF is considered. The inflation
needed to correct this effect can be important, even with
state-of-the-art models and especially immediately after an
analysis of dense observations with spatially uncorrelated
errors. This study thus suggests that in meteorological DA
using ensemble methods, the loss of variance that is nor-
mally attributed to the analysis step (i.e., sampling error
and inbreeding errors) may not be the only source of
underestimation, and that variance loss due to numerical
discretization could also play also a role.
## AUTHOR CONTRIBUTIONS

**Richard Ménard:** *Sergey Skachko:* Software; validation; visualization; writing Review Editing. **Olivier Pannekoucke:** investigation; validation; writing Review Editing.

## ACKNOWLEDGEMENTS

We would like to thank Michael Sitwell and Mateusz Reszka for proofreading the document. Also Quentin Errera and Simon Chabrillat of the Belgian Institute for Space Aeronomy for providing us with the 3D transport model. OP is supported by a KAP grant from LEFE INSU. The authors thank two anonymous reviewers for their constructive criticism.

## ORCID

**Richard Ménard** https://orcid.org/0000-0001-7073-5393

**Olivier Pannekoucke** https://orcid.org/0000-0002-3249-2818

## REFERENCES

Bates, J.R. and McDonald, A. (1982) Multiply-upstream, semi-Lagrangian advective schemes: analysis and application to a multi-level primitive equation model. *Monthly Weather Review*, 110, 1831–1842. http://dx.doi.org/10.1175/1520-0493(1982)110<1831:muslas>2.0.co;2.

Beljaars, A., Balsamo, G., Bechtold, P., Bozzo, A., Forbes, R., Hogan, R.J., Köhler, M., Morcrette, J.-J., Tompkins, A.M., Viterbo, P. and Wedi, N. (2018) The numerics of physical parametrization in the ECMWF model. *Frontiers in Earth Science*, 6, 137. http://dx.doi.org/10.3389/feart.2018.00137.

Berner, J., Shutts, G.J., Leutbecher, M. and Palmer, T.N. (2009) A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ECMWF ensemble prediction. *Journal of the Atmospheric Sciences*, 66, 603–626. http://dx.doi.org/10.1175/2008jas2677.1.

Bocquet, M. (2011) Ensemble Kalman filtering without the intrinsic need for inflation. *Nonlinear Processes in Geophysics*, 18(5), 735–750. http://dx.doi.org/10.5194/npg-18-735-2011.

Bonavita, M., Hólm, E., Isaksen, L., Fisher, M. (2016) The evolution of the ECMWF hybrid data assimilation system. *Quarterly Journal of the Royal Meteorological Society*, 142(694), 287–303. http://dx.doi.org/10.1002/qj.2652.

Buizza, R., Milleer, M., Palmer, T. N. (2007) Stochastic representation of model uncertainties in the ECMWF ensemble prediction system. *Quarterly Journal of the Royal Meteorological Society*, 125(560), 2887–2908. http://dx.doi.org/10.1002/qj.49712556006.

Cohn, S.E. (1993) Dynamics of Short-Term Univariate Forecast Error Covariances. *Monthly Weather Review*, 121(11), 3123–3149. http://dx.doi.org/10.1175/1520-0493(1993)121<3123:dostuf>2.0.co;2.

Daley, R. (1991) *Atmospheric Data Analysis*. Cambridge, UK: Cambridge University Press.

Daley, R. (1995) Estimating the Wind Field from Chemical Constituent Observations: Experiments with a One-Dimensional Extended Kalman Filter. *Monthly Weather Review*, 123(1), 181–198. http://dx.doi.org/10.1175/1520-0493(1995)123<0181:etwffc>2.0.co;2.

Dee, D.P., Uppala, S.M., Simmons, A.J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M.A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A.C.M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A.J., Haimberger, L., Healy, S.B., Hersbach, H., Hólm, E.V., Isaksen, L., Källberg, P., Köhler, M., Matricardi, M., McNally, A.P., Monge-Sanz, B.M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., Vitart, F. (2011) The ERA-Interim reanalysis: configuration and performance of the data assimilation system. *Quarterly Journal of the Royal Meteorological Society*, 137(656), 553–597. http://dx.doi.org/10.1002/qj.828.

Errera, Q., Daerden, F., Chabrillat, S., Lambert, J.C., Lahoz, W.A., Viscardy, S., Bonjean, S., Fonteyn, D. (2008) 4D-Var assimilation of MIPAS chemical observations: ozone and nitrogen dioxide analyses. *Atmospheric Chemistry and Physics*, 8(20), 6169–6187. http://dx.doi.org/10.5194/acp-8-6169-2008.

Evensen, G. (2003) The Ensemble Kalman Filter: theoretical formulation and practical implementation. *Ocean Dynamics*, 53(4), 343–367. http://dx.doi.org/10.1007/s10236-003-0036-9.

Grudzien, C., Bocquet, M. and Carrassi, A. (2020) On the numerical integration of the Lorenz-96 model, with scalar additive noise, for benchmark twin experiment. *Geoscientific Model Development*, 13, 1903–1924. http://dx.doi.org/10.5194/gmd-13-1903-2020.

Hamill, T.M. and Whitaker, J.S. (2005) Accounting for the error due to unresolved scales in ensemble data assimilation: a comparison of different approaches. *Monthly Weather Review*, 133, 3132–3147. http://dx.doi.org/10.1175/mwr3020.1.

Hatfield, S., Duben, P., Chantry, M., Kondo, K. and Miyoshi, T. (2018) Choosing the optimal numerical precision for data assimilation in the presence of model error. *Journal of Advances in Modeling Earth Systems*, 10, 2177–2191. http://dx.doi.org/10.1029/2018ms001341.

Houtekamer, P.L., Deng, X., Mitchell, H.L., Gagnon, N. and Baek, S.-J. (2014) Higher resolution in an operational ensemble Kalman filter. *Monthly Weather Review*, 142, 1143–1162. http://dx.doi.org/10.1175/mwr-d-13-00138.1.

Houtekamer, P.L., Mitchell Herschel, L. (1998) Data Assimilation Using an Ensemble Kalman Filter Technique. *Monthly Weather Review*, 126(3), 796–811. http://dx.doi.org/10.1175/1520-0493(1998)126<0796:dauaek>2.0.co;2.

Houtekamer, P.L., Zhang, F. (2016) Review of the Ensemble Kalman Filter for Atmospheric Data Assimilation. *Monthly Weather Review*, 144(12), 4489–4532. http://dx.doi.org/10.1175/mwr-d-15-0440.1.

Khattatov, B.V., Lamarque, J., Lyjak, L., Ménard, R., Levelt, P., Tie, X., Brasseur, G.P. and Gille, J.C. (2000) Assimilation of satellite observations of long-lived chemical species in global chemistry transport models. *Journal of Geophysical Research*, 105(D23), 29135–29144.

Lin, S.-J., Rood, R.B. (1996) Multidimensional Flux-Form Semi-Lagrangian-Transport Schemes. *Monthly Weather Review*, 124(9), 2046–2070. http://dx.doi.org/10.1175/1520-0493(1996)124<2046:mffslt>2.0.co;2.

Lyster, P., Cohn, S., Zhang, B., Chang, L.-P., Ménard, R., Olson, K. and Renka, R. (2004) A Lagrangian trajectory filter for constituent data assimilation. *Quarterly Journal of the Royal Meteorological Society*, 120(1), 181–198. http://dx.doi.org/10.1175/1520-0493(2004)120(1):181-198;2004.0181
<figure><img src="image_24.png" /></figure>

*Meteorological Society*, 130, 2315-2334. http://dx.doi.org/10.1256/qj.02.234.

Ménard, R., Chang, L.-P. (2000) Assimilation of Stratospheric Chemical Tracer Observations Using a Kalman Filter. Part II:χ²-Validated Results and Analysis of Variance and Correlation Dynamics. *Monthly Weather Review*, 128(8), 2672-2686.
http://dx.doi.org/10.1175/1520-0493(2000)128<2672aoscto>2.0.co;2.

Ménard, R., Cohn, S.E., Chang, L.-P. and Lyster, P.M. (2000) Assimilation of stratospheric tracer observations using a Kalman filter. Part I: formulation. *Monthly Weather Review*, 128, 2654-2671.
http://dx.doi.org/10.1175/1520-0493(2000)128<2654aoscto>2.0.co;2.

Ménard, R., Deshaies-Jacques, M., Gasset, N. (2016) A comparison of correlation-length estimation methods for the objective analysis of surface pollutants at Environment and Climate Change Canada. *Journal of the Air & Waste Management Association*, 66(9), 874-895. http://dx.doi.org/10.1080/10962247.2016.1177620.

Meng, Z., Zhang, F. (2007) Tests of an Ensemble Kalman Filter for Mesoscale and Regional-Scale Data Assimilation. Part II: Imperfect Model Experiments. *Monthly Weather Review*, 135(4), 1403-1423. http://dx.doi.org/10.1175/mwr3352.1.

Mitchell, H.L., Houtekamer, P.L. (2009) Ensemble Kalman Filter Configurations and Their Performance with the Logistic Map. *Monthly Weather Review*, 137(12), 4325-4343. http://dx.doi.org/10.1175/2009mwr2823.1.

Miyoshi, T. (2011) The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter. *Monthly Weather Review*, 139(5), 1519-1535. http://dx.doi.org/10.1175/2010mwr3570.1.

Pannekoucke, O., Berre, L., Desroziers, G. (2008) Background-error correlation length-scale estimates and their sampling statistics. *Quarterly Journal of the Royal Meteorological Society*, 134(631), 497-508. http://dx.doi.org/10.1002/qj.212.

Pannekoucke, O., Ricci, S., Barthelemy, S., Ménard, R., Thual, O. (2016) Parametric Kalman filter for chemical transport models. *Tellus A: Dynamic Meteorology and Oceanography*, 68(1), 31547. http://dx.doi.org/10.3402/tellusa.v68.31547.

Pannekoucke, O., Bocquet, M. and Ménard, R. (2018) Parametric covariance dynamics for the nonlinear diffusive burgers equation. *Nonlinear Processes in Geophysics*, 25, 481-495. https://doi.org/10.5194/npg-25-481-2018

Pannekoucke, O., Ménard, R., El Aabaribaoune, M., Plu, M. (2021) A methodology to obtain model-error covariances due to the discretization scheme from the parametric Kalman filter perspective. *Nonlinear Processes in Geophysics*, 28(1), 1-22. http://dx.doi.org/10.5194/npg-28-1-2021.

Rijshøgaard, L., Cohn, S., Li, Y. and Ménard, R. (1998) The use of spline interpolation in semi-lagrangian transport models. *Monthly Weather Review*, 126, 2008-2016.

Sacher, W., Bartello, P. (2008) Sampling Errors in Ensemble Kalman Filtering. Part I: Theory. *Monthly Weather Review*, 136(8), 3035-3049. http://dx.doi.org/10.1175/2007mwr2323.1.

Skachko, S., Errera, Q., Ménard, R., Christophe, Y., Chabriillat, S. (2014) Comparison of the ensemble Kalman filter and 4D-Var assimilation methods using a stratospheric tracer transport model. *Geoscientific Model Development*, 7(4), 1451-1465. http://dx.doi.org/10.5194/gmd-7-1451-2014.

Skachko, S., Ménard, R., Errera, Q., Christophe, Y., Chabriillat, S. (2016) EnKF and 4D-Var data assimilation with chemical transport model BASCOE (version 05.06). *Geoscientific Model Development*, 9(8), 2893-2908. http://dx.doi.org/10.5194/gmd-9-2893-2016.

Strang, G. (1968) On the construction and comparison of difference schemes. *SIAM Journal of Numerical Analysis*, 5, 506-517.
https://doi.org/10.1137/070541

Whitaker, J.S., Hamill, T.M. (2002) Ensemble Data Assimilation without Perturbed Observations. *Monthly Weather Review*, 130(7), 1913-1924. http://dx.doi.org/10.1175/1520-0493(2002)130<1913:edawpo>2.0.co;2.

Whitaker, J.S., Hamill, T.M. (2012) Evaluating Methods to Account for System Errors in Ensemble Data Assimilation. *Monthly Weather Review*, 140(9), 3078-3089. http://dx.doi.org/10.1175/mwr-d-11-00276.1.

Zhang, F., Snyder, C., Sun, J. (2004) Impacts of Initial Estimate and Observation Availability on Convective-Scale Data Assimilation with an Ensemble Kalman Filter. *Monthly Weather Review*, 132(5), 1238-1253. http://dx.doi.org/10.1175/1520-0493(2004)132<1238:ioieao>2.0.co;2.

**How to cite this article:** Ménard, R., Skachko, S. & Pannekoucke, O. (2021) Numerical discretization causing error variance loss and the need for inflation. *Quarterly Journal of the Royal Meteorological Society*, 147(740), 3498-3520. Available from: https://doi.org/10.1002/qj.4139

# APPENDICES

## A. CONTINUOUS AND DISCRETE SOLUTIONS BY OPERATOR SPLITTING

To prove that the system of Equations (6) and (7), solved in sequence, is a solution of Equation (4), and considering that $R(x_1, x_2, \tau; 0)$ is the solution operator of Equation (4) from time $t=0$ to $t=\tau$,

$$P(x_1, x_2, \tau) = R(x_1, x_2, \tau; 0) P(x_1, x_2, 0), \quad (A1)$$

we will show that it can be expressed (exactly) as a composition of the form,

$$R(x_1, x_2, \tau; 0) = R_2(x_2, \tau; 0) R_1(x_1, \tau; 0), \quad (A2)$$

where $R_1$ is the solution operator of Equation (6) and $R_2$ is the solution operator of Equation (7).

Indeed with operator splitting, the solution of Equation (6) expressed as

$$P^*(x_1, x_2, \tau; 0) = R_1(x_1, \tau; 0) P(x_1, x_2, 0) \quad (A3)$$
is used as the initial guess for solving Equation (7), yielding

$$P(x_1, x_2, \tau; 0) = R_2(x_1, \tau; 0) P^*(x_1, x_2, \tau; 0). \quad (A4)$$

Combining Equations (A3) and (A4) is equivalent to Equation (A2), which demonstrate that operator splitting approach Equations (6,7) is the solution of Equation (4).

It may be instructive to consider an example that will make this property obvious. Say that the wind $u$ does not depend on time. The solution of Equation (4) can be written as

$$P(x_1, x_2, t) = \exp \left\{ -t \left( u(x_1) \frac{\partial}{\partial x_1} + u(x_2) \frac{\partial}{\partial x_2} \right) \right\} P(x_1, x_2, 0), \quad (A5)$$

and the solutions of Equations (6) and (7) as

$$P(x_1, x_2, t) = \exp \left( -t u(x_1) \frac{\partial}{\partial x_1} \right) P(x_1, x_2, 0),$$

$$P(x_1, x_2, t) = \exp \left( -t u(x_2) \frac{\partial}{\partial x_2} \right) P(x_1, x_2, 0).$$

Combining these two solutions sequentially as in operator splitting gives the solution of the covariance propagation in $(x_1, x_2)$ (Equation (A5)).

We can relate the continuous formulation above with the discrete formulation as follows. In a discrete space formulation, a covariance function of the extended spatial space $(x_1, x_2)$ (Subsection 2.1) is now an entry of a covariance matrix $\mathbf{P}$ between pairs of grid points,

$$\mathbf{P} = \begin{pmatrix} P_{11} & P_{12} & \dots & P_{1n} \\ P_{21} & P_{22} & \dots & P_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ P_{n1} & P_{n2} & \dots & P_{nn} \end{pmatrix}, \quad (A6)$$

where

$$P_{ij} = \mathbf{P}\{x_1(i\Delta x), x_2(j\Delta x)\}$$

for regularly-spaced grid points $i, j = 1, \dots, n$.

In the covariance matrix (A6) different columns refer to different $x_2$ coordinates, while elements in a row refer to the $x_1$ coordinate. Thus, applying the solution operator $R_1$ with respect to $x_1$ on the covariance function $P(x_1, x_2)$ (A3) is equivalent to applying the model matrix $\mathbf{M}$ on columns of the covariance matrix $\mathbf{P}$. If $\mathbf{P}_{[c_i]}$ represents the $i$th column of $\mathbf{P}$, then continuous solution $P^*(x_1, x_2, t) = R_1 P(x_1, x_2, 0)$ is represented discretely as

$$\mathbf{P}^*_{[c_i]}(t) = \mathbf{M} \mathbf{P}_{[c_i]}(0)$$

for each column $c_i(i = 1, \dots, n)$. In matrix form (considering all the columns), it is then written as

$$\mathbf{P}^*(t) = \mathbf{M} \mathbf{P}(0).$$

Next, the discrete equivalent of $P(x_1, x_2, t) = R_2 P^*(x_1, x_2, t)$ is a model matrix operating on rows of the covariance $\mathbf{P}^*$, where we operate in sequence as in operator splitting. Let $\mathbf{P}_{[r_j]}$ represent the $j$th row of $\mathbf{P}$ so the continuous solution $P(x_1, x_2, t) = R_2 P^*(x_1, x_2, t)$ is represented discretely as

$$\mathbf{P}_{[r_j]}(t) = \mathbf{M} \mathbf{P}^*_{[r_j]}(t). \quad (A7)$$

To write Equation (A7) in matrix form, first note that in a matrix multiplication $\mathbf{MP}$ the matrix $\mathbf{M}$ normally operates on columns of $\mathbf{P}$, but this is not what is desired in Equation (A7). To operate on rows, take the transpose of the matrix and then operate on the columns of the transposed matrix. In matrix form, Equation (A7) can thus be written as

$$\mathbf{P}(t) = \mathbf{M} [\mathbf{P}^*(t)]^T. \quad (A8)$$

Combining Equation (1) with Equation (A8) we get

$$\mathbf{P}(t) = \mathbf{M}[\mathbf{M} \mathbf{P}(0)]^T = \mathbf{M} \mathbf{P}(0) \mathbf{M}^T, \quad (A9)$$

the standard KF form for a discrete space and discrete model.

Although both continuous and discrete forms can be solved by operator splitting, it will become apparent (as we shall see in this article) that they do not yield the same solutions.

## B. DISCRETIZATION IN LOCAL COVARIANCE COORDINATES

Consider a 1D periodic domain made up of $N$ grid points,

$$x(i) = (i - 1) \Delta x, \quad (i = 1, \dots, N),$$

where $\Delta x = L/N$ is the grid spacing and $L$ is the length of the domain. Consider a pair of points in the physical space with coordinates, $x_1 = x_1(i), x_2 = x_2(j)$ with indices $i$ and $j$ respectively. These two points have as local coordinates

$$\chi(i, j) = \frac{1}{2} [x_1(i) + x_2(j)] = \frac{1}{2}(i + j - 2)\Delta x,$$

$$\xi(i, j) = \frac{1}{2} [x_1(i) - x_2(j)] = \frac{1}{2}(i - j)\Delta x.$$
<figure><img src="image_25.png" /><figcaption>TABLE B1 Slicing in $\xi$</figcaption></figure>

<table><thead><tr><th>m = 0</th><th>m = 1</th><th>m = 2</th></tr></thead><tbody><tr><td>P(1,1); n = 0</td><td>P(2,1); n = 1</td><td>P(3,1); n = 0</td></tr><tr><td>P(2,2); n = 2</td><td>P(3,2); n = 3</td><td>P(4,2)=P(1,2); n = 2</td></tr><tr><td>P(3,3); n = 4</td><td>P(4,3)=P(1,3); n = 5</td><td>P(5,3)=P(2,3); n = 4</td></tr></tbody></table>

In this local coordinate system, we can also define a pair of indices $(n, m)$, as

$$ \chi(n) = n(i, j) \Delta x/2, \quad \xi(m) = m(i, j) \Delta x/2, \quad (B1) $$

where

$$ n(i, j) = i + j - 2, \quad m(i, j) = i - j. \quad (B2) $$

The propagation algorithm then consists of integrating Equation (31) a state vector in $\chi$ for each $\xi(m)$ values (i.e., parallel to the diagonal), followed by integrating Equation (32) a state vector in $\xi$ for each $\chi(n)$ values (i.e., perpendicular to the diagonal).

The key aspect of this scheme is the slicing of the covariance matrix, parallel and perpendicular to the diagonal, as well as defining vectors of equal length. Figure B1 shows the slicing of the original covariance matrix into vectors parallel to the diagonal using, for illustration, a reduced size matrix of $N = 3$. The red lines connected by arrows have all the same $\xi$ value or same $m$ index. On the main diagonal we have $m = 0$, then the next diagonal lower down is for $m = 1$ and the last diagonal is for $m = 2$. The element of the matrix contained in a slice $\xi$ (or $m$) is given by $P(k+m, k)$ for $k = 1, \dots, 3$. For the elements ending up outside the matrix domain, they can be mapped to elements of the matrix by using periodicity (dashed green lines). Table B1 shows the three $\xi$ slices that are used as initial condition for Equation (31). The matrix elements that are in the matrix domain are in bold, and the matrix elements extended by using periodicity are non-bold.

Figure B2 shows the slicing of the original covariance matrix into vectors perpendicular to the diagonal (or anti-diagonal). The blue lines connected by arrows have all the same $\chi$ or $n$ index. The rightmost line is for $n = 2$, then the next to the left is for $n = 1$ and the following on the left is for $n = 0$. The element of the matrix contained in a slice $\chi$ (or $n$) is given by $P(k, n+2-k)$ for $k = 3, \dots, 1$. Table B2 shows the three $\chi$ slices that are used as initial condition for Equation (32).

In general, for the $n, m$ indices, we note that an even (odd) index $n$ can only be associated with an even (odd) index $m$. This property can be obtained from the inverse

<figure><img src="image_26.png" /><figcaption>FIGURE B1 Slicing along $\chi$ for different value of $\xi$</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE B2 Slicing along $\xi$ for different value of $\chi$</figcaption></figure>

transformation of Equation (B2) where from $(n, m)$ we get $(i, j)$ with the knowledge that $i$ and $j$ are integers.

Once the slicing of the covariance matrix is done, the transport by the winds is simpler than may first appear, by noting that from Equations (B1) and (B2), we have

$$ \chi(n) + \xi(m) = (i - 1)\Delta x, \quad \chi(n) - \xi(m) = (j - 1)\Delta x. $$
<table><caption>TABLE B2 Slicing in $\chi$</caption><thead><tr><th>n = 0</th><th>n = 1</th><th>n = 2</th></tr></thead><tbody><tr><td>P(1,1); m = 0</td><td>P(1,2); m = 1</td><td>P(1,3); m = 0</td></tr><tr><td>P(2,0)=P(2,3); m = 2</td><td>P(2,1); m = 3</td><td>P(2,2); m = 2</td></tr><tr><td>P(3,-1)=P(3,2); m = 4</td><td>P(3,0)=P(3,3); m = 5</td><td>P(3,1); m = 4</td></tr></tbody></table>

The winds in the local covariance coordinates (Equations (29) and (30)) are then given by

$$ 
\begin{aligned}
u(\chi(n) + \xi(m)) + u(\chi(n) - \xi(m)) &= u(i-1) + u(j-1), \\
u(\chi(n) + \xi(m)) - u(\chi(n) - \xi(m)) &= u(i-1) - u(j-1),
\end{aligned}
$$

where we dropped the time index to simplify the notation.

### C. INFLATION FOR A FIRST-ORDER SCHEME

We will derive an inflation scheme that compensates for the variance loss as expressed in Equation (24). This derivation applies for first-order schemes, such as semi-Lagrangian with linear interpolation.

It is useful to separate the inflation associated with Equation (24) in two parts,

$$ 
\begin{aligned}
I_1 &= -HOT_1 = \alpha_i^n(1 - \alpha_i^n)[V^n(i) - P^n(i-1, i)], \\
I_2 &= -HOT_2 = \alpha_i^n(1 - \alpha_i^n)[V^n(i-1) - P^n(i, i-1)].
\end{aligned}
\quad (C1) $$

The variances and covariances can be approximated by taking an ensemble approach. Defining $\langle \rangle = \sum_{i=1}^{N_e} \langle \rangle/N_e$ as the ensemble mean, we have

$$ 
\begin{aligned}
V^n(i) &\approx \langle \tilde{c}_i^n \tilde{c}_i^n \rangle, \\
P^n(i-1, i) &\approx \langle \tilde{c}_{i-1}^n \tilde{c}_i^n \rangle, \\
P^n(i, i-1) &\approx \langle \tilde{c}_i^n \tilde{c}_{i-1}^n \rangle,
\end{aligned}
$$

where $N_e$ is the number of ensemble members, and $\langle \rangle$ indicates the departure from the ensemble mean. Note that

$$ 
\begin{aligned}
\langle \tilde{c}_i^n \tilde{c}_i^n \rangle - \langle \tilde{c}_{i-1}^n \tilde{c}_i^n \rangle &= \langle \tilde{c}_i^n - \tilde{c}_{i-1}^n \rangle \tilde{c}_i^n, \\
\langle \tilde{c}_{i-1}^n \tilde{c}_{i-1}^n \rangle - \langle \tilde{c}_i^n \tilde{c}_{i-1}^n \rangle &= \langle (\tilde{c}_i^n - \tilde{c}_{i-1}^n)(-\tilde{c}_{i-1}^n) \rangle,
\end{aligned}
$$

and adding those two expressions gives

$$ \langle (\tilde{c}_i^n - \tilde{c}_{i-1}^n)(\tilde{c}_i^n - \tilde{c}_{i-1}^n) \rangle = \langle (\tilde{c}_i^n - \tilde{c}_{i-1}^n)^2 \rangle, \quad (C2) $$

which is a positive quantity. Substituting Equation (C2) in the expression for inflation (C1), gives

$$ I^n(i) = I_1 + I_2 = \alpha_i^n(1 - \alpha_i^n) \langle (\tilde{c}_i^n - \tilde{c}_{i-1}^n)^2 \rangle. \quad (C3) $$

The next step consists of finding a stochastic version of this equation such that its covariance gives the expression (C3) above for the error covariance(variance), but is zero for the off-diagonal elements. To achieve a stochastic inflation term that is spatially uncorrelated and uncorrelated with the state, we first note that, if $X$ is a random vector independent of a random normally distributed vector $Y \sim N(0, I)$, then in computing the covariance of $X \circ Y$. where $\circ$ is the Schur product, we get

$$ 
\begin{aligned}
\langle (X \circ Y)_i (X \circ Y)_j \rangle &= \langle X_i Y_i X_j Y_j \rangle = \langle X_i X_j Y_i Y_j \rangle \\
&= \langle X_i X_j \rangle \langle Y_i Y_j \rangle = \langle X_i X_j \rangle \delta_{ij}.
\end{aligned}
\quad (C4) $$

Thus, a stochastic inflation of the form

$$ \sqrt{\alpha_i^n (1 - \alpha_i^n) (\tilde{c}_i^n - \tilde{c}_{i-1}^n)} Y_i $$

will generate a diagonal covariance whose variance is exactly the expression in Equation (C3). Thus we propose a stochastic inflated state of the form,

$$ \tilde{c}_i^{n+1} = \tilde{c}_i^n + \alpha_i^n (\tilde{c}_{i-1}^n - \tilde{c}_i^n) + \sqrt{\alpha_i^n (1 - \alpha_i^n)} (\tilde{c}_i^n - \tilde{c}_{i-1}^n) Y_i. \quad (C5) $$

### D. INFLATION FOR THE PARAMETRIC COVARIANCE FORMULATION

For the first-order semi-Lagrangian scheme with linear interpolation, it was established by Pannekoucke et al. (2021) (their equation (34a)) that the KF error variance of the discrete model obeys the equation

$$ \partial_t V + U \partial_x V = - \frac{V}{L_c^2} 2\kappa + \kappa \partial_x^2 V - \kappa \frac{(\partial_x V)^2}{2V}. \quad (D1) $$

The first term on the right-hand side corresponds, as we shall see, to the HOT or $L_n^D$ term discussed in Section 2.2. But there are more terms appearing in Equation (D1). We will briefly discuss how these terms could contribute to an inflation scheme.

The starting point is to discretize the variance propagation Equation (D1) using a semi-Lagrangian scheme with linear interpolation similar to Equations (23–24),

$$ \frac{V_{n+1}(i) - V_n(i)}{\Delta t} + u_i^n \frac{V_n(i) - V_n(i-1)}{\Delta x} = \frac{1}{\Delta t} L_n^C(i), \quad (D2) $$

where we used a constant and uniform wind, $U = u$ (Pannekoucke et al., 2021). The r.h.s. of Equation (D1) represents the tendency of the variance loss, which could be
<figure><img src="image_28.png" /></figure>

written as

$$L_n^C(i) = \left(-\frac{V}{L_c^2}(2\kappa) + \kappa \partial_x^2 V - \kappa \frac{(\partial_x V)^2}{2V}\right)_n(i) \\ = L_n^1(i) + L_n^2(i) + L_n^3(i), \qquad (D3)$$

where

$$L_n^1(i) = -\frac{V}{L_c^2}(2\kappa), \quad |_n(i), \qquad (D4)$$

$$L_n^2(i) = \kappa \partial_x^2 V \Big|_n(i), \qquad (D5)$$

$$L_n^3(i) = -\kappa \frac{(\partial_x V)^2}{2V} \Big|_n(i), \qquad (D6)$$

which still need to be discretized. But before we do this, it is actually useful rewrite the above equation, by multiplying Equation (D2) by $\Delta t$ and then rearranging terms so to yield

$$V_{n+1}(i) = (1 - \alpha_i^n) V_n(i) + \alpha_i^n V_n(i-1) + I_n^C(i), \qquad (D7)$$

an equation which has the same form as the discrete formulation Equation (23). The corresponding equation for the inflated variance $V^I$ should then be

$$V_{n+1}^I(i) = (1 - \alpha_i^n) V_n^I(i) + \alpha_i^n V_n^I(i-1) + I_n^C(i), \qquad (D8)$$

where the inflation term is simply the negative of variance loss term, as was done previously,

$$I_n^C(i) = -L_n^C(i). \qquad (D9)$$

Now let us derive discretized expressions for Equations (D4–D6). Note first that $\kappa$ is related to $\alpha(1 - \alpha)$ as follows

$$\kappa = \frac{u}{2}(\Delta x - u \Delta t) = \frac{\Delta x^2}{2 \Delta t} \alpha(1 - \alpha). \qquad (D10)$$

Then, the first term of the variance loss, Equation (D4), is equal to

$$L_n^1(i) = -\alpha_i^n(1 - \alpha_i^n)V_n(i)\frac{\Delta x^2}{L_c^2}. \qquad (D11)$$

In this study we do not have a separate equation for the correlation length. However, we note that the term $\Delta x^2/L_c^2$ actually appears in Equation (25) and it is simply related to the variance loss associated with the HOT (Equation (24)). Thus we can simply rewrite the equation above and obtain

<figure><img src="image_29.png" /><figcaption>FIGURE D1 Continuous error variance evolution with inflation. Results are after 12 hr of integration using all three terms of the inflation in Equations (D12)–(D14)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>FIGURE D2 As Figure D1, but after 3 days instead of 12 hr</figcaption></figure>

the corresponding inflation as

$$I_n^1(i) = -L_n^1(i) = \alpha_i^n(1 - \alpha_i^n) \\ \times \{V_n(i) - P_n(i, i-1) - P_n(i-1, i) + V_n(i-1)\}. \qquad (D12)$$

The other terms, Equations (D5–D6) depend on error variances and yield the following discretized inflation terms:

$$\begin{aligned} I_n^2(i) &= -L_n^2(i) \\ &= -\alpha_i^n(1 - \alpha_i^n) \frac{[V_n(i-1) - 2V_n(i) + V_n(i+1)]}{2}, \end{aligned} \qquad (D13)$$
<figure><img src="image_31.png" /></figure>

$$ 
\begin{aligned}
I_n^3(i) &= -L_n^3(i) \\
&= \alpha_i^n (1 - \alpha_i^n) \frac{[V_n(i+1) - V_n(i-1)]^2}{16V_n(i)}, \quad (D14)
\end{aligned}
$$

where we used a centred difference approximation for the gradient of variance in the third term. Using these inflation terms should approach the continuous solution, or the solution obtained by pure translation. We have

conducted experiments with the same parameters as in Subsection 6.3, and we display the error variance after 12 hr and 3 days in Figures D1 and D2. These results shows that the numerical scheme with the inflation scheme above is dispersive if applied as is. More research is needed to find the appropriate numerical scheme, but it is out of the scope of this article.