State-of-the-art stochastic data assimilation methods for high-dimensional non-Gaussian problems

Abstract

This paper compares several commonly used state-of-the-art ensemble-based data assimilation methods in a coherent mathematical notation. The study encompasses different methods that are applicable to high-dimensional geophysical systems, like ocean and atmosphere and provide an uncertainty estimate. Most variants of Ensemble Kalman Filters, Particle Filters and second-order exact methods are discussed, including Gaussian Mixture Filters, while methods that require an adjoint model or a tangent linear formulation of the model are excluded. The detailed description of all the methods in a mathematically coherent way provides both novices and experienced researchers with a unique overview and new insight in the workings and relative advantages of each method, theoretically and algorithmically, even leading to new filters. Furthermore, the practical implementation details of all ensemble and particle filter methods are discussed to show similarities and differences in the filters aiding the users in what to use when. Finally, pseudo-codes are provided for all of the methods presented in this paper.

Keywords:

• ensemble Kalman filter
• particle filter
• data assimilation
• high dimension
• non Gaussian

1. Introduction

Data assimilation (DA) is the science of combining observations of a system, including their uncertainty, with estimates of that system from a dynamical model, including its uncertainty, to obtain a new and more accurate description of the system including an uncertainty estimate of that description. The uncertainty estimates point to an efficient description in terms of probability density functions, and in this paper we discuss methods that perform DA using an ensemble of model states to represent these probability density functions.

Ensemble Kalman filters are currently highly popular DA methods that are applied to a wide range of dynamical models including oceanic, atmospheric, and land surface models. The increasing popularity of Kalman-Filter-based ensemble (EnKF) methods in these fields is due to the relative ease of the filter implementation, increasing computational power and the natural forecast error evolution in EnKF schemes with the dynamical model in time. However, due to technological and scientific advances, three significant developments have occurred in the last decade that force us to look beyond standard Ensemble Kalman Filtering, which is based on linear and/or Gaussian assumptions. Firstly, continuous increase in computational capability has recently allowed to run operational models at high resolutions so that the dynamical models have become increasingly non-linear due to the direct resolution of small-scale non-linear processes in these models, e.g. small-scale turbulence. Secondly, in several geoscientific applications, such as atmosphere, ocean, land surface, hydrology and sea – i.e. it is of interest to estimate variables or parameters that are bounded requiring DA methods that can deal with non-Gaussian distributions. Thirdly, the observational network around the world has increased manyfold for weather, ocean and land surface areas providing more information about the real system with greater accuracy and higher spatial and temporal resolution. Often the so-called observation operators that connect model states to observations of these new observations are non-linear, again asking for non-Gaussian DA methods. Thus, the research in non-linear DAmethods, which can be applied to high-resolution dynamical models and/or complex observation operators, has seen major developments in the last decade with the aim to understand how existing ensemble methods cope with non-linearity in the models, to develop new ensemble methods that are more suited to non-linear dynamical models, as well as to explore non-linear filters that are not limited to Gaussian distributions, such as particle filters or hybrids between particle and ensemble Kalman filters.

The origin of this paper lies in the EU-funded research project SANGOMA (Stochastic Assimilation for the Next Generation Ocean Model Applications). The project focused on generating a coherent and transparent database of the current ensemble-based data assimilation methods and development of data assimilation tools suitable for non-linear and high-dimensional systems concentrating on methods that do not require tangent linear approximations of the model or its adjoint. The methods described within this paper have been applied in operational oceanography, like the TOPAZ system (Sakov et al., 2012) or the FOAM system of the UK Met Office (see Blockley et al., 2014). While TOPAZ is already using an EnKF, FOAM applies an optimal interpolation scheme that takes in less dynamical information to estimate the error covariance matrix. That said this paper is aimed at a very broad audience and data assimilation methods discussed in this paper are not limited to applications to ocean or atmosphere models, hence, the methods are presented without the context of any specific dynamical model, allowing the reader to make the most of each technique for their specific application.

A number of reviews have been published recently, each collating parts of the development in data assimilation, e.g. Bannister (2017) gives a comprehensive review of operational methods of variational and ensemble-variational data assimilation, Houtekamer and Zhang (2016) review the ensemble Kalman filter with a focus on application to atmospheric data assimilation, Law and Stuart (2012) review variational and Kalman filter methods, and Bocquet et al. (2010) present a review of concepts and ideas of non-Gaussian data assimilation methods and discusses various sources of non-Gaussianity. The merits of this paper lie within:

• coherent mathematical description of the main methods that are used in the current data assimilation community for application to high-dimensional and non-Gaussian problems allowing the reader to easily oversee the differences between the methods and compare them;

• discussing ensemble Kalman Filters, particle filters, second-order exact filters and Gaussian Mixture Filters within the same paper using consistent notation;

• inclusion of practical application aspects of these methods, discussing computational cost, parallelising, localisation and inflation techniques;

• provision of pseudo-code algorithms for all of the presented methods;

along with inclusion of recent developments, such as the error-subspace transform Kalman filter (ESTKF) and recent particle filters, this paper goes beyond earlier reviews (e.g. Tippett et al., 2003; Hamill et al., 2006; van Leeuwen, 2009; Houtekamer and Zhang, 2016; Bannister, 2017).

The paper is organised as follows: in Section 2 the common ground through Bayes theorem is established, in Section 3 a historical overview is given for both ensemble Kalman and particle filter fields, and in Section 4 we define the basic problem solved by all of the methods presented in this paper. In Section 5 we discuss the most popular types of ensemble Kalman filter methods. Then, in Section 6, we discuss several particle filter methods that can be applied to high-dimensional problems. In Section 7, we describe ensemble filters with second-order accuracy, namely the particle filter with Gaussian resampling (PFGR), the non-linear ensemble transform filter (NETF), and the moment-matching ensemble filter. The Gaussian mixture filter is discussed in Section 8. The practical implementation of the filters including localisation, inflation, parallelisation and the computation cost as well as the aspect of non-linearity are discussed in Section 9. Finally, Appendix 1 provides pseudo codes for resampling techniques often used in particle filter methods, and Appendix 2 contains pseudo codes for all of the methods discussed in this paper.

We note that many of the filters discussed in this paper are available freely from Sangoma project website¹ along with many other tools valuable and/or necessary for data assimilation systems.

1.1. Notation

In the data assimilation community the currently most accepted notation is described in Ide et al. (1997). We adhere to this notation where possible while also making this paper acceptable and intuitive not only to data-assimilation experts but also to a wider audience including those who might like to explore data assimilation methods simply as tools for their specific needs. To this end, throughout this paper dimensions will always be described by capital N with an underscore indicating the space in question, that is

• $N_x$ - dimension of state space;

• $N_y$ - dimension of observation space;

• $N_e$ - dimension of ensemble/particle space.

Further, the time index is always denoted in parentheses in the upper right corner of the variables, i.e. ${\left(.\right)}^{(m)}$, except for operators such as $\mathcal{M}$ (dynamical model) and $\mathcal{H}$ (observation operator) where it is in the lower right corner. However, we will omit the time index when possible to ease the notation. We willrefer to each ensemble member (or each particle) by ${\mathsf{x}}_j$ where the index $j=1,\dots ,N_e$ and $N_e$ is the total number of the ensemble members (or particles).

When discussing Bayesian theory in Sections 2 and 6 purely random variables will be denoted by capital letters, and fixed or deterministic or observed quantities will be denoted by lowercase letters. Probability density functions will be denoted by $p\left(..\right)$, and q(..) and we will use lower case arguments in this context.

Throughout the paper, Greek letters will refer to various errors, e.g. observational or model errors. Finally, bold lowercase letters will denote vectors and bold uppercase letters will denote matrices.

2. Common ground through Bayes theorem

Various types of data assimilation methods, e.g. variational, ensemble Kalman filters, particle filters, etc., have originated from different fields and backgrounds, due to the needs of a particular community or application. However, all of these methods can be unified through Bayes theorem. In this section, we will give a summary of Bayes theorem showing how both ensemble Kalman filter (KF) methods and particle filter (PF) methods are linked in this context and what problems each of them solve. For an introduction to the Bayesian theory for data assimilation, the reader is referred to e.g. van Leeuwen and Evensen (1996) and Wikle and Berliner (2006).

Data assimilation is an approach for combining observations with model forecasts to obtain a more accurate estimate of the state and its uncertainty. In this context, we require

• data, that is observations $\mathsf{y}$ and a knowledge of their associated error distributions and

• a prior, that is a model forecast of the state, ${\mathsf{x}}^f$, and knowledge of the associated forecast and model errors;

to obtain the posterior, i.e. the analysis state ${\mathsf{x}}^a$, and its associated error. The posterior can be computed through Bayes theorem which states that(1) $$\begin{array}{c}\hfill p\left(\mathsf{x}|\mathsf{y}\right)=\frac{p\left(\mathsf{y}|\mathsf{x}\right)p\left(\mathsf{x}\right)}{p\left(\mathsf{y}\right)},\end{array}$$(1)

where $p\left(\mathsf{x}|\mathsf{y}\right)$ is the posterior or analysis probability density function, $p\left(\mathsf{y}|\mathsf{x}\right)$ is the observation probability density function or also called the likelihood, $p\left(\mathsf{x}\right)$ is the prior or forecast probability density function, and $p\left(\mathsf{y}\right)$ is the marginal probability density function of the observations, which can be thought of as a normalising constant. From now on, for the ease of the readability, we will abbreviate ‘probability density function’ with ‘pdf’.

Typically, data assimilation methods make the Markovian assumption for the dynamical model $\mathcal{M}$ and the conditional independence assumption of the observations. That is, we assume that the model state or the prior at time m, when conditioned on all previous states only depends on the state at the time $m-1$,(2) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(0:T)}\right)=p\left({\mathsf{x}}^{(0)}\right)\prod _{m=1}^{T}p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}^{(m-1)}\right).\end{array}$$(2)

Here, the superscript 0 : T is to be read for time indices from initial time to time T, which is typically called assimilation window in data assimilation. Further, observations are also usually assumed to be conditionally independent, i.e. they are assumed to be independent in time,(3) $$\begin{array}{c}\hfill p\left({\mathsf{y}}^{(1:T)}|{\mathsf{x}}^{(0:T)}\right)=\prod _{m=1}^{T}p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}^{(m)}\right).\end{array}$$(3)

Using Equations (2(2) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(0:T)}\right)=p\left({\mathsf{x}}^{(0)}\right)\prod _{m=1}^{T}p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}^{(m-1)}\right).\end{array}$$(2) ) and (3(3) $$\begin{array}{c}\hfill p\left({\mathsf{y}}^{(1:T)}|{\mathsf{x}}^{(0:T)}\right)=\prod _{m=1}^{T}p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}^{(m)}\right).\end{array}$$(3) ) we can rewrite Bayes theorem in Equation (1(1) $$\begin{array}{c}\hfill p\left(\mathsf{x}|\mathsf{y}\right)=\frac{p\left(\mathsf{y}|\mathsf{x}\right)p\left(\mathsf{x}\right)}{p\left(\mathsf{y}\right)},\end{array}$$(1) ) as(4) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(T)}|{\mathsf{y}}^{(T)}\right)\propto p({\mathsf{x}}^0)\prod _{m=1}^{T}p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}^{(m)}\right)p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}^{(m-1)}\right).\end{array}$$(4)

The Markovian assumption allows us to use new observations as they become available by updating the previous estimate of the state process without having to start the calculations from scratch. This is called sequential updating and the methods described in this paper all follow this approach.

Ensemble Kalman filter methods solve this problem using Gaussian assumptions for both prior and likelihood pdf’s. Multiplying two Gaussian pdf’s leads again to a Gaussian pdf, i.e. the posterior or analysis pdf will also be Gaussian. The posterior pdf will have only one global maximum, which will correspond to the ensemble mean (also mode and median since the pdf is Gaussian). In other words, the posterior pdf in ensemble Kalman filter methods described in Section 5 is found in terms of the first two moments (mean and covariance) of the prior and likelihood pdf’s. This is also true when ensemble Kalman filters are applied to non-linear dynamical models or observation operators, in which case the information from higher moments in an ensemble KF analysis update is ignored. This is a shortcoming in ensemble Kalman filters when applied to non-Gaussian problems. However, in general ensemble Kalman filter methods are robust when applied to non-linear models and catastrophic filter divergence, where the filter deviates strongly from the observations while producing unrealistically small error estimates, occurs mainly due to sparse or inaccurate observations (Verlaan and Heemink, 2001; Tong et al., 2016). It should, of course, be realised that in non-linear settings the estimates of the posterior mean and covariance might be off.

In particle filter methods, the posterior is obtained using the prior and likelihood pdf’s directly in Equation (1(1) $$\begin{array}{c}\hfill p\left(\mathsf{x}|\mathsf{y}\right)=\frac{p\left(\mathsf{y}|\mathsf{x}\right)p\left(\mathsf{x}\right)}{p\left(\mathsf{y}\right)},\end{array}$$(1) ) without restricting them to being Gaussian. If both prior and likelihood are Gaussian the resulting posterior or analysis pdf is also Gaussian and has a global maximum corresponding to the mean state. However, if either or both prior and likelihood pdf’s are non-Gaussian then the resulting posterior pdf will also not be Gaussian. In other words, if the dynamical model or mapping of the model variables to observation space are non-linear then particle filter methods will produce an analysis pdf which will provide knowledge of more than the first two statistical moments (mean and covariance), in contrast to ensemble Kalman filter methods. Thus, the analysis pdf could be skewed, multi-modal or of varying width in comparison to a Gaussian pdf. Hence, particle filters are, by design, able to produce analysis pdf’s for non-Gaussian problems. While standard particle filter methods suffer from filter divergence for large problems recently several particle filter variants have been developed that avoid this divergence.

In what follows, we will describe numerous filtering methods in Sections 5, 6, 7, and 8 and discuss how each method attempts to produce an analysis pdf for non-Gaussian and high-dimensional problems. However, firstly we provide an overview of the historical development of both ensemble Kalman filters and particle filter methods to show how these fields have evolved and what has given rise in the development of each of the methods.

3. History of filtering for data assimilation

Before we precede to the main point of our paper – describing in unified notation current state-of-the-art ensemble and particle filter methods for non-linear and non-Gaussian applications, their implementation, and practical application, a short summary is in order on the historical development in both ensemble Kalman filter and particle filter areas.

3.1. Development history of ensemble Kalman filters

Ensemble data assimilation (EnDA) started in 1994 with the introduction of the Ensemble Kalman filter (EnKF, Evensen (1994)). The use of perturbed observations was introduced a few years later simultaneously by Burgers et al. (1998) and Houtekamer and Mitchell (1998) to correct the previously too low spread of the analysis ensemble. This filter formulation defines today the basic ’Ensemble Kalman filter’, which we will denote as the Stochastic Ensemble Kalman Filter, with a slightly different interpretation and implementation, as will be described later. The first alternative variant of the original EnKF was introduced by Pham et al. (1998a) in the form of Singular ’Evolutive’ Interpolated Kalman (SEIK) filter. The SEIK filter formulates the analysis step in the space spanned by the ensemble and hence is computationally particularly efficient. In contrast to the EnKF, which was formulated as a Monte Carlo method, the SEIK filter was designed to find the analysis ensemble by writing each posterior member as a linear combination of prior members without using perturbed observations. Another ensemble Kalman filter that uses the space spanned by the ensemble was introduced with the Error-Subspace Statistical Estimation (ESSE) method (Lermusiaux and Robinson, 1999).

The filters mentioned above were all introduced for data assimilation in oceanographic problems. A new set of filter methods was introduced during the years 2001 and 2002 for meteorological applications. The Ensemble Transform Kalman Filter (ETKF, Bishop et al., 2001) was first introduced in the context of meteorological adaptive sampling. Further, the Ensemble Adjustment Kalman Filter (EAKF, Anderson, 2001) and the Ensemble Square Root Filter (EnSRF, Whitaker and Hamill, 2002) were introduced. The motivation for these three filters was to avoid the use of perturbed observations, which were found to introduce additional sampling error into the filter solution, with the meteorological community apparently being unaware of the development of the SEIK filter. The new filters were classified as ensemble square root Kalman filters and presented in a uniform notation by Tippett et al. (2003). Nerger et al. (2005a) further classified the EnKF and SEIK filters as error-subspace Kalman filters because the filters compute the correction in the error-subspace spanned by the ensemble. This likewise holds for the ETKF, EAKF, and EnSRF, however, these filters do not explicitly use a basis in the error subspace but use the ensemble to represent the space. When the EAKF and EnSRF formulations are used to assimilate all observation at once, these filters exhibit a much larger computational cost compared to the ETKF. To reduce the cost, the original study on the EnSRF (Whitaker and Hamill, 2002) already introduced a variant in which observations are assimilated sequentially, which assumes that the observation errors are uncorrelated. A similar serial formulation of the EAKF was introduced by Anderson (2003). This sequential assimilation of observations was assessed by Nerger (2015) and it was shown that this formulation can destabilise the filtering process in cases when the observations have a strong influence.

With regard to the classification as an ensemble square root Kalman filter, the SEIK filter is the first filter method that was clearly formulated in square root form. The original EnKF uses the square root form only implicitly but an explicit square root formulation of the EnKF was presented by Evensen (2003).

The methods above all solve the original equations of the Kalman filter but use the sample covariance matrix of the ensemble to represent the state error covariance matrix. An alternative was introduced with the Maximum-Likelihood Kalman Filter (MLEF, Zupanski, 2005). This filter represents the first variant of the class of hybrid filters that were introduced in later years. The filter computes the maximum-a posteriori solution (in contrast to the minimum-variance solution of the Kalman filter) by an iterative scheme.²

While the EnKFs were very successful in making the application of the Kalman filter feasible for the high-dimensional problems in oceanography and meteorology, the affordable ensemble size was always very limited. To counter the issue of sampling error in ensemble covariances (the ensemble-sampled covariance has a rank of not more than the ensemble size minus one while the applications were of very high state dimension) the method of covariance localisation was introduced by Houtekamer and Mitchell (1998) and Houtekamer and Mitchell (2001). Later, an alternative localisation was introduced for the ETKF (LETKF, Hunt et al., 2007) which uses a local analysis (also used previously, e.g. by Cohn et al. (1998)) where observations are down-weighted with increasing distance from the local analysis point through a tapering of the inverse observation covariances.

The relationship between the SEIK filter and the ETKF was investigated by Nerger et al. (2012a). The study leads to a new filter formulation, the Error-Subspace Transform Kalman Filter (ESTKF), which combined the advantages of both filter formulations.

The filters mentioned above represent main developments of the ensemble Kalman filters. However, there are many other developments, which are not included here. Some of them are discussed in the sections below, in particular with regard to localisation. Overall, while there are different reviews of selections of ensemble Kalman filters, a complete and coherent overview of the different methods is still missing.

3.2. Development history of particle filters

Particle filters, like ensemble Kalman filters, are variants of Monte Carlo methods in which the probability distribution of the model state given the observations is approximated by a number of particles; however, unlike ensemble Kalman filters, particle filters are fully non-linear data assimilation techniques. From a sampling point of view, Ensemble Kalman Filters draw samples directly from the posterior since the probability distribution function (pdf) is assumed to be a Gaussian. In a particle filter application, the shape of the posterior is not known, and hence one cannot sample directly from it. In its simplest form, samples are generated from the prior after which importance sampling is employed to turn them into samples from the posterior where each sample is weighted with its likelihood value.

Particle filters emerged before ensemble Kalman filters, and when Gordon et al. (1993) introduced resampling in the sequential scheme the method became mainstream in non-linear filtering. This basic scheme has been made more efficient for specific applications in numerous ways, like looking ahead, adding small perturbations to resampled particles to avoid that they are the same etc. (see Doucet et al., 2001 for a very useful review of the many methods available at that time). Attempts to apply the particle filter to geophysical systems are as old as 1996 (van Leeuwen and Evensen, 1996), with the first partially successful application by van Leeuwen (2003a). However, until recently, particle filters have been deemed to be computationally unfeasible for large-dimensional systems due to the filter degeneracy problem (Bengtsson et al., 2008; Snyder et al., 2008; van Leeuwen, 2009). This means that the likelihood weights vary substantially between the particles when the number of independent observations is large, such that one particle obtains a weight close to one, while all the others have weight very close to zero. New developments in the field generated particle filter variants that have been shown to work for large dimensional systems with a limited number of particles. These methods can be divided in two classes: those that use localisation (starting with van Leeuwen, 2003b; Bengtsson et al., 2003), followed more recently by local variants of the ensemble transform particle filter (ETPF, Reich, 2013) and the Local Particle Filter (Poterjoy, 2016a) and those that exploit the future observational information via proposal densities, such as the Implicit Particle Filter (Chorin and Tu, 2009), the Equivalent Weights Particle Filter (EWPF, van Leeuwen, 2010; van Leeuwen, 2011; Ades and van Leeuwen, 2013), and the Implicit Equal Weights Particle Filter (IEWPF, Zhu et al., 2016).

In another development, second-order exact filters have been developed that ensure that the first two moments of the posterior pdf are consistent with the particle filter, and higher-order moments are not considered. The first paper of this kind was the Particle Filter with Gaussian Resampling of Xiong et al. (2006), followed by the Merging Particle Filter (Nakano et al., 2007) and the Moment Matching Ensemble Filter (Lei and Bickel, 2011). All these filters seem to have been developed independently. The Non-linear Ensemble Transform Filter (Tödter and Ahrens, 2015) can be considered a local version of the filter by Xiong et al. (2006), ironically again developed independently.

A further approximation to particle filtering is the Gaussian Mixture Filter first introduced in the geosciences by Bengtsson et al. (2003), followed by the adaptive Gaussian mixture filter variants (Hoteit et al., 2008; Stordal et al., 2011). The advantage of these filters over the standard particle filter is that each particle is ’dressed’ by a Gaussian such that the likelihood weights are calculated using a covariance that is broader than the pure observational covariance, leading to better behaving weights at the cost of reducing the influence of the observations on the posterior pdf (see e.g. van Leeuwen, 2009).

4. The problem

Consider the following non-linear stochastic discrete-time dynamical system at a time when observations are available:(5) $$\begin{array}{cc}\hfill {\mathsf{x}}^{(m)}& ={\mathcal{M}}_m\left({\mathsf{x}}^{(m-1)}\right)+{\beta }^{(m)}\hfill \end{array}$$(5) (6) $$\begin{array}{cc}\hfill {\mathsf{y}}^{(m)}& ={\mathcal{H}}_m\left({\mathsf{x}}^{(m)}\right)+{\beta }_o^{(m)},\hfill \end{array}$$(6)

where ${\mathsf{x}}^{(m)}\in {\mathcal{R}}^{N_x}$ is the $N_x$ dimensional state vector, ${\mathsf{y}}^{(m)}\in {\mathcal{R}}^{N_y}$ is the observation vector of size $N_y\ll N_x$, ${\mathcal{M}}_m:{\mathcal{R}}^{N_x}\to {\mathcal{R}}^{N_x}$ is the forward model operator, ${\mathcal{H}}_m:{\mathcal{R}}^{N_x}\to {\mathcal{R}}^{N_y}$ is the observation operator, ${\beta }^{(m)}\in {\mathcal{R}}^{N_x}$ is the model noise (or error) distributed Gaussian with a covariance matrix ${\mathsf{Q}}^{(m)}$, and ${\beta }_o^{(m)}\in {\mathcal{R}}^{N_y}$ is the observation noise (or error) distributed Gaussian with covariance matrix ${\mathsf{R}}^{(m)}$.

Then we can define an ensemble of model forecasts obtained using Equation (5(5) $$\begin{array}{cc}\hfill {\mathsf{x}}^{(m)}& ={\mathcal{M}}_m\left({\mathsf{x}}^{(m-1)}\right)+{\beta }^{(m)}\hfill \end{array}$$(5) ) for each ensemble or particle member as follows,(7) $$\begin{array}{c}\hfill {\mathsf{X}}^{f,(m)}=\left[{\mathsf{x}}_1^{f,(m)},{\mathsf{x}}_2^{f,(m)},\dots ,{\mathsf{x}}_{N_e}^{f,(m)}\right]\in {\mathcal{R}}^{N_x\times N_e},\end{array}$$(7)

where superscript ${\left(.\right)}^f$ stands for forecast.

The aim of the stochastic data assimilation methods is to produce a posterior pdf or analysis distribution of the state, ${\mathsf{X}}^a$, at the time of the observations through combining the ensemble model forecast ${\mathsf{X}}^f$ with observations $\mathsf{y}$. In Section 5, we will discuss ensemble Kalman filter based methods and in Sections 6–8 we will discuss particle, second-order exact, and adaptive Gaussian mixture filter methods all achieving this aim through different approaches.

5. Ensemble Kalman filters

Given an initial ensemble ${\mathsf{X}}^{(0)}\in {\mathcal{R}}^{N_x\times N_e}$, the different proposed variants of the ensemble Kalman filter have the following steps in common:

• Forecast step: the ensemble members at each time step between the observations $0<k\le m$ are propagated using the full non-linear dynamical model: (8) $$\begin{array}{c}\hfill {\mathsf{x}}_j^{f,(k)}={\mathcal{M}}_k\left({\mathsf{x}}_j^{f,(k-1)}\right)+{\beta }_j^{(k)},\end{array}$$(8) starting at the previous analysis ensemble (if $k=1$, then this would be ${\mathsf{x}}_j^{f,(1)}={\mathcal{M}}_1\left({\mathsf{x}}_j^{(0)}\right)+{\beta }_j^{(1)}$), where $j=1,\dots ,N_e$ is the ensemble member index.

• Analysis step: at the observation time $k=m$ the ensemble forecast mean and covariance are updated using the available observations to obtain a new analysis ensemble.

The various ensemble methods differ in the analysis step. Here we will discuss current methods applicable for large-dimensional systems, namely, the original ensemble Kalman filter (EnKF) (Evensen, 1994) with stochastic innovations (Burgers et al., 1998; Houtekamer and Mitchell, 1998), the singular evolutive interpolated Kalman filter (SEIK) (Pham et al., 1998a), the error-subspace statistical estimation (ESSE) (Lermusiaux and Robinson, 1999; Lermusiaux et al., 2002; Lermusiaux, 2007), the ensemble transform Kalman filter (ETKF) (Bishop et al., 2001), the ensemble adjustment Kalman filter (EAKF) (Anderson, 2001), the original ensemble square root filter (EnSRF) (Whitaker and Hamill, 2002) with synchronous and serial observation treatment, the square root formulation of the EnKF (Evensen, 2003), the error subspace transform Kalman filter (ESTKF) (Nerger et al., 2012a), and the maximum likelihood ensemble filter (MLEF) (Zupanski, 2005; Zupanski et al., 2008). We will present these methods in the square root form and point out the different ways the analysis ensemble is obtained in each of the methods. Tippett et al. (2003) gives a uniform framework for EnSRFs, which we follow closely here. In the rest of this section for ease of notation we omit the time index ${(·)}^{(k)}$ since all of the analysis operations are done at time m.

