# Presentation: 573127

## NOAA/NCEP/EMC, Camp Springs, MD 3 November 2005The Maximum Likelihood Ensemble Filter development at the Colorado State UniversityMilija ZupanskiCooperative Institute for Research in the AtmosphereColorado State UniversityFort Collins, CO 80523-1375E-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)

**Notes:**
Outline slide

## Outline

**Motivation  **
**Hessian preconditioning**
**Maximum Likelihood Ensemble Filter (MLEF)**
**Preliminary results**
**Double-resolution MLEF**
**Future research directions**

**Notes:**
Outline slide

## 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

**Notes:**
Outline slide

## Hessian Preconditioning 

**Variational cost function **

**Hessian **

**Ideal preconditioning is a square-root of the Hessian matrix**

**Inverse Hessian **

**Perfect preconditioning: Hessian in minimization space is an identity matrix !**

**Hessian condition number in variational data assimilation ~ 60-100 !**

**Notes:**
Outline slide

## Hessian preconditioning 

**One-iteration minimization for quadratic cost function !**

**Notes:**
Outline slide

## Nonlinear observation operators

*** ******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*** solution by ***minimizing non-quadratic cost function***
- 		- Maximum Likelihood Ensemble Filter (MLEF) – *conditional mode*
- 		- Iterative KF – *conditional mean* **	**** **

**The KF, EKF, EnKF solution identical to the minimization of a quadratic cost-function**

**Notes:**
Outline slide

## 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

**Notes:**
Outline slide

## MLEF Framework

**Minimize cost function in the subspace *****spanned by ensemble perturbations***

**Use the forecast (prior) error covariance square-root**

**Similar to variational, however:**** **
**Non-differentiable iterative minimization with superior preconditioning**
**Solution in ensemble subspace (reduced rank)**
**Analysis uncertainty estimate**

**Ensemble size**

**State-space dimension**

**Notes:**
Outline slide

## Hessian preconditioning in MLEF

**Change of variable (Hessian preconditioning)**

*k* – iteration index

- Ensemble-size matrix

- Background (first guess)

**ETKF transformation utilized in Hessian preconditioning**

**Notes:**
Outline slide

## “Gradient” calculation in MLEF

**Not a gradient, rather a directional derivative in the direction of ensemble perturbation**
**Important for nonlinear/non-differentiable operators**

**Observation component of the gradient is **
**an innovation vector projection onto the ensemble perturbations**

**“****Gradient” calculation**

- Innovation vector

**Notes:**
Outline slide

## Analysis (posterior) error covariance

**Analysis error covariance estimated from minimization algorithm**
**At the minimum (*****x****min***=*****x****a***) 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**

- Analyisis (minimizer)

**Notes:**
Outline slide

## j

** ****Initial error covariance noisy, but quickly becomes spatially localized**
** ****No need to force error covariance localization **

**Cycle No. 4**

**Cycle No. 7**

**Cycle No. 10**

**Cycle No. 1**

*j*

*i*

**Analysis Error Covariance in KdVB model**

**Model dynamics forces adequate localization of uncertainties ! **

## MLEF with CSU global shallow-water model(Heikes and Randall 1995, MWR; Zupanski et. al. 2005, Tellus)

***Height analysis increment*** [***x****a*-***x****f*]

***Height RMS error*** [***x****a*-***x****t*]

**Initially noisy random perturbations quickly become smooth:**** **Consequence of error covariance localization by dynamics

**Notes:**
Outline slide

## Error covariance localization

**Lyapunov vector**

**Error covariance localization**
**(linear framework)**

**Possible explanation for error covariance localization: **
**Dynamic localization of Lyapunov vectors**

## 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, ********z****max**** = 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

**Notes:**
Outline slide

## LIQUID WATER CONTENT ASSIMILATION

- CONTROL

- EXP

- VERIF

**Better timing of maxima**

- G. Carrio, W. Cotton

**Vertically integrated observations:**
**LWP**

**Vertical structure of the analysis:**
**LWC**

## Independent observation

**ICE FORMING NUCLEI (IFN) CONCENTRATION**

- IFN above the inversion, as observed

- IFN below inversion as cloud forms

- G. Carrio, W. Cotton

- 22Z

## Double-resolution MLEF framework(THORPEX)

**Double-resolution MLEF framework****(THORPEX)**

** *****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:Forecast step

**Double-resolution MLEF:****Forecast step**

