Abstract

Real-size assimilation systems basically rely on estimation theory in which two sources of information, background and observations, are combined in an optimal way. However, the optimality of such large problems is not guaranteed, since they rely on different approximations. The observed values of the subparts of the cost function measuring the distance between the analysis and the two sources of information can be compared to their theoretical statistical expectations. Such a posteriori diagnostics can be used to optimize the statistics of observation and also background errors. Moreover, the expectations of the subparts of the cost function associated with observations are related to the weights of the different sources of observations in the analysis. It is shown that those theoretical statistical expectations are direct by-products of an ensemble of perturbed assimilations.

1. Introduction

Most operational meteorological centers are now using a variational scheme to produce global or limited-area analyses (Lewis and Derber 1985; Talagrand and Courtier 1987; Rabier et al. 2000). Variational schemes allow the use of a large spectrum of observations in a very efficient way since today both the size of the observation set and of the model state vector are quite large. These analysis schemes rely on the theory of linear estimation (Talagrand 1997). In the linear estimation theory, the different pieces of information (i.e., the observations and the background estimate of the state vector provided by a short-range forecast) are given weights that are inversely proportional to their error covariances. However, those error statistics are not perfectly known. This is especially the case for background error statistics. The determination and representation of those error statistics remain a major challenge in assimilation schemes.

A potential way to determine background error statistics is to rely on an ensemble of perturbed analyses. Originally introduced by Evensen (1994), the idea of an ensemble of assimilations has been for example applied by Houtekamer et al. (1996) in the framework of ensemble prediction and by Fisher (2003a) to improve the determination of background error statistics. In spite of their rather high numerical cost, large operational centers tend today to implement such Monte Carlo procedures to specify more realistic background error covariances (Berre et al. 2006; Belo Pereira and Berre 2006).

From a different point of view, Dee (1995) has defined a method based on the maximum likelihood to determine both observation and background error statistics.

Desroziers and Ivanov (2001) have proposed a method based on a consistency criterion originally proposed by Talagrand (1999). Chapnik et al. (2004) have further investigated the properties of the algorithm defined in Desroziers and Ivanov (2001) and have in particular shown that it was equivalent to a maximum likelihood approach. Chapnik et al. (2006) have applied this algorithm in an operational framework to determine observation error statistics. Such a tuning procedure has also been used by Buehner et al. (2005).

A specific feature of the previous algorithm is that it is based on a randomization approach that relies on a perturbation of observations. This is also the case in an ensemble of perturbed analyses, as it has been, for example, implemented by Fisher (2003a) or Berre et al. (2006). Such an ensemble of analyses is now operational at Météo-France in order to specify flow-dependent modulations of the background error variances in the operational four-dimensional variational data assimilation (4D-Var) assimilation scheme.

An objective of this paper is to show that an ensemble of perturbed analyses also allows the computation of consistency diagnostics and then the determination of more exact observation and background error variances (Berre et al. 2007).

Another point is that such consistency diagnostics are related to quantities that appear to be the respective sensitivities of the analysis to the different subsets of observations. A second objective of the paper is then to show that an ensemble of perturbed analyses also provides, without any significant additional numerical cost, the weights of the different subsets of observations in an assimilation cycle.

In the next section, the ensemble approach is introduced. In section 3, the consistency diagnostics and the way they can be computed in an ensemble of assimilations are given. Section 4 explains how the analysis sensitivities to the observations can also be seen as by-products of an ensemble of assimilations. The ensemble of perturbed analyses implemented at Météo-France in the global Action de Recherche Petite Echelle Grande Echelle (ARPEGE) model is described in section 5. The diagnostic of observation and background error statistics in the ARPEGE ensemble is shown in section 7. The analysis sensitivities to the different sets of observations are finally shown in the same operational framework. Conclusion and perspectives are given in section 9.

2. Ensemble formalism

In an assimilation cycle, the background xb is given by the evolution of the previous analysis xa by the forecast model M. The subsequent analyzed state xa is obtained as an optimal combination of the background and the observations yo. The two forecast and analysis steps are written as

 
formula
 
formula

where A stands for the possibly nonlinear analysis operator.

The estimate xa can be classically obtained as the solution of the minimization of the following cost function:

 
formula

where Jb(x) and Jo(x) are the background and observation terms, respectively. Matrices 𝗕 and 𝗥 stand for the background and observation error covariance matrices, respectively, and H is the possibly nonlinear observation operator including model integration in the 4D-Var formalism. Such a cost function can be solved by the minimization of a series of quadratic cost function with observation operators linearized around successive trajectories (Courtier et al. 1994).

It can be shown (see appendix A) that the following two equations stand for the evolution of the different errors involved in a forecast/analysis scheme, even in a slightly nonlinear analysis scheme such as 4D-Var:

 
formula
 
formula

where 𝗠 and 𝗛 are linearized versions of the model M and the observation operator H, respectively, and

 
formula

Vectors εb, εa, εm, and εo contain background, analysis, model, and observation errors, respectively, and εa− is the analysis error at the previous analysis step [such a linearization of the analysis step can also be found in Fisher and Andersson (2001), Zagar et al. (2005), or more recently in Daescu (2008)].

As underlined in Berre et al. (2006), this derivation of background and analysis errors holds even if the analysis is not perfect. The aim of an ensemble of perturbed forecasts–analyses is to mimic the evolution of the estimation errors appearing in the previous equations. For a given member l of an ensemble of L forecasts–analyses, the evolution of the perturbations applied to the unperturbed states satisfy the following equations:

 
formula
 
formula

where δlb and δla are the perturbations applied to the unperturbed background xb and analysis xa, respectively. The vector δlm is a vector of perturbations consistent with the covariance matrix 𝗤 of model errors:

 
formula

where ηlm is a vector of n random numbers (where n is the size of the state vector) with a standard Gaussian distribution (mean 0 and variance 1). Similarly, the vector δlo is a vector of perturbations consistent with the covariance matrix 𝗥 of observation errors:

 
formula

where ηlo is a vector of p random numbers (where p is the number of observations) with a standard Gaussian distribution.

Because the equations for the evolution of the perturbations are the same as those for the state estimation errors, the two equations for the perturbations, δlb and δla, will allow perturbations to evolve in a similar way as estimation errors. Then, the covariances of the perturbations will document the covariances of the estimation errors. In practice, this happens to be a powerful and relatively low-cost method when compared to the extended Kalman filter scheme, to obtain approximations of the background and analysis error covariance matrices.

Another point is that, although ensemble Kalman filter implementations most often rely on a relatively large ensemble (typically with 100 elements), some useful information on background error covariances can be obtained with small ensembles (5–10 elements). In that case a spatial filtering of the noisy raw variances has to be applied (Berre et al. 2007; Raynaud et al. 2008). Pannekoucke et al. (2007) also showed that the wavelet approach could be an adapted technique to filter the noisy correlations derived from a small ensemble. The wavelet approach to represent correlations in the 𝗕 matrix was first proposed by Fisher and Andersson (2001).

