620 likes | 819 Views
Examples of Four-Dimensional Data Assimilation in Oceanography. Ibrahim Hoteit. University of Maryland October 3, 2007. Outline. 4D Data Assimilation 4D-VAR and Kalman Filtering Application to Oceanography Examples in Oceanography 4D-VAR Assimilation Tropical Pacific, San Diego, …
E N D
Examples of Four-Dimensional Data Assimilation in Oceanography Ibrahim Hoteit University of Maryland October 3, 2007
Outline • 4D Data Assimilation • 4D-VAR and Kalman Filtering • Application to Oceanography • Examples in Oceanography • 4D-VAR Assimilation Tropical Pacific, San Diego, … • Filtering Methods Mediterranean Sea, Coupled models, Nonlinear filtering ... • Discussion and New Applications
Data Assimilation • Goal: Estimate the state of a dynamical system • Information: • Imperfect dynamical Model: • state vector, model error • transition operator form to • Sparse observations: • observation vector, observation error • observational operator • A priori Knowledge: and its uncertainties
Data Assimilation • Data assimilation: Use all available information to determine the best possible estimate of the system state • Observations show the real trajectory to the model • Model dynamically interpolates the observations • 3D assimilation:Determine an estimate of the state at a given time given an observation by minimizing • 4D assimilation:Determine given 4D-VAR and Kalman Filtering
4D-VAR Approach • Optimal Control:Look for the model trajectory that best fits the observations by adjusting a set of “control variables” minimize with the model as constraint: • is the control vector and may include any model parameter (IC, OB, bulk coefficients, etc) … and model errors • Use a gradient descent algorithm to minimize • Most efficient way to compute the gradients is to run the adjoint model backward in time
Kalman Filtering Approach • Bayesian estimation: Determine pdf of given • Minimum Variance (MV) estimate (minimum error on average) • Maximum a posteriori (MAP) estimate (most likely) • Kalman filter (KF): provides the MV (and MAP) estimate for linear-Gaussian Systems
The Kalman Filter (KF) Algorithm • Initialization Step: Analysis Step (observation) Forecast Step (model) • Kalman Gain • Analysis state • Analysis Error covariance • Forecast state • Forecast Error covariance
Application to Oceanography • 4D-VAR and the Kalman filter lead to the same estimate at the end of the assimilation window when the system is linear, Gaussian and perfect • Nonlinear system: • 4D-VAR cost function is non-convex multiple minima • Linearize the system suboptimal Extended KF (EKF) • System dimension ~ 108: • 4D-VAR control vector is huge • KF error covariance matrices are prohibitive • Errors statistics: • Poorly known • Non-Gaussian: KF is still the MV among linear estimators
4D Variational Assimilation • ECCO 1o Global Assimilation System • Eddy-Permitting 4D-VAR Assimilation • ECCO Assimilation Efforts at SIO • Tropical Pacific, San Diego, … In collaboration with the ECCO group, especially Armin Köhl*, Detlef Stammer*, Patrick Heimbach** *Universitat Hamburg/Germany, **MIT/USA
ECCO 1o Global Assimilation System • MITGCM (TAF-compiler enabled) • NCEP forcing and Levitus initial conditions • Model: • Data: • Assimilation scheme:4D-VAR with control of the initial conditions and the atmospheric forcing (with diagonal weights!!!) • ECCO reanalysis:1o global ocean state and atmospheric forcing from 1992 to 2004, …and from 1952 2001 (Stammer et al. …) • Altimetry (daily):SLA TOPEX, ERS • SST (monthly): TMI and Reynolds • Profiles (monthly) : XBTs, TAO, Drifters, SSS, ... • Climatology (Levitus S/T) and Geoid (Grace mission)
ECCO Solution Fit Equatorial Under Current (EUC) Johnson ECCO
ECCO Tropical Pacific Configuration • Regional: 26oS 26oN, 1/3o, 50 layers, ECCO O.B. • Data: • TOPEX, TMI SST, TAO, XBT, CTD, ARGO, Drifters; all at roughly daily frequency • Climatology: Levitus-T and S, Reynolds SST and GRACE • Control: Initial conditions , 2-daily forcing, and weekly O.B. • Smoothing: Smooth ctrl fields using Laplacian in the horizontal and first derivatives in the vertical and in time • First guess: Levitus (I.C.), NCEP (forcing), ECCO (O.B.) MITGCM Tropical Pacific OB = (U,V,S,T)
Eddy-Permitting 4D-VAR Assimilation • The variables of the adjoint model exponentially increase in time • Typical behavior for the adjoint of a nonlinear chaotic model • Indicate unpredictable events and multiple local minima • Correct gradients but wrong sensitivities • Invalidate the use of a gradient-based optimization algorithm • Assimilate over short periods (2 months) where the adjoint is stable • Replace the original unstable adjoint with the adjoint of a tangent linear model which has been modified to be stable (Köhl et al., Tellus-2002) • Exponentially increasing gradients were filtered out using larger viscosity and diffusivity terms in the adjoint model
HFL gradients after 45 days with increasing viscosities Visc = 1e11 & Diff = 4e2 10*Visc & 10*Diff 20*Visc & 20*Diff 30*Visc & 30*Diff
Initial temperature gradients after 1 year (2000) 10*Visc & 10*Diff 20*Visc & 20*Diff
Data Cost Function Terms 1/6;39 1/3;39 1;39 1;23
Control Cost Function Terms 1/6;39 1/3;39 1;39 1;23
Fit to Data 1/6;39 1/3;39 1;39 1;23
What Next … • Fit is quite good and assimilation solution is reasonable • Extend assimilation period over several years • Add new controls to enhance the controllability of the system and reduce errors in the controls • Improve control constraints … • Some references • Hoteit et al. (QJRMS-2006) • Hoteit et al. (JAOT-2007) • Hoteit et al. (???-2007)
Other MITGCM Assimilation Efforts at SIO • 1/10o CalCOFI 4D-VAR assimilation system • Predicting the loop current in the Gulf of Mexico … • San Diego high frequency CODAR assimilation • Assimilate hourly HF radar data and other data • Adjoint effectiveness at small scale • Information content of surface velocity data • MITGCM with 1km resolution and 40 layers • Control: I.C., hourly forcing and O.B. • First guess: one profile T, S and TAU (no U, V, S/H-FLUX) • Preliminary results: 1 week, no tides
Model Domain and Radars Coverage Time evolution of the normalized radar cost 1/6;39 1/3;39 1;39 1;23
Assimilation Solution: SSH / (U,V) & Wind Adj. 1/6;39 1/3;39 1;39 1;23
What Next … • Assimilation over longer periods • Include tidal forcing • Coupling with atmospheric model • Nesting into the CalCOFI model
Filtering Methods • Low-Rank Extended/Ensemble Kalman Filtering • SEEK/SEIK Filters • Application to Mediterranean Sea • Kalman Filtering for Coupled Models • Particle Kalman Filtering In collaboration with D.-T. Pham*, G. Triantafyllou**, G. Korres** *CNRS/France, **HCMR/Greece
Low-rank Extended/Ensemble Kalman Filtering • Reduced-order Kalman filters:Project on a low-dim subspace • Kalman correction along the directions of • Reduced error subspace Kalman filters: has low-rank • Ensemble Kalman filters:Monte Carlo approach to • Correction along the directions of
Singular Evolutive Kalman (SEEK) Filters • Low-rank (r) error subspace Kalman filters: Analysis Forecast • A “collection” of SEEK filters: • SEEK:Extended variant • SFEK: Fixed variant • SEIK:Ensemble variant with (r+1) members only! (~ETKF) Inflation and Localization
The Work Package WP12 in MFSTEP • EU project between several European institutes • Assimilate physical & biological observations into coupled ecosystem models of the Mediterranean Sea: • Develop coupled physical-biological model for regional and coastal areas of the Mediterranean Sea • Implement Kalman filtering techniques with the physical and biological model • … Investigate the capacity of surface observations (SSH, CHL) to improve the behavior of the coupled system
The Coupled POM – BFM Model One way coupled: Ecology does not affect the physics
A Model Snapshot 1/10o Eastern Mediterranean configuration 25 layers Elevation and Mean Velocity Mean CHL integrated 1-120m
Assimilation into POM • Model = 1/10o Mediterranean configuration with 25 layers • Observations = Altimetry, SST, Profiles T & S profiles, Argo data, and XBTs on a weekly basis • SEIK Filterwith rank 50 (51 members) • Initialization = EOFs computed from 3-days outputs of a 3-year model integration • Inflation factor = 0.5 • Localization = 400 Km
Mean Free-run RMS Error Assimilation into POM SSH RMS Misfits Mean Forecast RMS Error Forecast Free-Run Obs Error = 3cm Mean Analysis RMS Error Analysis
Salinity RMS Error Time SeriesFerryBox data at Rhone River SSH 07/12/2005 07/12/05: SATELLITE SSH FREE RUN FORECAST ANALYSIS 20/8/2014 ECOOP KICK-OFF
Assimilation into BFM • Model = 1/10o Eastern Mediterranean with 25 layers with perfect physics • Observations = SeaWiFS CHL every 8 days in 1999 • SFEK Filter = SEEK with invariant correction subspace • Correction subspace = 25 EOFs computed from 2-days outputs of a one year model integration • Inflation factor = 0.3 • Localization = 200 Km
Assimilation into BFM CHL RMS Misfits Free-Run Forecast Analysis
Ph Cross-Section at 28oE CHL Cross-Section at 34oN
Physical System • Ecological System Kalman Filtering for Coupled Models • MAP: Direct maximization of the joint conditional density standard Kalman filter estimation problem • Joint approach: strong coupling and same filter (rank) !!!
Decompose the joint density into marginal densities Dual Approach • Compute MAP estimators from each marginal density • Separate optimization leads to two Kalman filters … • Different degrees of simplification and ranks for each filter significant cost reduction • Same from the joint or the dual approach • The physical filter assimilate and
RRMS for state vectors • Twin-Experiments 1/10o Eastern Mediterranean (25 layers) • Joint:SEEK rank-50 • Dual:SEEK rank-50 for physics SFEK rank-20 for biology • . Physics Biology Ref Dual Joint Ref Dual Joint
What Next … • Joint/Dual Kalman filtering with real data • State/Parameter Kalman estimation • Better account for model errors • Some references • Hoteit et al. (JMS-2003), Triantafyllou et al. (JMS-2004), Hoteit et al. (NPG-2005), Hoteit et al. (AG-2005), (Hoteit et al., 2006), Korres et al. (OS-2007)
Nonlinear Filtering - Motivations • The EnKF is “semi-optimal”; it is analysis step is linear • The optimal solution can be obtained from the optimal nonlinear filter which provides the state pdf given previous data • Particle filter (PF) approximates the state pdf by mixture of Dirac functions but suffers from the collapse (degeneracy) of its particles (analysis step only update the weights ) • Surprisingly, recent results suggest that the EnKF is more stable than the PF for small ensembles because the Kalman correction attenuates the collapse of the ensemble
The Particle Kalman Filter (PKF) • The PKF uses a Kernel estimator to approximate the pdfs of the nonlinear filter by a mixture of Gaussian densities • The state pdfs can be always approximated by mixture of Gaussian densities of the same form: • Analysis Step: Kalman-type: EKF analysis to update and Particle-type: weight update (but using instead of ) • Forecast Step: EKF forecast step to propagate and • Resampling Step: …
Particle Kalman Filtering in Oceanography • It is an ensemble of extended Kalman filters with weights!! Particle Kalman Filtering • requires simplification of the particles error covariance matrices • The EnKF can be derived as a simplified PKF • Hoteit et al. (MWR-2007) successfully tested one low-rank PKF with twin experiments • What Next … • Derive and test several simplified variants of the PKF • Assess the relevance of a nonlinear analysis step: comparison with the EnKF • Assimilation of real data …
Discussion and New Applications • Advanced 4D data assimilation methods can be now applied to complex oceanic and atmospheric problems • More work is still needed for the estimation of the error covariance matrices, the assimilation into coupled models, and the implementation of the optimal nonlinear filter • New Applications: • ENSO prediction using neural models and Kalman filters • Hurricane reconstruction using 4D-VAR ocean assimilation! • Ensemble sensitivities and 4D-VAR • Optimization of Gliders trajectories in the Gulf of Mexico • … THANK YOU
4D-VAR or (Ensemble) Kalman Filter? 4D-VAR EnKF 4D-VAR or EnkF? …
Comparison with TAO-Array RMS RMS Zonal Velocity (m/s) 1/6;39 1/3;39 RMS Meridional Velocity (m/s) 1;39 1;23
San Diego HF Radar Currents Assimilation • Assimilate hourly HF radar data and other data • Goals: • Adjoint effectiveness at small scale • Information content of surface velocity data • Dispersion of larvae, nutrients, and pollutants • MITGCM with 1km resolution with 40 layers • Control: I.C., hourly forcing, and O.B. • First guess: one profile, no U and V, and no forcing! • Preliminary results: 1 week, no tides
Cost Function terms 1/6;39 1/3;39 1;39 1;23