A local particle filter and its Gaussian mixture extension implemented with minor modifications to the LETKF

A particle filter (PF) is an ensemble data assimilation method that does not assume Gaussian error distributions. Recent studies proposed local PFs (LPFs), which use localization, as in the ensemble Kalman filter, to apply the PF efficiently for high-dimensional dynamics. Among others, Penny and Miyoshi (2016) developed an LPF in the form of the ensemble transform matrix of the local ensemble transform Kalman filter (LETKF). The LETKF has been widely accepted for various geophysical systems, including numerical weather prediction (NWP) models. Therefore, implementing the LPF consistently with an existing LETKF code is useful.

This study develops a software platform for the LPF and its Gaussian mixture extension (LPFGM) by making slight modifications to the LETKF code with a simplified global climate model known as Simplified Parameterizations, Primitive Equation Dynamics (SPEEDY). A series of idealized twin experiments were accomplished under the ideal-model assumption. With large inflation by the relaxation to prior spread, the LPF showed stable filter performance with dense observations but became unstable with sparse observations. The LPFGM showed a more accurate and stable performance than the LPF with both dense and sparse observations. In addition to the relaxation parameter, regulating the resampling frequency and the amplitude of Gaussian kernels was important for the LPFGM. With a spatially inhomogeneous observing network, the LPFGM was superior to the LETKF in sparsely observed regions, where the background ensemble spread and non-Gaussianity were larger. The SPEEDY-based LETKF, LPF, and LPFGM systems are available as open-source software on GitHub (https://github.com/skotsuki/speedy-lpf, last access: 16 November 2022) and can be adapted to various models relatively easily, as in the case of the LETKF.

1 Introduction

Ensemble-based data assimilation (DA) has been broadly applied in geoscience fields such as weather and ocean prediction. The ensemble Kalman filter (EnKF) has been intensely investigated for the past two decades, such as for the perturbed observation method (Evensen, 1994; Burgers et al., 1998; Houtekamer and Mitchel, 1998; van Leeuwen, 2020) and for ensemble square root filters (e.g., Anderson, 2001; Bishop et al., 2001; Whitaker and Hamill, 2002; Hunt et al., 2007). The advantages of the EnKF are in its flow-dependent error estimates, represented by an ensemble, and in its relative ease of implementation in nonlinear dynamical systems, such as numerical weather prediction (NWP) models. The degrees of freedom of dynamical models (e.g., >O(10⁸) for NWP models) are much larger than the typically affordable ensemble size (<O(10³)). On the other hand, atmospheric and oceanic models show local low dimensionality (Patil et al., 2001; Oczkowski et al., 2005), and practical EnKF implementations use a localization technique that limits the impact of distant observations while also reducing the effective degrees of freedom. Among various kinds of EnKFs, the local ensemble transform Kalman filter (LETKF; Hunt et al., 2007) is widely utilized in operational NWP centers in the same manner, such as the Deutscher Wetterdienst (DWD) and the Japan Meteorological Agency (JMA). Analysis updates of the LETKF are performed by multiplying the ensemble transform matrix to the prior ensemble perturbation matrix, following the ensemble transform Kalman filter (ETKF, Bishop et al., 2001; Wang et al., 2004).

The particle filter (PF) is another ensemble DA method broadly applicable to nonlinear and non-Gaussian problems (see van Leeuwen, 2009, and van Leeuwen et al., 2019, for reviews on geoscience applications). The PF potentially solves some issues in the basic assumptions of the EnKF by permitting nonlinear observation operators and non-Gaussian likelihood functions (Penny and Miyoshi, 2016). For example, assimilating precipitation observations with a standard EnKF is difficult, partly because of their non-Gaussian errors (Lien et al., 2013, 2016; Kotsuki et al., 2017a). The PF can treat such variables with non-Gaussian errors properly. Several PF methods have been explored for low-dimensional problems in early studies (Gordon et al., 1993; van Leeuwen, 2009). However, applying the PF to high-dimensional dynamical systems is generally difficult, because the number of particles or the ensemble size must be increased exponentially with the system size to avoid a weight collapse in which very few particles occupy most of the weights (Snyder et al., 2008, 2015). If a weight collapse occurs, the PF loses the diversity of the particles after resampling. Therefore, we need to assure that the weights are similarly distributed among particles. Previous studies developed the equivalent-weights particle filter (EWPF; van Leeuwen, 2010; Ades and van Leeuwen, 2013, 2015; Zhu et al., 2016) to extend time until a weight collapse by using the proposal density to drive the particles toward the high-probability region of the posterior. However, even PFs with proposal densities are unable to prevent the weights from collapsing.

Alternatively, the local particle filter (LPF) uses localization to avoid a weight collapse by limiting the impact of observations within a local domain. Localization is a well-adopted method for the EnKF to treat sampling errors due to a limited ensemble size. Several LPFs have been proposed to apply the PF efficiently to high-dimensional systems (e.g., Bengtsson et al., 2003; van Leeuwen, 2009; Poterjoy, 2016; Penny and Miyoshi, 2016; Poterjoy and Anderson, 2016; Farchi and Bocquet, 2018; van Leeuwen et al., 2019). Among them, Penny and Miyoshi (2016) developed an LPF by means of the ensemble transform matrix of the LETKF. The LETKF has been commonly used for diverse geophysical systems, and consistent implementation with an available LETKF code would be useful. The DWD implemented the LPF in an operational global model based on the operational LETKF code and described stable performance in the operational setting (Potthast et al., 2019). Walter and Potthast (2022; hereafter WP22) extended the LPF to the Gaussian mixture for further improvements.

The goal of this study is to develop and provide a software platform to accelerate research on LPF and its Gaussian mixture extension (hereafter, LPFGM). Although many studies on LPF have used simple idealized chaotic models, such as the 40-variable Lorenz model (hereafter, “L96”; Lorenz, 1996; Lorenz and Emanuel, 1998), the present study uses a simplified atmospheric general circulation model, defined as the Simplified Parameterizations, Primitive Equation Dynamics (SPEEDY) model (Molteni, 2003). Miyoshi (2005) coupled the LETKF with the SPEEDY model for the first time. This study extends the SPEEDY-LETKF code and implements the LPF and LPFGM. As demonstrated below, the LPF and LPFGM can be implemented easily, with simple modifications to the existing LETKF code.

The remainder of this paper is organized as follows. Section 2 explains the mathematical formulation of the LETKF, LPF, and LPFGM, followed by a description of the specific modifications made to the existing LETKF code. Section 3 presents the experimental settings. Section 4 outlines the results and discussions. Lastly, Sect. 5 provides a summary.

2 Methodology

2.1 Local ensemble transform Kalman filter

Hunt et al. (2007) introduced the LETKF as a computationally efficient EnKF by combining the local ensemble Kalman filter (Ott et al., 2004) and ETKF (Bishop et al., 2001). Let X_t be a matrix composed of m ensemble state vectors. The ensemble mean vector and perturbation matrix of X_t are given by ${\stackrel{\mathrm{‾}}{\mathit{x}}}_t$ (∈ℝⁿ) and ${\mathbf{Z}}_t\equiv \left\{{\mathit{x}}_t^{\left(\mathrm{1}\right)}-{\stackrel{\mathrm{‾}}{\mathit{x}}}_t,\mathrm{\dots },{\mathit{x}}_t^{\left(m\right)}-{\stackrel{\mathrm{‾}}{\mathit{x}}}_t\right\}$ ($\in {\mathbb{R}}^{n\times m})$ respectively, where n is the system size. The subscript t indicates the time, and the superscript (i) denotes the ith ensemble member. The analysis equations of the LETKF are specified by

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathbf{X}}_{t,\mathrm{LETKF}}^{\mathrm{a}}={\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}+{\mathbf{Z}}_t^{\mathrm{b}}{\mathbf{T}}_{t,\mathrm{LETKF}},\text{(2)}& {\displaystyle }\begin{array}{rl}{\mathbf{T}}_{t,\mathrm{LETKF}}& ={\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}{\left({\mathbf{Y}}_t^{\mathrm{b}}\right)}^T{\mathbf{R}}_t^{-\mathrm{1}}\left({\mathit{y}}_t^{\mathrm{o}}-\stackrel{\mathrm{‾}}{H_t\left({\mathbf{X}}_t^{\mathrm{b}}\right)}\right)\cdot \mathbf{1}\\ & +{\left[\left(m-\mathrm{1}\right){\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}\right]}^{\mathrm{1}/\mathrm{2}},\end{array}\end{array}$$

where 1 is a row vector whose elements are all 1 (∈ℝ^m); T is the ensemble transform matrix ($\in {\mathbb{R}}^{m\times m}$); $\stackrel{\mathrm{̃}}{\mathbf{P}}$ is the error covariance matrix in the ensemble space ($\in {\mathbb{R}}^{m\times m}$); Y≡HZ is the ensemble perturbation matrix in the observation space ($\in {\mathbb{R}}^{p\times m}$); R is the observation error covariance matrix ($\in {\mathbb{R}}^{p\times p}$); y is the observation vector (∈ℝ^p); H is the linear observation operator matrix ($\in {\mathbb{R}}^{p\times n}$); and H is the observation operator that may be nonlinear. Here, p is the number of observations. The superscripts “o”, “b”, and “a” denote the observation, background (prior), and analysis (posterior), respectively. The matrix T_LETKF denotes the ensemble transform matrix of the LETKF. The first and second terms on the right-hand side of Eq. (2) correspond to the updates of ensemble mean and perturbation, respectively. The ETKF computes the analysis error covariance matrix in the ensemble space ${\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}$ by

$$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}={\left[{\displaystyle \frac{\left(m-\mathrm{1}\right)}{\mathit{\beta }}}\mathbf{I}+{\left({\mathbf{Y}}_t^{\mathrm{b}}\right)}^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{Y}}_t^{\mathrm{b}}\right]}^{-\mathrm{1}},\end{array}$$

where β is a multiplicative inflation factor. Equations (1) and (2) are derived from the Kalman filter equations, given by