3. A posteriori diagnostics

a. Analysis consistency diagnostics

As in Talagrand (1997) and Desroziers and Ivanov (2001), it is possible to introduce an extended vector of observations z, combining the proper observations yo, with dimension p, and the background vector xb, with the same dimension n as the unknown true state xt.

This is written as

 
formula

and

 
formula

Then, it is possible to write z under the form:

 
formula

where Γ is an extended observation operator and ε is the vector of background and observation errors with dimension n + p.

An important property pointed out by Talagrand (1999) is that if Ji stands for a subterm of J with pi elements, then, the statistical expectation of Ji(xa), in an optimal system, is

 
formula

where Γj and 𝗦i stand for the observation operator and the error covariance matrix, respectively, associated with those pi elements [a proof of Eq. (14) is given in appendix B for a slightly nonlinear assimilation scheme such as 4D-Var]. The matrix 𝗣a is here defined by

 
formula

In particular, for Γi = 𝗜n and 𝗦i = 𝗕, and knowing that 𝗞 = 𝗣a𝗛T𝗥−1, it follows that

 
formula

and for Γi = 𝗛 and 𝗦i = 𝗥,

 
formula

and then

 
formula

This means that, if the background and error statistics are correctly specified, then the expectation of the global cost function at its minimum should be equal to p. As pointed out by Bennett et al. (1993) and Talagrand (1999), this is a simple a posteriori consistency criterion of the analysis scheme. A deviation of the expectation of J(xa) from p indicates a misspecification in the statistics of the innovation vector.

b. Practical computation in an ensemble of assimilations

In practice, the statistical expectation of J(xa) often differs from p in an operational system, which means that there may be a misspecification in either the background or observation error covariances, or both, but in any case this simple criterion is not sufficient to indicate where the misspecification lies. To go further, one would like to compute the expectations E{Ji(xa)}, for any subpart of J, to examine where the deviation of the statistics of Ji(xa) from their theoretical values are the largest.

The direct use of Eq. (14) to derive this theoretical statistical expectation of Ji(xa) cannot be use in practice, since, in particular, 𝗣a is not explicitly computed in a large size variational analysis scheme.

However, as presented in Desroziers and Ivanov (2001), it can be shown that if the observation errors of a subset of observations i are uncorrelated with the remaining part of the observation errors, then

 
formula

where Πi is the projector that enables one to pass from the complete set of observations yo to the particular subset yio.

Similarly,

 
formula

where Πi is the projector that allows one to pass from the complete set of the elements of xb to a subset xib of elements, the errors of which are uncorrelated with the other errors of xb. Note that Eqs. (19) and (20) are still valid in a slightly nonlinear assimilation scheme such as 4D-Var [see appendix B for the proof of the generic expression in Eq. (14)].

Still following Desroziers and Ivanov (2001), it can be shown (see appendix C) that an estimation of E{Jio(xa)} is given by

 
formula
 
formula

where δlo = 𝗥1/2ηlo is a vector of perturbations of the complete set of observations; 𝗥i and 𝗛i are the error covariance matrix and the observation operator, respectively, associated with the observation subset i; and δlo,i = Πiδlo is the projection of δlo onto the subset of observations i.

Now, because δla = (𝗜 − 𝗞𝗛)δlb + 𝗞δlo in an ensemble of assimilations and also because the δlb and δlo perturbations are uncorrelated, a Monte Carlo estimation of E{Jio(xa)} is given by

 
formula

Similarly, a Monte Carlo estimation of a subpart of E{Jib(xa)} is given by

 
formula

where δb = 𝗕1/2ηb is a vector of perturbations of the complete set of elements of xb, 𝗕i is the background error covariance matrix associated with a subset i of elements of εb uncorrelated with the other elements of εb, and is the projection of δlb onto the subset i [see appendix C for a proof of Eqs. (23) and (24)].

Therefore, a key point of the present paper is that, since an ensemble of perturbed assimilations is based on perturbations of observations δlo and background δlb, then, such an ensemble should provide a way to evaluate the previous theoretical statistical expectations of subparts of the cost function.

Given the ensemble observation perturbations δlo and the corresponding δla on the analysis, it is easy to compute the E{Jio(xa)} since the 𝗥i matrix associated with the is known (by construction δlo = 𝗥1/2ηlo). Most often, 𝗥i is a diagonal matrix in operational systems, which also facilitates the computation.

Because the density of satellite observations is steadily increasing and some of these observations might be correlated, implementations of nondiagonal 𝗥i matrices are also under consideration (Bormann et al. 2003). In that case, 𝗥i is however built under an operator form, which would make the computation of Eq. (23) still possible.

The diagnostic of the E{Jib(xa)} through the use of Eq. (24) is more complicate, since in an ensemble of assimilations the perturbations on the background are not drawn from a given covariance matrix 𝗕, but are rather the result of the evolution of the analysis perturbations by the forecast model. While in the future, an ensemble estimate of may be considered for use, no application of this expression is shown in the present paper. However, the computation of E{Jib(xa)} for the whole Jb part, which is equal to Tr(𝗞𝗛) = Tr(𝗛𝗞), can be obtained indirectly by adding the Tr(Πi𝗛𝗞ΠiT) associated with the different subsets of observations (see section 7b).

c. Diagnostic of error variances

As discussed in Desroziers and Ivanov (2001), the possibility to compute the theoretical statistical expectation of the subparts of J(xa) also enables to optimize the error statistics associated to subsets of vector z. This can be done by rewriting the original cost function under the following form:

 
formula

where and are weighting parameters of background and observation errors, respectively, supposed to be homogeneous over a subset i. The rationale of the tuning procedure is to find the values of those weighting coefficients such as and become close to the respective expected statistical estimations E{Jib(xa)} and E{Jio(xa)} they would have if the system were optimal. This produces

 
formula

and

 
formula

Those equations are nonlinear in and , but their particular forms suggest the use of an iterative fixed-point method to solve them.

Desroziers and Ivanov (2001) and Chapnik et al. (2004) have shown that, in practice, the convergence of the method is very fast and that the values of the and obtained at the first iteration are very close to those obtained after several iterations.

Chapnik et al. (2004) have also shown that the and —verifying the previous relations—are the solutions of the application of a maximum likelihood method to estimate these parameters (Dee and da Silva 1999). Dee and da Silva (1999) and Chapnik et al. (2004) have demonstrated that the ability of the maximum likelihood method—and then of the previous method—to provide correct estimations of the error variances relies on the fact that a part of the innovation vector is spatially correlated and that the other part is not. This hypothesis is also assumed in Hollingsworth and Lönnberg (1986) to evaluate the background and observation variances. In practice, this hypothesis is most often valid, since background errors are spatially correlated and observation errors are not.

