Il Sistema di Assimilazione di Il Sistema di Assimilazione di Tipo Ensemble del CNMCA:Tipo Ensemble del CNMCA:

Descrizione e RisultatiDescrizione e Risultati

L.Torrisi, F.Marcucci, A.Troisi, A.Vocino e A.LancianiL.Torrisi, F.Marcucci, A.Troisi, A.Vocino e A.Lanciani 2° Servizio, CNMCA2° Servizio, CNMCA

CNMCA, Pratica di Mare, 08 Luglio 2011

Outline• Data assimilation (DA) in the NWP process• Probabilistic view of DA: Ensemble Data

Assimilation• The Ensemble DA implementation of CNMCA:

Local Ensemble Transform Kalman Filter (LETKF)

• Comparison with CNMCA 3DVar and 4DVar IFS analysis

• Tests with AMSU-A radiances• Future developments

NWP at CNMCA NWP at CNMCA

High resolution models should be initialized with high resolution, accurate initial conditions, both for deterministic and probabilistic forecasting.One of the strong points of the CNMCA NWP system has historically been the attention given to model initialization

14 km40 v.l.

7 km40 v.l.

2.8 km50 v.l.

- HRM hydrostatic model- parameterized convection

- compressible equations- parameterized convection

- compressible equations- explicit convection

CNMCA NWP SYSTEM until 31 May 11 CNMCA NWP SYSTEM until 31 May 11 3DVAR PSAS FGAT (T,u,v,qv,ps) every 3h using TEMP, PILOT, SYNOP, SHIP, BUOY, Wind Profiler, NOAA/METOP AMSUA rad., AMDAR-ACAR-AIREP, MSG/MODIS AMV, METOP/QUIKSCAT/ERS2 scatt. winds + Land SAF snow mask, IFS SST analysis once a day

Local Area Modelling: COSMO

Observation Monitor. and Data Assimilation:

COSMO-ME (7km) ITALIAN MET SERVICE

Post-ProcessingLM- -

DMOCOSMO-MECOSMO-ME

MFS (Oceanic Model)

Field dissemination (rotation / interp.)

Grid point correction (det./statis.)

• FFAA / NATO/ NURC• Civilian Protection Department• RAI3 Regional Forecast (TV)• EUMETSAT-HSAF • SRNWP-PEPS• COSMO verification WG• RAI Teletext• Intranet - prometeo.meteoam.it• Internet - www.meteoam.it• Other users

AWIAWI cloud type, ww, etc

temperature, humidity, cloudiness,wind, pressure, precipitation, etc

Text / Graphics

DMOCOSMO-ITCOSMO-IT

NETTUNONETTUNO3'3'

NETTUNONETTUNO1'1'

Field disseminationGraphics

wind

wind

BC fields

sign.heightmean dir, etc

sign.heightmean dir, etc

prec., etc

Data Assimilation in NWPData Assimilation in NWPThe basic idea of numerical weather prediction is to sample the state of the fluid at a given time (initial condition or analysis) and use the model equations to estimate the state of the fluid at some time in the future. Observations of the initial state of the atmosphere are combined with the results from a numerical model (the fist guess) to produce an analysis, which is considered as 'the best' estimate of the initial state of the system. Based on BLE theory the background and the observations are given a weight proportional to the inverse of their corresponding specified error covariances. As such the accuracy of a data assimilation scheme relies highly on the knowledge of the error statistics of both the model background and the observations

time

FC

OBS

OBS

True state (unknown)

AN

FC

Model Simulation

AN / FC cycle accumulates information from past observations

Data Assimilation in NWPData Assimilation in NWP

AN

Analysis Step

Data Assimilation at CNMCAThe first data assimilation cycle CNMCA was based on an OI scheme, then a 3DVar algorithm was implemented. 1.Can we improve on our currently operational 3DVar without the complications of 4DVar?2. Can we develop a system to consistently evaluate initial and forecast uncertainties?

Ensemble data assimilation has proved to be a viable and competitive alternative to variational methods (Houtekamer, 2005; Fertig et al., 2007; Miyoshi, 2007)Ensemble data assimilation provides “optimal” analysis errors estimates for ensemble forecasting and it overcomes the need for ad hoc inverse methods

Probabilistic view of DA

T=t0

pdf of current obs

T=t1

