# Hybrid Levenberg-Marquardt and weak constraint ensemble Kalman smoother method, Nonlinear Processes in Geophysics

Jan Mandel, El Houcine Bergou, Selime Gurol, Serge Gratton

► To cite this version:

Jan Mandel, El Houcine Bergou, Selime Gurol, Serge Gratton. Hybrid Levenberg-Marquardt and weak constraint ensemble Kalman smoother method, Nonlinear Processes in Geophysics. Nonlinear Processes in Geophysics, 2016, 23, pp.59-73. 10.5194/npgd-2-865-2015. hal-03165022

HAL Id: hal-03165022

https://hal.science/hal-03165022

Submitted on 11 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

This discussion paper is/has been under review for the journal *Nonlinear Processes in Geophysics* (NPG). Please refer to the corresponding final paper in NPG if available.

# Hybrid Levenberg–Marquardt and weak constraint ensemble Kalman smoother method

J. Mandel<fnref n="1" />, E. Bergou<fnref n="2" />,<fnref n="3" />, S. Gürol<fnref n="3" />, and S. Gratton<fnref n="3" />

Received: 21 April 2015 – Accepted: 23 April 2015 – Published: 26 May 2015

Correspondence to: J. Mandel (jan.mandel@gmail.com)

Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<footnote n="1">¹University of Colorado Denver, Denver, CO, USA</footnote>

<footnote n="2">²Institut de Mathématiques de Toulouse, Toulouse, France</footnote>

<footnote n="3">³INP-ENSEEIHT and CERFACS, Toulouse, France</footnote>
Abstract

We propose to use the ensemble Kalman smoother (EnKS) as the linear least squares solver in the Gauss–Newton method for the large nonlinear least squares in incremental 4DVAR. The ensemble approach is naturally parallel over the ensemble members and no tangent or adjoint operators are needed. Further, adding a regularization term results in replacing the Gauss–Newton method, which may diverge, by the Levenberg–Marquardt method, which is known to be convergent. The regularization is implemented efficiently as an additional observation in the EnKS. The method is illustrated on the Lorenz 63 and the two-level quasi-geostrophic model problems.

10

# 1 Introduction

Four dimensional variational data assimilation (4DVAR) is a dominant data assimilation method used in weather forecasting centers worldwide. 4DVAR attempts to reconcile model and data variationally, by solving a very large weighted nonlinear least squares problem. The unknown is a vector of system states over discrete points in time, when the data are given. The objective function minimized is the sum of the squares of the differences of the initial state from a known background state at the initial time and the differences of the values of observation operator and the data at every given time point. In the weak-constraint 4DVAR (Trémolet, 2007), considered here, the model error is accounted for by allowing the ending and starting states of the model at every given time point to be different, and adding to the objective function also the sums of the squares of those differences. The sums of the squares are weighted by the inverses of the appropriate error covariance matrices, and much of the work in the applications of 4DVAR goes into modeling those covariance matrices.

15

20

In the incremental approach (Courtier et al., 1994), the nonlinear least squares problem is solved iteratively by using a succession of linear least square solutions. The major cost in 4DVAR iterations is in evaluating the model, tangent and adjoint opera-

25

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>
tors, and solving large linear least squares. A significant software development effort is needed for the additional code to implement the tangent and adjoint operators to the model and the observation operators. Straightforward linearization, called the incremental approach (Courtier et al., 1994), leads to the Gauss–Newton method for nonlinear least squares (Bell, 1994; Tshimanga et al., 2008). However, Gauss–Newton iterations may not converge, not even locally. Finally, while the evaluation of the model operator is typically parallelized on modern computer architectures, there is a need to further parallelize the 4DVAR process itself.

The Kalman filter is a sequential Bayesian estimation of the gaussian state of a linear system at a sequence of discrete time points. At each of the time points, the use of the Bayes theorem results in an update of the state, represented by its mean and covariance. The Kalman smoother considers all states at all time points from the beginning to be a large composite state. Consequently, the Kalman smoother can be obtained from the Kalman filter by simply applying the same update as in the filter to the past states as well. However, historically, the focus was on efficient short recursions (Rauch et al., 1965; Strang and Borre, 1997), similar the sequential Kalman filter.

It is well known that weak constraint 4DVAR is equivalent to the Kalman smoother in the linear case. To apply the Kalman smoother in the nonlinear case, the problem needs to be linearized, leading to variants of the extended Kalman filter and the Gauss–Newton method. Use of the Kalman smoother to solve the linear least squares in the Gauss–Newton method is known as the iterated Kalman smoother, and considerable improvements can be obtained against running the Kalman smoother only once (Bell, 1994; Fisher et al., 2005).

The Kalman filter and smoother require maintaining the covariance of the state, which is not feasible for large systems, such as in numerical weather prediction. Hence, the ensemble Kalman filter (EnKF) and ensemble Kalman smoother (EnKS) (Evensen, 2009) use a Monte-Carlo approach for large systems, representing the state by an ensemble of simulations, and estimating the state covariance from the ensemble. The implementation of the EnKS in Stroud et al. (2010) uses the adjoint model explicitly,

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

with the short recursions and a forward and a backward pass, as in the KS. However, the implementations in Khare et al. (2008); Evensen (2009) do not depend on the adjoint model and simply apply EnKF algorithms to the composite state over multiple time points. We use the latter approach in the computations reported here.

In this paper, we propose to use the EnKS as a linear least squares solver in 4DVAR. The ensemble approach is naturally parallel over the ensemble members. The rest of the computational work is relatively cheap compared to the ensemble of simulations, and parallel dense linear algebra libraries can be used. The proposed approach uses finite differences from the ensemble, and no tangent or adjoint operators are needed. To stabilize the method and assure convergence, a Tikhonov regularization term is added to the linear least squares, and the Gauss–Newton method becomes the Levenberg–Marquardt method. The Tikhonov regularization is implemented within EnKS as a computationally cheap additional observation (Johns and Mandel, 2008). We call the resulting method EnKS-4DVAR. Theoretical convergence of the algorithm for large ensembles is studied in Bergou et al. (2014), where a rigorous proof is provided that the iterations of the EnKS-4DVAR method converge to those of the incremental 4DVAR for large ensembles in the $L^p$ norm, for any $p \in [1, \infty)$, in the large ensemble limit.

The EnKF has become a competitive method for data assimilation. Consequently, combinations of ensemble and variational approaches have become of considerable recent interest. Estimating the background covariance for 4DVAR from an ensemble was one of the first connections (Hamill and Snyder, 2000), and it is now standard and became operational (Wang, 2010). Gradient methods in the span of the ensemble for one analysis cycle (i.e., 3DVAR) include Zupanski (2005); Sakov et al. (2012) (with square root EnKF as a linear solver in Newton method), and Bocquet and Sakov (2012), who added regularization and use LETKF-like approach to minimize the nonlinear objective function over linear combinations of the ensemble. Bocquet and Sakov (2012) scale the ensemble to approximate the derivatives (tangent operators) as in Sakov et al. (2012). They call their approach the “bundle variant”, which is the same
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

as using finite differences to approximate derivatives. Here, we use a similar technique to approximate the derivatives. Liu et al. (2008, 2009); Liu and Xiao (2013) also minimize the (strong constraint) 4DVAR objective function over linear combinations of the ensemble by computations in the observation space. Their method, called Ens4DVAR, does not need tangent or adjoint operators also. Zhang et al. (2009) use a two-way connection between EnKF and 4DVAR, to obtain the covariance for 4DVAR, and 4DVAR to feed the mean analysis into EnKF. EnKF is operational at the National Centers for Environmental Prediction (NCEP) as part of its Global Forecast System Hybrid Variational Ensemble Data Assimilation System (GDAS), together with the Gridpoint Statistical Interpolation (GSI) variational data assimilation system (Developmental Testbed Center, 2015).

Additional work appeared after the first version of this paper was written (Mandel et al., 2013). Bocquet and Sakov (2014) extend the method of Bocquet and Sakov (2012) to 4DVAR and use finite difference approximations of the tangent operators, similarly as in Sakov et al. (2012) and here. However, Bocquet and Sakov (2014) nest the minimization loop for the 4DVAR objective function inside a square root version of the EnKS and minimize over the span of the ensemble, rather than nesting EnKS as a linear solver inside the 4DVAR minimization loop over the full state space as here. Their method is tied to the use of the sample covariance matrix of the state without localization of the covariance and to strong-constraint 4DVAR. The implementation of the EnKS in the computations reported here also uses sample covariance, and, consequently, the analysis ensemble consists of linear combinations of the forecast ensemble. However, limiting the EnKF to linear combinations only does not allow common approaches to localization (Sakov and Bertino, 2011).