Moreover, Chapnik et al. (2004) have shown that observation and background error variances provided by the algorithm are valid even if background error correlations are not correctly specified in matrix 𝗕.

However, in some cases, where observation errors are correlated, such tuning methods will be unable to provide valid estimations of observation error variances.

4. Degree of freedom for signal

An increased amount of research has been done recently in numerical weather prediction to assess the observation impact on analyses and subsequent short-range forecasts. Langland and Baker (2004) have for example proposed a procedure for estimating the impact of observations on a measure of short-range forecast error, using adjoint versions of the forecast model and the data assimilation procedure. A similar approach was followed by Zhu and Gelaro (2008). Desroziers et al. (2005b) have proposed a randomization procedure for evaluating the error variance reduction brought by observations on analyses and forecasts. Trémolet (2008) has also recently introduced a way of computing the adjoint of the assimilation scheme in an incremental variational formalism.

The analysis sensitivity to a particular subset i of observations can also be given by its degree of freedom for signal (DFS) introduced by Rodgers (2000). This quantity is defined by

 
formula

It measures the sensitivity of the analysis to a perturbation of a particular subset of observations and then the weight of these observations in the analysis. Cardinali et al. (2004) have used the DFS to measure the analysis sensitivity to observations in a real size assimilation system. Rabier et al. (2002) have also shown how to use such a diagnostic to select Infrared Atmospheric Sounding Interferometer (IASI) channels in order to extract the useful information from the very large amount of data provided by this instrument.

Fisher (2003b) has investigated the possibility to compute the total DFS brought by the complete set of observations in a real-size data assimilation. He has compared different methods to compute such a quantity. One of them is based on the estimation of the so-called influence matrix 𝗛𝗞 by a randomization procedure, as proposed by Girard (1987) and also by Hutchinson (1989). Such a randomization procedure has also been used by Wahba et al. (1995) to compute the generalized cross-validation criterion and inspired the randomized estimation of the E{Jio} proposed by Desroziers and Ivanov (2001).

Again, since an ensemble of perturbed analyses is based on explicit perturbations of observations and implicit perturbations of the background, it should provide estimations of the DFSi, with nearly no additional computational cost.

Hence, it can be shown (Chapnik et al. 2006) that the previous expression of the partial DFS can be rewritten as

 
formula

where i is a subset of observations.

One can recognize a part of Eq. (14) of the expectation of a subpart of the cost function. Talagrand (1999) has interpreted this expression as a measure of the contribution of the subset of observations i to the overall precision.

Relying on Eq. (8), the DFSi can be also rewritten, even in a nonlinear analysis,

 
formula

that can be estimated by a Monte Carlo procedure (see):

 
formula

Then, it appears that the analysis sensitivities to the different subsets of observations, through their DFSi, are available as direct by-products of an ensemble of perturbed analyses associated with a real-size assimilation scheme.

5. ARPEGE ensemble of perturbed assimilations

The French operational ARPEGE global model is a spectral model with 60 vertical levels and a stretched geometry that allows a higher resolution over France. With a T538 spectral resolution and a C2.4 stretching factor, the maximum spectral resolution is T1290, which corresponds to about 15 km in physical space.

The assimilation scheme associated with the ARPEGE model is based on an incremental 4D-Var formulation with two outer loops and with T107 and T149 spectral uniform resolution increments, respectively.

The variational formulation allows a wide range of observations to be used. The types of observations used in the ARPEGE assimilation are the following (listed by ARPEGE type number and code):

  1. Synop (surface observations),

  2. Aircraft,

  3. Satwind (satellite-derived winds),

  4. Driftbuoy,

  5. Radiosounding,

  6. Profiler,

  7. Brightness temperature (“Radiance”),

  8. Scatterometer, and

  9. GPSRO (GPS radio occultation).

The representation of matrix 𝗕 in the analysis scheme is based on a multivariate formulation (Derber and Bouttier 1999), that includes a normalization in physical space of the vorticity covariances by a field of standard deviations of vorticity background errors. This allows us to introduce some heterogeneity of the covariances, which especially induces a larger weight of observations in the areas where the vorticity error standard deviations are large, like in extratropical areas over oceans.

To make the fields of vorticity error standard deviations dependent on the meteorological situation and also on the data density, an ensemble of six perturbed analyses is now run operationally every 6 h, like for the operational 4D-Var. To reduce the cost of such an ensemble, the perturbed analyses are performed with a 3D-Var first guess at observation time (FGAT) scheme instead of 4D-Var and at a lower resolution than the unperturbed 4D-Var (T359 spectral uniform resolution for the nonlinear forecasts and T107 spectral resolution for the analysis increment). The 3D-Var FGAT scheme can be seen as a simplified implementation of 4D-Var, where the evolution of the analysis increment is neglected during the assimilation period (see e.g., Munro et al. 2004). Preliminary studies showed that the background error variances given by 4D-Var and 3D-Var FGAT were relatively close.

The set of observations used in the ensemble is the same as in the operational 4D-Var.

Since the perturbed backgrounds are different from ech other and also from the unperturbed operational background, the quality control of observations can lead to different rejections of observations, especially for satellite data, since this quality control is partially based on the amplitudes of the differences between the observations and the background. This is part of the nonlinearity of the forecasting/assimilation system that must be kept in order to simulate the dynamics of the background and analysis errors in a real-size system more closely.

For the time being, no model error perturbations are added. As a consequence the background errors deduced from the ensemble of assimilations are known to show too low variances and have to be amplified.

6. Estimation of Tr(Πi𝗛𝗞ΠiT) in the ARPEGE ensemble

As shown in the previous sections, the Tr(Πi𝗛𝗞ΠiT) are key quantities, since they appear in the E{Jio(xa)} and then in the tuning procedure of the observation error variances, but also in the measure of the analysis sensitivities to the different sets of observations since they are directly equal to the DFSi.

As discussed above, the Tr(Πi𝗛𝗞ΠiT), associated with a subset of observations i, can be evaluated by using the following randomized estimation:

 
formula

Because there is no unperturbed member in the ARPEGE 3D-Var FGAT ensemble of assimilations, the perturbations δla are calculated as the differences between two perturbed analyses xla and , as in Fisher (2003a). Accordingly, is the difference between two perturbed sets of observations and , so that the previous expression becomes

 
formula

with l′ = l + 1, for l = 1 to L – 1 and l′ = 1 for l = L.

