2.1 The simple smoother method

In Dong et al. (2021), a simple smoother method (hereafter referred to as the DHM smoother; it is an acronym made up of the initials of the authors' names) was presented for application in operational ocean reanalysis products, where the original analysis had been performed with a purely time-sequential approach, such as in forecasting situations for which future data are never available. This simple approach was designed to use the archive of increments to create a post hoc smoothing of the original reanalysis. Dong et al. (2021) showed the positive impact of this smoothing on a full ocean reanalysis and also on the low-dimensional Lorenz (1963, hereafter L63) system. The algorithm is as follows: let A_t be the forward-sequential (filtering) analysis at time t, and let I_t be the analysis increment field used to produce A_t. The smoother solution at time t is denoted S_t. The smoother algorithm is then written as follows. First, for S₀,

where $\mathrm{0}<{\mathit{\gamma }}_{\mathrm{a}}<\mathrm{1}$ is the increment decay rate per analysis time window, so that analysis increments from future analysis times decay in their influence on S₀. Similarly, for S₁,

By rearrangement, we obtain

with

where ${\mathrm{SI}}_t=S_t-A_t$ defines the smoother increment. These recursive relationships allow the smoother to be applied as a post-processing algorithm run backwards in time, starting with the final sequentially analysed time window and using the stored archive of filter increments. Later it will be convenient to define the increment decay per model time step, which we will write as γ, where γ^N=γ_a, and N is the number of model time steps between filter analyses. Then we will assume that each analysis window consists of one time step, with N=1 and γ=γ_a. The decay timescale τ associated with the smoothing is given, in time steps of δt, by

which is effectively a measure of the smoother lag. This equation can also be rearranged to be $\mathit{\gamma }=e^{-\frac{\mathit{\delta }t}{\mathit{\tau }}}$, implying an exponential increase in the forecast error. These smoother Eqs. (1) and (2) do not have a fixed-lag cutoff, unlike the fixed-lag smoothers discussed below.

We now discuss how this simple smoother is related to the conventional KS approach (a more formal proof of equivalence is given in the Appendix).

2.2 Extended Kalman filter and extended Kalman smoother

We start from the classical extended Kalman filter (ExtKF) and fixed-lag extended Kalman smoother (ExtKS) formulations, in which a tangent linear model is used for error covariance propagation when the model is nonlinear. Superscripts f, a, and s describe filter forecasts, filter analyses, and smoother analyses, respectively. The analysis of the Kalman filter at time k is given by

where the subscript represents the time step, x∈ℝⁿ is the n-dimensional state vector, y∈ℝ^m is the observations, and ℋ_k is the observation operator. In the ExtKF, the observation operator and the model can both be nonlinear, where the state vector evolves with a model ${\mathit{x}}_k=\mathcal{M}\left({\mathit{x}}_{k-\mathrm{1}}\right)$.

${\mathbf{K}}_k^{\mathrm{a}}\in {\mathbb{R}}^{n\times m}$ is the Kalman gain for the analysis, which is given by

with ${\mathbf{P}}_k^{\mathrm{f}}\in {\mathbb{R}}^{n\times n}$ being the forecast error covariance matrix, ${\mathbf{R}}_k\in {\mathbb{R}}^{m\times m}$ being the observation error covariance, and ^T being the transpose operator, all at time step k. Here, instead of the nonlinear observation operator, the tangent linear approximation, ${\mathbf{H}}_k\in {\mathbb{R}}^{m\times n}$, is used in the gain below, and the tangent linear model, ${\mathbf{M}}_k\in {\mathbb{R}}^{n\times n}$, is used for the propagation of the forecast error covariance matrix, ${\mathbf{P}}_k^{\mathrm{f}}={\mathbf{M}}_k{\mathbf{P}}_{k-\mathrm{1}}^{\mathrm{f}}{\mathbf{M}}_k^T$. Finally, the analysis error covariance can be derived from the forecast error covariance as follows:

For the fixed-lag ExtKS, Todling and Cohn (1996, hereafter TC96) derive backward-looking equations for the smoother, which run interleaved with every filter time step; however, here we will present forward-looking equations aimed at expressing the fully smoothed state, including contributions from multiple future filter steps, as presented for the simple smoother in Eqs. 1–5. The full equivalence between the TC96 notation and our notation is demonstrated in the Appendix.

The contributions from observations at time step k+ℓ to the smoother solution at time step k can be written in the same Kalman gain notation as follows:

We note that if the filter states are only stored at some assimilation frequency, e.g. once per day, then the index ℓ can be defined at the same frequency as these filter increments, as smoother states are only required at the same frequency as the stored filter states.

The full smoother solution for time step k is then obtained by the summation of smoother increments from all future time steps (here assumed to be truncated to the maximum lag L) as follows (see the Appendix):

The cross-time smoother gain matrix is simply a modified version of the standard ExtKF gain and can be written as follows:

There is a subtlety here because in the Appendix we will see that the cross-time error covariance ${\mathbf{P}}_{k,k+\mathrm{\ell }}$ is not independent of ${\mathbf{P}}_{k,k+\mathrm{\ell }-\mathrm{1}}$. However, this will not be the case in the simple smoother approximation as applied below.

To introduce the key simple smoother approximation, we rewrite the cross-time error covariance as a decay rate and consequently also neglect any interdependence of the smoother contributions from different times.

which is equivalent to assuming

This Eq. (13), when substituted into Eq. (10), clearly expresses the approximation being made to recover the simple smoother solution from the ExtKS equations.

When using Eq. (6), this gives

where ${\mathit{I}}_{k+\mathrm{\ell }}={\mathit{x}}_{k+\mathrm{\ell }}^{\mathrm{a}}-{\mathit{x}}_{k+\mathrm{\ell }}^{\mathrm{f}}$, thereby reproducing the simple smoother under the additional assumption that L≫τ in Eq. (5). The smoother is now defined entirely in terms of the sequential analysis increments that will therefore allow post-processing from an archive of increments from the sequential filter run. Another way to interpret this approximation is to say that the spatial and temporal error covariances in the KS are assumed to be separable, with the spatial (and cross-variable) error covariances being determined by the KF equations, but the temporal covariances (from times k+ℓ to k) being approximated by a simple decay. We will return to this description later when we seek to extend the approximations to the EnKS case.

It is also possible to make the equivalent approximations to the smoothed uncertainties. For each smoother increment introduced in Eq. (9), there will be a corresponding reduction in the smoother error covariance given by

so that the fully smoothed error covariance can be written as follows (see the Appendix):

Then, by making the simple smoother approximation, Eqs. (12) and (13) give

Now returning to use Eq. (8), we finally obtain

where ${\mathbf{IP}}_{k+\mathrm{\ell }}={\mathbf{P}}_{k+\mathrm{\ell }}^{\mathrm{f}}-{\mathbf{P}}_{k+\mathrm{\ell }}^{\mathrm{a}}$ are the filter error covariance increments, mirroring Eq. (15), forming the simple smoother equations for the increments. The smoothing Eqs. (15) and (19) could clearly both be written in a recursive format like that of Eq. (3) for easier post-processing. In the following sections, we investigate how well these approximations work through comparisons in the L63 system.

3.1 Assimilation set-up

A twin experiment using the L63 model was carried out to evaluate the smoother. The L63 uses a classical set-up with model equations.

where the standard model parameters are chosen as $\mathit{\sigma }=\mathrm{10},\mathit{\rho }=\mathrm{28},$ and $\mathit{\beta }=\frac{\mathrm{8}}{\mathrm{3}}$. The model parameter set-up is consistent with that in the DHM smoother (Dong et al., 2021). All experiments are run for 20 time units, with each time unit consisting of 100 time steps of 0.01. We performed a “truth” run first with x, y, and z values of 5 as the initial condition. Observations are assigned for x and y with a frequency of 5 and 20 time steps, respectively, and Gaussian errors added with standard deviations of 2. No observations are taken for z.

Dong et al. (2021) used a three-dimensional variational (3DVAR) method for assimilation into L63, with a fixed background error covariance prescribed as the time mean error covariance from a single separate L63 run. Here we ran the extended Kalman filter (ExtKF) and extended Kalman smoother (ExtKS) with 100 different initial conditions in which the background error covariance is modelled. However, with this assimilation frequency, we found that in order to avoid filter divergence, we needed a hybrid forecast error covariance that retained a 5 % weighting of the L63 climatological covariances in the background error. The smoother uses fixed-lag L=40 time steps, as we found that this lag gives the smallest errors in our experiments.

The simple smoother (DHM; Eqs. 1–5) was executed with γ_a=0.9 (equivalent to γ=0.9 in our experimental set-up), which is when the smoothing results have the smallest error compared to other values (experiments not shown). We also ran a modified Kalman smoother (hereafter MKS), using the approximated cross-time Kalman gain as in Eq. (13). This is implemented by directly substituting the approximation into the full KS equations described in the Appendix and is used to demonstrate DHM-equivalent results. The uncertainty estimates for the MKS smoother are also obtained in the same way.

3.2 ExtKS assimilation results

Across the 100 assimilation runs, we calculated the root mean square error (RMSE) time series against the truth for each smoothing method. Figure 1a, b, and c show a portion of the ($x,y,z$) RMSE time series, respectively, for the filter and the different smoothing methods in the solid lines. For most time steps, the KF errors are larger than the smoother errors. The full ExtKS has the smallest errors; however, the DHM and MKS are almost identical and lie in between those for the KF and KS. Also in Fig. 1, there are dashed lines representing the average of the smoothed standard deviation (SD) uncertainty estimates for the ExtKS runs (blue; Eq. 17) and MKS runs (green; Eq. 19), respectively. It should be emphasised that these uncertainty estimates are calculated entirely independently of the actual truth values themselves, which would not be known in a real assimilation problem. The level of agreement between these uncertainty estimates and the RMSE against truth is remarkable.

Figure 1The RMSE time series for the KF and KS, along with the MKS and the DHM smoother, in the L63 system for (a) x, (b) y, and (c) z. The RMSEs are averaged across 100 independent assimilation runs starting from different initial conditions and assimilating observations from the same 20 time unit truth run. Only two time units of the run are shown to allow the behaviour around the analysis times to be clearly seen. Dotted green lines and dotted blue lines (see the panel insets for details) are the posterior uncertainty SD (square root of error covariance P) estimates for the MKS and KS, respectively.

Download

The time mean RMSEs for $x,y,z$ are summarised in Table 1, along with the uncertainty SD, where calculated, over the entire 20 time units of the runs. Both DHM and MKS provide an improvement on the ExtKF results by 60 %–70 %, relative to the ExtKS improvements for x and y, although the RMSE for z is not reduced in DHM and MKS. This is perhaps because the instantaneous error covariances between z and the assimilated x,y variables, as used in DHM, are insufficient to improve z, whereas the full ExtKS allows some history of the x,y evolution to be used for deriving z smoothing increments. The RMSE numbers in parentheses are evaluated only at the filter update time steps, where observation data are assimilated. These errors are smaller than all-time-step RMSEs by ∼ 5 %, which is as a result of the data assimilation at these time steps. This is consistent with the RMSE time series in Fig. 1, where the Kalman filter (orange line) usually declines sharply when data are available. The smoother solutions, however, also yield much-improved analyses between observation time steps.

Figure 2As in Fig. 1 but the filter and smoother increments are shown for situations where the smoother increments, when applied, add to the filter increments.

Download

Figure 2a, b, and c show the actual $x,y,z$ increments, respectively, being introduced by the ExtKF and the different smoothers through the analysis times. In each case, the smoother increments are additional to the filter increments, which appear clearly as orange spikes for every five time steps when data are available. Smoother increments between analysis times can be seen, with the DHM and MKS increments being virtually identical and decaying backwards in time from each filter increment. The ExtKS increments are more complex, sometimes being similar to the DHM, but sometimes they can be considerably larger.

If we look at the mean ratio of the cross-time error covariances relative to the filter forecast error covariances in comparison to the simplified γ decay representation across time in Fig. 3, we can understand something about the performance of the smoothers. We do not expect these to be identical because the full cross-time smoother covariances for larger lags, ℓ, take account of the increments from intermediate times. For small lag values, the average cross-time error decay rates are fairly similar; however, for larger lag values, the model-derived cross-time error covariance on average takes the opposite sign. This happens on a similar timescale to the short oscillation period of x,y in L63 and is associated with the growing amplitude of these oscillations as they reach larger and smaller x,y values before the phase lobe transitions. This is a very model-specific behaviour, and the simple smoother's γ decay error covariances cannot represent this. This also explains why larger γ values make the simple smoother worse in L63 because the positive cross-time error covariances are used for larger lags when the negative cross-time error covariances are more appropriate.

The key point is that the simplified smoother DHM provides substantial improvement over the ExtKF, while incurring very little computational cost (no tangent linear model or TLM runs and no storage of cross-time error covariances) compared to the ExtKS. The DHM smoother can therefore be applied in its entirety through the post-processing of the filter results. While this was demonstrated in Dong et al. (2021), here we show more clearly how the equivalent MKS approximation is derived from the ExtKS equations, and we also show how the smoothed uncertainties can be cheaply post-processed and still give useful information.

In the next section, we extend the decay assumption for cross-time error covariances to apply to the ensemble Kalman filter and smoother equations, which are much more relevant to large nonlinear models for which the direct modelling of error covariances across time is infeasible.

Table 1Time mean RMSE against truth for the ExtKF and ExtKS, along with the modified MKS and the simplified DHM smoother for each variable in the L63 system (see also legend for Fig. 1). The numbers in parentheses are the mean RMSE values at observation time steps only (not independent data). The time mean of the standard deviations calculated for the uncertainties are also shown as SD. Time averaging is now over the entire 20 time units of each assimilation run.

Download Print Version | Download XLSX

Figure 3Cross-time error covariance decay rates for the ExtKS and the DHM smoother. For the ExtKS, the y axis is the smoother cross-time error covariance divided by the filter forecast error covariance and averaged for all times over the 20 time unit L63 run. This is compared to the decay of γ^ℓ used in the DHM smoother, as a test of Eq. (12).

Download

4.1 Approximating ensemble error covariances

In the ExtKF and ExtKS, a TLM propagates the flow-dependent error statistics, which are then used to calculate increments. However, the TLM reliability declines sharply with propagation time for a system as nonlinear as the L63 model. The ensemble Kalman filter (EnKF) can then give better results by estimating the error statistics with a finite ensemble of state realisations propagated by the full nonlinear model rather than by a TLM. This can then improve the quality of the forecast error covariance matrix. However, the update gains for the EnKF and ensemble Kalman smoother (EnKS) are defined identically to Eqs. (7) and (9), respectively, although the error covariances, ${\mathbf{P}}_k^{\mathrm{f}}$ (the error covariance at time step k) and ${\mathbf{P}}_{k,k+\mathrm{\ell }}$ (the error covariance between time steps k and k+ℓ) are calculated differently, since they are emulated from the limited ensemble of state vectors whose variability represents the uncertainty in the system. While ensemble filter methods have been adopted for larger environmental models, these have not generally added smoother steps because the cost to store, update, and apply posterior ensemble covariances still makes ensemble smoother methods generally infeasible.

However, these constraints can again be overcome by retaining the EnKF flow-dependent ensemble spreads to represent current errors, while making a simple decay approximation for the time shift error covariances, which is similar to our modified ExtKS in Sect. 2.2. For comparison purposes, we demonstrate this method in the general EnKS form in state space, keeping the same notation as we used previously, although the actual computation is performed in ensemble space, as in the ensemble transform Kalman filter (ETKF; e.g. Bishop et al., 2001; Zeng and Janjić, 2016).

where

is the normalised anomaly of the ith ensemble member from the mean in an N_e member ensemble of forecasts.

The first step in the ensemble filter is to update the ensemble mean,

and the second step is to update the uncertainty, using Eq. (8). Then, in order to regenerate the ensemble of analysis perturbations, the ETKF uses the following transformation:

where T_k is chosen to ensure that Eq. (8) is satisfied. Then, following the EnKF, the cross-time covariances in the EnKS can be expressed in state space as follows:

As explained in the Appendix, the full cross-time error covariances are calculated between the filter forecast (${\mathit{X}}_{k+\mathrm{\ell }}^{\mathrm{f}}$ in the EnKS) and previous partially smoothed states (X_k again in the EnKS), which require both past ensemble means and error covariances to be repeatedly smoothed. However, this is not necessary for the modified (simple) ensemble smoother, which we here call MEnKS. The ensemble mean smoothing using Eqs. 13 and 10 can be written as follows:

and then Eq. (18) can be used to obtain the smoothed uncertainties.

This is a great simplification because performing a full ensemble smoothing would require the whole past ensemble to be stored at all times and reprocessed. Equation (19) suggests that the error covariance increments must also be stored during the EnKF filter phase, which could still be a large storage requirement for a big model; however, in fact only the diagonal elements of P^s are likely to be of interest, i.e. the uncertainty variance of the state fields, or even just a subset of these, so only a smaller set of uncertainty increments may be of interest and worth storing for post-processing through Eq. (19). In the next subsection, we show the results of applying these approximations in L63.

4.2 EnKS assimilation results

Using the same L63 twin experiment as in Sect. 3, we solve a full smoother (EnKS) and use the modified (MEnKS) algorithm Eq. (28). To be consistent with the extended KS configuration, we use an ensemble size of 100 and a fixed lag of 40 time steps for smoothing. Figure 4a, b, and c show the $x,y,z$ RMSE and SD uncertainties, respectively, from these ensemble filter and smoother runs in the same format used in Fig. 1 for the extended Kalman filter.

These ensemble results are seen to produce lower RMSE values than the equivalent ExtKF or KS results (cf. Fig. 1), demonstrating the superiority of the ensemble assimilation method for dealing with the nonlinearity of the L63 model. Again, the EnKS substantially reduces the RMSE compared to the filter, and again, the approximated simple ensemble filter MEnKS gives intermediate RMSE results, with much less computational effort than the full EnKS. Although not as optimal as the EnKS, the simplified MEnKS shows a much smoother temporal evolution of the RMSE than the filter, which would be a significant improvement if, for example, applied to an ocean reanalysis. The post-processed uncertainty estimates also reproduce a reasonable estimate of the true RMSEs of the smoothed ensemble mean analyses.

Figure 5a, b, and c show the mean $x,y,z$ increment time series, respectively, for the ensemble filter and the two smoothers. The increments are smaller than those from the ExtKF/KS (Fig. 2), reflecting the improved assimilation approaches. Table 2 summarises the average RMSE and SD uncertainty results over the full 20 time units of the ensemble runs.

Figure 4As in Fig. 1 but for EnKF, EnKS, and MEnKS.

Download

Figure 5As in Fig. 2 but for EnKF, EnKS, and MEnKS.

Download

Table 2Time mean RMSE and SD uncertainties for each variable in EnKF, EnKS, and MEnKS in the L63 model, averaged over time units from 1–20. The numbers in parentheses are the RMSE values calculated at time steps with observations only (no independent data comparison). The lag is 40, and γ=0.9.

Download Print Version | Download XLSX