Data Assimilation by the Adjoint Method in Coastal and River Models

1. Data Assimilation by the Adjoint Method in Coastal and River Models by Graham Copeland and Igor Gejadze Workshop on Coastal Observatories. 17 – 19 October 2006, Proudman Oceanographic Laboratory. Department of Civil Engineering

2. Content Background to DA using Variational Methods:- Direct and Inverse Models The Variational Approach The Adjoint Method Examples of DA for the SWE and fsNSE Propagation of Adjoint Sensitivity Information Effective use of Data Model Coupling and Data Assimilation 1D St.Venant - 2D SWE

3. PROBLEM STATEMENT We require to either Reconstruct or Predict state variables (e.g. flow fields) based on:- OBSERVATIONS and GOVERNING EQUATIONS. This begs the questions:- What data are we likely to have? What data do we really need? How does the data control the solution? How do we make best use of data?

4. Limited Area Model Required Data – Model Controls Control Data Depth(x,y) Friction(x,y) Initial conditions: water levels & currents. Boundary conditions: currents or water levels . These values are needed to completely determine the solution but the boundary values are usually not all known from observations. Current v(0,y) Topography Depth(x,y) Friction coeff (x,y) Current v(L,y)

5. Direct and Inverse Models An initial or boundary value problem, known as a direct model, becomes inverse when:- some or all of its boundary or initial data is missing and is replaced by data located inside the domain. data data data data data data ill-posed ? Inverse also ill-posed

6. Direct and Inverse Models • So, since it most unlikely for all control data to be • available at the boundaries, all models are inverse. • However, in practice, reasonable assumptions are • made about the boundary controls and distributed • coefficients and so good solutions can be found. • But there are methods available that use available data • either • to recover the values of boundary controls • to recover the initial conditions • to recover values of distributed model parameters • to recover information about sources • and so lead to:- • assimilation of data into a model solution.

7. Variational Data Assimilation Methods Variational Methods have the following features:- They define an objective function J based on differences between the solution of the equations and observations i.e. a measure of ‘goodness of fit’. Using a variational method, they calculate a gradient measure that determines how to adjust the controls to minimise J and so improve the ‘goodness-of-fit’. They use an iterative procedure to move step by step towards a solution that minimises J.

8. The Adjoint Sensitivity Method This finds the gradient sensitivities of the type This sensitivity information is propagated away from the data locations back to the controls (e.g. at the boundary). The controls are then adjusted systematically and the direct model is recomputed. Eventually the values of controls are recovered that produce a solution that agrees with measured data. It carries out an Objective Calibration of the model. Useful for hindcasts and estimating initial conditions for forecasts

9. THE ADJOINT SENSITIVITY METHOD Shallow Water Equations Example 1- dimensional(x), unsteady(t) Integrate adjoint model backwards DATA x t initial conditions boundary conditions Integrate direct model forward Adjust initial and boundary conditions until Solution matches data

10. The Nonlinear shallow water equations: z Free surface t H (x,t) q Channel bed z(x,t) x

11. Measure data mismatch : by quadratic measuring function r Define the Objective function & the Lagrangian: Measure sensitivity to water level H(x0bs ,t0bs ) ≠ H0bs(x0bs ,t0bs) Or, if we have flow data qobs as well as Hobs r = 0.5 (H - Hobs)2(x – xobs) (t – tobs) + 0.5 (q - qobs)2(x – xobs) (t – tobs)

12. Take the first variation of J to get the adjoint equations in the Lagrange multipliers  and  (adjoint variables) Source terms

13. Sensitivities at Control Boundaries by Adjoint Method • Following the variational approach, the adjoint variables and  give us the sensitivities : ( m s2 ) the sensitivities of calculated water depths at x=xobs being different from an observed depthH0bsin terms of changes in inflow discharge “q” upstream. ( m2 s ) the sensitivities of water depths at x=xobs being different from an observed depth H0bs In terms of changes in water depth “H” downstream.

14. time x ‘Source’ for Adjoint (H-H0bs) Example: - Solution to forward equations Hydrograph amplitude = 2 m Trial Inflow hydrograph Calculated water depth (H) characteristic Discharge (q) Control (inflow) boundary

15. Source term & Adjoint Solution Source  t xobs x ‘sensitivity’ information how control boundaries affect levels at xobs ‘water level mismatch’ information creates a source that propagates backwards away from data towards boundaries

16. Solution to adjoint equations & Sensitivities on Control Boundaries Psi Phi   Results courtesy of H. Elhanafy, U of Strathclyde