Note that the δa perturbations could have been computed as the differences between the perturbed 3D-Var FGAT analyses and the control 4D-Var analysis. However, since the perturbed 3D-Var FGAT analyses are also performed at a lower resolution than the control 4D-Var, such a computation of the δa perturbations would not be consistent (another possibility, which is envisaged, would be to perform an additional unperturbed analysis in the same low-resolution 3D-Var FGAT mode).

Equation (33) includes an additional normalization factor of 2, because of the fact that it involves perturbation differences (instead of single perturbation values). The operation (1/σ1o)2 corresponds to a simple division by the variances of the observation errors used for the generation of the observation perturbations, since observation errors are, up to now, assumed to be uncorrelated in the ARPEGE assimilation. Here, Hi(xla) stands for the evaluation of the perturbed analysis xla for the subset of observations i by using the possibly nonlinear observation operator Hi associated with this particular subset.

In practice, the previous computation is rather easy to implement, because all the necessary elements, , Hi(xla), σio, stand in the archived perturbed observation databases. Furthermore, the scalar product appearing in the equation is simply the computation of the sum, over all the observations of subset i and all the perturbations l, of the product, observation by observation, of the difference between two observation perturbations and the difference between the two corresponding perturbed analyses evaluated at observation locations, normalized by the corresponding observation error variance.

Figure 1 shows the Tr(Πi𝗛𝗞ΠiT) quantities for a single date and three subsets of observations—Synop, Radiosounding, and satellite Radiance—and for the six perturbations l of the ARPEGE ensemble. To make those quantities comparable, they are divided by the corresponding numbers of observations pi.

Fig. 1.

Tr(Πi𝗛𝗞ΠiT)/pi for the different ARPEGE ensemble perturbations and a single date (0000 UTC 1 Jun 2008). (top left) Synop (dotted line with open circles), Radiosounding (dotted line with times signs), Radiance (dotted line with plus signs), and all observations (dashed line). (top right) Synop observations: geopotential Z (dotted line with open circles), zonal wind u (dotted line with times signs), and meridional wind υ (dotted line with plus signs). (bottom left) Radiosounding observations: temperature t (dotted line with open circles), zonal wind u (dotted line with times signs), meridional wind υ (dotted line with plus signs), and specific humidity (dotted line with small filled circles). (bottom right) HIRS Radiance channels: 11 (dotted line with open circles), 12 (dotted line with plus signs), 14 (dotted line with times signs), and 15 (dotted line with small filled circles).

Fig. 1.

Tr(Πi𝗛𝗞ΠiT)/pi for the different ARPEGE ensemble perturbations and a single date (0000 UTC 1 Jun 2008). (top left) Synop (dotted line with open circles), Radiosounding (dotted line with times signs), Radiance (dotted line with plus signs), and all observations (dashed line). (top right) Synop observations: geopotential Z (dotted line with open circles), zonal wind u (dotted line with times signs), and meridional wind υ (dotted line with plus signs). (bottom left) Radiosounding observations: temperature t (dotted line with open circles), zonal wind u (dotted line with times signs), meridional wind υ (dotted line with plus signs), and specific humidity (dotted line with small filled circles). (bottom right) HIRS Radiance channels: 11 (dotted line with open circles), 12 (dotted line with plus signs), 14 (dotted line with times signs), and 15 (dotted line with small filled circles).

The stability of the Tr(Πi𝗛𝗞ΠiT)/pi is high along the perturbations. This is especially the case for the whole set of observations or for the whole set of parameters of a type of observations. Of course, the stability of the randomized quantities is slightly diminished when they are associated with a more limited subset of observations, but it remains acceptable. This is likely due to the fact that the computations are done in a global system where the number of observations is quite large. Moreover, such quantities can be averaged over all the perturbations and also over different dates along a given period, that will still increase the robustness of their estimations.

Therefore, the possibility to compute the Tr(Πi𝗛𝗞ΠiT)/pi quantities over more than one perturbed analyses can be useful to improve their evaluation. The dispersion of the values obtained for the different perturbations can also indicate the confidence to give to these evaluations.

7. Error variance diagnostic in the ARPEGE ensemble

a. Tuning of observation error variances

As indicated in the previous sections, a first application of the estimation of the Tr(ΠiTΠi𝗛i𝗞) quantities is to use them to diagnose and optimize the observation error variances of the different subsets of observations.

Such a diagnostic has been systematically performed for all the observations used in the ARPEGE global assimilation. For some observations, such as Aircraft, Satwind, Radiosounding, and Profiler, the variances specified in matrix 𝗥 depend on the vertical.

Figure 2 shows the tuning coefficient sio diagnosed by the randomization procedure for wind and temperature measurements given by the Radiosounding type. The coefficients have been determined over an 18-day period, from 3 to 20 June 2008. The coefficients indicate that the wind error standard deviations should be kept unchanged in the lowest levels but multiplied by a factor around 0.6 in the highest levels. Note that the diagnosed normalization coefficients are similar for the zonal and meridional components of wind, which is rather expected, but also gives confidence in the ensemble randomized diagnostic. For temperature, the standard deviations should be slightly increased in the lowest and highest levels, but decreased by a factor around 0.6 near the 500-hPa level.

Fig. 2.

Normalization coefficients sio for the Radiosounding observation error standard deviations: (top left) zonal wind, (top right) meridional wind, and (bottom) temperature.

Fig. 2.

Normalization coefficients sio for the Radiosounding observation error standard deviations: (top left) zonal wind, (top right) meridional wind, and (bottom) temperature.

The normalization coefficients for the Aircraft observations (Fig. 3) indicate that the wind error standard deviations should be multiplied by values between 0.6 and 0.8, depending on the levels. Again, the normalization coefficients are very similar for the two components of wind. The values of the normalization coefficient are even lower for temperature than for wind, with values between 0.5 and 0.7.

Fig. 3.

As in Fig. 2, but for the Aircraft observation error standard deviations.

Fig. 3.

As in Fig. 2, but for the Aircraft observation error standard deviations.

The mean normalization coefficients by observation type (Fig. 4) show that there is a global overestimation of the observation error standard deviations in the ARPEGE assimilation. The needed normalization coefficients are around 0.6–0.8 for Synop, Aircraft, Driftbuoy, Radiosounding, Profiler, and GPSRO. They are close to 0.4–0.5 for Satwind, Radiance, and Scatterometer. For the three latter types of observations, the diagnosed large reductions of the observation error standard deviations may indicate a real overestimation of these standard deviations, but may be also due to the presence of error correlation for these types of observations. In that case, it is known that algorithms such as the one proposed by Desroziers and Ivanov (2001) will fail to recover the right level of error variance and will, on the contrary, interpret observation error correlation as background error correlation. A way to check that the normalization coefficients are correct is to introduce them in the analysis and to verify that the new diagnosed normalization coefficients are close to the previous ones.

