Outline

1. Ensemble Data AssimilationMassimo Bonavitamassimo.bonavita@ecmwf.int room 1014; ext. 2627ECMWFContributions from:Lars Isaksen, Elias Holm, Mike Fisher, Laure Raynaud

2. Outline • The Ensemble Data Assimilation method • What do we do with the EDA? • Hybrid DA • Use of EDA variances in 4D-Var • Use of EDA co-variances in 4D-Var • Summary

3. The EDA method For a linear system the data assimilation update is: Under the assumptions of statistically independent background (Pb), observation (R) and model errors (Q), the evolution of the system error covariances is given by:

4. The EDA method Consider now the evolution of the same system where we perturb the observations and the forecast model with random noise drawn from the respective error covariances: where ηk~N(0,R), ζk~N(0,Q). If we define the differences between the perturbed and unperturbed state and , their evolution is obtained by subtracting the unperturbed state evolution equations from the perturbed ones:

5. The EDA method i.e., the perturbations evolve with the same update equations of the state (Kalman gain and model operator). What about the errors? If we take the statistical expectation of the outer product of the perturbations:

6. The EDA method • These are the same equations for the evolution of the system error covariances: • provided that: • The applied perturbations ηk, ζk have the required covariances (R, Q); • At some stage in time

7. The EDA method What does all this mean in practice? We can use an ensemble of perturbed assimilation cycles to simulate the errors of our reference assimilation cycle; The ensemble of perturbed DAs should be as similar as possible to the reference DA (i.e., same or similar K matrix) The applied perturbations ηk, ζk must have the required error covariances (R, Q); There is no need to explicitly perturb the background xb

8. X1b(tk) X2b(tk) X1b(tk+1) X2b(tk+1) X1a(tk) X2a(tk) Analysis Analysis Forecast Forecast y+ε1o y+ε2o Boundary pert. 1 Boundary pert. 2 The EDA method ε2m ε1m ………………………………………………………….

9. The EDA method • 10 ensemble members using 4D-Var assimilations • T399 outer loop, T95/T159 inner loops. (Reference DA: T1279 outer loop, T159/T255/T255 inner loops) • Observations randomly perturbed • Observations error correlations taken into account • SST perturbed with realistically scaled structures • Model error represented by stochastic methods (SPPT, Leutbecher, 2009) • All 107 conventional and satellite observation types used

10. Outline • The Ensemble Data Assimilation method • What do we do with the EDA? • Hybrid DA • Use of EDA variances in 4D-Var • Use of EDA co-variances in 4D-Var • Summary

11. What do we do with the EDA? • The EDA simulates the error evolution of the 4DVar analysis cycle. As such it has two main applications: • Provide aflow-dependent estimate of analysis errors to initialize the ensemble prediction system (EPS) • Provide a flow-dependent estimate of background errors for use in 4D-Var assimilation

12. EVO-SVINI EDA-SVINI EVO-SVINI EDA-SVINI What do we do with the EDA? Improving Ensemble Prediction System by including EDA perturbations for initial uncertainty The Ensemble Prediction System (EPS) benefits from using EDA based perturbations. Replacing evolved singular vector perturbations by EDA based perturbations improve EPS spread, especially in the tropics. The Ensemble Mean has slightly lower error when EDA is used. N.-Hem. Tropics Ensemble spread and Ensemble mean RMSE for 850hPa T

13. What do we do with the EDA? • The EDA simulates the error evolution of the 4DVar analysis cycle. As such it has two main applications: • Provide a flow-dependent estimate of analysis errors to initialize the ensemble prediction system (EPS) • Provide a flow-dependent estimate of background errors for use in 4D-Var assimilation

14. What do we do with the EDA? For Gaussian error distributions and linear model and observation operators: xt= Mxt-1+ ηη~N(0,Q)(1) yt= Htxt + εtεt~N(0,R) (2) The Kalman Filter (KF) equations give the state forecast distribution xt|y1:t-1~N(xbt,Pbt): xbt= Mxat-1(3) Pbt= MPat-1M+Q(4) The analysis distribution xt|y1:t~N(xat,Pat) is given by: xat = xbt+ Kt(yt- Htxbt)(5) Pat= (I - KtHt)Pbt(6) Kt= PbtHTt(Rt+ HtPbtHTt)-1 (7)

