Title: Data Assimilation by the Adjoint Method in Coastal and River Models
1Data 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
2Content 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
3PROBLEM 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?
4Limited 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)
5Direct 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
6Direct 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.
7Variational 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.
8The 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
9THE 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
10The Nonlinear shallow water equations
z
Free surface
t
H (x,t)
q
Channel bed
z(x,t)
x
11Define 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)
12Take the first variation of J to get the
adjoint equations in the Lagrange multipliers ?
and ? (adjoint variables)
Source terms
13Sensitivities 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.
14Example - Solution to forward equations
Hydrograph amplitude 2 m
Trial Inflow hydrograph
Calculated water depth (H)
characteristic
Discharge (q)
Control (inflow) boundary
15Source 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
16Solution to adjoint equations Sensitivities on
Control Boundaries
Psi
Phi
?
?
Results courtesy of H. Elhanafy, U of Strathclyde
17Sinusoidal 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
18x
t
H(t) data
H data
?
?
Final H
Remainder of control not recovered
Control recovered
19Adjoint Method 2D section fsNSE Example problem
definition
p(x,y,t)
Adjoint variables S, u, v, p
2D vertical section
20Adjoint pressure p at mid-depth based on
velocity measurements
Antiphase. Discontinuity in p
Adjoint solution showing propagation of p
through model
21Adjoint pressure p at mid-depth based on tidal
elevation measurements
Same phase. Continuous p
Adjoint solution showing propagation of p
through model
22Wave 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
23u 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
24Separate 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
25inflow
open
Recovery of Boundary Controls in presence of
standing waves
26inflow
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
27Recovery 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
282 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
29Initial 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
30Flow model - 1
Surface water flow model essential part of the
full physical model used for eco-modelling.
31Flow 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 !
32Simple 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!
33Flow example with propagating dry/wet front
reference solution
b)
a)
c)
d)
34Difficulties 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
35Difficulties 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
36Characteristics 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.
37Characteristics 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!
38Characteristics 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.
39Difficulties 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
40Model 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
41Model 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.
42Model 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.
43Model 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.
44Numerical results on WFR coupling
45Generalized objective function
DA component
For the inlet coupling interface
For the outlet coupling interface
For the interfaces we use exactly what we know !
46Joint 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
47Joint 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
48Conclusions
- 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