Use of Incremental Analysis Updates in 4D-Var Data Assimilation

Received 13 February 2015; revised 9 June 2015; accepted 10 June 2015

A prototype incremental strong constraint four-dimensional variational (4D-Var) data assimilation system has been developed at the Goddard Space Flight Center (GSFC) Global Modeling and Assimilation Office (GMAO) (Tremolet and Todling, personal communication). This prototype uses the Grid-Point Statistical Interpolation developed in the National Centers for Environmental Prediction (Wu et al., 2002), the fifth generation of the Goddard Earth Observing System (GEOS-5) atmospheric model (Rienecker et al., 2008), and the GEOS-4 based tangent linear and adjoint models (Giering et al., 2003).

The 4D-Var prototype is superior to the NASA/GFSC 3D-Var system for a number of reasons. A key benefit is the use of the tangent linear and its adjoint model, which enables flow-dependent propagation of background error covariance over a given assimilation window. Preliminary experiments made with the prototype 4D-Var system showed improvements over 3D-Var. For instance, the forecast skill scores were improved, especially in the data-sparse Southern Hemisphere.

However, a verification study of the 4D-Var prototype found that large spikes exist in the time series of global mean precipitation (Fig. 1). While the exact cause of this has yet to be determined, it most likely arises from discontinuities at the beginning of data assimilation windows due to using initial condition increments as control variables. Such discontinuities can trigger spurious gravity waves that impact precipitation and other physics fields.

(Bloom et al., 1996) introduced an incremental analysis updating (IAU) technique to the NASA/GSFC 3D-Var analysis system to remove the same type of initial discontinuities in 3D-Var. (Polavarapu et al., 2000) also found a positive impact from using digital filter initialization to suppress gravity waves and to smooth analyses in a 4D-Var system, and (Polavarapu et al., 2004) later showed that IAU and incremental digital filtering are equivalent to each other for linear models using time-invariant coefficients. Because the NASA/GSFC 4D-Var prototype was based on the 3D-Var system with an IAU component embedded in the nonlinear forward model, it is quite straightforward to implement the IAU technique in the 4D-Var system.

The methodology underlying the current study has roots in advancements long before the development of operational 4D-Var. (Derber, 1989) proposed variational continuous assimilation (VCA) as a modification to the adjoint technique introduced by (Le Dimet and Talagrand, 1986). In the adjoint technique, the initial conditions are control variables that are adjusted to find the best fit to the data. Meanwhile, in VCA the control variables are the values that alter the model time derivatives over the assimilation interval. The operational implementation of 4D-Var in terms of increments was later proposed by (Courtier et al., 1994).

Figure 1.  Time series of globally averaged hourly total precipitation (units: mm d^-1). The red curve is for the conventional 4D-Var system, and the blue curve shows the IAU-based 4D-Var system.

In the current study the initial condition increments in conventional 4D-Var were replaced with time-integrated increments as control variables in the cost function, which is similar to the VCA technique. This new IAU-based 4D-Var system distributes analysis increments throughout the entire assimilation window both in the full model outer loops and in the tangent linear and the adjoint models. By working this way, a continuous analysis state is generated, and initial discontinuities are greatly reduced. The characteristics of the IAU-based 4D-Var solution could be quite different from those of the conventional 4D-Var without IAU since the latter incrementally adjusts the time derivatives or tendencies throughout the analysis window while the former adjusts only the initial condition at the beginning of the window. Theoretically, IAU-based 4D-Var is also superior to IAU-based 3D-Var, which does not take into account the propagation of the analysis increment in the analysis inner loop.

The outline of this paper is as follows: Section 2 describes both conventional 4D-Var and IAU-based 4D-Var. Solution convergence and computer costs are discussed in section 3. In section 4, some experimental results are presented, and the positive impacts of the IAU technique in the 4D-Var system are evaluated. A summary and discussion points are presented in section 5.

In this section, the algorithm employed by the NASA GMAO incremental strong constraint 4D-Var prototype is briefly described in section 2.1, and the generalized IAU-based 4D-Var system is presented in section 2.2. In section 2.3, two special cases of the generalized IAU-based 4D-Var system are discussed, one of which is simplified to the regular strong constraint 4D-Var, and the other is a strict IAU 4D-Var.