The present approach is not tied to the use of sample covariance. Rather, it allows an arbitrary implementation of the EnKS to be used, which makes localization possible, e.g., by tapering the sample covariance (Furrer and Bengtsson, 2007) or replacing the sample covariance by its diagonal in a spectral space (Kasanický et al., 2015). Implementations of the EnKS with localization exist, e.g., Butala (2012), which is based
on the Bryson–Frazier version of the classical formulation the KS, with a forward and a backward pass. Ensemble methods for the solution of the 4DVAR nonlinear least squares problem in the weak constraint 4DVAR, or ensemble methods for this problem which allow localization, do not seem to have been developed before.

The paper is organized as follows. In Sect. 2, we review the formulation of 4DVAR. The EnKS for the incremental linearized squares problem is reviewed in Sect. 3. The new method without tangent operators is introduced in Sect. 4. The modifications for the regularization and the Levenberg–Marquardt method are presented in Sect. 5. Section 6 contains the results of the computational experiments, and Sect. 7 is the conclusion.

## 2 Incremental 4DVAR and the Gauss–Newton method

We want to estimate $x_0, \dots, x_k$, where $x_i$ is the state at time $i$, from the background state, $x_0 \approx x_b$, the model, $x_i \approx M_i(x_{i-1})$, and the observations $H_i(x_i) \approx y_i$, where $M_i$ is the model operator, and $H_i$ is the observation operator. Quantifying the uncertainty by covariances, with $x_0 \approx x_b$ taken as $(x_0 - x_b)^T B^{-1} (x_0 - x_b) \approx 0$, etc., we get the nonlinear least squares problem

$$ \|x_0 - x_b\|_{B^{-1}}^2 + \sum_{i=1}^{k} \|x_i - M_i(x_{i-1})\|_{Q_i^{-1}}^2 + \sum_{i=1}^{k} \|y_i - H_i(x_i)\|_{R_i^{-1}}^2 \rightarrow \min_{x_{0:k}} \quad (1) $$

called weak-constraint 4DVAR (Trémolet, 2007). Originally, in 4DVAR, $x_i = M_i(x_{i-1})$; the weak constraint $x_i \approx M_i(x_{i-1})$ accounts for model error.

The least squares problem (Eq. 1) is solved iteratively by linearization,

$$ \begin{aligned} M_i(x_{i-1} + \delta x_{i-1}) &\approx M_i(x_{i-1}) + M'_i(x_{i-1}) \delta x_{i-1}, \\ H_i(x_i + \delta x_i) &\approx H_i(x_i) + H'_i(x_i) \delta x_i. \end{aligned} $$

Title Page

Abstract

Introduction

Conclusions

References

Tables

Figures

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

Back

Close

Full Screen / Esc

Printer-friendly Version

Interactive Discussion

5

10

15

20
For vectors $u_i$, $i = 1, \dots, k$, denote the composite (column) vector

$$u_{0:k} = \begin{bmatrix} u_0 \\ \vdots \\ u_k \end{bmatrix}.$$

In each iteration $x_{0:k} \leftarrow x_{0:k} + \delta x_{0:k}$, one solves the auxiliary linear least squares problem for the increments $\delta x_{0:k}$,

$$\begin{aligned} & \|x_0 + \delta x_0 - x_b\|_{\mathrm{B}-1}^2 + \sum_{i=1}^k \|x_i + \delta x_i - M_i(x_{i-1}) - M'_i(x_{i-1})\delta x_{i-1}\|_{\mathrm{Q}_i^{-1}}^2 \\ & \qquad + \sum_{i=1}^k \|y_i - H_i(x_i) - H'_i(x_i)\delta x_i\|_{\mathrm{R}_i^{-1}}^2 \to \min_{\delta x_{0:k}}. \end{aligned} \quad (2)$$

This is the Gauss–Newton method (Bell, 1994; Tshimanga et al., 2008) for nonlinear squares, known in 4DVAR as the incremental approach (Courtier et al., 1994). Denote

$$\begin{aligned} z_{0:k} &= \delta x_{0:k}, & z_b &= x_b - x_0, & m_i &= M_i(x_{i-1}) - x_i, & d_i &= y_i - H_i(x_i), \\ \mathbf{M}_i &= M'_i(x_{i-1}), & \mathbf{H}_i &= H'_i(x_i), \end{aligned} \quad (3)$$

and write the auxiliary linear least squares problem (Eq. 2) as

$$\begin{aligned} & \|z_0 - z_b\|_{\mathrm{B}-1}^2 + \sum_{i=1}^k \|z_i - \mathbf{M}_i z_{i-1} - m_i\|_{\mathrm{Q}_i^{-1}}^2 + \sum_{i=1}^k \|d_i - \mathbf{H}_i z_i\|_{\mathrm{R}_i^{-1}}^2 \to \min_{z_{0:k}} \\ & \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \min_{\delta x_{0:k}}. \end{aligned} \quad (4)$$

The function minimized in Eq. (4) is the same as the one minimized in the Kalman smoother (Bell, 1994). The Gauss–Newton method with the Kalman smoother as the linear least squares solver is known as the iterated Kalman smoother, and considerable improvements can be obtained against running the Kalman smoother, applied to the linearized problem, only once (Bell, 1994; Fisher et al., 2005).

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>
### 3 Ensemble Kalman filter and smoother

We present the EnKF and EnKS algorithms, essentially following Evensen (2009), in a form suitable for our purposes. We start with a formulation of the EnKF (Algorithm 1), in a notation useful for the extension to EnKS. The notation $v^\ell \sim N(m, A)$ means that $v^\ell$ is sampled from $N(m, A)$ independently of anything else. The ensemble of states of the linearized model at time $i$, conditioned on data up to time $j$ (that is, with the data up to time $j$ already ingested), is denoted by $Z_{i|j}^N = [z_{i|j}^1, \dots, z_{i|j}^N] = [z_{i|j}^\ell]$, where the ensemble member index $\ell$ always runs over $\ell = 1, \dots, N$, and similarly for other ensembles. Assume for the moment that the observation operator $H_i$ is linear, that is, $H_i(u) = H_i u$.

**Algorithm 1 EnKF**

Initialize $z_{0|0}^\ell \sim N(z_b, B)$, $\ell = 1, \dots, N$

For $i = 1, \dots, k$, advance in time

$$z_{i|t-1}^\ell = M_i(z_{i-1|t-1}^\ell) + v_i^\ell, \quad v_i^\ell \sim N(0, Q_i), \qquad (5)$$

followed by the analysis step

$$z_{i|t}^\ell = z_{i|t-1}^\ell - P_i^N H_i^T (H_i P_i^N H_i^T + R_i)^{-1} (H_i (z_{i|t-1}^\ell) - d_i - w_i^\ell), \quad w_i^\ell \sim N(0, R_i), \qquad (6)$$

where $P_i^N$ is the sample covariance computed from the ensemble $Z_{i|t-1}^N$.

We write the matrices in Eq. (6) as

$$P_i^N H_i^T = \frac{1}{N-1} A(H_i A)^T, \quad H_i P_i^N H_i^T = \frac{1}{N-1} H_i A(H_i A)^T, \qquad (7)$$
where **A** is the matrix of anomalies of the ensemble $\mathbf{Z}_{i|i-1}^N$,

$$ \mathbf{A} = [a^1, \dots, a^N] = [\mathbf{z}_{i|i-1}^1 - \bar{\mathbf{z}}_{i|i-1}, \dots, \mathbf{z}_{i|i-1}^N - \bar{\mathbf{z}}_{i|i-1}], \quad \bar{\mathbf{z}}_{i|i-1} = \frac{1}{N} \sum_{j=1}^{N} \mathbf{z}_{j|i-1}. \qquad (8) $$

In particular, $\mathbf{H}_i$ is used here only in the matrix-vector multiplications

$$ \mathbf{H}_i a^l = \mathbf{H}_i (\mathbf{z}_{i|i-1}^l - \bar{\mathbf{z}}_{i|i-1}) = \mathbf{H}_i \mathbf{z}_{i|i-1}^l - \frac{1}{N} \sum_{j=1}^{N} \mathbf{H}_i \mathbf{z}_{j|i-1}^j. \qquad (9) $$

Equation (9) allows the use of a nonlinear observation operator $\mathcal{H}_i$, which only needs to be evaluated on the members of the ensemble. This technique is commonly used for nonlinear observation operators, e.g., Chen and Snyder (2007); Mandel et al. (2009). The replacement of $\mathbf{H}_i$ by $\mathcal{H}_i$ in Eq. (6) is straightforward.