The ensemble methods discussed in this section are based on the Kalman filter (Kalman, 1960) where the updated ensemble mean follows the Kalman update for the state, given by(9) $$\begin{array}{c}\hfill {\overline{\mathsf{x}}}^a={\overline{\mathsf{x}}}^f+\mathsf{K}\left(\mathsf{y}-\mathcal{H}({\overline{\mathsf{x}}}^f)\right)={\overline{\mathsf{x}}}^f+\mathsf{K}\mathsf{d}\end{array}$$(9)

where $\mathsf{d}=\mathsf{y}-\mathcal{H}({\overline{\mathsf{x}}}^f)$ is the innovation. The ensemble covariance update follows the covariance update equation in the Kalman Filter, given by(10) $$\begin{array}{c}\hfill {\mathsf{P}}^a=(\mathsf{I}-\mathsf{K}\mathsf{H}){\mathsf{P}}^f,\end{array}$$(10)

where $\mathsf{K}$ is the Kalman gain given by(11) $$\begin{array}{c}\hfill \mathsf{K}={\mathsf{P}}^f{\mathsf{H}}^{\mathsf{T}}{(\mathsf{H}{\mathsf{P}}^f{\mathsf{H}}^{\mathsf{T}}+\mathsf{R})}^{-1}.\end{array}$$(11)

The matrix $\mathsf{H}$ is the linearised observation operator $\mathcal{H}(..)$ at the forecast mean ${\overline{\mathsf{x}}}^f$. Initially the Kalman filter was derived for a linear observation operator, but in the Extended Kalman Filter the non-linear observation operator is used as above.

Since for high-dimensional systems it is computationally not feasible to form the error covariance matrix $\mathsf{P}$, the analysis update of the covariance matrix in Equation (10(10) $$\begin{array}{c}\hfill {\mathsf{P}}^a=(\mathsf{I}-\mathsf{K}\mathsf{H}){\mathsf{P}}^f,\end{array}$$(10) ) is formulated in a square root form by computing a transform matrix and applying it to the ensemble perturbation matrix, which is a scaled square root of $\mathsf{P}$. That is, the analysis ensemble is then given by(12) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^a+{{\mathsf{X}}^{\prime }}^a,\end{array}$$(12)

where ${\overline{\mathsf{X}}}^a=({\overline{\mathsf{x}}}^a,\dots ,{\overline{\mathsf{x}}}^a)\in {\mathcal{R}}^{N_x\times N_e}$ is a matrix with the ensemble analysis mean in each column and the ensemble analysis perturbations are a scaled matrix square root of(13) $$\begin{array}{c}\hfill {\mathsf{P}}^a=\frac{{{\mathsf{X}}^{\prime }}^a{\left({{\mathsf{X}}^{\prime }}^a\right)}^{\mathsf{T}}}{N_e-1}.\end{array}$$(13)

To obtain the general square root form we write, using (10(10) $$\begin{array}{c}\hfill {\mathsf{P}}^a=(\mathsf{I}-\mathsf{K}\mathsf{H}){\mathsf{P}}^f,\end{array}$$(10) )(14) $$\begin{array}{cc}\hfill {{\mathsf{X}}^{\prime }}^a{({{\mathsf{X}}^{\prime }}^a)}^{\mathsf{T}}& =\left(\mathsf{I}-{\mathsf{P}}^f{\mathsf{H}}^{\mathsf{T}}{(\mathsf{H}{\mathsf{P}}^f{\mathsf{H}}^{\mathsf{T}}+\mathsf{R})}^{-1}\mathsf{H}\right){{\mathsf{X}}^{\prime }}^f{({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}\hfill \\ \hfill & ={{\mathsf{X}}^{\prime }}^f\left(\mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}^{-1}\mathsf{S}\right){({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}},\hfill \end{array}$$(14)

where $\mathsf{S}=\mathsf{H}{{\mathsf{X}}^{\prime }}^f$ is the ensemble perturbation matrix in observation space and(15) $$\begin{array}{c}\hfill \mathsf{F}=\mathsf{S}{\mathsf{S}}^{\mathsf{T}}+(N_e-1)\mathsf{R}\end{array}$$(15)

is the innovation covariance.

It is possible to use a slightly different way to calculate matrix $\mathsf{S}$ using the non-linear observation operator as $\mathsf{S}=\mathcal{H}({\mathsf{X}}^f)-\mathcal{H}({\overline{\mathsf{X}}}^f)$, in which, with a slight abuse of notation, $\mathcal{H}({\mathsf{X}}^f)=\left(\mathcal{H}({\mathsf{x}}_1^f),\cdots ,\mathcal{H}({\mathsf{x}}_{N_e}^f)\right)$, and similarly for $\mathcal{H}({\overline{\mathsf{X}}}^f)$. This can be used in any of the ensemble Kalman filters discussed below.

To find the updated ensemble analysis perturbations ${{\mathsf{X}}^{\prime }}^a$ we need to compute the square root $\mathsf{T}$ of the matrix(16) $$\begin{array}{c}\hfill \mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}^{-1}\mathsf{S}=\mathsf{T}{\mathsf{T}}^{\mathsf{T}},\end{array}$$(16)

where $\mathsf{T}$ is called a transform matrix. Different ways exist to compute the transform matrix $\mathsf{T}$ and here we will discuss the current methods applicable to large-dimensional systems.

For the ensemble-based Kalman filters presented in this paper we can write the analysis update as linear transformations using a weight vector $\overline{\mathsf{w}}$ for the ensemble mean and a weight matrix ${\mathsf{W}}^{\prime }$ for the ensemble perturbations as(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) (18) $$\begin{array}{cc}\hfill {{\mathsf{X}}^{\prime }}^a& ={{\mathsf{X}}^{\prime }}^f{\mathsf{W}}^{\prime }.\hfill \end{array}$$(18)

Notice, that the ensemble analysis perturbation matrix, ${{\mathsf{X}}^{\prime }}^a$, in Equation (18(18) $$\begin{array}{cc}\hfill {{\mathsf{X}}^{\prime }}^a& ={{\mathsf{X}}^{\prime }}^f{\mathsf{W}}^{\prime }.\hfill \end{array}$$(18) ) has a zero mean by construction. Further, we note that for most of the methods discussed in this section, matrix ${\mathsf{W}}^{\prime }$ is the transformation matrix $\mathsf{T}$ in Equation (16(16) $$\begin{array}{c}\hfill \mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}^{-1}\mathsf{S}=\mathsf{T}{\mathsf{T}}^{\mathsf{T}},\end{array}$$(16) ). However, this is not the case for EnKF, SEnKF and MLEF. Further, we can compute the analysis ensemble directly by(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19)

where $\overline{\mathsf{W}}=\left(\overline{\mathsf{w}},\dots ,\overline{\mathsf{w}}\right)$. In the sections below we will derive the weight matrices for each of the ensemble-based Kalman filter methods we discuss. The updated ensemble can then be obtained using Equation (19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ).

To aid simplicity in discussing the different methods we use the same letter for the variables with the same meaning, i.e. ${\mathsf{W}}^{\prime }$ is always the perturbation analysis transform matrix that transforms ${{\mathsf{X}}^{\prime }}^f$ into ${{\mathsf{X}}^{\prime }}^a$. Clearly, such variables do not necessarily have the same values for the various methods listed below. Thus, we subscript these variables common to all methods with a specific letter for each method. This letter is underlined in the title of each subsection that follows here, e.g. for EnKF we use ${{\mathsf{W}}^{\prime }}_N$. Note that some of the variables can have the same values for different methods, though. At the end of this section we will provide a table of the common variables with their dimensions and whether they are equal to the same variable in a different method.

5.1. The Stochastic Ensemble Kalman filter (E$\underline{n}$KF)

The Stochastic EnKF was introduced at the same time by Burgers et al. (1998) and Houtekamer and Mitchell (1998). It is a modified version of the original under-dispersive EnKF as introduced by Evensen (1994) by adding measurement noise to the innovations so that the filter maintains the correct spread in the analysis ensemble and prevents filter divergence.

Although the scheme was initially interpreted as perturbing observations, a more consistent interpretation is that the predicted observations are perturbed with the observation noise. The reason for this is that it doesn’t make sense to perturb observations since they already contain measurement noise (errors), e.g. from measuring instruments, and thus have already departed from the true state of the system. Also, Bayes Theorem, see Section 2 tells us that we need the probabilities of the states given this set of observations, not a perturbed set. The idea is that each ensemble member is statistically equivalent to the true state of the system, and the true observation is a perturbed measurement of the true state. So to compare that observation with the predicted observations the latter have to be perturbed with the measurement noise too to make this comparison meaningful. This reasoning is identical to that used in rank histograms in which observations are ranked in the perturbed predicted observations from the ensemble to be statistically equivalent.

Each ensemble member individually is explicitly corrected using the Kalman filter equations, and hence the square root form is implicit only as the transform matrix and its square root are never explicitly computed. In contrast to the other filters, the stochastic EnKF perturbs the predicted observations by forming a matrix(20) $$\begin{array}{c}\hfill {\mathsf{Y}}^f=\left(\mathcal{H}({\mathsf{x}}_1^f),\mathcal{H}({\mathsf{x}}_2^f),\dots ,\mathcal{H}({\mathsf{x}}_{N_e}^f)\right)+{\mathsf{Y}}^{\prime }\in {\mathcal{R}}^{N_y\times N_e},\end{array}$$(20)

where the observational noise (perturbation) matrix ${\mathsf{Y}}^{\prime }$ is given by:(21) $$\begin{array}{c}\hfill {\mathsf{Y}}^{\prime }=\left({ϵ}_1,{ϵ}_2,\dots ,{ϵ}_{N_e}\right)\in {\mathcal{R}}^{N_y\times N_e}\end{array}$$(21)

with the noise vectors ${ϵ}_j$ drawn from a Gaussian distribution with mean zero and covariance $\mathsf{R}$. We also introduce the observation matrix $\mathsf{Y}=\left(\mathsf{y},\mathsf{y},\dots ,\mathsf{y}\right)\in {\mathcal{R}}^{N_y\times N_e}$ consisting of $N_e$ identical copies of the observation vector.

The Stochastic EnKF uses the matrix $\mathsf{F}$ defined in Equation (15(15) $$\begin{array}{c}\hfill \mathsf{F}=\mathsf{S}{\mathsf{S}}^{\mathsf{T}}+(N_e-1)\mathsf{R}\end{array}$$(15) ) with prescribed matrix $\mathsf{R}$ and proceeds by transforming all ensemble members according to(22) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\mathsf{X}}^f+\frac{1}{N_e-1}{{\mathsf{X}}^{\prime }}^f{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}_N^{-1}\left(\mathsf{Y}-{\mathsf{Y}}^f\right),\end{array}$$(22)

Similar to the Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) )–(19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ), this can be written as(23) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\mathsf{X}}^f+{{\mathsf{X}}^{\prime }}^f{{\mathsf{W}}^{\prime }}_N\end{array}$$(23)

with(24) $$\begin{array}{c}\hfill {{\mathsf{W}}^{\prime }}_N=\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathbf{F}}_N^{-1}\left(\mathsf{Y}-{\mathsf{Y}}^f\right).\end{array}$$(24)

Due to the use of the observation ensemble $\mathsf{Y}$ no explicit transformation of the ensemble mean needs to be performed.

Algorithm  in Appendix 2 gives a pseudo-algorithm of the EnKF method.

We note that while the above description of the stochastic EnKF is widely accepted and implemented, it does produce the correct posterior covariance only in a statistical sense due to extra sampling errors while the ensemble mean is not affected by the sampling error by ensuring that observation noise matrix, ${\mathsf{Y}}^{\prime }$, has zero mean. However, in the limit of infinite ensemble size and when all sources of error (both observation and model) are correctly sampled, the stochastic EnKF does produce the correct posterior covariance (Whitaker and Hamill, 2002).

5.2. The singular evolutive interpolated Kalman filter (S$\underline{E}$IK)

The SEIK filter (Pham et al., 1998b Pham, 2001) was the first filter method that allowed for non-linear model evolution and that was explicitly formulated in square root form. The filter uses the Sherman–Morrison–Woodbury identity (Golub and Van Loan, 1996) to rewrite $\mathsf{T}{\mathsf{T}}^{\mathsf{T}}$ (Equation 16(16) $$\begin{array}{c}\hfill \mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}^{-1}\mathsf{S}=\mathsf{T}{\mathsf{T}}^{\mathsf{T}},\end{array}$$(16) ) as(25) $$\begin{array}{c}\hfill \mathsf{T}{\mathsf{T}}^{\mathsf{T}}=\mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}^{-1}\mathsf{S}={\left(\mathsf{I}+\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}.\end{array}$$(25)

Note, that the performance of this scheme depends on whether the product of the inverse of the observation error matrix, ${\mathsf{R}}^{-1}$, and a given vector can be efficiently computed, which is for instance the case when we assume that the observation errors are uncorrelated.

The SEIK filter computes the analysis step in the ensemble error subspace. This is achieved by defining a matrix(26) $$\begin{array}{c}\hfill {\mathsf{L}}_E={\mathsf{X}}^f{\mathsf{A}}_E,\end{array}$$(26)

where ${\mathsf{A}}_E\in {\mathcal{R}}^{N_e\times (N_e-1)}$ is a matrix with full rank and zero column sums. Commonly, matrix ${\mathsf{A}}_E$ is identified as(27) $$\begin{array}{c}\hfill {\mathsf{A}}_E=\left[\begin{array}{c}{\mathsf{I}}_{N_e-1\times N_e-1}\\ {\mathsf{0}}_{1\times N_e-1}\end{array}\right]-\frac{1}{N_e}\left[{\mathsf{1}}_{N_e\times N_e-1}\right],\end{array}$$(27)

where $\mathsf{0}$ is a matrix whose elements are equal to zero and $\mathsf{1}$ is a matrix whose elements are equal to one (Pham et al., 1998b). Matrix ${\mathsf{A}}_E$ implicitly subtracts the ensemble mean when the matrix $\mathsf{L}$ is computed. In addition, ${\mathsf{A}}_E$ removes the last column of ${{\mathsf{X}}^{\prime }}^f$. Thus, $\mathsf{L}$ is an $N_e\times N_e-1$ matrix that holds the first $N_e-1$ ensemble perturbations. The product of the square root matrices in the ensemble error space becomes now(28) $$\begin{array}{c}\hfill {\mathsf{T}}_E{\mathsf{T}}_E^{\mathsf{T}}={\left({\mathsf{A}}_E^{\mathsf{T}}{\mathsf{A}}_E+\frac{1}{N_e-1}{(\mathsf{H}{\mathsf{L}}_E)}^{\mathsf{T}}{\mathsf{R}}^{-1}(\mathsf{H}{\mathsf{L}}_E)\right)}^{-1}.\end{array}$$(28)

The matrix ${\mathsf{T}}_E{\mathsf{T}}_E^{\mathsf{T}}$ is of size $N_e-1\times N_e-1$. The square root ${\mathsf{T}}_E$ is obtained from the Cholesky decomposition of ${({\mathsf{T}}_E{\mathsf{T}}_E^{\mathsf{T}})}^{-1}$. Then, the ensemble transformation weight matrices in Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) )–(19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ) are given by(29) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_E& ={\mathsf{A}}_E{\mathsf{T}}_E\mathrm{\Omega },\hfill \end{array}$$(29) (30) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_E& =\frac{1}{\sqrt{N_e-1}}{\mathsf{A}}_E{\mathsf{T}}_E{\mathsf{T}}_E^{\mathsf{T}}{(\mathsf{H}{\mathsf{L}}_E)}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{d}.\hfill \end{array}$$(30)

Here, the columns of $\mathrm{\Omega }\in {\mathcal{R}}^{N_e-1\times N_e}$ are orthonormal and orthogonal to the vector ${(1,\dots ,1)}^{\mathsf{T}}$. $\mathrm{\Omega }$ can be either random or a deterministic rotation matrix. However, if a deterministic $\mathrm{\Omega }$ is used then Nerger et al. (2012a) shows that a symmetric square root of ${\mathsf{T}}_E{\mathsf{T}}_E^{\mathsf{T}}$ should be used for a more stable ensemble.

Algorithm  in Appendix 2 gives a pseudo-algorithm of the SEIK method.

5.3. The error-subspace statistical estimation (E$\underline{S}$SE)

The ESSE (Lermusiaux and Robinson, 1999 Lermusiaux et al., 2002 Lermusiaux, 2007) method is based on evolving an error subspace of variable size, that spans and tracks the scales and processes where the dominant errors occur (Lermusiaux et al., 2002). Here, we follow the formulation of Lermusiaux (2007) adapted to the unified notation used here.

The consideration of an evolving error subspace is analogous to the motivation of the SEIK filter. The main difference to other subspace filters mentioned here is how the ensemble matrix is truncated. That is, the full ensemble perturbation matrix ${{\mathsf{X}}^{\prime }}^f$ at the current analysis time with columns$$\begin{array}{c}\hfill {{\mathsf{x}}^{\prime }}_j^f={\mathsf{x}}_j^f-{\overline{\mathsf{x}}}^f,j=1,\dots ,N_e\end{array}$$

is approximated by the fastest growing singular vectors. The full ensemble perturbation matrix is decomposed using the reduced or thin singular value decomposition (SVD), (e.g. p. 72, Golub and Van Loan, 1996),(31) $$\begin{array}{c}\hfill {\mathsf{U}}_S{\mathrm{\Sigma }}_S{\mathsf{V}}_S^{\mathsf{T}}={{\mathsf{X}}^{\prime }}^f,\end{array}$$(31)

where ${\mathsf{U}}_S\in {\mathcal{R}}^{N_x\times N_e}$ is an orthogonal matrix of left singular vectors, ${\mathrm{\Sigma }}_S\in {\mathcal{R}}^{N_e\times N_e}$ is a diagonal matrix with singular values on the diagonal, and ${\mathsf{V}}_S^{\mathsf{T}}\in {\mathcal{R}}^{N_e\times N_e}$ is an orthogonal matrix of right singular vectors of ${{\mathsf{X}}^{\prime }}^f$. Next, normalised eigenvalues are computed via(32) $$\begin{array}{c}\hfill {\mathsf{E}}_S=\frac{1}{N_e-1}{\mathrm{\Sigma }}_S^2.\end{array}$$(32)

The matrices ${\mathsf{U}}_S$ and ${\mathsf{E}}_S$ are truncated to the leading eigenvalues. Using ${\tilde{\mathsf{U}}}_S$, ${\tilde{\mathsf{E}}}_S$ with rank ${\tilde{N}}_e\le N_e$ and ${\widehat{\mathsf{U}}}_S$, ${\widehat{\mathsf{E}}}_S$ with rank $\widehat{p}<{\tilde{N}}_e$ where the similarity coefficient $\rho$ is computed via(33) $$\begin{array}{c}\hfill \rho =\frac{\text{Tr}\left({\widehat{\mathsf{E}}}_S^{\frac{1}{2}}{\widehat{\mathsf{U}}}_S^{\mathsf{T}}{\tilde{\mathsf{U}}}_S{\tilde{\mathsf{E}}}_S^{\frac{1}{2}}\right)}{\text{Tr}\left({\tilde{\mathsf{E}}}_S\right)}\end{array}$$(33)

and $\text{Tr}(.)$ is the trace of a matrix. $\rho$ measures the similarity between two subspaces of different sizes. The process of reducing the subspace is repeated until $\rho$ is close to one, i.e. $\rho >\alpha$ where $1-ϵ\le \alpha \le 1$ is a user selected scalar limit.³ The dimension of the error subspace thus varies with time and in accord with model dynamics (Lermusiaux, 2007). Hence, in the following analysis update the reduced rank approximations(34) $$\begin{array}{ccc}\hfill {\tilde{\mathsf{U}}}_S& \approx \hfill & \hfill {\mathsf{U}}_S,\end{array}$$(34) (35) $$\begin{array}{ccc}\hfill {\tilde{\mathsf{E}}}_S& \approx \hfill & \hfill {\mathsf{E}}_S,\end{array}$$(35) (36) $$\begin{array}{ccc}\hfill {\tilde{\mathsf{V}}}_S& \approx \hfill & \hfill {\mathsf{V}}_S\end{array}$$(36)

are used where the right singular vector matrix ${\tilde{\mathsf{V}}}_S$ is also truncated to have size ${\tilde{N}}_e\times {\tilde{N}}_e$.

The product of the square root matrices, using Equation (14(14) $$\begin{array}{cc}\hfill {{\mathsf{X}}^{\prime }}^a{({{\mathsf{X}}^{\prime }}^a)}^{\mathsf{T}}& =\left(\mathsf{I}-{\mathsf{P}}^f{\mathsf{H}}^{\mathsf{T}}{(\mathsf{H}{\mathsf{P}}^f{\mathsf{H}}^{\mathsf{T}}+\mathsf{R})}^{-1}\mathsf{H}\right){{\mathsf{X}}^{\prime }}^f{({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}\hfill \\ \hfill & ={{\mathsf{X}}^{\prime }}^f\left(\mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}^{-1}\mathsf{S}\right){({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}},\hfill \end{array}$$(14) ), in the error subspace becomes(37) $$\begin{array}{c}\hfill {\mathsf{T}}_S{\mathsf{T}}_S^{\mathsf{T}}=\mathsf{I}-{\tilde{\mathsf{S}}}^{\mathsf{T}}{\tilde{\mathsf{F}}}^{-1}\tilde{\mathsf{S}}\end{array}$$(37)

where ensemble errors in observation space are given by $\tilde{\mathsf{S}}=\mathsf{H}{\tilde{\mathsf{U}}}_S{\tilde{\mathsf{E}}}_S^{\frac{1}{2}}{\tilde{\mathsf{V}}}_S^T$ and innovation covariances by $\tilde{\mathsf{F}}=\mathsf{H}{\tilde{\mathsf{U}}}_S\tilde{\mathsf{E}}{\tilde{\mathsf{U}}}^{\mathsf{T}}{\mathsf{H}}^{\mathsf{T}}+\mathsf{R}$.

The inverse of the $N_y\times N_y$-matrix $\tilde{\mathsf{F}}$ is obtained by performing the eigenvalue decomposition (EVD)(38) $$\begin{array}{c}\hfill \tilde{\mathsf{F}}=\mathrm{\Gamma }\mathrm{\Lambda }{\mathrm{\Gamma }}^{\mathsf{T}}\end{array}$$(38)

so that Equation (37(37) $$\begin{array}{c}\hfill {\mathsf{T}}_S{\mathsf{T}}_S^{\mathsf{T}}=\mathsf{I}-{\tilde{\mathsf{S}}}^{\mathsf{T}}{\tilde{\mathsf{F}}}^{-1}\tilde{\mathsf{S}}\end{array}$$(37) ) becomes(39) $$\begin{array}{c}\hfill {\mathsf{T}}_S{\mathsf{T}}_S^{\mathsf{T}}=\mathsf{I}-{\tilde{\mathsf{S}}}^{\mathsf{T}}\mathrm{\Gamma }{\mathrm{\Lambda }}^{-1}{\mathrm{\Gamma }}^{\mathsf{T}}\tilde{\mathsf{S}}.\end{array}$$(39)

Performing another EVD in Equation (39(39) $$\begin{array}{c}\hfill {\mathsf{T}}_S{\mathsf{T}}_S^{\mathsf{T}}=\mathsf{I}-{\tilde{\mathsf{S}}}^{\mathsf{T}}\mathrm{\Gamma }{\mathrm{\Lambda }}^{-1}{\mathrm{\Gamma }}^{\mathsf{T}}\tilde{\mathsf{S}}.\end{array}$$(39) ),(40) $$\begin{array}{c}\hfill {\mathsf{T}}_S{\mathsf{T}}_S^{\mathsf{T}}=\mathsf{Z}\mathrm{\Pi }{\mathsf{Z}}^{\mathsf{T}},\end{array}$$(40)

the symmetric square root becomes(41) $$\begin{array}{c}\hfill {\mathsf{T}}_S=\mathsf{Z}{\mathrm{\Pi }}^{\frac{1}{2}}{\mathsf{Z}}^T.\end{array}$$(41)

Hence, the ensemble transformation weight matrices needed to form the ensemble analysis mean and analysis perturbations in Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) )–(19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ) are given by(42) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_S^a& =\mathsf{Z}{\mathrm{\Pi }}^{\frac{1}{2}}{\mathsf{Z}}^{\mathsf{T}}\hfill \end{array}$$(42) (43) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_S^a& =\frac{1}{\sqrt{N_e-1}}\tilde{\mathsf{S}}{\tilde{\mathsf{F}}}^{-1}\mathsf{d}.\hfill \end{array}$$(43)

Note, that when computing the analysis ensemble mean and perturbations, the truncated ensemble perturbation matrix ${\widetilde{{\mathsf{X}}^{\prime }}}^f$ is used in the pseudo-algorithm  in Appendix 2. The truncation to the rank ${\tilde{N}}_e$ will results in a reduction of the ensemble size. To avoid that the ensemble size shrinks, Lermusiaux (2007) described an optional adaptive method to generate new ensemble members.

5.4. The ensemble transform Kalman filter (E$\underline{T}$KF)

The ETKF (Bishop et al., 2001) was derived to explicitly transform the ensemble in a way that results in the correct spread of the analysis ensemble. As the SEIK filter, the ETKF uses the Morrison-Woodbury identity to write(44) $$\begin{array}{c}\hfill {\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}={\left(\mathsf{I}+\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}.\end{array}$$(44)

In contrast to the SEIK filter, ${\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}$ is of size $N_e\times N_e$ and hence represents the error-subspace of dimension $N_e-1$ indirectly by the full ensemble.

Currently, the most widespread method to compute the update in the ETKF appears to be the formulation of the LETKF by Hunt et al. (2007), which we describe here. By performing the EVD of the symmetric matrix ${({\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}})}^{-1}={\mathsf{U}}_T{\mathrm{\Sigma }}_T{\mathsf{U}}_T^{\mathsf{T}}$ we obtain the symmetric square root(45) $$\begin{array}{c}\hfill {\mathsf{T}}_T={\mathsf{U}}_T{\mathrm{\Sigma }}_T^{-\frac{1}{2}}{\mathsf{U}}_T^{\mathsf{T}}.\end{array}$$(45)