In an incremental, strong-constraint 4D-Var system (Courtier et al., 1994), the model's trajectory and the departure of observations from the model are calculated by the use of a non-linear model in the outer-loop, \begin{equation} \bm{x}_i=\bm{M}_{i,i-1}\bm{x}_{i-1} , \end{equation} where x is the model state variable, M is the nonlinear forecast model, and the time step i is 0,1,2,…,n. The preconditioned Lanczos conjugate-gradient method (Fisher, 1998) is used to solve the minimization problem with the cost-function in the inner loop written as \begin{equation} J_0=\dfrac{1}{2}(\bm{\delta}\bm{x}_0+\bm{b})^{\rm T}\bm{B}^{-1}(\bm{\delta}\bm{x}_0+\bm{b})+\sum_{i=0}^n\bm{y}_i^{\rm T}\bm{R}_i^{-1}\bm{y}_i , \end{equation} where the control variable δx₀ is the correction from first guess x_g, and b=x_g-x_b is the difference between x_g and model state x_b, B is the background error covariance matrix, R _i is observation error, and y_i is the observation. The ( ) ^T notation denotes a transpose of a vector or a matrix. For i=0,1,…,n, $$ \bm{y}_i=\bm{H}_i\bm{M}_{i,0}\bm{\delta}\bm{x}_0+\bm{d}_i , $$ where H_i is the linearized observation operator, M_i,0 is the tangent linear model $$ \bm{M}_{i,0}=\Pi_{j=0}^{i-1}\bm{M}_{j+1,j}\equiv\bm{M}_{i,i-1}\bm{M}_{i-1,i-2}\cdots\bm{M}_{2,1}\bm{M}_{1,0} , $$ and the gradient of J₀ with respect to δx₀ is \begin{equation} \nabla J_0=\bm{B}^{-1}(\bm{\delta}\bm{x}_0+\bm{b})+\sum_{i=0}^n\bm{M}_{i,0}^{\rm T}\bm{H}_i^{\rm T}\bm{R}_i^{-1}\bm{y}_i . \end{equation}

In an actual implementation, the forward integrations are computed step by step for i=0,1,…,n, \begin{eqnarray} \bm{\delta}\bm{x}_i&=&\bm{M}_{i,i-1}\bm{\delta}\bm{x}_{i-1} ,\\ y_i&=&\bm{H}_i\bm{\delta}\bm{x}_i+d_{i} , \end{eqnarray} and the summation on the right-hand side of ∇ J₀ is derived backward through adjoint operations for i=n-1,…,1, \begin{eqnarray} \bm{\delta}\bm{x}_n^*&=&\bm{H}_n^{\rm T}\bm{R}_n^{-1}\bm{y}_n ,\\ \bm{\delta}\bm{x}_i^*&=&\bm{M}_{i+1,i}^{\rm T}\bm{\delta}\bm{x}_{i+1}^*+\bm{H}_i^{\rm T}\bm{R}_i^{-1}\bm{y}_i . \end{eqnarray} Then, at i=0, \begin{equation} \bm{\delta}\bm{x}_0^*=\bm{M}_{1,0}^{\rm T}\bm{\delta}\bm{x}_1^*+\bm{H}_0^{\rm T}\bm{R}_0^{-1}y_0\equiv\sum_{i=0}^n\bm{M}_{i,0}^{\rm T}\bm{H}_i^{\rm T}\bm{R}_i^{-1}\bm{y}_i . \end{equation}

In an IAU-based 4D-Var system, analysis increments are gradually inserted to the model at each time step instead of at the initial step. Noting the time-integrated increments \(\Delta\overlinex\) as the new control variables, and assuming the tendencies \(\alpha_i\Delta\overlinex\) are small fractions of the total integrated increments at time step i, we can write the nonlinear forward model in the outer loop as \begin{equation} \bm{x}_i=\bm{M}_{i,i-1}\bm{x}_{i-1}+\alpha_i\Delta\overline{\bm{x}} , \end{equation} with ∑_i=0ⁿα_i=1, which is implemented similarly to that of the IAU-based 3D-Var system.

