treatment in geophysical data assimilation - Some ideas - - PowerPoint PPT Presentation

Model error ✭✭✭✭✭✭

✭ ❤❤❤❤❤❤ ❤

treatment in geophysical data assimilation - Some ideas

Alberto Carrassi

Nansen Environmental and Remote Sensing Center, Norway Geophysical Institute, University of Bergen, Norway

With:

• P. Ailliot (Un. Brest, FR), M. Bocquet (ENPC-CEREA, FR), M. Lucini (Un. Reading, UK), L. Mitchell (Un. Adelaide, AUS),
• T. Miyoshi (RIKEN, JP), M. Pulido (Un. Reading, UK), P. Raanes (NORCE, NO), P. Tandeo (IMT, FR), S. Vannitsem

(RMI, BE), Y. Zhen (IMT, FR). 7th International Symposium on Data Assimilation - ISDA 2019 Carrassi et al. Model error in DA - ISDA 23rd January 2019 1 / 21

DA and model error The impact of model error

The impact of model error

◮ For years model error impacts on NWP predictions was considered small compared to the (growth of) i.c. error, and thus often neglected in DA. ◮ The amelioration of the i.c. & the increase of the forecast horizons (seasonal-to-interannual) led to a larger impact of the model error on prediction skill. ◮ In DA it often manifests as underestimation of the estimate state error co-variance ⇒ Inflation. ◮ Particularly on long timescales, model error becomes evident through the emergence of biases.

ECMWF IFS model coupled with NEMO ocean model. Sea surface forecast bias (Years 14–23). Figure from Magnusson et al., 2012

Carrassi et al. Model error in DA - ISDA 23rd January 2019 2 / 21

DA and model error The posing of the problem

Posing of the problem: Nonlinear Gaussian state-space model

It is usually assumed an HMM such as: xk = Mk:k−1(xk−1, λ) + ηk, yk = Hk(xk) + ǫk. (1) ◮ xk ∈ Rm and λ ∈ Rp are the model state and parameter vectors respectively. ◮ yk ∈ Rd are noisy observations related to the system’s state via the, generally nonlinear,

• bservation operator, H : Rm → Rd

◮ Mk:k−1 : Rm → Rm is usually a nonlinear, possibly chaotic, function from time tk−1 to tk. ◮ The model and the observational errors, ηk and ǫk, are usually assumed to be uncorrelated in time, mutually independent, and Gaussian distributed: ηk ∼ N(0, Qk) and ǫk ∼ N(0, Rk) Given the multiple sources of model error a stochastic approach is generally used. An accurate estimate of the model error covariance, Qk, is necessary.

Carrassi et al. Model error in DA - ISDA 23rd January 2019 3 / 21

DA and model error 1D illustration

The importance of a good Q - 1D illustration

Perfect Q Under-estimated Q Over-estimated Q

Univariate, linear case. xk =0.95xk−1 + ηk (2) yk =xk + ǫk (3) with ηk ∼ N(0, Qt) and ǫk ∼ N(0, Rt) ◮ The over-estimation is “safer” (i.e. “erring on the side of caution”, cf Raanes et al, 2019). ◮ Promote the use of inflation.

Tandeo et al, 2019 - Under review Carrassi et al. Model error in DA - ISDA 23rd January 2019 4 / 21

DA and model error 1D illustration

The importance of a good ||Q/R|| ratio - 1D illustration

◮ It is the ratio Q/R that matters for the accuracy of the state estimate. ◮ Good Q/R (no matter the individual estimates of Q and R) suffices to get good RMSE ◮ However it impacts differently the uncertainty quantification (i.e. coverage probability).

Carrassi et al. Model error in DA - ISDA 23rd January 2019 5 / 21

DA and model error 1D illustration

The importance of simultaneously estimating Q and R - 1D illustration

◮ Estimate Q or R with the Expectation Maximization (EM) (Shumway and Stoffer, 1982) ◮ Figure from Tandeo et al, 2019 - Under Review

It is not possible to fully compensate for the misrepresentation of Q/R by optimizing R/Q ⇒ The best is to estimate Q and R simultaneously.

Carrassi et al. Model error in DA - ISDA 23rd January 2019 6 / 21

DA and model error Estimating Q: key obstacles and objectives

Estimating Q: key obstacles and objectives

Large variety of possible error sources (incorrect parametrizations of physical processes, numerical discretizations, unresolved scales, etc..) The amount of available data insufficient to realistically describe the model error statistics, i.e. dim(y) = d ≪ dim(x) = m. Lack of a general framework for model error dynamics (as opposed to the dynamics of the i.c. error). Two separate but related objectives:

1 Is the white-noise assumption always a good one? 2 Can we efficiently estimate Qk along with the system state? Carrassi et al. Model error in DA - ISDA 23rd January 2019 7 / 21

DA and model error Time-correlated model error

Time-correlated model error - Formulation

Let assume to have the model: dx(t) dt = f(x, λ) used to describe the true process: dˆ x(t) dt = ˆ f(ˆ x, ˆ y, λ

′)

dˆ y(t) dt = ˆ h(ˆ x, ˆ y, λ

′)

◮ ˆ h(ˆ x, ˆ y, λ

′): unresolved scale; ∆λ = λ ′ − λ parametric error.

The evolution of the error covariance in the resolved scale: P(t) =< δx0δxT

0 > +

t

t0

dτ t

t0

dτ

′ < [f(x, λ) − ˆ

f(ˆ x, ˆ y, λ′)][f(x, λ) − ˆ f(ˆ x, ˆ y, λ′)] >T (4) ◮ The important factor controlling the evolution is the difference between the velocity fields, the tendencies f(x, λ) − ˆ f(ˆ x, ˆ y, λ)

Carrassi et al. Model error in DA - ISDA 23rd January 2019 8 / 21

DA and model error Time-correlated model error

Example: parametric error

◮ Assume the model resolves all scales ⇒ ˆ h = 0 and f = ˆ f. ◮ But it has error in the parameter δλ = 0. The (linearized) error evolution reads: δx(t) = x(t) − ˆ x(t) ≈ Mt:t0δx0 + t

t0

dτMt:τδµ(τ), δµ = ∂f ∂λ|λδλ (5) ◮ The model error acts as a deterministic process and the factor controlling its evolution is δµ(t) ◮ In view of the presence of the TLM propagator M, the flow instabilities are expected to influence the model error dynamics The model error covariance evolution reads: Q(t1, t2) = t1

t0

dτ t2

t0

dτ

′Mt1,τ < δµ(τ)δµ(τ ′)T > MT

t2,τ ′

(6)

Carrassi et al. Model error in DA - ISDA 23rd January 2019 9 / 21

DA and model error Time-correlated model error

Time-correlated model error - Formulation

◮ The evolution equation for the model error covariance cannot be implemented in high dimension. ◮ A suitable approximation can be obtained for short-time (e.g. the assimilation window). Q(t1, t2) [f(x, λ) − ˆ f(ˆ x, ˆ y, λ

′)][f(x, λ) − ˆ

f(ˆ x, ˆ y, λ

′)]T(t1 − t2)2 + O(3)

(7) ◮ The difference between the model and the nature tendencies, f(x, λ) − ˆ f(ˆ x, ˆ y, λ

′) is treated as

being correlated in time. ◮ The white-noise case would correspond to the terms f(x, λ) − ˆ f(ˆ x, ˆ y, λ

′) being delta-correlated

and the short-time evolution would be bound to be linear.

Carrassi et al. Model error in DA - ISDA 23rd January 2019 10 / 21

DA and model error Time-correlated model error

How to estimate the model-to-nature tendencies difference

Making use of the reanalysis ⇒ Qt ≈< (f − ˆ f)(f − ˆ f)T > t2 ◮ Needs to estimate the statistics of the velocity fields discrepancy. ◮ Use of the analysis increments from a reanalysis data-set assumed to be the “truth”: f − ˆ f = dx dt − dˆ x dt ≈ xf

r(t + τr) − xa r(t)

τr − xa

r(t + τr) − xa r(t)

τr = δxa

r

τr ⇒ Q(t) ≈< δxa

rδxa r T > τ 2

τ 2

r

with τr reanalysis assimilation interval and τ current assimilation interval.

Carrassi et al. Model error in DA - ISDA 23rd January 2019 11 / 21

DA and model error Time-correlated model error

Short-time EKF with re-analysis based Q. Comparison with inflated EKF

◮ L96 two scales. Neglect the fast scales in the model and observe 12/36 points on the coarse scale. ◮ (a): EKF with inflation Pf → (1 + ρ)Pf ◮ (b): Short-time EKF ST-EKF, with tuned re-analysis based Q → αQ ◮ (c): Analysis error for optimally tuned ST-EKF (α = 0.5 red line) and EKF (ρ = 0.09 black line)

0.1 0.3 0.5 0.05 0.1 0.15 0.2 0.25 ρ Normalized Error Variance 10