Using this decomposition, the ensemble transformation weight matrices needed to form the ensemble analysis mean and analysis perturbations in Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) )–(19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ) are given by(46) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_T& ={\mathsf{U}}_T{\mathrm{\Sigma }}_T^{-\frac{1}{2}}{\mathsf{U}}_T^{\mathsf{T}},\hfill \end{array}$$(46) (47) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_T& =\frac{1}{\sqrt{(N_e-1)}}{\mathsf{U}}_T{\mathrm{\Sigma }}_T^{-1}{\mathsf{U}}_T^{\mathsf{T}}{({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}{\mathsf{H}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{d}.\hfill \end{array}$$(47)

Using the symmetric square root produces a transform matrix which is closest to the identity matrix in the Frobenius norm (Hunt et al., 2007). Thus, the ETKF results in a minimum transform in the ensemble space, which is different from the notion of ’optimal transportation’ used in the ETPF (see Section 6.3).

The original publication introducing the ETKF (Bishop et al., 2001) did not specify the form of the matrix square root ${\mathsf{T}}_T$. There are different possibilities to compute it, and taking a simple single-sided square root could lead to implementations with a biased transformation, such that the transformation by ${\mathsf{W}}^{\prime }$ would not preserve the ensemble mean. However, using the symmetric square root approach this bias is avoided. Livings (2005) proposed another variant normalising first the forecast observation ensemble perturbation matrix so that the observations are dimensionless with standard deviation one(48) $$\begin{array}{c}\hfill \tilde{\mathsf{S}}=\frac{1}{\sqrt{N_e-1}}{\mathsf{R}}^{-\frac{1}{2}}\mathsf{S}.\end{array}$$(48)

Substituting (48(48) $$\begin{array}{c}\hfill \tilde{\mathsf{S}}=\frac{1}{\sqrt{N_e-1}}{\mathsf{R}}^{-\frac{1}{2}}\mathsf{S}.\end{array}$$(48) ) into (44(44) $$\begin{array}{c}\hfill {\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}={\left(\mathsf{I}+\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}.\end{array}$$(44) ) gives(49) $$\begin{array}{c}\hfill {\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}={\left(\mathsf{I}+{\tilde{\mathsf{S}}}^{\mathsf{T}}\tilde{\mathsf{S}}\right)}^{-1}.\end{array}$$(49)

To find the square root form next we perform the SVD(50) $$\begin{array}{c}\hfill {\tilde{\mathsf{S}}}^{\mathsf{T}}={\mathsf{U}}_T{\tilde{\mathrm{\Sigma }}}_T{\tilde{\mathsf{V}}}_T^{\mathsf{T}}.\end{array}$$(50)

In this case, the ensemble transformation weight matrices in Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) )–(19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ) become(51) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_T& ={\mathsf{U}}_T{\left(\mathsf{I}+{\tilde{\mathrm{\Sigma }}}_T{\tilde{\mathrm{\Sigma }}}_T^{\mathsf{T}}\right)}^{-\frac{1}{2}}{\mathsf{U}}_T^{\mathsf{T}},\hfill \end{array}$$(51) (52) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_T& =\frac{1}{\sqrt{N_e-1}}{\mathsf{U}}_T{(\mathsf{I}+{\tilde{\mathrm{\Sigma }}}^{\mathsf{T}}{\tilde{\mathrm{\Sigma }}}_T^{\mathsf{T}})}^{-1}{\tilde{\mathrm{\Sigma }}}_T{\tilde{\mathsf{V}}}_T^{\mathsf{T}}{\mathsf{R}}^{-\frac{1}{2}}\mathsf{d}.s\hfill \end{array}$$(52)

This formulation avoids the multiplication ${\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}$ and can hence prevent possible loss of accuracy due to rounding errors. However, this formulation also requires the computation of the square root of $\mathsf{R}$, which itself can result in rounding errors if $\mathsf{R}$ is not diagonal.

Algorithm  in Appendix 2 gives a pseudo-algorithm of the ETKF method.

5.5. The ensemble adjustment Kalman filter (E$\underline{A}$KF)

The EAKF was introduced by Anderson (2001). Similarly to the SEIK filter and the ETKF we require here that the matrix ${\mathsf{R}}^{-1}$ is readily available. Using scaled ensemble perturbations as discussed for the ETKF-formulation by Livings (2005) in Equations (48(48) $$\begin{array}{c}\hfill \tilde{\mathsf{S}}=\frac{1}{\sqrt{N_e-1}}{\mathsf{R}}^{-\frac{1}{2}}\mathsf{S}.\end{array}$$(48) )–(49(49) $$\begin{array}{c}\hfill {\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}={\left(\mathsf{I}+{\tilde{\mathsf{S}}}^{\mathsf{T}}\tilde{\mathsf{S}}\right)}^{-1}.\end{array}$$(49) ) we can write(53) $$\begin{array}{c}\hfill {\mathsf{T}}_A{\mathsf{T}}_A^{\mathsf{T}}={\left(\mathsf{I}+{\tilde{\mathsf{S}}}^{\mathsf{T}}\tilde{\mathsf{S}}\right)}^{-1}.\end{array}$$(53)

We perform the SVD on the scaled forecast ensemble observation perturbation matrix(54) $$\begin{array}{c}\hfill {\tilde{\mathsf{S}}}^{\mathsf{T}}={\mathsf{U}}_A{\mathrm{\Sigma }}_A{\mathsf{V}}_A^{\mathsf{T}}.\end{array}$$(54)

Note that ${\mathsf{U}}_A={\mathsf{U}}_T$, related to the similarity between the EAKF and the ETKF. We also use an EVD to obtain(55) $$\begin{array}{c}\hfill {\mathsf{P}}^f={\mathsf{Z}}_A{\mathrm{\Gamma }}_A{\mathsf{Z}}_A^{\mathsf{T}}.\end{array}$$(55)

The decomposition in Equation (55(55) $$\begin{array}{c}\hfill {\mathsf{P}}^f={\mathsf{Z}}_A{\mathrm{\Gamma }}_A{\mathsf{Z}}_A^{\mathsf{T}}.\end{array}$$(55) ) is usually performed as an SVD of the ensemble perturbation matrix ${{\mathsf{X}}^{\prime }}^f$, which approximates ${\mathsf{P}}^f$ using $N_e$ ensemble members.

Due to the ranks of the matrices decomposed in Equations (54(54) $$\begin{array}{c}\hfill {\tilde{\mathsf{S}}}^{\mathsf{T}}={\mathsf{U}}_A{\mathrm{\Sigma }}_A{\mathsf{V}}_A^{\mathsf{T}}.\end{array}$$(54) ) and (55(55) $$\begin{array}{c}\hfill {\mathsf{P}}^f={\mathsf{Z}}_A{\mathrm{\Gamma }}_A{\mathsf{Z}}_A^{\mathsf{T}}.\end{array}$$(55) ) there are at most $q=min(N_e-1,N_y)$ non-zero singular values in ${\mathrm{\Sigma }}_A$ and at most $N_e-1$ non-zero eigenvalues in ${\mathrm{\Gamma }}_A$. Thus, the matrices in the equations below can be truncated as follows: ${\mathsf{U}}_A\in {\mathcal{R}}^{N_e\times q}$, ${\mathrm{\Sigma }}_A\in {\mathcal{R}}^{q\times q}$, ${\mathsf{V}}_A^{\mathsf{T}}\in {\mathcal{R}}^{q\times N_y}$ and ${\mathrm{\Gamma }}_A\text{,}{\mathsf{Z}}_A\in {\mathcal{R}}^{N_e-1\times N_e-1}$. Then, the ensemble transformation weight matrices in Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) )–(19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ) are given by(56) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_A& ={\mathsf{U}}_A{\left(\mathsf{I}+{\mathrm{\Sigma }}_A{\mathrm{\Sigma }}_A^{\mathsf{T}}\right)}^{-\frac{1}{2}}{\mathrm{\Gamma }}_A^{-\frac{1}{2}}{\mathsf{Z}}_A^{\mathsf{T}}{{\mathsf{X}}^{\prime }}^f,\hfill \end{array}$$(56) (57) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_A& =\frac{1}{\sqrt{N_e-1}}{\mathsf{U}}_A{(\mathsf{I}+{\mathrm{\Sigma }}_A^{\mathsf{T}}{\mathrm{\Sigma }}_A)}^{-1}{\mathrm{\Sigma }}_T{\mathsf{V}}_A^{\mathsf{T}}{\mathsf{R}}^{-\frac{1}{2}}\mathsf{d}.\hfill \end{array}$$(57)

Note, that the EAKF perturbation weight matrix in Equation (56(56) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_A& ={\mathsf{U}}_A{\left(\mathsf{I}+{\mathrm{\Sigma }}_A{\mathrm{\Sigma }}_A^{\mathsf{T}}\right)}^{-\frac{1}{2}}{\mathrm{\Gamma }}_A^{-\frac{1}{2}}{\mathsf{Z}}_A^{\mathsf{T}}{{\mathsf{X}}^{\prime }}^f,\hfill \end{array}$$(56) 57(57) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_A& =\frac{1}{\sqrt{N_e-1}}{\mathsf{U}}_A{(\mathsf{I}+{\mathrm{\Sigma }}_A^{\mathsf{T}}{\mathrm{\Sigma }}_A)}^{-1}{\mathrm{\Sigma }}_T{\mathsf{V}}_A^{\mathsf{T}}{\mathsf{R}}^{-\frac{1}{2}}\mathsf{d}.\hfill \end{array}$$(57) ) is the same as applying the orthogonal matrix ${\mathrm{\Gamma }}_A^{-\frac{1}{2}}{\mathsf{Z}}_A^{\mathsf{T}}{{\mathsf{X}}^{\prime }}^f$ instead of the orthogonal matrix ${\mathsf{U}}_T$ in the ETKF perturbation transform matrix given by Equation (51(51) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_T& ={\mathsf{U}}_T{\left(\mathsf{I}+{\tilde{\mathrm{\Sigma }}}_T{\tilde{\mathrm{\Sigma }}}_T^{\mathsf{T}}\right)}^{-\frac{1}{2}}{\mathsf{U}}_T^{\mathsf{T}},\hfill \end{array}$$(51) 52(52) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_T& =\frac{1}{\sqrt{N_e-1}}{\mathsf{U}}_T{(\mathsf{I}+{\tilde{\mathrm{\Sigma }}}^{\mathsf{T}}{\tilde{\mathrm{\Sigma }}}_T^{\mathsf{T}})}^{-1}{\tilde{\mathrm{\Sigma }}}_T{\tilde{\mathsf{V}}}_T^{\mathsf{T}}{\mathsf{R}}^{-\frac{1}{2}}\mathsf{d}.s\hfill \end{array}$$(52) ) (Tippett et al., 2003).

The decomposition in Equation (55(55) $$\begin{array}{c}\hfill {\mathsf{P}}^f={\mathsf{Z}}_A{\mathrm{\Gamma }}_A{\mathsf{Z}}_A^{\mathsf{T}}.\end{array}$$(55) ) is costly due to the size of the matrix to be decomposed. For this reason, the EAKF is typically applied with serial observation processing as will be described for the EnSRF in Section 5.7.

Algorithm  in Appendix 2 gives a pseudo-algorithm of the EAKF method.

5.6. The ensemble square root filter (EnS$\underline{R}$F)

The EnSRF was introduced by Whitaker and Hamill (2002) to avoid the use of perturbed observations by a square root formulation. In the EnSRF the transform matrix is given by$${\mathsf{T}}_R{\mathsf{T}}_R^{\mathsf{T}}=\mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}^{-1}\mathsf{S}.$$

We first perform an EVD of $\mathsf{F}$ to obtain its inverse(58) $$\begin{array}{c}\hfill {\mathsf{F}}^{-1}={\mathrm{\Gamma }}_R{\mathrm{\Lambda }}_R^{-1}{\mathrm{\Gamma }}_R^{\mathsf{T}}.\end{array}$$(58)

