Maximizing Forecast Accuracy with Maximum Likelihood Ensemble Filter

1. NOAA/NCEP/EMC, Camp Springs, MD 3 November 2005The Maximum Likelihood Ensemble Filter development at the Colorado State UniversityMilija ZupanskiCooperative Institute for Research in the AtmosphereColorado State UniversityFort Collins, CO 80523-1375E-mail: ZupanskiM@CIRA.colostate.edu In collaboration with:Colorado State University: D. Zupanski, S. Fletcher, D. Randall, R. Heikes G. Carrio, W. Cotton Florida State University: I.M. Navon, B. Uzunoglu NOAA/NCEP: Zoltan Toth, Mozheng Wei, Yucheng Song Computational Support: NCEP IBM SP (frost) ,NCAR SCD (bluesky)

2. Outline • Motivation • Hessian preconditioning • Maximum Likelihood Ensemble Filter (MLEF) • Preliminary results • Double-resolution MLEF • Future research directions

3. Motivation • Uncertainties - Assign a degree of confidence in the produced analysis/forecast - Transport in time of the forecast/analysis state vector + uncertainty • Universality of assimilation/prediction - Same system can be used in wide range of applications - Portability • Single assimilation/prediction system - Complete feed-back between uncertainties - Easy to maintain and upgrade • Fewer assumptions/restrictions - Non-differentiable operators (discontinuity) - Highly nonlinear operators (microphysics, clouds) • User-friendly: Non-experts can use DA and EF - Allow more people to enjoy the benefit of new research

4. Hessian Preconditioning Variational cost function Hessian Inverse Hessian Ideal preconditioning is a square-root of the Hessian matrix Perfect preconditioning: Hessian in minimization space is an identity matrix ! Hessian condition number in variational data assimilation ~ 60-100 !

5. Hessian preconditioning One-iteration minimization for quadratic cost function !

6. Nonlinear observation operators • The KF, EKF, EnKF solution identical to the minimization of a quadratic cost-function • Two strategies for nonlinear observation operators: (1) Use linear KF solution, combined with nonlinear operators in covariance calculation - EnKF algorithms - EKF • Directly search for nonlinear solutionbyminimizing non-quadratic cost function • - Maximum Likelihood Ensemble Filter (MLEF) – conditional mode • - Iterative KF – conditional mean

7. Maximum Likelihood Ensemble Filter (MLEF)(Zupanski 2005, MWR; Zupanski and Zupanski 2005, MWR) • Estimate of the conditional mode of the posterior PDF • Ensembles used to estimate the uncertainty of the conditional mode • Non-differentiable minimization with Hessian preconditioning: • (Generalized conjugate-gradient and BFGS quasi-Newton algorithms) • Augmented control variable: initial conditions, model bias, empirical parameters, boundary conditions • Related to: (i) variational data assimilation, • (ii) Iterative Kalman filters, and • (iii) Ensemble Transform Kalman Filter – ETKF • Not sample based • Reduces to Kalman filter for linear operators and Gaussian PDF

8. MLEF Framework Use the forecast (prior) error covariance square-root State-space dimension Ensemble size Minimize cost function in the subspace spanned by ensemble perturbations Similar to variational, however: • Non-differentiable iterative minimization with superior preconditioning • Solution in ensemble subspace (reduced rank) • Analysis uncertainty estimate

9. Hessian preconditioning in MLEF Change of variable (Hessian preconditioning) Ensemble-size matrix Background (first guess) k – iteration index ETKF transformation utilized in Hessian preconditioning

10. “Gradient” calculation in MLEF • Not a gradient, rather a directional derivative in the direction of ensemble perturbation • Important for nonlinear/non-differentiable operators “Gradient” calculation Innovation vector Observation component of the gradient is an innovation vector projection onto the ensemble perturbations

11. Analysis (posterior) error covariance Analyisis (minimizer) • Analysis error covariance estimated from minimization algorithm • At the minimum (xmin=xa) use • Inverse Hessian = Analysis error covariance • Justification for the assumption • Accurate minimization implies small distance between the analysis and the truth • Good Hessian preconditioning allows efficient and accurate minimization • By monitoring minimization, assure the calculated solution is close to the true minimum