17. Sinusoidal solution to forward equations characteristic Wave of amplitude = 1 m driven by data at x = 0 at all times x time Data time series Control boundary Data amplitude = 0.8 m

18. x t H(t)data Hdata   Final H Remainder of control not recovered Control recovered

19. y Surface elevation data u(0,y,t) Unknown boundary condition Current data bed y = 0 x Adjoint Method 2D section fsNSE Example problem definition p(x,y,t) Adjoint variables S*, u*, v*, p* 2D vertical section

20. Adjoint pressure p* at mid-depth based on velocity measurements Antiphase. Discontinuity in p* Adjoint solution showing propagation of p* through model

21. Wave transport Advection u = 1 ms-1 c = (g d)0.5 20 ms-1 S* y=0 control boundary p*, u*,v* u*,v* y=40 m 800 m data ‘source’ Transport mechanisms in adjoint model • Wave transport affects all depths • Advective transport preserves profile but only within advective range

22. u sensitivity through vertical at x = 0 due to a single velocity source at mid-depth at x = 800 m u - sensitivity at the control boundary. Wave disturbance arrives first as p* followed by convective terms as uu* Wave celerity = 20x convective speed Vertical spread due relative vertical motion of source, not diffusion

23. Separate influences on boundary control by advection and by wave propagation in 2D SWE Recovered by plane SW wave at speed c=(g h)0.5 Inflow control boundary Recovered by advection at speed u Observed Current u

24. inflow open Recovery of Boundary Controls in presence of standing waves

25. inflow closed Recovery of Boundary Controls in presence of standing waves. Data from nodes does not recover controls at all. Only H or only u does not recover controls very well Co-located H and u recovers controls best. Must retain phase info’

26. Recovery of distributed parameters by Adjoint DA For example recover y of bathymetry in SWE model Fom observations of water level H and velocities q/H we can recover bathymetry z(x). This assumes that the boundary controls are known

27. 2 D example. Reference bathymetry used in an identical twin experiment Flow field computed over this bathymetry and used to provide ‘observations’ of H and velocity at sparse locations Recovered z(x) DA proceeds from an initially plane bed to recover bathymetry from sparse observations of H

28. Initial z(x) Recovery of bathymetry by DA of co-located observations of velocities u,v and H Recovered z(x) Demonstrates the importance of having both u,v and H data Results courtesy of M. Honnorat, INRIA Rhone-Alpes, Grenoble

29. Flow model - 1 Surface water flow model – essential part of the full physical model used for eco-modelling.

30. Flow model – 2 It does not look always possible to create a single (integral) surface flow model based on most complete flow eqns (3D fsNSE). It does not look elegant either. Can we use the existing models (blocks) to simulate integral phenomena? Assuming all blocks are based on different eqns, implemented using different methods, run by different users in different places, runs are not necessarily synchronized, etc. No essential modifications in existing software are allowed! NO DEFINITE ANSWER ! Can we assimilate data collected in different parts of the flow to model integral phenomena without creating the single (integral) model? Paradoxical, but more definitely YES !

31. Simple flow model: 1D St.Venant-2D SWE This simple set contains most properties of the general case. Interface between models is ‘liquid’ or ‘open’ boundary. 1) What must be specified at the model interface to provide correct transfer of information (in both directions!)? Exact answer exists only for 1D SWE 2) Keep in mind: models are different. 1D St.V model cannot provide sufficient information to specify well-posed 2D SWE model by definition. 3) How to arrange information exchange in time? Time step is very different ! Different models could be run in different places!

32. Flow example with propagating dry/wet front: reference solution b) a) c) d)

33. Difficulties in coupling flow models - 1 Interface between models is ‘liquid’ or ‘open’ boundary. Therefore, all difficulties related to ‘open’ boundaries are applied here. What information should be transferred between models? How to synchronize information transfer? Trivial case: two overlapping sub-domains, information exchange every time step, explicit numerical scheme Why does it work? Because both Derichlet and Neumann specify well-posed semi-discrete (frozen time) problem in sub-domain. However, not a global time problem!!! General rule to follow: Boundary conditions applied at the ‘open boundary’ model interface must specify well-posed flow problem for each model (sub-domain). Derichlet or Neumann

34. Difficulties in coupling flow models - 2 Can we afford information exchange every time step? Do we use explicit solvers for all models? Are the models perfectly synchronized? NO Waveform relaxation method, or Global Time Schwarz method (E. Lelarasmee, L.Halpern, M. Gander, …) Need for ‘global time’ boundary conditions at the model interface, which specify a well-posed flow problem for each model. Characteristic analysis : systematic (but still approximate) way

