Start Submission Become a Reviewer

Reading: Variational Assimilation of Surface Wave Data for Bathymetry Reconstruction. Part II: Second...

Download

A- A+
Alt. Display

Original Research Papers

Variational Assimilation of Surface Wave Data for Bathymetry Reconstruction. Part II: Second Order Adjoint Sensitivity Analysis

Authors:

R. A. Khan ,

Department of Meteorology (MISU), Stockholm University, Stockholm, SE; Department of Mathematics and Statistics, McMaster University, Hamilton, ON, CA
X close

N. K.-R. Kevlahan

Department of Mathematics and Statistics, McMaster University, Hamilton, ON, CA
X close

Abstract

Data assimilation methods have been proposed as a technique for reconstructing ocean bathymetry from observations of surface waves. To better understand this technique, we use second order adjoint (SOA) variational analysis to derive the sensitivity of the surface wave error to perturbations in the observations (such as their number, spacing and position relative to bathymetry profiles), given the reconstructed bathymetry. We apply SOA to the data assimilation scheme for the one-dimensional shallow water equations for bathymetry detection introduced in Khan and Kevlahan (2021). We derive the Hessian of a cost function J representing the error between forecast surface wave and the observations. We then use SOA to derive the sensitivity of the surface wave error given the reconstructed bathymetry to perturbations in the observations for both a compactly supported Gaussian bathymetry, and a sandbar profile bathymetry. We investigate the correlation between (i) low sensitivity of the surface wave given the reconstructed bathymetry, to the observations, and (ii) the error in the bathymetry reconstruction, as well as the sensitivity of the data assimilation scheme to perturbations of its parameters. Additionally, we determine whether the conclusions in Khan and Kevlahan (2021) for bathymetry reconstruction can be verified by the present sensitivity analysis. We observe that relatively large errors in the bathymetry reconstruction and large relative amplitudes of the Gaussian and sandbar bathymetry profiles are associated with higher sensitivity of the surface wave reconstruction error to the observations. However, sensitivity decreases when the observation network has a greater coverage of the bathymetry. Significantly, the sensitivity of the surface wave to the observations is orders of magnitude lower than the bathymetry reconstruction error itself. These results suggest optimal configurations of surface wave observations, help minimise costs for making observations, and could enhance the accuracy of tsunami models.

How to Cite: Khan, R.A. and Kevlahan, N.K.-R., 2022. Variational Assimilation of Surface Wave Data for Bathymetry Reconstruction. Part II: Second Order Adjoint Sensitivity Analysis. Tellus A: Dynamic Meteorology and Oceanography, 74(1), pp.187–203. DOI: http://doi.org/10.16993/tellusa.36
33
Views
4
Downloads
  Published on 14 Apr 2022
 Accepted on 18 Feb 2022            Submitted on 18 Feb 2022

1 Introduction

Data assimilation is integral to accurate climate, atmosphere and ocean modelling. Variational data assimilation algorithms such as 3D-VAR, 4D-VAR, and Kalman filtering techniques like the Ensemble Kalman filter (EnKF) are regularly used for numerical weather prediction and forecasts of climate trends. For example, data assimilation is used by the European Centre for Medium Range Weather Forecasting (ECMWF) for climate reanalysis, where archived observations are ‘reanalysed’, in order to create a comprehensive global data set describing the recent history of the earth’s climate, atmosphere and oceans. Data assimilation is used in tsunami forecast models, where observations of surface waves are used to reconstruct missing information such as initial conditions, and subsequently predict impact at coastlines (Nakamura et al., 2006). Tsunami modelling requires accurate knowledge of the bathymetry, and data assimilation has also been proposed as a technique to reconstruct bathymetry from measurements of surface waves. The goal is to reconstruct the bathymetry sufficiently accurately to generate precise predictions of the surface wave. This paper focuses on the sensitivity problem for bathymetry reconstruction.

In Part I (Khan and Kevlahan, 2021) of this study we developed and evaluated a variational data assimilation algorithm for the one-dimensional nonlinear shallow water equations. This algorithm estimates bathymetry using a finite set of observations of surface wave height. Our goal was to investigate whether variational data assimilation is a feasible method to estimate ocean bathymetry, in the shallow water model. The computational results presented in Part I indicated that convergence to the true bathymetry is improved by including more observation locations positioned ahead of the bathymetry. A necessary condition for convergence of the bathymetry reconstruction is that the amplitude of the initial conditions is less than 1% of the bathymetry height. Additionally, we investigated the effect of varying the amplitudes of the bathymetry and initial condition, and concluded that a necessary condition for data assimilation convergence was that the amplitude of the initial conditions be at least two orders of magnitude smaller than the normalised amplitude of the bathymetry. We then used density-based global sensitivity analysis (GSA) to quantify the sensitivity of the surface wave to the bathymetry reconstruction error. We concluded that reconstructing the bathymetry with a relative error of about 10% is sufficiently accurate for surface wave modelling in most cases.

Our objective for Part II of this study is to more rigorously quantify the sensitivity of the assimilation algorithm results to perturbations to the observations, such as to their number, spacing and positions relative to the bathymetry. To do this we use the second order adjoint (SOA) method developed by Shutyaev et al. (2018). The sensitivity analysis in the present study is complementary to the previous work since it quantifies how the results of the bathymetry assimilation change when the observations are perturbed.

In tsunami modelling we wish to choose the number and positions of the observation network such that the accuracy of the surface wave forecast is maximised, given some optimal reconstruction of the bathymetry. In order to solve this problem we need to quantify how sensitive the surface wave reconstruction is, given the reconstructed bathymetry, to the choice of observation configuration. To do this we define a response functional G as the least squares error in the surface wave height η(x, t) based on the reconstructed bathymetry β(b) compared to that produced by the true bathymetry β(t),

(1.1)
Gη,  u,  β(b)x=0Tηβ(t)ηβ(b)2dt,

where η(x, t) depends on its velocity u(x, t) through the shallow water equations. Let us define mj(t) to be observations of the surface wave at positions {xj} for j = 1, …, Nobs at continuous times t. (Recall that the reconstructed bathymetry β(b) is itself computed from observations of the surface wave using a variational data assimilation scheme.) Then,