−4

10 10

4

0.05 0.1 0.15 0.2 0.25 α 60 120 180 0.1 0.2 0.3 0.4 0.5 Day (a) (b) (c) Carrassi and Vannitsem, 2011 Carrassi et al. Model error in DA - ISDA 23rd January 2019 12 / 21

DA and model error Time-correlated model error

EnKF with short-time correlated model error

◮ L96 two scales. ETKF (Bishop et al, 2001) with “best tuned” multiplicative inflation and localization (red line). ◮ ETKF with model error matrix Q estimated using the short-time approximation and the re-analysis (ETKF-TC, green line). ◮ ETKF with time-varying model error, randomly sampled from the reanalysis-increment statistics (ETKF-TV blue line) such that xf

i = M(xa i ) + ηi τ τr

ηk ∼ N( ¯ δxa

r, Q)

i = 1, ..., N

Mitchell and Carrassi, 2015 Carrassi et al. Model error in DA - ISDA 23rd January 2019 13 / 21

DA and model error Time-correlated model error

4DVar with short-time correlated model error

Minimize the cost-function: 2J = τ τ (δxt1)TQ−1

t1t2(δxt2)dt1dt2+...

Model Lorenz 3-variables. Strong-constraint - Assume perfect model. Weak constraint 4DVar with uncorrelated model error: Qt = αB (blue) or Qt = Q(t)2 (red marks) Short-time weak constraint 4DVar - Q(t1, t2) ≈ Q0(t1)(t2)

mean quadratic error

10 3 10 2 10 1 100 101 102 103 0.05 0.1 0.15 0.2 0.25

/

tr = 10% Carrassi and Vannitsem, 2010 Carrassi et al. Model error in DA - ISDA 23rd January 2019 14 / 21

DA and model error Estimate Q using data

Time-batch estimated model error covariance

◮ The idea (Pulido et al, 2018) is to maximize the log-likelihood of the data (model evidence) as a function of the parameter θ l(θ) = ln

• p(xK:0, yK:1|θ)dxK:0

where θ can be λ, R or Q. ◮ Inserting an arbitrary PDF q(xK:0) and using the Jensen inequality we have l(θ) ≥

• q(xK:0) ln

p(xK:0, yK:1|θ) q(xK:0)

• dxK:0 ≡ Q(q, θ)

and the equality holds when q(xK:0) = p(xK:0|yK:1, θ) that is the PDF maximizing Q(q, θ) and a lower bound for l(θ). ◮ p(xK:0|yK:1, θ) can be obtained as the outcome of a DA procedure (e.g. EnKF, EnKS ...)

Carrassi et al. Model error in DA - ISDA 23rd January 2019 15 / 21

DA and model error Estimate Q using data

Time-batch estimated model error covariance Q

◮ This suggests a two-steps algorithms:

1 Expectation: Determine the distribution q that maximizes Q. This is given by

q∗ = p(xK:0|yK:1, θ

′). Note that p(xK:0|yK:1, θ ′) is the outcome (the posterior) of a data

assimilation algorithm for the HMM, evaluated at θ

′

2 Maximization: Determine the likelihood parameter θ∗ that maximizes Q(q∗, θ) over θ.

We have used the EnKF to estimate p(xK:0|yK:1, θ

′) in combination with:

the expectation–maximization, EnKF-EM the Newton–Raphson, EnKF-NR to maximize the likelihood associated to the parameters to be estimated.

Carrassi et al. Model error in DA - ISDA 23rd January 2019 16 / 21

DA and model error Estimate Q using data

Numeric with L96 model

20 40 60 80 100 EM iteration −6.8 −6.6 −6.4 −6.2 −6.0 −5.8 −5.6 Log-likelihood

(a)

20 40 60 80 100 EM iteration 10-4 10-3 10-2 10-1 ||Q−Qt ||F

(b)

K=100 K=500 K=1000 K=500, Ne=500 5 10 15 20 25 30 35 40 NR iteration −9.5 −9.0 −8.5 −8.0 −7.5 −7.0 −6.5 −6.0 −5.5 Log-likelihood

(a)

5 10 15 20 25 30 35 40 NR iteration 0.0 0.1 0.2 0.3 0.4 0.5 0.6 ||Q−Qt ||F

(b)

K=100 K=500 K=1000

