Using the (Iterative) Ensemble Kalman Smoother to Estimate the Time Correlation in Model Error

2.1 Problem formulation

Let $x\in \Re$ be the state variable of our system with initial conditions x₀ at t = 0. Consider τ independent model error jumps ϵ_t, for t = {1,…,τ} with zero mean. We denote the control variable as the column vector $z\in {\Re }^{\tau +1}$, with background distribution $N({\mu }_z^b,D_z)$:

where ${\mu }_z^b\in {\Re }^{\tau +1}$ and $D_z\in {\Re }^{(\tau +1)\times (\tau +1)}$ are the background mean vector and background error covariance matrix, respectively. The scalars b² and q² represent the background and model error variances. Note that z has a diagonal covariance matrix since we consider statistical independence amongst its elements. We use these elements to construct a simple linear system with time auto-correlated model errors in its evolution. The real linear evolution over one time step is prescribed by:

i.e. a simple auto-regressive component plus a model error realisation. This model error has distribution v(t) ~ N (0,q²), and the following structure in time:

Eq. (3) indicates that the model errors are auto-correlated in time, and this only depends on the lag |t–t|′ and a vector of N_θ parameters $\theta \in {\Re }^{N_{\theta }}$. For τ time steps, this yields an auto-correlation matrix $\Phi \in {\Re }^{\tau \times \tau }$. This symmetric Toeplitz has the following elements:

Being a positive-definite symmetric matrix, a Cholesky decomposition Φ is possible, i.e.:

where $L\in {\Re }^{\tau \times \tau }$ is lower triangular.

We now take Eq. (2) and write the time evolution of the system from t = 0 to t = τ in terms of the control vector z and the Cholesky factor L. Explicitly:

where $M_{0:\tau }\in {\Re }^{1\times (\tau +1)}$ is a row matrix with the model evolution from 0 to any time t. For this simple model, the elements are decreasing powers of α:

Let us define the composed evolution matrix $\widehat{M}\in {\Re }^{1\times (\tau +1)}$:

This includes the effect of the deterministic dynamics and the auto-correlation of the model error, and it will become useful in the next subsections. Separating the initial condition and the model errors, we write Eq. (6) as:

where $M_{1:\tau }\in {\Re }^{1\times \tau }$ and $ϵ\in {\Re }^{\tau }$ are obtained by removing the first element of M₀:_τ and z respectively. In fact, we can recover the value x_t at any time (not just the final time) using with the corresponding truncated elements in the second term.

We consider an observation of the truth at the end of the forecast window, i.e. at time t = τ. For simplicity, let the observation operator be the identity, so the observation equation is:

with the observation error: η ~ N(0,r²). Obtaining the analysis values for z was already discussed in Amezcua and van Leeuwen (2018), and Ren et al (2021). In this paper we discuss the solution of the joint state-variable estimation problem.

2.2 Bayesian solution for the joint state-parameter estimation

Considering both the control variable z and parameters θ to be random variables, the Bayesian solution of this problem is to obtain the posterior joint pdf of z and θ given the observation y. Namely,

The numerator is the joint pdf of z, θ and y. This is obtained as the product of the likelihood of y times the prior joint pdf of z and θ. If we consider these two to be statistically independent, then we have:

The marginal pdf of the observations is:

The prior distribution for the control variable and the likelihood are easy to characterise. Recall that we have:

From now on, we explicitly note that $\widehat{M}(\theta )$ depends on θ. Note that the joint estimation problem is complicated since we have the product $\widehat{M}(\theta )z$ in the likelihood, which limits the possibility of obtaining an analytical expression for the posterior p(x,θ|y). For this reason, we now discuss two solutions based in statistics of this pdf: the maximum-a-posteriori solution, and a popular approximation, the extended Kalman smoother solution.

2.3 The Maximum-a-posteriori solution

We can try maximising the joint pdf to obtain a maximum-a-posteriori (MAP) solution. This is equivalent of finding the minimum of the cost function:

Using the distributions in Eq. (14) and an arbitrary prior for θ, the minus logarithm of Eq. (12) is:

The minimisers {z*, θ*} of the cost-function can be obtained by taking the gradient of J(z, θ|y) with respect to both control variables and parameters and equating to zero:

