Data Assimilation in Spatio-Temporal Models with Non-Gaussian Initial States—The Selection Ensemble Kalman Model

Data assimilation of spatio-temporal models is a challenge in many fields of study, including, but not limited to, air pollution mapping, weather forecast, petroleum engineering, and ground water flow assessment. Over the years, methods have been developed to handle increasingly complex problems. It started with the Kalman filter as presented in the seminal publication [1]. The Kalman filter is based on a Gaussian initial model and Gauss-linear forward and observation models. It defined the foundation for data assimilation and is still used in many assimilation studies. The extended Kalman filter (EKF) [2] appeared as a natural methodological extension that allowed for nonlinearity in the Kalman filter framework by linearization. The ensemble Kalman filter (EnKF) [3,4] defined a Monte Carlo approach to the filter and it became popular as it allowed for nonlinearity in the forward and observation models without having to evaluate analytical gradients. The EnKF and its variants have proven to be efficient in solving high-dimensional and nonlinear problems, see [5,6]. In the EnKF, the initial ensemble members represent the initial state which may not have an analytical expression. The forward model then propagates the ensemble members forward in time. Pseudo observations are generated using the observation model. The conditioning of each ensemble member is made with the Kalman weights estimated from the ensemble to give the best linear update. In cases where the initial model is non-Gaussian, the distribution of the variable of interest conditioned on the data will tend toward Gaussianity as observations are assimilated due to the linear assimilation rule.

Non-Gaussian initial distributions may be conserved by using a univariate transform into Gaussian marginals while assuming multi-Gaussianity in the transformed space. A univariate back transform is then used to return to the original space. This approach has a long history in traditional statistics, geostatistics, and more recently in ensemble methods for data assimilation, which is referred to as copulas [7], normal score transform [8], and Gaussian anamorphosis [9], respectively. The latter has shown to improve the performance of the EnKF in many applications [10,11]. There are however some unresolved issues since Gaussian anamorphosis transforms the marginal distributions rather than the full distribution, and the effect on the resulting variables interdependence is uncertain.

The Ensemble Randomized Maximum Likelihood Filter (EnRML) [12] and its close relative the Iterative EnKF (IEnKF) [6] are primarily used to handle nonlinearities in the forward and observation models, but they will also retain certain non-Gaussian features in the filtering distribution. These filters require gradient evaluations to execute the update which can be complicated even if the adjoint state method is used. One alternative is to evaluate the gradient using the ensemble itself [13], but this approach introduces an approximation with unclear consequences, particularly in models with multimodal marginals.

Multimodality in the prior model can be represented using categorical auxiliary variables to construct Gaussian mixture prior models [14,15,16]. In a spatial setting, these models appear as a combination of Gaussian random fields whose parameters depend on the value taken by the categorical variable, but in order to retain spatial dependence, the categorical variable must also have a spatial dependence. This indicator spatial variable can be modeled as a Markov [17] or truncated pluri-Gaussian [18] random field. For both of these models, there are challenges related to temporal data assimilation, although some encouraging examples have been developed [19].

We define and study an alternative prior model, the selection-Gaussian random field [20,21], which may represent multimodality, skewness, and peakedness. This random field model is conjugate with respect to Gauss-linear forward and observation models, similarly to the Gaussian random field model. The posterior distribution is therefore analytically tractable under these assumptions [22]. For general forward and observation models, ensemble based algorithms along the lines of the EnKF can be designed. Such selection ensemble Kalman algorithms are the focus of this study, and they are evaluated on a couple of examples.

In Section 2, we introduce the selection ensemble Kalman model. It provides a framework for the use of the selection-Gaussian distribution as a prior in data assimilation. This framework is then used for ensemble filtering and smoothing through the selection EnKF (SEnKF) and the selection EnKS (SEnKS) algorithms. In Section 3, a synthetic case study of the diffusion equation, with two distinct test cases, showcases the ability of the proposed approaches to assess a parameter field and the initial state of a dynamic field. Results from the SEnKF and the SEnKS are compared to that of the traditional EnKF and the EnKS, respectively. In Section 4, potential shortcomings are discussed and the results are put into perspective with respect to applicability in more realistic applications. In Section 5, conclusions are presented.