In the inner loop of the IAU based 4D-Var system, the forward operations, for i=1,…,n, are \begin{eqnarray} \bm{\delta}\bm{x}_0&=&\alpha_0\Delta\overline{\bm{x}} ,\\ \bm{\delta}\bm{x}_i&=&\bm{M}_{i,i-1}\bm{\bm{\delta}}\bm{x}_{i-1}+\alpha_i\Delta\overline{\bm{x}} , \end{eqnarray} and the adjoint operations in a backward order for i=n-1,…,1,0, are \begin{eqnarray} \bm{\delta}\bm{x}_n^*&=&\bm{H}_n^{\rm T}\bm{R}_n^{-1}\bm{y}_n ,\quad \Delta\overline{\bm{x}}_n^*={\bf 0} ,\\ \bm{\delta}\bm{x}_i^*&=&\bm{M}_{i+1,i}^{\rm T}\bm{\delta}\bm{x}_{i+1}^*+\bm{H}_i^{\rm T}\bm{R}_i^{-1}\bm{y}_i ,\quad \Delta\overline{\bm{x}}_i^*\!=\!\alpha_{i+1}\bm{\delta x}_{i+1}^*+\Delta\overline{\bm{x}}_{i+1}^* ,\nonumber\\ \end{eqnarray} with the adjoint \(\Delta\overline\bmx^*\) of the control variables \(\Delta\overlinex\) being written as \begin{equation} \Delta\overline{\bm{x}}^*=\alpha_0\bm{\delta}\bm{x}_0^*+\Delta\overline{\bm{x}}_0^*=\sum_{i=0}^n\alpha_i\bm{\delta}\bm{x}_i^* , \end{equation}

The corresponding cost function in terms of \(\Delta\overlinex\) is \begin{equation} J_0=\dfrac{1}{2}(\Delta\overline{\bm{x}}+\bm{b})^{\rm T}\bm{B}^{-1}(\Delta\overline{\bm{x}}+\bm{b})+\sum_{i=0}^n\bm{y}_i^{\rm T}\bm{R}_i^{-1}\bm{y}_i , \end{equation} where for i=0,1,…,n \begin{equation} \bm{y}_i=\widetilde{\bm{H}}_i\widetilde{\bm{M}}_{i,0}\widetilde{\bm{I}}_0\Delta\overline{\bm{x}}+\bm{d}_i . \end{equation} The gradient of J₀ with respect to \(\Delta\overlinex\) is \begin{equation} \nabla J_0=\bm{B}^{-1}(\Delta\overline{\bm{x}}+\bm{b})+\widetilde{\bm{I}}_0^{\rm T}\sum_{i=0}^n\widetilde{\bm{M}}_{i,0}^{\rm T}\widetilde{\bm{H}}_i^{\rm T}\widetilde{\bm{R}}_i^{-1}\bm{y}_i , \end{equation} where all the symbols with tildes are redefined to reflect the presence of IAU vectors, \begin{eqnarray*} \widetilde{\bm{M}}_{i,0}&=&\prod_{j=0}^{i-1}\widetilde{\bm{M}}_{j+1,j}\equiv\widetilde{\bm{M}}_{i,i-1}\widetilde{\bm{M}}_{i-1,i-2}\cdots\widetilde{\bm{M}}_{2,1}\widetilde{\bm{M}}_{1,0} ,\\ \widetilde{\bm{M}}_{i,i-1}&=&\left[ \begin{array}{c@{\quad}c} \bm{M}_{i,i-1} & \alpha_i\\ 0, & 1 \end{array} \right] ,\quad \widetilde{\bm{H}}_i=[\bm{H}_i\quad 0] ,\quad \widetilde{\bm{I}}_0=\left[ \begin{array}{c} \alpha_0\\1 \end{array} \right] . \end{eqnarray*} Here, \(\bm\delta\widetilde\bmx_i\), the intermediate vectors with augmented IAU, and their adjoint vectors \(\bm\delta\widetilde\bmx_i^*\) at i=0,1,…,n, are written as \begin{eqnarray*} \bm{\delta}\widetilde{\bm{x}}_i=\left[ \begin{array}{c} \bm{\delta}\bm{x}_i\\ \Delta\overline{\bm{x}}_i \end{array} \right] ,\quad \bm{\delta}\widetilde{\bm{x}}_i^*=\left[ \begin{array}{c} \bm{\delta}\bm{x}_i^*\\ \Delta\overline{\bm{x}}_i^* \end{array} \right] .\\[-2mm] \end{eqnarray*} Thus, the forward tangent linear model operations for time i=1,…,n-1,n, are \begin{eqnarray} &&\hspace*{-10mm}\delta\widetilde{\bm{x}}_0\left[ \begin{array}{c} \alpha_0\\ 1 \end{array} \right]\Delta\overline{\bm{x}}=\widetilde{\bm{I}}_0\Delta\overline{\bm{x}} ,\\ \bm{\delta}\widetilde{\bm{x}}_i&=&\left[ \begin{array}{c@{\quad}c} \bm{M}_{i,i-1} & \alpha_i\\ 0 & 1 \end{array} \right]\bm{\delta}\widetilde{\bm{x}}_{i-1}=\widetilde{\bm{M}}_{i,i-1}\bm{\delta}\widetilde{\bm{x}}_{i-1}\nonumber\\ &=&\widetilde{\bm{M}}_{i,0} \bm{\delta}\widetilde{\bm{x}}_0=\widetilde{\bm{M}}_{i,0}\widetilde{\bm{I}}_0\Delta\widetilde{\bm{x}} . \end{eqnarray} The backward adjoint operations for i=n-1,…,1, are \begin{eqnarray} \bm{\delta}\widetilde{\bm{x}}_n^*&=&\widetilde{\bm{H}}_n^{\rm T}\widetilde{\bm{R}}_n^{-1}\bm{y}_n ,\\ \bm{\delta}\widetilde{\bm{x}}_i^*&=&\widetilde{\bm{M}}_{i+1,i}^{\rm T}\bm{\delta}\widetilde{\bm{x}}_{i+1}^*+\widetilde{\bm{H}}_i^{\rm T}\bm{R}_i^{-1}\bm{y}_i , \end{eqnarray} For i=0, \begin{equation} \bm{\delta}\widetilde{\bm{x}}_0^*=\widetilde{\bm{M}}_{1,0}^{\rm T}\bm{\delta}\widetilde{\bm{x}}_1^*+\widetilde{\bm{H}}_0^{\rm T}\bm{R}_0^{-1}\bm{y}_0 \equiv\sum_{i=0}^n\widetilde{\bm{M}}_{i,0}^{\rm T}\widetilde{\bm{H}}_i^{\rm T}\bm{R}_i^{-1}\bm{y}_i . \end{equation} Equations (18-22) for IAU-based 4D-Var are very similar to Eqs. (4) and (6-8), which allows for easy implementation of Eqs. (18-22).

