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
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
xxobs being different from an observed depth
H0bs in terms of changes in inflow discharge q
upstream.
( m2 s )
the sensitivities of water depths at xxobs
being different from an observed depth H0bs In
terms of changes in water depth H downstream.
14
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
H data
?
?
Final H
Remainder of control not recovered
Control recovered
19
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
Adjoint pressure p at mid-depth based on tidal
elevation measurements
Same phase. Continuous p
Adjoint solution showing propagation of p
through model
22
Wave transport
Advection u 1 ms-1
c (g d)0.5 ? 20 ms-1
S
y0
control boundary
p, u,v
u,v
y40 m
800 m
data source

• Transport mechanisms in adjoint model
• Wave transport affects all depths
• Advective transport preserves profile but only
  within advective range

23
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
24
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
25
inflow
open
Recovery of Boundary Controls in presence of
standing waves
26
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
27
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
28
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
29
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
30
Flow model - 1
Surface water flow model essential part of the
full physical model used for eco-modelling.
31
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 !
32
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!
33
Flow example with propagating dry/wet front
reference solution
b)
a)
c)
d)
34
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
35
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
36
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.
37
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!
38
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.

39
Difficulties in coupling flow models - 3

• The main problem is 1D-gt2D 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-gt2D 2D-gt1D
No information available! Trivial??
No problem with it. No point in controlling
40
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
41
Model coupling steady-state, non-homogenous,
optimal control 2
or non-overlapping case
A disjoint partition of
G
G
2
1
'
W
G
W
2
1
The objective function reflects our a-priori
natural guess about solutions at the coupling
interface.
42
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.
43
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.
44
Numerical results on WFR coupling
45
Generalized objective function
DA component
For the inlet coupling interface
For the outlet coupling interface
For the interfaces we use exactly what we know !
46
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
47
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
48
Conclusions

1. Coupling of flow models can be achieved together
   with DA within the same control loop.
2. 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.
3. 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.
4. Boundary conditions applied on the model
   interfaces should be based on characteristics or
   higher order characteristic invariants.
5. In the overlapping areas the finer model should
   adjust coarser model via defect correction