dG/dm is the sensitivity of the error in the surface wave to perturbations in the observations {mj(t)}, when the bathymetry is reconstructed from observations of the surface wave using the assimilation algorithm.

Our goal in this paper is to construct the response functional dG/dm, and then characterise its dependence on various types of perturbations of the observations {mj(t)}. These perturbations include changes in the positions of the observations {xj}, or perturbations in the shape of the true bathymetry (such as the amplitude or standard deviation of a Gaussian bathymetry profile) that result in a change in {mj(t)}. The choice of continuous or discrete time does not fundamentally change the results, and is made for mathematical convenience. In the discrete case (dG/dm)[x] is an Nobs × Nt matrix (a function of the space variable x), where Nt is the temporal resolution (e.g. the number of time steps), and the mi,j entry of dG/dm represents the sensitivity of the surface wave error (4.1) to the i-th observation at the j-th time step.

Our analysis answers the following three questions:

  1. How does the sensitivity of the surface wave error to perturbations in the measurements depend on the magnitude of the error of the reconstructed bathymetry?
  2. How does the sensitivity of the surface wave error depend on changes to parameters in the assimilation scheme (e.g. coverage of the bathymetry feature by the the observation points)? Is there an optimal spacing of observations that minimises the sensitivity?
  3. Are our qualitative computational observations from Khan and Kevlahan (2021) verified by the results of the sensitivity analysis?

In this study, we use the SOA algorithms outlined in Shutyaev et al. (2018) to analytically derive the Hessian of the cost function (2.4) minimised in the data assimilation scheme. We then take advantage of its properties to extract expressions for dG/dm. We present a numerical implementation of this algorithm for bathymetry data assimilation.

The variational data assimilation scheme is summarised in Section 2. Section 3 gives the analytical derivation of the Hessian and subsequent SOA sensitivity analysis for the bathymetry reconstruction data assimilation results. Section 4 presents the numerical implementation of the sensitivity analysis for the bathymetry data assimilation. Finally, Section 5 summarises the main results and suggests further considerations for future analyses.

2 Variational data assimilation for bathymetry reconstruction

We briefly summarise the data assimilation scheme for bathymetry reconstruction as implemented in Khan and Kevlahan (2021). The 1-D irrotational and incompressible nonlinear shallow water equations are,

Where β is the bathymetry to be reconstructed. We assume that the initial conditions ϕ(x) are compactly supported, and that the boundary conditions are periodic on a domain Ω = {x ∈ ℝ; –LxL}. The system is integrated in time for t ∈ [0 T], where T = 2L, and

(2.2)
Y=L2Ω  ×  0,T,
(2.3)
Yp=L2Ω

are the state space and the bathymetry function space respectively.

In order to find classical solutions of (2.1), we require the gradient J(β) to be in the Sobolev space H2(Ω), which imposes additional smoothness requirements on β and its derivative βx. In Khan and Kevlahan (2021), we implement a low pass filter, which removes higher frequencies in our reconstruction and effectively increases the regularity of our estimate of β at each iteration, raising it from the space L2(Ω) to H2(Ω).

We assimilate a set of measurements y(o)(t), which are observations of the true surface height perturbation η(x, t) at discrete positions {xj}, j = 1, …, Nobs and times t. This representation of the observation operator in time is convenient in the “optimise and then discretise” approach we have used for our variational assimilation scheme, and does not impose any smoothness requirements on the time dependence of the measurements. Discrete measurements at times tk can be represented easily as Dirac delta functions δ(ttk) in the continuous integrals representing the cost function (2.4).

We assume that we do not have complete information about β(x), and our objective is to minimise the error between the forecast solution η(x, t) given some guess for β, and the observations y(o)(t). We define this error in terms of a cost function,

(2.4)
Jβ=120Tj=1Mδxxjηx,  t;  βyj(o)t2dt,

where the Dirac delta functional is used to sample the state variable η at the same spatial locations as the observation points, and the summation from j = 1 to M represents the cumulative error at all observation points. The optimal bathymetry that minimises this cost function is the control variable satisfying

(2.5)
JL2β(b)=0,

where β(b) is the optimal approximation of the bathymetry β(x), and is the local minimizer of (2.4).

Since the focus of this work is the relationship between the error in the bathymetry reconstruction and the error in the resulting surface wave, we assume that the observations are free of noise. Actual applications would require stabilisation of noise present in the real-world measurements, which can be achieved by inclusion of a Tikhonov regularisation term in (2.4).

We note that the variational approach only works for continuous solutions, because the integration by parts used to extract ∇JL2(β) fails for discontinuous solutions. In principle, any discontinuities could be handled by including the position of any discontinuities as new variables, and then integrating by parts separately on each side of the discontinuity. This would require additional “boundary conditions” connecting the two solutions. However, in this paper we assume that the initial conditions and integration times are such that the solution remains continuous. In practice, there is always some diffusion in realistic simulations, either numerical or explicit horizontal Laplacian or hyperdiffusion, that ensures solutions remain smooth.

We formulate a Lagrangian constrained by (2.1) and some appropriately chosen adjoint variables (Lagrange multipliers) (η*, u*) that are solutions of

System (2.6) is chosen such that when the Lagrangian is integrated by parts in space and time, we are able to use the Riesz representation theorem and the Gâteaux derivative representation of J′(β; β′) given an arbitrary perturbation β′ to derive

(2.7)
L2J=0Tuη*xdt,

where

(2.8)
Jβ;β=Jβ,βL2(Ω),                      =LLL2Jββdx.

For a more detailed analysis of this derivation we refer the reader to Khan and Kevlahan (2021). Finally, the optimal reconstruction of the bathymetry β(b) is where

(2.9)
L2Jβ(o)=0Tuη*xdt=0.

In the assimilation scheme we use a gradient descent algorithm to iteratively find the optimal reconstruction of the initial conditions β(b) given some initial guess, such that (2.4) is minimised.

3 Second order adjoint sensitivity analysis

To derive the sensitivity the surface wave error to to perturbations in the observations, given the reconstructed bathymetry, using the methods outlined in Shutyaev et al. (2018) we need to formulate expressions for the Hessian of the cost function (2.4). Shutyaev et al. (2018) provide a general method based on properties of the Hessian, however they do not provide a derivation of the Hessian itself, which means we need to extend our variational adjoint analysis used to find the gradient of J, to find the Hessian of J.