Fig. 4.

Mean normalization coefficients sio of observation error standard deviations by observation type: 1 = Synop, 2 = Aircraft, 3 = Satwind, 4 = Driftbuoy, 5 = Radiosounding, 6 = Profiler, 7 = Brightness temperature, 9 = Scatterometer, and 10 = GPSRO.

Fig. 4.

Mean normalization coefficients sio of observation error standard deviations by observation type: 1 = Synop, 2 = Aircraft, 3 = Satwind, 4 = Driftbuoy, 5 = Radiosounding, 6 = Profiler, 7 = Brightness temperature, 9 = Scatterometer, and 10 = GPSRO.

b. Tuning of background error variances

As above-mentioned, the use of Eq. (24), to diagnose background variances, is difficult since the ensemble background perturbations are not drawn from a known covariance matrix, as it is, on the contrary, the case for observations.

However, if a single normalization factor of the 𝗕 matrix is sought, only the expectation E{Jb(xa)} is needed for the complete background vector. Moreover, as mentioned in section 3, E{Jb(xa)} should be equal to Tr(𝗛𝗞) in a slightly nonlinear assimilation. The estimation of Tr(𝗛𝗞) can be in turn deduced from the sum over all the observations of the Tr(Πi𝗛𝗞ΠiT) associated with the different subsets of observations.

Table 1 shows the Jb(xa) with unperturbed observations and backgrounds in comparison with the randomized estimation of Tr(𝗛𝗞). Values are given for a particular date for the four 6-h assimilations. The values of Jb(xa) are also shown for 3D-Var FGAT analyses performed with the same unperturbed observations and backgrounds used in the 4D-Var assimilations, as the ensemble of ARPEGE assimilations is also, up to now, performed in 3D-Var FGAT mode.

Table 1.

The Jb(xa) 4D-Var and 3D-Var FGAT, with unperturbed observations and backgrounds, compared with E{Jb(xa)} at 0000, 0600, 1200, and 1800 UTC 1 Jun 2008.

The Jb(xa) 4D-Var and 3D-Var FGAT, with unperturbed observations and backgrounds, compared with E{Jb(xa)} at 0000, 0600, 1200, and 1800 UTC 1 Jun 2008.
The Jb(xa) 4D-Var and 3D-Var FGAT, with unperturbed observations and backgrounds, compared with E{Jb(xa)} at 0000, 0600, 1200, and 1800 UTC 1 Jun 2008.

It can be seen that the values of Jb(xa) are systematically slightly larger for 3D-Var FGAT than for 4D-Var. This is expected since observations are known to have more weight in 3D-Var FGAT than in 4D-Var. Of course, this is not evidence that 3D-Var FGAT performs better than 4D-Var. On the contrary, this may be interpreted as a better filtering of the observation errors in 4D-Var.

The analysis sensitivity to observations is larger at 0000 and 1200 UTC than at 0600 and 1800 UTC, which is simply because radiosondes are mainly taken at 0000 UTC and that the required normalization coefficient sb can be computed with

 
formula

It varies from 0.61 to 0.71 with 4D-Var and from 0.64 to 0.76 with 3D-Var FGAT (Table 1), which means that the average background error standard deviations are too large in the ARPEGE 4D-Var. This says, in turn, that the background is given too less confidence in the assimilation.

c. Impact of the error variance tuning

The tuned variances for observation errors and the diagnosed normalization coefficient sb have been introduced in the 4D-Var scheme to measure the impact of such a tuning on the subsequent short-range forecasts.

Figure 5 shows an example of the root-mean-square (RMS) scores of the forecasts issued from the tuned analyses, compared to the operational scores. The scores are shown here for the 48-h wind forecasts against radiosondes at 250 hPa, for North America and 20°–90°N domains. The scores are relatively neutral globally, although the impact is slightly positive during the 12–17 November period. This relative neutrality is related to the fact that observation and background error variances tend to be reduced here, by chance, by the same amount of tuning.

Fig. 5.

RMS scores of 48-h wind forecasts vs radiosonde observations at 250 hPa: operational scores (solid line); scores with a tuning of observation error variances and a global tuning of background error variances (dashed lines). Scores have been computed from 1 Nov to 4 Dec 2008. They are shown for two domains: (top) North America and (bottom) 20°–90°N area.

Fig. 5.

RMS scores of 48-h wind forecasts vs radiosonde observations at 250 hPa: operational scores (solid line); scores with a tuning of observation error variances and a global tuning of background error variances (dashed lines). Scores have been computed from 1 Nov to 4 Dec 2008. They are shown for two domains: (top) North America and (bottom) 20°–90°N area.

An observed side effect of the reduction of both observation and background errors is that more observations are rejected, since a first part of the ARPEGE quality control is done by comparing the squared innovations with the corresponding sums of the background and observation error variances. Such a larger rejection of erroneous observations can explain a part of the improvement on the forecast scores seen in Fig. 5.

Satellite observations are among those for which the error variance has been proportionally the most reduced. This seems to be related to large positive impact observed at high levels (Fig. 6).

Fig. 6.

As in Fig. 5, but for RMS scores of 96-h geopotential forecasts against ECMWF analyses at 20 hPa. Scores are shown for two domains: (top) 20°–90°N area and (bottom) 20°–90°S area.

Fig. 6.

As in Fig. 5, but for RMS scores of 96-h geopotential forecasts against ECMWF analyses at 20 hPa. Scores are shown for two domains: (top) 20°–90°N area and (bottom) 20°–90°S area.

8. DFS in the ARPEGE ensemble

Another application of the estimation of the Tr(ΠiTΠi𝗛𝗞) is the evaluation of the DFSi associated with the different subsets of observations and variables.

Figure 7 (top panel) shows the DFS by observation type. Not surprisingly, the largest DFS are observed for Radiosounding and Brightness temperature types, followed by the DFS of Aircraft observations. Cardinali et al. (2004) have obtained similar values for the DFS per observation type, but with a larger value for satellite observations. This is likely due to the fact that the European Centre for Medium-Range Weather Forecasts (ECMWF) analysis relies on a larger number of satellite data.

Fig. 7.

Number and partial DFS of observations by observation type: 1 = Synop, 2 = Aircraft, 3 = Satwind, 4 = Driftbuoy, 5 = Radiosounding, 6 = Profiler, 7 = Brightness temperature, 9 = Scatterometer, and 10 = GPSRO. (top) Globe, (middle) 20°–90°N area, and (bottom) 20°–90°S area.

Fig. 7.

Number and partial DFS of observations by observation type: 1 = Synop, 2 = Aircraft, 3 = Satwind, 4 = Driftbuoy, 5 = Radiosounding, 6 = Profiler, 7 = Brightness temperature, 9 = Scatterometer, and 10 = GPSRO. (top) Globe, (middle) 20°–90°N area, and (bottom) 20°–90°S area.