with the gradients ${\nabla }_{z}J\in {\Re }^{\tau +1}$ and ${\nabla }_{\theta }J\in {\Re }^{N_{\theta }}$. If we also assume that the parameters follow a MVG –i.e. θ ~ (µ_θ, D_θ)– we are able to compute the gradients explicitly. This yields the following system of τ + 1 + N_θ equations:

with the Jacobian matrix $\frac{\partial \widehat{M}(\theta )}{\partial \theta }\in {\Re }^{(\tau +1)\times N_{\theta }}$ defined as:

One can solve z from Eq. (18a) to get:

with K and γ² defined as

and similarly from Eq. (18b):

Eq. (20) is the Kalman analysis equation for the posterior mean (Kalman, 1960; Kalman and Bucy, 1961). The complication, however, is in our case θ is an unknown. Hence, Eq. (20) needs to be solved in tandem with Eq. (22), which cannot be done analytically in general.

To actually calculate the derivative of the model with respect to the parameters, we note the following. For the j^th parameter θ_j, Eq. (19) can be readily computed using the Cholesky factor L defined in Eq. (5) in the following manner:

The derivative of the Cholesky matrix can be found using Theorem A.1 of Sarkka (2013):

where the matrix $T\in {\Re }^{\tau \times \tau }$ is defined as:

2.4 The extended Kalman Smoother solution

The extended Kalman Smoother solution can be derived directly from the MAP solution by centering the derivative of the model to the state and the parameters on the background values. Tracing back these derivatives we can rewrite Eq. (20) and Eq. (22) as:

Furthermore, in Eq. (27) we use a first-order Taylor expansion for $\widehat{M}(\theta )$ as:

We can solve then solve Eq. (26) and Eq. (27) by introducing a new augmented variable $\tilde{z}={(z^T,{\theta }^T)}^T$ and after some algebra we obtain:

where A_z is denoted as the posterior covariance in the extended variable space. In this case, the augmented gain $\tilde{K}\in {\Re }^{(\tau +1+N_{\theta })\times 1}$ and augmented residual variance ${\tilde{\gamma }}^2\in {\Re }_{\ge 0}$ are different from Eq. (21)

In the last step we have again assumed that state variables and parameters are statistically independent in the prior. Notice that even with a Gaussian prior p(θ) and a model that is linear in the state z the MAP and the extended KS solution can differ considerably. The MAP solution results in an implicit equation that needs to be solved iteratively. However, even if we manage to do this an uncertainty estimate is still lacking. It is possible, if the system is low-dimensional, to solve this problem with e.g. a particle filter. Unfortunately, for higher dimensions particle filters are degenerate, and a potential solution for that, localization, is difficult to use with global parameters. Instead, we resort to a popular approximate solution based on an iterative Ensemble Kalman smoother that can easily be scaled up to higher dimensions.

2.5 The (Iterative) Ensemble Kalman Smoother

The Iterative Ensemble Kalman Smoother (IEnKS) does solve the parameter estimation problem approximately, but it does introduce nonlinearity in the solution procedure. Hence it will provide a solution that is closer to the MAP solution than the Extended Kalman Smoother solution. In an IEnKS parameter estimation is typically performed via state augmentation, as in the extended Kalman Smoother. We define a variable $u\in {\Re }^{1+\tau +N_{\theta }}$ as:

where the augmented variable includes the state variable at all time steps and all the model error parameters, sometimes called the state formulation. Notice that this is slightly different than Eq. (1), where the augmented variable contains the initial conditions, the model error jumps, and the model error parameters, leading to the so-called forcing formulation. In a linear system, it is trivial to go from one formulation to the other, as explained in Tremolet (2006) and Amezcua and Van Leeuwen (2018). For ensemble implementations it is easier to work with state variables directly, so this is how we implement the experiments in the paper.

Ren, Amezcua, and Van Leeuwen (2021) studied the consequences of wrongly prescribing the model error decorrelation time scale, which was the only parameter. This is, they performed DA with a model error auto-correlation which was different from the one ‘guessed’ in the ensemble forecast. In the last part of the work, they used state augmentation in an EnKS –as formulated in Amezcua and Van Leeuwen (2018)– to estimate one time parameter in the model error. The authors did not, however, perform a systematic analysis with different types of model error structure. In this work we perform a deeper analysis and more extensive experiments.

