Data Assimilation in Meteorology Chris Budd Joint work with Chiara Piccolo, Mike Cullen (Met Office) Melina Freitag, Phil Browne, Emily Walsh, Nathan Smith and Sian.
Download ReportTranscript Data Assimilation in Meteorology Chris Budd Joint work with Chiara Piccolo, Mike Cullen (Met Office) Melina Freitag, Phil Browne, Emily Walsh, Nathan Smith and Sian.
Data Assimilation in Meteorology Chris Budd Joint work with Chiara Piccolo, Mike Cullen (Met Office) Melina Freitag, Phil Browne, Emily Walsh, Nathan Smith and Sian Jenkins (Bath) Understanding and forecasting the weather is essential to the future of planet earth and maths place a central role in doing this Accurate weather forecasting is a mixture of • Careful modelling of the complex physics of the ocean and atmosphere • Accurate computations on these models • Systematic collection of data • A fusion of data and computation Data assimilation is the optimal way of combining a complex model with uncertain data Integrated forecasting process Observations from space Upper-air observations Surface observations Weather radar Analysis Intervention NWP forecasts Fine-tuning Forecast products and guidance Modelling the global atmosphere and ocean Vertical exchange between layers of momentum, heat and moisture 15° W 60° N Horizontal exchange between columns of momentum, heat and moisture 3.75° 2.5° Vertical exchange between layers of momentum, heat and salts by diffusion, convection and upwelling 11.25° E 47.5° N Vertical exchange between layers by diffusion and advection Orography, vegetation and surface characteristics included at surface on each grid box Met Office Current Configurations UK North Atlantic and Europe (NAE) • Global 25km L70 model (was 40km L50) • Incremental 4D-Var Data Assimilation • 60km 24m ETKF ensemble (was 90km L38) • Regional NAE 12km L70 model • Incremental 4D-Var DA • 16km 24m L70 ETKF ensemble (was 24km L38) • UK 1.5km model (stretched) (was 4km) • Incremental 3D-Var DA Met. Office Global/Regional Ensemble Prediction System (MOGREPS) became fully operational in Sep 2008 after 3 years of trials. Focus on short-range out to 72hr. Data: Sources of observation Observation Volumes in 6 hours Category TEMPs PILOTs Wind Profiler Land Synops Ships Buoys Amdars Aireps GPS-RO Count % Category used 637 99% Satwinds: JMA 307 99% Satwinds: NESDIS 1355 39% Satwinds: EUMETSAT 16551 99% Scatwinds: Seawinds 3034 84% Scatwinds: ERS 8727 63% Scatwinds: ASCAT Count % use d 26103 4% 142478 3% 220957 1% 436566 1% 27075 2% 241626 4% 64147 23% SSMI/S 7144 12% SSMI 532140 1% 776 99% ATOVS 1127224 3% AIRS IASI 698048 1% 75824 6% 80280 3% Performance Improvements “Improved by about a day per decade” Met Office RMS surface pressure error over the N. Atlantic & W. Europe DA Introduced Andrew Lorenc What are the causes of improvements to NWP systems? 1. Model improvements, especially resolution and sub grid modelling 2. Better observations 3. Careful use of forecast & observations, allowing for their information content and errors. Achieved by variational data assimilation e.g. of satellite radiances. Basic Idea of Data Assimilation True state of the weather is xt Numerical Weather Prediction NWP calculation x gives a predicted state b with an error Make a series of observations y of some function H(x t ) of the true state Eg. Limited set of temperature measurements with error Now combine the prediction with the observations Both the NWP prediction and the data have errors. Can we optimally estimate the atmospheric state which is consistent with both the prediction and the data and estimate the resulting error? NOTE: Approximately 10^9 degrees of freedom 10^6 data points So significantly underdetermined problem Data y Best state estimate analysis xf xa NWP prediction xB xB Assume initially: Use x a to produce x forecast f 6 hours later variables 1. Errors are unbiased Gaussian 2. Data and NWP prediction errors are uncorrelated 3. H(x) is a linear operator Can estimate x a using Bayesian analyis: Maximum likelihood estimate of data y given Posterior xt P(x t y) P(y) P(y x t ) P(x t ) Prior Likelihood Best RMS unbiased estimate of the true state: BLUE Minimum error variance Assumptions about the error xB Data error: Gaussian, Covariance R Background (NWP) error: Gaussian, Covariance B BLUE: Find x a which minimises 1 1 T 1 J(xa ) xa xb B xa xb (Hxa y)T R1(Hxa y) 2 2 NOTE: P(x y)/P(x) eJ(x) If R and B are known then the best estimate of the analysis is T 1 x a x b K(y Hxb ), K BH (R HBH ) T Covariance of the analysis error A KRK (I KH)B(I KH) T Kalman filter: Continuously updates the forecast and its error given the incoming data. T Implementation: In the context of minimising the functional 1 1 T 1 T 1 J(xa ) xa xb B xa xb (Hxa y) R (Hxa y) 2 2 This is implemented as 3D-VAR (since 1999 in the Met Office) xB xa xf : Background, derived from 6 hour NWP forecast : Analysis : NWP forecast using xa as initial data Ensemble Kalman Filter EnKF This is a widely used Monte Carlo method that uses an ensemble of forecasts to estimate the terms in the Kalman filter Idea: Take a large number of initial states x i and estimate the resulting background states x B ,i xi xi xi xi xi Estimate 1 x xb,i,, N x B ,i x B ,i x B ,i x B ,i x B ,i 1 T B xB,i x xB,i x N 1 Basic Filtering Idea xB xa xa xB y Advantages: Works well with high dimensional systems Disadvantage: EnKF inaccurate with strong nonlinearity eg. Shear flow [Jones] 4D VAR … Preferred variational method Use window of several observations over 6 hours Obs. Jo x Previous forecast Jo xb Obs. Jb xa 9 UTC Jo Jo Obs. Corrected forecast Obs. 12 UTC Assimilation Window 15 UTC Time 4D-VAR idea: Evolutionary model M (nonlinear) Unknown initial state Times x0 t t 0 ,t1,t 2 , Leads to state estimates Data yi over window Find x 0 so that the estimates fit the data Smoothing Over a time window x1, x 2 ,... Minimise N 1 1 T 1 J(x 0 ) x 0 x b B x 0 x b (Hxi y i )T R1(Hxi y i ) 2 2 i 0 Subject to the strong model constraint x i1 M i (x i ) At present assume perfect model, but can also deal with certain types of model error (both random and systematic) by using a weak constraint instead Usually solved by introducing Lagrange multipliers 1 1 N T 1 J(x a ) x 0 x b B x 0 x b (Hxi y i )T R1(Hxi y i ) 2 2 i 0 N 1 i (x i1 M x i ) i 0 And solving the adjoint problems 0 Ji H T R1(Hxi y i ) i1 i mi x i 0 J0 B (x 0 x b ) H R (Hx0 y 0 ) 0 m0 x 0 1 T 1 Solution: Find x 0 to minimise nonlinear function J Need forward calculation to find x i and backward solve to find i expensive for high dimensional problems!!! Only VERY have limited time to dothe calculation (20 mins) Incremental 4D-Var: Cheaper! 1. Assume x 0 is close to x B 2. Linearise J about x B and minimise this function using an iterative method eg. BFGS (not usually) 3.Repeat if needed BUT: Relies on assumption of near linearity to work well Very effective method!! Met Office operational in 2004 [Lorenc, …. ] Used by many other centres Estimation of the background and covariance errors Good estimates of the covariance matrices R and B are important to the effectiveness of 4D-VAR 1. To get the physics correct 2. To avoid spurious correlations between parameters 3. To give well conditioned systems NOTE: B is a very large matrix, difficult to store and very difficult to update. Impractical to calculate using the Fokker-Plank equation R Different instrument error characteristics and errors of representativeness B Enormous: 10^8 x 10^8 Deduce structure from: 1 5 6 2 9 4 8 0 3 Historical data Known dynamical and physical structure of the atmosphere eg. Balance relationships [Bannister] Build meteorology into the calculation of B through Control Variable Transformations (CVTs) IDEA: Choose more ‘natural’ physical variables which have uncorrelated errors so that the transformed covariance matrix is block diagonal or even the identity Set x U U pUvUh , B UUT Reduces the complexity of the system AND gives better conditioning for the linear systems Up Uv Uh 1 1 1 Separates physical parameters into uncorrelated ones eg. temperature, wind, balanced and unbalanced Reduces vertical correlations by projecting onto empirical orthogonal vertical modes Reduces horizontal correlations by projecting onto spherical harmonics Effective, but errors arise due to lack of resolution of physical features leading to spurious correlations. Eg. Problems with stable boundary and inversion layers and assimilating radiosonde data Poor resolution leads to inaccurate predictions of fog and ice Solution one: increase global resolution VERY EXPENSIVE!!! Solution two: locally redistribute the computational mesh to resolve the features Cheap and effective! [Piccolo, Cullen, B,Browne, Walsh] Adjust the vertical coordinates to concentrate points close to the inversion layer and reduce correlations Introduce an extra transformation x U U pUaUvUh , U 1 a B UU T Adaptive mesh transformation applied to latitude-longitude coordinates Do this by using tools from adaptive mesh generation methods for PDES Set: z original height variable new ‘computational’ height variable Relate these via the equation dz M(z) d M called the ‘monitor function’ [B, Huang, Russell, Walsh] Take M large if there is active meteorology Eg. High potential vorticity 2 M 1 c z 2 Initially use background state estimate, then update Monitor function and the Adaptive Grid Piccolo&Cullen QJR Met Soc 2011 © Crown copyright Met Office First calculation UK4 domain: 3 Jan 2011 00z Updated calculation © Crown copyright Met Office Applied to the Met Office UK4 model Test case: 8th Feb 2010. Significant reduction in RMS error especially for temperatures Piccolo&Cullen, QJR Met Soc 2011 RMS T (K) RH (%) u (m/ s) v (m/s) Control 0.76 0.045 1.32 1.16 Test 0.64 0.045 1.29 1.16 Nobs 1011 901 819 819 Particularly effective for the 2m temperatures Used together with Met Office Open Road software to advise councils on road gritting over Christmas Adaptive mesh implemented operationally in November 2010. Now extending it to a fully three dimensional implementation [B,Browne,Piccolo] Other refinements to 4D-Var Change in the background norm [Freitag, B, Nichols] N 1 1 J(x 0 ) TV(x 0 x b ) (Hxi y i )T R1 (Hxi y i ) 2 2 i 0 Total variation: Gives significantly better resolution of fronts, shocks and other localised features But .. Hard to implement in high dimensions!! Dealing with model error If model has random errors with Covariance C can extend 4D-Var to find the minimiser x 0 , x1, x 2 , of N 1 1 T 1 J(x 0 ) x 0 x b B x 0 x b (Hxi y i )T R1(Hxi y i ) 2 2 i 0 1 T 1 x M (x ) C x i1 M i (x i ) i1 i i 2 However, most model errors eg. Numerical errors are systematic. Dealing with these and quantifying the uncertainty in DA is an area of active research [Jenkins, B, Smith, Freitag], [Cullen&Piccolo], [Stuart] Dealing with nonlinearity Lot of research into finding a compromise between dealing with the high dimensionality and nonlinearity in the system Better use of appropriate (eg. Lagrangian) data Tuning method to data [Jones, Stuart, Apte] Use of particle filters and MCMC methods [Peter Van Leeuwan] Conclusions Data assimilation is an optimal way of merging models with data Useful for model tuning, validation, evaluation, uncertainty quantification and reduction Very effective in meteorology Many other applications to Planet Earth eg. Climate change, oil reservoir modelling, geophysics, energy management and even crowd behaviour Problems with the simple Kalman Filter • Assumption of Gaussian random variables • Assumption of linearity • Assumption of known covariances • Covariance matrix B is VERY large for meteorological problems • Minimisation is a large complex problem Problems with 4D Var Accuracy • Estimation of the background covariance matrix • Ill conditioning of the linear systems • Reliance of near linearity • Inappropriate use of background data • Incorrect covariance between data • Unresolved random and systematic model error • Poor resolution in the model Improvements subject to much research