Representing model error in ensemble data assimilation

The paper investigates a method to represent model error in the ensemble data assimilation (EDA) system. The ECMWF operational EDA simulates the effect of both observations and model uncertainties. Observation errors are represented by perturbations with statistics characterized by the observation error covariance matrix whilst the model uncertainties are represented by stochastic perturbations added to the physical tendencies to simulate the effect of random errors in the physical parameterizations (STmethod). In this work an alternative method (XB-method) is proposed to simulate model uncertainties by adding perturbations to the model background field. In this way the error represented is not just restricted to model error in the usual sense but potentially extends to any form of background error. The perturbations have the same correlation as the background error covariance matrix and their magnitude is computed from comparing the high-resolution operational innovation variances with the ensemble variances when the ensemble is obtained by perturbing only the observations (OBS-method). The XB-method has been designed to represent the short range model error relevant for the data assimilation window. Spread diagnostic shows that the XBmethod generates a larger spread than the ST-method that is operationally used at ECMWF, in particular in the extratropics. Three-dimensional normal-mode diagnostics indicate thatXB-EDAspread projects more than the spread from the other EDAs onto the easterly inertia-gravity modes associated with equatorial Kelvin waves, tropical dynamics and, in general, model error sources. The background error statistics from the above described EDAs have been employed in the assimilation system. The assimilation system performance showed that the XBmethod background error statistics increase the observation influence in the analysis process. The other EDA background error statistics, when inflated by a global factor, generate analyses with 30–50 % smaller degree of freedom of signal. XB-EDAbackground error variances have not been inflated. The presented EDAs have been used to generate the initial perturbations of the ECMWF ensemble prediction system (EPS) of which theXB-EDA induces the largest EPS spread, also in the medium range, leading to a more reliable ensemble. Compared to ST-EDA, XB-EDA leads to a small improvement of the EPS ignorance skill score at day 3 and 7.