The comparison of the DFS with the corresponding numbers of observations shows that the DFS is not simply a question of number, since the weight of Radiosounding observations is proportionally larger than the weight of Brightness temperature observations. This can be interpreted as a lower confidence accorded to the last type of measurements, but it also reflects the fact that they are more integrated measures. Relative to their numbers, the DFS of Satwind and GPSRO measurements are high, which can be explained by the fact that they are given a large confidence and also because they are present in data-sparse areas like over the oceans.

In the Northern Hemisphere, the largest DFS are clearly for Radiosounding data, where they are quite numerous, followed by Aircraft and Brightness temperature data (Fig. 7, middle panel).

Despite the overall lower DFS per observation of Brightness temperature observations in the ARPEGE assimilation, they appear to be the main source of information in the Southern Hemisphere (Fig. 7, bottom panel), since there are few Radiosounding observations in this area of the globe.

Figure 8 (top panel) shows the DFS of the different observed variables given by all types of observations, except the Brightness temperature type. It can be seen that the largest DFS is for the wind components given by the upper-air observations (Aircraft, Satwind, Radiosounding, and Profiler) followed by the temperature measurements given by the same type of observations, except Profiler. For these variables, the DFS is nearly directly proportional to the number of measurements. The large DFS per observation of GPSRO observations again has to be noted. On the contrary, the DFS per observation is small for surface observations (Synop, Driftbuoy, and Scatterometer), except for the geopotential Z measurements given by Synop.

Fig. 8.

(left) Number and (right) partial DFS: (top) by variable, except brightness temperatures; (middle) by sounders; and (bottom) by channels (HIRS channels).

Fig. 8.

(left) Number and (right) partial DFS: (top) by variable, except brightness temperatures; (middle) by sounders; and (bottom) by channels (HIRS channels).

The DFS is a measure of the sum of the variations of the analysis brought by the variations of a given set of observations. Those variations of the analysis are clearly directly linked to the confidence given to the observations, through their observation error variances. The variations of the analysis are also constrained by the length scale of the background error correlations. If the length scale is relatively large in comparison with the density of observations, then the individual enabled variation of the analysis at an observation location will be small. On the contrary, if the length scale of the background error correlations is shorter, as for specific humidity, then the impact of observations will be larger. This can explain why the DFS per observation is high for specific humidity observations (Fig. 8, top panel). Of course, this does not say that humidity observations will have a main contribution on the reduction of forecast errors. The DFS has simply to be interpreted as a measure of the degrees of freedom for the analysis (the signal) to deviate from the background.

The DFS associated with different Brightness temperature sounders (Fig. 8, middle panel) are particularly interesting to examine, since they show that the DFS are not again simply functions of the numbers of observations. Note, in particular, the large values of the DFS per observation for High Resolution Infrared Radiation Sounder (HIRS), Advanced Microwave Sounding Unit-B (AMSU-B), Meteorological Satellite (Meteosat) clear-sky radiances (CSR), and, on the contrary, the low DFS/p for IASI observations.

Figure 8 (bottom panel) shows that the DFS per observation associated with HIRS channels 4–7 and 15 are quite low, when compared with those associated with channels 11–12, although they have comparable observation error standard deviation. A likely explanation of this behavior is that channels 4–7 and 15 are associated with much larger weighting functions than those of channels 11–12. As a consequence, the background error variances (not shown here), as diagnosed by a randomization procedure proposed by Andersson et al. (2000), show quite low values for channels 4–7 and 15 and larger values for channels 11–12. It means that the former channels will not produce a large deviation from the background. Another way to interpret this is that there is not a lot of information to get from a very integrated observation. A complementary explanation is that channels 11–12 are more sensitive to humidity than channels 4–7 and 15. Since background error correlations are shorter for humidity than for temperature, this will also allow channels 11–12 to produce larger local modifications to the background and then to be characterized by larger DFS.

9. Conclusions

Real-size assimilation systems rely on approximations that can make them suboptimal with respect to the combination of background and observation information. A posteriori diagnostics give a way to check the optimality of such systems and to improve them.

The use of an ensemble of assimilations is a way to optimize the description of background error covariance matrix 𝗕. Ensembles of assimilations can be based on an explicit perturbation of observations in order to provide implicit background perturbations, representative of background errors, from which useful statistics can be deduced. They are implemented, or planned to be implemented, in operational centers such as the ECMWF (Isaksen et al. 2007) or Météo-France (Berre et al. 2007).

It has been shown that such an ensemble can be used to evaluate the theoretical statistical estimations of subparts of the cost function measuring the distance between the analysis and the two sources of information, the background and the observations. The estimations appear to be robust in a real-size analysis system where the number of observations is large. Moreover, the robustness of the estimations can be increased by averaging them upon several dates.

The statistical estimations of subparts of the cost function can be used to tune the complete set of observation error standard deviations that appear to be overestimated in the ARPEGE 4D-Var assimilation. They can also allow the determination of a single normalization coefficient for the overall background error standard deviations; this error also happens to be (by chance) overestimated in the ARPEGE system.

The impact on the short-range forecasts of the reduction of both observation and background errors has been shown to be neutral to slightly positive at the low levels. The slight improvement on scores could be partly due to an improved rejection of erroneous observations against background, due to the better scaling of both observation and background error variances. A stronger positive impact has been observed in the stratosphere, where satellite radiance errors were overestimated to a larger extent than other error variances.

Since the ensemble of assimilations is based on perturbations of observations consistent with observation error variances, the right amplitude of those variances is also expected to have an influence on the correct dispersion of the ensemble. Although the ARPEGE analysis is based on the 4D-Var formalism, the ensemble of assimilations is run, for the time being, in a 3D-Var FGAT mode at Météo-France. 4D-Var ensembles have been occasionally run and actually show that the DFS given by the two types of ensembles are quite close. However, it is planned to more carefully investigate the differences between the background errors and DFS produced by 3D-Var FGAT and 4D-Var ensembles.

Desroziers et al. (2005a) have proposed other diagnostics, in observation space, to diagnose observation, background, and analysis errors. A first comparison, not shown here, between this method and the procedure described in the present paper shows a nice agreement between the two approaches. It is planned to systematically compare the error variances given by the two methods and to try to establish why the two methods compare so well although they rely on apparently different theoretical rationale.

Up to now, only a single normalization coefficient of matrix 𝗕 has been determined. There is a scope for investigating whether the a posteriori diagnostics described in this paper but also the observation space diagnostics proposed in Desroziers et al. (2005a) can be used to more finely optimize matrix 𝗕.

