Introduction to Data Assimilation NCEO Data-assimilation training days 11-12 July 2011

1. Introduction to Data AssimilationNCEO Data-assimilation training days11-12 July 2011 Peter Jan van Leeuwen Data Assimilation Research Center (DARC) University of Reading p.j.vanleeuwen@reading.ac.uk

2. Basic estimation theory Observation: T0 = T + e0 E{e02} = s02 E{em2} = sm2 E{e0em} = 0 E{e0} = 0 E{em} = 0 First guess: Tm = T + em Assume a linear best estimate: Tn = a T0 + b Tm with Tn = T + en . Find a and b such that 1) E{en} = 0 2) E{en2} minimal 1) Gives: E{en} = E{Tn-T} = E{aT0+bTm-T} = = E{ae0+bem+ (a+b-1) T} = (a+b-1) T = 0 Hence b=1-a.

3. Basic estimation theory 2) E{en2} minimal gives: E{en2} = E{(Tn –T)2} = E{(aT0 + bTm – T)2} = = E{(ae0+bem)2} = a2 E{e02} + b2E{em2} = = a2 s02 + (1-a)2 sm2 This has to be minimal, so the derivative wrt a has to be zero: 2 a s02 - 2(1-a) sm2 = 0, so (s02 + sm2)a – sm2 = 0, hence: sm2 s02 a = __________ and b = 1-a = __________ s02 + sm2 s02 + sm2 sm4 s02 + s04 sm2 s02 sm2 sn2 = E{en2} = ______________ = ________ (s02 + sm2)2 s02 + sm2

4. Basic estimation theory sm2 s02 Solution: Tn = _______ T0 + _______ Tm s02 + sm2 s02 + sm2 1 1 1 ___ = ___ + ___ and sn2 s02 sm2 Note: sn smaller than s0 and sm ! Best Linear Unbiased Estimate BLUE Just least squares!!!

5. Can we generalize this? • More dimensions • Nonlinear estimates (why linear?) • Observations that are not directly modeled • Biases • (Is this reference to the truth necessary?)

6. The basics: probability density functions P(u) 0.5 1.0 u (m/s)

7. The model pdf p (u(x1), u(x2), T(x3), … ) u(x2) u(x1) T(x3)

8. Observations • In situ observations: • irregular in space and time • e.g. sparse hydrographic observations, • Satellite observations: indirect • e.g. of the sea-surface

9. Data assimilation: general formulation Bayes theorem: The solution is a pdf!

10. Bayes’ Theorem Conditional pdf: Similarly: Combine: Even better:

11. Filters and smoothers Filter: solve 3D problem sequentially Smoother: solve 4D problem all at once Note: I assumed a dynamical (=time-varying) model here.

12. The Gaussian assumption I Most present-day data-assimilation methods assume that the model and observation pdf are Gaussian distributed. Advantage: only first two moments, mean and covariance, needed. Disadvantage: This is not always a good approximation. Gaussian distribution: Multivariate Gaussian for n-dimensional model vector (state or parameters or a combination):

13. The Gaussian assumption II The pdf of the observations is actually (and luckily) the pdf of the observations given the model state. Hence we need the model equivalent of the observations. We can express each observation as: in which d is an m-dimensional observation vector. The observation pdf is now given by:

14. (Ensemble) Kalman Filter I Use Gaussianity in Bayes at a specific time: Hence Complete the squares to find again a Gaussian (only for linear H !!!):

15. (Ensemble) Kalman Filter III • Two possibilities to find the expressions for the mean and covariance: • Completing the squares • Writing the new solution as a linear combination of model • and observations, assume unbiased and minimal trace of • error covariance. Both lead to the Kalman filter equations, which again are just the least squares solutions:

16. The error covariance:tells us how model variables co-vary Spatial correlation of SSHand SST in the Indian Ocean x x Likewise you can find the covariance of a model parameter with any other model variable. Haugen and Evensen, 2002