Introduction
Data assimilation systems combine observations and background state, usually a short-range forecast of 6 or 12 h.Alternative approaches to a deterministic initial condition system based on ensemble methodologies such as the ensemble Kalman filter or the ensemble variational analysis are nowadays widely used in numerical weather prediction (NWP).They have the advantage of providing information on flowdependent background error covariances.Nevertheless, ensemble data assimilation (EDA) systems are usually implemented in the context of a perfect model and the failure of representing model error affects for example the computation of the background error variances, which tends to be underestimated.Different strategies have been tested to take model error into account by rescaling the ensemble of analyses with constant, isotropic factors (inflation), the multi-model and multi-physical parameterization approach (Houtekamer et al., 1996(Houtekamer et al., , 2005)), the stochastic physical tendency (Buizza et al., 1999) and the backscatter stochastic kinetic energy schemes (Shutts, 2005;Berner et al., 2008).Comparison of different model error representations has been performed by Houtekamer et al. (2009), which indicates that the inflation approach has the largest contribution in their model error simulation.
Several operational NWP centres such as Météo-France, UK Met Office and Environment Canada have implemented an EDA system.Recently, Météo-France has implemented a six-member EDA where not only the observation uncertainties are represented but also the model uncertainties are explicitly taken into account.This is done by inflating the background field with a latitude-level varying factor in the ensemble (see Raynaud et al., 2009Raynaud et al., , 2012)).Their first operational EDA configuration (Berre et al., 2007) has been used in operations since July 2008.Flow-dependent background error variances for the operational 4D-Var assimilation system (Raynaud et al., 2011) were derived within a perfect model framework and the estimated variances were inflated "off-line" (i.e. after the ensemble has been completed) by using a posteriori diagnostics (Desroziers and Ivanov, 2001).
The inflation aims at representing model error contributions.Multiplicative inflation is in fact a simple and widely used technique to deal with unknown error sources.Because the off-line multiplicative inflation that was applied to the variances was not accounted for in the background perturbation update, the recently implemented operational EDA configuration includes the multiplicative inflation to enlarge the amplitude of forecast perturbations within the ensemble.
The Environment Canada ensemble system is based on a Kalman filter (Houtekamer and Mitchell, 2005;Houtekamer et al., 2005) and has been used operationally since January 2005.It provides an ensemble of initial conditions for the medium-range EPS (ensemble prediction system) and represents both observations and model sources of uncertainties.In particular the model error component has been extensively investigated.Initially a simplified and reduced amplitude form of the Canadian 3D-Var background error covariance was used to perturb the ensemble of background fields (Houtekamer et al., 2005).Then, different ways to determine the covariance for the additive model error component have been investigated (Hamill and Whitaker, 2005) and model error perturbations are added to the ensemble analysis rather than to the background ensemble (Houtekamer and Mitchell, 2005) to account for the data assimilation weakness.More recently, each member has a different model version (Meng and Zhang, 2007;Fujita et al., 2007) to represent uncertainties in model representation of physical processes.As already said above, Houtekamer and Mitchell (2005) concluded that the addition of isotropic model error perturbations to the ensemble of analyses is found to have the largest impact in terms of ensemble spread.
An EDA was operationally implemented at ECMWF in June 2010 (Isaksen et al., 2010).The EDA ensemble consists of ten independent members of lower resolution (with respect to the high-resolution operational 4D-Var system) 4D-Var data assimilation systems with perturbed observations and perturbed model tendencies.In particular, the observation uncertainties are represented by perturbing the observations and the model uncertainties by adding stochastic perturbations to the model tendencies during the first 12 h model evolution using the Stochastically Perturbed Parameterization Tendency scheme (SPPT, see Palmer et al., 2009, for a review).
The ECMWF EDA provides a flow-dependent or daily model background error covariance matrix that is supposed to improve the high resolution analysis system by better representing the daily dynamical synoptic features (Raynaud et al., 2008(Raynaud et al., , 2009(Raynaud et al., , 2011;;Buehner et al., 2010).Since its implementation, the EDA has been used together with the singular vectors to initialize the operational EPS (Molteni et al., 1996;Buizza et al., 2007) and to improve the simulation of initial uncertainties (Buizza et al., 2008), one of the fundamental aspects of the EPS design.
In this paper, a different way of representing model error in the operational ECMWF EDA is presented and compared to the standard SPPT method.The model uncertainties are represented by adding perturbations to the model background field.The magnitude of the perturbations varies with vertical level and with geographical latitude.They are estimated from a comparison between the innovation variance of the high resolution 4D-Var system, i.e. the difference between observation and background at the observation location, and the ensemble data assimilation variance (variance taken over an ensemble of assimilations over a 3-week period) in which only observation uncertainties are represented.The model error representation is therefore similar to the one introduced by Raynaud et al. (2012) at Météo-France that is referred to as the multiplicative perturbation method.The method described here is denoted as an additive perturbation method.The error represented is not restricted to the model error in the usual sense, i.e. the error that would be present in the forecast even if the initial condition were exact, but is related to any form of error, for example errors in the background covariance matrix coming from the operational ECMWF EDA.The present paper studies the EDA sensitivity to the different model error representations.The proposed method is compared to the operational one, which uses the stochastically perturbed parameterization tendency scheme to simulate model uncertainties, and to the EDA obtained by representing only observation uncertainties.A fourth EDA has also been designed to just quantify the impact of background cycling in the EDA where only observations are perturbed.Observation uncertainties are always equally represented in all EDAs examined.
Section 2 describes the methodology used to simulate model uncertainties.Section 3 analyses the spread characteristics of the investigated EDAs by using a variety of diagnostics and with an evaluation of the background error covariance matrix provided by the EDAs.The EPS performance sensitivity to the EDAs is also discussed.Conclusions are drawn in Sect. 4.

The XB-EDA
An ensemble of analyses attempts to generate a representative sample of possible states of a dynamical system.The samples are generated by the same assimilation system.From the optimal solution of the analysis problem, x a = f (x b , y), two input parameters can be identified: the observation vector y and the background vector x b obtained from a shortrange forecast, respectively.An ensemble of analyses can be generated by perturbing both input vectors.In particular, the observations uncertainties can be represented by perturbing vector y, whilst the model uncertainties (at least the short-range model error) can be represented by perturbing the model state vector x b .The perturbed analysis equation can be written as where ζ and η are perturbations defined as with f (λ, l, x) being a function of latitude (λ), model level (l) and model parameters (x).ζ and η are samples of vectors drawn from a multi-dimensional Gaussian distribution with zero mean and identity covariance matrix.To achieve that the final perturbations ζ and η have a covariance matrix specified by B and R, respectively, the square root of B and R is applied to the sequence of normally distributed vectors ζ and η.B and R are the estimated background and observation error covariance matrices; they are therefore only approximations of the true covariance matrices.When R is diagonal (i.e.uncorrelated observation errors) a simple multiplication by the observation error standard deviation σ o is applied (Eq. 1b).Only two sets of observations are perturbed with spatially correlated patterns.One is the atmospheric motion vector (AMV) observation (Bormann et al., 2003) and the other set is the sea-surface temperature field (Vialard et al., 2005).The magnitude of the final perturbation ζ is determined by f (λ, l, x) which is estimated by comparing the variance of the innovation vector d (over 3 weeks) with the ensemble data assimilation variance, Var(EDA) (estimated over an ensemble of assimilations over 3-week period), in the case the ensemble data assimilation is obtained by only perturbing the observations.The coefficient f (λ, l, x) is meant therefore to compensate for the discrepancy between the background error as obtained from the innovation d on the one hand, and the a priori background error covariance matrix B on the other.The innovation vector is the difference between the observation vector y and the background counterpart of the observation computed by using the nonlinear observation operator (H (x b )).Under the assumption of unbiasedness of the errors and de-correlation between the background and observation errors, the background error variance, as obtained from the innovation, is Var(d) − σ 2 o .The scalar function f (λ, l, x) is hence defined as where σ 2 o is the prescribed observation error variance.If Var(d) − σ 2 o is less or equal to the EDA variance it is imposed that f (λ, l, x) = 0.The perturbation amplitude modulation hence varies in the interval [0,1].The innovation variances have been computed for 10 hPa pressure layers for atmospheric measurements located between the surface and 50 hPa (wind observations), between surface and 5 hPa (temperature observations), and surface and 300 hPa when humidity observations are considered.Three latitude bands, namely Northern Hemisphere (20 • N, 90 • N), Southern Hemisphere (20 • S, 90 • S) and tropics (20 • S, 20 • N), and a 3-week data set have been considered.For the u and v component of the wind all conventional observations (radiosondes, pilots, synops, aircrafts and profilers), AMVs and scatterometer observations have been used to compute the innovation variance.For temperature conventional (radiosondes and aircrafts) and AMSU-A observations, for humidity radiosondes and All-Sky (SSM/I and TMI radiances) observations and, finally, for surface pressure all land and ocean stations have been used.The variation of f with latitude band, model level and for model parameters u, T , and q, is shown in Fig. 1. Figure 1a shows that for the u component f decreases in the troposphere, when pressure increases, down to zero in the tropics and down to 0.3 in the extratropics (level 1 at 0.01 hPa, identifies the top of the atmosphere).If observations are unavailable to estimate the innovation variance, the modulation factor f is kept constant, i.e. from model level 1-30 for wind, from model level 1-18 for temperature and model level 1-55 for humidity.Similar results for the modulation factor are obtained for the v component of the wind (not shown).In the lower troposphere close to the surface the modulation factor globally increases on average up to 0.6.For temperature (Fig. 1b) its magnitude increases with the increase of pressure on average for the three latitude bands from 0.3 to 0.5, the tropical modulation factor always being the smallest.
For humidity (Fig. 1c) f rapidly grows with the atmospheric pressure level up to 0.8 (Southern Hemisphere) towards the surface.Concerning surface pressure, the correction (1c) (not shown) is globally constant and around 0.4.During the cycling, f has been re-computed for retuning purposes using Eq.(1c) every 3 days and by using the past 3-day variance sample.However, the modulation factor has stabilized rather fast after 2 days of cycling.
The background is perturbed at the start of each assimilation window and the innovation vector computed along the trajectory starting from the perturbed background to correctly take into account the background changes (H is the non-linear model and observation operator) and to produce a balanced perturbed field.Figure 2 schematically represents the realization of the described XB-EDA ensemble.From a control (unperturbed) analysis the two set of perturbations, η for the observations and ζ for the background, are respectively added at the beginning of the 4D-Var assimilation window to create 10 different initial conditions (members).The two sets of perturbations are recomputed and added to the observations and the background field at every analysis cycle.
The ζ perturbation accounts for the short-range model error sources including the fraction of the analysis error that is due to the model error.Sources of error are therefore not only related to physical parameterizations but also to the dynamics, the spatial and temporal discretization, the linearized physical process and the misspecification of the probability distribution of errors in the observations and the background model.

The other EDA
In all EDAs presented here the observation uncertainties are represented by perturbing the observations whilst the model uncertainties are produced differently.Table 1 shows all the EDA configurations investigated.In particular OBS-EDA is the ensemble analysis where only observation uncertainties are represented (different observation values generate different background fields during the cycling process).The OBS-OBS ensemble is similar to the OBS-EDA one but the background fields are not cycled.At every cycle the 10-member background fields are replaced with the background field of the control, i.e. unperturbed analysis (Fig. 2).OBS-OBS is performed to evaluate the impact of the cycling by the background.The model error representation in ST-EDA is based on the SPPT scheme (Buizza et al., 1999;Palmer et al., 2009).Since November 2009, the operational EPS uses SPPT and the stochastic backscatter schemes (SPBS, Shutts, 2005;Berner et al., 2008) to simulate model uncertainties.The SPBS scheme simulates the inverse energy cascade due to the interaction between the unresolved and the resolved scales, and aims to compensate for the over-dissipation occurring in numerical models.The SPPT scheme, designed to simulate random model errors due to physical parameterizations, is still assumed to explain the largest source of model error in the EPS.The operational EDA configuration does not use the SPBS scheme but only the SPPT.In particular, SPPT model uncertainty is simulated using the 1-scale version of the stochastically perturbed parameterization tendency scheme, which perturbs the total parameterized tendency of physical processes.Since the 1-scale version use a timescale of 6 h, there is no need to cycle the model error perturbation across different data assimilation cycles.In the 1-scale version of the SPPT, the perturbations to the physical tendencies are defined to have a spatial correlation length of 500 km and a time correlation of 6 h, as in the original SPPT scheme (Buizza et al., 1999).
Because of the presumed ensemble analysis underdispersivity (to be confirmed later) a global inflation factor (= 1.4) has been applied to the background error standard deviation in the OBS, OBS-OBS and ST-EDAs methods to increase the ensemble spread and to penalize the model background further with respect to the observations in the assimilation process.The static background covariance matrix has, in fact, always been inflated in the ECMWF 4DVar assimilation system to avoid the excessive weight given to the background with respect to the observations.Indeed, studies on the observational influence in the analysis system have shown that globally and for a given assimilation cycle only 15 % of the information was provided by the observations while the remaining 85 % were due to the background (Cardinali et al., 2004;Cardinali, 2013).Unfortunately, the inflation is a constant that does not vary with respect to the parameters, with respect to the geographical location or weather situation and the resulting ensemble spread is simply globally amplified.

Results
In this section EDAs with different model error representations are compared and diagnosed.Each EDA includes 10 perturbed and 1 unperturbed 12 h 4D-Var assimilations (Rabier et al., 2000;Janisková et al., 2002;Tompkins and Janisková, 2004;Lopez and Moreau, 2005) at the resolution of T L 399L91 (spectral triangular truncation with 399 wave numbers and linear grid, and on 91 vertical levels) for the model forecast and T L 159L91 for the minimization calculation, respectively.

EDA spread
In Fig. 3 the averaged spread of the four data assimilation ensembles is compared for the zonal wind component.The average spread has been computed over the period 5 October-15 November 2008 (the first 5 days of the EDA computation have not been included in the evaluation to take into account "spin-up") from 6 h forecasts according to the expression where m i is the ith ensemble member, N = 10 and m is the ensemble mean.Expectation stands for averaging over longitude and over the selected period.
When only observation uncertainties are represented (u wind component OBS-EDA, Fig. 3a), the spread is mainly confined to the upper stratosphere (above model level 20, i.e. ∼ 10 hPa) and to the troposphere (below model level 40, i.e. ∼ 110 hPa) in the tropics.Very little spread is accomplished poleward of 40 • N or 40 • S. To understand how much of the spread is due to the cycling of OBS-EDA over successive assimilation windows, Fig. 3b shows the OBS-OBS spread.The OBS-OBS spread is mainly confined to the stratosphere and with a smaller amount to the tropical troposphere.When model errors are explicitly represented in the ensemble the spread increases according to the method applied.Compared to OBS-EDA, ST-EDA presents 10 % larger spread in the tropics and the mid-latitudes (Fig. 3c).The XB-EDA ensemble (Fig. 3d) shows the largest and more globally distributed spread.It amounts to ∼ 1 and 4 m s −1 larger than ST-EDA in the tropics and the stratosphere, respectively.It is also remarkably different in the high and medium extratropical troposphere and differences exceed 4-5 m s −1 in the midlatitude.Blank areas are values between 0 and 0.5.Figure 4a shows the difference of spreads in the OBS and OBS-OBS EDAs and illustrates the impact of cycling over successive assimilation windows.Figure 4b shows the difference of spreads in XB and OBS EDAs and illustrates the impact of the perturbation of the background (first part of Eq. 1b).Most of the former difference is located in the stratosphere and in the tropics while the second extends to a large part of the troposphere, especially in the Southern Hemisphere.The figures also suggest that the differences are everywhere positive (blank areas are values between 0 and 0.5).
Figure 5 shows the reduction of globally averaged spread of the OBS-, OBS-OBS-and ST-EDAs relative to the XB-EDA for each model level and the u component of the wind.The spread in XB-EDA is a function of the amplitude of the perturbation applied to the background field.The spread loss with respect to the XB ensemble is decreasing with the increase of model level for ST-and OBS-EDA whereas it is constant for the OBS-OBS EDAs.Close to the surface (i.e.model level 91) the first two EDAs lose 20 and 40 % spread and close to the top of the atmosphere 40-50 % on average at all latitudes, respectively.
In all EDAs the largest spread is located in the stratosphere where also the largest loss with respect to XB is observed.Similar results are obtained for the temperature field but the magnitude of the spread loss is 25 % smaller (not shown).

Spread case study
An example of spread differences among the OBS-, ST-and XB-EDAs in physical space is presented for 20-21 October 2008 in Fig. 6.In this period an intense baroclinic development took place over the northern West Pacific.In addition to a mature-stage cyclone moving toward the Gulf of Alaska, a deep cyclogenesis took place about 4000 km westward at about 50 • N, 180 • W (not shown).All methods produce spread associated with the two cyclones but differences exist in the structure and magnitude of the spread.The comparison of 6 h forecast vorticity spread at 850 hPa valid on 21 October at 12:00 UTC (Fig. 6) shows that the OBS-EDA spread associated with the western cyclone is not only smaller than in the other two EDAs but also located on the north and northeastern side of the system (Fig. 6a).The STand XB-EDAs both contain spread over a larger area.For the mature-stage cyclone in the eastern Northern Pacific, the spread in all three experiments has the typical comma shape (East Pacific, around 45 • N-150 • W) of frontal systems associated with baroclinic development.The ST-and XB-EDAs are the most similar although XB-EDA (Fig. 6c) has a larger maximal amplitude and occupies a larger area.It is worth noticing that XB-EDA is the only ensemble showing some degree of uncertainty in the polar region where only very few observations are available and consequently the analysis uncertainties should be larger.

Modal diagnosis of the ensemble spread
The discussion of the results in the previous sections is complemented hereafter with the diagnosis of the ensemble spread in terms of normal modes as described in Žagar et al. (2011).The 10 members of the four ensemble experiments are projected onto a set of three-dimensionally orthogonal vectors which are eigensolutions of the Navier-Stokes equations linearized about an horizontally homogeneous stable state of rest.This projection allows the attribution of the ensemble spread according to various horizontal and vertical scales as well as linearly balanced (quasi-geostrophic) and unbalanced (inertio-gravity, IG) parts of the flow.In Žagar et al. (2011), the method was applied to the analysis of the ensemble spread of the DART-CAM system whereas Žagar et al. (2013) used the method to study the balance properties of the ECMWF EDA system.In the ECMWF EDA for July 2007, it was found that about 50 % of the short-range forecast-error variance was associated with the IG modes and that the eastward-propagating IG component was dominant on all scales.Both results were associated with the majority of EDA variance being present in the tropics.On the other hand, the ensemble spread of the DART-CAM ensemble was characterized by a prevalence of the westward-moving IG modes which was found to be related to the covariance inflation.These studies suggest that the normal mode function (NMF) expansion is a useful diagnostic of EDA systems.
In the present study, we followed Žagar et al. (2013) to analyse model levels under 10 hPa (model levels 19-91, totalling 73 levels) in order to avoid very large spread in the mesosphere (see Fig. 3) that projects strongly on the leading vertical modes and can obscure the interpretation of the results.For the presented diagnostics, 6 h forecast starting at 18:00 UTC in the period 18 October-16 November 2008 are used (30 samples) for 10 ensemble members.The analysed data on the N64 Gaussian grid are projected onto 85 zonal wave numbers, 50 vertical modes and 40 meridional modes for each motion type, namely balanced, eastward inertio-gravity (EIG mode) and westward inertiogravity (WIG mode).
In modal space, the ensemble spread based on N ensemble members is defined as where χ ν,i is a non-dimensional complex projection coefficient for an ensemble member i while ν is a four-indices modal index which contains information about the zonal wave number, the meridional mode, the vertical mode and the wave type.The overbar stands for averaging over N ensemble members, i.e. χ υ is the ensemble mean, Each vertical mode is characterized by a value of "the equivalent depth" H which couples horizontal and vertical motions.The spread computed by Eq. (2b) applies at a single 6 h forecast range and the results are presented as time averages over 30 samples.For details of the projection procedure see Žagar et al. (2011) and references therein.One difference between the normal-mode diagnostics and other spread evaluation methods consists in analysing simultaneously the mass and wind fields making possible the physical interpretation of balance relations.Spectra of the balanced, EIG and WIG energy spread as a function of the zonal wave number are shown in Fig. 7.In agreement with what has been presented so far, the ensemble spread in the XB-EDA dominates over the ST-EDA and OBS-EDA at all scales and for all three motion types.The ST-EDA spread is closer to the OBS-EDA spread than to the XB-EDA spread.In all experiments, the EIG spread dominates over WIG at largest horizontal scales and in XB-EDA it is greater than the WIG spread on all scales due to the equatorial Kelvin modes (not shown).The XB-, ST-and OBS-EDAs have a smaller percentage of their spread in the largest scale with respect to OBS-OBS-EDA.Below zonal wave number 10, there is around 30 % of the total spread for all experiments, varying from 28 % for ST-EDA to 34 % for OBS-OBS-EDA (not shown).
On average, the total spread in XB-EDA is between 1.7 and 1.6 times greater than the OBS-EDA spread.The ST-EDA spread is around 1.25 times of the OBS-EDA spread.Both XB-and ST-EDA add relatively more spread in the IG part than in the balanced part with respect to OBS-EDA.Instead, the experiment without cycling (OBS-OBS-EDA) counts only for 40 % of the spread of OBS-EDA for all modes (not shown).When the spread is summed up across all scales, the percentages of ROT, EIG and WIG spread in the four experiments are 43, 28, 29 % for OBS-OBS-EDA, 41, 29, 30 % for OBS-EDA, 39, 30, 31 % for ST-EDA, and 40, 31, 29 % for XB-EDA, respectively.Overall, the XB-EDA is the only experiment with total EIG spread greater than the WIG spread.As can be seen in Fig. 8, which presents ratios between the balanced, EIG and WIG spread with the total spread as a function of the zonal scale, this applies for every zonal wave number.The dominance of the EIG spread in the XB-EDA is most likely associated with the larger tropical spread in this experiment (see also Fig. 3) and it is also in agreement with Žagar et al. (2013), who presented the same conclusion for the 3 and 12 h forecast-error variances in an earlier model cycle.As discussed there, easterly propagating tropical modes represent the most important variability and largest forecast error source in the tropics.
Contrary to XB-EDA, the OBS-OBS-EDA and ST-EDA experiments contain more WIG spread than EIG spread for zonal wave numbers greater than 6 and 11, respectively (Fig. 8).If the spread in each experiment is normalized by OBS-EDA, it is found that XB-EDA is characterized by an increased EIG spread on all scales with respect to OBS-EDA while ST-EDA increases more the WIG spread than the EIG spread with respect to OBS-EDA (not shown).This result may be a consequence of the variance inflation as found in Žagar et al. (2011) for an ensemble Kalman filter system DART/CAM.
Figure 8 also shows that in all experiments the IG spread is dominant over the balanced spread for the zonal wave numbers greater than 10.At the shortest analysed scales, the IG spread makes about 70 % of the total spread.Although differences in the EIG and WIG spread may seem small when expressed in percentages, they illustrate a sensitive balance affected either dynamically (larger growth of forecast errors in tropical easterly modes) or artificially (variance inflation) with potentially major impacts on subsequent forecasts.

Using EDAs to define background error covariance matrices for assimilation
Background error statistics, the static B covariance matrix, computed from the OBS-, ST-and XB-EDAs are provided to the assimilation system and three T L 399L91T L 255 resolution analyses are computed for the period 5 October-15 November 2008.Because at the time these experiments were performed, the operational configuration was still using the static B matrix, the use of a "flow dependent" B matrix was not possible.Description of the computation of the static B matrix from an ensemble analysis can be found in Fisher (2003); see also Derber and Bouttier, 1999, for covariance modelling.The background error covariance matrix is modelled using coordinate transformations and spherical wavelet techniques (Fisher, 2003).In addition, a nonlinear, analytical balance is included in the covariance model (Fisher, 2003).The OBS, ST and XB analyses use, respectively, OBS-, ST-and XB-EDAs estimated B matrices.Diagnostics have been performed to assess the background error covariance on the assimilation system.The analysis experiments have the same name of the EDA experiments but bold fonts are used instead.
The first diagnostic presented is based on the observation influence (OI) (Cardinali et al., 2004;Cardinali, 2013) which quantifies the observational leverage in the analysis.The mean OI is the degree of freedom for signal, DFS, or total observation influence (Tukey, 1972;Velleman and Welsch, 1981;Wahba et al., 1995;Purser and Huang, 1993) divided by the total number of observation N: H TL is the linear observation operator and K is the gain matrix.OI and DFS depend on the assigned accuracy of the observations and background as well as the model itself which is a space and time propagator.The DFS quantifies the number of statistically independent directions constrained by each observation.Differences in the OI or DFS in the three assimilation experiments reflect differences in the B matri- ces.The OI is proved to be bounded between 0 and 1; 0 influence indicates that an observation has not had influence on the estimate but only the background counted whilst OI = 1 means that an entire degree of freedom has been devoted to fit that observation point.The OI can be gathered e.g. by observation type; in Fig. 9 the OI in OBS, ST and XB analyses is shown for different satellite and conventional observation types.Results indicate that XB shows the largest OI.
In particular, the largest OI increase is noticed for wind reporting observations (0.3 OBS, 0.5 ST and 0.7 XB), GPS-RO (0.2 OBS, 0.3 ST and 0.7 XB), AMSU-B radiances (0.2 OBS, 0.3 ST and 0.4 XB) and All-Sky SSMI radiances (0.1 OBS, 0.2 ST and 0.3 XB).The OI diagnostic indicates that when model errors are under-represented in the ensemble analysis, the background error statistics are also underestimated and the observations have smaller leverage in the assimilation procedure.XB analysis provides better observations fit (not shown) in agreement with the higher OI.XB DFS is 50 % larger than OBS and 30 % larger than ST DFS.
A measure of the consistency of the assimilation system is provided by the diagnostics on the background-error statistics computed in observation space (Talagrand, 2002;Desroziers et al., 2005).If the K gain matrix is consistent with the "true" covariances for background and observation errors, the innovation d and the analysis errors should be decorrelated from a statistical point of view.It can be simply shown (Desroziers et al., 2005) that the covariance between the analysis increment in observation space (H x a −H x b ) and the innovation vector (d), quantities archived during the assimilation procedure, should satisfy The assigned background error variances, H TL BH T TL (in observation space), are also archived; therefore the difference between the assigned and estimated background variances can be computed and averaged over the period of interest.In the context of linear estimation theory, a consistent unbiased analysis should result in no difference between the estimated and assigned background error variance.The following variance consistency check (VCC): measures the difference between the background error variances estimated from the analysis residuals (Desroziers et al., 2005) and the background error variances assigned from the ensemble analysis normalized with respect to the estimated ones.
The VCC computed for the period 5 October-15 November 2008 for OBS, ST and XB analyses shows small but non-zero values.Figure 10 shows the VCC for AMSU-A and AMSU-B, HIRS, SSMI, SCAT, vertical profilers and AMV: XB VCC is smaller than OBS and smaller than or similar to ST.

EDA-only based ensemble prediction system (EPS)
forecasts Buizza et al. (2008) proposed to use EDA-based perturbations in the ECMWF operational ensemble prediction system (EPS), and June 2010 EDA-based perturbations have been introduced in the EPS to improve the simulation of initial uncertainties (Isaksen et al., 2010).The replacement in the EPS of the evolved singular vectors with EDA-based perturbations improved substantially the EPS spread over the tropics, with a detectable impact in the early forecast range also over the extratropics.A positive impact was also detected on the EPS skill.
The three EDAs discussed in this work can be used to assess the sensitivity of EPS forecasts to the EDA configurations.Two types of ensembles were run: the first type included EDA-only perturbations, and the second type also included singular vectors (this is the configuration of the operational EPS).Since the impact of the EDA is more evident in the EDA-only type, attention will be focused mainly on them.EDA-only EPSs have been run with the perturbations defined using the OBS-, ST and XB-EDAs.All three EPS configurations included 50 perturbed and 1 unperturbed members, with variable resolution T L 399L62 between forecast day 0 and 10, and T L 255L62 between forecast day 10 and 15 (in uncoupled mode).All forecasts have been run with both the SPPT and the SPBS stochastic scheme as in the operational EPS (in other words, even the OBS-EPS that starts from the OBS-EDA-based perturbations that did not use any stochastic model, included stochastic perturbations in each of the 50 perturbed members).Forecasts have been run for 18 cases, with initial conditions from 12 October to 14 November 2008 every other day (with 12:00 UTC as initial time).In all ensemble configurations, following the methodology used in the ECMWF EPS operational at the time when the experiments were conducted (for more details, see Isaksen et al. (2010) and references therein), the EDA-based component of the 50 EPS initial perturbations have been constructed by (a) defining 10 EDA-based perturbations by computing the difference between each of the 10 EDA perturbed members and the unperturbed (control) member and by (b) adding and subtracting these EDA-based perturbations from the unperturbed analysis, defined by the ECMWF operational highresolution 4D-Var system.Since this procedure provides only 20 perturbations, EPS members 21-40 have the same initial EDA-based perturbations as members 1-20, and EPS members 41-50 have the same as members 1-10.The fact that up to three EPS members can use the same EDA-based perturbations is not a problem in the ensembles run with initial perturbations generated using both EDA-based perturbations and singular vectors (as is the case for the operational EPS), since 25 different SV-based initial perturbations are also used to generate the 50 positive and negative SV-based perturbations.For the EDA-only ensembles, the EPS members starting with the same EDA-based perturbation diverge, albeit in a slower way than the ensembles initialized by blending EDA-and SV-based perturbations, since each EPS member is integrated with different model error perturbations generated by the stochastic physic schemes.
The performance of an EPS is usually measured by a range of metrics that compare, in a statistical sense, the forecast probability distribution function with the verification (either the analysis, or observations).Skill metrics that are routinely used include the ranked probability score and skill score, the Brier score (Brier, 1950), the area under a relative operating characteristic and the ignorance skill score (Roulston and Smith, 2002).The reader is referred to Wilks (1995) for a general overview, and e.g. to Palmer et al. (2007) for a review of metrics used to assess the skill of the ECMWF EPS.In this work, attention has been focused on three aspects: firstly the ensemble reliability, i.e. the consistency between forecast probabilities and observed frequencies of occurrence, measured by the agreement between ensemble spread and ensemble-mean error; secondly the error of the ensemble-mean forecast; and thirdly the skill of probabilistic forecasts measured by the continuous rank probability skill score (CRPSS) and the ignorance score.Considering the first aspect, in a reliable ensemble, on average the spread of the system measured by the standard deviation should be equal to the average error of the ensemble mean.This property follows from the fact that in such a system, one ensemble member can be considered as the verification (Buizza et al., 2005;Palmer et al., 2006).Figure 11 shows the EPS spread (measured by the standard deviation) and the root-mean-square error of the ensemble mean forecast for the zonal wind at 850 hPa (verified against the operational high resolution analysis), computed over the Northern Hemisphere (20-80 • N, Fig. 10a), over the Southern Hemisphere (20-80 • S, Fig. 10b) and over the tropics (20 • S-20 • N, Fig. 10c).
Figure 11 shows that for all configurations the ensembles are underdispersive, especially over the tropics, indicating that EDA-only perturbations, if used as generated by the EDA and not re-scaled, are not sufficient to produce reliable ensemble forecasts.Among the EDA configurations, the XB-EPS has the largest spread, with differences evident up to about forecast day 10.Considering the error of the ensemble mean, Fig. 11 shows that the ensemble mean forecasts are very similar, almost undistinguishable for most of the forecast times, with the XB-EPS showing the smallest error for the forecast times when the spread is closer to the ensemble-mean root-mean-square error level (e.g. between forecast day 5 and 8 over the SH and between forecast day 5 and 10 over the tropics).
Considering the skill of probabilistic forecasts, all the metrics mentioned above have been considered, and since results are all consistent, only CRPSS and ignorance skill scores will be shown.The CRPSS is the equivalent of the mean squared error for single forecasts, and give a measure of the average distance between the forecast and observed distributions; the corresponding skill score, the CRPSSs have been computed using a climatological probabilistic forecast as reference (thus a perfect probabilistic forecast would score 1, and a forecast as skillful as climatology 0).The ignorance skill score is a logarithmic score defined using information theory (Roulston and Smith, 2002), based on the information deficit (or ignorance) in the forecast.According to Benedetti (2010), for probabilistic forecast systems the ignorance score and skill scores are more fundamental scores than the Brier score and skill scores, given that these latter are second-order approximations of the former.A clear advantage of these two scores compared to the Brier score and skill score, or the area under a relative operating characteristic, is that both the CRPSS and the ignorance skill score consider the whole forecast probability distribution function of forecast states, and not simply some specific event (e.g. the probability of a variable exceeding a certain threshold).Thus, they provide a more complete assessment than these latter two (Wilks, 1995).
Figure 12a shows that in terms of CRPSS, the three EPSs have a similar performance.The same conclusions can be drawn from the ignorance skill score (Fig. 12b).However, it is worth noting that the small but consistent improvement at day 2 and 3 and day 6 as shown in Fig. 12a and b is observed for all parameters at different levels.Similar conclusions could be drawn by considering different performance metrics for these variables, e.g. the Brier skill score or the area under the relative operating characteristic for the probabilistic prediction of dichotomous events such as "wind anomalies above or below the climatological standard deviation" (not shown).Concluding, these results based on EDA-only ensembles indicate that the use of the XB-EDA in the EPS would lead to a slightly better match between spread and error, i.e. a better reliability, and to similar forecast skill as the use of the ST-and OBS-EDAs The EPS experiments based on EDAbased perturbations and singular vector perturbation indicate smaller spread differences, detectable only up to forecast day 5 instead of day 10, and practically no differences in skill (not shown).

Conclusions
In this paper, the representation of model error in the ECMWF 4DVar Ensemble Data Assimilation (EDA) system by additively perturbing the background field is shown.The method follows the idea used in the Meteorological Service of Canada (MSC) ensemble Kalman filter (Houtekamer et al., 2005) and in the ARPEGE 4D-Var system (Raynaud et al., 2012) to account for model uncertainties in the EDA by perturbing the model background field.In the MSC ensemble Kalman filter the perturbation magnitude is globally constant (but the representation of other errors is also considered) whilst at Météo-France the perturbation, as a function of latitude band and model level, is multiplied to the background field (multiplicative approach).Here, the additive approach is presented.
The idea behind it is that a large fraction of model error is represented by the short-range forecast error.Thus perturbing the 12 h model forecast would include error sources introduced by dynamics, spatial and temporal discretization and, in the case of assimilation, linearized physical processes and the misspecification of the probability distribution of errors in the observations and the numerical background model.The proposed methodology, contrary to the stochastically perturbed parameterization tendency scheme, does not require routine diagnostic and tuning.
Ensembles of data assimilation with different representations of model error have been compared.In particular, two more EDAs are examined, all with the same methodology to represent the observation uncertainties but different techniques to account for model uncertainties.One model error representation technique is based on the assumption that random model errors due to the physical process parameterizations are the main model error source (ST-EDA).In the ST-EDA, stochastic perturbations are added to the physical model tendency at each model time step.The other ensemble considered (OBS-EDA) only includes an observation error representation, which implicitly modifies the background fields in the assimilation cycling process.In the additive XBmethod, the magnitude of the perturbation is calculated by comparing the variance of the innovation vector of the high resolution analysis system with the ensemble data assimilation variance in which the ensemble data assimilation is produced by only perturbing the observations.The perturbation is a function of latitude band, vertical model level and the model parameter.
Results have shown that the EDA generated using XBmethod accounts for the largest spread, followed by the ST-EDA with stochastically perturbed parameterization tendencies and the OBS-EDA with observation-only perturbations.The increase of spread depends on the location: the stratosphere accounts for the largest increase of spread.In the tropics, the increase occurs mainly in the mid-troposphere, while in the extratropics spread enhancement is detected at all model levels and in the Southern Hemisphere in particular.The largest difference compared to other EDA configurations is found in the extratropics where all but the XB-EDA are not able to produce significant spread.However, the spread structure is quite similar among the different methods confirming that a large part of the model error sources are affecting the short-range forecast.Comparison with respect to the uncycled ensemble analysis (OBS-OBS EDA) shows that the background perturbation adds spread from high to low levels in the troposphere and into the extratropics in both Northern and Southern hemispheres.Normalmode diagnostics compare the percentage of balanced and unbalanced spread in the four ensembles XB-, ST-, OBS-OBS-and OBS-EDA.It was found that various experiments contain around 40 % of their spread in the balanced modes with ST-EDA being the least balanced.Both ST-and XB-EDA add both balanced and IG spread homogeneously at all scales when compared with OBS-EDA.An interesting difference between XB-and ST-EDA is found in the distribution of the added easterly IG spread relative to the westerly IG spread.XB-EDA adds more EIG spread on all scales whereas the ST-EDA increases the WIG more than the EIG spread when compared with OBS-EDA.Overall, XB-EDA is the only experiment which contains more spread in the easterly propagating IG modes than in the westerly IG component.It is believed that this is due to the tropical flow properties which are strongly influenced by the easterly propagating Kelvin waves.Although more investigation is needed to confirm this interpretation, the increase of spread in the westerly IG mode for the ST is likely to be related to the inflation of the model background error variances (as was found by Žagar et al., 2011, for the DART-CAM ensemble).In fact, it should be kept in mind that OBS-and ST-EDAs include a global inflation factor of the background error variances that penalizes the model background further with respect to the observations.The inflation is a spatial and temporal constant that does not depend on synoptic weather developments but only intends to introduce larger ensemble spread.
The covariance matrices produced by the three ensemble analyses have been provided to a higher resolution data assimilation system.The diagnostic performed on the three resulting analyses show that the largest observation influence (OI) is obtained from the assimilation system with XB-EDA background statistics.This is due to the larger background error variances.When the analysis residual diagnostic is applied to investigate the consistency of the assimilation systems considered, systematical smaller differences are found with XB showing closer agreement between the assigned and estimated background error variances.
Finally, the three EDAs have been used in ensemble prediction mode.Since June 2010, EDA-based perturbations have been used with singular vectors to simulate initial uncertainties in the ECMWF EPS (Buizza et al., 2008;Isaksen et al., 2010).Results have indicated that the use of the XB-EDA in the EPS would lead to the largest spread, with differences evident up to about forecast day 10 and with the smallest error for forecast times from day 5 to 8 over Southern Hemisphere and from day 5 to 10 over the tropics when the spread is closer to the ensemble mean root-mean-square error.In terms of EPS forecast skill, very small but consistent improvements up to day 7 have been detected when, in particular, the ignorance score is used.
In conclusion, the XB-method discussed in this work is shown to be a valuable alternative of the method used in the current ECMWF EDA to simulate model uncertainty in the ensemble analysis.It accounts for different sources of error coming from the dynamics, the parameterizations, the linearization and interpolation schemes and it is easier to tune and maintain.The tuning of the perturbation is performed automatically every 3 days from a 3-day sample of the high resolution operational innovations.A possible extension of the work presented would be to combine the background perturbation method with the stochastic physical tendency perturbation method or other methodologies that are considered to simulate longer range random model error sources.The estimation of the background perturbations magnitude should, in this case, be achieved by comparing the innovations of the high resolution analysis system with the EDA in which not only the observations are perturbed but also e.g.stochastic perturbations are added to the physical model tendency, that is, the ST-EDA.
The Supplement related to this article is available online at doi:10.5194/npg-21-971-2014-supplement.

Figure 1 .
Figure 1.Perturbation modulation factor f as a function of latitude band (North Hemisphere solid black line, Southern Hemisphere dotted line, and tropics solid grey line) and vertical model levels (model level 1 is 0.01 hPa and 91 is on average ∼ 1000 hPa) for (a) zonal wind, (b) temperature, and (c) humidity.The modulation factor is estimated over a 3-week period.

Figure 3 .
Figure 3. Zonally averaged cross section for the period 5 October to 15 November 2008 for the u component of the EDA spread.(a) OBS-EDA, (b) OBS-OBS-EDA, (c) ST-EDA and (d) XB-EDA.The vertical coordinate is model level.Contours are in m s −1 .

Figure 4 .Figure 5 .
Figure 4. Zonally averaged cross section from 5 October 2008 to 15 November 2008 for the u component of the EDA spread.(a) Difference in spreads between the OBS-EDA and the OBS-OBS-EDA spreads, (b) difference between the XB-EDA and the OBS-EDA spreads.Vertical coordinate is model level.Contours in m s −1 .

Figure 7 .
Figure 7. Time-averaged ensemble spread in the balanced and inertio-gravity modes for (a) XB-EDA, (b) ST-EDA, (c) OBS-and (d) OBS-OBS EDA.Black curves correspond to the spread ated with balanced modes (ROT), grey curves to spread due to easterly propagating inertia-gravity (EIG) modes whereas dotted curves represent spread due to westerly propagating inertia-gravity modes (WIG).

Figure 9 .
Figure9.Global average observation influence (OI) for the different observation types assimilated in the XB (black line), ST (dark grey line), OBS (light grey line) 4DVar analyses.AMSU-A and AMSU-B are microwave radiances, AIRS and IASI and HIRS infrared radiances, SSMI microwave imager radiance, GPS-RO satellite GPS radio occultation, OZONE retrieval, SCAT retrieved wind information from microwave scatterometer, atmospheric motion vector (AMV) from geostationary cloud imagery and vertical profiler consists of wind from radiosonde, pilot, aircraft and wind profiler observations.

Figure 11 .
Figure 11.Averaged spread measured by the standard deviation (solid lines) and error of the ensemble mean (lines with symbols) of the EPS run with EDA-only perturbations generated from OBS-EDA (light grey), ST-EDA (dark grey) and XB-EDA (black).Results refer to the zonal wind component at 850 hPa over the (a) Northern Hemisphere (20-80 • N), (b) Southern Hemisphere (20-80 • S) and (c) tropics (20 • S-20 • N).The average has been computed considering 18 cases, each with 51-member EPS forecast with initial condition from 12 to 14 November 2008, every other day (12:00 UTC only).

Figure 12 .
Figure12.Average CRPSS (a) and ignorance score (b) of EPS run with EDA-only perturbations generated from OBS-EDA (light grey), ST-EDA (dark grey) and XB-EDA (black) for temperature at 850 hPa over the Southern Hemisphere.The average has been computed considering 18 cases, 51-member EPS forecast with initial condition from 12 to 14 November 2008, every other day (12:00 UTC only).