Then, we can write the ensemble analysis covariance as(59) $$\begin{array}{cc}\hfill {{\mathsf{X}}^{\prime }}^a{({{\mathsf{X}}^{\prime }}^a)}^{\mathsf{T}}& ={{\mathsf{X}}^{\prime }}^f\left(\mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathrm{\Gamma }}_R{\mathrm{\Lambda }}_R^{-1}{\mathrm{\Gamma }}_R^{\mathsf{T}}\mathsf{S}\right){({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}\hfill \\ \hfill & ={{\mathsf{X}}^{\prime }}^f\left(\mathsf{I}-{\mathsf{G}}_R{\mathsf{G}}_R^{\mathsf{T}}\right){({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}\hfill \end{array}$$(59)

where ${\mathsf{G}}_R={\mathsf{S}}^{\mathsf{T}}{\mathrm{\Gamma }}_R{\mathrm{\Lambda }}_R^{-\frac{1}{2}}$. Decomposing ${\mathsf{G}}_R={\mathsf{U}}_R{\mathrm{\Sigma }}_R{\mathsf{V}}_R^{\mathsf{T}}$ using an SVD we obtain$$\begin{array}{cc}\hfill {{\mathsf{X}}^{\prime }}^a{({{\mathsf{X}}^{\prime }}^a)}^{\mathsf{T}}& ={{\mathsf{X}}^{\prime }}^f\left(\mathsf{I}-\left[{\mathsf{U}}_R{\mathrm{\Sigma }}_R{\mathsf{V}}_R^{\mathsf{T}}\right]{\left[{\mathsf{U}}_R{\mathrm{\Sigma }}_R{\mathsf{V}}_R^{\mathsf{T}}\right]}^{\mathsf{T}}\right){({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}\hfill \\ \hfill & ={{\mathsf{X}}^{\prime }}^f{\mathsf{U}}_R\left(\mathsf{I}-{\mathrm{\Sigma }}_R{\mathrm{\Sigma }}_R^{\mathsf{T}}\right){\mathsf{U}}_R^{\mathsf{T}}{({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}.\hfill \end{array}$$

The diagonal matrix holding the singular values is of dimension ${\mathrm{\Sigma }}_R\in {\mathcal{R}}^{N_e\times N_y}$ and has thus at most $min(N_e,N_y)$ nonzero singular values. To reduce the computational cost for the case of high dimensional models with $N_e\ll N_y$, we can truncate to get the much smaller matrix ${\mathrm{\Sigma }}_R\in {\mathcal{R}}^{N_e\times min(N_e,N_y)}$ (see Table 1). The square root form for the ensemble analysis perturbations is given by(60) $$\begin{array}{c}\hfill {{\mathsf{X}}^{\prime }}^a={{\mathsf{X}}^{\prime }}^f{\mathsf{U}}_R{\left(\mathsf{I}-{\mathrm{\Sigma }}_R{\mathrm{\Sigma }}_R^{\mathsf{T}}\right)}^{\frac{1}{2}},\end{array}$$(60)

and the ensemble transformation weight matrices needed to form the ensemble analysis mean and analysis perturbations in Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) )–(19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ) are given by(61) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_R& ={\mathsf{U}}_R{\left(\mathsf{I}-{\mathrm{\Sigma }}_R{\mathrm{\Sigma }}_R^{\mathsf{T}}\right)}^{\frac{1}{2}}{\mathsf{U}}_R^{\mathsf{T}},\hfill \end{array}$$(61) (62) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_R& ={\mathsf{S}}^{\mathsf{T}}{\mathrm{\Gamma }}_R{\mathrm{\Lambda }}_R^{-1}{\mathrm{\Gamma }}_R^{\mathsf{T}}\mathsf{d},\hfill \end{array}$$(62)

where in Equation (61(61) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_R& ={\mathsf{U}}_R{\left(\mathsf{I}-{\mathrm{\Sigma }}_R{\mathrm{\Sigma }}_R^{\mathsf{T}}\right)}^{\frac{1}{2}}{\mathsf{U}}_R^{\mathsf{T}},\hfill \end{array}$$(61) 62(62) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_R& ={\mathsf{S}}^{\mathsf{T}}{\mathrm{\Gamma }}_R{\mathrm{\Lambda }}_R^{-1}{\mathrm{\Gamma }}_R^{\mathsf{T}}\mathsf{d},\hfill \end{array}$$(62) ) we have post-multiplied the ensemble analysis perturbations by the orthogonal matrix of the left singular vectors ${\mathsf{U}}_R^{\mathsf{T}}$ to ensure that the analysis ensemble is unbiased (Livings et al., 2008; Sakov and Oke, 2008).

Algorithm  in Appendix 2 gives a pseudo-algorithm of the EnSRF method.

5.7. EnSRF with serial observation treatment

The serial observation treatment in the EnSRF was introduced by Whitaker and Hamill (2002) together with the EnSRF assimilating all observations at once. The serial treatment reduces the computing cost. Hence, the EnSRF is typically not applied with the bulk update described above, but with serial treatment of observations, which is possible if $\mathsf{R}$ is diagonal. In this case, each single observation can be assimilated separately. Thus, $\mathsf{F}$ reduces to the scalar F and $\mathsf{S}{\mathsf{S}}^{\mathsf{T}}$ to the scalar $S^2$. For a single observation ($N_y=1$), the matrix ${\mathsf{G}}_R$ becomes a vector given by(63) $$\begin{array}{c}\hfill {\mathsf{G}}_R=\frac{1}{\sqrt{F}}{\mathsf{S}}^{\mathsf{T}}.\end{array}$$(63)

All singular values of ${\mathsf{G}}_R$ are zero except the first, which is its norm,(64) $$\begin{array}{c}\hfill {\mathrm{\Sigma }}_R=\frac{S}{\sqrt{F}}\mathsf{e}\end{array}$$(64)

where $\mathsf{e}$ is a vector with $N_e$ zero elements except the first, which is one. The first column of ${\mathsf{U}}_R$ corresponds to the normalised vector ${\mathsf{S}}^{\mathsf{T}}$(65) $$\begin{array}{c}\hfill {\mathsf{U}}_R\mathsf{e}=\frac{1}{S}{\mathsf{S}}^{\mathsf{T}}.\end{array}$$(65)

The square root of the diagonal matrix in Equation (61(61) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_R& ={\mathsf{U}}_R{\left(\mathsf{I}-{\mathrm{\Sigma }}_R{\mathrm{\Sigma }}_R^{\mathsf{T}}\right)}^{\frac{1}{2}}{\mathsf{U}}_R^{\mathsf{T}},\hfill \end{array}$$(61) ) can be written as a sum of the identity matrix and a matrix proportional to $\mathsf{e}{\mathsf{e}}^{\mathsf{T}}$:(66) $$\begin{array}{c}\hfill {\left(\mathsf{I}-{\mathrm{\Sigma }}_R{\mathrm{\Sigma }}_R^{\mathsf{T}}\right)}^{\frac{1}{2}}=\mathsf{I}-(1-\sqrt{(N_e-1)R/F})\mathsf{e}{\mathsf{e}}^{\mathsf{T}}.\end{array}$$(66)

Using Equation (65(65) $$\begin{array}{c}\hfill {\mathsf{U}}_R\mathsf{e}=\frac{1}{S}{\mathsf{S}}^{\mathsf{T}}.\end{array}$$(65) ) and the fact that all columns of $\mathsf{U}$ are orthonormal, one obtains(67) $$\begin{array}{c}\hfill {{\mathsf{W}}^{\prime }}_R=\mathsf{I}-\frac{1-\sqrt{(N_e-1)R/F}}{S^2}{\mathsf{S}}^{\mathsf{T}}\mathsf{S}\end{array}$$(67)

and the weight vector for the update of the ensemble mean is(68) $$\begin{array}{c}\hfill {\overline{\mathsf{w}}}_R=\frac{1}{F}{\mathsf{S}}^{\mathsf{T}}\mathsf{d}.\end{array}$$(68)

The equations above are then applied in a series over each single observation. The equations are likewise valid when the EAKF is formulated with a serial observation treatment.

Algorithm  in Appendix 2 gives a pseudo-algorithm of the EnSRF method with serial observation treatment.

5.8. The square root formulation of the stochastic ensemble Kalman filter (SEnK$\underline{F}$)

The SEnKF was introduced by Evensen (2003) as a square root formulation of the stochastic EnKF. Defining ${\mathsf{Y}}^{\prime }$ as for the stochastic EnKF and using a matrix(69) $$\begin{array}{c}\hfill {\mathsf{F}}_F=\mathsf{S}{\mathsf{S}}^{\mathsf{T}}+{\mathsf{Y}}^{\prime }{{\mathsf{Y}}^{\prime }}^{\mathsf{T}}\end{array}$$(69)

we obtain the matrix(70) $$\begin{array}{c}\hfill {\mathsf{T}}_F{\mathsf{T}}_F^{\mathsf{T}}=\mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}_F^{-1}\mathsf{S}.\end{array}$$(70)

We could decompose ${\mathsf{F}}_F$ using an EVD but this is costly if $N_y\gg N_e$ (Evensen, 2003). Instead, we assume that forecast and observation errors are uncorrelated, i.e.(71) $$\begin{array}{c}\hfill \mathsf{S}{{\mathsf{Y}}^{\prime }}^{\mathsf{T}}\equiv 0,\end{array}$$(71)

so that(72) $$\begin{array}{c}\hfill {\mathsf{F}}_F=\mathsf{S}{\mathsf{S}}^{\mathsf{T}}+{\mathsf{Y}}^{\prime }{{\mathsf{Y}}^{\prime }}^{\mathsf{T}}=(\mathsf{S}+{\mathsf{Y}}^{\prime }){(\mathsf{S}+{\mathsf{Y}}^{\prime })}^{\mathsf{T}}.\end{array}$$(72)

Now we can use an SVD to decompose $\mathsf{S}+{\mathsf{Y}}^{\prime }={\mathsf{U}}_F{\mathrm{\Sigma }}_F{\mathsf{V}}_F^{\mathsf{T}}$, giving(73) $$\begin{array}{c}\hfill {\mathsf{F}}_F={\mathsf{U}}_F{\mathrm{\Sigma }}_F{\mathrm{\Sigma }}_F^{\mathsf{T}}{\mathsf{U}}_F^{\mathsf{T}},\end{array}$$(73)

which has a much smaller computational cost than decomposing ${\mathsf{F}}_F$ using an EVD when $N_y\gg N_e$.

The ensemble transformation is then computed according to Equation (23(23) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\mathsf{X}}^f+{{\mathsf{X}}^{\prime }}^f{{\mathsf{W}}^{\prime }}_N\end{array}$$(23) ) with the weight matrix given by(74) $$\begin{array}{c}\hfill {{\mathsf{W}}^{\prime }}_F={\mathsf{S}}^{\mathsf{T}}{\mathsf{U}}_F{\mathrm{\Sigma }}_F^{-1}{({\mathrm{\Sigma }}_F^{-1})}^{\mathsf{T}}{\mathsf{U}}_F^{\mathsf{T}}(\mathsf{Y}-{\mathsf{Y}}^f).\end{array}$$(74)

Algorithm  in Appendix 2 gives a pseudo-algorithm of the EnKF in square root form.

5.9. The error-subspace transform Kalman filter (EST$\underline{K}$F)

The ESTKF has been derived from the SEIK filter (Nerger et al., 2012a) by combining the advantages of the SEIK filter and the ETKF. The ESTKF exhibits better properties than the SEIK filter, like a minimum ensemble transformation as the ETKF. However, unlike the ETKF, the ESTKF computes the ensemble transformation in the error subspace spanned by the ensemble rather than using the ensemble representation of it. That is, the error subspace of the dimension $N_e-1$ is represented directly in the ESTKF (similarly to the SEIK filter) while in the ETKF the error subspace is represented indirectly using the full ensemble of size $N_e$.

Similar to the SEIK filter, a projection matrix ${\mathsf{A}}_K\in {\mathcal{R}}^{N_e\times N_e-1}$ is used whose elements are defined by(75) $$\begin{array}{c}\hfill {{\mathsf{A}}_K}_{\{i,j\}}:=\left\{\begin{array}{cc}\hfill 1-\frac{1}{N_e}\frac{1}{\frac{1}{\sqrt{N_e}}+1}& \mathrm{f}ori=j,i<N_e\hfill \\ \hfill -\frac{1}{N_e}\frac{1}{\frac{1}{\sqrt{N_e}}+1}& \mathrm{f}ori\ne j,i<N_e\hfill \\ \hfill -\frac{1}{\sqrt{N_e}}& \mathrm{f}ori=N_e\hfill \end{array}\right.\end{array}$$(75)

With this projection, the basis vectors for the error subspace are given by(76) $$\begin{array}{c}\hfill {\mathsf{L}}_K={\mathsf{X}}^f{\mathsf{A}}_K.\end{array}$$(76)

As for the matrix $\mathrm{\Omega }$ in the SEIK filter, the columns of matrix ${\mathsf{A}}_K$ are orthonormal and orthogonal to the vector ${(1,\cdots ,1)}^{\mathsf{T}}$. When the matrix ${\mathsf{L}}_K$ is computed, the multiplication with ${\mathsf{A}}_K$ implicitly subtracts the ensemble mean. Further, ${\mathsf{A}}_K$ subtracts a fraction of the last column of ${{\mathsf{X}}^{\prime }}^f$ from all other columns. In this way, the last column of ${{\mathsf{X}}^{\prime }}^f$ is not just dropped as in the SEIK filter, but its information is distributed over the other columns. The product of the square root matrices in the error subspace becomes now(77) $$\begin{array}{c}\hfill {\mathsf{T}}_K{\mathsf{T}}_K^{\mathsf{T}}={\left(\mathsf{I}+\frac{1}{N_e-1}{(\mathsf{H}{\mathsf{L}}_K)}^{\mathsf{T}}{\mathsf{R}}^{-1}(\mathsf{H}{\mathsf{L}}_K)\right)}^{-1}.\end{array}$$(77)

By performing the EVD of the symmetric matrix ${({\mathsf{T}}_K{\mathsf{T}}_K^{\mathsf{T}})}^{-1}={\mathsf{U}}_K{\mathrm{\Sigma }}_K{\mathsf{U}}_K^{\mathsf{T}}$ we obtain the symmetric square root(78) $$\begin{array}{c}\hfill {\mathsf{T}}_K={\mathsf{U}}_K{\mathrm{\Sigma }}_K^{-\frac{1}{2}}{\mathsf{U}}_K^{\mathsf{T}}.\end{array}$$(78)

Then, the ensemble transformation weight matrices needed to form the ensemble analysis mean and perturbations in Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) )–(19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ) are given by(79) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_K& ={\mathsf{A}}_K{\mathsf{T}}_K{\mathsf{A}}_K^{\mathsf{T}},\hfill \end{array}$$(79) (80) $$\begin{array}{cc}\hfill {\overline{\mathsf{w}}}_K& =\frac{1}{\sqrt{N_e-1}}{\mathsf{A}}_K{\mathsf{U}}_K{\mathrm{\Sigma }}_K^{-1}{\mathsf{U}}_K^{\mathsf{T}}{(\mathsf{H}{\mathsf{L}}_K)}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{d}.\hfill \end{array}$$(80)

Compared to the SEIK filter, both the matrices ${\mathsf{A}}_E$ and $\mathrm{\Omega }$ are replaced by ${\mathsf{A}}_K$ in the ESTKF. In addition, the ESTKF uses the symmetric square root of ${\mathsf{T}}_K{\mathsf{T}}_K^{\mathsf{T}}$. The use of ${\mathsf{A}}_K$ leads to a consistent projection onto the error subspace and back onto the state space, while the symmetric square root ensures that the minimum transformation is obtained. It is also possible to apply the ESTKF with a random ensemble transformation. For this case, the rightmost matrix ${\mathsf{A}}_K$ in Equation (79(79) $$\begin{array}{cc}\hfill {{\mathsf{W}}^{\prime }}_K& ={\mathsf{A}}_K{\mathsf{T}}_K{\mathsf{A}}_K^{\mathsf{T}},\hfill \end{array}$$(79) ) is replaced by a random matrix with the same properties as the deterministic ${\mathsf{A}}_K$.

Algorithm  in Appendix 2 gives a pseudo-algorithm of the ESTKF method.

5.10. The Maximum Likelihood Ensemble Filter ($\underline{M}$LEF)

The MLEF (Zupanski, 2005; Zupanski et al., 2008) calculates the state estimate as the maximum of the posterior probability density (pdf) function. This is in contrast to the ensemble Kalman filter methods described in this paper, which are based on the minimum variance approach, so targeting the mean. The maximum of the pdf is found by an iterative minimisation of the cost function using a generalised non-linear conjugate-gradient method.

The original MLEF filter (Zupanski, 2005) uses a second-order Taylor approximation to the analysis increments, which requires that the cost function is twice differentiable. However, this requirement is not necessarily satisfied in many real life non-linear applications, for example where the parameterisation of some processes is used in models or for strongly non-linear observation operators. Here, we present the revised MLEF by Zupanski et al. (2008) that avoids this requirement by using a non-differentiable minimisation algorithm.

In contrast to all other ensemble filters discussed above, the MLEF maintains the actual state estimate separately from the ensemble, which is used to provide the measurement of estimation error. Thus, the analysis perturbations in the MLEF are computed for each ensemble member in a square root form without re-centring them onto the analysis ensemble mean. Hence, this filter does not follow the same square root form as the filters described above and is presented last in this section.

In the MLEF, the ensemble analysis perturbations are defined using the difference between analysis and forecast for each ensemble member, ${{\mathsf{x}}^{\prime }}_j^a={\mathsf{x}}_j^a-{\mathsf{x}}_j^f$ and not between ensemble analysis states and the analysis mean. They are found using generalised Hessian preconditioning in state space. A change of variable is performed as follows(81) $$\begin{array}{cc}\hfill {\mathsf{x}}_j^a& ={\mathsf{x}}_j^f+{{\mathsf{x}}^{\prime }}_j^a,\hfill \end{array}$$(81) (82) $$\begin{array}{cc}\hfill {{\mathsf{x}}^{\prime }}_j^a& ={\mathsf{G}}^{\frac{1}{2}}{\xi }_j,\hfill \end{array}$$(82)

where the matrix ${\mathsf{G}}^{\frac{1}{2}}={{\mathsf{X}}^{\prime }}^f{\left(\mathsf{I}+\mathsf{C}\right)}^{-\frac{1}{2}}\in {\mathcal{R}}^{N_x\times N_e}$ represents the inverse square root of the generalised Hessian estimated at the initial point of minimisation, and ${\xi }_j$ is a control variable defined in ensemble subspace. Matrix $\mathsf{C}$ is the covariance matrix(83) $$\begin{array}{c}\hfill \mathsf{C}={({{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}{\mathsf{H}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{H}{{\mathsf{X}}^{\prime }}^f.\end{array}$$(83)

Equation (82(82) $$\begin{array}{cc}\hfill {{\mathsf{x}}^{\prime }}_j^a& ={\mathsf{G}}^{\frac{1}{2}}{\xi }_j,\hfill \end{array}$$(82) ) can be written as a transformation of ensemble perturbations by(84) $$\begin{array}{c}\hfill {{\mathsf{x}}^{\prime }}_j^a={{\mathsf{X}}^{\prime }}^f{\mathsf{w}}_{M,j}\end{array}$$(84)

where the elements of the weight vector ${\mathsf{w}}_{M,j}\in {\mathcal{R}}^{N_e}$ for ensemble member j are given by(85) $$\begin{array}{c}\hfill {\mathsf{w}}_{M,j}={\left(\mathsf{I}+\mathsf{C}\right)}^{-\frac{1}{2}}{\xi }_j.\end{array}$$(85)

Now we use an EVD of $\mathsf{C}={\mathrm{\Gamma }}_M{\mathrm{\Lambda }}_M{\mathrm{\Gamma }}_M^{\mathsf{T}}$ to write Equation (85(85) $$\begin{array}{c}\hfill {\mathsf{w}}_{M,j}={\left(\mathsf{I}+\mathsf{C}\right)}^{-\frac{1}{2}}{\xi }_j.\end{array}$$(85) ) as(86) $$\begin{array}{c}\hfill {\mathsf{w}}_{M,j}={\mathrm{\Gamma }}_M{\left(\mathsf{I}+{\mathrm{\Lambda }}_M\right)}^{-\frac{1}{2}}{\mathrm{\Gamma }}_M^{\mathsf{T}}{\xi }_j.\end{array}$$(86)

We note, that in a linear case matrix ${\mathsf{G}}^{\frac{1}{2}}$ is a square root of ${\mathsf{P}}^a$. Indeed the same decomposition and inversion was used to find the square root analysis perturbations for the ETKF, see Equation (45(45) $$\begin{array}{c}\hfill {\mathsf{T}}_T={\mathsf{U}}_T{\mathrm{\Sigma }}_T^{-\frac{1}{2}}{\mathsf{U}}_T^{\mathsf{T}}.\end{array}$$(45) ).

After successfully accomplishing the Hessian preconditioning, the next step in the iterative minimisation is to calculate the gradient in the ensemble-spanned subspace. The preconditioned generalised gradient at the k-th minimisation iteration is obtained by(87) $$\begin{array}{c}\hfill {\mathrm{\nabla }}_{G}J({\mathsf{x}}_k)={\left[\mathsf{I}+\mathsf{C}\right]}^{-1}{\xi }_k-\mathsf{Z}{({\mathsf{x}}_k)}^{\mathsf{T}}{\mathsf{R}}^{-\frac{1}{2}}\left[\mathsf{y}-\mathcal{H}({\mathsf{x}}_k)\right],\end{array}$$(87)

where(88) $$\begin{array}{cc}\hfill \mathsf{Z}(\mathsf{x})& =[{\mathsf{z}}_1(\mathsf{x}),{\mathsf{z}}_2(\mathsf{x}),\dots ,{\mathsf{z}}_{N_e}(\mathsf{x})]\hfill \end{array}$$(88) (89) $$\begin{array}{cc}\hfill {\mathsf{z}}_j(\mathsf{x})& ={\mathsf{R}}^{-\frac{1}{2}}[\mathcal{H}(\mathsf{x}+{{\mathsf{x}}^{\prime }}_j^f)-\mathcal{H}(\mathsf{x})].\hfill \end{array}$$(89)

Upon convergence we have obtained an optimal state analysis$$\begin{array}{c}\hfill {\mathsf{x}}^a={\mathsf{x}}_k.\end{array}$$

To complete the non-differential formulation of the MLEF, ensemble analysis perturbations are computed as follows(90) $$\begin{array}{c}\hfill {{\mathsf{X}}^{\prime }}^a={{\mathsf{X}}^{\prime }}^f{\left[\mathsf{I}+{\left(\mathsf{Z}({\mathsf{x}}^a)\right)}^{\mathsf{T}}\mathsf{Z}({\mathsf{x}}^a)\right]}^{-\frac{1}{2}}\end{array}$$(90)

where $\mathsf{Z}({\mathsf{x}}^a)$ is obtained by substituting $\mathsf{x}={\mathsf{x}}^a$ into Equation (89(89) $$\begin{array}{cc}\hfill {\mathsf{z}}_j(\mathsf{x})& ={\mathsf{R}}^{-\frac{1}{2}}[\mathcal{H}(\mathsf{x}+{{\mathsf{x}}^{\prime }}_j^f)-\mathcal{H}(\mathsf{x})].\hfill \end{array}$$(89) ).

Algorithm  in Appendix 2 gives a pseudo-algorithm of the MLEF method.

5.11. Summary of ensemble Kalman Filter methods

In this section, we have described ten most popular ensemble Kalman filter methods that are applicable to high-dimensional non-Gaussian problems.

This collection of methods could be categorised in different ways, for example in deterministic ensemble filters, where the analysis is found through explicit mathematical transformations (SEIK, ETKF, EAKF, EnSRF, ESTKF, MLEF), and stochastic ensemble filters, where perturbed forecasted observations are used (EnKF, SEnKF). Burgers et al. (1998) and Houtekamer and Mitchell (1998) showed that in order to maintain sufficient spread in the ensemble and prevent filter divergence, the observations should be treated as random variables, i.e. perturbed, while our interpretation is slightly different, as described above. This stochasticity, of course, leads to extra sampling noise in the filters. On the other hand, Lawson and Hansen (2004) showed that for large ensemble sizes, stochastic filters can handle non-linearity better than the deterministic filters. This is due to the additional Gaussian observation spread normalising the ensemble update in the stochastic filter, which tends to erase the non-Gaussian higher moments non-linear error growth has generated. However, current computational power restricts us to small ensemble sizes for high-dimensional problems, in which case stochastic filters add another source of sampling error thus underestimating the analysis update (Whitaker and Hamill, 2002).

While all ensemble Kalman filter methods use low-rank approximations of the state error covariance matrix, some of the methods in this section are referred to as error-subspace ensemble filters because they directly operate in the error subspace spanned by ensemble rather than using the ensemble representation of it. Such filters are SEIK (see Section 5.2), ESTKF (see Section 5.9), and ESSE (see Section 5.3). Nerger et al. (2005b) compares the stochastic EnKF with the SEIK filter in an idealised high-dimensional shallow water model with non-linear evolution, showing that the main difference between the filters lies in the efficiency of the representation of the covariance matrix $\mathsf{P}$. In general, the EnKF filter will require a larger ensemble size $N_e$ to achieve the same performance as the SEIK filter. The relation of the ETKF and SEIK methods has been studied by Nerger et al. (2012a), where also the ESTKF has been derived. Apart from computing the ensemble transformation in the error subspace in case of SEIK and ESKTF, the three filters are essentially equivalent. However, for the SEIK filter it has been found that the application of the Cholesky decomposition can lead to unevenly distributed variance in the ensemble members. To this end, the ESTKF and ETKF method are preferable unless a random matrix $\mathrm{\Omega }$ is used in the SEIK filter.

The methods described in this section each have nuances in which they differ one from another as well as underlying common ground. Writing these methods using unified mathematics notation allows us to see these algorithmic differences and commonalities more readily. Many filters described above have several common variables and while for some methods the variables have different sizes, others can not only have the same size but actually the same value, too. Table 1 summarises the sizes of the common variables between the methods and below we comment on whether they have the same value.

Table 1. Overview of the sizes of matrices that are used in different filter methods.

Variable	EnKF	SEIK	ESSE	ETKF	EAKF	EnSRF	SEnKF	ESTKF
$\mathsf{T}$	$N_e\times N_e$	$N_e-1\times N_e-1$	${\tilde{N}}_e\times {\tilde{N}}_e$	$N_e\times N_e$	$N_e\times N_e$	$N_e\times N_e$	$N_e\times N_e$	$N_e-1\times N_e-1$
$\mathsf{U}$	–	–	$N_x\times N_e$	$N_e\times N_e$	$N_e\times N_e$	$N_e\times N_e$	$N_e\times N_e$	$N_e-1\times N_e-1$
$\mathrm{\Sigma }$	–	–	$N_e\times N_e$	$N_e\times N_e$	$N_e\times N_y$	$N_e\times min(N_e,N_y)$	$N_e\times N_y$	$N_e-1\times N_e-1$
$\mathsf{V}$	$N_y\times N_y$	–	$N_e\times N_e$	–	$N_y\times N_y$	$min(N_e,N_y)\times N_y$	–	–
$\mathsf{L}$	–	$N_x\times N_e-1$	–	–	–	–	–	$N_x\times N_e-1$
$\mathsf{A}$	–	$N_e\times N_e-1$	–	–	–	–	–	$N_e\times N_e-1$

Apart from the matrices listed in Table 1, there are the finally resulting weight matrices ${\mathsf{W}}^{\prime }$ and $\overline{\mathsf{W}}$. The matrix $\overline{\mathsf{W}}$ is identical for all filters. Thus, for all filters except EnKF, SEnKF and MLEF, which don’t use this matrix, the mean of the analysis ensemble is identical. The EnKF and SEnKF are an exception because they do not explicitly transform the ensemble mean and introduce sampling error due to the perturbed observations. The maximum likelihood approach of the MLEF also results in a different analysis state estimate if the ensemble distribution is not Gaussian and the observation operator is non-linear. Further the square ${\mathsf{W}}^{\prime }{{\mathsf{W}}^{\prime }}^T$ is identical for all filters except the EnKF, SEnKF and MLEF. Thus, the analysis covariance matrix ${\mathsf{P}}^a$ will be identical.

In contrast to the equality of the matrix $\overline{\mathsf{W}}$, the matrix ${\mathsf{W}}^{\prime }$ is different for almost all methods. Thus, while many methods yield the same analysis ensemble covariance matrix, their ensemble perturbations are distinct. In the ETKF and ESTKF methods, the dimensions of the matrices $\mathsf{T}$, $\mathsf{U}$, and $\mathrm{\Sigma }$ have distinct dimensions. However, the ensemble transformation weight matrices ${\mathsf{W}}^{\prime }$ of both methods are identical (Nerger et al., 2012a).

In general, the choice of using a particular ensemble method depends on a number of the system’s components: the dynamical model at hand, model error, number and types of observations, ensemble size. Given these degrees of freedom it is not possible to attribute one data assimilation method to be better suited for a general situation. In practise operational meteorological applications most widely use the LETKF and the serial formulations of the EnSRF and EAKF, while in oceanography there are many applications of the stochastic EnKF, the SEIK filter, and the ESTKF. There are less applications of the ESSE and MLEF, despite the fact that these filters are algorithmically interesting because of the ensemble-size adaptivity of ESSE and the maximum-likelihood solution MLEF. From the algorithmic viewpoint, the stochastic EnKF will be useful if the stochasticity can be an advantage and if large ensembles can be used. Further, the filters SEIK, ETKF, ESTKF differ from the EnKF, EnSRF and EAKF also in the application of distinct localisation methods (see Section 9.1 for the discussion of localisation). EnKF, EnSRF and EAKF allow for localisation in the state space, which could be advantageous for some observation types (Campbell et al., 2010). The serial formulation of the EnSRF and EAKF requires that the observation error covariance matrix is diagonal. Thus, these filters cannot directly be applied if the observation errors are correlated. A transformation into variables that are uncorrelated is possible in theory, but it is most likely not practical for large sets of observations.

6. Particle filters

In this section we will consider the standard particle filter followed by three efficient variants of the particle filters: the Equivalent Weights Particle Filter (EWPF, van Leeuwen, 2010;van Leeuwen, 2011; Ades and van Leeuwen, 2013), the Implicit Equal-Weights Particle Filter (Zhu et al., 2016) and the Ensemble Transform Particle Filter (Reich, 2013). Other variants of local particle filters are discussed in Section 9 on practical implementation. Another interesting particle filter for high-dimensional systems, the so called implicit particle filter (Chorin and Tu, 2009; Chorin et al., 2010; Morzfeld et al., 2012; van Leeuwen et al., 2015), is not discussed here as it needs a 4D-Var-like minimisation for each particle. The Multivariate Rank Histogram Filter (MRHF, Metref et al., 2014b), based on the Rank-Histogram Filter of Anderson (2010) that performs well in highly non-Gaussian regimes, has been recently developed in the European project Sangoma.⁴ However, it is still under development for high-dimensional systems and its idea is only shortly described in Section 9.1.5. Often particle filters are defined as providing approximations of $p\left({\mathsf{x}}^{(0:m)}|{\mathsf{y}}^{(1:m)}\right)$, but we restrict ourselves to particle filters that are approximations of the marginal posterior pdf $p\left({\mathsf{x}}^{(m)}|{\mathsf{y}}^{(1:m)}\right)$ as there are at present no efficient algorithms for the former for high-dimensional geophysical systems, and we have forecasting in mind. Furthermore, for ease of presentation we take all earlier observations for granted, leading to the marginal posterior at time m being denoted as $p\left({\mathsf{x}}^{(m)}|{\mathsf{y}}^{(m)}\right)$.

6.1. The standard particle filter

This particle filter is also known as the bootstrap filter or Sequential Importance Resampling (SIR). The probability distribution function (pdf) in particle filtering, represented by $N_e$ particles or ensemble members at time k, is given by(91) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(m)}\right)=\frac{1}{N_e}\sum _{j=1}^{N_e}\delta \left({\mathsf{x}}^{(m)}-{\mathsf{x}}_j^{(m)}\right),\end{array}$$(91)

where ${\mathsf{x}}^{(m)}\in {\mathcal{R}}^{N_x}$ is the $N_x$-dimensional state of the system that has been integrated forward in time using the stochastic forward model and $\delta (x)$ is a Dirac-delta function. We let time m to be the time of a current set of observations with the previous observation set at time 0. Then the stochastic forward model for times $0<k<m$ for each particle $j=1,\dots ,N_e$ is given by(92) $$\begin{array}{c}\hfill {\mathsf{x}}_j^{(k)}={\mathcal{M}}_k\left({\mathsf{x}}_j^{(k-1)}\right)+{\beta }_j^{(k)},\end{array}$$(92)

where ${\beta }_j^{(k)}\in {\mathcal{R}}^{N_x}$ are random terms representing the Gaussian distributed model errors with mean zero and covariance matrix $\mathsf{Q}$, and ${\mathcal{M}}_k:{\mathcal{R}}^{N_x}\to {\mathcal{R}}^{N_x}$ is the deterministic model from time $k-1$ to k. Thus, the model state transition from time $k-1$ to k is fully described by the transition density given by(93) $$\begin{array}{c}\hfill p\left({\mathsf{x}}_j^{(k)}|{\mathsf{x}}_j^{(k-1)}\right)=\mathcal{N}\left({\mathcal{M}}_k\left({\mathsf{x}}^{(k-1)}\right),\mathsf{Q}\right),\end{array}$$(93)

which will be of later use.

Using Bayes theorem(94) $$\begin{array}{c}\hfill p\left({\mathsf{x}}_j^{(m)}|{\mathsf{y}}^{(m)}\right)=\frac{p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m)}\right)}{p\left({\mathsf{y}}^{(m)}\right)}p\left({\mathsf{x}}_j^{(m)}\right)\end{array}$$(94)

and the Markovian property of the model, the full posterior at observation time m is written as(95) $$\begin{array}{c}\hfill p\left({\mathsf{x}}_j^{(m)}|{\mathsf{y}}^{(m)}\right)=\sum _{j=1}^{N_e}{w}_j^{(m)}\delta \left({\mathsf{x}}^{(m)}-{\mathsf{x}}_j^{(m)}\right)\end{array}$$(95)

where the weights $w_j^{(m)}$ are given by(96) $$\begin{array}{c}\hfill w_j^{(m)}\propto p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m)}\right)p\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)w_j^{(m-1)}\end{array}$$(96)

and each $w_j^{(m-1)}$ is the product of all the weights from all time steps $0<k\le m-1$. The conditional pdf $p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}^{(m)}\right)$ is the pdf of the observations given the model state ${\mathsf{x}}^{(m)}$ which is often taken to be Gaussian(97) $$\begin{array}{cc}\hfill p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}^{(m)}\right)& \propto exp\left[-\frac{1}{2}{\left({\mathsf{y}}^{(m)}-{\mathcal{H}}_m\left({\mathsf{x}}^{(m)}\right)\right)}^{\mathsf{T}}\right.\hfill \\ \hfill & \left.{\mathsf{R}}^{-1}\left({\mathsf{y}}^{(m)}-{\mathcal{H}}_m\left({\mathsf{x}}^{(m)}\right)\right)\right].\hfill \end{array}$$(97)

To obtain equal-weight posterior particles one applies resampling, in which particles with high weights are duplicated, while particles with low weights are abandoned. Several schemes have been developed to perform resampling, and three of the most-used schemes are presented in Appendix 1.

The problem in high-dimensional spaces with a large number of independent observations is that these weights vary enormously over the particles, with one particle obtaining a weight close to one, while all the others have a weight very close to zero. This is the so-called degeneracy problem related to the ‘curse of dimensionality’: any resampling scheme will produce $N_e$ copies of the particle with the highest weight, and all variation in the ensemble has disappeared.

Hence, as mentioned at the beginning of this section, to apply a particle filter to a high-dimensional system additional information is needed to limit the search space of the filter. One option is to use localisation directly on the standard particle filter. Local particle filters, like the so-called Local Particle Filter (Poterjoy, 2016a) will be discussed in the section on localisation in particle filters. We next discuss the proposal-density particle filters since this technique could be applied to all filters and permits us to achieve equal-weights for the particles in a different way.

6.2. Proposal-density particle filters

To avoid that the ensemble degenerates we aim at ensuring that equally significant particles are picked from the posterior density. To do this we have to ensure that all particles end up in the high-probability area of the posterior pdf, and that they have very similar, or even equal, weights. For the former we can use a scheme that pulls the particles towards the observations. Several methods can be used for this, including traditional methods like 4DVar, a variational method, and ensemble Kalman filters and smoothers. However, the main ingredient in efficient particle filters is the step that ensures that the weights of the different particles are close before any resampling step.

We start by writing the prior at time m as follows:(98) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(m)}\right)=\int p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}^{(m-1)}\right)p\left({\mathsf{x}}^{(m-1)}\right)d{\mathsf{x}}^{(m-1)}.\end{array}$$(98)

Without loss of generality but for simplicity we assume that the particle weights in the ensemble at the previous time step $m-1$ are equal, so(99) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(m-1)}\right)=\frac{1}{N_e}\sum _{j=1}^{N_e}\delta \left({\mathsf{x}}^{(m-1)}-{\mathsf{x}}_j^{(m-1)}\right).\end{array}$$(99)

Using Equation (99(99) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(m-1)}\right)=\frac{1}{N_e}\sum _{j=1}^{N_e}\delta \left({\mathsf{x}}^{(m-1)}-{\mathsf{x}}_j^{(m-1)}\right).\end{array}$$(99) ) in Equation (98(98) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(m)}\right)=\int p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}^{(m-1)}\right)p\left({\mathsf{x}}^{(m-1)}\right)d{\mathsf{x}}^{(m-1)}.\end{array}$$(98) ) leads directly to:(100) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(m)}\right)=\frac{1}{N_e}\sum _{j=1}^{N_e}p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right),\end{array}$$(100)

hence, from Equation (93(93) $$\begin{array}{c}\hfill p\left({\mathsf{x}}_j^{(k)}|{\mathsf{x}}_j^{(k-1)}\right)=\mathcal{N}\left({\mathcal{M}}_k\left({\mathsf{x}}^{(k-1)}\right),\mathsf{Q}\right),\end{array}$$(93) ) the prior can be seen as a mixture density, with each density centred around one of the forecast particles.

One can now multiply the numerator and denominator of Equation (100(100) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(m)}\right)=\frac{1}{N_e}\sum _{j=1}^{N_e}p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right),\end{array}$$(100) ) by the same factor $q\left(x^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)$, in which ${\mathsf{x}}_{1:N_e}^{(m-1)}$ is defined as the collection of all particles at time $m-1$, and the conditioning on j denotes that each particle does in general has a different parent to start from. This leads to(101) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^{(m)}\right)=\frac{1}{N_e}\sum _{j=1}^{N_e}\frac{p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)}{q\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)}q\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)\end{array}$$(101)

where $q\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)$ is the so-called proposal transition density, or proposal for short, whose support should be equal to or larger than that of $p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}^{(m-1)}\right)$. Note that the proposal density as formulated here is slightly more general than the usual $q\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_j^{(m-1)},{\mathsf{y}}^{(m)}\right)$ through allowing for the explicit dependence on all particles at time $m-1$.

Drawing from this density we find for the posterior:(102) $$\begin{array}{cc}\hfill p\left({\mathsf{x}}^{(m)}|{\mathsf{y}}^{(m)}\right)& =\frac{p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m)}\right)}{p\left({\mathsf{y}}^{(m)}\right)}p\left({\mathsf{x}}^{(m)}\right)\hfill \\ \hfill & =\frac{1}{N_e}\sum _{j=1}^{N_e}{w}_j\delta \left({\mathsf{x}}^{(m)}-{\mathsf{x}}_j^{(m)}\right)\hfill \end{array}$$(102)

where $w_j$ are the particle weights given by(103) $$\begin{array}{c}\hfill w_j=\frac{p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m)}\right)}{p\left({\mathsf{y}}^{(m)}\right)}\frac{p\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)}{q\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)}.\end{array}$$(103)

Using Bayes’ theorem, the numerator in the expression for the weights can be expressed as(104) $$\begin{array}{cc}& p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}^{(m)}\right)p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)\hfill \\ \hfill & =p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_j^{(m-1)},{\mathsf{y}}^{(m)}\right)p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)\hfill \end{array}$$(104)

Therefore, the particle weight of ensemble member j can be written as:(105) $$\begin{array}{c}\hfill w_j=\frac{p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)}{p\left({\mathsf{y}}^{(m)}\right)}\frac{p\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_j^{(m-1)},{\mathsf{y}}^{(m)}\right)}{q\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)}.\end{array}$$(105)

In the so-called optimal proposal density (Doucet et al., 2000) one chooses$$\begin{array}{c}\hfill q\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)=p\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_j^{(m-1)},{\mathsf{y}}^{(m)}\right),\end{array}$$

leading to weights $w_j\propto p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)$. For systems with a large number of independent observations these weights are again degenerate (see, e.g. Snyder et al., 2008; Ades and van Leeuwen, 2013; Snyder et al., 2015).

Several efficient particle filter schemes have been developed utilising the proposal density to avoid this degeneracy. Here we discuss the Equivalent-Weights Particle Filter (EWPF) and the Implicit Equal-Weights Particle Filter (IEWPF). As mentioned, the Implicit Particle Filter (Chorin et al., 2010), which allows for an extension of the one-time-step optimal proposal particle filter to a full time window explores a 4DVar-like method on each particle. Since it needs an adjoint of the underlying model, it is not discussed in this paper.