$$\begin{array}{}\text{(4)}& {\displaystyle }{\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{a}}={\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}+{\mathbf{K}}_t\left({\mathit{y}}_t^{\mathrm{o}}-\stackrel{\mathrm{‾}}{H_t\left({\mathbf{X}}_t^{\mathrm{b}}\right)}\right),\text{(5)}& {\displaystyle }{\mathbf{K}}_t={\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T{\left({\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T+{\mathbf{R}}_t\right)}^{-\mathrm{1}}={\mathbf{P}}_t^{\mathrm{a}}{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}},\text{(6)}& {\displaystyle }\begin{array}{rl}{\mathbf{P}}_t^{\mathrm{a}}& =\left(\mathbf{I}-{\mathbf{K}}_t{\mathbf{H}}_t\right){\mathbf{P}}_t^{\mathrm{b}}\\ & ⟺{\left({\mathbf{P}}_t^{\mathrm{a}}\right)}^{-\mathrm{1}}={\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t,\end{array}\end{array}$$

where K is the Kalman gain ($\in {\mathbb{R}}^{n\times p}$), and P is the error covariance matrix in the model space ($\in {\mathbb{R}}^{n\times n}$). Derivations of Eqs. (5) and (6) are detailed in Appendix A. The EnKF approximates the error covariance matrix by $\mathbf{P}\approx \frac{\mathrm{1}}{m-\mathrm{1}}\mathbf{Z}{\left(\mathbf{Z}\right)}^T$. Hunt et al. (2007) provide more details on deriving Eqs. (1) and (2) from the Kalman filter equations. For nonlinear observation operators, the following approximation is used:

$$\begin{array}{}\text{(7)}& \mathbf{Y}\equiv \mathbf{HZ}\approx H\left(\mathbf{X}\right)-\stackrel{\mathrm{‾}}{H\left(\mathbf{X}\right)}\cdot \mathbf{1}.\end{array}$$

The LETKF provides ${\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}$ and ${\left({\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}\right)}^{\mathrm{1}/\mathrm{2}}$ of Eq. (2) by solving the eigenvalue decomposition of ${\left({\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}\right)}^{-\mathrm{1}}=\mathbf{U}\mathbf{\Lambda }{\mathbf{U}}^T$, where U is a square matrix ($\in {\mathbb{R}}^{m\times m}$) composed of eigenvectors, and Λ is a diagonal matrix ($\in {\mathbb{R}}^{m\times m}$) composed of the corresponding eigenvalues. The eigenvalue decomposition leads to ${\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}=\mathbf{U}{\mathbf{\Lambda }}^{-\mathrm{1}}{\mathbf{U}}^T$ and ${\left({\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}\right)}^{\mathrm{1}/\mathrm{2}}=\mathbf{U}{\mathbf{\Lambda }}^{-\mathrm{1}/\mathrm{2}}{\mathbf{U}}^T$.

The localization is practically important for mitigating sampling errors in the ensemble-based error covariance with a limited ensemble size (Houtekamer and Zhang, 2016). With localization, the LETKF computes the transform matrix T_LETKF at every model grid point independently by assimilating a subset of observations within the localization cut-off radius. The LETKF employs localization by inflating the observation error variance so that observations distant from the analysis model grid point have fewer impacts (Hunt et al., 2007; Miyoshi and Yamane, 2007).

For local analysis schemes, a spatially smooth transition of the transform matrix T_LETKF is essential to prevent abrupt changes in the analyses of neighboring grid points. The LETKF realizes a smooth transition of the transform matrix by using the symmetric square root of ${\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}$ (Hunt et al., 2007). The symmetric square root matrix minimizes the mean-square distance between identity matrix I and ${\left[\left(m-\mathrm{1}\right){\stackrel{\mathrm{̃}}{\mathbf{P}}}_t^{\mathrm{a}}\right]}^{\mathrm{1}/\mathrm{2}}$; therefore, the analysis ensemble perturbations can be closer to the background ensemble perturbation.

The EnKF generally underestimates the error variance, mainly because of model errors, nonlinear dynamics, and limited ensemble size. Therefore, a covariance inflation technique is used to inflate the underestimated error variance (Houtekamer and Mitchell, 1998). Among several kinds of covariance inflations methods, the present study considers multiplicative inflation (Anderson and Anderson, 1999) and relaxation to the prior scheme (Zhang et al., 2004), as implemented by Whitaker and Hamill (2012). In multiplicative inflation, the ensemble-based covariance is multiplied by a factor β (P^b→βP^b). This multiplicative inflation is employed in Eq. (3) in the LETKF. In this study, the LETKF experiments use the approach of Miyoshi (2011), which adaptively estimates a spatially varying inflation factor on the basis of the innovation statistics of Desroziers et al. (2005). In realistic problems, covariance relaxation methods are often used to inflate the posterior perturbation (e.g., Terasaki et al., 2019; Kotsuki et al., 2019b). This study utilizes the relaxation to prior spread (RTPS; Whitaker and Hamill, 2012), given by

$$\begin{array}{}\text{(8)}& {\mathbf{Z}}_{t\left(k\right)}^{\mathrm{a}}\leftarrow \left(\left(\mathrm{1}-\mathit{\alpha }\right)+\mathit{\alpha }{\displaystyle \frac{{\mathit{\sigma }}_{t\left(k\right)}^{\mathrm{b}}}{{\mathit{\sigma }}_{t\left(k\right)}^{\mathrm{a}}}}\right){\mathbf{Z}}_{t\left(k\right)}^{\mathrm{a}},\end{array}$$

where α is the RTPS parameter, and σ is the ensemble spread. The subscript (k) denotes the kth model variable or the kth component of the state vector. Although this study only uses multiplicative inflation for the LETKF experiments, the posterior error perturbation is inflated by RTPS for LPF and LPFGM.

2.2 Local particle filter with ensemble transform matrix

Here, we describe the LPF in the form of the ensemble transform matrix of the LETKF (Reich, 2013; Penny and Miyoshi, 2016; Potthast et al., 2019). The PF is a direct Monte Carlo realization of Bayes' theorem, given by

$$\begin{array}{}\text{(9)}& \mathit{\pi }\left({\mathit{x}}_t\mathrm{|}{\mathit{y}}_{\mathrm{1}:t}^{\mathrm{o}}\right)={\displaystyle \frac{\mathit{\pi }\left({\mathit{y}}_t^{\mathrm{o}}\mathrm{|}{\mathit{x}}_t\right)\mathit{\pi }\left({\mathit{x}}_t\mathrm{|}{\mathit{y}}_{\mathrm{1}:t-\mathrm{1}}^{\mathrm{o}}\right)}{\mathit{\pi }\left({\mathit{y}}_t^{\mathrm{o}}\mathrm{|}{\mathit{y}}_{\mathrm{1}:t-\mathrm{1}}^{\mathrm{o}}\right)}},\end{array}$$

where $\mathit{\pi }\left({\mathit{x}}_t\mathrm{|}{\mathit{y}}_{\mathrm{1}:t}^{\mathrm{o}}\right)$ is the probability of state x given all observations y^o up to time t. The PF approximates the prior probability density function (PDF), which appears in the numerator of Eq. (9), using an ensemble forecast

$$\begin{array}{}\text{(10)}& {\mathit{\pi }}_{\mathrm{LPF}}\left({\mathit{x}}_t\mathrm{|}{\mathit{y}}_{\mathrm{1}:t-\mathrm{1}}^{\mathrm{o}}\right)\approx {\sum }_{i=\mathrm{1}}^m{\mathit{w}}_t^{\mathrm{b}\left(i\right)}\mathit{\delta }\left({\mathit{x}}_t-{\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right),\end{array}$$

where δ is the Dirac's delta function, and w^b(i) is the prior weight of the ith particle (i.e., ensemble member) at the previous analysis time. If particles are resampled at the previous analysis time, all particles have the same weight ${\mathit{w}}^{\mathrm{b}\left(i\right)}=\mathrm{1}/m\left(i=\mathrm{1},\mathrm{\dots },m\right)$. The sum of the weights ${\sum }_{i=\mathrm{1}}^m{\mathit{w}}_t^{\mathrm{a},b\left(i\right)}$ is always 1. The denominator on the right-hand side of Eq. (9) can be estimated by

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathit{\pi }\left({\mathit{y}}_t^{\mathrm{o}}\mathrm{|}{\mathit{y}}_{\mathrm{1}:t-\mathrm{1}}^{\mathrm{o}}\right)& =\int \mathit{\pi }\left({\mathit{y}}_t^{\mathrm{o}}\mathrm{|}{\mathit{x}}_t\right)\mathit{\pi }\left({\mathit{x}}_t\mathrm{|}{\mathit{y}}_{\mathrm{1}:t-\mathrm{1}}^{\mathrm{o}}\right)\mathrm{d}{\mathit{x}}_t\\ & \approx {\sum }_{i=\mathrm{1}}^m{\mathit{w}}_t^{\mathrm{b}\left(i\right)}\mathit{\pi }\left({\mathit{y}}_t^{\mathrm{o}}\mathrm{|}{\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right).\end{array}.\end{array}$$

This study assumes a Gaussian likelihood function, given by

$$\begin{array}{}\text{(12)}& \mathit{\pi }\left({\mathit{y}}_t^{\mathrm{o}}\mathrm{|}{\mathit{x}}_t\right)\propto \exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{y}}_t^{\mathrm{o}}-H_t\left({\mathit{x}}_t\right)\right)}^T{\mathbf{R}}_t^{-\mathrm{1}}\left({\mathit{y}}_t^{\mathrm{o}}-H_t\left({\mathit{x}}_t\right)\right)\right].\end{array}$$

This assumption means the observation error distribution is assumed to be Gaussian. Then, the posterior can be expressed by

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathit{\pi }}_{\mathrm{LPF}}\left({\mathit{x}}_t\mathrm{|}{\mathit{y}}_{\mathrm{1}:t}^{\mathrm{o}}\right)\approx {\sum }_{i=\mathrm{1}}^m{\mathit{w}}_t^{\mathrm{a}\left(i\right)}\mathit{\delta }\left({\mathit{x}}_t-{\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right),\text{(14)}& {\displaystyle }\begin{array}{rl}{\mathit{w}}_t^{\mathrm{a}\left(i\right)}& ={\mathit{w}}_t^{\mathrm{b}\left(i\right)}{\mathit{\pi }}_{\mathrm{LPF}}\left({\mathit{y}}_{\mathrm{1}:t}^{\mathrm{o}}\mathrm{|}{\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right)/\\ & \left\{{\sum }_{k=\mathrm{1}}^m{\mathit{w}}_t^{\mathrm{b}\left(k\right)}{\mathit{\pi }}_{\mathrm{LPF}}\left({\mathit{y}}_{\mathrm{1}:t}^{\mathrm{o}}\mathrm{|}{\mathit{x}}_t^{\mathrm{b}\left(k\right)}\right)\right\}\\ & ={\mathit{w}}_t^{\mathrm{b}\left(i\right)}{\mathit{q}}_t^{\left(i\right)}/\left\{{\sum }_{k=\mathrm{1}}^m{\mathit{w}}_t^{\mathrm{b}\left(k\right)}{\mathit{q}}_t^{\left(k\right)}\right\},\end{array}\text{(15)}& {\displaystyle }{\mathit{q}}_t^{\left(i\right)}=\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{y}}_t^{\mathrm{o}}-H_t\left({\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right)\right)}^T{\mathbf{R}}_t^{-\mathrm{1}}\left({\mathit{y}}_t^{\mathrm{o}}-H_t\left({\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right)\right)\right],\end{array}$$