In this work we perform state augmentation experiments using two DA methods, the stochastic EnKS and its iterative variant (IEnKS), as formulated in Evensen (2018). This formulation aims to find the minimum of a cost function using ensemble members, either in one step (EnKS) or multiple steps (IEnKS). In the IEnKS a cost function is defined for each ensemble member j by rewriting the cost function Eq. (16) in terms of the extended variable u as:

in which ${\mu }_u^j$ is the prior ensemble member j, and

is a perturbed observation with the perturbation drawn from the observation error pdf, so η^j ~ N(0,r²). In Eq. (32) it is clear that the product $\widehat{M}(u)u$ is responsible for the difficulty of this problem. Evensen (2018) proposed a simple Gauss-Newton iteration for each ensemble member as:

where i is the iteration index, $C\in {\Re }^{(1+\tau +N_{\theta })\times (1+\tau +N_{\theta })}$ is a symmetric approximation of the Hessian of the cost function in which the second derivative of $\widehat{M}(u)u$ is neglected, so

and $\nabla J^j(u)\in {\Re }^{1+\tau +N_{\theta }}$ is the gradient of the cost function for ensemble member j. The iteration step length δ is tunable. The iteration is started at $u_0^j={\mu }_u^j$. The simple descent algorithm can be used with the ensemble gradient and the ensemble approximated Hessian information, as shown in equation (36) of Evensen (2018). For efficiency, we use the particular implementation described in Evensen et al. (2019). Tuning factor δ is an art. In our implementation we reduce the step size the closer we get to the minimum via:

with 0<c<1. The stochastic EnKS can be obtained from this formulation by using only one iteration and setting δ = 1.

3.1 Model error formulations

The first type of model error was used both in Amezcua and Van Leeuwen (2018) and Ren, Amezcua, and Van Leeuwen (2021). The correlation of model error between two time steps is:

In this case, the auto-correlation decays exponentially, and it only has one parameter – the decay time scale ω (ω>0). When ω tends to 0 the model error becomes a white-noise process. When ω tends to infinity, the model error becomes fixed (a bias). Summarizing:

The second model error structure we explore contains both decaying and oscillatory elements:

The first term is a geometric memory term, with base –1 < φ < 1, which can be considered a decay parameter. The second term is oscillatory with frequency f. To avoid instability of the system, the decay parameter φ is bounded as 0 ≤ φ ≤ 1, and the frequency f is bounded by 0 < f < 0.5. When the frequency is 0.5, the covariance becomes purely decaying with time and it becomes purely oscillatory when the decaying parameter tends to 1.0. Summarising:

3.2 Evolution models

To illustrate the issues that arise when estimating model error parameters that are related to temporal correlations, we restrict our experimental set up to two zero dimensional models. The first is the simple linear scalar model used for the analysis in section 2. The true model is:

Note that this true model is a stochastic model and we assume v ∼ N(0,q²). The true model uses a real decay timescale ω_r, while the forecast model assumes a guessed time auto-correlation parameter ω_g. This value needs to be updated in the DA cycle.

The second set of experiments is much more ambitious since we use a non-linear map as the real system. In particular, we use the well-known logistic map (Robert, 1976):

where γ is the control coefficient that gives rise to stable solutions when 0 ≤ γ ≤ 4; by stable we mean solutions that do not leave the interval [0,1]. In our experiments we choose γ = 4 which puts the model in a chaotic regime. Note that in these experiments, the true model is deterministic. This seemingly simple non-linear map is non-invertible, with initialisation errors growing and saturating quickly. It has therefore been used before to test the performance of different configurations of the EnKF (Mitchell and Houtekamer, 2009). The true autocorrelation of the model error (computed offline as described in the introduction) has both positive and negative values, while decreasing in absolute magnitude as a function of the lag. Therefore, a 2-parameter memory is more appropriate for this case, although it increases the difficulty of the problem by having to estimate an extra parameter.

As a forecast model, we propose a stochastic linear model