12. Analysis Error Covariance in KdVB model Cycle No. 1 Cycle No. 4 Cycle No. 7 Cycle No. 10 j i • Initial error covariance noisy, but quickly becomes spatially localized • No need to force error covariance localization Model dynamics forces adequate localization of uncertainties !

13. MLEF with CSU global shallow-water model(Heikes and Randall 1995, MWR; Zupanski et. al. 2005, Tellus) Height analysis increment [xa-xf] Height RMS error [xa-xt] Initially noisy random perturbations quickly become smooth:Consequence of error covariance localization by dynamics

14. Error covariance localization (linear framework) Lyapunov vector Possible explanation for error covariance localization: Dynamic localization of Lyapunov vectors

15. Assimilation of real boundary-layer cloud observations using the LES RAMS model • 23 2-h DA cycles: 18UTC 2 May 1998 – 00 UTC 5 May 1998 (Mixed phase Arctic boundary layer cloud at Sheba site) • Experiments initialized with typical clean aerosol concentrations • May 4 was abnormal: high IFN and CCN above the inversion • x= 50m, zmax= 30m (2d domain: 50col, 40lev), t=2s, Nens=48 • Sophisticated microphysics in RAMS/LES, Prognostic IFN, CCN • Control variables: _il, u, v, w, N_x, R_x (8 species), IFN, CCN (dim= 22 variables x 50 columns x 40 levels = 44,000) • Radar/lidar/aircraft observations (retrievals) of IWP, LWP G. Carrio, W. Cotton

16. LIQUID WATER CONTENT ASSIMILATION CONTROL Vertical structure of the analysis: LWC EXP Vertically integrated observations: LWP VERIF Better timing of maxima G. Carrio, W. Cotton

17. ICE FORMING NUCLEI (IFN) CONCENTRATION IFN below inversion as cloud forms IFN above the inversion, as observed 22Z Independent observation G. Carrio, W. Cotton

18. Operational resolution forecast model for the control • Low resolution forecast model for the ensembles • Cost function defined in operational resolution • Minimization in ensemble subspace • Motivation - 3DVAR/4DVAR operational systems • - Computational savings • - Fewer number of ensembles required than in the full operational setup • - More adequate number of degrees of freedom in the ensembles Double-resolution MLEF framework(THORPEX)

19. Double-resolution MLEF:Forecast step Interpolation operator from the operational (high) to coarse (low) resolution Forecast error covariance column vector (low resolution): MC – low resolution forecast model M – operational resolution (control) forecast model xC – low resolution state vector x – operational resolution state vector

20. Operational resolution cost-function Double-resolution MLEF:Analysis step Change of variable (low resolution) Low resolution observation operator Interpolated from operational resolution Directional gradient (ensemble subspace) For double-resolution MLEF, need an interpolation operatorC, low resolution modelMC, and low resolution observation operatorHC

21. THORPEX related development of the MLEF • MLEF with T6228 GFS and SSI • Code development completed and debugged • GFS model + SSI (interpolation from model to observations) • Currently tested PREPBUFR observations, will include satellite/radar/lidar Compare MLEF with other EnKF methods • Evaluate if dynamical localization holds • Ensemble size, robustness of the algorithm • Code efficiency: script driven algorithm, exploits the NCEP code structure • Model bias and parameter estimation • Capability included in the current MLEF/GFS version • Reduce the large number of degrees of freedom (by a projection operator) • Double-resolution MLEF • Test this capability, using two GFS resolutions • Evaluate the ensemble size issue

22. Other Research and Future • Development of a fully non-Gaussian algorithm • Allow for non-Gaussian state variable errors (initial conditions, empirical parameters) • Generalized algorithm with a list of PDFs CloudSat assimilation • Observation information content • Relative information content from various observation types/groups • Microscale models • Boundary layer, 50m-500m horizontal resolution • Probabilistic transfer and interaction between scales • Carbon data assimilation • Exploit assimilation of new measurements (OCO-Orbiting Carbon Observatory) • MLEF with super-parameterization • Assimilation of clouds and precipitation observations – MMF (Multiscale Modeling Framework) • NASA GEOS + super-parameterization • Climate models and predictability