35. Characteristics and higher order invariants – 1 Characteristic analysis of the 2D SWE: Eigenvalues: Invariants: Incoming invariants to be prescribed !, outgoing – extrapolated from the interior Characteristic analysis is approximate (even assuming small bathymetry and friction slopes) since x is considered as dominant incidence and flow direction. This results to the partial reflection (when using w as bc), which increases when incidence direction deviates from normal. Problem: flow may support many wave packages coming by different incidence angle.

36. Characteristics and higher order invariants – 2 Outgoing quantities from = Incoming quantities to Vice versa Using characteristic invariants at the interface between sub-domains is eventually equivalent to the waveform relaxation method. The solution can be obtained by successively solving problems in sub-domains (global time Schwarz ). In theory (for the 1D SWE), the process should converge in 2 iterations!

37. Characteristics and higher order invariants – 3 More complex invariants can be build based on the theory of absorbing BC (B. Engquist, A. Majda) Higher order invariants used as bc allow reflections of waves coming to the boundary by angles deviating from the normal to be essentially reduced. Example: second order invariants (E. Mazauric-Nourtier, build for linearized SWE ) (in) (in) (out) Remarks: • Still (x) remains the dominant flow and incidence direction. Improves treatment of waves in a certain cone around the normal. 2. Numerical implementation of higher order invariants (more than 2nd) could be unstable.

38. Difficulties in coupling flow models - 3 The main problem is 1D->2D information transfer. The 1D model cannot provide full set of boundary conditions for the 2D model : • For the 2D model all invariants are distributed in the tangential direction. • There is no tangential invariant in the 1D formulation exists. How to distribute along interface? Need distribution rule! Uniform?? 1D->2D 2D->1D No information available! Trivial?? No problem with it. No point in controlling

39. Model coupling: steady-state, non-homogenous, optimal control - 1 Virtual control methods: J-L. Lions, O. Pironneau, A. Quarteroni convection convection-diffusion - virtual controls Coupling is considered as the minimization problem:

40. Model coupling: steady-state, non-homogenous, optimal control – 2 or non-overlapping case A disjoint partition of G G 2 1 ' G W W 2 1 The objective function reflects our a-priori ‘natural’ guess about solutions at the coupling interface.

41. Model coupling: global time, non-homogenous, optimal control – 3 Open bc for 1D: Coupling is considered as the minimization problem for the objective functional: - fluxes Incoming invariant for 2D problem is the unknown control. Solution of the control problem gives result which tends to the solution by Schwarz iterations under trivial assumptions. Thus, by itself optimal control approach is not enough to get better results. Its advantage is that it may allow additional information to be used.

42. Model coupling & data assimilation If data is available and has to be assimilated, then coupling and assimilation can be performed in a single optimization loop. Extended objective functional: Could contain other controls, i.e. Common DA objective sensors This approach should be preferred for the non-homogenous model because the ‘rigid’ coupling at the interface could be eventually harmful. Otherwise, solution in 2D area is dominated by data, which fills in gaps in our knowledge about bc at the coupling interface.

43. Numerical results on WFR coupling

44. Generalized objective function DA component For the inlet coupling interface For the outlet coupling interface For the interfaces we use exactly what we know !

45. Joint DA-coupling algorithm – 1 (consistent 1D/2D meshes) Results of data assimilation: recovering the unknown inlet BC for the 1D model (together with BC for the 2D local model, which are not presented). Both elevation and discharge data are used. The convergence history: for the generalized objective function (in bold line) and for its components

46. Joint DA-coupling algorithm – 2 (consistent 1D/2D meshes) Results of data assimilation: recovering the unknown inlet BC for the 1D model (together with BC for the 2D local model, which are not presented). Only elevation data are used. The convergence history: for the generalized objective function (in bold line) and for its components

47. Conclusions • Coupling of flow models can be achieved together with DA within the same control loop. • The models and their adjoints run independently exchanging information between global time runs. No need to create single integrated model or generate its single integrated adjoint. • This arrangement should be preferred when incompatible (due to different reasons) models are involved into joint consideration. Measured data fills gaps in a-priori data needed in ‘coarse-to-fine’ information transfer. Thus, coupling together with DA is simpler task than without. • Boundary conditions applied on the model interfaces should be based on characteristics or higher order characteristic invariants. • In the overlapping areas the ‘finer’ model should adjust ‘coarser’ model via ‘defect correction’