A- A+
Alt. Display

# A non-linear least squares enhanced POD-4DVar algorithm for data assimilation

## Abstract

This paper presents a novel non-linear least squares enhanced proper orthogonal decomposition (POD)-based 4DVar algorithm (referred as NLS-4DVar) for the non-linear ensemble-based 4DVar. In the algorithm, the Gauss–Newton iterative method is employed to handle the non-quadratic non-linearity of the 4DVar cost function while the overall structure of the algorithm still resembles the original POD-4DVar algorithm. It is proved that the original POD-4DVar algorithm is a special case of the proposed NLS-4DVar algorithm under the assumption of the linear relationship between the model perturbations (MPs) and the simulated observation perturbations (OPs). Under the assumption it is also shown that the solution of POD-4DVar algorithm coincides with the solution of the proposed NLS-4DVar algorithm. On the contrary, if the linear relationship assumption is dropped, the solution of the POD-4DVar algorithm is only the first iteration of the proposed NLS-4DVar algorithm. As a result, our analysis provides an explanation for the degraded and inaccurate performance of the POD-4DVar algorithm when the underlying forecast model or (and) the observation operator is strongly non-linear. The potential merits and advantages of the proposed NLS-4DVar are demonstrated by a group of Observing System Simulation Experiments with Advanced Research WRF (ARW) using accumulated rainfall-observations.

Keywords:
How to Cite: Tian, X. and Feng, X., 2015. A non-linear least squares enhanced POD-4DVar algorithm for data assimilation. Tellus A: Dynamic Meteorology and Oceanography, 67(1), p.25340. DOI: http://doi.org/10.3402/tellusa.v67.25340
Published on 01 Dec 2015
Accepted on 1 Dec 2014            Submitted on 1 Jul 2014

## 1. Introduction

It is well known that four-dimensional variational data assimilation (4DVar) (e.g. Lewis and Derber, 1985; Le Dimet and Talagrand, 1986; Courtier and Talagrand, 1987; Courtier et al., 1994; Bormann and Thepaut, 2004; Park and Zou, 2004; Caya et al., 2005; Bauer et al., 2006; Rosmond and Xu, 2006; Gauthier et al., 2007) has provided a useful and robust tool for numerical weather prediction (NWP), especially as more and more observation data becomes available via remote sensing techniques. 4DVar has the following two attractive features: (1) the physical model provides a temporal smoothness constraint; (2) it can simultaneously assimilate the observations at multiple time points in an assimilation window. Generally speaking, the successful use of the incremental method (Courtier et al., 1994) and adjoint technique (e.g. Lewis and Derber, 1985; Le Dimet and Talagrand, 1986; Courtier and Talagrand, 1987) has opened a new era in the numerical weather forecast and they are widely used at many operational numerical weather forecast centres (e.g. Bormann and Thepaut, 2004; Park and Zou, 2004; Caya et al., 2005; Bauer et al., 2006; Rosmond and Xu, 2006; Gauthier et al., 2007). Unfortunately, due to its strong dependence on the adjoint model of the forecast model (Bormann and Thepaut, 2004; Park and Zou, 2004; Caya et al., 2005; Bauer et al., 2006; Rosmond and Xu, 2006; Gauthier et al., 2007), computer implementation of 4DVar is often expensive and not easy to program because of the need to maintain and update the adjoint model of the forecast model. The situation is even worse when the forecast model is highly non-linear and/or the model physics involves discontinuous parameters (Xu, 1996). Furthermore, the applications of the static, thus flow-independent, background-error covariance in the traditional 4DVar often causes its poor assimilation performance (Beck and Ehrendorfer, 2005; Zhang et al., 2009; Cheng et al., 2010).

On the contrary, ensemble Kalman filter (EnKF; cf. Evensen, 1994, 2004; Houtekamer and Mitchell, 1998; Keppenne, 2000; Anderson, 2001; Houtekamer and Mitchell, 2001; Whitaker and Hamill, 2002; Ott et al., 2004; Zhang et al., 2004; Hunt et al., 2007) provides another advanced method for data assimilation. EnKF becomes increasingly attractive due to its simple conceptual formulation and relative ease of implementation. Moreover, it has the capability to provide flow-dependent error estimates from the statistics of the ensemble samples. However, unlike 4DVar, EnKF cannot guarantee the dynamic balance in its results since it is essentially a sequential data assimilation method.