Assume that a pdf of the state of the atmosphere is available. The first step is to assimilate recent observations, thereby sharpening the pdf. The second step is to propagate the pdf forward in time until new observations are available. If the pdf is initially sharp (i.e., the distribution is relatively narrow), then chaotic dynamics and model uncertainty will usually broaden the probability distribution. The update and forecast steps are then repeated.

Full BayesianFull Bayesianor Extended or Extended Kalman FilterKalman FilterApproachApproach

Prior pdf of state of atmospheregiven past obs

Posterior pdfgiven all obs

Adapted from Miyoshi

Probabilistic view of DA

T=t0AN w/ errors FC w/ errors

P

OBS w/ errorsAN w/ errors

T=t1

Problem ：Degrees of freedom of system and of the error are too large to express explicitly (computationally prohibitive), even assuming normality of error statistics and linearity of error growth (Kalman Filter hypothesis)

Full BayesianFull Bayesianor Extended or Extended Kalman FilterKalman FilterApproachApproach

Ensemble Data AssimilationObs.Analysis Ens. mean

T=t0 T=t1

Generate ensemble members best representing the analysis errors

T=t2

EnDA = ensemble fc. + ensemble update

AN initial condition with errors

FC Ens. mean

P

MontecarloMontecarloApproachApproach

randomly sample the initial pdf

Ensemble Data Assimilation

time BG Forecasts

timeInitial conditions

BG Forecasts

Initial conditions

Analysis Step

Forecast Step

• Ensemble-based assimilation algorithms use Monte-Carlo techniques. As an approximation, the the initial pdf is randomly sampled. Thereafter, the effects of chaos, model error, and observations are simulated.

• With a large enough random sample, probabilities estimated from the ensemble relative frequency will converge to the probabilities that would be calculated explicitly (process computationally feasible and robust).

MontecarloMontecarloApproachApproach

Flow Dependent Errors

Adapted from Schraff (2010)

Advanced DA Methodology• EnDA considers flow-dependent error structures, or the “errors of

the day”. Unlike data assimilation algorithms such as 3D-Var, the background-error covariances can vary in time and space. If this error covariance model is relatively accurate (sampling and model errors have to be taken into account), it will thus provide a better adjustment to the observations. No adjoint operator needed.

Natural linkage between DA and ensemble forecasting• Ensemble forecasts are designed to estimate the flow-dependent

uncertainty of the forecast, while data assimilation techniques require accurate estimates of forecast uncertainty in order to optimally blend the prior forecast(s) with new observations. Ensemble based assimilation methods integrate the two steps analysing the analysis errors in addition to analysis itself. EnDA provides “ideal” ensemble perturbations

Ensemble Data Assimilation

Ensemble Kalman Filter• Ensemble Kalman Filter (EnKF) is an approximation to the KF in

that background-error covariance Pb is estimated from a finite ensemble of forecasts xb.

• Kalman Filter has an update (analysis) step where the state estimate xb and an estimate of the forecast uncertainty Pb are adjusted to new observations y, and a forecast step, where the updated state and the uncertainty estimate are propagated forward to the time when the next set of observations become available.

( )

( ) BGnumber ofm perturbBG xxXariance BGm

e meanBG ensemblxmt ensembleBG forecasxxxx

bbbTbbb

m

j

bj

bbm

bbb

. cov XX1

1P

1x ,...,,1

21

=−=−

=

== ∑=

Analysis step

Forecast step

AnalysisEnsemble

Analysis Error Covariance

Background FCEnsemble Background Error Covar.

( )Tt

bt

b

m11 XX

11

++−

+ q( ) Tqq

mQ

11−

=

H = forward obs operator mappingmodel variables to obs space R = observation err. covariance

M = numerical model

model error covariance (unbiased)

Kalman Gain Matrix

• There are two general classes of EnKF, depending on the kind of update step, stochastic (fully Monte Carlo, perturbed obs methods) and deterministic (square root filters). Both classes propagate the ensemble of analyses with nonlinear forecast models.

• In stochastic EnKF, random noise is applied during the update step to simulate observation uncertainty (ensemble of DA systems).

• In deterministic EnKF, the obs are assimilated to update the ensemble mean and the analysis ensemble perturbations can be formed by transforming the forecast ensemble perturbation through a transform matrix T.

Ensemble Kalman Filter

i=timek=memberAnalysis

Ensemble

AnalysisEnsemble Mean

Analysis Ensemble Perturb.

Analysis Ensemble

At CNMCA the LETKF (Hunt et al. 2007) formulation was chosen (deterministic approach), because algorithmically simple to code, intrinsically parallel, etc. The basic idea is: do analysis in the space of the ensemble perturbations:

– ensemble mean analysis is the linear combination of forecast ensemble states which best fits the observational dataset– analysis ensemble members are locally linear combinationsof background ensemble members– separate analysis are done at every grid point, selecting only obs in a vicinity. This explicit localization combats the problem of spurious correlations between distant locations due to limited sample size provided by an ensemble

aba

a

aba

w Xxx WXX

w Xxx

b

ba

b

+==

+=

]w, ..... ,w[Ww]P~1)-[(mW

)])H(x-)(H(x, ..... ,)H(x-)[(H(xY

]YRY1)I-[(mP~ ))xH(-(yRYP~w

aaaa

aa

bbm

bb1

b

1b-1bTa

b-1bTaa

+==

=

+=

=−

Ensemble Kalman Filter (LETKF)Ensemble Kalman Filter (LETKF)

Analysis Ensemble

Analysis Ensemble MeanAnalysis Ensemble Perturb.

Ensemble Kalman Filter (LETKF)Ensemble Kalman Filter (LETKF)

Diagonal elements are much larger than off-diagonal elements (small corrections)Good forecast members get larger weight in all analysis members

kj

m

jjkk

b

,a

1

bba

aba

w Xxx

kmember afor w Xxx

∑=

+=

+=

wa

Adapted from Schraff (2010)

Model Error• Forecast errors derive not only from errors in the initial conditions but also from errors due to the model deficiencies (called model error). The sources of model error can be due to lack of resolution, approximate parameterizations of physical processes, etc.

• The EnKF tends to underestimate the uncertainty leading to overconfidence in the background state estimate (the analysis underweights the observations). If the discrepancy becomes too large over time, the observations are essentially ignored by the analysis, and the dynamics of the data assimilation system become decoupled from the truth. Generally this tendency is countered by an ad hoc procedure (with at least one tunable parameter) that inflates either the background covariance or the analysis covariance during each data assimilation cycle.

“Multiplicative inflation” instead multiplies the background covariance matrix (or equivalently, the perturbations of the background ensemble members from their mean) by a constant factor larger than one

“Additive inflation” adds random perturbations with a certain

covariance structure to the analysis covariance during each cycle

- State Dependent Multiplicative Inflaction according to Whitaker et al (2010)

- Climatological Additive Noise

Covariance Inflaction

Scale factor

randomly selected, 24h forecast differences

> 0

In the current implementation of LETKF model errors and sampling errors are taken into account using:

40 member ensemble at 0.09° (~10Km) grid spacing (HRM model), 40 hybrid p-sigma vertical levels (top at 10 hPa)

6-hourly assimilation cycle run and (T,u,v,qv,ps) as a set of control variables

Observations: RAOB, SYNOP, SHIP, BUOY, AIREP, AMDAR, ACAR, AMV (MSG, MODIS), WindPROF, SCAT(METOP,ERS2), AMSU-A (METOP,NOAA) radiances (very soon)

Horizontal localization with 800 Km circular local patches (obs weight smoothly decay ∝ r-1)

Vertical localization to layers whose depth increases from 0.2 scale heights at the lowest model levels to 2 scale heights at the model top

Lateral boundary condition perturbation using EPS Climatological perturbed SST Daily blending of the mean upper level analysis with the IFS analysis

to compensate the limited satellite data usage

CNMCA LETKF Implementation CNMCA LETKF Implementation

7 km40 v.l.

2.8 km50 v.l.

- compressible equations- parameterized convection

- compressible equations- explicit convection

CNMCA NWP SYSTEM since 1 June 11CNMCA NWP SYSTEM since 1 June 11LETKF analysis ensemble (40 members) every 6h using TEMP, PILOT, SYNOP, SHIP, BUOY, Wind Profiler, NOAA/METOP AMSUA rad., AMDAR-ACAR-AIREP, MSG AMV, METOP/ERS2 scatt. winds + Land SAF snow mask, IFS SST analysis once a day

Local Area Modelling: COSMO

Ensemble Data Assimilation:

COSMO-ME (7km) ITALIAN MET SERVICE

10 km40 v.l.

- HRM hydrostatic model- parameterized convection

Mean/Control StateAnalysis

OPERATIONAL NWP SYSTEM

Forecast 00

Observations (±3h)Observations (±3h)

IFS B.C.IFS B.C.

LETKF AnalysisEnsemble

00 06 12

Forecast 12

18120600

COSMO-ME (7 km)

00 06 12 18

HRM F.G.Ensemble

18120600

18IFS Analysis

Forecast 00COSMO-IT (2.8 km)

Dat

a A

ssim

ilatio

n Sy

stem

Local Area Modelling

12BlendedMean

18-00Nudging

18-00

18COSMO-MEBC

Obs

Forecast 12 06-12Nudging

06-12

06COSMO-MEBC

Obs

00SN mask

SST

EPSPert.

(since 1 June 2011)

Long Run from LETKF

BG Forecasts

timeInitial conditions BG Forecasts

Analysis Step

BG Forecast Step

Deterministic Long RunfromAnalysis Ensemble Mean

( )

( )

( ) ( ) 2

1

2

1

2

1

1

1

1

1

1

iv

iv

n

ii

ui

u

n

iii

n

iii

n

iii

ii

yxyxn

ErrorSquareMeanRootVectorWind

yxn

ErrorSquareMeanRoot

yxn

ErrorAbsoluteMean

YXyxn

ErrorMean

sizesamplennsobservatioyforecastsx

−+−=

−=

−=

−=−=

===

∑

∑

∑

∑

=

=

=

=

Objective Verification

RMSE measures the magnitude of the error retaining the variable unit. The smaller, the better

ME is the average of errors, it indicates the average error direction (positive bias = over-forecast, negative bias = underforecast)

MAE Average of magnitude of errors (not indicative of the error direction). The smaller, the better

3DVar/LETKF comparison

3DVar/LETKF comparison

3DVar/LETKF comparison

3DVar/LETKF comparison

hits false alarms

misses correct negatives

Precipitation prediction is an example of a yes/no forecast. Some scores (FB, TS, ETS, FAR, etc.) are usually computed from the contingency table.FB=(hits+false alarms)/(hits+misses) The Frequency Bias (FB) measures the ratio of the frequency of forecast events to the frequency of observed events and indicates whether the forecast system has a tendency to underforecast (FB<1) or overforecast (FB>1) events (precipitation frequency).“Good” value: 1TS=hits/(hits+false alarms+misses)The Threat Score (TS) measures the fraction of observed and/or forecast events that were correctly predicted. The TS is somewhat sensitive to the climatology of the event, tending to give poorer scores for rare events. The Equitable Threat Score is designed to help offset this tendency. ETS=(hits-hits expected by chance)/(hits+false alarms+misses–hits expected by chance)where hits expected by chance = (total forecasts of the event) * (total observations of the event) / (sample size)Good value: the larger, the better

Objective Verification 

Contingency Table

forecast

observedyes no

yes

no

3DVar/LETKF comparisonApril 2011

T+6h

T+12h

3DVar/LETKF comparisonMay 2011

T+6h

T+12h

3DVar/LETKF comparisonEuropean stationsApril and May 2011

3DVar/LETKF comparisonEuropean stationsApril and May 2011

3DVar/LETKF comparisonEuropean stationsApril and May 2011

3DVar/LETKF comparisonEuropean stationsApril and May 2011

Deterministic Run from LETKF

DeterministicLong Run fromControl State

BG Forecasts

timeInitial conditions

BG Forecasts

Analysis Step

BG Forecast Step

Control initial condition

DeterministicLong Run fromAnalysisEnsemble Mean

Along with standard LETKF analysis, a control state LETKF analysis is computed. It uses a control run, instead of the BG ensemble mean

)) H(xK (yxx

xHyK

cbc

ac

b

−+=

−+= ))(( xx ba

Control Forecast

time

Control/Mean LETKF comparison

COSMO-ME IC from LETKF Analysis Ensemble Mean

COSMO-ME IC from Control State LETKF

Control/Mean LETKF comparisonEuropean stations

Control/Mean LETKF comparisonEuropean stations

Control/Mean LETKF comparisonEuropean stations

Control/Mean LETKF comparisonEuropean stations

Comparison of LETKF/IFS analysis

The deterministic COSMO-ME runs, initialized by the 00 UTC LETKF analysis ensemble mean and by the global IFS 4DVAR system and driven by IFS boundary conditions, are objectively verified against the European observation network. In this evaluation you have to take into account the reduced number of observation types (no radiances in this experiment) used in the CNMCA LETKF system.

IFS/LETKF analysis comparisonEuropean stations

IFS/LETKF analysis comparisonEuropean stations

IFS/LETKF analysis comparisonEuropean stations

IFS/LETKF analysis comparisonEuropean stations

Operational EnDA • The CanadianMeteorological Centre (CMC) has used the EnKF to provide the initial conditions for the operational ensemble prediction system (EPS) since January 2005 • The Met Office started using an ETKF operationally to initialise the EPS in the summer of 2005.• At European Centre for Medium range Weather Forecast (ECMWF) a hybrid approach is used:

• In June 2010 differences between members of an ensemble of data assimilations (EDA) are used instead of the evolved singular vectors to create initial spread between EPS forecast members

• In May 2011 background error variances are computed from the ensemble of data assimilations (EDA) used by the deterministic 4D-Var assimilation

• The National Center for Environmental Prediction (USA), the Japan Meteorological Agency and the DWD are working on hybrid formulations between EnDA and variational systems. • The CNMCA is the first meteorological centre which uses operationally a pure ensemble data assimilation (LETKF) to initialize a deterministic NWP model (COSMO-ME).

COSMO KENDA

Selecting which satellite radiances to assimilate is complicated by the fact that they not have a single well-defined vertical locationThe weighting function at a particular model point indicates the sensitivity of that observation to the state at that model grid point “MAXIMUM-BASED SELECTIONMAXIMUM-BASED SELECTION”” METHODMETHOD (Fertig et al. 2007) Assign radiance observations to the model level for which the magnitude of the weighting function (wf) is largest. Use this location to select observations within the local region in model space AMSU-A are treated as “single-level” observations Use also the wf shape as vertical covariance localization function

AMSU-A Assimilation

Radiances Treatment

Weighting function (transmittance vert. derivative)

p(hPa)

Radiances Radiances from from RTTOVRTTOV

Radiances Treatment in LETKF

ny

Mean FG time = n

Observation Increments with

Mean FG and QC

Bias CorrectionQuality Control

Observation Increments

ˆ ( )bnH x

a a an n nx x X= +

Mean FG Mean FG Radiances from Radiances from

RTTOVRTTOV Observation Increments with

Mean FG

LETKF Analysis time = n

Radiances Conv. Obs.

FG Ensemble

bnx

ˆ ˆ( ) ( )b bn nH x H x′

AN Ensemble

AMSU-A selection from NOAAxx and MetOp:- Only sea obs from 4-10- Rain check on CH 4- Grody LWP check

AMSU-A Assimilation

Use Ensemble Mean or Cont. as reference for BC and QC

Impact of AMSUA rad.Relative difference (%) in RMSE computed against IFS analysis for 00 UTC EuroHRM runs (28km) from 11-10-2009 to 10-11-2009

Impact of AMSUA rad.Relative difference (%) in RMSE computed against IFS analysis for 00 UTC EuroHRM runs (28km) from 11-10-2009 to 10-11-2009

Operational assimilation of satellite observations

BG radiances are simulated with RTTOV, NWP-SAF’s radiative transfer model, from atmospheric profiles of temperature, variable gas concentrations, cloud and surface properties.

The assimilation of the huge number of brightness temperatures and the tangent linear calculations are computationally expensive.

Thus the following actions are ongoing:

• Optimization and parallelization of the RTTOV / CNMCA DA interface, now fully compatible with COSMO adapters; • Update to the last version of RTTOV (version 10);• New improvements on

1. Bias correction algorithms,2. Other sensors (AMSU-B, MHS…)

Use of IASI Retrievals• The level2 data from hyperspectral IR

sounders are considered crucial to improve the performance of the current and future CNMCA Data Assimilation System (DAS).

• IASI retrievals from EUMETSAT operational dissemination (EUMETCAST level2 product) are currently used in passive monitoring mode to continuously assess the deviations from first-guess analysis fields.

• Dedicated Observing System Simulation Experiments (OSSEs) to evaluate the impact on both 3D-Var and LETKF DAS are on-going.

• Long-lasting experience with IASI represents the starting point for future challenges addressed by the next generation hyperspectral sensors of EUMETSAT programmes.

Conclusions and Future Developments

• LETKF, pre-operational from the end of December 2010 at CNMCA, has substitute 3D-Var for COSMO-ME the 1 June 2011• The use of AMSU-A should furtherly improve the LETKF DA. • Assimilation of IASI retrievals, AMSU-B/MHS will be investigated.• Non-linearities and balancing are issues to address• Tests with COSMO model and shorter assimilation window• Further tuning of model error representation (tuning of covariance localization, evolved additive noise, bias correction, etc.)• Implement a Short-Range EPS based on LETKF

Grazie per l’attenzione!