6.2.1. The equivalent-weights particle filter

The EWPF works as follows:

(1)	Determine the optimal proposal weight $w_j\propto p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)$ for each particle. Note that these weights vary enormously in high-dimensional systems.
(2)	Choose a target weight $w_{\text{target}}$ based on these weights that a certain percentage of particles can reach. For instance, if the target weight is set to the lowest of these weights we keep 100% of the particles. A choice of 50% will mean that the target weight is set to the medium value of these weights.
(3)	Calculate the position in state space of each particle such that it has a weight exactly equal to the target weight. This is where the proposal density comes in. Note that some of the particles cannot reach this target weight no matter how we move them, and these are brought back into the ensemble via the resampling step in point 5.
(4)	Add a small random perturbation to each particle and recalculate its weight.
(5)	Resample all particles such that their weights are equal again.

It is in step 3 that we use the fact that the proposal density is dependent on all previous particles, and not just particle j. This step is the main reason for the efficiency of the filter.

As an example, when the error in the model equations is additive Gaussian and the observation operator is linear an analytical solution can be found for the maximum weight for each particle j, or actually, the minimum of minus the $log$ of that weight called ${\varphi }_j$:(106) $$\begin{array}{cc}\hfill {\varphi }_j& ={\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T{\left(\mathsf{H}Q{H}^T+\mathsf{R}\right)}^{-1}\hfill \\ \hfill & \times \left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right).\hfill \end{array}$$(106)

Then a target weight is set from these ${\varphi }_j$’s. The target weight splits the ensemble of particles in two groups: those particles that have a higher optimal proposal weight, and those with a lower optimal proposal weight. The latter are abandoned at this point, and will be regenerated in the resampling step 5.

For the retained particles, there is an infinite number of ways to move a particle in state space such that it reaches the target weight. In the EWPF that problem is solved by assuming(107) $$\begin{array}{c}\hfill {\widehat{\mathsf{x}}}_j^{(m)}=\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)+{\alpha }_j\mathrm{{\rm Y}}\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)\end{array}$$(107)

in which ${\alpha }_j$ is a scalar, and $\mathrm{{\rm Y}}$ is defined as(108) $$\begin{array}{c}\hfill \mathrm{{\rm Y}}=\mathsf{Q}{\mathsf{H}}^T{\left(\mathsf{H}\mathsf{Q}{\mathsf{H}}^T+\mathsf{R}\right)}^{-1}.\end{array}$$(108)

Under this assumption the number of solutions is reduced to two, and the two values for ${\alpha }_j$ are given by(109) $$\begin{array}{c}\hfill {\alpha }_j=1\pm \sqrt{1-b_j/a_j}\end{array}$$(109)

in which(110) $$\begin{array}{cc}\hfill a_j& =\frac{1}{2}{\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T\hfill \\ \hfill & \times {\mathsf{R}}^{-1}{\mathsf{H}}^T\mathrm{{\rm Y}}\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)\hfill \end{array}$$(110)

and(111) $$\begin{array}{cc}\hfill b_j& =\frac{1}{2}{\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T\hfill \\ \hfill & \times {\mathsf{R}}^{-1}\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)+\hfill \\ \hfill & log{w}_{\text{target}}-log{w}_j^{(m-1)},\hfill \end{array}$$(111)

in which $w_j^{(m-1)}$ is the weight of particle j accumulated over previous time steps, included here for completeness. Note that $w_{\text{target}}$ is the target weight selected from $\varphi$’s in Equation (106(106) $$\begin{array}{cc}\hfill {\varphi }_j& ={\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T{\left(\mathsf{H}Q{H}^T+\mathsf{R}\right)}^{-1}\hfill \\ \hfill & \times \left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right).\hfill \end{array}$$(106) ) (e.g. if we choose to keep 80% particles $-log(w_{\text{target}})={\{{\tilde{\varphi }}_j\}}_{j=⌈N_e\ast 0.8⌉}$ where ${\{{\tilde{\varphi }}_j\}}_{j=1,\dots ,N_e}$ is a sorted list of optimal-proposal weight of each particle) and that ${\alpha }_j=1$ pushes the particle to its optimal-proposal weight position. The solution resembles the optimal proposal solution in which the deterministic part of the proposal is scaled to ensure equal weights. Also note the resemblance of the deterministic part with the shape of that used in a Kalman filter when we replace Q with the ensemble covariance of the state.

When the number of independent observations is large the optimal proposal density particle filter is degenerate, meaning that one particle gets a much larger weight than all the others. The EWPF is not degenerate because a set percentage of all particles has a similar weight (before the resampling step). The EWPF does not, however, converge for large $N_e$ to the posterior pdf because of this equivalent-weights construction, in which high-weight particles are moved such that their weight becomes lower, equal to the target weight. So the scheme is biased. However, the large $N_e$ limit is not that relevant in practise as the affordable number of particles will be low, below say 10,000, and typically of $O(20-100)$. In that setting, the Monte-Carlo error will be substantial, and the bias should be measured against the Monte-Carlo error. As long as the latter is larger than the former the scheme is a valid alternative in high-dimensional systems.

Algorithm  in Appendix 2 gives a pseudo-algorithm for the EWPF.

6.2.2. The implicit equal-weights particle filter

This scheme is very similar to that of the EWPF:

(1)	Determine the optimal proposal weight $w_j\propto p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)$ for each particle. Note that these weights vary enormously in high-dimensional systems.
(2)	Choose a target weight based on these weights that a certain percentage of particles can reach. Typically the target weight is chosen as the minimum of the maximal weights, so that all particles are kept.
(3)	Draw a random perturbation vector for each particle, and add this to the particle position that leads to maximal weight. So far the scheme is the same as that used in the optimal proposal density.
(4)	Scale each random vector such that each particle will reach the target weight.
(5)	Resample the particles such that their weights are equal in case the kept percentage is lower than 100%.

The main difference between this scheme and the EWPF is that in the EWPF we scale the deterministic part of the optimal proposal to reach a target weight, while here we scale the random part of the optimal proposal.

The implicit part of our scheme follows from drawing samples implicitly from a standard Gaussian distributed proposal density $q\left(\xi \right)$ instead of the original one $q\left({\mathsf{x}}^{(m)}|{\mathsf{x}}^{(m-1)},{\mathsf{y}}^{(m)}\right)$, as in (Chorin and Tu, 2009). These two pdfs are related by:(112) $$\begin{array}{c}\hfill q\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)=\frac{q\left({\xi }_j\right)}{\left||{\frac{\mathrm{d}\mathsf{x}}{\mathrm{d}\xi }}_j\right||}\end{array}$$(112)

where $\left||\frac{\mathrm{d}\mathsf{x}}{\mathrm{d}{\xi }_j}\right||$ denotes the absolute value of the determinant of the Jacobian matrix of the ${\mathcal{R}}^{N_x}\to {\mathcal{R}}^{N_x}$ transformation ${\mathsf{x}}_j=g_j({\xi }_j)$. The transformation $g_j(.)$ is now defined via the following implicit relation between variable ${\mathsf{x}}_j^{(m)}$ and $\xi$ as(113) $$\begin{array}{c}\hfill {\mathsf{x}}_j^{(m)}={\mathsf{x}}_j^a+{\alpha }_j^{1/2}{\mathsf{P}}^{1/2}{\xi }_j^{(m)}\end{array}$$(113)

where ${\mathsf{x}}_j^a$ is the mode of $p\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_j^{(m-1)},{\mathsf{y}}^{(m)}\right)$, $\mathsf{P}$ a measure of the width of that pdf, and ${\alpha }_j$ a scalar that depends on ${\xi }_j^{(m)}$.

The ${\alpha }_j$ are now chosen such that all particles get the same weight $w_{\text{target}}$, so the scalar ${\alpha }_j$ is determined for each particle from:(114) $$\begin{array}{cc}\hfill w_j& =\frac{p\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_j^{(m-1)},{\mathsf{y}}^{(m)}\right)p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)}{q\left({\xi }_j\right)}\left||\frac{\mathrm{d}\mathsf{x}}{\mathrm{d}{\xi }_j}\right||w_j^{(m-1)}\hfill \\ \hfill & =w_{\text{target}}\hfill \end{array}$$(114)

This ensures that the filter is not degenerate in systems with arbitrary dimensions and an arbitrary number of independent observations. Because of the target-weight construction the filter does not converge to the correct posterior pdf, and the same discussion as for the EWPF applies here, namely that as long as this bias is smaller than the Monte-Carlo error this filter is a valid candidate for high-dimensional non-linear filtering.

As an example we assume now that observation errors and model errors are Gaussian and that the observation operator $\mathsf{H}\in {\mathcal{R}}^{N_y\times N_x}$ is linear. Then we find that(115) $$\begin{array}{cc}& p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}^{(m)}\right)p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)\hfill \\ \hfill & =\frac{1}{A}exp\left[-\frac{1}{2}{\left({\mathsf{y}}^{(m)}-\mathsf{H}{\mathsf{x}}^{(m)}\right)}^T{\mathsf{R}}^{-1}\left({\mathsf{y}}^{(m)}-\mathsf{H}{\mathsf{x}}^{(m)}\right)\right.\hfill \\ \hfill & -\frac{1}{2}\left.{\left({\mathsf{x}}^{(m)}-\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T{\mathsf{Q}}^{-1}\left({\mathsf{x}}^{(m)}-\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)\right]\hfill \\ \hfill & =\frac{1}{A}exp\left[-\frac{1}{2}{\left({\mathsf{x}}^{(m)}-{\mathsf{x}}_j^a\right)}^T{\mathsf{P}}^{-1}\left({\mathsf{x}}^{(m)}-{\mathsf{x}}_j^a\right)\right]exp\left(-\frac{1}{2}{\varphi }_j\right)\hfill \\ \hfill & =p\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_j^{(m-1)},{\mathsf{y}}^{(m)}\right)p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)\hfill \end{array}$$(115)

where(116) $$\begin{array}{cc}\hfill \mathsf{P}& ={\left({\mathsf{Q}}^{-1}+{\mathsf{H}}^T{\mathsf{R}}^{-1}\mathsf{H}\right)}^{-1},\hfill \end{array}$$(116) (117) $$\begin{array}{cc}\hfill {\mathsf{x}}_j^a& =\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)+\mathrm{{\rm Y}}\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right),\hfill \end{array}$$(117)

and(118) $$\begin{array}{cc}\hfill {\varphi }_j& ={\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T{\left(\mathsf{H}\mathsf{Q}{\mathsf{H}}^T+\mathsf{R}\right)}^{-1}\hfill \\ \hfill & \times \left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right).\hfill \end{array}$$(118)

This leads to a complicated non-linear differential equation for ${\alpha }_j$ that involves the determinant of $\mathsf{P}$. Since we are interested in high-dimensional problems we consider this equation in the limit of large state dimension $N_x$. In that limit it turns out that we can integrate this equation, leading to the much simpler equation (see Appendix in Zhu et al. (2016)):(119) $$\begin{array}{c}\hfill ({\alpha }_j-1){\gamma }_j-N_{x}log({\alpha }_j)+{\varphi }_j-log{w}_j^{(m-1)}=log{w}_{\text{target}}.\end{array}$$(119)

in which ${\gamma }_j={\xi }_j^T{\xi }_j$. This equation could be approximated by using numerical methods, such as the Newton method, etc., but analytical solutions based on the so-called Lambert W function do exist. We do not elaborate on these here.

Algorithm  in Appendix 2 gives a pseudo-algorithm for the IEWPF.

6.2.3. Between observations: relaxation steps

If the system is not observed at every time step, the schemes mentioned above can be used over the time window between observations. No analytical solutions can be obtained in this case so that the solution has to be found iteratively. However, this procedure is rather expensive as it typically involves solving a problem similar to a 4DVar on each particle.⁵ Thus, typically simpler schemes are employed between observation times. These schemes will be less efficient, although we can ensure that Bayes’ theorem is fulfilled exactly for each particle.

In the following, we demonstrate the use of relaxation between observation times. We use the future observations to relax the particles at time k towards observations at next time $m>k$ by using instead of Equation (92(92) $$\begin{array}{c}\hfill {\mathsf{x}}_j^{(k)}={\mathcal{M}}_k\left({\mathsf{x}}_j^{(k-1)}\right)+{\beta }_j^{(k)},\end{array}$$(92) ) the modified forward model(120) $$\begin{array}{cc}\hfill {\mathsf{x}}_j^{(k)}& ={\mathcal{M}}_k\left({\mathsf{x}}_j^{(k-1)}\right)+{\tilde{\beta }}_j^{(k)}+\hfill \\ \hfill & +\widetilde{\mathrm{{\rm Y}}}\left[{\mathsf{y}}^{(m)}-{\mathcal{H}}_k\left({\mathsf{x}}_j^{(k-1)}\right)\right],\hfill \end{array}$$(120)

where ${\tilde{\beta }}_j^{(k)}\in {\mathcal{R}}^{N_x}$ are random terms representing the model error distributed according to a given covariance matrix $\tilde{\mathsf{Q}}$,⁶${\mathcal{M}}_k$ is the same deterministic model as in Equation (92(92) $$\begin{array}{c}\hfill {\mathsf{x}}_j^{(k)}={\mathcal{M}}_k\left({\mathsf{x}}_j^{(k-1)}\right)+{\beta }_j^{(k)},\end{array}$$(92) ), $\widetilde{\mathrm{{\rm Y}}}$ is a relaxation matrix given by(121) $$\begin{array}{c}\hfill \widetilde{\mathrm{{\rm Y}}}=\tau (k)\mathsf{Q}{\mathsf{H}}^T{\mathsf{R}}^{-1}.\end{array}$$(121)

Here, $\tau (k)$ is a time dependent scalar that determines the strength of the relaxation. ${\mathsf{y}}^{(m)}\in {\mathcal{R}}^{N_y}$ is the vector of $N_y$ observations at time m and ${\mathcal{H}}_k:{\mathcal{R}}^{N_x}\to {\mathcal{R}}^{N_y}$ is the observation operator mapping model space into observation space. Note that the observations ${\mathsf{y}}^{(m)}$ exist at the later time $m>k$. The modified transition density is now given by(122) $$\begin{array}{cc}\hfill q\left({\mathsf{x}}_j^{(k)}|{\mathsf{x}}_j^{(k-1)},{\mathsf{y}}^{(m)}\right)=& \mathcal{N}\left({\mathcal{M}}_k\left({\mathsf{x}}^{(k-1)}\right)\right.\hfill \\ \hfill & \left.+\widetilde{\mathrm{{\rm Y}}}\left[{\mathsf{y}}^{(m)}-{\mathcal{H}}_k\left({\mathsf{x}}^{(k-1)}\right)\right],\mathsf{Q}\right),\hfill \end{array}$$(122)

and the modified weights $w_j^{(k)}$ are accumulated as(123) $$\begin{array}{cc}\hfill w_j^{(k)}& \propto \frac{p\left({\mathsf{x}}_j^{(k)}|{\mathsf{x}}_j^{(k-1)}\right)}{q\left({\mathsf{x}}_j^{(k)}|{\mathsf{x}}_j^{(k-1)},{\mathsf{y}}^{(m)}\right)}{w}_j^{(k-1)}\hfill \\ \hfill & \propto \prod _{t=1}^k\frac{p\left({\mathsf{x}}_j^{(t)}|{\mathsf{x}}_j^{(t-1)}\right)}{q\left({\mathsf{x}}_j^{(t)}|{\mathsf{x}}_j^{(t-1)},{\mathsf{y}}^{(m)}\right)}.\hfill \end{array}$$(123)

This simple modification of the forward model to include information about future observations using a relaxation term is only consistent with Bayes Theorem when the weights that are introduced by this modification are properly taken into account, and it leads to efficient schemes if it is used in combination with an equal-weight scheme, like the EWPF or the IEWPF. Algorithm  in Appendix 2 gives a pseudo-algorithm of the relaxation step used in the EWPF and the IEWPF.

Note that it would also be possible to use other methods like ensemble smoothers or ensemble 4Dvar-like methods to move particles between observations, but we will not elaborate on those here.

6.3. The Ensemble Transform Particle Filter (ETPF)

The idea of the Ensemble Transport Particle Filter (Reich, 2013) is to avoid resampling by finding a linear transportation map between the prior and the posterior ensemble such that the prior particles are minimally modified, while ensuring that the posterior particles have equal weight. We write each posterior particle as a linear combination of the prior particles as(124) $$\begin{array}{c}\hfill {\mathsf{x}}_j^a=N_e\sum _{i=1}^{N_e}{\mathsf{x}}_i^f{t}_{\mathit{ij}}\end{array}$$(124)

in which we ensure that the particles have the correct mean via(125) $$\begin{array}{c}\hfill \sum _{i=1}^{N_e}{t}_{\mathit{ij}}=\frac{1}{N_e},\sum _{j=1}^{N_e}{t}_{\mathit{ij}}=w_i.\end{array}$$(125)

This still leads to $N_e^2-2$ undetermined elements $t_{\mathit{ij}}$. These are found by minimising the movement from old to new particles, by minimising(126) $$\begin{array}{c}\hfill J(\mathsf{T})=\sum _{i,j}^{N_e}{t}_{\mathit{ij}}\left||{\mathsf{x}}_i^f-{\mathsf{x}}_j^f\right|{|}^2\end{array}$$(126)

under the condition that $t_{\mathit{ij}}\ge 0$. The above formulation is an example of an optimal transportation algorithm, see e.g. the review by Chen and Reich in van Leeuwen et al. (2015). This scheme can be combined with any proposal density discussed in the previous section.

If the dynamical model is deterministic one needs to add some small random noise to the particles to avoid ensemble collapse. Typically this noise is assumed to be Gaussian with zero mean and covariance ${\mathsf{P}}^a=h^2{\mathsf{P}}^f$ with $0<h<1$ a free parameter. This term is an ad-hoc addition related to inflation in Ensemble Kalman Filters.

Algorithm  in Appendix 2 gives a pseudo-algorithm of the ETPF.

7. Second-order exact ensemble Kalman filters

Several extensions to ensemble Kalman filters have been proposed to overcome the linearity or Gaussianity assumptions. A large number of filters exists that try to bridge an ensemble Kalman filter and particle filter by defining smoothly varying parameters that move the filter between these two extremes based on the degeneracy of the particle filter. In high-dimensional systems, however, all of these filters become ensemble Kalman filters as any particle filter contribution results in complete degeneracy. These filters (not discussed here) will become useful when localisation is applied.

In a non-linear, non-Gaussian case the ensemble Kalman filters will necessarily produce an analysis where the mean and covariance are biased due to the assumption of a Gaussian prior pdf (Lei and Bickel, 2011). Here, we will discuss four ensemble filters that concentrate on getting the first two moments of the posterior distribution correct in non-linear situations. These are the Particle Filter with Gaussian Resampling of Xiong et al. (2006), the Non-linear Ensemble Transform Filter (Tödter and Ahrens, 2015), the Moment-Matching Ensemble Filter (Lei and Bickel, 2011), and the Merging Particle Filter (Nakano et al., 2007).

7.1. Particle Filter with Gaussian Resampling (PFGR) and Non-linear Ensemble Transform Kalman Filter (NETF)

The Particle Filter with Gaussian Resampling (PFGR, Xiong et al. (2006)) introduced an explicit ensemble transformation matching the mean and covariance matrix. The Non-linear Ensemble Transform Filter (NETF, Tödter and Ahrens, 2015) is a recent reinvention of this algorithm formulated to obtain an ensemble transformation that is analogous to that of the ETKF. In addition, the NETF was introduced with localisation, so that the filter can be applied to high-dimensional systems (see Section 9.1). The presentation here follows the more modern formulation of the NETF in analogy to the ETKF presented before. As a novel feature, the presented formulation avoids the explicit computation of the analysis state, that is given by the weighted ensemble mean.

The PFGR and the NETF are designed to exactly match the first two moments of the posterior pdf in Bayes theorem without assuming that the prior or likelihood are normally distributed. The forecast ensemble is transformed into an analysis ensemble by applying a weight vector to obtain the analysis mean state and a transform matrix to obtain analysis ensemble perturbations, analogous in form to a square root filter (Equations (17(17) $$\begin{array}{cc}\hfill {\overline{\mathsf{x}}}^a& ={\overline{\mathsf{x}}}^f+{{\mathsf{X}}^{\prime }}^f\overline{\mathsf{w}},\hfill \end{array}$$(17) ) to (19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) )).

As in most particle filters, the likelihood weights that arise from Bayes’ theorem(127) $$\begin{array}{c}\hfill w_j=\frac{p\left(\mathsf{y}|{\mathsf{x}}_j\right)}{{\sum }_{k=1}^{N_e}p\left(\mathsf{y}|{\mathsf{x}}_k\right)}\end{array}$$(127)

are used. For normally distributed observation errors, the weight of each member is at first given by(128) $$\begin{array}{c}\hfill w_j\propto exp\left[-\frac{1}{2}{\left(\mathsf{y}-\mathcal{H}\left({\mathsf{x}}_j^f\right)\right)}^T{\mathbf{R}}^{-1}\left(\mathsf{y}-\mathcal{H}\left({\mathsf{x}}_j^f\right)\right)\right]\end{array}$$(128)

and then normalised so that the weights sum up to one. Before the weights are computed, the ensemble perturbations should be inflated by an inflation factor $\gamma >1$ as in the ensemble-based Kalman filters (for inflation see Section 9.2). Using the weight vector $\mathsf{w}={\left(w_1,\dots ,w_{N_e}\right)}^T$ the transform matrix is(129) $$\begin{array}{cc}\hfill \mathsf{T}{\mathsf{T}}^T& =N_e\left[\mathrm{d}iag(\mathsf{w})-\mathsf{w}{\mathsf{w}}^T\right].\hfill \end{array}$$(129)

Here, $\mathrm{d}iag(\mathsf{w})$ is a diagonal matrix that contains the weights $w_j$ on the diagonal. The factor $N_e$ was not present in the formulation by Xiong et al. (2006). It was introduced by Tödter and Ahrens (2015) to ensure that the ensemble has the correct analysis variance. As in the ensemble Kalman filters, the eigenvalue decomposition of $\mathsf{T}{\mathsf{T}}^T=\mathsf{U}\mathrm{\Sigma }{\mathsf{U}}^T$ yields the ensemble transformation(130) $$\begin{array}{cc}\hfill \mathsf{T}& =\mathsf{U}{\mathrm{\Sigma }}^{1/2}{\mathsf{U}}^T.\hfill \end{array}$$(130)

Combining the weight vector and transform matrix as in Equation (19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ), the analysis ensemble is given by(131) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\mathsf{X}}^f\left(\mathsf{T}\mathrm{\Lambda }+[\mathsf{w},\dots ,\mathsf{w}]\right).\end{array}$$(131)

Here, $\mathrm{\Lambda }$ is an random matrix. Xiong et al. (2006) use a random matrix sampled from a normal distribution with mean zero and standard deviation one. They use $\mathrm{\Lambda }$ because in Equation (130(130) $$\begin{array}{cc}\hfill \mathsf{T}& =\mathsf{U}{\mathrm{\Sigma }}^{1/2}{\mathsf{U}}^T.\hfill \end{array}$$(130) ) they omit all eigenvalues that are very close to zero and need to restore an ensemble of full size. In contrast, Tödter and Ahrens (2015) use a mean-preserving orthogonal matrix (see Pham, 2001) analogous to that used in the SEIK filter. They motivate the use of $\mathrm{\Lambda }$ also by the reduction of ensemble outliers and showed experimentally that the random transformation with mean preserving properties leads to a more stable data assimilation process.

Note, that the transformation in Equation (131(131) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\mathsf{X}}^f\left(\mathsf{T}\mathrm{\Lambda }+[\mathsf{w},\dots ,\mathsf{w}]\right).\end{array}$$(131) ) is applied to the ensemble matrix ${\mathsf{X}}^f$ instead of the ensemble perturbation matrix ${\mathsf{X}}^{{}^{\prime }f}$ without subsequent addition of the analysis mean state (see e.g. Equation 19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ). This is possible because of the property of $\mathsf{T}$ to implicitly subtract the ensemble mean, while the multiplication of ${\mathsf{X}}^f$ with the weight vector array adds the analysis mean state.

For high-dimensional systems, a localisation of the analysis step is required. It was introduced by Tödter and Ahrens (2015) in analogy to the localisation of the ETKF and SEIK filters (see Sec. 9.1). Algorithm  in Appendix 2 gives a pseudo-algorithm of the PFGR and NETF and Algorithm shows the computation of the weights for Gaussian observation errors.

7.2. Moment-Matching Ensemble Filter (MMEF)

A stochastic algorithm that has second-order correct statistics was developed by Lei and Bickel (2011). In this moment-matching ensemble filter (MMEF) we generate an ensemble of perturbed pseudo-observations, ${\mathsf{Y}}^f$, as in the SEnKF (see Equations 20(20) $$\begin{array}{c}\hfill {\mathsf{Y}}^f=\left(\mathcal{H}({\mathsf{x}}_1^f),\mathcal{H}({\mathsf{x}}_2^f),\dots ,\mathcal{H}({\mathsf{x}}_{N_e}^f)\right)+{\mathsf{Y}}^{\prime }\in {\mathcal{R}}^{N_y\times N_e},\end{array}$$(20) and 21(21) $$\begin{array}{c}\hfill {\mathsf{Y}}^{\prime }=\left({ϵ}_1,{ϵ}_2,\dots ,{ϵ}_{N_e}\right)\in {\mathcal{R}}^{N_y\times N_e}\end{array}$$(21) )(132) $$\begin{array}{c}\hfill {\mathsf{Y}}_j^f\sim p_{\mathcal{H}({\mathsf{x}}_{\mathcal{j}})}\left(\mathsf{y}|{\mathsf{x}}_j\right)\end{array}$$(132)

using $\mathcal{H}\left({\mathsf{x}}_j\right)$ as variable in the density, so $\mathsf{y}$ is fixed. Then the analysis mean for each particle is generated using a corresponding pseudo-observation as follows(133) $$\begin{array}{c}\hfill {\overline{\mathsf{x}}}^a\left({\mathsf{y}}_j^f\right)=\sum _{k=1}^{N_e}{w}_k\left({\mathsf{y}}_j^f\right){\mathsf{x}}_k={\mathsf{X}}^f\mathbf{w}\left({\mathsf{y}}_j^f\right)\end{array}$$(133)

in which $w_k\left({\mathsf{y}}_j^f\right)$ is given by(134) $$\begin{array}{c}\hfill w_k\left({\mathsf{y}}_j^f\right)=\frac{p\left({\mathsf{y}}_j^f|{\mathsf{x}}_k\right)}{{\sum }_{l=1}^{N_e}p\left({\mathsf{y}}_j^f|{\mathsf{x}}_l\right)}.\end{array}$$(134)