where γ_g is a prescribed damping coefficient for the forecast model, with its value within 0 < γ_g < 1 to keep the model stable. The model error has the same properties as in the experiments where the real model is linear. In section 1, we described an off-line process to compute the statistics of this model error from a long-time collection of differences between perfect and imperfect forecasts. What we look for, instead, is to estimate the time-related parameters online, along with the state estimation.

3.3 Implementation details

Our experiments follow the fraternal (non-identical) twin set up. A real initial condition is prescribed and the true model is run for a long time. In the case of the true stochastic linear model shown in Eq. (2), a particular realisation of the model error is used for this true run. For this model we use α = 0.8 for all experiments. Synthetic observations are obtained at different times by applying Eq. (10) with r² = (0.01)². We generate synthetic observations every model time step, although we use different subsets in different experiments. The ensemble forecast is generated by selecting N_e initial conditions coming from N(0,b²), with b² = (0.1)². The model error realisations are generating using Eq. (2.1). This is, uncorrelated model errors are generated from N(0,q²), with q² = (0.1)². The time auto-correlation is inserted by computing the desired Cholesky factor L(θ) which depends in the parameter vector θ, performing the operation indicated in Eq. (2.1). For model error parameters we have the following:

• For the one-parameter model error, the prior distribution is ${\omega }^b~N({\mu }_{\omega }^b,{\sigma }_{\omega }^2)$ with ${\mu }_{\omega }^b=0.3$ and ${\sigma }_{\omega }^2={(0.5)}^2$.

• For the two-parameter model error, we consider the two parameters to be uncorrelated. Their marginal background distributions are: ${\varphi }^b~N({\mu }_{\varphi }^b,{\sigma }_{\varphi }^2)$ with ${\mu }_{\varphi }^b=0.3$ and ${\sigma }_{\varphi }^2={(0.2)}^2$, and $f^b~N({\mu }_f^b,{\sigma }_f^2)$ with ${\mu }_f^b=0.3$ and ${\sigma }_f^2={(0.2)}^2$.

Both the EnKS and IEnKS use stochastic implementations, (Van Leeuwen, 2020) i.e. one needs to add observational noise to observation values predicted by each ensemble member. In this way the real observation and the predicted observations are generated in the same way. In all our experiments, we use an assimilation window of τ = 20 time steps plus the initial condition. When we assimilate 1 observation per window, this is located at the end time. When we assimilate 20 observations per window, these correspond to all time steps except that of the initial conditions. In order to avoid the introduction of sampling errors –which would be difficult to disentangle from the other effects we study–, we use a relatively large ensemble with N_e = 200. This size is kept constant for all experiments.

For the IEnKS, we vary the number of iterations and the step length δ. We consider the IEnKS to have converged to a solution when the difference |u_i+1–u_i| < 0.01, which is approximately 1% of the range of the parameters.

4.1 Illustration of state estimation

In the first set of experiments we aim to illustrate the performance of the EnKS for state estimation. To this end, we use the 1-parameter memory setting and choose background values very close to the real parameter, both the linear and non-linear models. The results of these experiments are displayed in Figure 3. In all panels, the truth is shown with a black thick line, randomly selected background ensemble members in blue, and randomly selected analysis members in red. For the analysis, the thick red line corresponds to the mean. We do not show all members to avoid visual cluttering. For this experiment, we only study a single assimilation window.

Figure 3

The top row of this figure shows the results with the linear model. In this case, the EnKS is capable of recovering an analysis mean quite close to the truth. The availability of more observations makes an important difference, as noted in the results shown in the two columns (left shows 1 observation per window, right shows 20 observations per window). In this simple linear case, even with very sparse observations, the EnKS can still provide a fair estimation of the state.

The bottom row shows the results with the logistic map. In the case with a single observation at the end of the simulation window, the EnKS was unsuccessful, which is expected considering the length of the assimilation window, and that joint estimation is being performed. To obtain a better results, more information from the observations is needed. With dense observations, observing every time-step for instance, the EnKS provides a correct estimation of all the state values at all time steps.

4.2 Parameter estimation in the linear model

4.2.1 One-parameter estimation