As shown in the paper, a direct by-product of an ensemble of perturbed analyses is the evaluation of the weights of observations through their DFS. Since the associated computations are costless, they can be implemented easily to provide an additional routine monitoring tool in a real-size system.

Acknowledgments

We thank the two referees for their constructive comments on the first version of the paper.

REFERENCES

REFERENCES
Andersson
,
E.
,
M.
Fisher
,
R.
Munro
, and
A.
McNally
,
2000
:
Diagnosis of background errors for radiances and other observable quantities in a variational data assimilation scheme, and the explanation of a case of poor convergence.
Quart. J. Roy. Meteor. Soc.
,
126
,
1455
1472
.
Belo Pereira
,
M.
, and
L.
Berre
,
2006
:
The use of an ensemble approach to study the background error covariances in a global NWP model.
Mon. Wea. Rev.
,
134
,
2466
2489
.
Bennett
,
A.
,
L.
Leslie
,
C.
Hagelberg
, and
P.
Powers
,
1993
:
Tropical cyclone prediction using a barotropic model initialized by a generalized inverse method.
Mon. Wea. Rev.
,
121
,
1714
1729
.
Berre
,
L.
,
S.
Stefanescu
, and
M.
Belo Pereira
,
2006
:
The representation of analysis effect in three error simulation techniques.
Tellus
,
58A
,
196
209
.
Berre
,
L.
,
O.
Pannekoucke
,
G.
Desroziers
,
S. E.
Stefanescu
,
B.
Chapnik
, and
L.
Raynaud
,
2007
:
A variational assimilation ensemble and the spatial filtering of its error covariances: Increase of sample size by local spatial averaging.
Proc. ECMWF Workshop on Recent Developments in Data Assimilation for Atmosphere and Ocean, Reading, United Kingdom, ECMWF, 151–168
.
Bormann
,
N.
,
S.
Saarinen
,
G.
Kelly
, and
J-N.
Thépaut
,
2003
:
The spatial structure of observation errors in atmospheric motion vector from geostationary satellite data.
Mon. Wea. Rev.
,
131
,
706
718
.
Buehner
,
M.
,
P.
Gauthier
, and
Z.
Liu
,
2005
:
Evaluation of new estimates of background- and observation-error covariances for variational assimilation.
Quart. J. Roy. Meteor. Soc.
,
131
,
3373
3383
.
Cardinali
,
C.
,
S.
Pezzuli
, and
E.
Andersson
,
2004
:
Influence matrix diagnostic of a data assimilation system.
Quart. J. Roy. Meteor. Soc.
,
130
,
2767
2786
.
Chapnik
,
B.
,
2005
:
Réglage des statistiques d’erreur en assimilation variationnelle (Tuning of error statistics in variational assimilation).
Thèse de doctorat de l’Université Toulouse III, 196 pp
.
Chapnik
,
B.
,
G.
Desroziers
,
F.
Rabier
, and
O.
Talagrand
,
2004
:
Properties and first application of an error statistics tuning method in variational assimilation.
Quart. J. Roy. Meteor. Soc.
,
130
,
2253
2275
.
Chapnik
,
B.
,
G.
Desroziers
,
F.
Rabier
, and
O.
Talagrand
,
2006
:
Diagnosis and tuning of observational error statistics in a quasi operational data assimilation setting.
Quart. J. Roy. Meteor. Soc.
,
132
,
543
565
.
Courtier
,
P.
,
J-N.
Thépaut
, and
A.
Hollingsworth
,
1994
:
A strategy for operational implementation of 4D-Var, using an incremental approach.
Quart. J. Roy. Meteor. Soc.
,
120
,
1367
1387
.
Daescu
,
D.
,
2008
:
On the sensitivity equations of four-dimensional variational (4DVAR) data assimilation.
Mon. Wea. Rev.
,
136
,
3050
3065
.
Dee
,
D. P.
,
1995
:
On-line estimation of error covariance parameters for atmospheric data assimilation.
Mon. Wea. Rev.
,
123
,
1128
1145
.
Dee
,
D. P.
, and
A.
da Silva
,
1999
:
Maximum-likelihood estimation of forecast and observation error covariance parameters. Part I: Methodology.
Mon. Wea. Rev.
,
127
,
1822
1834
.
Derber
,
J. C.
, and
F.
Bouttier
,
1999
:
A reformulation of the background error covariances in the ECMWF global data assimilation system.
Tellus
,
51A
,
195
221
.
Desroziers
,
G.
, and
S.
Ivanov
,
2001
:
Diagnosis and adaptive tuning of observation-error parameters in a variational assimilation.
Quart. J. Roy. Meteor. Soc.
,
127
,
1433
1452
.
Desroziers
,
G.
,
L.
Berre
,
B.
Chapnik
, and
P.
Poli
,
2005a
:
Diagnosis of observation, background and analysis-error statistics in observation space.
Quart. J. Roy. Meteor. Soc.
,
131
,
3385
3396
.
Desroziers
,
G.
,
P.
Brousseau
, and
B.
Chapnik
,
2005b
:
Use of randomization to diagnose the impact of observations on analyses and forecasts.
Quart. J. Roy. Meteor. Soc.
,
131
,
2821
2837
.
Evensen
,
G.
,
1994
:
Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics.
J. Geophys. Res.
,
99
, (
C5
).
3104
3127
.
Fisher
,
M.
,
2003a
:
Background error covariance modelling.
Proc. ECMWF Workshop on Recent Developments in Data Assimilation for Atmosphere and Ocean, Reading, United Kingdom, ECMWF, 97–123
.
Fisher
,
M.
,
2003b
:
Estimation of entropy reduction and degrees of freedom for signal for large variational analysis systems.
Tech. Rep. 397, ECMWF, 20 pp
.
Fisher
,
M.
, and
E.
Andersson
,
2001
:
Developments in 4D-Var and Kalman filtering.
Tech. Rep. 347, ECMWF, 38 pp
.
Girard
,
D.
,
1987
:
A fast Monte Carlo cross-validation procedure for large least squares problems with noisy data.
Tech. Rep. 687-M, IMAG, Grenoble, France, 22 pp
.
Hollingsworth
,
A.
, and
P.
Lönnberg
,
1986
:
The statistical structure of short-range forecast errors as determined from radiosonde data. Part I: The wind field.
Tellus
,
38A
,
111
136
.
Houtekamer
,
P.
,
L.
Lefaivre
,
J.
Derome
,
H.
Ritchie
, and
H.
Mitchell
,
1996
:
A system simulation approach to ensemble prediction.
Mon. Wea. Rev.
,
124
,
1225
1242
.
Hutchinson
,
M.
,
1989
:
A stochastic estimator for the trace of the influence matrix for Laplacian smoothing splines.
Commun. Stat. Simul.
,
18
,
1059
1076
.
Isaksen
,
L.
,
M.
Fisher
, and
J.
Berner
,
2007
:
Use of analysis ensembles in estimating flow-dependent background error variance.
Proc. ECMWF Workshop on Flow-Dependent Aspects of Data Assimilation, Reading, United Kingdom, ECMWF, 65–86
.
Langland
,
R.
, and
N.
Baker
,
2004
:
Estimation of observation impact using the NRL atmospheric variational data assimilation adjoint system.
Tellus
,
56A
,
189
201
.
Lewis
,
J. M.
, and
J. C.
Derber
,
1985
:
The use of adjoint equations to solve a variational adjustment problem with advective constraints.
Tellus
,
37A
,
309
322
.
Munro
,
R.
,
C.
Köpken
,
G.
Kelly
,
J-N.
Thépaut
, and
R.
Saunders
,
2004
:
Assimilation of Meteosat radiance data within the 4D-Var system at ECMWF: Data quality monitoring, bias correction and single-cycle experiments.
Quart. J. Roy. Meteor. Soc.
,
130
,
2293
2313
.
Pannekoucke
,
O.
,
L.
Berre
, and
G.
Desroziers
,
2007
:
Filtering properties of wavelets for the local background error correlations.
Quart. J. Roy. Meteor. Soc.
,
133
,
363
379
.
Rabier
,
F.
,
H.
Järvinen
,
E.
Klinker
,
J-F.
Mahfouf
, and
A.
Simmons
,
2000
:
The ECMWF operational implementation of four-dimensional variational assimilation. I: Experimental results with simplified physics.
Quart. J. Roy. Meteor. Soc.
,
126
,
1143
1170
.
Rabier
,
F.
,
N.
Fourri
,
D.
Chafa
, and
P.
Prunet
,
2002
:
Channel selection method for infrared atmospheric sounding interferometer radiances.
Quart. J. Roy. Meteor. Soc.
,
128
,
1011
1032
.
Raynaud
,
L.
,
L.
Berre
, and
G.
Desroziers
,
2008
:
Spatial averaging of ensemble-based background-error variances.
Quart. J. Roy. Meteor. Soc.
,
134
,
1003
1014
.
Rodgers
,
C.
,
2000
:
Inverse Methods for Atmospheric Sounding Theory and Practice.
World Scientific Publishing, 256 pp
.
Talagrand
,
O.
,
1997
:
Assimilation of observations, an introduction.
J. Meteor. Soc. Japan
,
75
,
191
209
.
Talagrand
,
O.
,
1999
:
A posteriori verification of analysis and assimilation algorithms.
Proc. ECMWF Workshop on Diagnosis of Data Assimilation Systems, Reading, United Kingdom, ECMWF, 17–28
.
Talagrand
,
O.
, and
P.
Courtier
,
1987
:
Variational assimilation of meteorological observations with the adjoint vorticity equation. I: Theory.
Quart. J. Roy. Meteor. Soc.
,
113
,
1311
1328
.
Trémolet
,
Y.
,
2008
:
Computation of observation sensitivity and observation impact in incremental variational data assimilation.
Tellus
,
60A
,
964
978
.
Wahba
,
G.
,
D.
Johnson
,
F.
Gao
, and
J.
Gong
,
1995
:
Adaptative tuning of numerical weather prediction models: Randomized GCV in three- and four-dimensional data assimilation.
Mon. Wea. Rev.
,
123
,
3358
3369
.
Zagar
,
N.
,
E.
Andersson
, and
M.
Fisher
,
2005
:
Balanced tropical data assimilation based on a study of equatorial waves in ECMWF short-range forecast errors.
Quart. J. Roy. Meteor. Soc.
,
131
,
987
1011
.
Zhu
,
Y.
, and
R.
Gelaro
,
2008
:
Observation sensitivity calculations using the adjoint of the gridpoint statistical interpolation (GSI) analysis system.
Mon. Wea. Rev.
,
136
,
335
351
.