Similarly, the analysis mean for actual observations is computed via(135) $$\begin{array}{c}\hfill {\overline{\mathsf{x}}}^a\left(\mathsf{y}\right)=\sum _{k=1}^{N_e}{w}_k\left(\mathsf{y}\right){\mathsf{x}}_k={\mathsf{X}}^f\mathbf{w}\left(\mathsf{y}\right).\end{array}$$(135)

Furthermore, we calculate equivalent expressions for covariances for perturbed and actual observations as follows(136) $$\begin{array}{cc}\hfill {\mathsf{P}}^a\left({\mathsf{y}}_j^f\right)& =\sum _{k=1}^{N_e}{w}_k\left({\mathsf{y}}_j^f\right)\left({\mathsf{x}}_k-{\overline{\mathsf{x}}}^a\left({\mathsf{y}}_j^f\right)\right)\hfill \\ \hfill & \times {\left({\mathsf{x}}_k-{\overline{\mathsf{x}}}^a\left({\mathsf{y}}_j^f\right)\right)}^T\hfill \end{array}$$(136) (137) $$\begin{array}{cc}\hfill {\mathsf{P}}^a\left(\mathsf{y}\right)& =\sum _{k=1}^{N_e}{w}_k\left(\mathsf{y}\right)\left({\mathsf{x}}_k-{\overline{\mathsf{x}}}^a\left(\mathsf{y}\right)\right)\hfill \\ \hfill & \times {\left({\mathsf{x}}_k-{\overline{\mathsf{x}}}^a\left(\mathsf{y}\right)\right)}^T.\hfill \end{array}$$(137)

Then each of the ensemble members or particles is updated via(138) $$\begin{array}{c}\hfill {\mathsf{x}}_j^a={\overline{\mathsf{x}}}^a\left(\mathsf{y}\right)+{\mathsf{P}}^a{\left(\mathsf{y}\right)}^{1/2}{\mathsf{P}}^a{\left({\mathsf{y}}_j^f\right)}^{-1/2}\left({\mathsf{x}}_j-{\overline{\mathsf{x}}}^a\left({\mathsf{y}}_j^f\right)\right).\end{array}$$(138)

This filter gives the correct posterior mean and covariance in the large-ensemble limit (Lei and Bickel, 2011). To see this, note that ${\mathsf{x}}_j-{\overline{\mathsf{x}}}^a\left({\mathsf{y}}_j^f\right)$ is distributed according to $N\left(0,{\mathsf{P}}^a\left({\mathsf{y}}_j^f\right)\right)$, so ${\mathsf{P}}^a{\left({\mathsf{y}}_j^f\right)}^{-1/2}\left({\mathsf{x}}_j-{\overline{\mathsf{x}}}^a\left({\mathsf{y}}_j^f\right)\right)$ is distributed $\mathcal{N}\left(0,\mathsf{I}\right)$, and hence the distribution of ${\mathsf{x}}_j^a$ is $\mathcal{N}\left({\overline{\mathsf{x}}}^a\left(\mathsf{y}\right),{\mathsf{P}}^a\left(\mathsf{y}\right)\right)$.

This filter cannot be used in high-dimensional systems, even when localisation is applied, because it needs the evaluation of several full covariance matrices. However, we can explore ensemble perturbations that are used to calculate these covariances, as in all ensemble Kalman filter schemes. The following was not discussed by Lei and Bickel (2011), but is a practical way to make the filter useful in high-dimensional systems.

We can express each covariance matrix ${\mathsf{P}}^a\left({\mathsf{y}}_j^f\right)$ directly in terms of the forecast ensemble as(139) $$\begin{array}{c}\hfill {\mathsf{P}}^a\left({\mathsf{y}}_j^f\right)={\mathsf{X}}^f\mathsf{T}\left({\mathsf{y}}_j^f\right){{\mathsf{X}}^f}^T\end{array}$$(139)

where matrix $\mathsf{T}\left({\mathsf{y}}_j^f\right)$ is given by(140) $$\begin{array}{c}\hfill \mathsf{T}\left({\mathsf{y}}_j^f\right)=\text{diag}\left(\mathbf{w}({\mathsf{y}}_j^f)\right)-\mathbf{w}\left({\mathsf{y}}_j^f\right){\mathbf{w}}^T\left({\mathsf{y}}_j^f\right).\end{array}$$(140)

The square root of this matrix is(141) $$\begin{array}{c}\hfill {\mathsf{P}}^a{\left({\mathsf{y}}_j^f\right)}^{1/2}={\mathsf{X}}^f\mathsf{T}{\left({\mathsf{y}}_j^f\right)}^{1/2}.\end{array}$$(141)

To find the inverse of this matrix we perform an SVD on the prior ensemble matrix(142) $$\begin{array}{c}\hfill {{\mathsf{X}}^{\prime }}^f=\mathsf{U}\mathrm{\Lambda }{\mathsf{V}}^T\end{array}$$(142)

and compute also the EVD on the much smaller square matrices $\mathsf{T}\left({\mathsf{y}}_j^f\right)$(143) $$\begin{array}{c}\hfill \mathsf{T}\left({\mathsf{y}}_j^f\right)={\tilde{\mathsf{U}}}_j{\tilde{\mathrm{\Lambda }}}_j{\tilde{\mathsf{U}}}_j^T.\end{array}$$(143)

Using Equations (142(142) $$\begin{array}{c}\hfill {{\mathsf{X}}^{\prime }}^f=\mathsf{U}\mathrm{\Lambda }{\mathsf{V}}^T\end{array}$$(142) ) and (143(143) $$\begin{array}{c}\hfill \mathsf{T}\left({\mathsf{y}}_j^f\right)={\tilde{\mathsf{U}}}_j{\tilde{\mathrm{\Lambda }}}_j{\tilde{\mathsf{U}}}_j^T.\end{array}$$(143) ) we find(144) $$\begin{array}{c}\hfill {\mathsf{P}}^a{\left({\mathsf{y}}_j^f\right)}^{-\frac{1}{2}}={\tilde{\mathsf{U}}}_j^T{\tilde{\mathrm{\Lambda }}}_j^{1/2}\mathsf{V}{\mathrm{\Lambda }}^{-1}{\mathsf{U}}^T.\end{array}$$(144)

Hence, we can write the update equation of the MMEF as(145) $$\begin{array}{c}\hfill {\mathsf{x}}_j^a={\overline{\mathsf{x}}}^a+{{\mathsf{X}}^{\prime }}^f{\mathsf{T}}^{1/2}{\tilde{\mathsf{U}}}_j^T{\tilde{\mathrm{\Lambda }}}_j^{1/2}\mathsf{V}{\mathrm{\Lambda }}^{-1}{\mathsf{U}}^T\left({\mathsf{x}}_j-{\overline{\mathsf{x}}}^a\left({\mathsf{y}}_j^f\right)\right).\end{array}$$(145)

This expression is suitable for high-dimensional applications when the matrices $\mathsf{T}({\mathsf{y}}_j^f)$ are computed with localisation.

7.3. Merging Particle Filter (MPF)

The merging particle filter generates several sets of posterior ensembles and merges them via a weighted average to obtain a new set of particles that has the correct mean and covariance but is more robust than the standard particle filter. Specifically, the method draws a set of q ensembles each of size $N_e$ from the weighted prior ensemble at the resampling step. Denote each ensemble member as ${\mathsf{x}}_{j,i}$ for ensemble member j in ensemble i. Then new merged ensemble members are generated via(146) $$\begin{array}{c}\hfill {\mathsf{x}}_j^a=\sum _{i=1}^q{\alpha }_i{\mathsf{x}}_{j,i}.\end{array}$$(146)

To ensure that the new ensemble has the correct mean and covariance, the coefficients ${\alpha }_j$ need to fulfil the two conditions(147) $$\begin{array}{c}\hfill \sum _{j=1}^q{\alpha }_j=1;\sum _{j=1}^q{\alpha }_j^2=1,\end{array}$$(147)

where each ${\alpha }_j$ also has to be a real number.

When $q>3$ there is no unique solution for the $\alpha$’s, while for $q=3$ we find:(148) $$\begin{array}{cc}\hfill {\alpha }_1& =\frac{3}{4}\hfill \\ \hfill {\alpha }_2& =\frac{\sqrt{13}+1}{8}\hfill \\ \hfill {\alpha }_3& =-\frac{\sqrt{13}-1}{8}.\hfill \end{array}$$(148)

Although not discussed by Nakano et al. (2007) this scheme will be degenerate for high-dimensional problems. However, we can make the $\alpha$’s space-dependent when $q>3$ and then apply localisation.

8. Adaptive Gaussian mixture filter

Both ensemble Kalman and Monte Carlo-based techniques discussed in Sections 5 and 6, respectively, have their drawbacks. The Gaussian mixture filter (Anderson and Anderson, 1999; Bengtsson et al., 2003; Hoteit et al., 2008) attempts to avoid these by approximating an arbitrary form of the prior by combining multiple Gaussian priors. This gives it the advantage that both the local Kalman filter type correction step as well as the weighting and resampling step of a particle filter can be applied. This possibility also makes it applicable to highly non-linear and high dimensional systems. In this paper, we discuss the adaptive Gaussian mixture filter developed by Stordal et al. (2011) as a representative scheme, out of all Gaussian mixture filters that have been proposed.

In the Gaussian mixture filter, the prior distribution is approximated by a mixture density (Silverman, 1986) where each ensemble member forms the centre of a Gaussian density function(149) $$\begin{array}{c}\hfill p({\mathsf{x}}^f)=\sum _{j=1}^{N_e}\frac{1}{N_e}\mathcal{N}\left({\mathsf{x}}_j^f,{\tilde{\mathsf{P}}}^f\right),\end{array}$$(149)

where $\mathcal{N}({\mathsf{x}}_j,\tilde{\mathsf{P}})$ denotes a multivariate Gaussian kernel density with ensemble member ${\mathsf{x}}_j$ as mean and covariance matrix ${\tilde{\mathsf{P}}}^f=h^2{\mathsf{P}}^f$, in which ${\mathsf{P}}^f$ is the covariance of the whole forecast ensemble and h is a bandwidth parameter. Stordal et al. (2011) discuss that the optimal choice of the bandwidth h is $h_{\mathit{opt}}\sim N_e^{-1/5}$ if we are only interested in the marginal properties of the individual components of $\mathsf{x}$, but that it might be beneficial to choose $h>h_{\mathit{opt}}$ to reduce the risk of filter divergence, since the choice of the bandwidth parameter determines the magnitude of the Kalman filter update step. Thus, the parameter h is treated as the design parameter and is defined by the user. Note that each particle represents the mean of a Gaussian kernel and that the uncertainty associated with each particle is given by the covariance of that Gaussian kernel (Stordal et al., 2011).

If the likelihood is Gaussian, the posterior pdf is again a Gaussian mixture, now with pdf(150) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^a|\mathsf{y}\right)=\sum _{j=1}^{N_e}{w}_j\mathcal{N}\left({\mathsf{x}}_j^a,{\tilde{\mathsf{P}}}^a\right).\end{array}$$(150)

Here, the weights $w_j$ are propotional to $\mathcal{N}\left(\mathsf{y}-\mathsf{H}{\mathsf{x}}_j^f,{\mathsf{R}}^a\right)$ with ${\sum }_j{w}_j=1$ and ${\mathsf{R}}^a=\mathsf{H}{\tilde{\mathsf{P}}}^f{\mathsf{H}}^T+\mathsf{R}$. So, compared to the particle filter the covariance used in the weights is inflated with a term $\mathsf{H}{\tilde{\mathsf{P}}}^f{\mathsf{H}}^T$, leading to more equal weights. Each mean ${\mathsf{x}}_j^a$ and the covariance matrix ${\tilde{\mathsf{P}}}^a$ are obtained using one of the EnKF variants.