23. Thank you !

24. ISSUE A: Hessian preconditioning vs. low-rank Starting point 2 Starting point 1 xmin J=const • For high Hessian condition numbers ~50-100, minimization success is unpredictable • Low-rank assumption in ensemble DA may miss important perturbations (directions)

25. ISSUE B: Error covariances Dynamical localization vs. prescribed structure • Addressed by advanced methods • Computational efficiency • Moisture related variables (microphysics) Weak correlations across frontal zone L x1 Strong correlations along frontal zone x2

26. Forecast (prior) error covariance • Control vector is the most likely forecast • Square-root used in the algorithm (full covariance can be calculated) • No sampling of error covariance • Provides dynamic continuity between the analysis and forecast

27. Atmospheric observations have non-Gaussian statistics: • - precipitation • - specific humidity (moisture) • - ozone • - cloud droplet concentration • - microphysical variables • Atmospheric state variables have non-Gaussian statistics: - specific humidity (moisture) - ozone • - microphysical variables • - concentrations (clouds, aerosols) Non-Gaussian MLEF framework Gaussian data assimilation framework is generally used Need to evaluate the impact of this assumption

28. MLEF approach • Define non-Gaussian cost function from the conditional PDFs • Minimize such defined cost function – calculate the conditional mode • Algorithmically, define a list of PDFs, follow the appropriate code branching • Start with relatively well known Lognormal PDF • Examine Gaussian assumptions using a non-Gaussian mathematical framework • Fletcher and Zupanski (2005a-SIAM J.Appl.Math.; 2005b-J.Roy.Statist.Soc.B), Non-Gaussian MLEF framework mode median mean • Mode, Mean, Median are identical in Gaussian (or any symmetric) PDF • Mode, Mean, Median are all different in Lognormal (or any skewed) PDF

29. Log-Normal PDF Gaussian PDF Non-Gaussian MLEF framework Lognormal errors are multiplicative: Gaussian errorsare additive:

30. Non-Gaussian MLEF experiment with SWM:Lognormal height observation errors • Assume: • Gaussian prior PDF • Lognormal observation PDF (height) • Lognormal observation operator H(x)=exp[a(x-b)] Minimize mixed Normal-Lognormal cost function: Additional Lognormal observation term Gaussian prior PDF Gaussian-like lognormal observation term [for ln(x)] Higher nonlinearity of the cost function compared to the pure Gaussian

31. sLogn= 1.0 ; m= 1.27 sGauss= 1.5 m sLogn= 1.0 ; m = 3.15  sGauss= 2.5 m Impact of Lognormal observation errors:Analysis RMS errors Success of the Gaussian MLEF depends on the observation statistics

32. sLogn= 1.0 ; m= 1.27 sGauss= 1.5 m sLogn= 1.0 ; m = 3.15  sGauss= 2.5 m Impact of Lognormal observation errors:Innovation histogram N(0,1) Innovation statistics is significantly impacted by the PDF framework

33. MLEF Analysis Step Use the forecast (prior) error covariance square-root Minimize cost function in the subspace spanned by ensemble perturbationspif Similar to variational, however: • Non-differentiable iterative minimization with preconditioning • No differentiability assumption: works for all bounded operators • Generalized gradient, generalized Hessian • Perfect preconditioning for quadratic cost function • Solution in ensemble subspace • Reduced dimensions of the analysis correction subspace • Focus on unstable, growing perturbations in the analysis • Search for the vectors (ensemble perturbations) that span the attractor subspace

34. Information Content AnalysisNASA GEOS-5 Single Column Model ds measures effective DOF of an ensemble-based data assimilation system (e.g., MLEF). Useful for addressing DOF of the model error. Dusanka Zupanski, CIRA/CSU Zupanski@CIRA.colostate.edu

35. A : Reduction of uncertainty (0-) MLEF 450 ens Bayesian MLEF (full rank) successfully reproduces Bayesian inversion results