While works such as Wang et al. (1992) provide a derivation of the Hessian vector product for the initial conditions assimilation, their derivation is for the finite dimensional case, and assumes a vector form for both the state variables and the control variable. In our case the derivation of the first order adjoint is for the infinite dimensional case in the space L2(Ω) over some domain Ω, where we used the L2 inner product and the Riesz representation theorem to extract our gradient ∇L2J(ϕ), and subsequently used a low pass filter in the numerical implementation to raise it from the space L2(Ω) to H2(Ω). For that reason, a derivation of the Hessian in the same functional space as for our first order adjoint is appropriate, and in doing so we aim to extract the ‘Gâteaux Hessian’ for Hilbert spaces with the following definition:

If f is twice Gâteaux differentiable at x, we can identify D2f(x) with the operator ∇2f(x) ∈ B(H) in the sense that

(3.1)
yHzH,D2fxyz=z,2fxyB(H),

where we call ∇2f(x) the (Gâteaux) Hessian of f at x, B(H) is the space of continuous linear functionals in H (in our case L2(Ω), and Df(x)y is the Gâteaux derivative of f in the direction y.

In the remainder of this section we present the derivation of the Hessian of J(β), and subsequently the sensitivity analysis for the optimal bathymetry reconstructed using data assimilation of surface wave measurements, using methods outlined in (Shutyaev et al., 2018).

3.1 Hessian of J(β) for bathymetry reconstruction

We know that the Gâteaux derivative of J with respect to the bathymetry and some perturbation direction β′ is

(3.2)
Jβ;β=LL0Tuη*xdt  βdx.

Consider a second perturbation of J′(β; β′), β^ where we have ββ+εβ^. Then the second order Gâteaux derivative of J is

Jβ;β;β^=ddεJβ+εβ^;β|ε=0                            =ddεLL0Tux,  t;β+εβ^η*xx,  t;β+εβ^dtβdxε=0

We consider a regular perturbation expansion of the integrand, approximating it by the series f0 + f1ϵ + 𝒪(ϵ2). We can see this is equivalent to a Taylor expansion about ϵ = 0. Then

    ddεJβ+εβ^;β|ε=0=ddεLL0Tfo+f1εdtβdx|ε=0=LL0Tf1dtβdx.

In order to find the term f1, we recognise that this is equivalent to the coefficient of the linear term in the Taylor approximation, and f1=ddεux,  t;β+εβ^xη*x,  t;β+εβ^|ε=0. Thus, we require u and adjoint variable η* given the perturbed bathymetry β+εβ^. To find the resulting forward and adjoint systems given the perturbation, we assume this perturbation β^ in the bathymetry produces the following perturbations to our state and adjoint variables;

  • uu + û and ηη + η^ for the shallow water system.
  • u* → u* + ū and η* → η* + η¯ for the adjoint system.

The resulting perturbed model for the state variables û, η¯ is

and the second order adjoint (SOA) model is

Giving us

(3.5)
2Jβ;β^=0Tu^η*x+uη¯xdt.

We define the Hessian H acting on the perturbation β^ as the successive solution of the perturbed and SOA models such that

(3.6)
Hβ^=0Tu^η*x+uη¯xdt.

This form of the Hessian of the cost function can be verified using a modified form of the kappa test analysis outlined in Kevlahan et al. (2019),

(3.7)
κε=limε01εJβ+εβ^;  βJβ;  ηLL0Tu^η*x+uη¯xdtβdx.

If the derivations are correct, κ(ϵ) should converge to 1 as ϵ → 0. The results for the kappa test for the Hessian H are presented in Section 4 in Figure 1(d).

a–c: differing shapes for initial condition and bathymetry. d–f: numerical convergence and verification of algorithm 1
Figure 1 

(a)–(c): The three test cases for bathymetry β(x) (dashed line) and surface wave initial conditions ϕ(x) (solid line) for the data assimilation. The surface wave initial conditions η^, bathymetry β^, and average depth H are not to scale in these diagrams, as η^ was restricted to 1% of β^ across most of the numerical analyses, and β = 0.1. (d): The log error in the convergence of the kappa test in (3.1), to verify the numerical calculation of the Hessian. (e) and (f):Hν (red) and F (blue) for cases I and II where ν is the solution of Hν = F found using the matlab linear solver bicgstabl.

A key observation is that we have only derived the action of the Hessian on some perturbation β^, whereas in Khan and Kevlahan (2021) we were able to derive the gradient of J(ϕ) for any arbitrary perturbation η′. It does not seem possible to proceed as before and find the Hessian for any arbitrary perturbation using variational methods. This is because the gradient ∇J is an element of the Hilbert space L2(Ω), whereas the Hessian H is an operator defined on this space. Indeed, Wang et al. (1992) were also limited to derivation of a ‘Hessian vector product’ in the finite dimensional case. Despite this, the current derivation is sufficient for the following sensitivity analysis.

3.2 Sensitivity analysis for bathymetry reconstruction

Define the optimality system for bathymetry assimilation as the successive solution of

Where λ(x) is the optimal reconstruction of the bathymetry β(x) giving:

(3.10)
0Tuη*xdt=0,

and m(t) are the observations taken at positions {xj} at for j = 1, …, Nobs at continuous times t. We impose periodic boundary conditions in space, where the domain is Ω = {x ∈ ℝ; x ∈ [–L, L]}.

Let us consider some arbitrary response functional G(η, u, λ). Then, by the chain rule, the sensitivity of G to perturbations in the observations m can be defined as

(3.11)
dGdm=dGdηdηdm+dGdududm+dGdλdλdm.

Let us now consider a perturbation in the observations mm + m^ giving

  • uu + û; and ηη + η^ for the shallow water system (3.8).
  • u* → u* + u˜* and η* → η* + η˜* for the adjoint system (3.9).
  • λ → β + λ^ for the optimal bathymetry.

The perturbed system becomes

(3.14)
0Tu^η*x+uη˜*xdt=0.

Then we can say

(3.15)
dGdm,m^=Gx,η^Y+Gu,u^Y+Gλ,λ^YP

where Y and Yp are the state space, and the bathymetry function space as defined in (2.2), and (2.3) respectively. Now, let us introduce some adjoint variables Pi, i = 1, …, 5, where PiY, i = 1, …, 4, and P5Yp. Then if we take the inner product of P1 and P2 with the systems (3.12), inner product of P3 and P4 with (3.13), and P5 with (3.14), we get the following duality relation

(3.16)
0TLLP1η^t+xuη^+x(η+1λ)u^xλ^uxλu^                    +P2u^t+u^ux+η^x                      +P3η˜*t+u^η*x+uη˜*x+u˜*xj=1Mη^xj,t;  λ^m^jtδxxj                      +P4u˜*t+η+1λη˜*x+η^η*x+uu˜*x+u^u*xλ^η*x                      +P5u^η*x+uη˜*xdtdx=0.

Integrating (3.16) by parts in space and time, we are able to transfer the derivatives onto the adjoint variables Pi, i = 1, …, 5 instead of on û, η^, u˜*, u˜*. If we pick the following systems for Pi,

(3.19)
0TP4η*x+uP1xdt=Gλ,

then as a result of integration by parts and the choice of systems for Pi, (3.16) reduces to

(3.20)
j=1MP3xj,t,m^=Gη,η^Y+Gu,u^Y+Gλ,λ^Yp.

By the Riesz representation theorem and equivalence of inner products in (3.15) and (3.20), we define the sensitivity of the response functional G(η, u, λ) to perturbations in the observations m as

(3.21)
Gm=j=1MP3xj,t.

We note that, because of the Dirac delta functional, the integrals represented by the inner products on the right-hand side of (3.20) are well defined, even with discrete observations in space. As we are only perturbing the locations of the measurements (i.e. the parameters xj of the Dirac function), this should not affect the smoothness required in the inner product on the left-hand side, which logically should be in the intersection of the spaces Y and Yp.

Equations (3.18) and (3.17) are a coupled system of four variables, with two initial time conditions and two final time conditions, making it challenging to solve. However, we observe that (3.17) is equivalent to the perturbed system for the Hessian ∇2J(ϕ) (3.1) with P5 = β^(x), and (3.18) is equivalent to the second order adjoint (SOA) system (3.1) with forcing term (–G/∂η, –G/∂u)T. Let us replace P5 with the auxiliary variable ν. Then Shutyaev et al. (2018) show that the solutions to the adjoint systems (3.17) and (3.18) are equivalent to solving

(3.22)
Hν=F,

where F is defined as

(3.23)
F=Gλ0Tuγxdt,

and γ is the solution of the forced first order adjoint system

Under the assumption H is positive definite (as our cost function J is convex), we can find a unique ν for every F such that Hν = F Shutyaev et al., (2018). We are then able to find the sensitivity of the surface wave to perturbations in the measurements, G/∂m, by the steps outlined in Algorithm 1.

Algorithm 1

Calculation of Second Order Adjoint Sensitivity mG for Bathymetry Assimilation.


1: Define λG0T  uxγdt, where γ is the solution of (3.24).
2: Solve H ν = F for ν, where H is the Hessian operator acting on ν, and F is the forcing term defined by (3.2) in step 1.
3: Solve the system (3.17) by substituting the control variable P5(x) with ν (as found in step 2) to find the adjoint variable P3(x, t).
4: Define mG=j=1MP3x,t, where P3 has been sampled at the locations of the observation points {xj}.

4 Results

The results of this analysis are a significant extension and complement to the work undertaken in Khan and Kevlahan (2021). In this section, we present the numerical implementation of the second order adjoint sensitivity analysis for the bathymetry assimilation with a travelling surface wave over a Gaussian bathymetry and a sandbar profile bathymetry.

We recall that the response functional G is the least squares error between the surface wave given the true bathymetry, and the surface wave resulting from the optimally reconstructed bathymetry,

(4.1)
Gη,  u,  β(b)x=0Tηβ(t)ηβ(b)2dt,

and that G depends implicitly on the measurements mj(t) (taken at positions {xj} for j = 1, …, Nobs at times t, which may be either continuous or discrete), through the reconstructed bathymetry η(β(b)). We now use the Algorithm 1 to find G/∂m, the sensitivity of the error in the surface wave to perturbations in the observations {mj(t)}, when the bathymetry has been reconstructed using the assimilation algorithm outlined in Khan and Kevlahan (2021).

As outlined in Section 1, our analysis aims to answer the following three questions:

  1. How does the sensitivity of the surface wave error to perturbations in the measurements depend on the magnitude of the error of the reconstructed bathymetry?
  2. How does the sensitivity of the surface wave error depend on changes to parameters in the assimilation scheme (e.g. coverage of the bathymetry feature by the the observation points)? Is there an optimal spacing of observations that minimises the sensitivity?
  3. Can our qualitative observations from Khan and Kevlahan (2021) be verified by the results of the sensitivity analysis?

We examine these questions one by one, beginning with Question I.

4.1 Question I

We analyse how the sensitivity of the surface wave error to perturbations in the measurements depends on the magnitude of the error of the reconstructed bathymetry.

If the error in the bathymetry reconstruction is defined as

(4.2)
β(t)β(b)L2β(t)L2,

(where β(t) is the true bathymetry and β(b) is the optimal reconstruction), then our objective is to understand whether a relatively higher error for the bathymetry reconstruction (4.2) is correlated with increased sensitivity G/∂m where the response functional G is given by (4.1). This question is significant because high sensitivity G/∂m is undesirable, as it implies that variations in the position and number of observations have a large effect on the accuracy of surface wave predictions

To answer this question, we compare the error in convergence to the exact bathymetry of the optimal reconstruction as seen in Khan and Kevlahan (2021) and the sensitivity G/∂m for the cases considered. These results are summarised in Table 1, and shown in Figure 1.

Table 1

Cases considered for data assimilation algorithm, and comparison of the relative L2 reconstruction error (4.2) in the bathymetry as shown in Figure 4 of Khan and Kevlahan (2021), and the time integrated sensitivity 0TG/mdt of the surface wave error to the observations.


CASE BATHYMETRY INITIAL CONDITIONS ERROR SENSITIVITY

I Gaussian Gaussian 𝒪(10–3) 𝒪(10–9)

II Sandbar Gaussian 𝒪(10–2) 𝒪(10–5)

We recall that these cases were chosen to analyse convergence in scenarios where the support of ϕ(x) and the support of β(t)(x) overlap or are disjoint, and to evaluate the effect of a surface wave with compactly supported initial conditions (Cases I and II) or periodic initial conditions (Case III). We consider Gaussian and sandbar profiles for the bathymetry as a 1-D approximation for peaks and ridges characterising ocean bathymetry. The results of the kappa test (3.1) for each case are given in Figure 1(d). Using these data assimilation results, we implement Algorithm 1 to find the sensitivity of the cost function G to the observations m.

Recall that our analysis here aims to help explain the convergence results presented in Part I. (Khan and Kevlahan, 2021). In Part I the best convergence was for Case I, with a minimum relative bathymetry reconstruction error of 𝒪(10–3) (Figure 4 in Khan and Kevlahan (2021)). For Cases II and III, the optimal results did not accurately recover the exact bathymetry, and the error in the reconstructed convergence error was correspondingly larger, 𝒪(10–2) and 𝒪(10–1) respectively.

To find the optimal ν such that Hν = F as in (3.22), we used the Bi-Conjugate Gradient Stabilised Method (BICGSTAB) to solve the system for ν. Figure 1(e) and (f) illustrate the convergence of Hν to the right hand side F for the optimal ν. The BICGSTAB method is a Krylov linear solver designed for non-symmetric linear operators, and was found to prove better convergence in this analysis than with comparable solvers such as GMRES (Generalised Minimum Residual Method). We observe that the error in convergence in Case I (Figure 1(e)) is due to noise in the right hand side F, and that Hν in Figure 1(f) shows some overfitting. A reasonable tolerance for convergence of BICGSTAB is a relative residual error less than 10–4, however due to the small-scale noise shown in Figure 1(e) and (f), this tolerance was not met. Convergence is improved for BICGSTAB when the system is preconditioned (Barrett et al., 1994). Consequently, errors may be due to the fact that preconditioners for the operator H were not easily computable in the algorithm, and hence were not used. As the error in Figure 1(e) and (f) seems to be due to a small amount of noise, for the purposes of this study we consider the convergence to be satisfactory.

Thus, assuming that we have found the optimal value for ν (representing the adjoint variable P5) given F and the operator H, we can derive the sensitivity G/∂m = δ(x – xj)P3(x, t) as the solution of the system (3.17) using ν as the bathymetry perturbation λ, where δ(x – xj) is the observation operator mapping the state variables onto the observation space.

The sensitivities for Cases I and II are presented in Figure 2(a) and (b) as a function of time at three observation points: 1 (first), 23 (centre) and 45 (last). We computed the sensitivities for Nobs = 45 observation points, as this value produced optimal convergence in Khan and Kevlahan (2021). The observation points are spaced equidistantly with Δx = 0.06 and with the first point at 0.1L. We recall that the sensitivity dG/dm is an Nobs × Nt matrix.

Sensitivity over time at selected positions, and time–integrated sensitivity at all positions
Figure 2 

(a, b) The sensitivity dG/dm as a function of time (with final time t = T), for assimilation results for Case I. There are Nobs = 45 observations, equidistantly spaced with Δx = 0.06 and with the first point at 0.1L. Results show dG/dm at three distinct observation points mj, where j = 1 (first observation), Nobs2 (the median observation), and Nobs (the last observation). (c, d) The time integrated sensitivity 0TG/mdt at each observation point.

Figure 2(a) shows that for Case 1 the sensitivity at each observation point is small and similar in magnitude, 𝒪(10–8). Comparing the time dependence at the different observation points, we notice that the oscillations are largest at the times when the surface wave is observed by the relevant observation point even though the amplitude of the sensitivity remains low. Since the relative convergence error we observed in Part I for the Case I bathymetry reconstruction was small, 𝒪(10–4), our initial hypothesis is that a low error in the reconstructed bathymetry is correlated with low sensitivity of the surface wave error to perturbations in the observations for this particular reconstructed bathymetry.

To further substantiate this observation, we compare the results for Cases I and II in Figure 2(a) and (b). We see that the sensitivity is much higher for Case I, 𝒪(10–4) at the same three observation points. The oscillations indicating the passing wave are present for Case II as well, and have a relatively lower frequency, possibly due to the flatter curvature of the sandbar bathymetry compared to a Gaussian bathymetry, resulting in a more gradual effect on the surface wave. We note that the bathymetry reconstruction error (4.2) for Case II is 𝒪(10–2), an order of magnitude higher than for Case I. Therefore, the increased sensitivity shown in Figure 2(b) for Case II compared to Case I is correlated with a larger reconstruction error in the bathymetry. The fact that the sensitivity is at least three orders of magnitude higher than Case I further strengthens our initial hypothesis that a higher error in bathymetry reconstruction leads to higher sensitivity of the surface wave error to perturbations in the observations.

Figure 2(c) and (d) shows the time integrated sensitivity across all 45 observation points for Cases I and II. As we can see, the hypothesis that large bathymetry reconstruction errors lead to higher sensitivity of the surface wave error to perturbations in the observations is confirmed across the the complete set of observation points.

We summarise the results comparing sensitivity of the surface wave error to perturbations in the observations and relative error in bathymetry reconstruction in Table 1. We now consider Question II, by observing the trend when we perturb influential parameters in the algorithm, such as the placement and number of observation points.

4.2 Question II

We analyse the sensitivity of the surface wave error to perturbations in the in the locations of the of observation points, as well as on the shape of the bathymetry. Our objective is to determine if there is an optimal spacing of the observations that minimises the sensitivity. Results are presented for Cases I and II, where we vary the spacing of observation points, and shape of the bathymetry respectively. The latter consideration is motivated by the significant difference in sensitivity we observed between Case I (Gaussian bathymetry with relatively low standard deviation) and Case II, (a sandbar bathymetry that is less localized).

First, we consider the sensitivity of the surface wave error to the spacing of the observation points, keeping the total number of observation points fixed. Our objective is to find the spacing of a fixed number of uniformly spaced observation points that minimises the sensitivity of the surface wave error to perturbations in the observations. We choose Nobs = 45 for these experiments. We run several iterations of the data assimilation algorithm, where in each iteration the position of the last observation point x45 is fixed, and as Δx is varied, the the first observation point is shifted back (as in Figure 3d, e and f).

The time–integrated sensitivity as location of the first observation point is varied
Figure 3 

Case I: The time integrated sensitivity of the surface wave error 0TG/mdt as the location of the first observation point is varied such that the observation points cover a greater proportion of the domain and the initial conditions support.

We wish to observe how the resulting sensitivity of the surface wave error to the observations changes, and thus answer the question: does increasing the spacing of the observations (and consequently increasing coverage of the domain by observation points) always lead to a decrease in the sensitivity dG/dm? If so, the implication for tsunami models is that, given a fixed number of evenly spaced observation points, it is best to choose a spacing that covers as much of the bathymetry feature as possible.

However, we also need to ensure that the spacing is small enough that the smallest scales in the surface wave resulting from the bathymetry can be observed. To verify this, in Figure 4 we show the spectrum of the surface wave with and without bathymetry (for both Cases I and II). We observe that for both cases, the bathymetry does not introduce scales smaller than those of the initial condition into the surface wave solution. In other words, it is sufficient that the observation points are spaced closely enough to resolve the initial conditions (i.e. one half the minimum effective length scale of the initial conditions).

Shape and spectrum of the surface wave at t= 1.95 given a flat bathymetry, and non–zero bathymetry for Cases I and II respectively
Figure 4 

(a, c) The surface wave at t = 1.95 given a flat bathymetry (red), and non-zero bathymetry (blue). The amplitudes of bathymetry and initial condition are not to scale, however the location is accurately represented. (b, d) Spectrum of the surface wave given a flat bathymetry and non-zero bathymetry for Cases I and II respectively.

Figure 3 shows the time integrated sensitivity 0TG/mdt of the surface wave error for Case I, corresponding to different observation point spacings Δx. The values of Δx and the resulting observation configurations are given in the bottom panels of Figure 3. We observe that for Case I, despite a large variation in the sensitivity between the first and last observation points, the amplitude of the sensitivity dG/dm is small for each choice of Δx, at most 𝒪(10–8). This suggests that for Gaussian bathymetry and initial conditions (as seen in Figure 3(d)–(f)), larger spacing of observations (or more coverage of the bathymetry) does not significantly increase the sensitivity of the error in the surface wave produced by the reconstructed bathymetry, to the observations mj(t).

Results for Case II are shown in Figure 5. For easy comparison across results we give the semi-log graph of the absolute time integrated sensitivity |0TG/m  dt|. Increasing the spacing Δx dramatically improves |0TG/m  dt|. We note that this may be due to either the increased coverage of the bathymetry support by the observation points, or the larger spacing. We observe that configurations of observations for Case II in Figure 5(b) and (c) show larger coverage of the bathymetry feature by observations than in 5(a), where only half of the bathymetry is observed by the right propagating surface wave. Based on this result, we suggest that while increasing the spacing Δx does not have a significant impact on |0TG/m  dt|, increasing the coverage of the observation network such that it covers the support of the bathymetry, results in lower sensitivity of the response functional G to the observations mj(t).

The absolute time integrated sensitivity as observations cover a greater proportion of the domain
Figure 5 

Case II: the absolute time integrated sensitivity of the surface wave error |0TG/m  dt| as the location of the first observation point is varied such that the observation points cover a greater proportion of the domain and the initial conditions support.

There is a clear difference between Cases I and II, where we have the same initial conditions (Gaussian) but different bathymetry. To better understand the reasons for this difference, in our next analysis we run trials of the data assimilation scheme where we begin with a Gaussian, as in Case I, and iteratively increase the standard deviation of the bathymetry such as it becomes closer in shape qualitatively to a sandbar type bathymetry, with flatter curvature. We consider two different positions for the initial condition, given in Figure 6(a)–(f) and 6(g)–(l) respectively. This is because for Case I the bathymetry is to the right of the support of the initial condition, and so we replicate this positioning in Figure 6(a)–(f). However, in Case II the support for both the bathymetry and the initial conditions is centred at x = 0, and so we replicate this for the set of results in Figure 6(g)–(l).

The absolute time integrated sensitivity as the standard deviation of the bathymetry Gaussian is increased
Figure 6 

The absolute time integrated sensitivity |0TG/m  dt| as the standard deviation of the bathymetry Gaussian is increased. (a–f) show results with the initial condition to the right of the bathymetry, like Case I. (g–l) show results with both the initial conditions and bathymetry centred at x = 0, like Case II.

For both initial conditions, it is clear that the flatter the bathymetry shape (i.e. the longer the bathymetry feature), the higher the sensitivity of the response functional to the locations of the observations. This suggests that the position of the bathymetry relative to the initial conditions support does not influence sensitivity, but that extent of the bathymetry does. Each trial was conducted with Nobs = 45 and Δx = 0.1, with the position of the first observation at x1 = –L (complete coverage of the domain by the observations).

Having gained some insight into the influence of spacing and number of points on sensitivity, our ultimate objective is to gauge whether these conclusions regarding the sensitivity of the surface wave error produced by the bathymetry reconstruction to observations confirm and help us understand the main conclusions of Khan and Kevlahan (2021).

4.3 Question III

In Khan and Kevlahan (2021) we presented results for optimally reconstructing bathymetry using surface wave observations in a variational data assimilation algorithm. We concluded that the data assimilation results are improved by (i) increasing the number of observation points, and (ii) maintaining an optimal ratio between the amplitude of the initial conditions and bathymetry relative to the average depth. We also observed that the error in the surface wave produced by the reconstructed bathymetry was orders of magnitude smaller than the error in the bathymetry reconstruction, suggesting low sensitivity of the surface wave to bathymetry reconstruction error.

In terms of the current sensitivity analysis, the question we wish to address is:

Does satisfying the aforementioned conditions (i) and (ii) also result in low sensitivity of the surface wave error (4.1) to perturbations in the observations {mj(t)}, and is this sensitivity also orders of magnitude lower than the reconstruction error in the bathymetry?

We provide an answer by verifying the following conclusions using the sensitivity analysis results.

  • C.i  Does a sub-optimal ratio of the amplitude of the initial conditions to bathymetry (i.e. those that were found to produce non-convergent results in Khan and Kevlahan (2021)) result in greater sensitivity of the surface wave error to observations?
  • C.ii  If the reconstruction error in the bathymetry is high (i.e. at least 10%), does that imply that the sensitivity of the resulting error in the surface wave to observations is also proportionately high?

(C.i) is based on results in Khan and Kevlahan (2021), where we analysed the affect of varying η^/β^ and β^/H on the resulting reconstruction error in the bathymetry. Here represents the amplitude of the initial condition, β^ is the amplitude of the bathymetry, and H is the average sea depth (normalised to H = 1). We observed that convergence was sub-optimal when η^/β^ was greater than 𝒪(10–2) and bathymetry amplitudes were large (over 10% of the average sea depth H).

Since satisfying η^/β^ ≤ 𝒪(10–2) is necessary for convergence of the data assimilation, we do not explicitly violate this condition in our analysis. Instead, we analyse the effect on the sensitivity of the surface wave error to observations, as β^ is varied (and η^ is fixed at 1% of β^). The relationship between the ratio β^/H and the sensitivity G/∂m can be observed in Figures 7 and 8. The figures correspond to Cases I and II respectively. We present the time-integrated sensitivity as the bathymetry height β^ is varied from 0.01 to c, where c is the approximate cut-off value beyond which the surface wave error was larger than 0.1% (as shown in Khan and Kevlahan (2021)). In Figure 7(a) and 7(c), results are shown for Case I with β^/H = 1%, 15%, and 30% respectively. We observe that the sensitivity increases from 10–9 to 10–7 as the bathymetry height increases from 1% of the depth to 30% of the depth.

The absolute time integrated sensitivity as the relative amplitude of the bathymetry is increased for Case I
Figure 7 

Case I: Absolute time integrated sensitivity |0TG/m  dt| as the relative amplitude of the bathymetry is increased.

The absolute time integrated sensitivity as the relative amplitude of the bathymetry is increased for Case 2
Figure 8 

Case II: absolute time integrated sensitivity |0TG/m  dt| as the relative amplitude of the bathymetry is increased.

Similarly, in Figure 8, we consider the sensitivity for Case I with β^/H = 2%, 8%, and 20%. We see that as β^ increases, the sensitivity increases from 𝒪(10–8) by at at least two orders of magnitude for each increasing value of β^, indicating a clear correlation between the normalised height of the bathymetry β^ and the sensitivity of the surface wave error to observations.

To summarise, in each of the cases sensitivity of errors in the free surface wave to perturbations in the measurements increases with bathymetry height β^/H. We note that for Case I especially, this increase occurs when β^/H increases from 1% to 15%. For Case I sensitivity did not vary significantly when β^/H doubled from 15% to 30%. (see Figure 7(b) and Figure 7(c)). Based on these observations, we conclude that a lower relative bathymetry height β^/H decreases the sensitivity of G to changes in the observations {mj(t)}.

Perhaps the most significant question we aim to address with this analysis is whether a large error in the reconstructed bathymetry is associated with a high sensitivity of the error in the surface wave to perturbations in the observations, as stated in (C.ii). In Khan and Kevlahan (2021) we found that the error in the surface wave was orders of magnitude lower than the error in the bathymetry, suggesting low sensitivity of the surface wave to bathymetry reconstruction error. This implies that relatively large tolerance levels for bathymetry reconstruction error in tsunami models may be acceptable, at least in the one-dimensional case. We now wish to make this conclusion more rigorous, and verify whether the SOA sensitivity of the surface wave error to the observations G/∂m also exhibits the same behaviour.

If G/∂m were the same order of magnitude as bathymetry reconstruction error, this would imply that if the bathymetry reconstruction is sub-optimal, the accuracy of surface wave prediction will be sensitive to perturbations to the observations. This is undesirable from a forecasting perspective as it indicates predictions may vary greatly based on small changes in the observations (such as noise or small changes in position). Due to the ill-posedness of inverse problems of this kind, errors in measurements can be amplified and subsequently have a large affect on the surface wave error.

Fortunately, the sensitivity analysis presented here allows us to conclude that, in all cases, considered the sensitivity G/∂m is orders of magnitude smaller than the optimal values for the bathymetry reconstruction error presented in Table 1. This relationship does not change even for the worst sensitivity results (observed for Case II in Figure 8(c), where the sensitivity is 𝒪(10–3). This is still an order of magnitude smaller than the lowest error observed for Case II, which was 𝒪(10–2).

Additionally, we observed that even when the bathymetry reconstruction error is sub-optimal the sensitivity can be relatively far smaller, as demonstrated in Figures 7(a) and 8(a) where the sensitivity is very small, 𝒪(10–9) and 𝒪(10–8) respectively, while the bathymetry reconstruction error is large, 𝒪(10–1). Thus, although we found that the sensitivity of surface wave error to observations increases with the error in the reconstructed bathymetry, it is nevertheless many orders of magnitude smaller.

5 Conclusion

Variational data assimilation has been proposed as a technique for reconstructing ocean bathymetry from observations of surface waves. We have used the second order adjoint technique to quantify the sensitivity of the surface wave to perturbations in the observations (such as their number, spacing and position relative to bathymetry profiles), given approximately reconstructed bathymetry. Our goal is to better understand the sensitivity of the assimilation process to the observations, and also to quantify the acceptable level of bathymetry error for accurate modelling of surface waves.

We first derived the Hessian product Hν given some arbitrary perturbation of the control variable ν, and then demonstrate that deriving the sensitivity G/∂m involves solving the forced equation Hν = F with the right hand side F dependent on the optimal assimilation results. For the present study, we chose the response G to be the relative L2 error in the surface wave, produced by the reconstructed bathymetry.

We numerically implemented the sensitivity algorithm for the bathymetry assimilation case, in order to further investigate and confirm the conclusions about the accuracy of variational bathymetry assimilation from Khan and Kevlahan (2021). We then used this algorithm to answer three questions arising from our bathymetry assimilation algorithm.

Question I asks whether there is a link between the sensitivity G/∂m and the convergence of the reconstructed bathymetry to the exact form. The results for three bathymetry and initial conditions configurations outlined in Table 1 show that a higher error in bathymetry reconstruction is associated with higher sensitivity of the surface wave error to perturbations in the observations.

Question II asks how changing parameters in the data assimilation scheme, such as the spacing of observation points, affects sensitivity of the surface wave. We analysed the sensitivity G/∂m as we varied the observation point placement and number of points. The main conclusion from varying the spacing of the observation points is that covering all of the bathymetry feature with the observation points dramatically decreases sensitivity (provided the spacing is fine enough to resolve all significant scales of surface wave).

To gain further insight into the differences between cases with localized or extended bathymetry features, we conducted trials where we incrementally increased the standard deviation of a Gaussian bathymetry feature, until its curvature qualitatively resembled a sandbar. We observed that the higher standard deviation led to an increase in the time-integrated sensitivity, regardless of the position of the bathymetry relative to the initial condition.

Finally in Question III, we investigated whether our qualitative conclusions from Part I (Khan and Kevlahan, 2021) could be confirmed by the sensitivity analysis. Our first conclusion addressed the relationship between the bathymetry and initial conditions amplitudes β^ and η^, and the average depth H. In Part I we observed that a necessary condition for convergence of the reconstructed bathymetry was that the relative height of the initial conditions compared to the bathymetry height η^/β^ be less than or equal to 𝒪(10–2) when β^/H was larger than 10% of the depth. In the sensitivity analysis we varied the relative bathymetry amplitude β^/H (keeping η^ at 1% of β^, to investigate the effect of the relative bathymetry amplitude on the sensitivity of the surface wave response G to the observations. For both Cases I and II, we saw a general increase in sensitivity with the relative bathymetry height β^/H.

Secondly, we verified the observation from Part I that the error in the surface wave is orders of magnitude lower than the error in the reconstructed bathymetry, suggesting low sensitivity of the surface wave to reconstruction error. We therefore investigated whether the sensitivity G/∂m exhibited the same behaviour, i.e. whether the sensitivity of the surface wave error to observations was orders of magnitude lower than the bathymetry reconstruction error. We clearly observed that in all cases considered, the sensitivity G/∂m was orders of magnitude lower than the optimal values of the reconstruction error presented in Table 1.

In conclusion, the analyses in this chapter confirm the results observed for the data assimilation in Part I (Khan and Kevlahan, 2021). We have shown that the necessary conditions for convergence of the bathymetry reconstruction error, also correspond to low sensitivity of the surface wave error to the observations. By showing that this sensitivity is orders of magnitude lower than the bathymetry reconstruction error, we conclude that even if the reconstruction of the bathymetry is sub-optimal, the forecast of surface wave exhibits low sensitivity to changes in observations. High sensitivity of the surface wave error to observations would imply that predictions of surface waves may vary greatly due to small changes in the observations. As errors in measurements can be amplified and subsequently have a large affect on the surface wave error, the low sensitivity observed in these analyses is encouraging from a forecasting perspective, especially in situations where the relative bathymetry reconstruction error is sub-optimal (larger than 10%).

Improvements to the current analysis could be made by ensuring better convergence of Hν = F, such that the residual error is decreased further. Equivalent results for the initial conditions assimilation may also shed more light on the sensitivity of surface wave propagation to observations.

We note that these results are for an idealised 1-D case. The next step would be verify the conclusions of this analysis for the 2-D data assimilation, and use realistic ocean bathymetry and observation data. In Khan and Kevlahan (2021) we provided a concise overview on the importance of bathymetry for accurate tsunami modelling. Consequently, the results from the current analysis help us to quantify the accuracy required for the reconstructed bathymetry, such that the sensitivity of the the error in the surface wave to the observations remains low. We also gain a better understanding of how large bathymetry features impact the surface wave accuracy, and the effects of the number and placement of observations. In particular, our results suggest that efforts to improve bathymetry mapping should focus on features with very high relative heights (e.g. seamounts). Therefore, the results observed here provide an encouraging first step towards a more realistic implementation for tsunami models, and may serve as a benchmark for future analyses.

The analyses in this work quantify the sensitivity of the surface wave error given the reconstructed bathymetry to observations. They complement the analysis of sensitivity of the surface wave error to parameters in the models (like the bathymetry and initial conditions amplitudes β^ and η^) in Khan and Kevlahan (2021), where we use Global Sensitivity Analysis (GSA) techniques to derive sensitivity indices quantifying the variation in the surface wave error resulting from the respective parameters.

Acknowledgements

We would like to thank Dr. Arthur Vidard and Dr. Eric Blayo at the Laboratoire Jean Kuntzmann for their suggestion that we consider second order adjoint sensitivity analysis to further investigate our results, and for their valuable feedback on this research project.

Funding information

This project was funded by the NSERC Discovery Grant of Dr. Nicholas Kevlahan, and research exchange at the Laboratoire Jean Kuntzmann was facilitated by the McMaster-CNRS Fund travel grant.

Competing interests

The authors have no competing interests to declare.

References

  1. Barrett, R, et al. 1994. Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods. Society for Industrial and Applied Mathematics. DOI: https://doi.org/10.1137/1.9781611971538 

  2. Kevlahan, N, Khan, R and Protas, B. 2019. On the convergence of data assimilation for the one-dimensional shallow water equations with sparse observations. In: Adv Comput Math, 45: 3195–3216. DOI: https://doi.org/10.1007/s10444-019-09733-6 

  3. Khan, RA and Kevlahan, NK-R. 2021. Variational assimilation of surface wave data for bathymetry reconstruction. Part I: algorithm and test cases. DOI: https://doi.org/10.1080/16000870.2021.1976907 

  4. Nakamura, K, et al. 2006. Sequential Data Assimilation: Information Fusion of a Numerical Simulation and Large Scale Observation Data. In: Journal of Universal Computing, 12(6): 608–626. 

  5. Shutyaev, V, Le Dimet, F-X and Parmuzin, E. 2018. Sensitivity analysis with respect to observations in variational data assimilation for parameter estimation. In: Nonlinear Processes in Geophysics, 25(2): 429–439. DOI: https://doi.org/10.5194/npg-25-429-2018 

  6. Wang, Z, Navon, I, Le Dimet, F-X and Zou, X. Jan. 1992. The second order adjoint analysis: Theory and applications. In: Meteorology and Atmospheric Physics, 50: 3–20. DOI: https://doi.org/10.1007/BF01025501 

comments powered by Disqus