The differences between the conventional incremental, strong-constraint 4D-Var and IAU-based 4D-Var are illustrated in Fig. 2. In conventional 4D-Var, analysis increments are added at the initial step, which creates discontinuities at the beginning of an assimilation window. On the other hand, IAU-based 4D-Var gradually adds small fractions of analysis increments to a nonlinear forecast model in the outer loops at each time step and acts like a continuous assimilation method. The IAU-based 4D-Var system also gradually adds small fractions of increments to the tangent linear model in the inner loops.

Figure 2.  Schematic diagram of both conventional 4D-Var and the IAU-based 4D-Var system. In conventional 4D-Var, analysis increments are added at initial time steps, where the discontinuities exist. In the IAU-based 4D-Var, the analysis increments are gradually added to the nonlinear forecast model in the outer loops. The IAU-based 4D-VAR system also gradually adds increments to the tangent linear model and its adjoint in the inner loops.

In data assimilation a weak constraint refers to corrections made to the forecast model rather than the initial condition, so IAU-based 4D-Var is actually a weak-constraint scheme. However, the algorithm is implemented exactly as in the incremental strong-constraint algorithm. The only difference is that time-integrated increments replace initial condition increments as control variables. As for computational overhead, the increase is negligible for IAU-based 4D-Var when compared with conventional 4D-Var.

The generalized IAU-based 4D-Var can be simplified to the conventional incremental 4D-Var when $$ \alpha_0=1,\quad {\rm and} \alpha_i=0 {\rm for} i=1,2,\cdots,n . $$ For another extreme case, the generalized algorithm reduces to a strict IAU-based 4D-Var when $$ \alpha_0=0,\quad {\rm and} \alpha_i=\dfrac{1}{n} {\rm for} i=1,2,\cdots,n. $$ In this study, we only run experiments for the above two special cases to evaluate the impacts of the IAU in a 4D-Var system.