In this paper, $f(\mathit{y})$ denotes the probability density function (pdf) of a random variable $\mathit{y}$, ${\phi }_n(\mathit{y};\mathit{\mu },\Sigma )$ denotes the pdf of the Gaussian n-vector $\mathit{y}$ with expectation n-vector $\mathit{\mu }$ and covariance $(n\times n)$-matrix $\Sigma$. Furthermore, ${\Phi }_n(A;\mathit{\mu },\Sigma )$ denotes the probability of the aforementioned Gaussian n-vector $\mathit{y}$ to be in $A\subset {\mathbb{R}}^n$. We also use ${\mathit{i}}_n$ to denote the all-ones n-vector, ${\mathit{I}}_n$ to denote the identity $(n\times n)$-matrix and $\mathbf{1}(S)$ to denote the indicator function that equals 1 when S is true and 0 otherwise. We consider log-diffusivity to be an adimensional quantity and it will therefore not be given a unit.

Consider the unknown temporal n-vector ${\mathit{r}}_t$ for $t\in {\mathcal{T}}_r:\{0,1,\dots ,T,T+1\}$. Let $\mathit{r}=\{{\mathit{r}}_0,{\mathit{r}}_1,\dots ,{\mathit{r}}_T,{\mathit{r}}_{T+1}\}$ denote the variable of interest and let ${\mathit{r}}_{i:j}$ denote $\{{\mathit{r}}_i,{\mathit{r}}_{i+1},\dots ,{\mathit{r}}_j\},\forall (i,j)\in {\mathcal{T}}_r^2,i\le j$. Assume that the temporal m-vectors of observations ${\mathit{d}}_t$ for $t\in {\mathcal{T}}_d:\{0,1,\dots ,T\}$ are available, and define $\mathit{d}=\{{\mathit{d}}_0,{\mathit{d}}_1,\dots ,{\mathit{d}}_T\}$ and ${\mathit{d}}_{i:j}=({\mathit{d}}_i,\dots ,{\mathit{d}}_j\}$ accordingly. The model specified hereafter defines a hidden Markov (HM) model [23] as displayed in Figure 1.

Prior model: The prior model on $\mathit{r}$ consists of an initial and a forward model, where $f({\mathit{r}}_0)$ is the pdf of the initial state and $f({\mathit{r}}_{1:T+1}|{\mathit{r}}_0)$ defines the forward model.

(a) Initial distribution: The distribution for the initial state $f({\mathit{r}}_0)$ is assumed to be in the class of selection-Gaussian distributions [20,21]. Consider a Gaussian $(n+n)$-vector $[\tilde{\mathit{r}},\mathit{\nu }]$, with n-vectors ${\mathit{\mu }}_{\tilde{r}}$ and ${\mathit{\mu }}_{\nu }$, $(n\times n)$-matrix ${\Gamma }_{\nu |\tilde{r}}$, and where ${\Sigma }_{\tilde{r}}$, ${\Sigma }_{\nu }$, and ${\Sigma }_{\nu |\tilde{r}}$ are all three covariance $(n\times n)$-matrices with ${\Sigma }_{\nu }={\Gamma }_{\nu |\tilde{r}}{\Sigma }_{\tilde{r}}{\Gamma }_{\nu |\tilde{r}}^T+{\Sigma }_{\nu |\tilde{r}}$. Define a selection set $A\subset {\mathbb{R}}^n$ of dimension n and let ${\mathit{r}}_0=[\tilde{\mathit{r}}|\mathit{\nu }\in A]$; then, ${\mathit{r}}_0$ is in the class of selection-Gaussian distribution and its pdf is,

Note that the class of Gaussian distributions constitutes a subset of the class of selection-Gaussian distributions with ${\Gamma }_{\nu |\tilde{r}}=0\times {\mathit{I}}_n$. The dependence in $[\tilde{\mathit{r}},\mathit{\nu }]$ represented by ${\Gamma }_{\nu |\tilde{r}}$ and the selection subset A are crucial user-defined parameters with the latter being temporally constant. The selection-Gaussian model may represent multimodal, skewed, and/or peaked marginal distributions, see [21]. In this study, the initial distribution is defined to be a discretized stationary selection-Gaussian random field with parametrization,