From Eq. (7), it follows that the analysis ensemble $\mathbf{Z}_{i|i}^N$ consists of linear combinations of the forecast ensemble, which can be written as multiplying the forecast ensemble by a suitable transformation matrix $\mathbf{T}_i^N$,

$$ \mathbf{Z}_{i|i}^N = \mathbf{Z}_{i|i-1}^N \mathbf{T}_i^N, \quad \mathbf{T}_i^N \in \mathbb{R}^{N \times N}. \qquad (10) $$

The EnKS is obtained by applying the same analysis step (Eq. 6) as in the EnKF to the ensemble $\mathbf{Z}_{0:i|i-1}$ of composite states from time 0 to $i$, conditioned on data up to time $i-1$,

$$ \mathbf{Z}_{0:i|i-1}^N = \begin{bmatrix} \mathbf{z}_{0|i-1}^N \\ \vdots \\ \mathbf{z}_{i|i-1}^N \end{bmatrix}, $$
<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

in the place of $\mathbf{Z}_{i|j-1}$, with the observation matrix $\tilde{\mathbf{H}}_{0:i} = [0, \dots, \mathbf{H}_i]$. Then, Eq. (6) becomes

$$z_{0:i|j}^{\ell} = x_{0:i|j-1}^{N} - \mathbf{P}_{0:i,0:i}^{N} \tilde{\mathbf{H}}_{0:i}^{\mathrm{T}} (\tilde{\mathbf{H}}_{0:i} \mathbf{P}_{0:i,0:i} \tilde{\mathbf{H}}_{0:i}^{\mathrm{T}} + \mathbf{R}_i)^{-1} (\tilde{\mathbf{H}}_{0:i} z_{0:i|j-1}^{\ell} - \mathbf{d}_i^{\ell} - \mathbf{w}_i^{\ell}),$$

where $\mathbf{P}_{0:i,0:i}^{N}$ is the sample covariance matrix of $\mathbf{Z}_{0:i|j-1}^{N}$. Fortunately, the matrix-vector and matrix-matrix products can be simplified,

$$\tilde{\mathbf{H}}_{0:i} z_{0:i|j-1}^{\ell} = [0, \dots, 0, \mathbf{H}_i] z_{0:i|j-1}^{\ell} = \mathbf{H}_i z_{0:i|j-1}^{\ell} \quad (11)$$

$$\mathbf{P}_{0:i,0:i}^{N} \tilde{\mathbf{H}}_{0:i}^{\mathrm{T}} = \mathbf{P}_{0:i,i}^{N} \mathbf{H}_{i}^{\mathrm{T}}, \quad \tilde{\mathbf{H}}_{0:i} \mathbf{P}_{0:i,0:i} \tilde{\mathbf{H}}_{0:i}^{\mathrm{T}} = \mathbf{H}_{i} \mathbf{P}_{i,i}^{N} \mathbf{H}_{i}^{\mathrm{T}}, \quad (12)$$

which gives the analysis step (Eq. 14) in the EnKS (Algorithm 2).

**Algorithm 2 EnKS**

Initialize $z_{0|0}^{\ell} \sim N(z_b, \mathbf{B})$, $\ell = 1, \dots, N$.

For $i = 1, \dots, k$, advance in time,

$$z_{i|j-1}^{\ell} = \mathbf{M}_i z_{i-1|j-1}^{\ell} + \mathbf{m}_i + \mathbf{v}_i^{\ell}, \quad \mathbf{v}_i^{\ell} \sim N(0, \mathbf{Q}_i), \quad (13)$$

followed by the analysis step

$$z_{0:i|j}^{\ell} = x_{0:i|j-1}^{N} - \mathbf{P}_{0:i,i}^{N} \mathbf{H}_{i}^{\mathrm{T}} (\mathbf{H}_{i} \mathbf{P}_{i,i}^{N} \mathbf{H}_{i}^{\mathrm{T}} + \mathbf{R}_{i})^{-1} (\mathbf{H}_{i} z_{0:i|j-1}^{\ell} - \mathbf{d}_{i}^{\ell}), \quad \mathbf{d}_{i}^{\ell} \sim N(0, \mathbf{R}_{i}), \quad (14)$$

where $\tilde{\mathbf{H}}_{0:i} = [0, \dots, \mathbf{H}_i]$, and $\mathbf{P}_{i,i}^{N}$ is the sample covariance matrix of $\mathbf{Z}_{i|j-1}^{N}$.

The EnKS can be implemented in a straightforward manner by applying the same transformation as in the EnKF to the composite state from times 0 to $i$, $\mathbf{Z}_{0:i|j}^{N} = \mathbf{Z}_{0:i|j-1}^{N} \mathbf{T}_{i}^{N}$, where $\mathbf{T}_{i}^{N}$ is the transformation matrix in Eq. (10) (Brusdal et al., 2003, Eq. 20).
## 4 Nonlinear EnKS-4DVAR method

The formulation of the auxiliary least squares problem (Eq. 2) relies on the linearized (i.e., tangent) model operators $M_i$ and $H_i$ and their adjoints. The linearized model $M_i = M'_i(x_{i-1})$ occurs only in advancing the time as action on the ensemble members $\delta x^\ell = z^\ell$,

$$M_i z_{i-1}^\ell + m_i = M'_i(x_{i-1}) z_{i-1}^\ell + M_i(x_{i-1}) - x_i$$

Approximating by finite differences based at $x_{i-1}$ with step $\tau > 0$, we get

$$M_i z_{i-1}^\ell + m_i \approx \frac{M_i(x_{i-1} + \tau z_{i-1}^\ell) - M_i(x_{i-1})}{\tau} + M_i(x_{i-1}) - x_i. \quad (15)$$

Thus, advancing the linearized model in time requires $N + 1$ evaluations of $M_i$, at $x_{i-1}$ and $x_{i-1} + \tau \delta x_{i-1}^\ell$. The observation matrix $H_i$ occurs only in the action on the deviations of the ensemble of increments, $H_i(z_i^\ell - \bar{z}_i)$, $\bar{z}_i = \frac{1}{N} \sum_{j=1}^N z_j^\ell$. Approximating by finite differences based at $x_i$, with step $\tau > 0$, we have in the multiplication (Eq. 7) of the ensemble covariance,

$$H_i(z_i^\ell - \bar{z}_i) \approx \frac{H_i(x_i + \tau z_i^\ell) - H_i(x_i)}{\tau} - \frac{1}{N} \sum_{j=1}^{N} \frac{H_i(x_i + \tau z_j^\ell) - H_i(x_i)}{\tau} \quad (16)$$

Thus, evaluating the action of the linearized observation operator on the ensemble requires $N + 1$ evaluations of $H_i$, at $x_i$ and $x_i + \tau z_i^\ell$. We call the resulting method EnKS-4DVAR (Algorithm 3).
First, initialize

$$x_0 = x_b, \quad x_i = M_i(x_{i-1}), i = 1, \dots, k,$$

if not given already. One iteration (Eq. 2) of the incremental 4DVAR is then implemented as follows.

Given $x_0, \dots, x_k$, initialize $z_{0|0}^\ell \sim N(z_b, B)$ with $z_b = x_b - x_0$.

For $i = 1, \dots, k$, advance $z^\ell$ in time following Eq. (13), with the linearized operator approximated from Eq. (15),

$$z_{i|i-1}^{\ell} = \frac{M_i(\mathbf{x}_{i-1} + \tau z_{i-1|i-1}^{\ell}) - M_i(\mathbf{x}_{i-1})}{\tau} + M_i(\mathbf{x}_{i-1}) - \mathbf{x}_i + \mathbf{v}_i^{\ell}, \quad (17)$$

$\mathbf{v}_i^\ell \sim N(0, Q_i)$, followed by the analysis step (Eq. 14), with the multiplication by the matrix $H_i$ approximated by linearizing the observation operator from Eq. (16).

For $i = 1, \dots, k$, update $x_i \leftarrow x_i + \frac{1}{N} \sum_{\ell=1}^{N} z_{i|k}^\ell$.

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>
**Algorithm 4 Nonlinear EnKS**

Initialize $\mathbf{x}_{0|0}^k \sim N(\mathbf{x}_b, \mathbf{B})$.

For $i = 1, \dots, k$, advance in time

$$ \mathbf{x}_{i|t-1}^\ell = \mathcal{M}_i(\mathbf{x}_{i-1|t-1}^\ell) + \mathbf{v}_i^\ell, \quad \mathbf{v}_i^\ell \sim N(0, \mathbf{Q}_i) \quad (18) $$

followed by the analysis step