where q is the likelihood. Equation (14) results in posterior weights that satisfy ${\sum }_{i=\mathrm{1}}^m{\mathit{w}}_t^{\mathrm{a}\left(i\right)}=\mathrm{1}$.

To mitigate weight collapse, the local PF (LPF) solves the PF equations by assimilating local observations surrounding the analysis grid point. This study uses Penny and Miyoshi's (2016) approach that computes the analysis at every model grid point independently, with the observation error covariance R being inflated by the inverse of a localization function, as in the LETKF (Penny and Miyoshi, 2016).

Sampling importance resampling (SIR) is a technique for applying the PF for high-dimensional dynamics with a limited amount of particles. The resampling process rearranges the particles to effectively represent the densest areas of the posterior PDF. After resampling, each particle has equal weights, and the posterior PDF is given by

$$\begin{array}{}\text{(16)}& {\mathit{\pi }}_{\mathrm{LPF}}\left({\mathit{x}}_t\mathrm{|}{\mathit{y}}_{\mathrm{1}:t}^{\mathrm{o}}\right)\approx {\displaystyle \frac{\mathrm{1}}{m}}{\sum }_{i=\mathrm{1}}^m\mathit{\delta }\left({\mathit{x}}_t-{\mathit{x}}_t^{\mathrm{a}\left(i\right)}\right).\end{array}$$

The resampling process can be expressed as equally valid to the ensemble transform matrix of the ETKF (Reich, 2013; Penny and Miyoshi, 2016), given by

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathbf{X}}_{t,\mathrm{LPF}}^{\mathrm{a}}& ={\mathbf{X}}_t^{\mathrm{b}}{\mathbf{T}}_{t,\mathrm{LPF}}=\left({\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}+{\mathbf{Z}}_t^{\mathrm{b}}\right){\mathbf{T}}_{t,\mathrm{LPF}}\\ & ={\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}+{\mathbf{Z}}_t^{\mathrm{b}}{\mathbf{T}}_{t,\mathrm{LPF}},\end{array}\end{array}$$

where T_LPF denotes the ensemble transform matrix of the LPF. Here, we applied $\left({\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}\right){\mathbf{T}}_{t,\mathrm{LPF}}={\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}$ based on the following necessary condition of the ensemble transform matrix for the LPF:

$$\begin{array}{}\text{(18)}& {\sum }_{i=\mathrm{1}}^m{\mathbf{T}}_{t,\mathrm{LPF}}^{(i,j)}=\mathrm{1}j=\mathrm{1},\mathrm{\dots },m.\end{array}$$

In addition to Eq. (18), the ideal resampling matrix satisfies the following two conditions:

• $$\begin{array}{}\text{(19)}& \begin{array}{rl}{\sum }_{j=\mathrm{1}}^m{\mathbf{T}}_{t,\mathrm{LPF}}^{(i,j)}& =m\cdot {\mathit{w}}_t^{\mathrm{a}\left(i\right)}\text{ for }i=\mathrm{1},\mathrm{\dots }m\\ & \iff {\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}={\sum }_{i=\mathrm{1}}^m{\mathit{x}}_t^{\mathrm{b}\left(i\right)}\cdot {\mathit{w}}_t^{\mathrm{a}\left(i\right)},\end{array}\end{array}$$

• spatially smooth transition of T_t,LPF,

where m is the ensemble size, and T^(i,j) indicates the ith-row and jth-column element of the matrix T.

The resampling matrix significantly affects the filter performance (Farchi and Bocquet, 2018). The present study constructs resampling matrices based on the Monte Carlo approach. This resampling method uses m random numbers r⁽ⁱ⁾ ($i=\mathrm{1},\mathrm{\dots },m$) drawn from uniform distribution U ([0,1]) and accumulated weights w_t,acc:

$$\begin{array}{}\text{(20)}& {\mathit{w}}_{t,\mathrm{acc}}^{\left(\mathrm{0}\right)}=\mathrm{0},{\mathit{w}}_{t,\mathrm{acc}}^{\left(i\right)}={\mathit{w}}_{t,\mathrm{acc}}^{(i-\mathrm{1})}+{\mathit{w}}_t^{\mathrm{a}\left(i\right)}i=\mathrm{1},\mathrm{\dots },m.\end{array}$$

By definition, ${\mathit{w}}_{t,\mathrm{acc}}^{\left(m\right)}$ is 1. After sorting r to be the ascending order r_sorted, we generate the resampling matrix using Algorithm 1. This procedure is similar to the resampling of Potthast et al. (2019), but it employs an additional treatment so that the resampling matrix is close to the identity matrix (line 10 of Algorithm 1; Fig. 1b and c).

Algorithm 1 Generation of the resampling matrix T.

Require:

r_sorted and w_t,acc

1:

All components of T and z are initialized to be 0 # z: an indicator of resampled particles

2:

do j=1 to m # index of posterior particle

3:

do i=1 to m # index of prior particle

4:

if   $\mathit{\{}{\mathit{w}}_{t,\mathrm{acc}}^{(i-\mathrm{1})}<r_{\mathrm{sored}}^{\left(j\right)}\le {\mathit{w}}_{t,\mathrm{acc}}^{\left(i\right)}\mathit{\}}$ then

5:

z^(j)=i # resampled particle

6:

end if

7:

end do

8:

end do

9:

do j=1 to m # (loop for posterior particle)

10:

k=z^(j)

11:

if   $\mathit{\{}{\mathbf{T}}^{(k,k)}=\mathrm{0}\mathit{\}}$ then

12:

${\mathbf{T}}^{(k,k)}=\mathrm{1}$ # diagonal components

13:

z^(j)=0

14:

end if

15:

end do

16:

do j=1 to m # (loop for posterior particle)

17:

k=z^(j)

18:

if {k≠0} then

19:

do l=1 to m # (loop for columns of T)

20:

if { sum(${\mathbf{T}}^{(:,l)})=\mathrm{0}\mathit{\}}$ then

21:

${\mathbf{T}}^{(k,l)}=\mathrm{1}$ # off-diagonal components

22:

z^(j)=0

23:

exit

24:

end if

25:

end do

26:

end if

27:

end do

Figure 1 Examples of the resampling matrix in the case of 40 particles with a given weight generated by uniform random numbers. (a) Weight for 40 particles. (b, c) Examples of sampled resampling matrices by Algorithm 1 using different random numbers r. (d, e) Resampling matrices using the Monte Carlo stochastic approach with 200 and 10 000 sampled matrices, respectively.

Classical resampling can produce the same particles, resulting in the loss of diversity of posterior particles. Therefore, additional treatments are required to generate slightly different posterior particles. Potthast et al. (2019) proposed adding Gaussian random noise to T_t,LPF to prevent the generation of the same particles. As an alternative, the present study uses the Monte Carlo approach that repeats Algorithm 1 many times with different random numbers and takes the average of the generated resampling matrices (Fig. 1). Using different random numbers r, the resampling matrices differ (Fig. 1b and c), even if the same weight is used (Fig. 1a). By averaging the resampling matrices generated by using different random numbers, we get matrices that have higher weights in the diagonal components and lower weights in the off-diagonal components (Fig. 1d and e). The Monte Carlo approach makes the resampling matrix closer to the identity matrix, which is beneficial for a smooth transition in space. In addition, this stochastic approach approximates Eq. (19) using the Monte Carlo approach. The generated transform matrix with 200 samples (Fig. 1c) is close to that with 10 000 samples (Fig. 1d) in the case of 40 particles. Therefore, the transform matrices are generated by averaging 200 sampled matrices in subsequent LPF experiments. Hereafter, we call this resampling method “MC resampling”, which is used in the following experiments. The number of required samples for MC resampling is briefly investigated in Appendix B.

The effective particle (or ensemble) size N_eff (Kong et al., 1994) is useful to measure the diversity of particles in LPF:

$$\begin{array}{}\text{(21)}& N_{t,\mathrm{eff}}=\mathrm{1}/{\sum }_{i=\mathrm{1}}^m{\left({\mathit{w}}_t^{\mathrm{a}\left(i\right)}\right)}^{\mathrm{2}}\in \left[\mathrm{1},m\right].\end{array}$$

If N_eff is sufficiently large (N_eff≈m), no resampling is needed. This study considers a tunable parameter N₀ as a criterion for resampling:

$$\begin{array}{}\text{(22)}& {\mathbf{T}}_{t,\mathrm{LPF}}\leftarrow \left\{\begin{array}{ll}\mathbf{I}& \text{ if }{N}_{t,\mathrm{eff}}>N_{\mathrm{0}}\\ {\mathbf{T}}_{t,\mathrm{LPF}}& \text{ else}\end{array}\right.\end{array}$$

Without resampling ($N_{t,\mathrm{eff}}>N_{\mathrm{0}}$), the posterior weight is succeeded for the subsequent forecast as follows:

$$\begin{array}{}\text{(23)}& {\mathit{w}}_{t+\mathrm{1}}^{\mathrm{b}}\leftarrow \left(\mathrm{1}-\mathit{\tau }\right)\cdot {\mathit{w}}_t^{\mathrm{a}}+\mathit{\tau }/m,\end{array}$$

where τ is the tunable forgetting factor $\mathit{\tau }\in [\mathrm{0}:\mathrm{1}]$. Here, τ=1 means the subsequent prior weight ${\mathit{w}}_{t+\mathrm{1}}^{\mathrm{b}}$ has the same weights $\mathrm{1}/m$, whereas the posterior weight ${\mathit{w}}_t^{\mathrm{a}}$ is completely succeeded to ${\mathit{w}}_{t+\mathrm{1}}^{\mathrm{b}}$ when τ=0. This weight succession (Eq. 23) can be interpreted as a temporal localization that reduces the impact of observations temporally distant from the assimilation time. The weight succession in local PF is not trivial, because the weight (or likelihood) of particles would move with the dynamical flow. Similar discussions can be found for the advection of localization functions (e.g., Ota et al., 2013; Kotsuki et al., 2019a). This study assumes no flow motion of the weights and simply uses Eq. (23) at each grid point independently.

Although the LETKF applies inflation to the prior error covariance (i.e., P^b→βP^b), this inflation is suboptimal to the LPF, partly because the weight collapse occurs more easily if multiplicative inflation is applied to the prior perturbation (see Eq. 15). Therefore, the LPF usually applies inflation to the posterior particles (e.g., Penny and Miyoshi, 2016; Farchi and Bocquet, 2018). The present study uses the RTPS to inflate posterior perturbations for the LPF. We do not use the RTPS for the LETKF experiments, since the Miyoshi (2011) adaptive multiplicative inflation is known to outperform the manually tuned RTPS for idealized twin experiments with SPEEDY, based on the authors' preliminary experiments and Yoichiro Ota (personal communication, 2021).

Preliminary experiments showed that the RTPS outperformed the relaxation to prior perturbation (RTPP; Zhang et al., 2004) for the LPF. Having some noise in the transform matrix is important for the LPF to maintain the diversity of posterior particles. The RTPP makes the transform matrix closer to the identity matrix, resulting in less diverged posterior particles. Therefore, the RTPS would be a more suitable relaxation method than the RTPP for LPF.

2.3 LPF with Gaussian mixture extension

For the LPF to avoid weight collapse, assimilating too many independent observations in the local area is not desirable (van Leeuwen et al., 2019). To solve this problem, hybrid algorithms of EnKF and PF have been explored to efficiently assimilate massive observation data. The Gaussian mixture extension of the LPF is one such hybrid algorithm (Hoteit et al., 2008; Stordal et al., 2011; WP22). In the Gaussian mixture extension, the prior PDF is approximated by a combination of Gaussian distributions centered at the values of the particles, given by:

$$\begin{array}{}\text{(24)}& {\mathit{\pi }}_{\mathrm{LPFGM}}\left({\mathit{x}}_t\mathrm{|}{\mathit{y}}_{\mathrm{1}:t-\mathrm{1}}^{\mathrm{o}}\right)\approx {\sum }_{i=\mathrm{1}}^m{\mathit{w}}_t^{\mathrm{b}\left(i\right)}N\left({\mathit{x}}_t^{\mathrm{b}\left(i\right)},{\widehat{\mathbf{P}}}_t^{\mathrm{b}}\right),\end{array}$$

where $N\left({\mathit{x}}_t^{\mathrm{b}\left(i\right)},{\widehat{\mathbf{P}}}_t^{\mathrm{b}}\right)$ is the Gaussian kernel, with mean ${\mathit{x}}_t^{\mathrm{b}\left(i\right)}$ and covariance ${\widehat{\mathbf{P}}}_t^{\mathrm{b}}$. Here, hat indicates matrices for the Gaussian kernels (e.g., prior error covariance ${\widehat{\mathbf{P}}}^{\mathrm{b}}$, observation error covariance $\widehat{\mathbf{R}}$, and Kalman gain $\widehat{\mathbf{K}})$. The covariance of the Gaussian kernel uses the sampled covariance matrix ${\mathbf{P}}_t^{\mathrm{b}}$, such that

$$\begin{array}{}\text{(25)}& {\widehat{\mathbf{P}}}_t^{\mathrm{b}}=\mathit{\gamma }{\mathbf{P}}_t^{\mathrm{b}}\approx {\displaystyle \frac{\mathit{\gamma }}{m-\mathrm{1}}}{\mathbf{Z}}_t^{\mathrm{b}}{\left({\mathbf{Z}}_t^{\mathrm{b}}\right)}^T,\end{array}$$

where γ (>0) is a tunable parameter that regulates the amplitude and width of the Gaussian kernel (i.e., uncertainty of particles). For example, larger γ reduces amplitude and increases the width of the Gaussian kernel. Because the kernel is supposed to have a Gaussian distribution, increasing the Gaussian kernel's width results in a decrease in amplitude. The LPFGM results in the same analysis of the LPF when γ is 0 (γ→0).

The Gaussian mixture performs a two-step update to obtain the posterior particles. The first update moves the center of the Gaussian kernel with observations. Since kernels are Gaussian, we can use the Kalman filter for the first update:

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\mathit{x}}_t^{\mathrm{a}\left(i\right)}& ={\mathit{x}}_t^{\mathrm{b}\left(i\right)}+{\widehat{\mathbf{K}}}_t\left({\mathit{y}}_t^{\mathrm{o}}-H\left({\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right)\right)(i=\mathrm{1},\mathrm{\dots },m)\\ & \iff {\mathbf{X}}_{t,\mathrm{GM}}^{\mathrm{a}}={\mathbf{X}}_t^{\mathrm{b}}+{\widehat{\mathbf{K}}}_t\left({\mathit{y}}_t^{\mathrm{o}}\cdot \mathbf{1}-H_t\left({\mathbf{X}}_t^{\mathrm{b}}\right)\right),\end{array}\end{array}$$

where $\widehat{\mathbf{K}}={\widehat{\mathbf{P}}}^{\mathrm{a}}{\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}$, as in Eq. (5). In the LETKF, the Kalman gain K is computed by $\mathbf{K}={\mathbf{P}}^{\mathrm{a}}{\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}\approx {\mathbf{Z}}^{\mathrm{b}}{\widehat{\mathbf{P}}}^{\mathrm{a}}{\left({\mathbf{Z}}^{\mathrm{b}}\right)}^{\mathbf{T}}{\mathbf{H}}^T{\mathbf{R}}^{-\mathrm{1}}$. Since the Gaussian kernel uses the ensemble-based error covariance (Eq. 25), we can apply the exact algorithms of the LETKF to compute ${\widehat{\mathbf{P}}}^{\mathrm{a}}$ by replacing β of Eq. (3) with γ (i.e., ${\widehat{\mathbf{P}}}^{\mathrm{a}}={\left[\frac{\left(m-\mathrm{1}\right)}{\mathit{\gamma }}\mathbf{I}+{\left({\mathbf{Y}}_t^{\mathrm{b}}\right)}^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{Y}}_t^{\mathrm{b}}\right]}^{-\mathrm{1}})$. Here, the same gain matrix $\widehat{\mathbf{K}}$ is applied to update each particle independently, while the gain matrix $\widehat{\mathbf{K}}$ is based on the forecast error covariance estimated from the entire ensemble. Consequently, Eq. (26) is equivalent to

$$\begin{array}{}\text{(27)}& {\mathbf{X}}_{t,\mathrm{GM}}^{\mathrm{a}}={\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}+{\mathbf{Z}}_t^{\mathrm{b}}{\mathbf{T}}_{t,\mathrm{GM}},\end{array}$$

where

$$\begin{array}{}\text{(28)}& {\mathbf{S}}_t={\widehat{\mathbf{P}}}_t^{\mathrm{a}}{\left({\mathbf{Y}}_t^{\mathrm{b}}\right)}^T{\mathbf{R}}_t^{-\mathrm{1}}\left({\mathit{y}}_t^{\mathrm{o}}\cdot \mathbf{1}-H_t\left({\mathbf{X}}_t^{\mathrm{b}}\right)\right),\end{array}$$

and

$$\begin{array}{}\text{(29)}& {\mathbf{T}}_{t,\mathrm{GM}}={\mathbf{S}}_t+\mathbf{I}.\end{array}$$

Here, T_t,GM denotes the ensemble transform matrix of the Gaussian mixture.

The second update resamples the particles based on the likelihood of the posterior kernel, given by

$$\begin{array}{}\text{(30)}& {\mathit{q}}_t^{\left(i\right)}=\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\mathit{y}}_t^{\mathrm{o}}-H_t\left({\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right)\right)}^T{\widehat{\mathbf{R}}}_t^{-\mathrm{1}}\left({\mathit{y}}_t^{\mathrm{o}}-H_t\left({\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right)\right)\right],\end{array}$$

where

$$\begin{array}{}\text{(31)}& {\widehat{\mathbf{R}}}_t={\mathbf{R}}_t+{\mathbf{H}}_t{\widehat{\mathbf{P}}}_t^{\mathrm{b}}{\mathbf{H}}_t^T.\end{array}$$

Hoteit et al. (2008) and Stordal et al. (2011) used Eq. (30) for computing the likelihood of posterior kernels. Alternatively, WP22 suggested that using Eq. (15) instead of Eq. (30) is a reasonable approximation in the case of a smaller variance of ${\widehat{\mathbf{P}}}_t^{\mathrm{b}}$ compared to the observation departure ${\mathit{y}}_t^{\mathrm{o}}-H_t\left({\mathit{x}}_t^{\mathrm{b}\left(i\right)}\right)$. This study follows the WP22's approximation, because computing the inverse of Eq. (31) is computationally much more expensive than computing the inverse of diagonal R. Using WP22's approximation, the solution of the LPFGM is given by

$$\begin{array}{}\text{(32)}& \begin{array}{rl}{\mathbf{X}}_{t,\mathrm{LPFGM}}^{\mathrm{a}}& ={\mathbf{X}}_{t,\mathrm{GM}}^{\mathrm{a}}{\mathbf{T}}_{t,\mathrm{LPF}}\\ & =\left({\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}+{\mathbf{Z}}_t^{\mathrm{b}}{\mathbf{T}}_{t,\mathrm{GM}}\right){\mathbf{T}}_{t,\mathrm{LPF}}\\ & ={\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}+{\mathbf{Z}}_t^{\mathrm{b}}{\mathbf{T}}_{t,\mathrm{GM}}{\mathbf{T}}_{t,\mathrm{LPF}}.\end{array}\end{array}$$

Here, we used $\left({\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}\right){\mathbf{T}}_{t,\mathrm{LPF}}={\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}\cdot \mathbf{1}$ (cf. Eq. 18). Consequently, the transform matrix of the LPFGM (T_LPFGM) is given by

$$\begin{array}{}\text{(33)}& {\mathbf{T}}_{t,\mathrm{LPFGM}}={\mathbf{T}}_{t,\mathrm{GM}}{\mathbf{T}}_{t,\mathrm{LPF}}.\end{array}$$