15. What do we do with the EDA? • The KF propagates information from past observations through the background state xbtand its error covariance Pbt • 4D-Var does not cycle the state error covariance information, only the state estimate • 4D-Var is equivalent to a Kalman filter in which the background error covariance matrix is reset to some static matrix Bat every analysis update (12h in the current ECMWF implement.) • Two possible remedies: • Run 4DVar over a long enough window so that the influence of the initial state and errors on the final state analysis is negligible (Long window Weak constraint 4D-Var, talk by Y. Tremolet)

16. Analysis experiment started with satellite data reintroduced on 15th August 2005 Memory of the initial state disappears after approx. 3days From Mike Fisher

17. What do we do with the EDA? • An analysis window of ≥3 days would allow 4D-Var analysis at final time to be independent of initial state and error estimates • However: • An effective model errorparameterization must be applied to reconcile model and measurements over a long analysis time window (i.e., weak constraint 4D-Var) • We would still lack an estimate of analysis errors

18. What do we do with the EDA? • Two possible ways around this: • Run 4DVar over a long enough window so that the influence of the initial state and errors on the final state analysis is negligible • Use a sequential method (EDA, EnKF) to cycle error covariance information.

19. Outline • The Ensemble Data Assimilation method • What do we do with the EDA? • Hybrid DA • Use of EDA variances in 4D-Var • Use of EDA co-variances in 4D-Var • Summary

20. Hybrid DA Hybrid approx.: Use flow-dependent state error estimates (from an EnKF/EDA system) in a 3/4D-Var analysis algorithm This solution would: • Integrate flow-dependent state error covariance information into the variational analysis • Keep the full rank representation of B and its implicit evolution inside the assimilation window • More robust than pure EnKF for limited ensemble sizes and large model errors • Allow (eventual) localization of ensemble perturbations to be performed in state space; • Allow for flow-dependent QC of observations

21. Hybrid DA Hybrid systems • Ensemble perturbations can be used in a 3-4DVar analysisin a number of different ways: • Use ensemble covariances to estimate B matrix at the start of a variational analysis window (EDA method; alpha control variable method, Lorenc 2003) • Use of ensemble covariances over the all 4DVar analysis window(En4DVAR, E4DVAR, 4D-EN-VAR) in lieu of TL/ADJ model

22. Hybrids: α control variable Conceptually add a flow-dependent term to the climatological B matrix: Bc is the static, climatological covariance Pe○ Cloc is the localised ensemble covariance In practice this is done through augmentation of control variable: and introducing an additional term in the cost function: • Alpha control variable method (Met Office, NCEP/GMAO)

23. Hybrids: EDA In the ECMWF 4D-Var, the B matrix is defined implicitly in terms of a transformation from the background departure (x-xb) to a control variable χ: (x-xb) = Lχ So that the implied B=LLT. In the current wavelet formulation (Fisher, 2003), the variable transform can be written as: T is the balance operator Σb is the gridpointvariance of background errors Cj(λ,φ) is the vertical covariance matrix for wavelet index j ψj are the set of radial basis function that define the wavelet transform.

24. Hybrids: EDA Cj(λ,φ)are full vertical covariance matrices, function of (λ,φ). They determine both the horizontal and vertical background error correlation structures; In standard 4D-Var T and Cjare computed off-line using a climatology of EDA perturbations. Σb is computed by random sampling of the static B matrix (randomization procedure, Fisher and Courtier, 1995) How do we make this error covariance model flow-dependent? We look for flow-dependent EDA estimates of Σband Cj(λ,φ)

25. Outline • The Ensemble Data Assimilation method • What do we do with the EDA? • Hybrid DA • Use of EDA variances in 4D-Var • Use of EDA co-variances in 4D-Var • Summary

26. Use of EDA variances in 4DVar • We want to use EDA perturbations to simulate 4DVar flow-dependent background error covariance evolution • We start with the diagonal of the B matrix, Σb

27. Use of EDA variances in 4DVar What do raw ensemble variances look like? Spread of Vorticity FG t+9h 500hPa

28. Use of EDA variances in 4DVar Define: Veda = EDA sampled variance V* = True error variance Then the estimation error can be decomposed: Veda-V* = E[Veda-V*] + (Veda-E[Veda]) systematic random Similarly for mean square of the estimation error: E[(Veda-V*)2] = (E[Veda-V*])2 + E[(Veda-E[Veda])2] systematic random We should try to minimize both!