Our first experiments with the linear model consist of trying to estimate the parameter ω via the EnKS, using the one-parameter autocorrelation shown in Eq. (37). The results are shown in Figure 4. This figure has four panels. Panels (a), (b) and (c) have the same format, the difference amongst the panels is the number of observations per assimilation window. In each panel, the black vertical line shows the true value of the parameter. Different pdfs are displayed, computed as histograms using the ensemble members. The blue line shows the background pdf for the estimated parameter, whereas the orange-red lines correspond to the analysis pdfs after different number of assimilation windows. We see that the DA system works well for one-parameter estimation, especially with multiple observations. In fact, even with a single observation at the end of each simulation window, the analysis pdf moves towards the true ω, and the variance of the posterior pdfs is smaller than that of the prior, and it gets smaller as the number of consecutive assimilation windows increases. The effect of using more observations in each simulation window is to accelerate the convergence of the posterior pdf towards a stationary distribution. Panel (d) illustrates more clearly the evolution of the analysis estimator for sequential assimilation windows (in the horizontal axis). The blue horizontal line shows the value of the background mean ${\overline{\omega }}^b\approx 1$, whereas the horizontal black line shows the true value of the parameter ω^r = 0.5. In this case, we only track the analysis mean (solid colour lines) and the analysis standard deviation (shaded colour areas). With one observation at the end of the assimilation window, the convergence of the analysis mean is slow, taking about 50 assimilation windows before the analysis mean and standard deviation stabilise. There is, however, a considerable bias with respect to the true parameter. Increasing the number of observations in the assimilation window accelerates the convergence of the mean estimator. The cases with 5 and 20 observations per assimilation window approach their final value between 10 and 20 assimilation windows. The final estimate with 5 observations is still slightly larger than the true parameter. The estimate with 20 observations satisfactorily converges toward the real parameter from below, although at the end of the 80 assimilation windows a small gap still remains, much larger then the estimated uncertainty in the ensemble.

Figure 4

We perform similar experiments with the IEnKS. Recall that in this case the number of iterations can influence the performance of the smoother. We analyse the effect of the number of iterations in Figure 5. This figure has two similar panels corresponding to 1 observation (left) and 20 observations (right) per window. The blue horizontal line denotes the background mean, and the black line the real parameter value. In each panel we plot several lines corresponding to analysis mean after different number of assimilation windows, as well as a shaded area corresponding to the analysis standard deviation. In the horizontal axis we have the number of iterations in each window. These iterations use a fixed step length δ. Perhaps the most important message of this figure is that the IEnKS results are independent of the number of assimilation windows since there is not a noticeable difference between the different lines. For one observation at the end of the assimilation window, the estimator has not fully converged after 20 iterations. It seems that at least 5 iterations per window are necessary for the IEnKS to show reasonable performance. In the case of 20 observations per window, it takes between 15 to 20 iterations for the IEnKS to converge to the true value of the parameter. Again, the number of assimilation windows does not matter. Based on these results we fix the number of iterations in the IEnKS to N_iter = 10 for the rest of our experiments. In practice, one aims to perform as few iterations as possible to minimise computational expenses. We implement the IEnKS with a different numbers of observations and simulation windows. These results shown in Figure 6. The overall format of this figure is the same as that of Figure 4. Panels (a) to (c) reveal that the convergence to the analysis pdf is faster with more observations per window, and the resulting consecutive analysis pdf’s have smaller variance. Panel (d) reveals interesting results. For the three observation frequencies, the major changes to the analysis mean occur in the first 20 assimilation windows. For the cases with 1 and 5 observations per window, the bias in the final analysis mean is reduced considerably with respect to the estimators obtained by the EnKS. With 20 observations per window, the estimation is remarkably accurate, with the analysis mean oscillating around the real parameter value after 50 windows and remaining there.

Figure 5

Figure 6

In this subsection we have shown that EnKS converges to the real parameter value but a small bias remains, while its uncertainty estimate tends to contain that real value. The results show significant improvement increasing in the number of observations and the EnKS benefits from more assimilation windows. Next step is to extend our experiments to a more complex model error setting with two unknown parameters.

4.2.2 Two-parameter estimation