In high-dimensional systems, the covariance matrices are never formed explicitly, and the algorithm in Stordal et al. (2011) cannot be used. Hoteit et al. (2008) used an update based in the SEIK filter (see Section 5.2). For a more modern formulation, we provide here an algorithm based on Stordal et al. (2011) but explore an ETKF to avoid the explicit computation of $\tilde{\mathsf{P}}$. First, the matrix(151) $$\begin{array}{c}\hfill {\mathsf{T}}_{\mathit{GM}}{\mathsf{T}}_{\mathit{GM}}^{\mathsf{T}}={\left(\mathsf{I}+\frac{h^2}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}\end{array}$$(151)

is generated with $\mathsf{S}=\mathsf{H}{{\mathsf{X}}^{\prime }}^f$ similar to Equation (44(44) $$\begin{array}{c}\hfill {\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}={\left(\mathsf{I}+\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}.\end{array}$$(44) ), but including the factor $h^2$. Then, we perform an EVD of the symmetric matrix ${({\mathsf{T}}_{\mathit{GM}}{\mathsf{T}}_{\mathit{GM}}^{\mathsf{T}})}^{-1}={\mathsf{U}}_{\mathit{GM}}{\mathrm{\Sigma }}_{\mathit{GM}}{\mathsf{U}}_{\mathit{GM}}^{\mathsf{T}}$ and obtain the symmetric square root(152) $$\begin{array}{c}\hfill {\mathsf{T}}_{\mathit{GM}}={\mathsf{U}}_{\mathit{GM}}{\mathrm{\Sigma }}_{\mathit{GM}}^{-1/2}{\mathsf{U}}_{\mathit{GM}}^{\mathsf{T}}.\end{array}$$(152)

This is used to update the mean of each Gaussian kernel by calculating the ETKF update on each of the prior particles as(153) $$\begin{array}{c}\hfill {\overline{\mathsf{w}}}_{j,GM}=\frac{1}{\sqrt{(N_e-1)}}{\mathsf{U}}_{\mathit{GM}}{\mathrm{\Sigma }}_{\mathit{GM}}^{-1}{\mathsf{U}}_{\mathit{GM}}^{\mathsf{T}}{(h{{\mathsf{X}}^{\prime }}^f)}^{\mathsf{T}}{\mathsf{H}}^{\mathsf{T}}{\mathsf{R}}^{-1}{\mathsf{d}}_j\end{array}$$(153)

in which ${\mathsf{d}}_j=\mathsf{y}-\mathsf{H}{\mathsf{x}}_j^f$. The new centres of the Gaussian mixture densities are now found as(154) $$\begin{array}{c}\hfill {\mathsf{x}}_j^a={\mathsf{x}}_j^f+h{{\mathsf{X}}^{\prime }}^f{\overline{\mathsf{w}}}_{j,GM}.\end{array}$$(154)

A square root of the posterior covariance of each Gaussian mixture density is found by(155) $$\begin{array}{c}\hfill \mathsf{Z}=h{{\mathsf{X}}^{\prime }}^f{{\mathsf{W}}^{\prime }}_{\mathit{GM}}\end{array}$$(155)

in which(156) $$\begin{array}{c}\hfill {{\mathsf{W}}^{\prime }}_{\mathit{GM}}={\mathsf{U}}_{\mathit{GM}}{\mathrm{\Sigma }}_{\mathit{GM}}^{-1/2}{\mathsf{U}}_{\mathit{GM}}^{\mathsf{T}}.\end{array}$$(156)

Thus, for Equation (150(150) $$\begin{array}{c}\hfill p\left({\mathsf{x}}^a|\mathsf{y}\right)=\sum _{j=1}^{N_e}{w}_j\mathcal{N}\left({\mathsf{x}}_j^a,{\tilde{\mathsf{P}}}^a\right).\end{array}$$(150) ) we have ${\tilde{\mathsf{P}}}^a={(N_e-1)}^{-1}\mathsf{Z}{\mathsf{Z}}^T$, but to evaluate the equation one can use the square root $\mathsf{Z}$, so it is not required to compute ${\tilde{\mathsf{P}}}^a$ explicitly.

Until this point, the algorithm is the standard Gaussian mixture filter. The adaptive part of the filter was introduced by Stordal et al. (2011) and has been demonstrated to avoid filter divergence due to ensemble degeneration. Further, the adaptivity allows us to choose smaller values of the bandwidth parameter h. To stabilise the Gaussian mixture filter, we interpolate the original analysis weights with a uniform weight as(157) $$\begin{array}{c}\hfill w_j^{\alpha }=\alpha w_j+(1-\alpha )N_e^{-1}.\end{array}$$(157)

For the adaptivity, $\alpha$ is chosen to be(158) $$\begin{array}{c}\hfill \alpha =N_{\text{eff}}{N}_e^{-1},\end{array}$$(158)

where $N_{\text{eff}}=1/({\sum }_{l=1}^{N_e}{w}_l^2)$ is the effective ensemble size. To avoid ensemble degeneration one can further add a resampling step as in particle filters. It is performed if $N_{\text{eff}}<N_c$, with $N_c$ a value that can be chosen freely, for instance $N_c=0.5N_e$. The full scheme then becomes:

(a)

When $N_{\text{eff}}\ge N_c$ no resampling is needed, so the weights are calculated as above and transported with each particle to the next set of observations.

(b)

When $N_{\text{eff}}<N_c$ we will resample according to any of the resampling schemes in Appendix 1. This leads to a new set of states for the centres of the Gaussian mixtures denoted ${\mathsf{x}}_{j,(i)}$ in which j denotes the index of the state for resampling, and i its index after resampling. Note that several of the new states will coincide. To avoid identical samples we draw our final new ensemble from the Gaussian mixtures, as follows (159) $$\begin{array}{ccc}\hfill {\xi }_i& \sim \hfill & \hfill \mathcal{N}(0,\mathsf{I}),\end{array}$$(159) (160) $$\begin{array}{cc}\hfill {\mathsf{x}}_i^a& ={\mathsf{x}}_{j,(i)}+\mathsf{Z}{\xi }_i,\hfill \end{array}$$(160) (161) $$\begin{array}{cc}\hfill w_i& =\frac{1}{N_e}.\hfill \end{array}$$(161) Further, we set ${\tilde{\mathsf{P}}}^a={(N_e-1)}^{-1}{\mathsf{X}}^a{({\mathsf{X}}^a)}^T$, but use this only in factorised form.

Note that in this scheme we never calculate a full state covariance.

It is important to realise what the adaptive part does. Indeed, by construction, the filter is not degenerate, but at the expense of strongly reducing the influence of the observations when $\alpha$ is small. In high dimensional systems with a large number of independent observations localisation is essential to avoid using the scheme as a sum of ensemble Kalman filters only.

In the scheme by Bengtsson et al. (2003), the mean of each Gaussian pdf is chosen at random from the ensemble, and the covariance in each Gaussian pdf is estimated from the ensemble members which are local in state space, including a localisation and smoothing step. Since the scheme has not been applied to high-dimensional systems it will not be discussed here.

9. Practical implementation of the ensemble methods

This section is devoted to issues related to the practical implementation of the ensemble methods. In particular, we address the need for localisation and inflation in some of the presented ensemble methods to counteract the issues arising from ensemble undersampling in large scale problems such as ocean and atmosphere prediction. We also discuss the computational cost of each method as presented and the parallelisation of ensemble data assimilation methods. We will conclude this section with a discussion on the suitability of the ensemble data assimilation methods applied to non-linear dynamical models.

9.1. Localisation in EnDA

The success of the EnDA methods is highly dependent on the size of the ensemble being adequate for the system we apply these methods to. Thus, for large scale problems, where the number of state variables is many magnitudes larger than the number of ensemble members, ensemble undersampling can cause major problems in EnDA methods: underestimated ensemble variance, filter divergence, and errors in estimated correlations, in particular spurious long-range correlations. In such cases, spatial localisation is a necessary tool to minimise the effect of undersampling.

Localisation damps long-range correlations, e.g. in the ensemble covariance matrix (’covariance localisation’, see Section 9.1.2). This damping can be applied to the extend to keep only correlations over limited distances and erase long-range correlations in the analysis step. Thus, localisation decouples the analysis update at distant locations in a model grid. The underlying assumption of localisation is that the assimilation problem has in fact a local structure. This means, that correlation length scales are much shorter than the extent of the model grid so that only correlations over short distances are relevant while for long distances the sampling error in the ensemble-estimated covariance matrices dominates (see, e.g. Morzfeld et al., 2017). This seems to be fulfilled for many oceanic and atmospheric applications. For example, Patil et al. (2001) described a locally low dimension for atmospheric dynamics. The success of localised filters in oceanic and atmospheric data assimilation applications also shows that this condition is dominantly fulfilled, even though it is known that long-range correlations (teleconnections) exist in the atmosphere and ocean. However, if a modelling problem does not have a local structure or if too little observations are available or the observations only represent long-range properties of the system, localisation cannot be applied.

Localisation is usually applied either explicitly by considering only observations from a region surrounding the location of the analysis or implicitly by modifying $\mathsf{P}$ or $\mathsf{R}$ so that observations from beyond a certain distance cannot affect the analysis state. The way in which such localisation is applied is still an active field of research and many variants of localisation schemes have emerged over the last decade. There are two main types of spatial localisation techniques (or simply localisation) that are widely used in ensemble data assimilation: Covariance localisation (also termed P- or B-localisation) and observation localisation (also denoted R-localisation). Both methods will be discussed here together with domain localisation, which is required for the application of observation localisation. In addition, a number of adaptive localisation schemes was developed over the recent years. A selection of these schemes is discussed in Section 9.1.4.

In general, all localisation schemes are empirical. While they improve the estimations by ensemble filters, they can disturb balances in the model state (Lorenc, 2003; Kepert, 2009). Further, the interaction of localisation with the serial observation processing usually applied with the EnSRF and EAKF methods can reduce the stability of these filters (Nerger, 2015).

9.1.1. Domain localisation

Domain localisation or local analysis is the oldest localisation technique. For ensemble Kalman filters it was first applied by Houtekamer and Mitchell (1998), but the method was also applied in earlier schemes of optimal interpolation (see Cohn et al., 1998). In domain localisation we only use the ensemble perturbations that belong to the domain $D_{\gamma }$ in which the analysis correction of the state estimate is computed. For example, this domain can be a vertical column of grid points or a single grid point. Thus, we use a linear transformation ${\mathsf{D}}_{\gamma }$ to obtain(162) $$\begin{array}{c}\hfill {{\mathsf{x}}^{\prime }}_{j,\gamma }={\mathsf{D}}_{\gamma }{{\mathsf{x}}^{\prime }}_j^f,\end{array}$$(162)

where $j=1,\dots ,N_e$ and $\gamma =1,\dots ,\mathrm{\Gamma }$ with $\mathrm{\Gamma }$ being the total number of subdomains. To localise, we now only use observations within a specified distance – the localization radius – around the local domain $D_{\gamma }$. This defines a local observation domain ${\widehat{D}}_{\gamma }$. Using the corresponding linear transformation ${\widehat{\mathsf{D}}}_{\gamma }$ we can transform the observation error covariance $\mathsf{R}$, the global observation vector $\mathsf{y}$, and the global observation operator $\mathsf{H}$ analogously to Equation (162(162) $$\begin{array}{c}\hfill {{\mathsf{x}}^{\prime }}_{j,\gamma }={\mathsf{D}}_{\gamma }{{\mathsf{x}}^{\prime }}_j^f,\end{array}$$(162) ) to their local parts(163) $$\begin{array}{cc}\hfill {\mathsf{y}}_{\gamma }& ={\widehat{\mathsf{D}}}_{\gamma }\mathsf{y},\hfill \end{array}$$(163) (164) $$\begin{array}{cc}\hfill {\mathsf{H}}_{\gamma }& ={\widehat{\mathsf{D}}}_{\gamma }\mathsf{H},\hfill \end{array}$$(164) (165) $$\begin{array}{cc}\hfill {\mathsf{R}}_{\gamma }& ={\widehat{\mathsf{D}}}_{\gamma }\mathsf{R}{\mathsf{D}}_{\gamma }^{\mathsf{T}}.\hfill \end{array}$$(165)

Thus, we neglect observations that are outside of the domain ${\widehat{D}}_{\gamma }$. Then a general local analysis state is given by(166) $$\begin{array}{c}\hfill {\mathsf{X}}_{\gamma }^a={\overline{\mathsf{X}}}_{\gamma }^f+{{\mathsf{X}}^{\prime }}_{\gamma }^f\left({\overline{\mathsf{W}}}_{\gamma }+{{\mathsf{W}}^{\prime }}_{\gamma }\right),\end{array}$$(166)

where ${\overline{\mathsf{W}}}_{\gamma }$ and ${{\mathsf{W}}^{\prime }}_{\gamma }$ are computed using local ensemble forecast perturbations and local observations from the domain ${\widehat{D}}_{\gamma }$. For a complete analysis update, a loop over all local analysis domains has to be performed with a local analysis update for each domain.

Applying domain localisation allows significant savings in computing time since solving for the analysis update is not performed globally but on much smaller local domains. Accordingly, updates on the smaller scale domains can be done independently and therefore parallel (Nerger et al., 2006) even if the observation domains overlap. In ensemble-based Kalman filters, domain localisation was used predominantly with filters that use the analysis error covariance matrix for the calculation of the gain like SEIK, ETKF, ESTKF, all discussed in detail in Section 5. In these algorithms, the forecast error covariance matrix is never explicitly computed. Examples of the application of domain localisation can be found, e.g. in Brusdal et al. (2003) and Testut et al. (2003).

Blindly using domain localisation can result in boxed analysis fields if neighbouring local domains are updated using significantly different observation sets. Thus, great care needs to be taken to choose domains so that they overlap sufficiently to produce smooth global analysis fields with minimal increase in computational cost. Today, domain localisation is typically applied with observation localisation (Hunt et al., 2007), which is discussed in Section 9.1.3.

9.1.2. Covariance localisation

Covariance localisation (also termed P-localisation or B-localisation, depending on whether the background covariance matrix is denoted $\mathsf{P}$ or $\mathsf{B}$ as in variational assimilation schemes) is a localisation method that is directly applied to the ensemble covariance matrix. The ensemble undersampling causes spurious cross-correlations between state variables. As realistic long-range correlations are typically small, the sampling errors are particularly pronounced for long distances. The direct covariance localisation can be used to reduce the long-range correlations in the forecast error covariance and hence damp the spurious correlations. In addition, the rank of the ensemble covariance is increased, giving more degrees of freedom to the analysis update (Hamill et al., 2001; Whitaker and Hamill, 2002).

Typically, covariance localisation is applied by first forming a correlation matrix $\mathsf{C}$ and then taking a Schur product (an element by element matrix multiplication) of this correlation matrix and the forecast error covariance. Thus, given some ${\mathsf{P}}^f$, our localised forecast error covariance will be(167) $$\begin{array}{c}\hfill {\mathsf{P}}_L^f=\mathsf{C}\circ {\mathsf{P}}^f.\end{array}$$(167)

The localization matrix $\mathsf{C}$ is usually formed of correlation functions with compact support similar in shape to a Gaussian function [e.g.][]Gaspari1999. Practically, the computation of the covariance matrix ${\mathsf{P}}^f$ can be avoided by applying the localisation matrix to the matrices ${\mathsf{P}}^f{\mathsf{H}}^T$ and $\mathsf{H}{\mathsf{P}}^f{\mathsf{H}}^T$ (see Equation (11(11) $$\begin{array}{c}\hfill \mathsf{K}={\mathsf{P}}^f{\mathsf{H}}^{\mathsf{T}}{(\mathsf{H}{\mathsf{P}}^f{\mathsf{H}}^{\mathsf{T}}+\mathsf{R})}^{-1}.\end{array}$$(11) )).

We note that, from all the ensemble-based Kalman filter methods presented in Section 5, covariance localisation can only be applied to the EnSRF and EAKF, since for these methods observations can be processed serially, and in the stochastic EnKF.

9.1.3. Observation localisation

In the case of square root filters, presented in Section 5, the full covariance matrix is never formed. Instead, only the ensemble perturbation matrix ${{\mathbf{X}}^{\prime }}^f$ is calculated at each analysis step. Petrie and Dance (2010) showed that covariance-localisation for square root filters cannot be approximately decomposed into a square root of the correlation matrix $\rho$,(168) $$\begin{array}{c}\hfill \left(\rho {\rho }^{\mathsf{T}}\right)\circ \left({\mathsf{X}}^{\prime }{\mathsf{X}}^{\prime \mathsf{T}}\right)\ne \left(\rho \circ {\mathsf{X}}^{\prime }\right){\left(\rho \circ {\mathsf{X}}^{\prime }\right)}^{\mathsf{T}}\end{array}$$(168)

thus covariance localisation cannot be applied. For such filters, e.g. SEIK, ETKF and ESTKF, the observation localisation is a more natural choice and is currently used instead of covariance localisation (Hunt et al., 2007; Miyoshi and Yamane, 2007; Janjić et al., 2011).

Observation localisation is applied by modifying the observation error covariance matrix $\mathsf{R}$. More specifically, one modifies its inverse ${\mathsf{R}}^{-1}$ so that the inverse observation variance decreases to zero with the distance of an observation from an analysis grid point. To be able to define the distance, it is necessary to perform the analysis with the domain localization method as described in Section 9.1.1. An abrupt cutoff could be obtained by setting observation variances to zero beyond a given distance. This would be equivalent to the simple domain localisation of Section 9.1.1 and could result in non-smooth analysis updates. For a smooth analysis, e.g. Brankart et al. (2003) described to increase the observation error variance with increasing distance from the analysis grid point. Hunt et al. (2007) proposed to use a gradual observation localisation in the LETKF acting on ${\mathsf{R}}^{-1}$, which is likewise applicable with the SEIK filter and the ESTKF. In this case, elements of ${\mathsf{R}}^{-1}$ are multiplied by a smoothly decreasing function of distance from the analysis grid point. This modification smoothly reduces the observation influence and excludes observations outside a defined radius by prescribing their error to be infinitely large. As for covariance-localisation, the method uses a Schur product as(169) $$\begin{array}{c}\hfill \tilde{\mathsf{R}}=\tilde{\mathsf{C}}\circ \mathsf{R}.\end{array}$$(169)

Here, the same correlation function (Gaspari and Cohn, 1999) as for covariance localisation can be used to construct the localisation matrix. However, in contrast to covariance localisation, $\tilde{\mathsf{C}}$ is not a correlation matrix as the values on the diagonal of this matrix vary with the distance between the observation and the local analysis domain. Then, the analysis update is computed as in the case of domain localisation, but using the weight-localised matrix $\tilde{\mathsf{R}}$. For computational savings we would in practise also discard any observations with zero weight from the analysis computations.

Both observation and covariance localisation can lead to similar assimilation results. In general, the optimal localisation has been found to be a bit larger for covariance localisation than for observation localisation (Greybush et al., 2011). The reason for this difference lies in the different effect of the localisations in the Kalman gain as was explained by Nerger et al. (2012b).

9.1.4. Adaptive localisation schemes

The localisation methods described above are widely used and can be applied without much additional computing cost. However, the optimal localisation radius is a priori unknown and needs to be tuned in numerical experiments. For the tuning one performs several data assimilation experiments with different localisation radii, perhaps over shorter time periods, and selects the radius that results in the smallest estimation errors. Regarding the theoretical understanding of localisation, Kirchgessner et al. (2014) showed for the case of observation localisation when each grid point is observed that the optimal localisation radius should be reached when the sum over the observation weights equals the ensemble size. This finding allows for a simple form of adaptivity or a starting point for further tuning. Further, Perianez et al. (2014) showed that both the sampling error in the ensemble covariance matrix and the observation error influence the optimal localisation radius. As the sampling error has a largest influence when the true correlations are small, the dynamically generated correlations also influence the optimal localisation radius (Zhen and Zhang, 2014; Flowerdew, 2015).

To avoid the need for numerical tuning and to better adapt the localisation to the dynamically created correlation structure, several adaptive localisation methods have been developed, which we shortly mention here. A common approach is to damp the spurious correlations that are caused by sampling errors due to the small ensemble size. Anderson (2007) developed a hierarchical localisation method, in which the ensemble is partitioned into sub-ensembles. Then, the sub-ensembles are used to estimate the sampling errors. Bishop and Hodyss (2009) proposed an adaptive localisation method that uses a power of the correlations to damp small correlations and pronounce those correlations that are significant. This method can find correlations even at longer distances. Further, methods have been developed to find empirical localisation functions. In these methods, one attempts to find for a single observation the weight factor that minimises the deviation from a true solution (Anderson, 2012; Lei and Anderson, 2014; Flowerdew, 2015). These methods are typically tuned once based on observation system simulation experiments (OSSEs), in which one knows the true state. When the OSSEs are configured realistically, the obtained localisation functions should be applicable for the assimilation of real observations after the tuning.

The major advantage of the methods proposed so far is that they are able to adaptively specify the localisation function or radius according to the dynamically generated covariance structure. However, the methods still need tuning, which can be even more costly than for the fixed covariance and observation localisation methods. For example, the method by Bishop and Hodyss (2009) requires the specification of an envelope function around the locations that are found by powering the correlations and the number of powers that are computed. Lei and Anderson (2014) also showed that the localisation function can change when it is applied iteratively such that a sufficient number of iterations have to be computed.

Apart from the adaptive localisation methods, further methods like spectral localisation (Buehner and Charron, 2007) and localisation in different variables (i.e. stream function, velocity potential, Kepert, 2006) have been developed. However, none of these methods are yet a standard for operational centres.

9.1.5. Localisation in particle filtering

Several variants of the Particle Filter that explore localisation have been developed recently, following its success in Ensemble Kalman Filters. An issue with directly localising $\mathsf{R}$ or using domain localisation is that the weight of each particle is a global property of the filter (van Leeuwen, 2009). That is, the same particle could have a high weight in one area and low weight in another making it ambiguous whether this particle should be resampled or not. Keeping parts of a number of particles that all perform well in a certain area of the domain and parts of other particles in other areas of the domain would lead to balance problems between variables and sharp gradients in the fields. In contrast, when performing parameter estimation a smooth variation of parameter values is less likely to cause imbalances in the model variables, and localisation is straightforward, as pioneered by Vossepoel and Van Leeuwen (2006).

Particle filters that use a proposal density, such as the EWPF discussed in Section 6.2.1 indirectly use localisation through the model error covariance matrix $\mathsf{Q}$. This localisation does not explicitly work on the weights but on how the states are updated, because a natural choice is to pre-multiply each update of a particle with that matrix. Since the model error covariance matrix will mainly contain short length-scale correlations related to missing or inaccurate physics at the model grid scale, each point in the state space is only influenced by observations within the radius set by that covariance matrix. In fact, as noted in Section 6.2.1, we do have the freedom to choose this matrix differently from $\mathsf{Q}$, so other choices closer to our needs are possible. This is because the effects of this choice will be taken into account in the computation of the weights of each particle. This has not been explored in any detail in the literature.

Of the full particle filters, the ETPF (Reich, 2013) can easily be localised by taking for each grid point only observations close to that grid point into account and making the transformation matrix space-dependent to ensure smooth transitions between different regions. This can for example be achieved by calculating the transformation matrix in a limited number of grid points and interpolate that matrix between grid points. This would also reduce the number of computations, which would otherwise be prohibitive (see Section 9.4 on computational costs).

The PFGR and the NETF perform an ensemble transformation similar to the ETKF, but with a transform matrix $\mathsf{T}$ computed from particle filter weights. Accordingly, observation localisation can be applied to the NETF (Tödter and Ahrens, 2015) by smoothing the weight matrix over space. This can also be applied to the MMPF in the high-dimensional implementation. Also the MPF can be localised by making the weights local and using a systematic resampling method like Stochastic Universal Resampling (see Appendix 1). In practise, more might be needed, e.g. the extra averaging as advocated by Penny and Miyoshi (2016) described below.

Several localisation schemes have been proposed and discussed in the review van Leeuwen (2009) and those will not be repeated here. The most obvious thing to do is to weight and resample locally, and somehow glue the resampled particles together via averaging at the edges between resampled local particles (van Leeuwen, 2003b). Recently, Penny and Miyoshi (2016) used this idea with more extensive averaging, and their scheme runs as follows. First, for each grid point j the observations close to that grid point are found and the weight of each particle i is calculated based on the likelihood of only those observations:(170) $$\begin{array}{c}\hfill w_{i,j}=\frac{p({\mathsf{y}}_j|{\mathsf{x}}_{i,j})}{{\sum }_{k=1}^{N_e}p({\mathsf{y}}_j|{\mathsf{x}}_{i,j})}\end{array}$$(170)

in which ${\mathsf{y}}_j$ denotes the set of observations within the localisation area. This is followed by resampling via Stochastic Universal Resampling to obtain ensemble members $x_{i,j}^a$ with $i=1,\dots ,N_e$ for each grid point j. As mentioned before, the issue is that two neighbouring grid points can have different sets of particles, and smoothing is needed to ensure that the posterior ensemble consists of smooth particles. This smoothing is performed for each grid point j for each particle i by averaging over the $N_p$ neighbouring points within the localisation area around grid point j:(171) $$\begin{array}{c}\hfill x_{i,j}^a=\frac{1}{2}{x}_{i,j}^a+\frac{1}{N_p}\sum _{k=1}^{2N_p}{x}_{i,j_k}^a\end{array}$$(171)

in which $j_k$ for $k=1,\dots ,N_p$ denotes the grid point index for those points in the localisation area around grid point j. The resampling via Stochastic Universal Resampling is done such that the weights are sorted before resampling, so that high-weight particles are joined up to reduce spurious gradients.

While this scheme does solve the degeneracy problem in simple one-dimensional systems it is unclear if it will work well in complex systems such as the atmosphere in which fronts can easily be smoothed out, and non-linear balances broken, see e.g. the discussion in van Leeuwen (2009).

A new scheme has recently been proposed in Poterjoy (2016a), which involves a very careful process of ensuring smooth posterior particles and retaining non-linear relations. The filter processes each observation sequentially, as follows. First, adapted weights are calculated for the first element ${\mathsf{y}}_1$ of the observation vector, as(172) $$\begin{array}{c}\hfill {\tilde{w}}_i=\alpha p({\mathsf{y}}_1|{\mathsf{x}}_i)+1-\alpha \end{array}$$(172)

These weights are then normalised by their sum $\tilde{W}$. Then we resample the ensemble according to these normalised weights to form particles ${\mathsf{x}}_{k_i}$.

Here, $\alpha$ is an important parameter in this scheme, with $\alpha =1$ leading to standard weighting, and $\alpha =0$ leading to all weights being equal to 1. Its importance lies in the fact that the weights are always larger than $1-\alpha$, so even a value close to 1, say $\alpha =0.99$, leads to a minimum weight of 0.01 that might seem small, but it means that particles that are more than 1.7 observational standard deviations away from the observations have their weights cut off to something close to $1-\alpha$. This seriously limits the influence the observation can have on the ensemble. Furthermore, the influence of $\alpha$ does depend on the size of the observational error, which is perhaps not what one would like. It is included to avoid losing any particle.

Now, we do the following for each grid point j. For each member i we calculate a weight(173) $$\begin{array}{c}\hfill {\tilde{w}}_i=\alpha \rho ({\mathsf{y}}_1,{\mathsf{x}}_j,r)p({\mathsf{y}}_1|{\mathsf{x}}_i)+1-\alpha \rho ({\mathsf{y}}_1,{\mathsf{x}}_j,r)\end{array}$$(173)

in which $\rho (.)$ is the localisation function with localisation radius r. The normalised weights for this grid point, $w_i$, are obtained by dividing ${\tilde{w}}_i$ by the summed weights over all the particles. Note, again, the role played by $\alpha$. Then, the posterior mean for this observation at this grid point is calculated as(174) $$\begin{array}{c}\hfill {\overline{\mathsf{x}}}_j=\sum _{i=1}^{N_e}{w}_i{\mathsf{x}}_{i,j}\end{array}$$(174)

in which ${\mathsf{x}}_{i,j}$ is grid point j of particle i. Next, a number of scalars are calculated that ensure smooth posterior fields (Poterjoy, 2016a):(175) $$\begin{array}{cc}\hfill {\sigma }_j^2& =\sum _{i=1}^{N_e}{\omega }_i{({\mathsf{x}}_{i,j}-{\overline{\mathsf{x}}}_j)}^2\hfill \\ \hfill c_j& =\frac{N_e(1-\alpha \rho ({\mathsf{x}}_j,{\mathsf{y}}_1,r))}{\alpha \rho ({\mathsf{x}}_j,{\mathsf{y}}_1,r)\tilde{W}}\hfill \\ \hfill r_{1,j}& =\sqrt{\frac{{\sigma }_j^2}{\frac{1}{N_e-1}\sum _{i=1}^{N_e}{({\mathsf{x}}_{k_i,j}-{\overline{\mathsf{x}}}_j+c_j({\mathsf{x}}_{i,j}-{\overline{\mathsf{x}}}_j))}^2}}\hfill \\ \hfill r_{2,j}& =c_j{r}_{1j}\hfill \end{array}$$(175)

so that the final estimate becomes:(176) $$\begin{array}{c}\hfill {\mathsf{x}}_{i,j}^a={\overline{\mathsf{x}}}_j+r_{1,j}({\mathsf{x}}_{k_i,j}-{\overline{\mathsf{x}}}_j)+r_{2,j}({\mathsf{x}}_{i,j}-{\overline{\mathsf{x}}}_j).\end{array}$$(176)

This procedure is followed for each grid point so that at the end we have an updated set of particles that have incorporated the first observation. As a next step the whole process is repeated for the next observation, with a small change that ${\tilde{w}}_i$ is multiplied by ${\tilde{w}}_i$ from the previous observation, until all observations have been assimilated. In this way, the full weight of all observations is accumulated in the algorithm. Now the importance of $\alpha$ comes to full light: without $\alpha$ the ensemble would collapse because the $\tilde{w}$’s would be degenerate when observations are accumulated.

The final estimate shows that each particle at grid point j is the posterior mean at that point plus a contribution from the deviation of the posterior resampled particle from that mean and a contribution from the deviation of the prior particle from that mean. So each particle is a mixture of posterior and prior particles, and departures from the prior are suppressed. When $\alpha =1$, so for a full particle filter, we find for grid points at the observation locations that $c_j=0$ because it is $\rho ({\mathsf{y}}_1,{\mathsf{x}}_j,r)=1$ here. Accordingly, it is $r_{2,j}=0$ and $r_{1,j}\approx 1$ and indeed the scheme gives back the full particle filter.

Between observation locations it can be shown that the particles have the correct first and second order moments, but higher-order moments are not conserved. To remedy this a probabilistic correction is applied at each grid point as follows. The prior particles are dressed by Gaussians with width 1 and weighted by the likelihood weights to generate the correct posterior pdf. The posterior particles are dressed in the same way, each with weight $1/N_e$. Then the cumulative distribution functions (cdf’s) for the two densities are calculated using a trapezoidal rule integration. A cubic spline is used to find the prior cdf values at each prior particle i, denoted by cdf(i). Then a cubic spline is fitted to the other cdf, and the posterior particle i is found as the inverse of its cdf at value cdf(i). See Poterjoy (2016a) for details. The result of this procedure is that higher order moments are brought back into the ensemble between observed points.

This scheme, although rather complicated, is the only local particle filter scheme that has been applied to high-dimensional geophysical systems based on primitive equations in Poterjoy and Anderson (2016b). (van Leeuwen, 2003b applied a local particle filter to a high-dimensional quasi-geostrophic system, but that system is quite robust to sharp gradients as it does not allow gravity waves.)

Another interesting local particle filter is the Multivariate Rank Histogram Filter (Metref et al., 2014a). The idea is to write the posterior pdf in terms of an observed marginal multiplied by a set of conditional pdfs. For example, for a 3-dimensional system in which variable $x_1$ is observed we have:(177) $$\begin{array}{cc}\hfill p(x_1,x_2,x_3|y)& =\frac{p(y|x_1)}{p(y)}p(x_1,x_2,x_3)\hfill \\ \hfill & =\frac{p(y|x_1)}{p(y)}p(x_1,x_2)p(x_3|x_1,x_2)\hfill \\ \hfill & =\frac{p(y|x_1)}{p(y)}p(x_1)p(x_2|x_1)p(x_3|x_1,x_2).\hfill \end{array}$$(177)

The filter now uses the rank-histogram idea of Anderson (2010) on each component, resulting in a fully non-Gaussian update of each component. Localisation can be easily applied directly in this algorithm as it is a transformation algorithm and the transformation can be made local. Unfortunately, this procedure becomes too expensive when the system is high dimensional. However, via a so-called mean-field approximation we suppress the conditioning on non-observed variables, so that we find:(178) $$\begin{array}{c}\hfill p(x_1,x_2,x_3|y)\approx \frac{p(y|x_1)}{p(y)}p(x_1)p(x_2|x_1)p(x_3|x_1).\end{array}$$(178)

This will make the algorithm parallelisable and suitable for high-dimensional applications, although that has not been explored yet.

9.2. Ensemble covariance inflation

In practice, an ensemble Kalman filter can diverge from the truth due to systematic underestimation of the error variances in the filter, possibly caused by model errors or ensemble undersampling as discussed in Section 9.1. In particular, estimating a too large amount of long range correlation will reduce the estimated variance too strongly. Regardless of the cause, underestimating the uncertainty leads to a filter that is overly confident in the state estimate. Thus, the analysis step of the filter puts increasingly more weight on the ensemble background estimate than on the observations and, at some point, it disregards observations completely. Localisation is one method to reduce the undersampling. However, for high-dimensional systems, localisation alone is not sufficient to ensure a stable assimilation process and covariance inflation is applied to further increase the sampled variance and thus stabilise the filter. In addition, the inflation can partly account for model error in case of an imperfect model (Pham et al., 1998b; Hamill, 2001; Anderson, 2001; Whitaker and Hamill, 2002; Hunt et al., 2007).

Most common is a fixed multiplicative covariance inflation (Anderson and Anderson, 1999). The method uses the inflation factor r to perform a multiplicative inflation for each ensemble member ${\mathsf{x}}_j^{a,f}$. With $j=1,\dots ,N_e$ being ensemble member indices, it is given by(179) $$\begin{array}{c}\hfill {\mathsf{x}}_j^{a,f}=r\left({\mathsf{x}}_j^{a,f}-{\overline{\mathsf{x}}}^{a,f}\right)+{\overline{\mathsf{x}}}^{a,f}\end{array}$$(179)

where r usually is chosen to be slightly greater than one. The specification of an optimal inflation factor may vary according to the size of the ensemble ((Hamill, 2001); Whitaker and Hamill, 2002) and the choice of r will depend on various factors, such as dynamics of the model, type of the ensemble filter used as well as the length scale of covariance localisation.

Related to covariance inflation is the so-called ’forgetting-factor’ $\rho$ introduced by Pham et al. (1998b). The forgetting factor is usually chosen to be slightly lower than one and is typically applied in the square root filters like the ETKF, SEIK and ESTKF. For example, in the ETKF it is applied to $\mathsf{T}{\mathsf{T}}^T$, e.g. in Equation (44(44) $$\begin{array}{c}\hfill {\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}={\left(\mathsf{I}+\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}.\end{array}$$(44) ) as(180) $$\begin{array}{c}\hfill {\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}={\left(\rho \mathsf{I}+\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}.\end{array}$$(180)

In this way, the inflation and forgetting factors are related as $\rho =r^{-2}$. Equation (180(180) $$\begin{array}{c}\hfill {\mathsf{T}}_T{\mathsf{T}}_T^{\mathsf{T}}={\left(\rho \mathsf{I}+\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}.\end{array}$$(180) ) allows one to apply inflation in a computationally very efficient way because $\mathsf{T}{\mathsf{T}}^T$ is much smaller than the ensemble states to which the inflation is applied in Equation (179(179) $$\begin{array}{c}\hfill {\mathsf{x}}_j^{a,f}=r\left({\mathsf{x}}_j^{a,f}-{\overline{\mathsf{x}}}^{a,f}\right)+{\overline{\mathsf{x}}}^{a,f}\end{array}$$(179) ).

Next to the multiplicative inflation, an additive inflation has been proposed. The multiplicative inflation leads to an inflation that is relative to the variance level. Thus, large variances will be inflated much more than small variances. This behaviour can be avoided with additive inflation (Ott et al., 2004), which can also be applied in combination with the multiplicative inflation. In additive inflation, all variances are inflated by the same amount, rather than a relative factor. This difference can be useful if the variances vary strongly as in this case the additive inflation acts stronger on the very small variances.

The optimal strength of the inflation is usually determined by tuning experiments, i.e. running experiments with different inflation values and analysing which value results in the smallest estimation errors. Usually a single fixed value of r or $\rho$ is chosen for all grid points. This situation is mainly motivated by the fact that a manual tuning of spatially varying inflations is not feasible for high-dimensional models. To avoid the tuning, several adaptive inflation methods have been proposed. Brankart et al. (2003) proposed to use the relation(181) $$\begin{array}{c}\hfill \mathrm{t}r({\rho }^{-1}\mathsf{S}{\mathsf{S}}^T+\mathsf{R})={(\mathsf{y}-\mathcal{H}({\overline{\mathsf{x}}}^f))}^T(\mathsf{y}-\mathcal{H}({\overline{\mathsf{x}}}^f))\end{array}$$(181)

to estimate a temporally variable forgetting factor $\rho$ for multiplicative inflation. This equation is one of the statistical consistency relations in observation space that Kalman filters should fulfil (Desroziers et al., 2005). Further, Anderson (2009) proposed a method to adaptively estimate spatially and temporally varying inflation factors. This method also aims to fulfil Equation (181(181) $$\begin{array}{c}\hfill \mathrm{t}r({\rho }^{-1}\mathsf{S}{\mathsf{S}}^T+\mathsf{R})={(\mathsf{y}-\mathcal{H}({\overline{\mathsf{x}}}^f))}^T(\mathsf{y}-\mathcal{H}({\overline{\mathsf{x}}}^f))\end{array}$$(181) ) but uses Bayesian estimation to obtain the inflation values. All of these adaptive methods do assume that we have a very good knowledge of the error covariance of the observations. Apart from adaptively inflating the ensemble spread, adaptive inflation of observation errors has been proposed by Minamide and Zhang (2017) for assimilating all-sky satellite brightness temperatures.

An alternative to the inflation can be to explicitly account for the sampling error caused by the finite ensemble size as is done in the finite-size ensemble transform Kalman filter (Bocquet, 2011). This method, while still denoted ’Kalman filter’, requires the iterative minimisation of a cost functional and is hence distinct from the Ensemble Kalman filter variants in Section 5, which compute a one-step analysis update.

9.3. Parallelisation of EnDA

The need to integrate an ensemble of model states leads to large computational costs, because instead of computing a single model integration as in normal modelling applications an ensemble of $\mathcal{O}$(10–100) members has to be propagated. To reduce the time to perform the costly computations one can apply parallelisation of the data assimilation program and then use high-performance computers with a large number of processors to perform the computations. The ensemble integrations as the most costly part of the computations can be easily parallelised. In fact, the integration of each ensemble state is independent from the other states. Thus, this step could be parallelised by simply starting the numerical model $N_e$ times. Each model state has to be initialised from a different restart file and one has to store the final state of each model integration to keep the information on the forecast ensemble. Subsequently to the ensemble forecasts, one starts the data assimilation program, which reads the ensemble information from the files, computes the analysis step, and writes a set of new restart files to prepare the next forecast phase. The computations of the analysis step can also be parallelised as is outlined below. This implementation scheme of data assimilation can be termed ’offline coupling’ (Nerger and Hiller, 2013). While being flexible, the frequent writing and reading of the large files holding the ensemble states can take a significant amount of time.

A more sophisticated parallelisation of the ensemble data assimilation problem with a high-dimensional ocean model was discussed by Keppenne and Rienecker (2002) and Keppenne and Rienecker (2003). This method applied a domain-decomposition to the model and then integrated several ensemble states concurrently. The forecast ensemble was then collected by the use of the parallelisation technique SHMEM, which was also used for exchanging data in between processors during the analysis step of the EnKF applied in this study. Keeping the analysis step and the ensemble forecasts within one program reduced the overall computing time because the writing and reading of model state files is reduced.

The analysis step of the ensemble filters can also be parallelised using parallelisation methods like the Message Passing Interface (MPI, Gropp et al., 1994). The parallelization differs depending on whether localisation is used and on which of the filters is used. For the filter methods that assimilate all observations at once (in contrast to the serial observation processing of the EAKF and EnSRF) using the domain-decomposition of a model was found to be more efficient than using ensembles which are distributed over several processors because the amount of data that has to be exchanged using MPI is smaller for domain-decomposition (Nerger et al., 2005a).

For the ensemble Kalman filters with domain localisation, the local analysis update is independent for each local domain. Thus, this part is naturally parallel and can be distributed with MPI, the shared-memory standard OpenMP, or a combination of both. However, because the observation domains have a larger spatial extent they can reach into the grid domain held by neighbouring processors. The local analysis step needs the difference (innovation) between the observation and the corresponding part of the observed state vector. These differences need to be first computed by the processor that holds the sub-domain and then exchanged in between the different processors computing the analysis step. This computation of the observation innovations and their exchange using MPI is only required once before the loop over all local analysis domains can be computed in parallel (Nerger and Hiller, 2013). The cost for these operations depends on the total number of observations and on their distribution over the model grid. For many observations this can limit the parallel speedup of the analysis update as was shown for the localised SEIK filter by Nerger and Hiller (2013).

The EAKF and EnSRF are typically applied with serial observation processing and covariance localisation. In this case, the parallelisation of the analysis step has to take into account that for each assimilated observation the full model state has to be updated. Hence, also the innovation differences between the not yet assimilated observations and the corresponding observed model state change after each update. Anderson (2007) proposed to let each processor separately update the innovations so that the required parallel communication is limited. This parallelisation does not take the localisation into account. Taking into account that the localisation results in a limited reach of the observation influence, Wang et al. (2013) proposed another parallelisation strategy.

The analysis step of the ensemble Kalman filters requires only the model states. This allows for a generic coupling between the model and the analysis step. In particular, one can implement filter algorithms such that they can be coupled in the same way with different models. This allows one to build generic frameworks for ensemble data assimilation (Nerger et al., 2005a; Nerger and Hiller, 2013; Browne and Wilson, 2015). In the generic form, the ensemble forecast can still be computed by concurrent parallel model forecasts. The transfer of the forecast state information can then be performed either directly in memory by subroutine calls (Nerger and Hiller, 2013) or by parallel communication using MPI (Nerger, 2004; Browne and Wilson, 2015). These strategies allow a tight ’online’ coupling of the model and the data assimilation code that computes the analysis updates. The coupling can be achieved with minimal changes in the model code.

For the implementation of the EWPF and IEWPF different parallelisation schemes are applicable for the computations at each nudging step in between observations (Equations (120(120) $$\begin{array}{cc}\hfill {\mathsf{x}}_j^{(k)}& ={\mathcal{M}}_k\left({\mathsf{x}}_j^{(k-1)}\right)+{\tilde{\beta }}_j^{(k)}+\hfill \\ \hfill & +\widetilde{\mathrm{{\rm Y}}}\left[{\mathsf{y}}^{(m)}-{\mathcal{H}}_k\left({\mathsf{x}}_j^{(k-1)}\right)\right],\hfill \end{array}$$(120) ) and (123(123) $$\begin{array}{cc}\hfill w_j^{(k)}& \propto \frac{p\left({\mathsf{x}}_j^{(k)}|{\mathsf{x}}_j^{(k-1)}\right)}{q\left({\mathsf{x}}_j^{(k)}|{\mathsf{x}}_j^{(k-1)},{\mathsf{y}}^{(m)}\right)}{w}_j^{(k-1)}\hfill \\ \hfill & \propto \prod _{t=1}^k\frac{p\left({\mathsf{x}}_j^{(t)}|{\mathsf{x}}_j^{(t-1)}\right)}{q\left({\mathsf{x}}_j^{(t)}|{\mathsf{x}}_j^{(t-1)},{\mathsf{y}}^{(m)}\right)}.\hfill \end{array}$$(123) )) and at observation time for the EWPF between Equations (106(106) $$\begin{array}{cc}\hfill {\varphi }_j& ={\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T{\left(\mathsf{H}Q{H}^T+\mathsf{R}\right)}^{-1}\hfill \\ \hfill & \times \left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right).\hfill \end{array}$$(106) ) to (111(111) $$\begin{array}{cc}\hfill b_j& =\frac{1}{2}{\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T\hfill \\ \hfill & \times {\mathsf{R}}^{-1}\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)+\hfill \\ \hfill & log{w}_{\text{target}}-log{w}_j^{(m-1)},\hfill \end{array}$$(111) ) and for the IEWPF between Equations (112(112) $$\begin{array}{c}\hfill q\left({\mathsf{x}}^{(m)}|{\mathsf{x}}_{1:N_e}^{(m-1)},{\mathsf{y}}^{(m)}\right)=\frac{q\left({\xi }_j\right)}{\left||{\frac{\mathrm{d}\mathsf{x}}{\mathrm{d}\xi }}_j\right||}\end{array}$$(112) ) to (114(114) $$\begin{array}{cc}\hfill w_j& =\frac{p\left({\mathsf{x}}_j^{(m)}|{\mathsf{x}}_j^{(m-1)},{\mathsf{y}}^{(m)}\right)p\left({\mathsf{y}}^{(m)}|{\mathsf{x}}_j^{(m-1)}\right)}{q\left({\xi }_j\right)}\left||\frac{\mathrm{d}\mathsf{x}}{\mathrm{d}{\xi }_j}\right||w_j^{(m-1)}\hfill \\ \hfill & =w_{\text{target}}\hfill \end{array}$$(114) ). Before the observation time, the computations for the random forcing ${\tilde{\beta }}_j^{(m)}$ in Equation (120(120) $$\begin{array}{cc}\hfill {\mathsf{x}}_j^{(k)}& ={\mathcal{M}}_k\left({\mathsf{x}}_j^{(k-1)}\right)+{\tilde{\beta }}_j^{(k)}+\hfill \\ \hfill & +\widetilde{\mathrm{{\rm Y}}}\left[{\mathsf{y}}^{(m)}-{\mathcal{H}}_k\left({\mathsf{x}}_j^{(k-1)}\right)\right],\hfill \end{array}$$(120) ) are independent for each particle since a different forcing is drawn from the covariance for each of them. Similar, the nudging term in Equation (120(120) $$\begin{array}{cc}\hfill {\mathsf{x}}_j^{(k)}& ={\mathcal{M}}_k\left({\mathsf{x}}_j^{(k-1)}\right)+{\tilde{\beta }}_j^{(k)}+\hfill \\ \hfill & +\widetilde{\mathrm{{\rm Y}}}\left[{\mathsf{y}}^{(m)}-{\mathcal{H}}_k\left({\mathsf{x}}_j^{(k-1)}\right)\right],\hfill \end{array}$$(120) ) and the update of the weights are independent for each particle. Thus, these operations can be performed in parallel and there is no need to gather all particles on a single process. The computation of the matrix $\mathrm{{\rm Y}}$ in Equations (108(108) $$\begin{array}{c}\hfill \mathrm{{\rm Y}}=\mathsf{Q}{\mathsf{H}}^T{\left(\mathsf{H}\mathsf{Q}{\mathsf{H}}^T+\mathsf{R}\right)}^{-1}.\end{array}$$(108) ) and (117(117) $$\begin{array}{cc}\hfill {\mathsf{x}}_j^a& =\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)+\mathrm{{\rm Y}}\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right),\hfill \end{array}$$(117) ) is computationally the most expensive part. When the observation operator does not change over time, this matrix can be precomputed before beginning the assimilation. The downside of this approach is that this matrix can be huge and requires a lot of memory if the state dimension and number of observations are large. Otherwise, since the same matrix is used by all particles, it is possible to distribute the computation to all processes allocated for the particles, e.g. using a parallel matrix solver. At observation time, most of the computations are again independent for all particles. Only the maximum weight obtained from Equations (106(106) $$\begin{array}{cc}\hfill {\varphi }_j& ={\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T{\left(\mathsf{H}Q{H}^T+\mathsf{R}\right)}^{-1}\hfill \\ \hfill & \times \left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right).\hfill \end{array}$$(106) ) and (118(118) $$\begin{array}{cc}\hfill {\varphi }_j& ={\left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right)}^T{\left(\mathsf{H}\mathsf{Q}{\mathsf{H}}^T+\mathsf{R}\right)}^{-1}\hfill \\ \hfill & \times \left({\mathsf{y}}^{(m)}-\mathsf{H}\mathcal{M}\left({\mathsf{x}}_j^{(m-1)}\right)\right).\hfill \end{array}$$(118) ) for EWPF and IEWPF, respectively, must be exchanged over all processes holding particles, so that the target weight $w_{\text{target}}$ can be computed. Further parallelisation, e.g. to use the domain decomposition of the model, might also be possible. However, the matrix $\mathsf{Q}$ is frequently implemented in form of operators. As the parallelisation is always dependent on the particular implementation of the matrix $\mathsf{Q}$ it cannot be generalised for all models.

9.4. Computational cost

In Section 5, we presented various ensemble-base Kalman filter methods in a clean mathematical way for ease of comparison and clarity. In Appendix 2, we give a practical and precise pseudo-algorithms on how to implement each method. Providing detailed operation counts for all the ensemble methods presented in this paper would be too lengthy but more importantly the actual operation count would depend on many details such as operators $\mathcal{H}$ and $\mathsf{R}$, which are case specific for the model and observations. The operation counts provided here have been obtained by counting them in the pseudo codes in Appendix 2.

Generally, the leading order of operation counts in the different filters are those that scale with third order in any of the dimensions $N_x$, $N_y$, and $N_e$. For the SEIK, ETKF, ESTKF, and the EAKF methods, the leading order of operation count is $O\left(N_y{N}_e^2+N_e^3+N_x{N}_e^2\right)$ if the observation error covariance matrix is diagonal. The main cost is the update of the ensemble by multiplying with the weight matrix in Equation (19(19) $$\begin{array}{c}\hfill {\mathsf{X}}^a={\overline{\mathsf{X}}}^f+{{\mathsf{X}}^{\prime }}^f\left(\overline{\mathsf{W}}+{\mathsf{W}}^{\prime }\right),\end{array}$$(19) ) which has a complexity of $O\left(N_x{N}_e^2\right)$. Computing the matrix $\mathbf{T}{\mathbf{T}}^T$, e.g. in Equation (25(25) $$\begin{array}{c}\hfill \mathsf{T}{\mathsf{T}}^{\mathsf{T}}=\mathsf{I}-{\mathsf{S}}^{\mathsf{T}}{\mathsf{F}}^{-1}\mathsf{S}={\left(\mathsf{I}+\frac{1}{N_e-1}{\mathsf{S}}^{\mathsf{T}}{\mathsf{R}}^{-1}\mathsf{S}\right)}^{-1}.\end{array}$$(25) ) involves multiplications of matrix $\mathbf{S}$ with ${\mathbf{R}}^{-1}$ which has a complexity of $O\left(N_y{N}_e^2\right)$. Finally computing the transform matrix $\mathbf{T}$ by a Cholesky decomposition or an EVD has a computational complexity of $O\left(N_y{N}_e^2\right)$. While the leading order of operation counts is identical for all four filters, the SEIK and ESTKF are in general computationally faster than the ETKF, or the bulk formulation of the EAKF, despite equal leading operation counts due to details in the algorithms. For the EAKF, computing the SVDs of the matrices ${\mathbf{X}}^{\prime f}$ and $\tilde{\mathbf{S}}$, whose costs scale with $O\left(N_x{N}_e^2\right)$ and $O\left(N_y{N}_e^2\right)$, respectively, increases the computing time without changing the leading order of operation counts. Thus, the leading order of operation count does not reflect the computing speed.

The serial observation handling that is usually applied in the EAKF and EnSRF leads to an operation count of $O\left(N_y{N}_x{N}_e\right)$ in the leading order. Because only the ensemble updates are of third order complexity in the serial update, it can be faster than the bulk updates that assimilate all observations at once. This is even the case when localisation is used. However, in combination with localisation, the stability of the serial formulations can be deteriorated (Nerger, 2015).

The leading order operation count of the stochastic EnKF with perturbed forecasted observations is $O(N_y{N}_e^2+N_y^2{N}_e+N_y^3+N_x{N}_e^2)$. Here again the ensemble update, which scales as $O\left(N_x{N}_e^2\right)$, is usually the most costly operation. However, the EnKF is usually more costly than the filters mentioned before because of the inversion of the $N_y\times N_y$ matrix ${\mathsf{F}}_F$ (Equation (69(69) $$\begin{array}{c}\hfill {\mathsf{F}}_F=\mathsf{S}{\mathsf{S}}^{\mathsf{T}}+{\mathsf{Y}}^{\prime }{{\mathsf{Y}}^{\prime }}^{\mathsf{T}}\end{array}$$(69) )), which has a complexity of $O\left(N_y^3\right)$. Parallelising this inversion can help to reduce the computing time. A computing cost $O\left(N_y^3\right)$ also occurs for the bulk formulation of the EnSRF due to the EVD computed in Equation (58(58) $$\begin{array}{c}\hfill {\mathsf{F}}^{-1}={\mathrm{\Gamma }}_R{\mathrm{\Lambda }}_R^{-1}{\mathrm{\Gamma }}_R^{\mathsf{T}}.\end{array}$$(58) ).

When localisation is used, the change in the cost compared to the global formulation depends on the localisation method used. For covariance localisation (Section 9.1.2), the cost for computing the weight in matrix $\mathsf{C}$ and to apply it to ${\mathsf{P}}^f$ or the matrices ${\mathsf{P}}^f{\mathsf{H}}^T$ and $\mathsf{H}{\mathsf{P}}^f{\mathsf{H}}^T$ is added. For observation localisation (Section 9.1.3), the cost to compute the analysis for a local analysis domain with the bulk update methods is $O(N_{y,\gamma }{N}_e^2+N_e^3+N_{x,\gamma }{N}_e^2)$, where $N_{y,\gamma }$ is the number of local observations and $N_{x,\gamma }$ is the size of a local state vector that is corrected. Because both $N_{y,\gamma }$ and $N_{x,\gamma }$ are usually much smaller than the global dimensions $N_y$ and $N_x$, a single local analysis update is cheaper than the global update. However, the local analysis update has to be computed for each local analysis domain. Thus, the cost for the analysis with observation localisation is usually significantly higher than the global analysis. However, the local analysis can be easily parallelised to reduce the computing time as was described in Section 9.3.

The computing cost in the ETKF can be reduced using a projection matrix $\mathsf{A}$ analogous to the SEIK and ESKTF methods. For the ETKF, this projection is square-matrix with diagonal entries of $1-1/N_e$ and off-diagonal entries of $-1/N_e$. The advantage of using $\mathsf{A}$ is that one can avoid the explicit computation of ${\mathsf{X}}^{\prime }$ in favour of applying $\mathsf{A}$ to smaller matrices when evaluating the analysis equations (Nerger et al., 2012a).

Table 2. Overview of filter methods available from Sangoma project website.

Filter	Section	Comment
EnKF	5.1	includes covariance localisation
EnsRF	5.6	includes covariance localisation
ETKF/LETKF	5.4	ETKF is without localisation; LETKF includes domain/observation localisation
ESTKF/LESTKF	5.9	ESTKF is without localisation; LESTKF includes domain/observation localisation
NETF	7.1	without localisation

The computational cost for the particle-based non-linear filter NETF (Section 7.1) is similar to that of the ETKF since the analysis is performed in the $N_e$-dimensional subspace spanned by the ensemble members. In addition, the NETF does not compute an inverse matrix thus avoiding computational instabilities caused by small singular values, which are sometimes neglected in ETKF implementations for that reason (Sakov et al., 2012). If localisation is applied to the NETF, the local analysis computations are independent and can be evaluated in parallel as for the ETKF. The generation of random rotation matrices consumes additional resources; however, it is possible to resort to a collection of pre-calculated random matrices since they only depend on ensemble size $N_e$.

9.5. Ensemble data assimilation and non-linearity

The original EnKF was developed to overcome stability problems of the extended Kalman filter (see Jazwinski, 1970) that were discovered with ocean data assimilation applications (Evensen, 1993). Due to its use of an ensemble to propagate the state error covariance matrix, the EnKF is suited for non-linear models in this phase. However, the analysis step is based on the Kalman filter and is only optimal for Gaussian distributions. Obviously, a non-linear model forecast will transform a Gaussian distribution into a non-Gaussian distribution. Hence, the optimality of the Kalman filters is no longer preserved and the estimated analysis state and the error estimates will be suboptimal. This is a common issue for all ensemble filters whose analysis step is based on the equations of the Kalman filter. Nonetheless, the many existing data assimilation studies with non-linear models, e.g. of the ocean or atmosphere, with different formulations of the ensemble Kalman Filters show that these filters are rather stable with regard to non-linearity.

Second-order accurate ensemble filters, like the NETF, MMEF and MPF in Section 7 as well as the adaptive Gaussian mixture filter described in Section 8 avoid the assumption that the forecast ensemble has a Gaussian distribution. Thus, they should be better suited for non-linear systems. When the methods are applied with localisation, they can also be applied with large systems (e.g. the NETF in Tödter et al., 2016). However, filters like the NETF are still approximations to the full non-linear analysis that is performed by particle filters.

Particle filters do not rely on any assumption on the error distribution of the state estimate. However, the observation errors are frequently assumed to be Gaussian as is for the particle filters presented in Section 6. Additionally, while the EWPF and IEWPF do not require knowledge of forecast errors, they both require good knowledge of model errors, i.e. $\mathsf{Q}$. Of course, good knowledge of model errors is always beneficial to forecasting irrespective of the data assimilation method used, but for the application EWPF and IEWPF model errors are essential. While standard particle filters suffer the curse of dimensionality when applied to large systems, the EWPF and IEWPF by construction are designed to work for high-dimensional models, including those which are highly non-linear, with a small number of particles, e.g. Ades and van Leeuwen (2013) and Zhu et al. (2016). However, when applying the relaxation scheme between observations in an EWPF it is important to keep in mind that one has to choose the relaxation term $\widetilde{\mathrm{{\rm Y}}}$ in Equation (120(120) $$\begin{array}{cc}\hfill {\mathsf{x}}_j^{(k)}& ={\mathcal{M}}_k\left({\mathsf{x}}_j^{(k-1)}\right)+{\tilde{\beta }}_j^{(k)}+\hfill \\ \hfill & +\widetilde{\mathrm{{\rm Y}}}\left[{\mathsf{y}}^{(m)}-{\mathcal{H}}_k\left({\mathsf{x}}_j^{(k-1)}\right)\right],\hfill \end{array}$$(120) ) very carefully. We can choose this term to suit the needs of our model and indeed we need to do so carefully by selecting an appropriate relaxation strength function $\rho$ and covariance matrix. The relaxation term can be chosen to be constant between observation times, but that would not be a good idea if the system experiences oscillations between observations. In that case, the strength term can be chosen to be linearly increasing with the time lag to the next observation we are nudging particles to, or non-linear with maximum strength close to the observations.

In all local particle filters that we discussed the posterior particles are linear combinations of the prior particles. This has the potential to break non-linear balances between variables in the model. However, the linear combinations are typically formed such that only prior particles are added that are close to each other in state space, and hence quite similar. So this is not necessarily a disadvantage.

10. Summary and conclusion

This overview paper provides a coherent algorithmic summary and highlights differences between many currently used ensemble data assimilation methods that can be applied to high dimensional and non-linear problems such as ocean or weather prediction including well-known ensemble-based Kalman filters as well as recently developed particle filter methods and the Gaussian mixture filter.

We have presented these methods in a mathematically coherent way allowing the reader to compare many methods easily. In particular, we have presented all ensemble-based Kalman filter methods in form of a square root filter. In addition, we have included practical pseudo-algorithms for all methods since for computational reasons many of them would not be implemented in the form they are mathematically described. For some of the particle filters and for the Gaussian mixture filter we have presented the theory along with the step by step algorithm.

Finally, we have discussed important issues for practical implementation of the ensemble methods including various methods of localisation, inflation, parallelisation, computational cost and ensemble applicability to non-linear problems.

Concluding, a wealth of ensemble-based data-assimilation methods have been developed, and although they seem quite different in theory, the numerical implementations are quite similar. The implementations turn out to be quite similar to those for the particle filters, even those that explore a proposal density, where the state covariances that play an essential role in the Kalman filters are replaced by the covariance of the model errors. The main difference is that the state covariances are evolving over time and are always of low rank, while the model error covariance is given and of full rank but sparse. This means that different numerical algorithms need to be used to solve the equations when the system of interest has a high dimension.

11. Code availability

We note that many of the algorithms here have been efficiently implemented as a part of the Sangoma project and are freely available to everyone on the project website http://sourceforge.net/p/sangoma/ along with many other tools useful to data assimilation. Table 2 provides a list of available filters. Please note that the filter implementation was done independently from this paper so that not all filters described here are available. For simplicity these filters have been implemented without parallelisation and are hence only usable for moderately large problems with a state dimension of O(${10}^5$).

Further, all of these analysis methods have been implemented in at least one of the toolboxes connected to the Sangoma project, these are: EMPIRE, OAK, SESAM, OPENDA, BELUGA/SEQUOIA, NERSC and PDAF. For example, the set of filters listed in Table 2, plus the SEIK filter (5.2) with localization, are available in a parallelised implementation for high-dimensional problems in the freely available data assimilation framework PDAF (Nerger et al., 2005a; Nerger and Hiller, 2013). Further, the EWPF (Section 6.2.1) and IEWPF (Section 6.2.2) are available in EMPIRE (Browne and Wilson, 2015).

Acknowledgements

PJvL thanks the European Research Council (ERC) for funding of the CUNDA project under the European Unions Horizon 2020 research and innovation programme.

Additional information

Funding

This work was supported by the SANGOMA EU Project [grant number FP7-SPACE-2011-1-CT-283580-621 SANGOMA].

Notes

No potential conflict of interest was reported by the authors.

1 www.data-assimilation.net/.

2 The discussion of the increasingly growing developments in hybrid data assimilation methods is beyond the scope of this paper, instead we refer the reader to a very recent review article by Bannister (2017), and papers by Frei and Künsch (2013) and Chustagulprom et al. (2016) aiming to bridge particle and ensemble Kalman filter methods.

3 Note that $\rho =1$ if $\widehat{p}={\tilde{N}}_e$.

4 Many of the analysis methods discussed in this paper including MRHF have been implemented in Sangoma and are available for free to download from www.data-assimilation.net, as well as many other data assimilation tools for diagnostics, utilities etc..

5 Interestingly, the ECMWF is using an ensemble of 4DVars for their weather forecasting scheme, and it is relatively easy to turn this into a set of particles using 4DVar as proposal (see e.g. van Leeuwen et al., 2015).

6 The model error covariance matrices are usually assumed to be equal, i.e. $\tilde{\mathsf{Q}}=\mathsf{Q}$. .

Appendix 1

Resampling methods

In this section, we give descriptions of a number of resampling techniques that can be applied to the particle filter and Gaussian mixture filter methods to turn weighted particles into equal-weight particles. The resampling techniques included here are probabilistic resampling, stochastic universal resampling and residual resampling. However, they are by no means exclusive and other techniques could be used.

Probabilistic resampling (PR)

The probabilistic resampling or the basic random resampling is the most straightforward to implement as we sample directly from the density given by the weights.

Given the weights ${\{w_j\}}_{j=1}^{N_e}$ associated with the ensemble of particles, where the sum of weights is equal to one, the total number of particles $N_e$ and the number of particles to be generated ${\tilde{N}}_e$, we generate an index of the sampled particles using the Algorithm .

The required input for the PR is: $w\in {\mathcal{R}}^{N_e}$ a vector of particle weights, $N_e$ the total number of particles in the filter, and ${\tilde{N}}_e$ the number of particles to be sampled and the method returns an index $I\in {\mathcal{R}}^{{\tilde{N}}_e}$, which can then be used to select the sampled particles ${\mathsf{x}}_j^{\ast }={\mathsf{x}}_{I(j)}$ for $j=1:{\tilde{N}}_e$.

Note that this scheme introduces sampling noise by drawing ${\tilde{N}}_e$ times from a uniform distribution.

Stochastic universal resampling (SUR)

Stochastic universal resampling is also known as systematic resampling. It performs resampling in the same way as the basic random resampling algorithm except instead of drawing each $u_j$ independently from $\mathcal{U}(0,1)$ for $j=1,\dots ,N_e$, it uses a uniform random number u according to $u\sim \mathcal{U}[0,1/N_e]$ and $u_j=u+(j-1)/N_e$ (Bolić et al., 2003).

Given the weights ${\{w_j\}}_{j=1}^{N_e}$ associated with the ensemble of particles, where the sum of weights is equal to one, the total number of particles $N_e$ and the number of particles to be generated ${\tilde{N}}_e$, we generate an index of the sampled particles using the Algorithm .

The required input for the SUR is: $w\in {\mathcal{R}}^{N_e}$ a vector of particle weights, $N_e$ the total number of particles in the filter, and ${\tilde{N}}_e$ the number of particles to be sampled and the method returns an index $I\in {\mathcal{R}}^{{\tilde{N}}_e}$ which can then be used to select the sampled particles ${\mathsf{x}}_j^{\ast }={\mathsf{x}}_{I(j)}$ for $j=1:{\tilde{N}}_e$.

Note, that this method has a lower sampling noise than probabilistic resampling since only one random variable is drawn.

Residual resampling (RR)

The RR algorithm samples the particles in two parts. In the first part the number of replications of particles is calculated, but since the method does not guarantee that the number of resampled particles is $N_e$, the residual $N_r$ is computed. The second step requires resampling, which produces $N_r$ of the final ${\tilde{N}}_e$ particles. In Algorithm  this is done by PR, but other resampling technique can be used.

The required input for the RR is: $w\in {\mathcal{R}}^{N_e}$ a vector of particle weights, $N_e$ the total number of particles in the filter, and ${\tilde{N}}_e$ the number of particles to be sampled and the method returns an index $I\in {\mathcal{R}}^{{\tilde{N}}_e}$ which can then be used to select the sampled particles ${\mathsf{x}}_j^{\ast }={\mathsf{x}}_{I(j)}$ for $j=1:{\tilde{N}}_e$. Note, that we used the PR method to obtain an array $IR\in {\mathcal{R}}^{N_r}$ with the indices of the additional sampled particles, which we then stored in the remaining empty cells of the index array $I\in {\mathcal{R}}^{N_e}$.

Note, that this method reduces the sampling noise, but not as much as the SUR method.

Appendix 2

Filter algorithms for practical implementation

This Appendix contains practical pseudo-algorithms of all the ensemble filter methods presented in Sections 5–7. To discuss the computational cost of each method in Section 9.4 we used the algorithms presented in this appendix because for some filter methods they are more computationally efficient or numerically stable than mathematically elegant versions given in Section 5. The algorithms are written in the way one would implement them efficiently in Fortran. For compactness, the algorithms don’t show that the final step of the ensemble filters can usually be written in a blocked form, so that only the allocation one large ensemble array $\mathsf{X}$ is required. This is different for the MLEF, where two arrays of size $N_x\times N_e$ are required. If indices are given for matrices, the notation follows Fortran in that the first index defines the row, while the second index specifies the column.

Note, that in all the algorithms that follow in this appendix any variable with a number subscript is a temporary variable used to reduce the computational time and storage space needed for the algorithm. Further, for ease of reading these algorithms, we use abbreviations: SVD for singular value decomposition and EVD for eigenvalue decomposition. The values in the right column of each algorithm give the dimension of the resulting array, which helps to determine the computational cost of the operations.

Below, we use the notation $\mathcal{H}\left({\mathsf{X}}^f\right)$ as a shorthand notation for applying the possibly non-linear observation operator $\mathcal{H}$ individually to each ensemble state in ${\mathsf{X}}^f$.