29. Use of EDA variances in 4DVar • Noise level is due to sampling errors: 10 memberensemble • EDA is a stochastic system: variance of variance estimator ~ 1/Nens • We need a system to effectively filter out noise from first guess ensemble forecast variances: Reduce the random component of the estimation error

30. Use of EDA variances in 4DVar “Mallat et al.: 1998, Annals of Statistics, 26,1-47” Define Ge(i) as the random component of the sampling error in the estimated ensemble variance at gridpointi: Then the covariance of the sampling noise can be shown to be a simple function of the expectation of the ensemble-based covariance matrix: (1)

31. Use of EDA variances in 4DVar A consequence of (1) is that: (2) i.e., sampling noise is smaller scale than background error. If the variance field varies on larger scales then the background error we can use a spectral filter to disentangle noise error from the sampled variance field

32. Use of EDA variances in 4DVar

33. Use of EDA variances in 4DVar • There is indeed a scale separation between signal and sampling noise! • Truncation wavenumber is determined by maximizingsignal-to-noiseratio of filtered variances (details in Raynaud et al., 2009; Bonavita et al., 2011) • Optimal truncation wavenumberdepends on parameter and model level

34. Use of EDA variances in 4DVar Raw Ensemble StDev VO ml64 Filtered Ensemble StDev VO ml64

35. Use of EDA variances in 4DVar Raw Ensemble StDev Spec. Hum. ml64 Filtered Ensemble StDev Spec. Hum. ml64

36. Use of EDA variances in 4DVar Operational StDev Random. Method (Fisher & Courtier, 1995) Sample of static B errors Filtered Ensemble StDev VO ml64

37. Use of EDA variances in 4DVar Should we also do something about the systematic component of the estimation error? E[(Veda-V*)2] = (E[Veda-V*])2 + E[(Veda-[Veda])2] systematic random A statistically consistent ensemble satisfies: (1-1/Nens)-1<Ens_Var> = (1+1/Nens)-1<Ens_Mean_Square_Error>

38. Use of EDA variances in 4DVar Vorticity ml 78 (~850hPa) Ensemble Error Ensemble Spread Spread - Error

39. Use of EDA variances in 4DVar • To get statistically consistent EDA variances we need to perform an online calibration (Ensemble Variance Calibration; Kolczynsky et al., 2009, 2011; Bonavita et al., 2011) • Calibration factors depend onlatitude bands and parameter • Calibration factors are also flow-dependent, i.e. depend on the size of the expected error

40. Analysis Forecast Analysis Forecast SST+εiSST EDA Cycle xb+εib i=1,2,…,10 xb xa xf+εif xf xa+εia y+εio EDA scaled Var Variance post-process Variance Rescaling Variance Filtering EDA scaled variances εifraw variances 4DVar Cycle

41. Use of EDA variances in 4DVar • How does the flow-dependent structure of the EDA variances affect the 4D-Var analysis? • Single Observation experiments

42. Use of EDA variances in 4DVar Z500 Geopotential (shaded) and MSLP valid 30-09-2010 at 21Z

43. Use of EDA variances in 4DVar Vorticity spread ml=78 (850hPa) EDA “Randomization method”

44. Use of EDA variances in 4DVar Single obs. experiment: Tobs-Tfg=+1K, (34N,74W), 900 hPa EDA “Randomization method”

45. Use of EDA variances in 4DVar Single obs. experiment: Tobs-Tfg=+1K, (34N,74W), 900 hPa 4D-Var+EDA Standard 4D-Var Tana-Tfg 900 hPa VOana-VOfg 850 hPa

46. Use of EDA variances in 4DVar • The observation weight in the analysis is increased in the area of large background uncertainty: EDA Randomiz. ΔT ΔVo • The EDA analysis increments show a degree of flow-dependency

47. Use of EDA variances in 4DVar What is the impact of EDA variances on scores? • CY36R4, T1279L91 • ffg8 20100111 - 20100331 (control: fezj): WINTER • ffge20100802 – 20101030 (control: 0051): SUMMER

48. Geopotential RMSE reduction winter summer Blue=☺

49. Temperature RMSE reduction winter summer Blue=☺

50. Outline • The Ensemble Data Assimilation method • What do we do with the EDA? • Hybrid DA • Use of EDA variances in 4D-Var • Use of EDA co-variances in 4D-Var • Summary