◮ The EnKF-EM requires the optimal value in the maximization step to be computed analytically which limits the range of its applications ⇒ Ok in a Gaussian framework, an iterative minimization in nonlinear cases. ◮ In the EnKF-NR one makes use of approximate formulae for the model evidence. ◮ Convergence of the NR and EM maximization as a function of the iterations for different evidencing window lengths (K = 100, 500, 1000). ◮ (a) Log-likelihood function. ◮ (b) Frobenius norm of the model noise estimation error. ◮ In about 10 iterations, they converge to a good estimation.

Pulido et al, 2018 Carrassi et al. Model error in DA - ISDA 23rd January 2019 17 / 21

DA and model error Adaptive inflation

However always use inflation... (better if) adaptively

◮ Even with a good Q, you “always” need inflation due to sampling error and non-linearity/non-Gaussianity. ◮ Can avoid tuning by adaptive inflation; e.g. EAKF-adaptive by Anderson, 2007 or ETKF-adaptive by Miyoshi, 2011 ◮ A survey of existing methods in Raanes et al, 2019 ◮ Raanes et al, 2019 hybridized the “finite-size” EnKF-N (Bocquet, 2011) and the ETKF-adaptive ⇒ EnKF-N-hybrid targets explicit both sampling and model error

Carrassi et al. Model error in DA - ISDA 23rd January 2019 18 / 21

DA and model error Adaptive inflation

However always use inflation... (better if) adaptively

◮ L96 two scales with only the coarse scale observed, N = 20 members. ◮ EnKF-N-hybrid yields best filter accuracy, but only by slight margin. ◮ Results done and reproducible using the DA Python platform DAPPER. Codes hosted here. ◮ Figure from Raanes et al, 2019

5 10 15 20 25 30

Forcing (F) — both for truth and DA

0.0 0.2 0.4 0.6 0.8 1.0

RMSE ETKF excessive ETKF tuned EAKF adaptive ETKF adaptive EnKF-N hybrid

◮ Further sophistication of single-factor adaptive inflation estimation schemes is unlikely to yield significant, further improvements.

Carrassi et al. Model error in DA - ISDA 23rd January 2019 19 / 21

DA and model error Conclusion

Conclusion

Treating model error as a stochastic term is very convenient and coherent with a Bayesian formulation of the problem. But in many real problems (and particularly in climate science) model error is actually time-correlated and its impact grows with the prediction horizon. This is very much relevant when using smoother. Estimating the model error covariance matrix is extremely difficult in high-dimension. Classical state-augmentation approaches do not work well because the model error component of the error covariance is bound to monotonically decrease with time. A new method, based on the computation of the observational likelihood (the the model evidence) as a function of the model error covariance is introduced. The method requires the computation of the posterior that can be obtained (under Gaussian hypothesis) using EnKF, EnKS. Inflation is always needed to cope with non-Gaussianity and sampling error, but also for not-optimal Q. An extension of the EnKF-N originally devised for sampling error has been introduced to simultaneously deal with sampling and model error.

Carrassi et al. Model error in DA - ISDA 23rd January 2019 20 / 21

DA and model error Bibliography

Bibliography

Carrassi, A. and S. Vannitsem, 2010. Accounting for model error in variational data assimilation. A Deterministic Formulation. Mon. Weather. Rev., 138, 3369-3386 Carrassi A. and S. Vannitsem, 2011. State and parameter estimation with the extended Kalman filter using a deterministic approach. Q. J. Roy. Meteor. Soc., 137, 435-451 Mitchell, L. and A. Carrassi, 2015. Accounting for model error due to unresolved scales within ensemble Kalman filtering. Q. J. Roy. Meteor. Soc., 141, 1417–1428 Pulido, M., P. Tandeo, M. Bocquet, A. Carrassi and M. Lucini, 2018. Stochastic parametrization identification using ensemble Kalman filtering combined with expectation-minimization and Newton-Raphson maximum likelihood methods. Tellus, 70, 1442099 Raanes, P., M. Bocquet, and A. Carrassi, 2019. Adaptive covariance inflation in the ensemble Kalman filter by Gaussian scale mixtures. Q. J. R. Meteorol Soc., 2019, 1–23. https://doi.org/10.1002/qj.3386 Tandeo, P., P. Ailliot, M. Bocquet, A. Carrassi, T. Miyoshi, M. Pulido and Y. Zhen, 2019. Joint Estimation

• f Model and Observation Error Covariance Matrices in Data Assimilation: a Review. Submitted. Available

at https://arxiv.org/pdf/1807.11221.pdf

Carrassi et al. Model error in DA - ISDA 23rd January 2019 21 / 21