The foregoing analysis shows that in order to maximise their full potential 4DVar and EnKF should be coupled to form a more efficient composite method, which seems a natural idea to compensate the shortcomings of both methods by combining their strength (e.g. Lorenc, 2003; Harlim and Hunt, 2007; Tian et al., 2008, 2011; Wang et al., 2010; Zhang et al., 2009; Cheng et al., 2010; Zhang and Zhang, 2012; Wang et al., 2013). According to the classification suggested by Andrew C. Lorenc (http://www.wcrp-climate.org/WGNE/BlueBook/2013/individual-articles/01_Lorenc_Andrew_EnVar_nomenclature.pdf), the existing composite methods in the literature can be roughly divided into three groups, namely, the ‘hybrid’ (e.g. Wang et al., 2008a, 2008b; Clayton et al., 2013), the ‘En4DVar’ (e.g. Zhang et al., 2009; Zhang and Zhang, 2012) and the ‘4DEnVar’ assimilation methods (Qiu et al., 2007; Liu et al., 2008; Tian et al., 2008; Wang et al., 2010; Tian et al., 2011; Tian and Xie, 2012). The ‘hybrid 4DVar’ denotes 4DVar methods using a combination of climatological and ensemble covariances. In the ‘En4DVar’, the EnKF and 4DVar runs are implemented in parallel and the information is exchanged complementarily (Zhang and Zhang, 2012): The resulted data assimilation system is based on an ensemble forecasting framework in which EnKF is used to estimate the ensemble background covariance and to update the analysis perturbations, while the 4DVar analysis is performed afterward and is taken as the final analysis for the combined 4DVar and EnKF method. Here 4DVar is implemented in its traditional way except the ensemble (flow-dependent) background covariances being used in place of its original static ones as inputs to the 4DVar minimisation. It should be pointed out that the adjoint model is still indispensable in the ‘hybrid’ and ‘En4DVar’ assimilation methods because it is used in the traditional 4DVar implementation. Similar to ‘En4DVar’, the ‘4DEnVar’ assimilation methods are also based on the ensemble forecasting framework. They start from the 4DVar cost function and are derived under the assumption that the 4DVar analysis increments belong to the linear space expanded by the ensemble-based simulated samples. Therefore, the analysis increments could be expressed as linear combinations of the model perturbations (MPs). Substituting such a linear combination into the 4DVar cost function, the control (state) variables, which are the coefficients in the linear combination, can be computed explicitly under the assumption of linear dependence between the MPs and observation perturbations (Ops). As a result, the adjoint model is no longer needed and the assimilation process is significantly simplified. However the linear dependence assumption between the MPs and OPs is likely to be challenged by the forecast model and/or by the observation operator which are highly non-linear. To a large extent, the 4DEnVar assimilation methods are very similar to the Ensemble Kalman Smoother (Evensen and van Leeuwen, 2000; Harlim and Hunt, 2007; Milewski and Bourqui, 2013). As the theoretical foundation for many 4DEnVar assimilation methods, the assumption of the linear dependence between the MPs and OPs has often been questioned. It is of great theoretical and practical importance to understand the impact of the linear dependence assumption to the assimilation performance of 4DEnVar, especially, when the forecast model and/or the observation operator are highly non-linear. For example, the proper orthogonal decomposition (POD)-based ensemble 4/3DVar (referred as POD-4/3DVar, Tian et al., 2011) has been successfully applied to real carbon cycle data assimilation (Tian et al., 2014) and radar assimilation (Pan et al., 2012), which shows a robust performance in the atmospheric transport data assimilation. Furthermore, we have also built a POD-4DVar based land data assimilation system on the community land model platform (Tian et al., 2009, 2010). Very recently results from the real data assimilation experiments with the Weather Research and Forecasting (WRF) model show that significant improvements in predicting the convective system and thus precipitation are achieved due to improved initial conditions for the storm's dynamics and microphysics through POD-4DVar-based assimilation of the radial velocity and reflectivity data (Zhang et al., 2014). Nevertheless, those experiments also demonstrate an appropriate assimilation window has to be chosen for POD-4DVar to guarantee its linear assumption between the OPs and MPs and to incorporate the observations as much as possible to provide more observation information when applying it to conduct real radar data assimilation for a heavy convective-rainfall case study (Zhang et al., 2014). One of main goals of this paper is to provide an explanation to such a fundamental issue.

To address the above issue, we first convert the ensemble-based 4DVar problem into an equivalent non-linear least squares problem and then propose a Gauss–Newton type iterative algorithm to approximate the solution of the non-linear 4DVar problem, which has to be terminated after finite number of iterations in practice according to some pre-determined computationally acceptable stopping criterion. Remarkably, we show by a simple analysis that the original POD-based 4DVar approach (Tian et al., 2008, 2011) is a special case of the proposed non-linear least squares approach which is amount to taking the first iteration of our Gauss–Newton iterative algorithm as the solution of the 4DEnVar. Therefore, it is expected that such a crude approach in general gives a fairly rough approximation to the NLS-4DVar optimal solution when the forecast model and/or the observation operator are highly non-linear, although it does provide an accurate solution under the assumption of the linear dependence between the MPs and OPs.

The rest of the paper is organised as follows. In Section 2, we first derive an equivalent non-linear least squares formulation for the non-linear ensemble-based 4DVar and then propose a Gauss–Newton type iterative algorithm for numerically solving the resulted non-linear least squares problem. In Section 3, we present a group of observing system simulation experiments (referred to as OSSEs) with Advanced Research WRF (ARW) using accumulated rainfall-observations to gauge the performance of our NLS-4DVar approach. Finally, the paper is concluded with a summary and a few concluding remarks given in Section 4.

## 2. Methodology

By minimising the following incremental form of the 4DVar cost function, one can obtain an optimal increment of the initial condition (IC), ${\mathbf{x}}_{a}^{\text{'}}$, at the initial time t0

(1a )
$\text{J}\left({\mathbf{x}}^{\text{'}}\right)=\frac{1}{2}{\left({\mathbf{x}}^{\text{'}}\right)}^{T}{\mathbf{B}}^{-1}\left({\mathbf{x}}^{\text{'}}\right)+\frac{1}{2}\sum \text{k}=1\text{S}{\left[{\text{H}}_{\text{k}}{\text{M}}_{t\to t\text{k}}\left({\mathbf{x}}_{b}+{\mathbf{x}}^{\text{'}}\right)-{\mathbf{y}}_{obs}\right]}^{\mathbf{T}}{\mathbf{R}}_{\text{k}}^{-1}\left[{\text{H}}_{\text{k}}{\text{M}}_{t\to t\text{k}}\left({\mathbf{x}}_{b}+{\mathbf{x}}^{\text{'}}\right)-{\mathbf{y}}_{obs}\right].$

Equation (1a) can be rewritten as follows

(1b )
$\text{J}\left({\mathbf{x}}^{\text{'}}\right)=\frac{1}{2}{\left({\mathbf{x}}^{\text{'}}\right)}^{T}{\mathbf{B}}^{-1}\left({\mathbf{x}}^{\text{'}}\right)+\frac{1}{2}{\left[{\text{L}}^{\text{'}}\left({\mathbf{x}}^{\text{'}}\right)-{\mathbf{y}}_{obs}^{\text{'}}\right]}^{\mathbf{T}}{\mathbf{R}}^{-1}\left[{\text{L}}^{\text{'}}\left({\mathbf{x}}^{\text{'}}\right)-{\mathbf{y}}_{obs}^{\text{'}}\right],$

where x=x−xb is the perturbation of the background field xb at t0,

(2 )
${\mathbf{y}}_{\text{obs}}^{\text{'}}=\left[\begin{array}{c}{\mathbf{y}}_{\text{obs},1}^{\text{'}}\\ {\mathbf{y}}_{\text{obs},2}^{\text{'}}\\ ⋮\\ {\mathbf{y}}_{\text{obs},\text{S}}^{\text{'}}\end{array}\right],$
(3 )
${L}^{\text{'}}=\left(\begin{array}{c}{L}_{1}^{\text{'}}\\ {L}_{2}^{\text{'}}\\ ⋮\\ {L}_{\text{S}}^{\text{'}}\end{array}\right),$
(4 )
${L}_{k}^{\text{'}}\left({\mathbf{x}}^{\text{'}}\right)={L}_{k}\left({\mathbf{x}}_{b}+{\mathbf{x}}^{\text{'}}\right)-{L}_{k}\left({\mathbf{x}}_{b}\right),$
(5 )
${\mathbf{y}}_{\text{obs},k}^{\text{'}}={\mathbf{y}}_{\text{obs},k}-{L}_{k}\left({\mathbf{x}}_{b}\right),$
(6 )
${L}_{k}={H}_{k}{M}_{t0\to {t}_{k}},$

and

(7 )
$\mathbf{R}=\left[\begin{array}{cccc}{\mathbf{R}}_{1}& 0& \cdots & 0\\ 0& {\mathbf{R}}_{2}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\mathbf{R}}_{\text{S}}\end{array}\right].$

Here the superscript T stands for matrix transpose, b is the background value, tk denotes the kth observation time point, S is the total number of observation time points in the assimilation window, Hk is the observation operator, and matrices B and Rk are the background and observation error covariances, respectively.

Similar to the original POD-4DVar approach (Tian et al., 2008, 2011), the NLS-4DVar approach also starts from an ensemble of N initial fields xj(j=1,…,N) and further assumes that their analysis solution ${\mathbf{x}}_{a}^{\text{'}}$ is expressed by the linear combinations of the MPs (${x}_{\text{j}}^{\text{'}}={x}_{\text{j}}-{x}_{b}$, j=1,…, N) as follows

(8 )
${\mathbf{x}}_{\text{a}}^{\text{'}}={\mathbf{P}}_{\text{x}}\beta ,$

where ${\mathbf{P}}_{\text{x}}=\left({\mathbf{x}}_{1}^{\text{'}},{\mathbf{x}}_{2}^{\text{'}},\cdots ,{\mathbf{x}}_{\text{N}}^{\text{'}}\right)$ and $\beta ={\left({\beta }_{1},{\beta }_{2},\cdots ,{\beta }_{\text{N}}\right)}^{\mathbf{T}}$.

Likewise, we set

(9 )
${\mathbf{P}}_{\text{y}}=\left({\mathbf{y}}_{1}^{\text{'}},{\mathbf{y}}_{2}^{\text{'}},\cdots ,{\mathbf{y}}_{\text{N}}^{\text{'}}\right),$

where ${\mathbf{y}}_{\text{j}}^{\text{'}}={\text{L}}^{\text{'}}\left({\mathbf{x}}_{\text{j}}^{\text{'}}\right)$. As in Tian et al. (2010), the POD (Ly and Tran, 2001, 2002) transformations to the OP and MP matrixes (Py and Px) are also conducted in our NLS-4DVar formulations to guarantee the orthogonality of $\left({\mathbf{y}}_{1}^{\text{'}},{\mathbf{y}}_{2}^{\text{'}},\cdots ,{\mathbf{y}}_{\text{N}}^{\text{'}}\right)$. For notation brevity, we still use Px and Py to represent the POD-transformed OP and MP matrixes.

Substituting eq. (8) into eq. (1) and expressing the cost function in terms of a new control variable β yield

(10 )
$J\left(\beta \right)=\frac{1}{2}{\beta }^{\mathbf{T}}{\left({\mathbf{P}}_{\text{x}}\right)}^{\mathbf{T}}{\mathbf{B}}^{-1}\left({\mathbf{P}}_{\text{x}}\right)\beta +\frac{1}{2}{\left[{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)-{\mathbf{y}}_{\text{obs}}^{\text{'}}\right]}^{\mathbf{T}}{\mathbf{R}}_{}^{-1}\left[{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)-{\mathbf{y}}_{\text{obs}}^{\text{'}}\right],$

In addition, substituting the ensemble background covariance $\mathbf{B}=\frac{\left({\mathbf{P}}_{\text{x}}\right){\left({\mathbf{P}}_{\text{x}}\right)}^{\mathbf{T}}}{N-1}$ (Evensen, 2004; Tian et al., 2011), eq. (10) is further reduced into the following form:

(11 )
$J\left(\beta \right)=\frac{1}{2}\left(N-1\right)\cdot {\beta }^{\mathbf{T}}{\left({\mathbf{P}}_{\text{x}}\right)}^{\mathbf{T}}{\left[{\left({\mathbf{P}}_{\text{x}}\right)}^{\mathbf{T}}\right]}^{-1}{\left({\mathbf{P}}_{\text{x}}\right)}^{-1}\left({\mathbf{P}}_{\text{x}}\right)\beta +\frac{1}{2}{\left[{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)-{\mathbf{y}}_{\text{obs}}^{\text{'}}\right]}^{\mathbf{T}}{\mathbf{R}}_{}^{-1}\left[{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)-{\mathbf{y}}_{\text{obs}}^{\text{'}}\right]=\frac{1}{2}\left(N-1\right)\cdot {\beta }^{\mathbf{T}}\beta +\frac{1}{2}{\left[{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)-{\mathbf{y}}_{\text{obs}}^{\text{'}}\right]}^{\mathbf{T}}{\mathbf{R}}_{}^{-1}\left[{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)-{\mathbf{y}}_{\text{obs}}^{\text{'}}\right]$

where the inverses of (Px)T and (Px) should be understood as the Moore–Penrose or pseudo-inverse (Evensen, 2004).

Furthermore, since R is symmetric, it has the Cholesky factorisation

(12 )
$\mathbf{R}={\mathbf{R}}_{+}^{1/2}{\left({\mathbf{R}}_{+}^{1/2}\right)}^{\mathbf{T}},$

hence,

(13 )
$J\left(\beta \right)=\frac{1}{2}\left(N-1\right){\beta }^{\mathbf{T}}\beta +\frac{1}{2}{\left[{\mathbf{R}}_{+}^{-1/2}{L}^{\text{'}}\left({\mathbf{P}}_{x}\beta \right)-{\mathbf{R}}_{+}^{-1/2}{\mathbf{y}}_{\text{obs}}^{\text{'}}\right]}^{\mathbf{T}}\left[{\mathbf{R}}_{+}^{-1/2}{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)-{\mathbf{R}}_{+}^{-1/2}{\mathbf{y}}_{\text{obs}}^{\text{'}}\right].$

Introducing

(14 )
$Q\left(\beta \right)=\left(\begin{array}{c}{\mathbf{R}}_{+}^{-1/2}\left[{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)-{\mathbf{y}}_{\text{obs}}^{\text{'}}\right]\\ \sqrt{N-1}\beta \end{array}\right),$

then eq. (11) can be rewritten as follows

(15 )
$J\left(\beta \right)=\frac{1}{2}Q{\left(\beta \right)}^{\mathbf{T}}Q\left(\beta \right).$

Apparently, the original 4DVar problem (1–7) has been transformed into the above non-linear least squares formulation (Dennis and Schnabel, 1996) through the steps given in eqs. (8–14).

Both approximations

(16 )
${L}^{\text{'}}\left({\mathbf{x}}_{\text{j}}^{\text{'}}\right)={\mathbf{y}}_{\text{j}}^{\text{'}}\approx {\mathbf{H}}_{}^{\text{'}}{\mathbf{M}}_{}^{\text{'}}\left({\mathbf{x}}_{\text{j}}^{\text{'}}\right),$

and

(17 )
${L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)\approx {\mathbf{H}}^{\text{'}}{\mathbf{M}}^{\text{'}}\left({\mathbf{P}}_{\text{x}}\beta \right)\approx \sum _{\text{j}=1}^{N}{\mathbf{y}}_{\text{j}}^{\text{'}}{\beta }_{\text{j}}={\mathbf{P}}_{\text{y}}\beta .$

are widely used in the literature. Here the tangent linear operator ${\mathbf{H}}_{}^{\text{'}}{\mathbf{M}}_{}^{\text{'}}$ is given by

(18 )
${\mathbf{H}}_{}^{\text{'}}{\mathbf{M}}_{}^{\text{'}}=\left[\begin{array}{cccc}{L}_{1}^{\text{'}}& 0& \cdots & 0\\ 0& {L}_{2}^{\text{'}}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {L}_{\text{S}}^{\text{'}}\end{array}\right],$

and

(19 )
${L}_{k}^{\text{'}}\left({\mathbf{x}}^{\text{'}}\right)={H}_{k}{M}_{t0\to {t}_{k}}\left({\mathbf{x}}_{b}+{\mathbf{x}}^{\text{'}}\right)-{H}_{k}{M}_{{t}_{0}\to {t}_{k}}\left({\mathbf{x}}_{b}\right),\text{}k=1,\cdots ,S,$

Consequently, the first-derivate matrix (or Jacobian matrix) JacQ(β) of Q(β) can be computed approximately as follows

(20 )
${J}_{\text{ac}}Q\left(\beta \right)=\frac{\partial Q\left(\beta \right)}{\partial \beta }\approx \left(\begin{array}{c}{\mathbf{R}}_{+}^{-1/2}{\mathbf{P}}_{\text{y}}^{}\\ \sqrt{N-1}\mathbf{I}\end{array}\right),$

where I denotes the N×N identity matrix. Therefore, the Gauss–Newton iteration for the non-linear least squares problem (15) is defined by (Dennis and Schnabel, 1996)

(21 )
${\beta }^{i+1}={\beta }^{i}-{\left[{\left({J}_{\text{ac}}Q\left({\beta }^{i}\right)\right)}^{\mathbf{T}}\left({J}_{\text{ac}}Q\left({\beta }^{i}\right)\right)\right]}^{-1}{\left({J}_{\text{ac}}Q\left({\beta }^{i}\right)\right)}^{\mathbf{T}}Q\left({\beta }^{i}\right)$

for i=0,1,2,…, Imax (Imax is the maximum iteration number).

Substituting (14) and (20) into (21) and rearranging terms we get

(22 )
${\beta }^{i+1}={\beta }^{i}-{\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}_{}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}\left(N-1\right){\beta }^{i}+{\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left[{\mathbf{y}}_{\text{obs}}^{\text{'}}-{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}{\beta }^{i}\right)\right].$

Substituting (8) into (22) yields

(23 )
${\beta }^{i+1}={\beta }^{i}-{\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}_{}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}\left(N-1\right){\beta }^{i}+{\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left[{\mathbf{y}}_{\text{obs}}^{\text{'}}-{L}^{\text{'}}\left({\mathbf{x}}_{\text{a}}^{\text{'},i}\right)\right].$

Using eq. (17) as well as the orthogonality of $\left({\mathbf{y}}_{1}^{\text{'}},{\mathbf{y}}_{2}^{\text{'}},\cdots ,{\mathbf{y}}_{\text{N}}^{\text{'}}\right)$, βi can be expressed as

(24 )
${\beta }^{i}={\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{P}}_{\text{y}}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{L}^{\text{'}}\left({\mathbf{P}}_{\text{x}}{\beta }^{i}\right)={\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{P}}_{\text{y}}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{L}^{\text{'}}\left({\mathbf{x}}_{\text{a}}^{\text{'}i}\right).$

Equation (23) can be further rewritten as

(25 )
${\beta }^{i+1}={\beta }^{i}-{\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}_{}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}\left(N-1\right){\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{P}}_{\text{y}}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{L}^{\text{'}}\left({\mathbf{x}}_{\text{a}}^{\text{'},i}\right)+{\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left[{\mathbf{y}}_{\text{obs}}^{\text{'}}-{L}^{\text{'}}\left({\mathbf{x}}_{\text{a}}^{\text{'},i}\right)\right],$
(26 )
${\mathbf{P}}_{\text{x}}{\beta }^{i+1}={\mathbf{P}}_{\text{x}}{\beta }^{i}-{\mathbf{P}}_{\text{x}}{\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}_{}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}\left(N-1\right){\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{P}}_{\text{y}}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{L}^{\text{'}}\left({\mathbf{x}}_{\text{a}}^{\text{'},i}\right)+{\mathbf{P}}_{\text{x}}{\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left[{\mathbf{y}}_{\text{obs}}^{\text{'}}-{L}^{\text{'}}\left({\mathbf{x}}_{\text{a}}^{\text{'}i}\right)\right].$

Hence

(27 )
${\mathbf{x}}_{\text{a}}^{\text{'},i+1}={\mathbf{x}}_{\text{a}}^{\text{'},i}+{\mathbf{P}}_{\text{x}}\stackrel{}{{\mathbf{P}}_{\text{y}}}{L}^{\text{'}}\left({\mathbf{x}}_{\text{a}}^{\text{'},i}\right)+{\mathbf{P}}_{\text{x}}\stackrel{}{{\mathbf{P}}_{\text{y}}}\left[{\mathbf{y}}_{\text{obs}}^{\text{'}}-{L}^{\text{'}}\left({\mathbf{x}}_{\text{a}}^{\text{'},i}\right)\right],$

where

(28 )
(29 )
$\stackrel{}{{\mathbf{P}}_{\text{y}}}={\left[{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}\left({\mathbf{P}}_{\text{y}}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}{\left({\mathbf{P}}_{\text{y}}\right)}^{\mathbf{T}}{\mathbf{R}}^{-1}.$

Introducing the time index k(k=1,…,S), eqs. (27–29) can be rewritten as

(27b )
(28b )

and

(29b )
$\stackrel{}{{\mathbf{P}}_{\text{y},k}}={\left[\sum _{k=1}^{S}{\left({\mathbf{P}}_{\text{y},k}\right)}^{\mathbf{T}}\left({\mathbf{P}}_{\text{y},k}\right)+\left(N-1\right)\mathbf{I}\right]}^{-1}{\left({\mathbf{P}}_{\text{y},k}\right)}^{\mathbf{T}}{\mathbf{R}}_{k}^{-1}.$

In order to filter out the spurious long-range correlations (Houtekamer and Mitchell, 2001) between the initial perturbations and the corresponding prediction increments on model grids, we replace the matrixes ${\mathbf{P}}_{x}\stackrel{}{{\mathbf{P}}_{y}}$ and y by $\rho \circ {\mathbf{P}}_{x}\stackrel{}{{\mathbf{P}}_{y}}$ and , respectively, where $\mathbf{B}\circ \mathbf{C}$ stands for the Schur product of matrices B and C of the same dimension which is a matrix whose (i,j) entry is given by bi,j·ci,j. Thus we obtain

(30 )

which can be rewritten explicitly as

(30b )

The elements of the matrix ρ in eq. (30) can be calculated according to the formula

(31 )
${\rho }_{\text{i},\text{j}}={C}_{0}\left({d}_{\text{h},\text{i},\text{j}}/{d}_{\text{h},0}\right)\cdot C\left({d}_{\text{v},\text{i},\text{j}}/{d}_{\text{v},0}\right)0,$

where the filtering function C0 is defined by (Gaspari and Cohn, 1999)

(32 )
${C}_{0}\left(r\right)=\left\{\begin{array}{ll}-\frac{1}{4}{r}^{5}+\frac{1}{2}{r}^{4}+\frac{5}{8}{r}^{3}-\frac{5}{3}{r}^{2}+1,\hfill & 0\le r\le 1,\hfill \\ \frac{1}{12}{r}^{5}-\frac{1}{2}{r}^{4}+\frac{5}{8}{r}^{3}+\frac{5}{3}{r}^{2}-5r+4-\frac{2}{3}{r}^{-1},\hfill & 1

and dh,0 and dv,0 are the horizontal and vertical covariance localisation Schur radii, respectively.

It is easy to see that the background field Xb corresponds to the first Gauss–Newton iterate using the trivial starting field ${\mathbf{x}}_{\text{a}}^{\text{'},0}=0$, which results in

(33 )
${\mathbf{x}}_{\text{a}}^{\text{'},1}=\left(\rho \circ {\mathbf{P}}_{\text{x}}\stackrel{}{{\mathbf{P}}_{\text{y}}}\right){\mathbf{y}}_{\text{obs}}^{\text{'}}.$

In this case (${\mathbf{x}}_{a}^{\text{'},0}=0$), the 4DVar analysis of Xa is computed by

(34 )
${\mathbf{x}}_{\text{a}}^{}={\mathbf{x}}_{\text{b}}^{}+\left(\rho \circ {\mathbf{P}}_{\text{x}}\stackrel{}{{\mathbf{P}}_{\text{y}}}\right){\mathbf{y}}_{\text{obs}}^{\text{'}}.$

Equation (34) is exactly the analysis solution [eq. (38) in Tian et al., 2011] in the POD-based ensemble 4DVar (Tian et al., 2011; Tian and Xie, 2012) approach. Consequently, the POD-4DVar analysis solution is just the first iteration of the proposed NLS-4DVar sequence generated by eq. (27). Moreover, most of the existing ensemble-based 4DVar methods also fall into such a special case of one iteration of the NLS-4DVar approach applied to the ensemble 4DVar optimisation problem (10).

Furthermore, if both Hk and ${M}_{{t}_{0}\to {t}_{k}}$are linear, it is easy to check that ${L}^{\text{'}}={\mathbf{H}}^{\text{'}}{\mathbf{M}}^{\text{'}}={H}_{k}{M}_{{t}_{0}\to {t}_{k}}$ and eq. (27) reduces to

(35 )
${\mathbf{x}}_{\text{a}}^{\text{'}}=\left(\rho \circ {\mathbf{P}}_{\text{x}}\stackrel{}{{\mathbf{P}}_{\text{y}}}\right){\mathbf{y}}_{\text{obs}}^{\text{'}}.$

Equation (35) shows that the NLS-4DVar sequence generated by eq. (27) becomes a constant sequence which coincides with the original POD-4DVar solution under the assumption that both Hk and ${M}_{{t}_{0}\to {t}_{k}}$ are linear.

From its definition, it is not hard to see that the implementation of the NLS-4DVar requires the following steps:

Prepare xb,xj, ${\mathbf{R}}_{k}^{-1}$ and ;

Compute Px,${\mathbf{P}}_{y}^{k}$ and ${\mathbf{y}}_{\text{obs},k}^{\text{'}}$;

Compute $\left(\rho \circ {\mathbf{P}}_{\text{x}}\stackrel{}{{\mathbf{P}}_{\text{y}}}\right)$ and ;

${\mathbf{x}}_{\text{a}}^{\text{'},1}=\left(\rho \circ {\mathbf{P}}_{\text{x}}\stackrel{}{{\mathbf{P}}_{\text{y}}}\right){\mathbf{y}}_{\text{obs}}^{\text{'}}$

do i=1,…,Imax (where Imax is the maximum iteration number)

Compute $L\left({\mathbf{x}}_{a}^{\text{'},i}\right)$ using the non-linear models H and M with the analysis value ${\mathbf{x}}_{a}^{\text{'},i}$ at the ith step;

Compute .

It should be noted that in the above algorithm the value of ${L}^{\text{'}}\left({\mathbf{x}}_{a}^{\text{'},i}\right)$ is calculated using the non-linear forecast model M and the non-linear observation operator H directly, rather than using their corresponding tangent linear operators. Therefore, the NLS-4DVar approach is designed to handle non-linear (and linear) data assimilation.

## 3. Evaluations through OSSEs

An OSSE is considered as one of the best benchmark tests to evaluate a data assimilation methodology since it can provide both the ‘true’ state and the corresponding ‘observation’. On the contrary, real data assimilation is often not an effective mean to assess a new methodology. It is well known that the observation errors in space and time are inevitable in any available observations, therefore, the exact underlying true state is hardly known. Once the observations are assimilated into the IC, there is no basis for evaluating the assimilation results.

In this section, a group of OSSEs are carefully designed to examine the performance of NLS-4DVar in comparison to POD-4DVar. As stated in Wang et al. (2010), rainfall observations through 4DVar may significantly improve rainfall-forecasting skill (Zou and Kuo, 1996) because the precipitation process depends upon many model variables (e.g. wind, temperature, and water vapour mixing ratio) through model dynamics and thermodynamics. Obviously, the observation operator that links the model states and the rainfall observations is highly non-linear and often cannot be explicitly expressed. Naturally, it is an appropriate choice to evaluate the proposed NLS-4DVar method using the rainfall observations. In our experiments, we choose a 6-hour assimilation window with accumulated rainfall observation data available at 3 and 6 hours after the analysis time.

To establish a more realistic framework, we use the WRF model version 3 (Skamarock et al., 2008; http://www.wrf-model.org/index.php, the Advanced Research WRF version or ARW is used in our experiments) as the forecast model for the OSSEs. The main physical components of the WRF model used in our experiments include the Rapid Radiative Transfer Model (RRTM) based on Mlawer et al. (1997) for long-wave radiation, the Dudhia (1989) shortwave radiation scheme, the Yonsei University (YSU) PBL scheme (Hong et al., 2006), the Purde Lin explicit cloud microphysics parameterisation (Lin et al., 1983; Rutledge and Hobbs, 1984; Chen and Sun, 2002), and the Noah LSM land scheme (Chen and Dudhia, 2001). The following experiments (including the NLS-4DVar and POD-4DVar data assimilation and control simulation runs) are done over a mesh of 120×100 (longitude×latitude) grid points with grid spacing of 30 km. The domain covers the whole of China around a region from 18°N~42°N to 95°E~125°E. In the vertical direction, 27 layers from σ=0 to σ=1 are used.

The assumed ‘true’ atmospheric states over 30 hours are produced as follows. A 12-hour spinning-up run is conducted from 1200 UTC 07 June to 0000 UTC 08 June, 2010 with 1°×1° National Centers for Environmental Prediction (NCEP) Final (FNL) Global Tropospheric Analyses, to obtain the ‘true’ initial fields at the analysis time (i.e. 0000 UTC 08 June, 2010), which is used to initiate the 30-hour ‘true’ simulations. The observations are thus generated at 3 and 6 hours after the analysis time by adding Gauss white random noises (with a given constant error standard deviation ~0.2 mm) in space and time to the interpolated ‘true’ accumulated rainfall fields at 565 stations (Fig. 1) in the observation region.

Fig. 1

The 565 observational stations used in this study.

The background fields used in our experiments are prepared very similarly with the ‘true’ case but they are produced from a 24-hour forecast initialised by the NCEP/FNL data at 24 hours prior to the analysis time. A sizeable difference is thus found between the ‘true’ state and the background fields. From the background fields at the analysis time, a 30-hour integration regarded as the control run (referred to as ‘CTL’) is obtained through the assimilation window to produce the basic states for assimilation and evaluate the results from both NLS-4DVar and POD-4DVar assimilation experiments.

A 4D moving sampling strategy (Wang et al., 2010; Tian et al., 2014) is adopted to produce the MPs (Px) and their corresponding simulated OPs (Py) from 78-hour forecasts with one output per hour and the background simulations. The POD-transformed OP and MP matrixes also can easily be obtained. The two 78-hour forecasts are initialised with 1°×1° NCEP/FNL global analysis data, respectively, at 24 and 48 hours prior to the analysis time. Each 78-hour forecast includes 73 six-hour moving windows and the ensemble size is 146 (=73×2) as a result.

In the following experiments, the performance of the NLS-4DVar method is examined in comparison to the POD-4DVar method. To have a fair comparison, all the assimilation settings (including the assimilated variables, the analysis time, the assimilation window, the observations, the sampling strategy, the ensemble size and the localization Schur radii) are the same. It should be noted that, like most 4DEnVar methods, we have to determine the localisation Schur radii for NLS-4DVar either by experience or through sensitivity experiments since only a static correlation matrix (31) is used in the current algorithm, which will definitely bring some uncertainty. This phenomenon is indeed observed in our real data assimilation experiments (e.g. Zhang et al., 2014). To assess the overall performance of the NLS-4DVar method (marked as ‘NLS’), especially in contrast to the POD-4DVar method (marked as ‘POD’), we first compare the root mean square errors (RMSEs) of 30-hour forecast of accumulated rainfall (from 1~30 hours after the analysis time) with the ICs from the ‘NLS’, ‘POD’ and ‘CTL’, respectively. Figure 2 shows the RMSEs of 30-hour forecast of accumulated rainfall with the ICs from the assimilations, both ‘NLS’ and ‘POD’ produce smaller values than those from ‘CTL’, which indicates both ‘NLS’ and ‘POD’ can steadily improve the rainfall forecast over 30 hours. Meanwhile, Fig. 2 also indicates that NLS-4DVar outperforms POD-4DVar moderately, due to the iterative strategy adopted in NLS-4DVar. Additionally, sensitivity experiments show that the NLS-4DVar solution is already convergent only after two iterations (i.e. Imax=2).

Fig. 2

RMSEs of 30-hour forecast of accumulated rainfall (from 1~30 hours) with the ICs from ‘CTL’ (solid line), ‘POD’ (long-dashed line) and ‘NLS’ (short-dashed line).

To further investigate the potential of NLS-4DVar, we also compare the difference of accumulated rainfall at 18 hours after the analysis time with the ICs from ‘NLS’, ‘POD’, ‘CTL’ and ‘truth’ in Fig. 3, respectively. Generally, ‘CTL’ and ‘POD’ share a stronger positive bias region around (26°~30°N, 122°~126°E). Correspondingly, NLS-4DVar only has a little weaker positive bias in the same region, which demonstrates the NLS-4DVar method can not only reduce the overall forecast errors but also ameliorate the detailed patterns of forecasted rainfall substantially. The results also suggest that NLS-4DVar works fairly effectively, as it does help to reduce the errors in both the immediate vicinity and outside of the observations even though the observational error is not negligible and the distribution of observations is sporadic and only at the land surface. Certainly, the superior performance of NLS-4DVar is attributed to its ability to provide a better description of the atmosphere and is closer to the ‘true’ state at the analysis by incorporating the observation information into the background. Figure 4 shows that the difference of temperature at the σ=0.512 level at the analysis time between ‘NLS’ and ‘truth’ is much smaller than those between the ‘CTL’, ‘POD’ and ‘truth’ as a whole. A significant decrease of the difference in the region (24°~41°N, 117°~128°E) appears in ‘NLS’ compared with ‘POD’ and ‘CTL’ (see Fig. 4a, b and c). Meanwhile, Fig. 4b and c also show that POD-4DVar can also improve the accuracy of the IC significantly, which leads to a general decrease (but less than NLS-4DVar) in the difference of temperature at the selected layer. Furthermore, to investigate the effectiveness of NLS-4DVar, we also compare the vertical profiles of RMSEs of some basic model variables from ‘NLS’ and ‘POD’ at the beginning and the end of the assimilation window, calculated at all horizontal model grid points of each layer. The same conclusion can be drawn that NLS-4DVar produces more accurate atmospheric states than POD-4DVar does. As seen in Figs. 5 and 6, NLS-4DVar yields a smaller error than POD-4DVar on most model layers. This implies that the iterative strategy adopted in NLS-4DVar really contributes to its superior performance.

Fig. 3

Accumulated rainfall difference (mm) between (a) ‘NLS’, (b) ‘POD’, (c) ‘CTL’ and the ‘truth’ at 18 hours after the analysis time.

Fig. 4

Temperature difference (K) on the σ=0.512 level between (a) ‘NLS’, (b) ‘POD’, and (c) ‘CTL’ and the ‘true’ state at the analysis time.

Fig. 5

Vertical profiles of RMSEs of (a) zonal wind (m s−1), (b) meridional wind (m s−1), (c) temperature (°C), and (d) water vapour mixing ratio (g kg−1) of ‘NLS’ (solid curve) and ‘POD’ (dashed line) at the start of the assimilation window.

Fig. 6

Same as Fig. 5, but for at the end of the assimilation window.

## 4. Summary and concluding remarks

The existing 4DEnVar methods (e.g. Tian et al., 2008, 2011) are commonly based on the linear assumption between the MPs and OPs as their theoretical foundation. When the forecast model or observation operator is highly non-linear (frequently occurring in the actual assimilation implementations), the above linear assumption becomes problematic. As a result, its theoretical foundation collapses. To address this issue, a non-linear least squares iterative algorithm is proposed to iteratively approximate the optimal non-linear 4DVar solution in a computationally acceptable way, which defines the NLS-4DVar approach. Particularly, the iterative algorithm can produce satisfactory results with a cost of a very few iterations (only three iterations are used in our OSSEs), consequently, it can be implemented more efficiently than the standard adjoint-based 4DVar (Wang et al., 2010).

The effectiveness, robustness and potential merits of the proposed NLS-4DVar method are demonstrated by several sets of OSSEs based on the platform of ARW using accumulated rainfall-observations. For comparison, the original POD-4DVar method is simultaneously used to carry out the same OSSEs under the same conditions. Our main conclusions are summarised as follows:

• NLS-4DVar can be easily implemented as the original POD-4DVar method by utilizing an explicit iterative equation rather than resorting to any existing optimization algorithm;
• The solution of the original POD-4DVar method is only the first iteration of the NLS-4DVar iterative method, which largely explains inferior performance of POD-4DVar.
• A group of OSSEs in a nearly realistic framework demonstrates that the iterative strategy employed in NLS-4DVar is able to alleviate the non-linear problems arising from data assimilation, the resulting NLS-4DVar method not only can reduce the overall forecast errors but also can substantially ameliorate the detailed patterns of forecasted rainfall.
• Applications of NLS-4DVar in the meso- and micro-scale systems are expected to be favourable because of their high non-linearities.
• As an ensemble-based method, NLS-4DVar also needs a moderation function to filter out the spurious long range correlations but only adopts a non-adaptive correlation matrix (which is usually determined either by experience or through sensitivity experiments) in current formulation. Extra efforts need be made to explore a flow-adaptive moderation of spurious ensemble correlations and its use in NLS-4DVar.

## 5. Acknowledgements

This work was partially supported by the National High Technology Research and Development Program of China (Grant No. 2013AA122002), the Knowledge Innovation Program of the Chinese Academy of Sciences (Grant No. KZCX2-EW-QN207) and the Special Fund for Meteorological Scientific Research in Public Interest (GYHY201306045). The authors would also like to thank two anonymous reviewers for their valuable comments and suggestions that are helpful for improving the paper.

## References

1. Anderson J. L . An ensemble adjustment Kalman filter for data assimilation . Mon. Weather Rev . 2001 ; 129 : 2884 – 2903 .

2. Bauer P. , Lopez P. , Benedetti A. , Salmond D. , Saarinen S. , co-authors . Implementation of 1D+4D-Var assimilation of precipitation affected microwave radiances at ECMWF II: 4D-Var . Q. J. Roy. Meteorol. Soc . 2006 ; 132 : 2307 – 2332 .

3. Beck A. , Ehrendorfer M . Singular-vector-based covariance propagation in a quasigeostrophic assimilation system . Mon. Weather Rev . 2005 ; 133 : 1295 – 1310 .

4. Bormann N. , Thepaut J. N . Impact of MODIS polar winds in ECMWF's 4DVAR data assimilation system . Mon. Weather Rev . 2004 ; 132 : 929 – 940 .

5. Caya A. , Sun J. , Snyder C . A comparison between the 4DVAR and the ensemble Kalman filter techniques for radar data assimilation . Mon. Weather Rev . 2005 ; 133 : 3081 – 3094 .

6. Chen F. , Dudhia J . Coupling an advanced land-surface/hydrology model with the Penn State/NCAR MM5 modeling system. Part I: Model description and implementation . Mon. Weather Rev . 2001 ; 129 : 569 – 585 .

7. Chen S.-H. , Sun W.-Y . A one-dimensional time dependent cloud model . J. Meteorol. Soc. Jpn . 2002 ; 80 : 99 – 118 .

8. Cheng H. , Jardak M. , Alexe M. , Sandu A . A hybrid approach to estimating error covariances in variational data assimilation . Tellus A . 2010 ; 62 : 288 – 297 .

9. Clayton A. M. , Lorenc A. C. , Barker D. M . Operational implementation of a hybrid ensemble/4D-Var global data assimilation system at the Met Office . Q. J. Roy. Meteorol. Soc . 2013 ; 139 : 1445 – 1461 .

10. Courtier P. , Talagrand O . Variational assimilation of meteorological observations with the adjoint vorticity equation II: numerical results . Q. J. Roy. Meteorol. Soc . 1987 ; 113 : 1329 – 1347 .

11. Courtier P. , Thepaut J. N. , Hollingsworth A . A strategy for operational implementation of 4DVar using an incremental approach . Q. J. Roy. Meteorol. Soc . 1994 ; 120 : 1367 – 1387 .

12. Dennis J. E. Jr. , Schnabel R. B . Numerical Methods for Unconstrained Optimization and Nonlinear Equations (Classics in Applied Mathematics) .

13. Dudhia J . Numerical study of convection observed during the winter monsoon experiment using a mesoscale two-dimensional model . J. Atmos. Sci . 1989 ; 46 : 3077 – 3107 .

14. Evensen G . Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics . J. Geophys. Res . 1994 ; 99 ( C5 ): 10143 – 10162 .

15. Evensen G . Sampling strategies and square root analysis schemes for the EnKF . Ocean Dynam . 2004 ; 54 : 539 – 560 .

16. Evensen G. , van Leeuwen P. J . An ensemble Kalman smoother for nonlinear dynamics . Mon. Weather Rev . 2000 ; 128 : 1852 – 1867 .

17. Gaspari G. , Cohn S. E . Construction of correlation functions in two and three dimensions . Q. J. Roy. Meteorol. Soc . 1999 ; 125 : 723 – 757 .

18. Gauthier P. , Tanguay M. , Laroche S. , Pellerin S. , Morneau J . Extension of 3DVAR to 4DVAR: implementation of 4DVAR at the meteorological service of Canada . Mon. Weather Rev . 2007 ; 135 ( 6 ): 2339 – 2354 .

19. Harlim J. , Hunt B. R . Four-dimensional local ensemble transform Kalman filter: numerical experiments with a global circulation model . Tellus A . 2007 ; 59 : 731 – 748 .

20. Hong S.-Y. , Noh Y. , Dudhia J . A new vertical diffusion package with an explicit treatment of entrainment processes . Mon. Weather Rev . 2006 ; 134 : 2318 – 2341 .

21. Houtekamer P. L. , Mitchell H. L . Data assimilation using an ensemble Kalman filter technique . Mon. Weather Rev . 1998 ; 126 : 796 – 811 .

22. Houtekamer P. L. , Mitchell H. L . A sequential ensemble Kalman filter for atmospheric data assimilation . Mon. Weather Rev . 2001 ; 129 : 123 – 137 .

23. Hunt B. R. , Kostelich E. J. , Ott E. , Szunyogh I . Efficient data assimilation spatiotemporal chaos: a local ensemble transform Kalman filter . Physica D . 2007 ; 230 : 112 – 126 .

24. Keppenne C. L . Data assimilation into a primitive-equation model with a parallel ensemble Kalman filter . Mon. Weather Rev . 2000 ; 128 : 1971 – 1981 .

25. Le Dimet F. X. , Talagrand O . Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects . Tellus A . 1986 ; 38 : 97 – 110 .

26. Lewis J. M. , Derber J. C . The use of the adjoint equation to solve a variational adjustment problem with advective constraints . Tellus A . 1985 ; 37 : 309 – 322 .

27. Lin Y.-L. , Farley R. D. , Orville H. D . Bulk parameterization of the snow field in a cloud model . J. Clim. Appl. Meteorol . 1983 ; 22 : 1065 – 1092 .

28. Liu C. , Xiao Q. , Wang B . An ensemble-based four-dimensional variational data assimilation scheme. Part I: Technique formulation and preliminary test . Mon. Weather Rev . 2008 ; 136 : 3363 – 3373 .

29. Lorenc A . The potential of the ensemble Kalman Filter for NWP: a comparison with 4Dvar . Q. J. Roy. Meteorol. Soc . 2003 ; 129 : 3183 – 3203 .

30. Ly H. V. , Tran H. T . Modeling and control of physical processes using proper orthogonal decomposition . Math. Comput. Model . 2001 ; 33 : 223 – 236 .

31. Ly H. V. , Tran H. T . Proper orthogonal decomposition for flow calculations and optimal control in a horizontal CVD reactor . Q. Appl. Math . 2002 ; 60 : 631 – 656 .

32. Milewski T. , Bourqui M. S . Potential of an ensemble Kalman smoother for stratospheric chemical-dynamical data assimilation . Tellus A . 2013 ; 65 : 18541 .

33. Mlawer E. J. , Taubman S. J. , Brown P. D. , Iacono M. J. , Clough S. A . Radiative transfer for inhomogeneous atmosphere: RRTM, a validated correlated-k model for the longwave . J. Geophys. Res . 1997 ; 102 ( D14 ): 16663 – 16682 .

34. Ott E. , Hunt B. R. , Szunyogh I. , Zimin A. V. , Kostelich E. J. , co-authors . A local ensemble Kalman filter for atmospheric data assimilation . Tellus A . 2004 ; 56 : 415 – 428 .

35. Pan X. , Tian X. , Li X. , Xie Z. , Shao A. , co-authors . Assimilating Doppler radar radial velocity and reflectivity observations in the weather research and forecasting model by a proper orthogonal decomposition based ensemble three-dimensional variational assimilation method . J. Geophys. Res . 2012 ; 117 : D17113 .

36. Park K. , Zou X . Toward developing an objective 4DVAR BDA scheme for hurricane initialization based on TPC observed parameters . Mon. Weather Rev . 2004 ; 132 ( 8 ): 2054 – 2069 .

37. Qiu C. , Shao A. , Xu Q. , Wei L . Fitting model fields to observations by using singular value decomposition: an ensemble-based 4DVar approach . J. Geophys. Res . 2007 ; 112 : D11105 .

38. Rosmond T. , Xu L . Development of NAVDAS-AR: nonlinear formulation and outer loop tests . Tellus A . 2006 ; 58 ( 1 ): 45 – 58 .

39. Rutledge S. A. , Hobbs P. V . The mesoscale and microscale structure and organization of clouds and precipitation in midlatitude cyclones. XII: a diagnostic modeling study of precipitation development in narrow cloud-frontal rain bands . J. Atmos. Sci . 1984 ; 20 : 2949 – 2972 .

40. Skamarock W. C. , Klemp J. B. , Dudhia J. , Gill D. O. , Barker D. M. , co-authors . A Description of the Advanced Research WRF Version 3 .

41. Tian X. , Xie Z . Implementations of a square-root ensemble analysis and a hybrid localisation into the POD-based ensemble 4DVar . Tellus A . 2012 ; 64 18735 .

42. Tian X. , Xie Z. , Dai A . An ensemble-based explicit four-dimensional variational assimilation method . J. Geophys. Res . 2008 ; 113 : D21124 .

43. Tian X. , Xie Z. , Dai A. , Jia B. , Shi C . A microwave land data assimilation system: scheme and preliminary evaluation over China . J. Geophys. Res . 2010 ; 115 : D21113 .

44. Tian X. , Xie Z. , Dai A. , Shi C. , Jia B. , co-authors . A dual-pass variational data assimilation framework for estimating soil moisture profiles from AMSR-E microwave brightness temperature . J. Geophys. Res . 2009 ; 114 : D16102 .

45. Tian X. , Xie Z. , Liu Y. , Cai Z. , Fu Y. , co-authors . A joint data assimilation system (Tan-Tracker) to simultaneously estimate surface CO2 fluxes and 3-D atmospheric CO2 concentrations from observations . Atmos. Chem. Phys . 2014 ; 14 : 13281 – 13293 .

46. Tian X. , Xie Z. , Sun Q . A POD-based ensemble four-dimensional variational assimilation method . Tellus A . 2011 ; 63 : 805 – 816 .

47. Wang B. , Liu J. , Wang S. , Cheng W. , Juan L. , co-authors . An economical approach to four-dimensional variational data assimilation . Adv. Atmos. Sci . 2010 ; 27 ( 4 ): 715 – 727 .

48. Wang X. , Barker D. M. , Snyder C. , Hamill T. M . A hybrid ETKF-3DVAR data assimilation scheme for the wrf model. Part I: Observing system simulation experiment . Mon. Weather Rev . 2008a ; 136 : 5116 – 5131 .

49. Wang X. , Barker D. M. , Snyder C. , Hamill T. M . A hybrid ETKF-3DVAR data assimilation scheme for the WRF model. Part II: Real observation experiments . Mon. Weather Rev . 2008b ; 136 : 5132 – 5147 .

50. Wang X. , Parrish D. , Kleist D. , Whitaker J. S . GSI 3DVar-based ensemble-variational hybrid data assimilation for NCEP global forecast system: single resolution experiments . Mon. Weather Rev . 2013 ; 141 : 4098 – 4117 .

51. Whitaker J. S. , Hamill T. M . Ensemble data assimilation without perturbed observations . Mon. Weather Rev . 2002 ; 130 : 1913 – 1924 .

52. Xu Q . Generalized adjoint for physical processes with parameterized discontinuities. Part I: Basis issues and heuristic examples . J. Atmos. Sci . 1996 ; 53 ( 8 ): 1123 – 1142 .

53. Zhang B. , Tian X. , Sun J. , co-authors . PODEn4DVar-based radar data assimilation scheme: formulation and preliminary results from real-data experiments with advanced research WRF (ARW) . Tellus A . 2014

54. Zhang F. , Snyder C. , Sun J . Tests of an ensemble Kalman filter for convective-scale data assimilation: impact of initial estimate and observations . Mon. Weather Rev . 2004 ; 132 : 1238 – 1253 .

55. Zhang F. Q. , Zhang M. , Hansen J. A . Coupling ensemble Kalman filter with four dimensional variational data assimilation . Adv. Atmos. Sci . 2009 ; 26 ( 1 ): 1 – 8 .

56. Zhang M. , Zhang F. Q . E4DVar: coupling an ensemble Kalman filter with four-dimensional variational data assimilation in a limited-area weather prediction model . Mon. Weather Rev . 2012 ; 140 : 587 – 600 .

57. Zou X. , Kuo Y.-H . Rainfall assimilation through an optimal control of initial and boundary conditions in a limited-area mesoscale model . Mon. Weather Rev . 1996 ; 124 : 2859 – 2882 .