Namely, the LPFGM can be described as the ensemble transform matrix form. Representing the LPFGM with only one transform matrix T_LPFGM is practically beneficial if one aims to reduce computational costs by the weight interpolation for the LETKF, in which the transform matrices at higher-resolution model grid points are interpolated by transform matrices at coarser model grid points (Yang et al., 2009; Kotsuki et al., 2020). The weight interpolation is also useful for ensuring a spatially smooth transition of the transform matrix for LPFs (Potthast et al., 2019).

Figure 2 A simple scalar example of observation as well as prior and posterior PDFs for (a) LPF, (b) fifth particle of LPFGM, and (c) LPFGM. Green, blue, red, purple, and dashed circles represent an observation and prior particles, moved particles, posterior particles, and replaced particles, respectively.

The two-step update of the LPFGM may appear to use the same observations twice, but this is not true. To understand the principles here, we consider a simple scalar example with H=1.0, with illustrations in Fig. 2. Let $\mathit{\pi }\left(y_t^{\mathrm{o}}\right)$ be an observation PDF (Fig. 2, top row), and the prior and posterior PDFs of the LPF are given by ${\mathit{\pi }}_{\mathrm{LPF}}\left(x_t\right)\approx \frac{\mathrm{1}}{m}{\sum }_{i=\mathrm{1}}^m\mathit{\delta }\left(x_t-x_t^{\mathrm{b}\left(i\right)}\right)$ and ${\mathit{\pi }}_{\mathrm{LPF}}\left(x_t\mathrm{|}{y}_t^{\mathrm{o}}\right)$ ,respectively (Fig. 2a, second and third rows). The LPF employs resampling by approximating the posterior PDF as a combination of prior particles, such that ${\mathit{\pi }}_{\mathrm{LPF}}\left(x_t\mathrm{|}{y}_t^{\mathrm{o}}\right)\approx {\sum }_{i=\mathrm{1}}^m{w}_t^{\mathrm{a}\left(i\right)}\mathit{\delta }\left(x_t-x_t^{\mathrm{b}\left(i\right)}\right)$ (Fig. 2a, bottom row). Next, we focus on the fifth particle of the LPFGM (Fig. 2b). Prior and posterior PDFs of the fifth particle are given by $\mathit{\pi }\left(x_t^{\left(\mathrm{5}\right)}\right)=N\left(x_t^{\mathrm{b}\left(\mathrm{5}\right)},{\widehat{P}}_t^{\mathrm{b}}\right)$ and $\mathit{\pi }\left(x_t^{\left(\mathrm{5}\right)}\mathrm{|}{y}_t^{\mathrm{o}}\right)=N\left(x_t^{\mathrm{b}\left(\mathrm{5}\right)},{\widehat{P}}_t^{\mathrm{b}}\right)\cdot N\left(y_t^{\mathrm{o}},R_t\right)$, respectively (Fig. 2b, second and third rows). Since the Gaussian kernels are assumed for the prior particles, the center of the posterior kernel moves such that $\mathit{\pi }\left(x_t^{\left(\mathrm{5}\right)}\mathrm{|}{y}_t^{\mathrm{o}}\right)\propto N\left(x_t^{\mathrm{a}\left(\mathrm{5}\right)},{\widehat{P}}_t^{\mathrm{a}}\right)$, where $x_t^{\mathrm{a}\left(\mathrm{5}\right)}$ and ${\widehat{P}}_t^{\mathrm{a}}$ can be computed by the Kalman filter (from the blue circle to the red circle in Fig. 2b). Since the LPFGM moves all particles, the posterior PDF of the LPFGM is given by ${\mathit{\pi }}_{\mathrm{LPFGM}}\left(x_t\mathrm{|}{y}_t^{\mathrm{o}}\right)\approx \frac{\mathrm{1}}{m}{\sum }_{i=\mathrm{1}}^m\mathit{\pi }\left(x_t^{\left(i\right)}\mathrm{|}{y}_t^{\mathrm{o}}\right)$ (Fig. 2c, third row). These movements correspond to the first update of the LPFGM. In contrast to the LPF, the posterior PDF of the LPFGM is approximated by a combination of the posterior Gaussian kernels (red circles of Fig. 2c, third row). The LPFGM employs the resampling based on the moved particles such that ${\mathit{\pi }}_{\mathrm{LPFGM}}\left(x_t\mathrm{|}{y}_t^{\mathrm{o}}\right)\approx {\sum }_{i=\mathrm{1}}^m{w}_t^{\mathrm{a}\left(i\right)}N\left(x_t^{\mathrm{a}\left(i\right)},{\widehat{P}}_t^{\mathrm{a}}\right)$ (Fig. 2c, bottom row). This resampling corresponds to the second update of the LPFGM. As seen in this example, the two-step update of the LPFGM does not use the same observations twice.

Figure 3 Workflows of LETKF, LPFGM, and LPF. Gray, purple, green, and red boxes represent the components of LETKF, a similar component as in the LETKF, the components of the LPF, and a unique component of the LPFGM, respectively. Subroutines and source files are also described in boxes.

2.4 Implementation and computational complexity