To validate the algorithm, the first step is to examine solution convergence. We ran the outer loop 10 times and fixed the number of inner loop iterations at 70 for both conventional 4D-Var and strict IAU-based 4D-Var. Both the conventional 4D-Var experiment and IAU-based algorithm were carried out with version fdda_b1p3-das_215-3 of the NASA/GFSC 4D-Var prototype. The default configuration was adopted for setting up the control experiment, except that the inner loop iteration numbers and outer loop numbers were changed. Figure 3 shows the evolution of the cost per observation (Jo/p) as a function of the accumulated number of inner loop iterations for 10 outer loops. The convergence of Jo/p is obtained after only three outer loop iterations, though the convergence values are quite different between the two experiments. Meanwhile, Fig. 4 shows the reduction of gradient norms as a function of inner-loop iteration number for the first three outer loops. The gradient norm from IAU-based 4D-Var decreases much more quickly than that of the conventional 4D-Var, which shows that the IAU-based algorithm converges faster. However, IAU-based 4D-Var converges to a bigger minimum value (0.470) than the conventional 4D-Var system (0.389). This implies that conventional 4D-Var has less model constraint and more strongly fits observations, but this could be a less balanced solution. IAU-based 4D-Var has more constraint from the model and a more balanced solution, but it does not fit observations as well.

Figure 3.  The evolution of cost per observation (Jo/p) as a function of accumulated inner loop iteration numbers. The solid lines represent Jo/p as computed in the nonlinear trajectory runs, while the dashed lines represent Jo as computed in the inner loop. The red color is for conventional 4D-VAR, and the blue is for IAU-based 4D-VAR.

Figure 4.  The J gradient reduction as a function of inner loop iteration number for the first three outer loops. The red curve is for conventional 4D-VAR, and the blue curve is for the IAU-based 4D-VAR system.

Figure 5.  Initial condition increments of temperature for conventional 4D-VAR at 2100 UTC 8 September 2007 at the 51st model level (a), and on the model surface (c), and time-integrated increments for the 6-h assimilation window starting from 2100 UTC 8 September 2007 at the 51st model level (b), and on the model surface (d).

In this section experiments were designed to diagnose the impacts of IAU on 4D-Var. An assimilation experiment was run for both conventional and IAU-based 4D-Var starting from the same initial condition at 2100 UTC 8 September 2007 and ending at 2100 UTC 31 October 2007. The experiments used three outer loops based on the aforementioned convergence study, but the number of inner loop iterations was cut from 70 to 35 for the second and third outer loops to reduce the computational cost.

Figure 5 shows the temperature control variable as an example comparison between the two systems. Shown in the left panels are the initial condition increments from conventional 4D-Var at the 51st model level (≈5.6 hPa; Fig. 5a) and on the model surface (Fig. 5c) at 2100 8 September 2007. The right-hand panels show time-integrated temperature increments from the IAU-based 4D-Var system in a 6 h assimilation window starting from 2100 8 September 2007 at the 51st model level (Fig. 5b), and on the model surface (Fig. 5d). The difference in temperature increments between conventional 4D-Var and IAU-based 4D-Var is quite small. However, for IAU-based 4D-Var the time-integrated increment at level 51 has more detailed small structure over the Southern Hemisphere. In addition, the biggest increments in conventional 4D-Var are near locations with many rawinsonde observations, while IAU-based 4D-Var has big increments farther away.

The time series of hourly total precipitation averaged over the globe from both conventional 4D-Var and IAU-based 4D-Var are first analyzed to see if introducing IAU into 4D-Var can reduce discontinuities. As shown in Fig. 1, the discontinuities found at the beginning of assimilation windows in the conventional 4D-Var (red curve) are completely smoothed out by the IAU implementation (blue line).

Verification of other state variables is also conducted to check the impacts of IAU in 4D-Var. Figure 6 shows an example time series of background and analysis surface temperature at (30^°N, 120^°E) from conventional 4D-Var and the IAU-based 4D-var systems. For this example, the conventional 4D-Var analysis itself is of high quality, and only small discontinuities appeared at initial time of each assimilation window. The evolutions of background and analysis surface temperature from IAU based 4D-Var are quite close to the conventional 4D-Var.