APPENDIX A

Proof of Eq. (5)

The incremental implementation of 4D-Var, with the linearization of observation operators around successive updated nonlinear trajectories, attempts to solve the original problem of minimizing the following cost function, with nonlinear observation operators, including model integration:

 
formula

The analysis xa is the vector that makes the gradient of J vanish. That can be written as

 
formula
 
formula
 
formula
 
formula

where 𝗛 is the linearized observation operator around xa.

It follows that

 
formula
 
formula

with

 
formula
 
formula

APPENDIX B

Proof of Eq. (14)

Chapnik (2005) has given a proof of Eq. (14) with a linear observation operator. This proof is generalized here for the nonlinear case.

The difference between the analysis xa and the extended observation vector z (merging the background xb and the proper observation vector yo) can be written as

 
formula
 
formula
 
formula

where Γ is the linearized generalized observation operator around xa.

Following Eq. (5), the analysis error εa can be written as

 
formula
 
formula
 
formula

where Γ+ is the left inverse of Γ such as Γ+Γ = 𝗜.

Using Eqs. (B3) and (B6) and the linearity property of the expectation operator, it follows that

 
formula
 
formula
 
formula
 
formula

It is easy to check that (ΓΓ + − 𝗜)𝗦(ΓΓ + −𝗜)T = 𝗦 − Γ𝗣aΓT and thus

 
formula
 
formula

APPENDIX C

Proof of Eqs. (23) and (24)

Using Eq. (14), the expectation of a subpart of Jo(xa) can be written as

 
formula
 
formula

If the observation errors associated with subset i are independent from the other observation errors, then and 𝗥−1 commutes with ΠiTΠi so that

 
formula
 
formula
 
formula

On another hand, the trace of a matrix 𝗔 can be approximated with

 
formula
 
formula

where the ηl are a set of L vectors of random numbers with standard Gaussian distribution (mean 0 and variance 1; Girard 1987).

Equation (C5) can be rewritten as

 
formula
 
formula

Applying the randomized estimation of the trace of a matrix to this last equation gives

 
formula
 
formula
 
formula

where δlo = 𝗥1/2ηl is a vector of perturbations on the complete set of observations and is the restriction of δlo to the subset of observations i. In an ensemble of assimilations, the perturbations on the background, the observations, and the analysis are linked by the relation δa = (𝗜 − 𝗞𝗛)δb + 𝗞δo. Furthermore, the perturbations δo and δb are uncorrelated, so that , which proves Eq. (23).

The proof of Eq. (24) can be obtained with the same kind of derivations.

Footnotes

Corresponding author address: Gérald Desroziers, GAME/CNRM, Météo-France, CNRS, 42, av. G. Coriolis, Toulouse, 31057 CEDEX, France. Email: gerald.desroziers@meteo.fr

This article included in the Mathematical Advances in Data Assimilation (MADA) special collection.