17. (Ensemble) Kalman Filter IV Interestingly sm2 s02 Old solution: Tn = _______ T0 + _______ Tm sm2 + s02 sm2 + s02 Hence we find the same solution, but now generalized to more dimensions and partial or indirect observations (H operator)!

18. (Ensemble) Kalman Filter V • Operational steps: • Determine mean and covariance from the prior pdf • Use Kalman filter equations to update the model variables • Propagate the new state forward in time using the model equations • Propagate the new covariance in time using the linearized model • equations • Back to 2) • In the ensemble Kalman filter the mean and covariance are found from an ensemble of model states/runs.

19. And parameter estimation? Augment the state vector with the parameters (i.e. add parameters to the state vector). This makes the parameter estimates time variable. How do we find the covariances? 1) Prescribe the covariances between parameters and observed variables, or 2) Find best fits using Monte-Carlo runs (lots of tries, can be expensive).

20. Nonlinear measurement operators H (but observations still Gaussian distributed) • 1) Use the nonlinear H in the Kalman filter equations • (no rigorous derivation possible). • 2) Augment the state vector x with H(x) • (still assuming Gausian pdf for new model vector). • Solve the nonlinear problem. For large-scale problems • variational methods tend to be efficient: 3DVar (and 4DVar). • (Sometimes also used when H is linear.)

21. Variational methods A variational method looks for the most probable state, which is the maximum of this posterior pdf also called the mode. Instead of looking for the maximum one solves for the minimum of a so-called costfunction. The pdf can be rewritten as in which Find min J from variational derivative: J is costfunction or penalty function

22. Gradient descent methods J 1’ 3 4 6 5 2 1 model variable

23. 3DVar • 3DVar solves for the minimum of the costfunction, which is the • maximum of the posterior pdf. • Advantages: • When preconditioned appropriately it can handle very • large dimensional systems (e.g. 10^8 variables) • It uses an adjoint, which tells us the linear sensitivity • of each model variable to all other model variables. • This is useful for e.g. observation impact. • Disadvantages: • The adjoint has to be coded • It doesn’t automatically provide an error estimate. For weakly • nonlinear systems the inverse of the 2nd variational derivative, • the Hessian, is equal to the covariance, but it comes at great cost.

24. 4DVar There is an interesting extension to this formulation to a smoother. Filter: Solve a 3D problem at each observation time. 4DVar: Solve a 3D problem at the beginning of the time window using all observations. Note we get a new full 4D solution!

25. 4DVar: The observation pdf The part of the costfunction related to the observations can again be written as: But now we have sets of observations at different times tobs, and H is a complex operator that includes the forward model. If observation errors at different times are uncorrelated we find:

26. 4DVar: the dynamical model The dynamical model is denoted by M: Using the model operator twice brings us to the next time step: And some short-hand notation:

27. 4DVar The total costfunction that we have to minimize now becomes: in which the measurement operator Hi contains the forward model: defined as

28. Weak-constraint 4DVar So far, we have considered errors in observations and the model initial condition. What about the forward model? If the forward model has errors, i.e. if there are errors present in the model equations the model pdf has to be extended. This leads to extra terms in the costfunction:

29. Parameter estimation Bayes: Looks simple, but we don’t observe model parameters…. We observe model fields, so: in which H has to be found from model integrations.

30. And what if the distributions are strongly non-Gaussian? Fully non-Gaussian data-assimilation methods exist, but none has proven to work for large-scale systems. A few do have potential, among which the Particle Filter.

31. Summary and outlook • We know how to formulate the data assimilation problem • using Bayses Theorem. • We have derived the Kalman Filter and shown that it is • the best linear unbiased estimator (BLUE). • We derived 3D and 4DVar and discussed some of their • properties. • This forms the basis for what is to come the rest of the day! • Notation: when we talk about the state vector and its • observations notation will be changed: • becomes x, and d becomes y. Beware !!!