In contrast, the time series of surface temperature at (30^°N, 117.5^°W) in Fig. 7 demonstrates poor performance of the conventional strong constraint 4D-Var. The model background temperatures are much warmer than observations in many assimilation windows, and the analysis overcorrected the systematic warm bias, which resulted in an unrealistic saw-tooth pattern. The IAU-based 4D-Var distributed the increment over the entire assimilation window without over-correction and with a much smoother analysis, although the warm bias still exists.

Figure 6.  Time series of surface temperature background and analysis at (70^°N, 137.5^°E). Closed red circles and the red solid line are the background and analysis from conventional 4D-Var without IAU (upper panel). Closed blue circles and the blue solid line are the background and analysis from IAU-based 4D-Var (lower panel).

The time series in Fig. 8, which shows temperature at the 51st model level at (30^°N, 120^°E), demonstrates noisier gravity wave effects in conventional 4D-Var than in the IAU-based 4D-Var system. Figure 9 shows a time series of globally averaged hourly temperature analyses from conventional 4D-Var (red) and the IAU-based 4D-Var system (blue) at the 51st model level (top panel) and at the model surface (bottom panel). Spurious high-frequency waves are observed in the conventional 4D-Var, while they are absent from the IAU-based 4D-Var.

The time series of globally averaged 6-hourly temperature analyses from IAU-based 3D-Var (green), conventional 4D-Var (red), and IAU-based 4D-Var (blue) are given in Fig. 10. The IAU-based 4D-Var analysis is even smoother than that of the IAU-based 3D-Var, but the difference between 4D-Var analyses with and without IAU is much smaller than the difference between 3D-Var and 4D-Var.

Figure 7.  As in Fig. 6 but for surface temperature at (30^°N, 117.5^°W).

Figure 8.  As in Fig. 6 but for temperature at the 51st model level (≈5.6 hPa) at (30^°N, 120^°E).

Figure 9.  Time series of globally-averaged hourly temperature analyses from conventional 4D-Var (red) and IAU-based 4D-Var (blue). The upper panel shows the time series at model level 51 (≈5.6 hPa), and the lower panel shows the time series at the model surface.

Figure 10.  Time series of global average 6-hourly temperature analyses for IAU-based 3D-Var (green), conventional 4D-Var without IAU (red), and IAU-based 4D-Var (blue). The upper panel shows the time series at the 51st model level (≈5.6 hPa ), and the lower panel shows the time series at the model surface.

Figure 11.  October 2007 means (solid) and standard deviations (dashed) of the observed-minus-background residuals (top panels) and observed-minus-analysis (bottom panels) for atmospheric water vapor mixing ratio in four different regions for conventional 4D-Var without IAU (red), and IAU based 4D-Var.

To verify IAU-based 4D-Var, the statistics of the observation-minus-background (OmB) departures and observation-minus-analysis (OmA) were calculated for 4D-Var experiments with and without IAU. Giving atmospheric water vapor mixing ratio as an example, Fig. 11 shows the monthly mean bias (solid) and standard deviations (dashed) of OmB residuals (top panels) and OmA residuals (bottom panels) for four different regions. Although the OmA standard deviation from 4D-Var with IAU is slightly larger than 4D-Var without IAU, the OmB mean biases and standard deviations from both conventional 4D-Var and IAU-based 4D-Var are very similar. The IAU-based analysis slightly degrades in the OmA standard deviation statistics because of more constraint by the model. For other variables such as temperature, zonal and meridional winds, the results are similar (figures not shown here).

In this paper we introduced IAU into 4D-Var and compared it to conventional strong constraint 4D-Var without IAU. The implementation of IAU in 4D-Var generates a continuous analysis state, and large spikes in global mean precipitation present in conventional 4D-Var were greatly reduced. The easy conversion to the strong constraint formula from IAU is another outstanding characteristic and is why the IAU-based 4D-Var does not have additional computational cost. The experiment in IAU-based 4D-Var used the same background error from conventional 4D-Var without any tuning. This might also be one of the reasons that the OmA statistics of the IAU-based 4D-Var are slightly larger than for conventional 4D-Var. Tuning of the background error will be carried out in the future to improve the OmA statistics and to obtain better forecast skill scores.