- Interpolation operator from the operational (high) to coarse (low) resolution

- Forecast error covariance column vector (low resolution):

*M**C* 	– low resolution forecast model
*M* 	– operational resolution (control) forecast model
***x***C 	– low resolution state vector
***x*** 	– operational resolution state vector

## Double-resolution MLEF:Analysis step

**Double-resolution MLEF:****Analysis step**

- Operational resolution cost-function

- Change of variable (low resolution)

- Directional gradient (ensemble subspace)

- Interpolated from operational resolution

- Low resolution observation operator

**For double-resolution MLEF, need an interpolation operator** *C*,
**low resolution model** *M**C*, **and low resolution observation operator** *H**C*

## 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

**Notes:**
Outline slide

## 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

**Notes:**
Outline slide

## Thank you !

**Notes:**
Outline slide

## ISSUE A: 

***x***min

*J*=*const*

- Starting point 1

- Starting point 2

-  For high Hessian condition numbers ~50-100, minimization success is unpredictable
-  Low-rank assumption in ensemble DA may miss important perturbations (directions)

- ISSUE A:
- Hessian preconditioning vs. low-rank

**Notes:**
Outline slide

## ISSUE B: Error covariances 

***x***1

***x***2

**L**

***Strong*** correlations
- along frontal zone

***Weak*** correlations
- across frontal zone

- ISSUE B: Error covariances
- Dynamical localization vs. prescribed structure

-  Addressed by advanced methods
-  Computational efficiency
-  Moisture related variables (microphysics)

**Notes:**
Outline slide

## 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**

**Notes:**
Outline slide

## Non-Gaussian MLEF framework

**Non-Gaussian MLEF framework**

** ****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)

**Gaussian data assimilation framework is generally used**
**Need to evaluate the impact of this assumption**

## Non-Gaussian MLEF framework

**Non-Gaussian MLEF framework**

** ****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 (2005*a*-SIAM *J.Appl.Math*.; 2005*b*-*J.Roy.Statist.Soc.B*),

** ****Mode, Mean, Median are *****identical***** in Gaussian (or any symmetric) PDF**
** ****Mode, Mean, Median are *****all different***** in Lognormal (or any skewed) PDF**

*mode*

*median*

*mean*

## Non-Gaussian MLEF framework

**Non-Gaussian MLEF framework**

***Lognormal errors are multiplicative:***

**Gaussian PDF**

**Log-Normal PDF**

***Gaussian errors***** *****are additive:***

## Non-Gaussian MLEF experiment with SWM: Lognormal height observation errors

**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: ***

***Higher nonlinearity of the cost function compared to the pure Gaussian****** ***

**Gaussian-like** lognormal observation term [for *ln*(***x***)]

**Gaussian** prior PDF

- Additional** Lognormal **observation term

## Impact of Lognormal observation errors: Analysis RMS errors

**Impact of Lognormal observation errors:**** ****Analysis RMS errors**

*******Logn***= 1.0 ; **********= 3.15 ******** *******Gauss***= 2.5 m**

**Success of the Gaussian MLEF depends on the observation statistics **

*******Logn***= 1.0 ; *********** *****= 1.27** ****** *******Gauss***= 1.5 m**

## Impact of Lognormal observation errors: Innovation histogram

**Impact of Lognormal observation errors:**** ****I****nnovation histogram**

*******Logn***= 1.0 ; **********= 3.15 ******** *******Gauss***= 2.5 m**

**Innovation statistics is significantly impacted by the PDF framework **

*******Logn***= 1.0 ; *********** *****= 1.27** ****** *******Gauss***= 1.5 m**

*N*(0,1)

## MLEF Analysis Step

**Minimize cost function in the subspace *****spanned by ensemble perturbations***** *****p****i**f*

**Use the forecast (prior) error covariance square-root**

**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

**Notes:**
Outline slide

## Information Content AnalysisNASA GEOS-5 Single Column Model

- Zupanski@CIRA.colostate.edu

**Information Content Analysis****NASA GEOS-5 Single Column Model**

***d******s***** measures effective DOF of an ensemble-based data assimilation system (e.g., MLEF). Useful for addressing DOF of the model error.**

**Notes:**
Outline slide

## A : Reduction of uncertainty (0-) 

**Bayesian**

**MLEF**
- 450 ens

- A : Reduction of uncertainty (0-)

- MLEF (full rank) successfully reproduces Bayesian inversion results