For a given spatial correlation $(n\times n)$-matrix ${\Sigma }_{\tilde{r}}^{\rho }$, a stationary selection-Gaussian random field is fully parametrized by ${\mathbf{\Theta }}^{SG}=({\mu }_{\tilde{r}},{\mu }_{\nu },{\sigma }_{\tilde{r}},{\Sigma }_{\tilde{r}}^{\rho },\gamma ,A)$. Similarly, a stationary Gaussian random field is parametrized by ${\mathbf{\Theta }}^G=({\mu }_r,{\sigma }_r,{\Sigma }_r^{\rho })$.

(b) Forward model: The forward model given the initial state $[{\mathit{r}}_{1:T+1}|{\mathit{r}}_0]$ is defined as with where ${\omega }_t(·,·)\in {\mathbb{R}}^n$ is the forward model with random n-vector ${\mathit{ϵ}}_t^r$, independent and identically distributed (iid) for each t. This forward model may be nonlinear, but, since it only involves the variable at the previous time step ${\mathit{r}}_t$, it defines a first-order Markov chain. Note that $f({\mathit{r}}_{t+1}|{\mathit{r}}_t)$ cannot generally be written in closed form.

Likelihood model: The likelihood model for $[\mathit{d}|\mathit{r}]$ is defined as conditional independent with single-site response, with where ${\psi }_t(·,·)\in {\mathbb{R}}^m$ is the likelihood function with random m-vector ${\mathit{ϵ}}_t^d$, iid for each t. Note that $f({\mathit{d}}_t|{\mathit{r}}_t)$ cannot generally be written in closed form.

Posterior model: The posterior model for the HM model in Figure 1 is given by and is also a Markov chain, see [23,24]. This model is denoted the selection ensemble Kalman model. If the forward and likelihood models are Gauss-linear, the posterior model is also selection-Gaussian and analytically tractable, see [22]. When the forward and/or likelihood models are nonlinear, however, approximate or sampling based assessment of the posterior model must be made. For this purpose, we introduce the selection ensemble Kalman filter (SEnKF) and smoother (SEnKS) in the spirit of the traditional ensemble Kalman model [3].

The traditional EnKF algorithm aims at assessing the forecast pdf $f({\mathit{r}}_{T+1}|{\mathit{d}}_{0:T})$, and it is justified by general HM model recursions, see [23]. The algorithm is initiated by and utilizes the recursion for $t=1,\dots ,T$,

The expressions are represented by an ensemble of realizations, which in each recursion is conditioned using a linearized approximation with Kalman weights estimated from the ensemble. Thereafter, the ensemble is forwarded to the next time step. The SEnKF introduced in this study relies on the same relations as above, but it operates on the augmented $(n+n)$-vector $[{\tilde{\mathit{r}}}_{·},\mathit{\nu }]$, see Equation (2). Hence, the forward model is defined as where the auxiliary n-vector ${\nu }_t$ is temporally constant.

The likelihood model is defined as

The SEnKF algorithm provides an ensemble representation of and, based on this ensemble, empirical sampling based inference, see [21], is used to obtain the forecast of interest:

The SEnKF algorithm is specified in Algorithm A1 in Appendix A.

The traditional EnKS algorithm aims at evaluating the interpolation pdf $f({\mathit{r}}_{0:T}|{\mathit{d}}_{0:T})$ with corresponding HM model recursions, see [23]. The algorithm is initiated by and the recursions for $t=1,\dots ,T$,

The expressions are represented by an ensemble of realizations. Forwarding is made on the ensemble and the conditioning is empirically linearized. Note that the dimension of the model increases very fast, one may therefore only store the interpolation pdf $f({\mathit{r}}_s|{\mathit{d}}_{0:T})$ at the time point s of interest. The SEnKS introduced in this study relies on the relations defined above and uses an extended $(n+n)$-vector $[\tilde{\mathit{r}},\mathit{\nu }]$ as defined in Equation (2). The forward and likelihood models are identical to those defined for the filter. The SEnKS algorithm provides an ensemble representation of and by using empirical sampling based inference, see [21], the interpolation of interest is assessed,