We implemented the LPFGM based on the available LETKF code from Miyoshi (2005) and from follow-on studies (Miyoshi and Yamane, 2007; Kondo and Miyoshi, 2019; Kotsuki et al., 2020; https://github.com/takemasa-miyoshi/letkf, last access: 16 November 2022). Figure 3 compares the workflows of the LETKF, LPF, and LPFGM. All DA methods execute the same first four steps (steps A–D). After step D, the LETKF involves four additional steps (steps E–H) to compute T_LETKF. The LPF computes the weights of particles (step J), followed by the generation of the transform matrix T_LPF (step K). The LPFGM first executes steps E–G, as in the LETKF, and then executes step I to compute T_GM, followed by the LPF algorithms to compute T_LPF (steps J and K). Finally, the LPFGM multiplies T_GM and T_LPF to compute T_LPFGM (step L). At the end of the process, the transform matrix is applied to the prior perturbation matrix Z^b to produce the analysis ensemble (step M) followed by the RTPS (step N).

If the LETKF code is available, the LPF in the transform matrix form can be developed by coding two steps, J and K. The LPFGM can also be developed easily by coding two more steps, I and L, if the LETKF and LPF codes are available.

Table 1 Computational complexities of LETKF, LPF, and LPFGM. Each step corresponds to the steps in Fig. 3. Cross marks represent the steps used for LETKF, LPF, and LPFGM.

1: This computation depends on the localization scale
2: These computations assume the diagonal R (i.e., uncorrelated observation error)
n: system size
m: ensemble size
p: the number of observations
p_L: the number of local observations
N_MC: the number of times for generating resampling matrices by Algorithm 1

Here, we compare the computational complexities of the LETKF, LPF, and LPFGM algorithms (Table 1). The total cost (C_T) of a DA cycle is identical to the overhead of the assimilation system (C_H) plus n times the average local analysis cost (C_L) and m times the cost of a one-member model forecast (C_M):

$$\begin{array}{}\text{(34)}& C_{\mathrm{T}}\left(n,p,m\right)=C_{\mathrm{H}}\left(n,p,m\right)+n\cdot C_{\mathrm{L}}\left(p_{\mathrm{L}},m\right)+m\cdot C_{\mathrm{M}},\end{array}$$

where p_L is the number of local observations within the localization cut-off radius. We assume that the overhead and model costs are equivalent among the three methods. In addition, the total computational cost of DA usually depends on the local analysis cost (C_L):

$$\begin{array}{}\text{(35)}& {\displaystyle }{C}_{\mathrm{L}}^{\mathrm{LETKF}}=O\left(m^{\mathrm{3}}+m^{\mathrm{2}}{p}_{\mathrm{L}}\right),\text{(36)}& {\displaystyle }{C}_{\mathrm{L}}^{\mathrm{LPF}}=O\left(m^{\mathrm{2}}{N}_{\mathrm{MC}}\right),\text{(37)}& {\displaystyle }{C}_{\mathrm{L}}^{\mathrm{LPFGM}}=O\left(m^{\mathrm{3}}+m^{\mathrm{2}}{p}_{\mathrm{L}}+m^{\mathrm{2}}{N}_{\mathrm{MC}}\right),\end{array}$$

where N_MC is the number of times the resampling matrices are generated by Algorithm 1. The number of local observations p_L is usually much greater than the ensemble size m and N_MC. In this case, O(m²p_L) is dominant for the LETKF and the LPFGM.

Here, the computational cost of LPF is more expensive than that with a simpler resampling algorithm, such as the stochastic uniform resampling (O(m); Penny and Miyoshi 2016; Farchi and Bocquet, 2018), due to the relatively complex Algorithm 1 (O(m²N_MC)). We could not use such simpler approaches in this study, since they do not yield ideal resampling matrices that satisfy the two conditions (Eq. 19 and spatially smooth transition of T_t,LPF). The computational cost O(m²N_MC) of LPF is still much smaller than the LETKF if N_MC≪p_L.

3 Experimental settings

3.1 SPEEDY model

This study used the intermediate global atmospheric general circulation model SPEEDY (Molteni, 2003) to compare the LETKF, LPF, and LPFGM. The SPEEDY model is a computationally inexpensive hydrostatic model with fundamental physical parameterization schemes such as surface flux, radiation, convection, cloud, and condensation. The SPEEDY model has 96×48 grid points in the horizontal plane (T30 ∼3.75^∘ × 3.75^∘) and seven sigma-coordinate seven vertical layers. The SPEEDY model consists of five prognostic variables: temperature (T), specific humidity (Q), zonal wind (U), meridional wind (V) at seven layers, and surface pressure (P_s). The SPEEDY model coupled with the LETKF (SPEEDY-LETKF) has been widely used in DA studies (Miyoshi, 2005, and many follow-on studies, e.g., Miyoshi, 2011; Kondo et al., 2013; Kondo and Miyoshi, 2016, 2019; Kotsuki and Bishop, 2022). We implemented the LPF and LPFGM based on the existing SPEEDY−LETKF code, following the procedures described in Sect. 2.4.

Figure 4 Observing networks for (a) REG2 and (c) RAOB experiments. Small black dots and red crosses represent model grid points and observing points, respectively.

3.2 Experimental design

The experimental settings follow the previous SPEEDY−LETKF experiments (Miyoshi, 2005, and follow-on studies, e.g., Kotsuki et al., 2020). A series of idealized and identical experiments (also known as observing system simulation experiments) were conducted without model errors. We first performed a spin-up run for one year, initialized by the standard atmosphere during rest, followed by the nature run started at 00:00 UTC on 1 January of the second year. We assumed diagonal observation error covariance R (i.e., uncorrelated observation error). Gaussian noise was added to the nature run to produce observation data at 6 h intervals. The standard deviations of the Gaussian noise are 1.0 K for T, 1.0 m s⁻¹ for U and V, 0.1 g kg⁻¹ for Q, and 1 hPa for surface pressure. This study considered two observing networks (Fig. 4a and b): the regularly distributed network (hereafter REG2) and the radiosonde-like inhomogeneous network (hereafter RAOB). We observe T, U, and V at all seven layers, whereas Q was observed at the first–fourth layers. The ensemble size is 40, and their initial conditions were taken from an independent single deterministic SPEEDY forecast with sufficient spin-up simulation.

Table 2 List of experiments.

ρ_h: horizontal localization scale (Eq. 38)
α: relaxation parameter of RTPS (Eq. 8)
γ: amplitude of Gaussian kernel of the LPFGM (Eq. 25)
N₀: parameter controls resampling frequency (Eq. 22)
τ: forgetting factor of weight (Eq. 23)
1: Instead of RTPS, the LETKF uses Miyoshi (2011)'s adaptive multiplicative inflation.

Table 2 summarizes the settings of the LETKF, LPF, and LPFGM experiments. The observation error variance was inflated for the localization by using the Gaussian-based function given by

$$\begin{array}{}\text{(38)}& l=\left\{\begin{array}{l}\exp \left[-\frac{\mathrm{1}}{\mathrm{2}}\left\{{\left(d_{\mathrm{h}}/{\mathit{\rho }}_{\mathrm{h}}\right)}^{\mathrm{2}}+{\left(d_{\mathrm{v}}/{\mathit{\rho }}_{\mathrm{v}}\right)}^{\mathrm{2}}\right\}\right]\\ \text{ if }{d}_{\mathrm{h}}<\mathrm{2}\sqrt{\mathrm{10}/\mathrm{3}}{\mathit{\rho }}_{\mathrm{h}}\text{ and }{d}_{\mathrm{v}}<\mathrm{2}\sqrt{\mathrm{10}/\mathrm{3}}{\mathit{\rho }}_{\mathrm{v}}\\ \mathrm{0}\text{ else}\end{array}\right.\end{array}$$

where l is the localization function, and its inverse l⁻¹ is multiplied and used to inflate R. d_h and d_v are horizontal and vertical distances (km and log(Pa)) between the observation and analysis grid point. ρ_h and ρ_v are tunable horizontal and vertical localization scales (km and log(Pa)). Subscripts “h” and “v” represent horizontal and vertical respectively. This study set the vertical localization scale ρ_v to be 0.1 log(Pa) following Greybush et al. (2011) and Kondo et al. (2013). The horizontal localization scales of the LETKF were tuned manually before the experiments to minimize the first-guess root mean square error (RMSE) of the fourth layer (∼500 hPa) temperature of the SPEEDY-LETKF, since this is among the most important variables for medium-range NWP.

The LETKF experiments used the adaptive multiplicative inflation method of Miyoshi (2011), in which the inflation factor β was estimated adaptively. On the basis of sensitivity experiments for γ, we chose γ=1.5. Sensitivity to γ is discussed in Sect. 4.3.3.

We first performed a one-year (from January to December) SPEEDY-LETKF experiment over the second year following the one-year spin up. We then performed LETKF, LPF, and LPFGM experiments from January to April in the third year, initialized by the first-guess ensemble of the LETKF at 00:00 UTC on 1 January of the third year. The results from the last three months, i.e., February to April, were used for verification. We assessed the 6 h forecast or background RMSE for T at the fourth model level (∼500 hPa). While this study mainly discusses the RMSE and ensemble spread for T, similar results are observed for other variables. Here, this study will concentrate on the fourth-level T as an important variable for NWPs, since the major goal of this study is to investigate the stabilities of LPF and LPFGM compared with the LETKF. Humidity verification or accounting for nonlinear observation operators and non-Gaussian observation errors are still crucial studies to investigate the advantages of the LPF and LPFGM with respect to the LETKF.

Figure 5 Time-mean background RMSEs (solid lines) and ensemble spreads (dashed lines) for T (K) at the fourth model level (∼500 hPa) as a function of RTPS parameter α averaged over three months of the third year (February–April), with REG2 observations. Magenta, red, blue, and green lines are LPF experiments with horizontal localization scales of 700, 600, 500, and 400 km, respectively. Dashed lines represent RMSE of LETKF (0.1257 K), with adaptive multiplicative inflation instead of RTPS.

4 Results and discussion

4.1 Experiments with a regularly distributed observing network

We first compare the LETKF, LPF, and LPFGM with REG2, where the manually tuned horizontal localization scale for the LETKF is 600 km. First, sensitivities to the horizontal localization scale (ρ_h) and RTPS parameter (α) are investigated for the LPF. Figure 5 indicates the time-mean background RMSEs and ensemble spreads for LETKF and LPF. The LPF requires large inflation (α≥1.0) to avoid filter divergence. With large inflation (α≥1.0), the LPF experiment with ρ_h=600 km resulted in the smallest RMSEs among the three LPF experiments. The best-performing localization scale of LPF was found to be similar to that of the LETKF. The four LPF experiments exhibited larger RMSEs than the LETKF, as demonstrated by the previous studies with a low-dimensional L96 model (Poterjoy, 2016; Penny and Miyoshi, 2016; Farichi and Bouquet, 2018). The LPF experiment with ρ_h=700 km shows filter divergence when α=1.0. The ensemble spreads of LPF are smaller than the RMSEs when α=1.0. This under-dispersive ensemble is a typical behavior of LPF (e.g., Poterjoy and Anderson, 2016).

Figure 6 Time series of globally averaged background (a) RMSEs (solid lines) and ensemble spreads (dashed lines) for T (K) and (b) effective particle size N_eff at the fourth model level (∼500 hPa), with REG2 observations. Black lines show the LETKF. Red, blue, and green lines are the LPF experiments, with localization scales of 600, 500, and 400 km, respectively. RTPS parameter α is set to 1.00 in all three LPF experiments. The abscissa shows month and day.

Second, we compare the time series of the background RMSEs, ensemble spreads, and effective particle sizes N_eff (Fig. 6). The RMSE and ensemble spread are consistent in the LETKF. However, the LPF generally shows smaller ensemble spreads than the RMSEs. Since the beginning of the experiments on 1 January, the three LPF experiments showed rapid increases in the ensemble spread (Fig. 6a) and a rapid decrease in the effective particle size (Fig. 6b) within two weeks. After the rapid changes, the ensemble spreads and effective particle sizes were stabilized. The three LPF experiments increased RMSEs until the beginning of March. The LPF with ρ_h=500 km showed filter divergence in April. The LPF with ρ_h=400 km and ρ_h=600 km also showed increasing trends in RMSE until the end of April, while their ensemble spreads and effective particle sizes were stable.

Figure 7 Similar to Fig. 5 but for the LPFGM with REG2 and 600 km localization scale. Magenta, red, blue, and green lines show the cases with N₀=40, 10, 5, and 2, respectively. Dashed lines represent RMSE of LETKF (0.1257 K). Panels show the LPFGM with forgetting factors for (a) τ=0.0, (b) τ=0.1, (c) τ=0.2, (d) τ=0.3, (e) τ=0.5, and (f) τ=1.0, respectively.

Next, sensitivities to the RTPS parameter (α), resampling frequency (N₀), and forgetting factor (τ) are investigated for the LPFGM. The best performing localization scale of LPF was identical to that of the LETKF (Fig. 5). Therefore, we choose the same 600 km localization scale for the LPFGM. Figure 7 compares the time-mean background RMSEs and ensemble spreads for LETKF and LPFGM. In the LPFGM experiments, we investigate six experimental settings for τ=0.0, 0.1, 0.2, 0.3, 0.5, and 1.0. With τ=0.0 (Fig. 7a), the LPFGM shows RMSEs similar to those of LETKF, with best-performing parameters for α and N₀. Regulating the resampling frequency (N₀) would be needed for LPFGM, because the LPFGM shows filter divergence when resampling is employed at all assimilation steps, excluding the case of α=1.0 (magenta line in Fig. 7a). The LPFGM employs two ensemble updates: the Gaussian mixture and resampling. Owing to the first update by the Gaussian mixture, the LPFGM would not require resampling at all assimilation steps. With τ=0.0, the best-performing relaxation parameter α is approximately 0.6–0.9, excluding the case employing resampling at all assimilation steps.

Increasing τ leads to the LPFGM being less sensitive to the resampling parameter N₀ and RTPS parameter α (Fig. 7b–f), implying that the LPFGM is more stable when the LPFGM forgets weights (τ≥0.1). Additionally, the LPFGM requires less inflation (i.e., smaller α) when the forgetting factor τ becomes larger. Owing to the first update by the Gaussian mixture, succeeding weights completely may be unnecessary for the LPFGM. In addition, the method of succeeding weights in the local PF is not trivial, because the weights would move with the dynamical flow. Because the results of the LPFGM experiments without weight succession (τ=1.0) were superior to those with weight succession (τ=0.0), the remainder of this paper focuses on the LPFGM experiments without weight succession (τ=0.0).

.

Figure 8 Similar to Fig. 6, but showing the LPFGM experiments with REG2 and 600 km localization scale. Magenta, red, blue, and green lines are LPFGM experiments with N₀=40, 10, 5, and 2, respectively. Forgetting factor τ is 1.0. The RTPS parameter is set to 1.00 in the experiment with N₀=40 (i.e., all-time resampling) and to 0.50 for in other experiments.

The time series of background RMSEs and the effective particle sizes for LPFGM experiments show that the LPFGM with every-time resampling (N₀=40) exhibits large RMSEs and the smallest effective particle sizes (Fig. 8, magenta). In contrast, the LPFGM experiments with infrequent resampling (N₀=10, 5, and 2) show small RMSEs similar to that of the LETKF and maintain larger effective particle sizes than the LPF (Fig. 8, red, blue, and green lines).

Figure 9 Spatial patterns of background RMSE for T (K) at the fourth model level (∼500 hPa) for (a) LETKF, (b) LPF, and (c) LPFGM, with best-performing localization scale and RTPS parameter, averaged over February–April with REG2 and 600 km localization scale. Panels (d) and (e) show the differences between LETKF-LPF and LETKF-LPFGM, respectively. The warm (cold) color indicates that the LPF or LPFGM is better (worse) than the LETKF.

Finally, we compare the spatial patterns of the RMSEs at a fourth model level for the LETKF, LPF, and LPFGM, with best-performing parameter settings (Fig. 9). The LETKF and LPFGM show larger RMSEs in the tropics and polar regions (Fig. 9a and c), possibly because of uncertainties from convective dynamics in the tropics and sparser observations in the polar regions. The LPF shows a larger RMSE than the LETKF globally (Fig. 9b and d). In contrast, slight improvements were observed globally in the LPFGM relative to the LETKF, as indicated by generally warm colors in Fig. 9e, especially around the North Pole.

4.2 Experiments with a realistic observing network

Here, we compare the LETKF, LPF, and LPFGM with the realistic observing network RAOB. With this spatially inhomogeneous observing network, all LPF experiments showed filter divergence, even with a broad range of the localization and RTPS parameters. Since there are fewer assimilated observations in RAOB than in REG2, the weight collapse would not be the primary cause of the LPF filter divergence. Because the LPF creates posterior particles by linearly combining prior particles, it is preferable for the LPF to have observations within the range of the prior particles. Filter divergence must be prevented by maintaining synchronization between the LPF and the observations. The LPF would require more observations for the synchronization than the LETKF, which was also shown in the authors' initial tests with L96. In addition, the LPFGM was unstable with weight succession (τ=0.0). Therefore, this section focuses on the comparison of the LETKF and LPFGM without weight succession (τ=1.0) only.

Figure 10 Similar to Fig. 5 but showing the LPFGM experiments with the RAOB and 1100 km localization scale. Magenta, red, blue, and green lines are the LPFGM experiments with N₀=40, 10, 5, and 2, respectively. Forgetting factor τ is fixed at 1.0. Dashed black lines represent the RMSE of the LETKF (0.4762 K); (b) enlarges (a) for the range of the RMSE of 0.25–0.45.

First, sensitivities to the RTPS parameter (α) and resampling frequency (N₀) are investigated for the LPFGM. Figure 10 compares the time-mean background RMSEs and ensemble spreads for the LETKF and LPFGM. The best-performing horizontal localization scale for the LETKF is 1100 km, and the same localization scale is used for the LPFGM. The LPFGM with every-time resampling (N₀=40) shows the largest RMSEs (Fig. 10a, magenta). By contrast, some LPFGM experiments show smaller RMSEs than the LETKF (Fig. 10b). In this experimental setting, the resampling parameter of N₀=2.0 shows the most stable performance. With N₀=2.0, the LPFGM outperforms the LETKF slightly for α=0.5–0.7 (green line of Fig. 10b).

Figure 11 Similar to Fig. 8 but showing the LPFGM experiments with the RAOB and 1100 km localization scale. Magenta, red, blue, and green lines are LPFGM experiments with N₀=40, 10, 5, and 2, respectively. Forgetting factor τ is 1.0. The RTPS parameter is set to 1.00 in the experiment with N₀=40 (i.e., all-time resampling) and to 0.60 for in other experiments.

The time series of the background RMSEs and effective particle sizes show that the LPFGM with every-time resampling (N₀=40; Fig. 11a, magenta) and relatively frequent resampling (N₀=10; Fig. 11a, red) caused filter divergences. The times when the filters diverged seem to correspond to the times when the effective particle sizes decreased (Fig. 11b, magenta and red). This reduction in effective particle size is a typical sign of filter divergence for PF (Snyder et al., 2008). One may speculate that the assumption of ${\widehat{\mathbf{R}}}_t\approx {\mathbf{R}}_t$ in the LPFGM (Sect. 2.3) would be a reason for the filter divergence. However, adopting the appropriate norm (${\widehat{\mathbf{R}}}_t={\mathbf{R}}_t+{\mathbf{H}}_t{\widehat{\mathbf{P}}}_t^{\mathrm{b}}{\mathbf{H}}_t^T$) for the LPFGM did not eliminate the filter divergence (not shown). An effective solution for avoiding the filter divergence with every-time resampling is still unknown, but a possible reason is the method of creating a transform matrix using the Monte Carlo method that approximates two ideal conditions (Eq. 19 and a spatially smooth transition). For deeper understanding, using other resampling techniques such as the optimal transport (see Sect. 4.3.2) would be useful. In contrast, the LPFGM experiments with infrequent resampling (N₀=5 and 2) show stable RMSEs and maintain a larger effective particle size. These stable LPFGM experiments maintain the amplitude of the ensemble spreads similarly to that of the RMSEs.

Figure 12 Spatial patterns of background RMSE for T (K) at the fourth model level (∼500 hPa) for (a) LETKF and (b) LPFGM, averaged over February–April, with RAOB and 1100 km localization scale. Panel (c) shows the difference between LETKF and LPFGM. Warm (cold) color indicates that the LPFGM is better (worse) than the LETKF. Black dashed rectangles show the region (120–60^∘ W and 70–30^∘ S) where the differences between LETKF and LPFGM are investigated in Fig. 14.

We further explore regions where the LPFGM outperforms the LETKF. Figure 12 compares the space-based patterns of the background RMSEs for T (K) at the fourth model level for LETKF and LPFGM. Errors in sparsely recognized regions – such as the South Pacific Ocean, Indian Ocean, and polar regions – tend to be larger (Fig. 12a and b). The LPFGM outperforms the LETKF in such sparsely observed regions (Fig. 12c). The regions showing improvements correspond to the regions where the ensemble spreads and first-guess non-Gaussianity are larger (Fig. 13). Here, we measure the non-Gaussianity by means of Kullback–Leibler divergence (Kullback and Leibler, 1951).

Figure 13 Spatial distributions of the background (a, b) Kullback–Leibler divergence and (c, d) ensemble spread for T (K) at the fourth model level (∼500 hPa), averaged over February–April, with RAOB and 1100 km localization scale; (a, c) and (b, d) show the LETKF and LPFGM, respectively. Black dashed rectangles show the region (120–60^∘ W and 70–30^∘ S) where the differences between LETKF and LPFGM are investigated in Fig. 14.

The first update of the LPFGM (the Gaussian mixture) is similar to the LETKF update, as discussed. Therefore, the improvements of the LPFGM with respect to the LETKF would be owing to the resampling step. Previous studies demonstrated that the LPF outperformed the LETKF for non-Gaussian DA with non-Gaussian observation and prior PDFs (Poterjoy and Anderson, 2016; Poterjoy, 2016; Penny and Miyoshi, 2016). Therefore, the improvements of the LPFGM would come from the consideration of non-Gaussian prior PDF in sparsely observed regions.

Figure 14 Time series of 3 d-mean (a) RMSE (K), (b) ensemble spread (K), (c) Kullback–Leibler divergence, and (d) effective particle size N_eff at fourth model level temperature T, averaged over the region indicated by rectangles in Figs. 12 and 13 (120–60^∘ W and 70–30^∘ S). The black line represents the LETKF. Green and purple lines are the LPFGM whose resampling frequencies N₀ are 2 and 0, respectively. Vertical red lines represent the cases when the LETKF shows large RMSE greater than 0.8 K. In (d), the effective particle size is also computed for the LETKF using the first-guess ensemble.

Finally, we investigate the region where the difference in RMSE between the LETKF and LPFGM is large (120–60^∘ W and 70–30^∘ S; indicated by black dashed rectangles in Figs. 12 and 13). Figure 14 compares the time series of RMSE, ensemble spread, effective particle size, and Kullback–Leibler divergence averaged over the region for the LETKF and LPFGM. Here, we conducted an additional LPFGM experiment that employs no resampling (i.e., N₀=0) to investigate the importance of resampling. Vertical red lines in Fig. 14 represent the cases when the LETKF has large RMSEs greater than 0.8 K. Figure 14a indicates that the LPFGM tends to mitigate large RMSEs, which is in contrast to the LETKF. Since no significant increase is seen in the Kullback–Leibler divergence at the four cases (Fig. 14c), the large RMSEs of LETKF would not be caused by analyses with highly non-Gaussian first-guess ensembles. For the last three cases, the LPFGM has significantly smaller RMSEs than the LETKF. The improvement of the LPFGM would be partially led by larger ensemble spread than the LETKF (Fig. 14b). However, the LPFGM with resampling (N₀=2) outperforms the LPFGM without resampling (N₀=0), indicating that resampling improves the RMSE. The LETKF and LPFGM exhibit reduced effective particle size at the three later cases (Fig. 14d). Because of the smaller effective particle size, fewer ensembles predict states that are closer to observations, and the LPFGM employed more resampling in the region. In addition to the particle shift, the LPFGM would reduce RMSE further by the resampling.

The LPFGM occasionally resulted in larger RMSE than the LETKF (e.g., beginning of March in Fig. 14a). However, the LPFGM has the potential to outperform the LETKF overall by reducing the RMSE in sparsely observed regions. While operational NWP systems assimilate massive observations, the LPFGM would be useful when spatially and temporally sparse observations are used – e.g., in the twentieth-century reanalysis projects (e.g.,  Compo et al., 2011; Laloyaux et al., 2018) and paleoclimate reconstructions (Acevedo et al., 2017; Okazaki and Yoshimura, 2017).

4.3 Factors requiring further investigation

Finally, we discuss other issues for potential further improvement of the LPF and LPFGM in future studies.

4.3.1 Inflation

This study used the RTPS to inflate the posterior perturbation for the LPF and LPFGM. As shown in Figs. 5, 7, and 10, the RMSEs of the LPF and LPFGM are sensitive to the RTPS parameter. There is a need to investigate methods that estimate the RTPS parameter adaptively, as in the EnKF (Ying and Zhang, 2015; Kotsuki et al., 2017b). However, the adaptive relaxation methods used in the EnKF cannot be applied directly to the LPF, because they use the innovation statistics (Desroziers et al., 2005) that assume analysis updates by the Kalman gain. Therefore, substantially different adaptive relaxation methods should be investigated for the LPF and LPFGM.

An alternative method of inflation is adding random noise to the transform matrix (Potthast et al., 2019). Regulating the amplitude of the random noise was not trivial in the authors' preliminary experiments with L96 (not shown). Random noise that is too small results in a loss of diversity of posterior particles, whereas excessively large random noise results in an overly dispersive ensemble. Potthast et al. (2019) also suggested estimating the amplitudes of the random noise based on the innovation statistics. However, determining the minimum and maximum values of the amplitudes was not trivial from the authors' experience with L96. Therefore, a method to determine the optimum amplitude of random noise should be investigated as well. Moreover, additive inflation methods, as used in the EnKF, can be beneficial for producing diverse prior particles (Mitchell and Houtekamer, 2000; Corazza et al., 2003). For example, Penny and Miyoshi's (2016) LPF applied an additive inflation Gaussian noise whose variance is scaled to a magnitude of local analysis error variance. Investigating better and adaptive inflation methods is very important to stabilize LPFs.

4.3.2 Transform matrix for LPF

The transform matrix significantly affects the filter performance. Farchi and Bocquet (2018) demonstrated that the optimal transport method of Reich (2013) resulted in lower RMSE than the commonly used stochastic uniform resampling approach with the L96 model. However, it is not clear if optimal transport is beneficial for high-dimensional systems such as the SPEEDY model.

In this study, all-time resampling was detrimental for the LPFGM experiments (Figs. 7 and 10).This may be due to the MC resampling that does not exactly but only approximately satisfies Eq. (19) (${\stackrel{\mathrm{‾}}{\mathit{x}}}_t^{\mathrm{b}}={\sum }_{i=\mathrm{1}}^m{\mathit{x}}_t^{\mathrm{b}\left(i\right)}\cdot {\mathit{w}}_t^{\mathrm{a}\left(i\right)}$). With a small effective particle size, the MC resampling almost satisfies Eq. (19), since the resampling is not affected by the sampling noise in the selection of particles. In contrast, satisfying Eq. (19) is more difficult when the effective particle size is larger due to the sampling noise. A possible solution to this problem is to use the optimal transport method (Reich, 2013) that satisfies Eq. (19) exactly. Another essential property for local PFs is the spatially smooth transition of the transform matrix. The optimal transport method would also be useful for ensuring the spatially smooth transition of the transform matrix, because the posterior weights in nearby grids are typically similar. It is necessary to investigate better methods for generating resampling matrices.

4.3.3 Tunable parameters of LPF and LPFGM

The tunable parameters of the LPF and LPFGM should be investigated further. For example, parameter N₀, which controls the resampling frequency, significantly affected the filter accuracy and stability. Adaptive determination of this parameter can prevent time-consuming parameter tuning.

Figure 15 Time-mean background RMSEs for T (K) at the fourth model level (∼500 hPa) as a function of RTPS parameter α, averaged over three months of the third year (February–April); (a) shows REG2 with 600 km localization scale, and (b) shows RAOB with 1100 km localization scale. Red, blue, and green lines are LPFGM experiments where γ=2.5, 1.5, and 0.5, respectively. Black dashed lines represent RMSE of LETKF (0.1257 K) with adaptive multiplicative inflation instead of RTPS. Other parameters are set to be the best-performing parameters: N₀=2 and τ=1.0.

The sensitivity to parameter γ should also be investigated for the LPFGM. WP22 proposed using γ=2.5. Since γ controls the amplitude and width of Gaussian kernels, the filter accuracy and stability of the LPFGM would be sensitive to this parameter. We briefly examined the sensitivity to this parameter with REG2 and RAOB (Fig. 15; γ=2.5, 1.5, and 0.5 for red, blue, and green lines). The results indicate that the filter accuracy of the LPFGM is sensitive to this parameter, especially for a spatially inhomogeneous RAOB observing network. For our experimental settings, the best-performing γ was 1.5 for both REG2 and RAOB. However, the optimal parameter would differ for different models, observing networks, and observing frequency. Also, there might be a relation between the optimal γ and ensemble size. Therefore, adaptive determination of this parameter is helpful.

The LPFGM without weight succession (τ=1.0) resulted in lower RMSE than that with weight succession (τ=0.0). However, the optimal parameter τ would be somewhere between 0.0 and 1.0. In addition, the weight succession at each model grid point would be suboptimal. Hence, the method for weight succession should also be explored further.

5 Summary

This study aims to develop a software platform for the LETKF, LPF, and LPFGM with the intermediate global atmospheric model SPEEDY. The main results of this investigation are briefly listed as follows:

1. The LPF and LPFGM were developed with only minor modifications to the existing LETKF system with the SPEEDY model.

2. With dense observations (REG2), the LPF showed stable filter performance with large inflation by the RTPS. The best-performing localization scale of the LPF was identical to that of the LETKF. The LPF forecast was less accurate than the LETKF forecast. With sparse observations (RAOB), the LPF did not work.

3. The LPFGM showed more stability and lower forecast RMSEs than the LPF. In addition to the RTPS parameter, regulating the resampling frequency and the amplitude of Gaussian kernels was important for the LPFGM. The LPFGM without weight succession resulted in more stability and lower RMSEs than that with weight succession. With RAOB, the LPFGM forecast was more accurate than the LETKF forecast in sparsely observed regions, where the background ensemble spread and non-Gaussianity are larger.

As discussed in Sect. 4.3, there is much room for improvement in the LPF and LPFGM. While the LPFGM potentially provides a more accurate forecast than the LETKF, the LPFGM has more tunable parameters than the LETKF. Manually tuning these parameters by trying numerous experiments is computationally expensive. Therefore, adaptive methods for determining such tuneable parameters need to be explored. Also, it is important to investigate computationally more efficient methods to generate resampling matrices for the LPF and LPFGM.

The SPEEDY-based LETKF, LPF, and LPFGM used in this study are available as open-source software on GitHub (https://github.com/skotsuki/speedy-lpf, last access: 16 November 2022). This can act as a useful platform to investigate the LPF and LPFGM further in comparison with the well-known LETKF.

Appendix A: Derivations of Kalman gain and analysis error covariance

Here, we describe the derivations of Kalman gain and analysis error covariance (Eqs. 5 and 6). The Kalman gain ${\mathbf{K}}_t={\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T{\left({\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T+{\mathbf{R}}_t\right)}^{-\mathrm{1}}$ can be changed as follows:

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}{\mathbf{K}}_t& ={\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}\\ & \cdot \left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]{\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T{\left({\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T+{\mathbf{R}}_t\right)}^{-\mathrm{1}}\end{array}\text{(A2)}& {\displaystyle }\begin{array}{rl}& ={\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}\\ & \cdot \left[{\mathbf{H}}_t^T+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T\right]{\left({\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T+{\mathbf{R}}_t\right)}^{-\mathrm{1}}\end{array}\text{(A3)}& {\displaystyle }\begin{array}{rl}& ={\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}\left[{\mathbf{R}}_t+{\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T\right]\\ & \cdot {\left({\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}{\mathbf{H}}_t^T+{\mathbf{R}}_t\right)}^{-\mathrm{1}}\end{array}\text{(A4)}& {\displaystyle }={\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}.\end{array}$$

Using ${\mathbf{K}}_t={\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}$, analysis error covariance ${\mathbf{P}}_t^{\mathrm{a}}=\left(\mathbf{I}-{\mathbf{K}}_t{\mathbf{H}}_t\right){\mathbf{P}}_t^{\mathrm{b}}$ can be changed as follows:

$$\begin{array}{}\text{(A5)}& {\displaystyle }{\mathbf{P}}_t^{\mathrm{a}}=\left(\mathbf{I}-{\mathbf{K}}_t{\mathbf{H}}_t\right){\mathbf{P}}_t^{\mathrm{b}}={\mathbf{P}}_t^{\mathrm{b}}-{\mathbf{K}}_t{\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}\text{(A6)}& {\displaystyle }={\mathbf{P}}_t^{\mathrm{b}}-{\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t{\mathbf{P}}_t^{\mathrm{b}}\text{(A7)}& {\displaystyle }\begin{array}{rl}& ={\mathbf{P}}_t^{\mathrm{b}}-{\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}\\ & \cdot \left[\left({\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right)-{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}\right]{\mathbf{P}}_t^{\mathrm{b}}\end{array}\text{(A8)}& {\displaystyle }={\mathbf{P}}_t^{\mathrm{b}}-\left\{\mathbf{I}-{\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}\right\}{\mathbf{P}}_t^{\mathrm{b}}\text{(A9)}& {\displaystyle }={\left[{\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t\right]}^{-\mathrm{1}}.\end{array}$$

Consequently, we can derivate ${\left({\mathbf{P}}_t^{\mathrm{a}}\right)}^{-\mathrm{1}}={\left({\mathbf{P}}_t^{\mathrm{b}}\right)}^{-\mathrm{1}}+{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}{\mathbf{H}}_t$ from Eq. (A9) and ${\mathbf{K}}_t={\mathbf{P}}_t^{\mathrm{a}}{\mathbf{H}}_t^T{\mathbf{R}}_t^{-\mathrm{1}}$ from Eq. (A4).

Appendix B: The number of required samples for the MC resampling

This appendix investigates the number of required samples for the MC resampling to obtain accurate transform matrices stably. Here, we assume that 10 000 samples are sufficient to obtain an accurate transform matrix (${\mathbf{T}}_{N_{\mathrm{MC}}=\mathrm{10}\mathrm{000}}$). An absolute error (AE) of the transform matrix with G samples (${\mathbf{T}}_{N_{\mathrm{MC}}=G}$) is defined to measure the accuracy of the transform matrix:

$$\begin{array}{}\text{(B1)}& \text{AE}={\sum }_{i=\mathrm{1}}^m{\sum }_{j=\mathrm{1}}^m\mathrm{|}{\mathbf{T}}_{N_{\mathrm{MC}}=G}^{\left(i,j\right)}-{\mathbf{T}}_{N_{\mathrm{MC}}=\mathrm{10}\mathrm{000}}^{\left(i,j\right)}\mathrm{|}.\end{array}$$

Figure B1 shows the AE as a function of the number of samples for ensemble members 10, 20, 40, and 80. For this investigation, we obtained 1000 independent weights generated by uniform random numbers, and average, minimum, and maximum AEs of 1000 cases are shown. For four ensembles, AEs decrease rapidly for the first 1000 samples, followed by a gradual decrease until AE = 0. To reach 10 % of initial error at N_MC=1, MC resampling requires about 70, 300, 400, and 600 samples for ensemble members 10, 20, 40, and 80. It suggests that more ensemble requires more samples to obtain an accurate transform matrix stably by means of MC resampling.

Figure B1 Absolute errors of the transform matrix of MC resampling as a function of the number of samples (N_MC) verified against the transform matrix with 10 000 samples. Absolute errors are computed for 1000 independent cases that have different weights generated by uniform random numbers. Black bold lines show the average of 1000 cases, and gray shades represent minimum and maximum errors of the 1000 cases. Panels (a)–(d) show experiments with 10, 20, 40, and 80 ensembles, respectively. Dashed line shows the 10 % of initial error at N_MC=1.