For this part of the experiments, we implement the model error with formulated in Eq. (39), which contains both decaying and oscillatory parameters. We start by experimenting with the EnKS. This method, however, fails to find the correct values for both parameters even in the case with observations every time step. This case is shown in Figure 7. This figure shows the background (blue line) and analysis (orange-red lines) pdf’s for the oscillatory (left) and decaying (right) parameters. The vertical black lines correspond the true parameter values. The figure illustrates how the mean of the analysis pdf converges towards wrong values of the parameters, with noticeable variance reduction after successive assimilation windows. This behaviour suggests that the minimisation process in the IEnKS is converging to an incorrect local minimum. To explore the failure of the EnKS in this case, it is useful to display the cost function of the problem. This requires computing Eq. (32) using Eq. (39) and Eq. (41). For simplicity, we do this for the case of a single operation at the end of the assimilation window. Before explaining the result, we need to recall some aspects of this function:

• The input of this function is the vector u and the output is a scalar. The function maps 1 + τ + N_θ values into a single one.

• The background elements needed by the cost function are the background mean ${\mu }_u^b$ and and covariance D_u.

• The observation elements needed by the cost function are the actual observation y and its variance r².

Figure 7

In Figure 8 we start with a simple setting. We use the exact mean ${\mu }_u^b$ and covariance D_u, as well as the observational variance r² from section 3.3. Since we are not interested in the behaviour of the cost function with respect to the state variables, we set the state variables to a fixed value. and we only let the values θ = [f,φ]^T vary freely. We plot a cross-section of the cost-function in the two-dimensional parameter space. We choose x₀ = 0 and ε₁ = … = ε_τ = 0.1 –note that we have to choose a value different from zero, otherwise the cost function would be blind to the variation in model error parameters.

Figure 8

Figure 8 has nine panels. Each one of the columns correspond to a different length of the assimilation window (recall that the observation is taken at the end), and each one of the columns correspond to a different given observation. We already know that the resulting cost function is not quadratic, and that the deformation from a quadratic cost function grows as the number of time step grows, and as the difference between the given observation and the predicted observation is larger. This is exactly what we observe in the figure. Both the first and third columns, when the observations are y 4.0 respectively, show the largest deformation of the cost function and this grows for the rows from top to bottom. The middle panel (y = 0) sees little deformation from a paraboloid. For each panel, the blue dot represents the background mean values of the parameters, and the pink dot corresponds to the global minimum of the cost function. In all cases, the cost function is still concave and there is a unique minimum to be found.

The simple case shown before is not what the EnKS encounters in practice. First of all, we do not know the real background covariance matrix, so this comes from a sample estimator. Second, the cost function is applied to each ensemble member separately, and therefore the background mean is just the background value for each member. The final analysis estimator is the mean of the N_e estimated minima. We explore this in Figure 9. First of all, we fix given observation to y = –4. Each column represents a different ensemble member (3 members out of the 200 we used to compute background statistics), and each row represents a different length of the assimilation window. To produce these cost functions, the background values are random realisations of the distributions for the initial conditions, model errors and parameter values shown in section 3.3. Again, we are interested in the behaviour in the two-dimensional parameter space, so we fix all the state variables to the same constants as before. What we see is very different from Figure 8. Figure 9 shows cost functions which are not convex, in fact they have very intricate topographic features such as narrow elongated ridges and valleys, and even some local mimima and maxima. The global cost minimum is difficult to converge to in a single application of the EnKS. In fact, this does not happen. We see that as the assimilation window length increases, the complicated features of the cost function increase. This agrees with the unsuccessful results we had seen in Figure 7. Therefore, it is necessary to apply a number of smaller steps via iterations in an IEnKS.

Figure 9

We now apply the IEnKS to avoid getting stuck in a valley or a local minimum in the minimisation process. The results for the two-parameter estimation using IEnKS are shown in Figures 10 and 11. We apply the IEnKS with 10 iterations and a fixed step length δ = 0.3.

Figure 10

Figure 11