The SEnKS algorithm is specified in Algorithm A2 in the Appendix A. Both algorithms, SEnKF and SEnKS, contain empirically linearized conditioning and asymptotic results, when the ensemble size goes to infinity, and are consistent only for Gauss-linear forward and likelihood models. Under these assumptions, the model is analytically tractable; however, see [22]. In spite of this lack of asymptotic consistency for general HM models, the ensemble Kalman scheme has proven surprisingly reliable for high-dimensional, weakly nonlinear models even with very modest ensemble sizes [25].

We consider two test cases to illustrate the relevance of the selection ensemble Kalman algorithms presented in Section 2. The model, common to both test cases, is based on the diffusion equation. The test cases are designed such that it will be opportune to consider bi-modal initial distributions. In the first test case, we compare the SEnKF to the traditional EnKF with a focus on predicting the diffusivity field that contains a high diffusivity channel. In the second test case, we compare the SEnKS to the traditional EnKS with a focus on evaluating the initial temperature field that is divided into two distinct areas where the initial temperature is substantially higher in one than in the other.

The traditional EnKF and EnKS algorithms provide an ensemble that directly represents the posterior models $f({\mathit{r}}_{T+1}|{\mathit{d}}_{0:T})$ and $f({\mathit{r}}_{0:T}|{\mathit{d}}_{0:T})$, respectively. Hence, the posterior models can be assessed by displaying statistics based on these ensembles. Reliable assessment of the posterior model in the two test cases can be obtained with approximately 1000 ensemble members. The SEnKF and SEnKS algorithms under study provide an ensemble of the augmented posterior models, $f({\tilde{\mathit{r}}}_{T+1},\mathit{\nu }|{\mathit{d}}_{0:T})$ and $f({\tilde{\mathit{r}}}_{0:T},{\mathit{\nu }}_{0:T}|{\mathit{d}}_{0:T})$, respectively. In order to obtain the posterior models of interest,$f({\mathit{r}}_{T+1}|{\mathit{d}}_{0:T})$ and $f({\mathit{r}}_{0:T}|{\mathit{d}}_{0:T})$, the conditioning on $\mathit{\nu }\in A$ must be made by empirical sampling based inference, see [21,22]. This inference requires the estimation of the expectation vector ${\mathit{\mu }}_{\tilde{r}\nu }$ and covariance matrix ${\Sigma }_{\tilde{r}\nu }$. The two test cases are defined on a $(21\times 21)$-grid for both ${\tilde{\mathit{r}}}_t$ and $\mathit{\nu }$—hence in dimension 882. The expectation and covariance will have 882 and 389,403 unique entries, respectively. Our experience from this study is that approximately 10,000 ensemble members are required to obtain reliable assessment of the posterior models of interest. To reduce the ensemble size, we have tested various localization approaches [29], without notable success, and leave the subject for further research.

Data assimilation of spatio-temporal variables with multimodal spatial histograms is challenging. Traditional ensemble Kalman algorithms enforce a regression towards the mean due to the linearized conditioning on observations, hence the multimodality is averaged out. We introduce the selection ensemble Kalman algorithms, termed SEnKF and SEnKS. These algorithms are based on recursive expressions similar to the ones justifying the traditional ensemble Kalman algorithms, but they are defined in an augmented space including the selection variable. From the two case studies, we conclude that multimodality is much better represented by the selection ensemble Kalman algorithms than by the traditional ones. We obtain RMSE reductions in the range of 20 to 30%.

The traditional ensemble Kalman algorithms provide an ensemble representation of the posterior model of interest hence making assessment of the posterior pdf simple. The selection ensemble Kalman algorithms are defined in an augmented space and conditioning on the selection variable must be made a posteriori. For this conditioning to be reliable, the ensemble size needs to be much larger than for the traditional algorithms. Hence, there is a trade-off between improved reproduction of multimodal characteristics of the phenomenon under study and the computational demands. In our case study, the ensemble size needed to be increased by approximately a factor of ten.

We have not fully explored the possibilities of robust estimation of model parameters in the conditioning of the selection variable. This robustification may reduce the ensemble size requirements. Note that parallelization in forwarding of the ensemble is possible and it will reduce the computer demands.