$$ \mathbf{X}_{0:i|t}^N = \mathbf{X}_{0:i|t-1}^N - \mathbf{A}\mathbf{G}^T(\mathbf{G}\mathbf{G}^T + \mathbf{R}_i)^{-1} \left[ \mathcal{H}_i(\mathbf{x}_{i|t-1}^\ell) - \mathbf{y}_i - \mathbf{w}_i^\ell \right]_{\ell=1,N'} \quad (19) $$
$$ \mathbf{w}_i^\ell \sim N(0, \mathbf{R}_i) $$

$$ \mathbf{A} = [a^1, \dots, a^N], \quad a_i^\ell = \mathbf{x}_{i|t-1}^\ell - \frac{1}{N} \sum_{j=1}^N \mathbf{x}_{i|j-1} $$

$$ \mathbf{G} = [g_i^1, \dots, g_i^N], \quad g_i^\ell = \mathcal{H}_i(\mathbf{x}_{i|t-1}^\ell) - \frac{1}{N} \sum_{j=1}^N \mathcal{H}_i(\mathbf{x}_{i|j-1}^j) \quad (20) $$

Note that for small $\tau$, the resulting method is asymptotically equivalent to the method with the derivatives (Bergou et al., 2014). Surprisingly, it turns out that in the case when $\tau = 1$, we recover the standard EnKS applied directly to the nonlinear problems, that is, with the linearized advance in time (Eq. 5) replaced by application of the original, nonlinear operator $\mathcal{M}_i$ (Algorithm 4, obtained by outting together Eq. (12) with Eqs. (7)–(9). We write the nonlinear EnKS as operating on the original ensemble of the states $\mathbf{X}^N = [\mathbf{x}^\ell]_{\ell=1}^N$ rather than on the increments $z^\ell = \delta\mathbf{x}^\ell$.) In particular, the incremental 4DVAR does not converge unless it is already at a stationary point, because each iteration delivers the same result, up to the randomness of the EnKS.

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>
**Theorem 1** If $\tau = 1$, then one step of EnKS-4DVAR (Algorithm 3) becomes the non-linear EnKS (Algorithm 4). In particular, in that case, the result of the step does not depend on the previous iterate.

Proof. Indeed, Eq. (17) becomes

$$ 
\begin{aligned}
z_{i|j-1}^{\ell} &= \frac{\mathcal{M}_i(\mathbf{x}_{i-1} + z_{i-1|j-1}^{\ell}) - \mathcal{M}_i(\mathbf{x}_{i-1})}{1} + \mathcal{M}_i(\mathbf{x}_{i-1}) - \mathbf{x}_i + \mathbf{v}_i^{\ell} \\
&= \mathcal{M}_i(\mathbf{x}_{i-1} + z_{i-1|j-1}^{\ell}) - \mathbf{x}_i + \mathbf{v}_i^{\ell},
\end{aligned}
$$

hence, $\mathbf{x}_i + z_{i|j-1}^{\ell} = \mathcal{M}_i(\mathbf{x}_{i-1} + z_{i-1|j-1}^{\ell}) + \mathbf{v}_i^{\ell}$, which is exactly the same as advancing the ensemble member $\ell$ following Eq. (18) with $\mathbf{x}_{i-1|j-1}^{\ell} = \mathbf{x}_{i-1} + z_{i-1|j-1}^{\ell}$. Similarly, Eq. (16) becomes with $\tau = 1$,

$$ 
\begin{aligned}
a_i^{\ell} &= \frac{\mathcal{H}_i(\mathbf{x}_i + z_{i|j-1}^{\ell}) - \mathcal{H}_i(\mathbf{x}_i)}{1} - \frac{1}{N} \sum_{j=1}^{N} \frac{\mathcal{H}_i(\mathbf{x}_i + z_{i|j-1}^{j}) - \mathcal{H}_i(\mathbf{x}_i)}{1} \quad (21)
\end{aligned}
$$

$$ = \mathcal{H}_i(\mathbf{x}_i + z_{i|j-1}^{\ell}) - \frac{1}{N} \sum_{j=1}^{N} \mathcal{H}_i(\mathbf{x}_i + z_{i|j-1}^{j}), \quad (22) $$

which is exactly the same as Eq. (20) with $\mathbf{x}_{i|j-1}^{\ell} = \mathbf{x}_i + z_{i|j-1}^{\ell}$. Finally, the innovation term in Eq. (14) becomes using Eq. (3),

$$ \mathbf{H}_i z_{i|j-1}^{\ell} - d_i = \frac{\mathcal{H}_i(\mathbf{x}_i + z_{i|j-1}^{\ell}) - \mathcal{H}_i(\mathbf{x}_i)}{1} - [\mathbf{y}_i - \mathcal{H}_i(\mathbf{x}_i)] = \mathcal{H}_i(\mathbf{x}_{i-1}^{\ell}) - \mathbf{y}_i, $$

which is exactly the same as in Eq. (19). □
## 5 Tikhonov regularization and the Levenberg-Marquardt method

The Gauss-Newton method may diverge, but convergence to a stationary point of (Eq. 1) can be recovered by a control of the step $\delta x$. Adding a constraint of the form $||\delta x_i|| \le \varepsilon$ leads to globally convergent trust region methods (Gratton et al., 2013). Here, we add $\delta x_i$ in a Tikhonov regularization term of the form $\gamma||\delta x_i||_{S_i^{-1}}^2$, which controls the step size as well as rotates the step direction towards the steepest descent, and obtain the Levenberg-Marquardt method (Levenberg, 1944; Marquardt, 1963) $x_{0:k} \leftarrow x_{0:k} + \delta x_{0:k}$, where

$$ ||\delta x_0 - z_b||_{B^{-1}}^2 + \sum_{i=1}^k ||\delta x_i - M_i \delta x_{i-1} - m_i||_{Q_i^{-1}}^2 + \sum_{i=1}^k ||d_i - H_i \delta x_i||_{R_i^{-1}}^2 + \gamma \sum_{i=0}^k ||\delta x_i||_{S_i^{-1}}^2 \rightarrow \min_{\delta x_{0:k}} \quad (23) $$

Under suitable technical assumptions, the Levenberg-Marquardt method is guaranteed to converge globally if the regularization parameter $\gamma \ge 0$ is large enough (Gill and Murray, 1978; Osborne, 1976). Estimates for the convergence of the Levenberg-Marquardt method in the case when the linear system is solved only approximately exist (Wright and Holt, 1985).

Similarly as in Johns and Mandel (2008), we interpret the regularization terms $\gamma||\delta x_i||_{S_i^{-1}}^2$ in Eq. (23) as arising from additional independent observations $\delta x_i \sim N(0, \gamma^{-1}S_i)$, which can be assimilated separately, resulting in a mathematically equivalent but often more efficient two-stage method – simply run the EnKF analysis (Eq. 6) twice. With the choice of $S_i$ as identity or, more generally a diagonal matrix, the implementation of these large observations can be made efficient (Mandel et al., 2009).

Note that unlike in Johns and Mandel (2008), where the regularization was applied to a nonlinear problem and thus the sequential data assimilation was only approximate, here the EnKS is run on the auxiliary linearized problem (Eq. 23), so all distributions

<figure><img src="image_44.png" /></figure>

<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>
are gaussian and the equivalence of solving Eq. (23) at once and assimilating the observations sequentially is statistically exact.

# 6 Computational results

In this section, we investigate the performance of EnKS-4DVAR method, described in this paper, by solving the nonlinear least-squares problem (Eq. 1) in which the dynamical models are chosen either the Lorenz 63 system (Lorenz, 1963) or the two-level quasi-geostrophic model (Fandry and Leslie, 1984).

We first consider experiments where the regularisation is not necessary to guarantee the convergence (i.e., $γ = 0$). Lorenz 63 equations are used as a forecast model for these experiments. Section 6.1 describes the Lorenz 63 model and presents numerical results on the convergence. Using the same model, in Sect. 6.2, we investigate the impact of the finite differences parameter $τ$, used to approximate the derivatives of the model and observation operators, along the iterations.

Experiments where the regularisation is necessary to guarantee the convergence are shown in Sect. 6.3, and we analyse the impact of the regularisation parameter $γ$ on the application to the two-level quasi-geostrophic model.

Note that for the experiments presented here, we do not use localization, hence we choose large ensemble sizes. In all experiments, the regularization covariance $S_j = I$.

## 6.1 Numerical experiments using Lorenz 63 model

The Lorenz 63 equations (Lorenz, 1963) are given by the nonlinear system

$$ \frac{dx}{dt} = -\sigma(x - y), \quad \frac{dy}{dt} = \rho x - y - xz, \quad \frac{dz}{dt} = xy - \beta z, $$

where $x = x(t)$, $y = y(t)$, $z = z(t)$ and $σ, ρ, β$ are parameters, whose values are chosen as 10, 28 and 8/3 respectively for the experiments described in this paper. These

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>
values result in a chaotic behaviour with two regimes as illustrated in Fig. 1. This figure shows the Lorenz attractor, which has two lobes connected near the origin, and the trajectories of the system in this saddle region are particularly sensitive to perturbations. Hence, slight perturbations can alter the subsequent path from one lobe to the other.

The system is discretized using the fourth-order Runge-Kutta method. The state at time t is denoted by $X_t = [x(t), y(t), z(t)]^T$, $X_t \in \mathbb{R}^3$.

To evaluate the performance of EnKS-4DVAR method (Algorithm 3), we will test it using the classical twin experiment technique, which consists on fixing an initial true state, denoted by truth₀, and then integrating the initial truth in time using the model to obtain the true state truthᵢ = M(truthᵢ₋₁) at each time i. We then build the data yᵢ by applying the observation operator Hᵢ to the truth at time i and by adding a Gaussian perturbation N(0, Rᵢ). Similarly, the background xᵣ is sampled from the Gaussian distribution with the mean truth₀ and the covariance matrix B. Then, we try to recover the truth using the observations and the background.

We perform numerical experiments without model error. The initial truth is set to truth₀ = [1, 1, 1]ᵀ and the background covariance is chosen as the identity matrix of order three, i.e. **B** = **I**₃. The model is advanced in time by using the Runge-Kutta method with a time step of 0.1 time unit. The time window length is k = 50 time steps (5 time units). The observation operator is defined as $H_i(x, y, z) = (x^2, y^2, z^2)$. At each time i, the observations are constructed as follows: $y_i = H_i(truth_i) + v_i$, where vᵢ is sampled from N(0, R) with R = **I**₃. Observations are taken for each time step (i = 1, ..., 50). The ensemble size is fixed to N = 100.

Figure 2 shows the estimator of the state vector $X_i$, i = 1, ..., 10, for the first five iterations of Algorithm 3. Figure 3 shows the root square error (RSE) for the same iterations shown in Fig. 2. RSE is defined as

$$RSE_i^{(j)} = \sqrt{\frac{1}{n}(\text{truth}_i - \mathbf{x}_i^{(j)})^T(\text{truth}_i - \mathbf{x}_i^{(j)})}, \quad j = 1, \dots, 5, \quad (24)$$
where truth<sub>i</sub> is the true vector state at time i, x<sub>i</sub><sup>(j)</sup> is the jth iterate of Algorithm 3 at time i and n is the length of x<sub>i</sub>. Table 1 shows the root mean square error (RMSE) for each iterate given by

$$
\text{RMSE}^{(j)} = \frac{1}{k} \sum_{i=0}^{k} \text{RSE}_{i}^{(j)}, \quad j = 1, \dots, 5, \qquad (25)
$$

where *k* is the number of time steps.

From Table 1 and Figs. 2 and 3, it can be seen that the iterates of Algorithm 3 con-
verge to the solution (without using a regularization). For these experiments, we ob-
serve that RMSE is reduced significantly in five iterations. Note that the error does not
converge to zero, because of the approximation and variability inherent in the ensemble
approach.

6.2 The impact of the finite difference parameter

In this section, we investigate the influence of the finite differences parameter $\tau$ used to approximate the derivatives of the model and observation operators. We use the same experimental set-up as described in the previous section. The numerical results are based on 30 runs of Algorithm 3 with eight iterations for Lorenz 63 problem, with the following choices for the parameter $\tau$: $1, 10^{-1}, 10^{-2}, 10^{-3}, 10^{-4}, 10^{-5}$ and $10^{-6}$.

The mean of the objective function values are shown in Table 2. The box plots of the objective function for the first four iterations are shown in Fig. 4 and for the last four iterations are shown in Fig. 5.

These figures and table show the impact of the parameter $\tau$ on the objective function minimization along the iterations. For $\tau = 1$ (when we use the classical non linear EnKS), the results are almost the same after the first iteration, in this case performing more iterations do not improve the results. However, when $\tau \le 10^{-1}$ the objective function is decreasing along iterations. For $\tau = 10^{-1}$, more iterations are needed to reduce the objective function significantly. In the case of small values of $\tau$, for instance

<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /></figure>
$\tau \le 10^{-2}$, few iterations are enough to reduce the objective function significantly. When $\tau = 10^{-3}$, the results are slightly different from the results with smaller values of $\tau$.

As a conclusion for these experiments, it is better to choose $\tau \le 10^{-3}$, such that the results will be less sensitive to the value of $\tau$. Note that this value is problem dependent. In practice, to avoid divergence due to the finite difference approximations, it is better to choose $\tau$ as small as possible, and since the computers use finite-precision arithmetic, we need to be careful to the effects of computer rounding.

We can also observe that for the first iteration, the best decrease in objective function is obtained when $\tau = 1$, and the worst decrease is obtained for $\tau = 10^{-6}$. Moreover, for the first four iterations the bigger $\tau$ is, the better results will be, but for the last four iterations, the smaller $\tau$ is, the better results will be. Hence, an adaptive $\tau$ can be a better choice than a fixed $\tau$ over iterations. For instance, for these experiments we can start with $\tau = 1$, then decrease its value along iterations. Exploration of the best strategy to choose $\tau$ over iterations will be studied in the future works.

## 6.3 Numerical tests using a two-layer Quasi Geostrophic model (QG)

The EnKS-4DVAR algorithm has been implemented into Object Oriented Prediction System (OOPS) (Trémolet, 2013), which is a data assimilation framework developed by European Centre for Medium-Range Weather Forecasts (ECMWF). Numerical experiments are performed by using the simple two-layer quasi-geostrophic model of OOPS platform. The details for the model and the data assimilation system are given in Sects. 6.3.1 and 6.3.2 respectively. Numerical experiments are performed to solve the weak-constraint data assimilation problem (Eq. 1) by using EnKS-4DVAR with regularization. Numerical results are presented in Sect. 6.3.3.

### 6.3.1 A two-layer quasi-geostrophic model

The two-layer quasi-geostrophic channel model is widely used in theoretical atmospheric studies, since it is simple enough for numerical calculations and it adequately

<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /></figure>

<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>
captures an important aspect of large-scale dynamics in the atmosphere: in the horizontal direction of the atmospheric flow, the Coriolis force caused by the rotation of the Earth, and the pressure gradient force are in approximate balance.

The two-layer quasi-geostrophic model equations are based on the non-dimensional quasi-geostrophic potential vorticity, whose evolution represents large scale circulations of the atmosphere. The quasi-geostrophic potential vorticity on the first (upper) and second (lower) layers can be written respectively as

$$q_1 = \nabla^2 \psi_1 - \frac{f_0^2 L^2}{g'H_1} (\psi_1 - \psi_2) + \beta y, \quad q_2 = \nabla^2 \psi_2 - \frac{f_0^2 L^2}{g'H_2} (\psi_2 - \psi_1) + \beta y + R_s, \quad (26)$$

where $\psi$ is the stream function, $\nabla^2$ is the two-dimensional Laplacian, $R_s$ represents orography or heating, $\beta$ is the (non-dimensionalised) northward variation of the Coriolis parameter at the fixed latitude $y$, $f_0$ is the Coriolis parameter at the southern boundary of the domain. $L$ is the typical length scale of the motion we wish to describe, $H_1$ and $H_2$ are the depths of the two layers, $g' = g\Delta\theta/\bar{\theta}$ is the reduced gravity where $\bar{\theta}$ is the mean potential temperature, and $\Delta\theta$ is the difference in potential temperature across the layer interface. Details of the derivation of these non-dimensional equations can be found in Fandry and Leslie (1984); Pedlosky (1979).

Potential vorticity in each layer is conserved and thus is described by

$$\frac{D_i q_i}{Dt} = 0, \quad i = 1, 2. \qquad (27)$$

where $D_i/Dt$, is the total derivative, defined by

$$\frac{D_i}{Dt} = \frac{\partial}{\partial t} + u_i \frac{\partial}{\partial x} + v_i \frac{\partial}{\partial y} \qquad (28)$$

and

$$u_i = -\frac{\partial \psi_i}{\partial y}, \quad v_i = \frac{\partial \psi_i}{\partial x}, \qquad (29)$$

<figure><img src="image_60.png" /></figure>

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /></figure>
are the horizontal velocity components in each layer. Therefore, the potential vorticity at each time step is determined by using the conservation of potential vorticity given by Eq. (27). In this process, time stepping consists of a simple first order semi-Lagrangian advection of potential vorticity.

Given the potential vorticity at a fixed time, Eq. (26) can be solved for the stream function at each gridpoint and then the velocity fields obtained through Eq. (29). The equations are solved by using periodic boundary conditions in the west–east direction and Dirichlet boundary condition in the north–south direction. For the experiments in this paper, we choose $L = 10^6$ m, $H_1 = 6000$ m, $H_2 = 4000$ m, $f_0 = 10^{-4} \text{s}^{-1}$, $\beta = 1.5$. For more details on the model and its solution, we refer the reader to Fisher et al. (2011).

The domain for the experiments is 12 000 km by 6300 km for both layers. The horizontal discretization consists of 40 × 20 points, so that the east–west and the north–south resolution is approximately 300 km. The dimension of the state vector of the model is then 1600. Note that the state vector is defined only in terms of the stream function.

### 6.3.2 Experimental setup

The performance of EnKS-4DVAR with regularization is analyzed by using twin experiments (Sect. 6.1).

The truth is generated from a model with layer depths of $D_1 = 6000$ m and $D_2 = 4000$ m, and the time step is set to 300 s, whereas the assimilating model has layer depths of $D_1 = 5500$ m and $D_2 = 4500$ m, and the time step is set to 3600 s. These differences in the layer depths and the time step provide a source of model error.

For all the experiments presented here, observations of non-dimensional stream function, vector wind and wind speed were taken from a truth of the model at 100 points randomly distributed over both levels. Observations were taken every 12 hours. We note that the number of observations is much smaller than the dimension of the state vector. Observation errors were assumed to be independent from each others and uncorrelated in time. The standard deviations (SD) were chosen to be equal to 0.4 for stream function observation error, 0.6 for vector wind and 1.2 for wind speed.

5

10

15

20

25

Title Page

Abstract

Introduction

Conclusions

References

Tables

Figures

<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /></figure>

Back

Close

Full Screen / Esc

Printer-friendly Version

Interactive Discussion

<figure><img src="image_68.png" /></figure>
The observation operator is the bi-linear interpolation of the model fields to horizontal observation locations.

The background error covariance matrix (B matrix) and the model error covariances (matrices $Q_k$) used in these experiments correspond to vertical and horizontal correlations. The vertical and horizontal structures are assumed to be separable and the vertical correlation function is assumed to be constant over the grid. In the horizontal plane covariance matrices correspond to isotropic, homogeneous correlations of stream function with Gaussian spatial structure obtained from a Fast Fourier Transform approach (Dietrich and Newsam, 1997; Nowak et al., 2003). For the background covariance matrix B, the SD and the horizontal correlation length scale in this experiments was set to 0.8 and $10^6$ m respectively. The vertical correlation function value was taken as 0.2. For the model error covariance matrices $Q_j$, the SD and the horizontal correlation length scale was set to 0.2 and $2 \times 10^6$ m respectively. The vertical correlation function value was taken as 0.5.

The window length is set to 10 days, with two sub-windows of 5 days ($k = 2$). The ensemble size is chosen to be $N = 3000$.

### 6.3.3 Numerical results

Figure 6 shows the objective function values along iterations, when using incremental 4DVAR method. The objective function is oscillating along the iterations, therefore incremental 4DVAR method without regularization is diverging.

Figures 7 and 8 show the objective function values along iterations for eight different choices of regularization parameter $\gamma$. It can be seen from these two figures that when $\gamma = 0$, Algorithm 3 is diverging as expected (since we do not use regularization and we only approximate the linearized subproblem using ensembles). For small values of $\gamma$ (for instance $\gamma \le 10^{-1}$), the objective function is not monotonically decreasing, hence Algorithm 3 is still diverging even if we use the regularization. Therefore, small values of $\gamma$ can not guarantee the convergence. For large values of $\gamma$ (for instance $\gamma \ge 10$),

<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

<figure><img src="image_71.png" /></figure>

<figure><img src="image_72.png" /></figure>
we can observe the decrease on the objective function along iterations. Moreover, the fastest decrease on the objective function is obtained for $γ = 10$.

In conclusion, when the regularization is used, the choice of the regularization parameter $γ$ is crucial to ensure the convergence. For instance, for small values of $γ$, the method can still diverge, and for large values of $γ$, the objective function decreases, but slowly (and many iterations may be needed to attain some predefined decrease). Therefore the regularization parameter should be neither “very small” nor “very large”. An adaptive $γ$ over iterations can be a better compromise which will be explored in future studies.

## 7 Conclusions

We have proposed a stochastic solver for the incremental 4DVAR weak constraint method. The regularization term added to the Gauss–Newton method, resulting in a globally convergent Levenberg–Marquardt method, maintains the structure of the linearized least squares subproblem, enabling us to use ensemble Kalman smoother as linear solver while simultaneously controlling the convergence. We have formulated the EnKS-4DVAR method (Algorithm 3) and have shown that it is capable of handling strongly nonlinear problems. We have demonstrated that the randomness of the EnKS version used (with perturbed data) eventually limits the convergence to a minimum, bit a sufficiently large decrease of the objective function can be achieved for successful data assimilation. On the contrary, we suspect that the randomization may help to increase the supply of the search directions over the iterations, as opposed to deterministic methods locked into one low-dimensional subspace, such as the span of a one given ensemble.

We have numerically illustrated the new method on the Lorenz 63 model and the two-level quasi-geostrophic model. We have analyzed the impact of the finite differences parameter $τ$ used to approximate the derivatives of the model and observation operators. We have shown that for $τ = 1$, the iterates obtained from EnKS-4DVAR are

Title Page

Abstract

Introduction

Conclusions

References

Tables

Figures

<figure><img src="image_73.png" /></figure>

<figure><img src="image_74.png" /></figure>

<figure><img src="image_75.png" /></figure>

<figure><img src="image_76.png" /></figure>

Back

Close

Full Screen / Esc

Printer-friendly Version

Interactive Discussion

<figure><img src="image_77.png" /></figure>
<figure><img src="image_78.png" /></figure>

<figure><img src="image_79.png" /></figure>

<figure><img src="image_80.png" /></figure>

<figure><img src="image_81.png" /></figure>

equivalent to those of obtained from the standard nonlinear EnKS. Based on computational experiments, we have concluded that in the first iteration, it is better to use the classical ensemble Kalman smoother as proposed in Evensen and van Leeuwen (2000) (i.e., $τ = 1$, Algorithm 4 here), and then to decrease $τ$ over the iterations.

For the second part of the experiments, we have shown the performance of the EnKS-4DVAR method with regularization on the two-level quasi-geostropic problem, one of the standard model problems for atmospheric circulation. We have observed that the incremental 4DVAR method is not converging for a long time window length, and that the regularization is necessary to guarantee convergence. We have concluded that the choice of the regularization parameter is crucial to ensure the convergence and different choices of this parameter can change the rate of decrease in the objective function. As a summary, an adaptive regularization parameter can be a better compromise to achieve the approximate solution in a reasonable number of iterations.

The choice of the parameters used in our approach is of crucial importance for the computational cost of the algorithm, for instance the number of iterations to obtain some desired reduction. The exploration in more detail of the best strategies to adapt these parameters course of the iterations will be studied elsewhere.

The base method, used in the computational experiments here, is using sample covariance. However, there is nothing to prevent the use of more sophisticated variances of EnKS with localization and the covariance inflation, and square root filters instead of EnKS with data perturbation. These issues, as well as, the performance on larger and realistic problems, will be studied elsewhere.

**Acknowledgements.** This research was partially supported by the Fondation STAE project ADTAO, the Czech Science Foundation under the grant GA13-34856S, and the US National Science Foundation under the grant DMS-1216481. A part of this work was done when Mandel was visiting INP-ENSEEIHT and CERFACS, and when Bergou and Gratton were visiting the University of Colorado Denver.
## References

Bell, B.: The iterated Kalman smoother as a Gauss-Newton method, SIAM J. Optimiz., 4, 626-636, doi:10.1137/0804035, 1994. 867, 871

Bergou, E., Gratton, S., and Mandel, J.: On the convergence of a non-linear ensemble Kalman smoother, arXiv:1411.4608, submitted, 2014. 868, 877

Bocquet, M. and Sakov, P.: Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems, Nonlin. Processes Geophys., 19, 383-399, doi:10.5194/npg-19-383-2012, 2012. 868, 869

Bocquet, M. and Sakov, P.: An iterative ensemble Kalman smoother, Q. J. Roy. Meteor. Soc., 140, 1521-1535, doi:10.1002/qj.2236, 2014. 869

Brusdal, K., Brankart, J. M., Halberstadt, G., Evensen, G., Brasseur, P., van Leeuwen, P. J., Dombrowsky, E., and Verron, J.: A demonstration of ensemble based assimilation methods with a layered OGCM from the perspective of operational ocean forecasting systems, J. Marine Syst., 40-41, 253-289, doi:10.1016/S0924-7963(03)00021-6, 2003. 874

Butala, M. D.: A localized ensemble Kalman smoother, in: 2012 IEEE Statistical Signal Processing Workshop (SSP), August 5-8, 2012, Ann Arbor, MI, USA, 21-24, IEEE, doi:10.1109/SSP.2012.6319665, 2012. 869

Chen, Y. and Snyder, C.: Assimilating vortex position with an ensemble Kalman filter, Mon. Weather Rev., 135, 1828-1845, doi:10.1175/MWR3351.1, 2007. 873

Courtier, P., Thépaut, J.-N., and Hollingsworth, A.: A strategy for operational implementation of 4D-Var, using an incremental approach, Q. J. Roy. Meteor. Soc., 120, 1367-1387, doi:10.1002/qj.49712051912, 1994. 866, 867, 871

Developmental Testbed Center: NOAA Ensemble Kalman Filter Beta Release v1.0, available at: http://www.dtcenter.org/com-GSI/users/docs (last access: March 2015), 2015. 869

Dietrich, C. R. and Newsam, G. N.: Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix, SIAM J. Sci. Comput., 18, 1088-1107, 1997. 886

Evensen, G.: Data Assimilation: the Ensemble Kalman Filter, 2nd edn., Springer, Berlin, 307 pp., doi:10.1007/978-3-642-03711-5, 2009. 867, 868, 872

Evensen, G. and van Leeuwen, P. J.: An ensemble Kalman smoother for non-linear dynamics, Mon. Weather Rev., 128, 1852-1867, doi:10.1175/1520-0493(2000)128:1852:AEKSFN>2.0.CO;2, 2000. 888

<figure><img src="image_82.png" /></figure>

<figure><img src="image_83.png" /></figure>

<figure><img src="image_84.png" /></figure>

<figure><img src="image_85.png" /></figure>
Fandry, C. and Leslie, L.: A two-layer quasi-geostrophic model of summer trough formation in the Australian subtropical easterlies, J. Atmos. Sci., 41, 807–817, 1984. 880, 884

Fisher, M., Leutbecher, M., and Kelly, G. A.: On the equivalence between Kalman smoothing and weak-constraint four-dimensional variational data assimilation, Q. J. Roy. Meteor. Soc., 131, 3235–3246, doi:10.1256/qj.04.142, 2005. 867, 871

Fisher, M., Trémolet, Y., Auvinen, H., Tan, D., and Poli, P.: Weak-Constraint and Long-Window 4D-Var, Tech. rep., European Centre for Medium-Range Weather Forecasts, ECMWF Technical Memorandum 655, Reading, UK, 2011. 885

Furrer, R. and Bengtsson, T.: Estimation of high-dimensional prior and posterior covariance matrices in Kalman filter variants, J. Multivariate Anal., 98, 227–255, doi:10.1016/j.jmva.2006.08.003, 2007. 869

Gill, P. E. and Murray, W.: Algorithms for the solution of the nonlinear least-squares problem, SIAM J. Numer. Anal., 15, 977–992, doi:10.1137/0715063, 1978. 879

Gratton, S., Gürol, S., and Toint, P.: Preconditioning and globalizing conjugate gradients in dual space for quadratically penalized nonlinear-least squares problems, Comput. Optim. Appl., 54, 1–25, doi:10.1007/s10589-012-9478-7, 2013. 879

Hamill, T. M. and Snyder, C.: A hybrid ensemble Kalman filter–3D variational analysis scheme, Mon. Weather Rev., 128, 2905–2919, doi:10.1175/1520-0493(2000)128<2905:AHEKFV>2.0.CO;2, 2000. 868

Johns, C. J. and Mandel, J.: A two-stage ensemble Kalman filter for smooth data assimilation, Environ. Ecol. Stat., 15, 101–110, doi:10.1007/s10651-007-0033-0, 2008. 868, 879

Kasanický, I., Mandel, J., and Vejmelka, M.: Spectral diagonal ensemble Kalman filters, Nonlin. Processes Geophys. Discuss., 2, 115–143, doi:10.5194/npgd-2-115-2015, 2015. 869

Khare, S. P., Anderson, J. L., Hoar, T. J., and Nychka, D.: An investigation into the application of an ensemble Kalman smoother to high-dimensional geophysical systems, Tellus A, 60, 97–112, doi:10.1111/j.1600-0870.2007.00281.x, 2008. 868

Levenberg, K.: A method for the solution of certain non-linear problems in least squares, Q. Appl. Math., 2, 164–168, 1944. 879

Liu, C. and Xiao, Q.: An ensemble-based four-dimensional variational data assimilation scheme. Part III: Antarctic applications with advanced research WRF using real data, Mon. Weather Rev., 141, 2721–2739, doi:10.1175/MWR-D-12-00130.1, 2013. 869

5

10

15

20

25

30

<figure><img src="image_86.png" /></figure>

<figure><img src="image_87.png" /></figure>

<figure><img src="image_88.png" /></figure>

<figure><img src="image_89.png" /></figure>
<figure><img src="image_90.png" /></figure>

<figure><img src="image_91.png" /></figure>

<figure><img src="image_92.png" /></figure>

<figure><img src="image_93.png" /></figure>

<figure><img src="image_94.png" /></figure>

Liu, C., Xiao, Q., and Wang, B.: An ensemble-based four-dimensional variational data assimilation scheme. Part I: Technical formulation and preliminary test, Mon. Weather Rev., 136, 3363–3373, doi:10.1175/2008MWR2312.1, 2008. 869

5

Liu, C., Xiao, Q., and Wang, B.: An ensemble-based four-dimensional variational data assimilation scheme. Part II: Observing system simulation experiments with Advanced Research WRF (ARW), Mon. Weather Rev., 137, 1687–1704, doi:10.1175/2008MWR2699.1, 2009. 869

10

Lorenz, E. N.: Deterministic nonperiodic flow, J. Atmos. Sci., 20, 130–141, doi:10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2, 1963. 880

Mandel, J., Beezley, J. D., Coen, J. L., and Kim, M.: Data assimilation for wildland fires: ensemble Kalman filters in coupled atmosphere–surface models, IEEE Contr. Syst. Mag., 29, 47–65, doi:10.1109/MCS.2009.932224, 2009. 873, 879

15

Mandel, J., Bergou, E., and Gratton, S.: 4D-VAR by ensemble Kalman smoother, arxiv:1304.5271, 2013. 869

Marquardt, D. W.: An algorithm for least-squares estimation of nonlinear parameters, J. Soc. Ind. Appl. Math., 11, 431–441, doi:10.1137/0111030, 1963. 879

Nowak, W., Tenkleve, S., and Cirpka, O.: Efficient computation of linearized cross-covariance and auto-covariance matrices of interdependent quantities, Math. Geol., 35, 53–66, 2003. 886

20

Osborne, M. R.: Nonlinear least squares – the Levenberg algorithm revisited, J. Aust. Math. Soc. B, 19, 343–357, doi:10.1017/S033427000000120X, 1976. 879

Pedlosky, J.: Geophysical Fluid Dynamics, Springer, New York Heidelberg Berlin, 715 pp., 1979. 884

25

Rauch, H. E., Tung, F., and Striebel, C. T.: Maximum likelihood estimates of linear dynamic systems, AIAA J., 3, 1445–1450, 1965. 867

Sakov, P. and Bertino, L.: Relation between two common localisation methods for the EnKF, Computat. Geosci., 10, 225–237, doi:10.1007/s10596-010-9202-6, 2011. 869

Sakov, P., Oliver, D. S., and Bertino, L.: An iterative EnKF for strongly nonlinear systems, Mon. Weather Rev., 140, 1988–2004, doi:10.1175/MWR-D-11-00176.1, 2012. 868, 869

30

Strang, G. and Borre, K.: Linear Algebra, Geodesy, and GPS, Wellesley-Cambridge Press, Wellesley, MA, 490 pp., 1997. 867
Stroud, J. R., Stein, M. L., Lesht, B. M., Schwab, D. J., and Beletsky, D.: An ensemble Kalman filter and smoother for satellite data assimilation, J. Am. Stat. Assoc., 105, 978–990, doi:10.1198/jasa.2010.ap07636, 2010. 867

Trémolet, Y.: Model-error estimation in 4D-Var, Q. J. Roy. Meteor. Soc., 133, 1267–1280, doi:10.1002/qj.94, 2007. 866, 870

Trémolet, Y.: Object-Oriented Prediction System, available at: http://www.data-assimilation.net/Events/Year3/OOPS.pdf (last access: 15 May 2015), 2013. 883

Tshimanga, J., Gratton, S., Weaver, A. T., and Sartenaer, A.: Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation, Q. J. Roy. Meteor. Soc., 134, 751–769, doi:10.1002/qj.228, 2008. 867, 871

Wang, X.: Incorporating ensemble covariance in the gridpoint statistical interpolation variational minimization: a mathematical framework, Mon. Weather Rev., 138, 2990–2995, doi:10.1175/2010MWR3245.1, 2010. 868

Wright, S. J. and Holt, J. N.: An inexact Levenberg–Marquardt method for large sparse nonlinear least squares, J. Aust. Math. Soc. B, 26, 387–403, doi:10.1017/S0334270000004604, 1985. 879

Zhang, F., Zhang, M., and Hansen, J.: Coupling ensemble Kalman filter with four-dimensional variational data assimilation, Adv. Atmos. Sci., 26, 1–8, doi:10.1007/s00376-009-0001-8, 2009. 869

Zupanski, M.: Maximum likelihood ensemble filter: theoretical aspects, Mon. Weather Rev., 133, 1710–1726, doi:10.1175/MWR2946.1, 2005. 868

<figure><img src="image_95.png" /></figure>

<figure><img src="image_96.png" /></figure>

<figure><img src="image_97.png" /></figure>

<figure><img src="image_98.png" /></figure>
<table><caption>**Table 1.** The root mean square error given by Eq. (25) for the first six Gauss–Newton iterations of Algorithm 3, for Lorenz 63 problem. The whole state is observed. Ensemble size is 100. The time window length is 50 time steps. Finite differences parameter is $10^{-3}$.</caption><thead><tr><th>Iteration</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>RMSE</td><td>20.16</td><td>15.37</td><td>3.73</td><td>2.53</td><td>0.09</td><td>0.09</td></tr></tbody></table>

<figure><img src="image_99.png" /></figure>

<figure><img src="image_100.png" /></figure>

<figure><img src="image_101.png" /></figure>

<figure><img src="image_102.png" /></figure>
<table><caption>**Table 2.** Mean of the objective function from 30 runs of the EnKS-4DVAR algorithm (Algorithm 3) for the Lorenz 63 problem and for different values of $\tau$ (finite differences parameter). The whole state is observed. Ensemble size is 50. The time window length is 50 time steps.</caption><thead><tr><th>Iter.</th><th>τ = 1</th><th>τ = 10<sup>-1</sup></th><th>τ = 10<sup>-2</sup></th><th>τ = 10<sup>-3</sup></th><th>τ = 10<sup>-4</sup></th><th>τ = 10<sup>-5</sup></th><th>τ = 10<sup>-6</sup></th></tr></thead><tbody><tr><td>1</td><td>1.02e+06</td><td>1.39e+09</td><td>3.21e+09</td><td>3.54e+09</td><td>3.58e+09</td><td>3.58e+09</td><td>3.58e+09</td></tr><tr><td>2</td><td>1.39e+06</td><td>5.27e+07</td><td>1.70e+08</td><td>1.93e+08</td><td>1.96e+08</td><td>1.96e+08</td><td>1.96e+08</td></tr><tr><td>3</td><td>1.32e+06</td><td>4.14e+06</td><td>2.99e+06</td><td>3.69e+06</td><td>3.76e+06</td><td>3.77e+06</td><td>3.77e+06</td></tr><tr><td>4</td><td>1.38e+06</td><td>5699</td><td>3266</td><td>4431</td><td>4581.31</td><td>4594</td><td>4598</td></tr><tr><td>5</td><td>1.55e+06</td><td>1299</td><td>89.22</td><td>65.69</td><td>65.4442</td><td>65.41</td><td>65.26</td></tr><tr><td>6</td><td>1.34e+06</td><td>830.1</td><td>17.08</td><td>6.933</td><td>6.844</td><td>6.856</td><td>6.923</td></tr><tr><td>7</td><td>2.05e+06</td><td>826.8</td><td>10.75</td><td>1.885</td><td>1.89082</td><td>1.8</td><td>1.721</td></tr><tr><td>8</td><td>1.47e+06</td><td>847.4</td><td>10.82</td><td>1.68</td><td>1.63813</td><td>1.547</td><td>1.641</td></tr></tbody></table>

Title Page

Abstract

Introduction

Conclusions

References

Tables

Figures

<figure><img src="image_103.png" /></figure>

<figure><img src="image_104.png" /></figure>

<figure><img src="image_105.png" /></figure>

<figure><img src="image_106.png" /></figure>

Back

Close

Full Screen / Esc

Printer-friendly Version

Interactive Discussion
<figure><img src="image_107.png" /><figcaption>**Figure 1.** The Lorenz attractor, initial values $x(0) = 1$, $y(0) = 1$, and $z(0) = 1$, discretization time step is $dt = 0.1$ time unit.</figcaption></figure>

<figure><img src="image_108.png" /></figure>

<figure><img src="image_109.png" /></figure>

<figure><img src="image_110.png" /></figure>

<figure><img src="image_111.png" /></figure>

<figure><img src="image_112.png" /></figure>
<figure><img src="image_113.png" /><figcaption>**Figure 2.** The three components *x*, *y*, *z* of the truth and the first five Gauss-Newton iterations from Lorenz 63 problem, for the first 10 time steps. The initial conditions for the truth are *x*(0) = 1, *y*(0) = 1, and *z*(0) = 1. Time step is *dt* = 0.1 time unit. Observations are the full state at each time step. Ensemble size is 100. The time window length is 50 time steps. Finite differences parameter is 10⁻³.</figcaption></figure>

<figure><img src="image_114.png" /></figure>

<figure><img src="image_115.png" /></figure>

<figure><img src="image_116.png" /></figure>

<figure><img src="image_117.png" /></figure>

<figure><img src="image_118.png" /></figure>
<figure><img src="image_119.png" /><figcaption>**Figure 3.** Root square error given by Eq. (24) for the first five Gauss–Newton iterations from Lorenz 63 problem. The problem setting is the same as in Fig. 2.</figcaption></figure>

<figure><img src="image_120.png" /></figure>

<figure><img src="image_121.png" /></figure>

<figure><img src="image_122.png" /></figure>

<figure><img src="image_123.png" /></figure>

<figure><img src="image_124.png" /></figure>
<figure><img src="image_125.png" /></figure>

<figure><img src="image_126.png" /></figure>

<figure><img src="image_127.png" /></figure>

<figure><img src="image_128.png" /></figure>

<figure><img src="image_129.png" /></figure>

<figure><img src="image_130.png" /><figcaption>**Figure 4.** Box plots of objective function values for Lorenz 63 problem. From the left to the right and from the top to the bottom, the figures correspond to the results of the first, the second, the third and the fourth iteration respectively. The whole state is observed. Ensemble size is 50. The time window length is 50 time steps. In each box, the central line presents the median (red line), the edges are the 25th and 75th percentiles (blue line), the whiskers extend to the most extreme data points the plot algorithm considers to be not outliers (black line), and the outliers are plotted individually (red dots).</figcaption></figure>
<figure><img src="image_131.png" /></figure>

<figure><img src="image_132.png" /></figure>

<figure><img src="image_133.png" /></figure>

<figure><img src="image_134.png" /></figure>

<figure><img src="image_135.png" /><figcaption>**Figure 5.** Same as Fig. 4, but for the fifth, the sixth, the seventh and the eighth iteration respectively.</figcaption></figure>
<figure><img src="image_136.png" /><figcaption>**Figure 6.** Objective function values along incremental 4DVAR iterations, for two-level quasi-geostrophic problem from Sect. 6.3.2.</figcaption></figure>

<figure><img src="image_137.png" /></figure>

<figure><img src="image_138.png" /></figure>

<figure><img src="image_139.png" /></figure>

<figure><img src="image_140.png" /></figure>
<figure><img src="image_141.png" /></figure>

<figure><img src="image_142.png" /></figure>

<figure><img src="image_143.png" /></figure>

<figure><img src="image_144.png" /></figure>

<figure><img src="image_145.png" /><figcaption>**Figure 7.** Objective function values along EnKS-4DVAR with regularization iterations for two-level quasi-geostrophic problem (Sect. 6.3.2). From the left to the right and from the top to the bottom: $\gamma = 0, \gamma = 0.001, \gamma = 0.1, \gamma = 1$.</figcaption></figure>
<figure><img src="image_146.png" /></figure>

<figure><img src="image_147.png" /></figure>

<figure><img src="image_148.png" /></figure>

<figure><img src="image_149.png" /></figure>

<figure><img src="image_150.png" /></figure>

<figure><img src="image_151.png" /><figcaption>**Figure 8.** As Fig. 7, but for $\gamma = 10, \gamma = 100, \gamma = 500, \gamma = 1000$, respectively.</figcaption></figure>