Figure 10 shows the ensemble background and analysis pdfs (after different number of assimilation windows) for the oscillatory (top row) and decaying (bottom row) parameters, for 1 observation (top row) and 20 observations (bottom row). This figure reveals that the IEnKS works fairly well in this case. Similar to the results with the single parameter estimation, the two-parameter estimation results improve with the increase in the number of observations and assimilation windows in the simulation period. The IEnKS estimation converges towards the correct values of both parameters. In the case of only 1 observation in the assimilation window, the pdf stops short of the real value (black line), while in the case of 20 observations the estimation is better. Figure 11 tracks the evolution of the analysis mean and standard deviation for the two parameters (one for each panel) as the number of assimilation windows increases, and for three total number of observations: 1, 5 and 20. The true and background mean values are shown with horizontal lines, black and blue correspondingly. Note the improvement in the estimation as the number of observations increase. Even for observations every time step, there is a bias in the final estimation of the frequency parameter, while the decay parameter is captured in the ensemble uncertainty. Compared with the results from the one-parameter estimation, the two-parameter estimation problem seems to be much more complicated, and requires more observations to make the estimation work.

In closing the experiments for this subsection, we want to see how the IEnKS works with different initial guess for the two parameters. Given the complicated shape of the cost function, it is conceivable that we may get stuck in local minimums, even when using this sophisticated iterative method. We divide the two-dimensional parameter space into four quadrants, and we choose starting points in each one of the quadrants. These results are shown in Figure 12. This figure shows results in the case of one observation (left) and 20 observation (right). In both cases, the true values of the parameters are at the centre of the quadrants (denoted with a black dot). We see that the position of the initial condition can have serious impact on the estimation results. When we only observe at the end of the window most initial conditions do not lead to a value close to the true value. In fact, we see that the DA system cannot distinguish between parameter values lying roughly on a straight line defined by the red dots. On the other hand, when we have observations at every time-step, the posterior seems to get fairly close to the true values of the parameters, but many local minima exist around the true values. The results show the importance of the initial guess of the parameters on the parameter estimation results.

Figure 12

After experimenting with the linear model, our next step is extending the experiments of two-parameter model error autocorrelation to the nonlinear model, the logistic map.

4.3 Parameter estimation in a non-linear model

The last experiments we perform are also the most challenging. In this case our true model is the logistic map with the damping coefficient γ = 3.99, and the proxy model is linear model with auto-correlated model error. The true model error autocorrelation in this case has both an oscillatory and decaying behaviour (this can be diagnosed offline as described in section 3), so we directly try two-parameter estimation. Since the EnKS failed on the two-parameter estimation even with the linear model, we decide to only apply the IEnKS with 10 iterations, over times windows of 10 time steps long.

In our first experiments we try to estimate f and φ using the IEnKS with a fixed iteration step length δ = 0.3 as used in the linear experiments. The results of these experiments are shown in Figure 13 for both parameters (rows) and for 1 and 20 observations in the window (columns). With 1 observation at the end of the window the smoother does not seem to converge, and with 20 observations convergence seems to occur, but to incorrect values of the parameters. Increasing the number of iterations pr window or the number of windows did not improve results.

Figure 13

Since the IEnKS with a fixed step length δ fails to estimate the parameters directly, we explored three further experimental settings. In the first we transform the estimation problem to a more linear setting by estimating lag-1 and lag-2 autocorrelations of the model error (labelled as AC(1) and AC(2) respectively), and then transform the results to f and φ as the following:

In the second setting we use a decaying step length δ while using the untransformed parameters. However, Figure 14 shows that neither method is successful, suggesting that the problem is highly nonlinear and a more careful tuning of the minimization is required.

Figure 14

In the third setting we both transform problem to first estimating the autocorrelations and applying a decaying step length. However, with only 1 observation in the assimilation window the IEnKS failed, no matter what we tried. For a fully observed system with 10 observations per assimilation window successful results can be achieved. We show results for reducing the step length δ by 2% after each iteration (i.e. c = 0.98 in Eq. (36)). From the two top plots in Figure 15, we can see that the results are similar to the successful results for the linear model shown in Figure 10. The posterior shows improvement with more simulation windows, and the variance decreases with more windows as well. Figures 15 (c) and (d) show that the posterior mean moves towards the right values for both parameters, with slowly decaying uncertainty.

Figure 15

Even though the minimization is successful the convergence is extremely slow, even in this simple zero-dimensional model. This illustrates that estimating parameters in model errors that are related to temporal correlations is a very hard problem.