# Variational data assimilation for the wave equation in heterogeneous media

Erik Burman, Janosch Preuss, Tim van Beeck

<figure><img src="image_1.png" /></figure>

## To cite this version:

Erik Burman, Janosch Preuss, Tim van Beeck. Variational data assimilation for the wave equation in heterogeneous media. 2025. hal-05444361v2

HAL Id: hal-05444361

https://hal.science/hal-05444361v2

Preprint submitted on 14 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Variational data assimilation for the wave

equation in heterogeneous media

Numerical investigation of stability

Erik Burman<fnref n="1" />, Janosch Preuss<fnref n="2" />, Tim van Beeck<fnref n="3" />

Contributing authors: e.burman@ucl.ac.uk; janosch.preuss@inria.fr;
t.beeck@math.uni-goettingen.de;

Abstract

In recent years, several numerical methods for solving the unique continuation problem for the wave equation in a homogeneous medium with given data on the lateral boundary of the space-time cylinder have been proposed. This problem enjoys Lipschitz stability if the geometric control condition is fulfilled, which allows devising optimally convergent numerical methods. In this article, we investigate whether these results carry over to the case in which the medium exhibits a jump discontinuity. Our numerical experiments suggest a positive answer. However, we also observe that the presence of discontinuities in the medium renders the computations far more demanding than in the homogeneous case.

**Keywords:** Unique continuation, data assimilation, wave equation, discontinuous coefficients, geometric control condition

<footnote n="1">¹Department of Mathematics, University College London, Gower Street,
London, WC1E 6BT, United Kingdom.</footnote>

<footnote n="2">²Inria Project-Team Makutu, Université de Pau et des Pays de l'Adour,
Avenue de l'Université - BP 1155, Pau, 64013, France.</footnote>

<footnote n="3">³Institute for Numerical and Applied Mathematics, University of
Göttingen, Lotzestr. 16-18, Göttingen, 37083, Germany.</footnote>
# 1 Introduction

We consider the discretization of a data assimilation problem for the wave equation in heterogeneous media. In the following, let $\Omega \subset \mathbb{R}^d$, $d \in \{1, 2, 3\}$, be a bounded domain with smooth boundary $\partial\Omega$. We assume that $\Omega$ is split by a smooth $(d-1)$-dimensional interface $\Gamma$ into two disjoint components such that $\Omega \setminus \Gamma = \Omega_1 \cup \Omega_2$ and $\Omega_1 \cap \Omega_2 = \emptyset$.

We consider a scalar wave speed

$$c(x) := \begin{cases} c_1(x), & x \in \Omega_1, \\ c_2(x), & x \in \Omega_2, \end{cases} \quad \text{with } c_i \in C^\infty(\Omega_i), i = 1, 2,$$

which may be discontinuous across the interface $\Gamma$, but is constant in time. For a final time $T > 0$, we consider the space-time domain $Q = (0, T) \times \Omega$ with lateral boundary $\Sigma := (0, T) \times \partial\Omega$, and for any nonempty subset $\omega \subset \Omega$, we set $\omega_T := (0, T) \times \omega$.

Let $\Box_c$ be the wave operator with discontinuous wave-speed defined as

$$\Box_c u := \partial_t^2 u - \operatorname{div}(c^2 \nabla u).$$

In particular, we write $\Box_1$ if $c_i = 1$ for $i = 1, 2$, corresponding to the case of a homogeneous medium with unit wave speed.

Then, we consider the problem: find $u : Q \to \mathbb{R}$ such that

$$(\Box_c u, u) = (0, 0) \quad \text{in } Q \times \Sigma, \qquad (1.1)$$

along with the natural interface conditions

$$[u]_{\Gamma} = 0 \text{ on } \Sigma, \quad [c^2 \nabla u]_{\Gamma} \cdot n_{\Gamma} = 0 \text{ on } \Sigma, \qquad (1.2)$$

where $n_\Gamma$ is the normal vector on $\Gamma$ pointing into $\Omega_+$ and the jump across $\Gamma$ is defined<fnref n="1" /> as $[u]_\Gamma = u_+ - u_-$ with $u_\pm$ being the traces of $u|_{\Omega_\pm}$ on $\Gamma$.

To ensure that the problem is well posed, one usually prescribes initial data

$$(u, \partial_t u)|_{t=0} = (u_0, u_1) \text{ in } \Omega, \qquad \text{(IVP)}$$

<footnote n="1">¹And analogously for vector-valued functions.</footnote>
where $u_0 \in H^1(Q)$ and $u_1 \in L^2(Q)$ are given. For a detailed analysis of the well-posedness of the problem (1.1)-(1.2)+(IVP) in heterogeneous media, we refer to [1].

In this work, we consider a data assimilation problem instead, where the initial data is unknown, and we assume that there are given measurements $u_\omega \in L^2(\omega_T)$ prescribed in a data domain $\omega_T$, with $\omega \subset \Omega$, $\omega \neq \emptyset$:

$$u = u_{\omega} \text{ in } \omega_T. \qquad (\text{DA})$$

In a medium without an interface, the problem (1.1)+(DA) is well-studied; see for example [2–7]. In particular, if the boundary $\partial\Omega$ is strictly convex, the problem is Lipschitz stable provided that the set $\omega_T \subset Q$ fulfills the geometric control condition (GCC) [8]. Roughly speaking, a set $\omega \subset \Omega$ and a final time $T$ fulfill the GCC if all geometric optic rays traveling with unit speed<fnref n="2" /> in $\Omega$, accounting for possible reflections on the boundary $\partial\Omega$, intersect the set $(0, T) \times \omega$. For a more precise definition, we refer to [8]. Assuming that the GCC is satisfied, we have the following result.

**Theorem 1.1.** Assume that $\partial\Omega$ is strictly convex and $c_i = 1$ for $i = 1, 2$. If $\omega_T \subset Q$ fulfills the GCC, then there exists a constant $C > 0$ such that for any $\phi \in H^1(Q)$, we have the following estimates:

$$\begin{aligned} \|\phi\|_{L^{\infty}(0,T;L^2(\Omega))} + \|\partial_t\phi\|_{L^2(0,T;H^{-1}(\Omega))} &\le C (\|\phi\|_{L^2(\omega_T)} + \|\phi\|_{L^2(\Sigma)} + \|\Box_1\phi\|_{H^{-1}(Q)}), \\ \|\phi|_{t=0}\|_{L^2(\Omega)} + \|\partial_t\phi|_{t=0}\|_{H^{-1}(\Omega)} &\le C (\|\phi\|_{L^2(\omega_T)} + \|\phi\|_{L^2(\Sigma)} + \|\Box_1\phi\|_{H^{-1}(Q)}). \end{aligned}$$

*Proof.* See [9, Thm. A.4] and [9, Rem. A.5]. □

In the heterogeneous case, where $c_1 \neq c_2$, the problem is less studied. In the absence of the GCC, we have to assume that the final time $T$ is large enough to ensure that problem (1.1)-(1.2)+(DA) is well-posed. To be precise, Thm. 5.19 of [10] (restated below as Theorem 1.2) guarantees that solutions of (1.1)-(1.2)+(DA) are unique if

$$T > 2 \sup_{x \in \Omega} \mathrm{dist}_c(x, \omega), \qquad (1.3)$$

<footnote n="2">²We consider the metric to be adapted to the wave speed c.</footnote>
where $\mathrm{dist}_c(x, \omega) := \inf_{y \in \omega} d_c(x, y)$. Here, $d_c$ is the distance function defined as the length of the shortest continuous path in $\Omega$ connecting $x$ and $y$ with respect to the metric adapted to the discontinuous wave speed $c$. We remark that the condition on the time in (1.3) is identical to the one assumed in the classical unique continuation result of Tataru for a smooth coefficient $c$. While the condition (1.3) guarantees unique solvability, the stability of problem (1.1)-(1.2)+(DA) can be poor. However, as shown in [10, Thm. 5.19], Assumption (1.3) grants the following stability result holds.

**Theorem 1.2** (Thm. 5.19 of [10]). Let $\omega \subset \Omega$ be nonempty and let $T > 0$ be s.t. (1.3) holds. Then, there exists $C, \kappa$ such that for any initial data $(u_0, u_1) \in H_0^1(\Omega) \times L^2(\Omega)$ and $u$ solving (1.2)+(IVP) and

$$ \Box_c u \in L^2(Q), \quad u = 0 \text{ on } \Sigma, $$

one has for any $\mu > 0$ that

$$ \|(u_0, u_1)\|_{L^2 \times H^{-1}} \le C e^{\kappa \mu} \left( \|u\|_{L^2(\omega_T)} + \|\Box_c u\|_{L^2(Q)} \right) + \frac{C}{\mu} \|(u_0, u_1)\|_{H^1 \times L^2}. \quad (1.4) $$

In particular, if $(u_0, u_1) \ne (0, 0)$ we have that

$$ \|(u_0, u_1)\|_{L^2 \times H^{-1}} \le C \frac{\|(u_0, u_1)\|_{H^1 \times L^2}}{\log \left( 1 + \frac{\|(u_0, u_1)\|_{H^1 \times L^2}}{\|u\|_{L^2(\omega_T)} + \|\Box_c u\|_{L^2(Q)}} \right)}. \quad (1.5) $$

In this article, we consider a discontinuous-in-time finite element method to solve the data assimilation problem (1.1)-(1.2)+(DA) in the heterogeneous case. This method was originally introduced for the homogeneous case in [11]. The focus of this article is not on the numerical analysis of the proposed discretization, but rather on investigating the numerically observed stability under the assumption of the GCC. Specifically, we want to investigate numerically whether the following assumption – motivated by Theorem 1.1 for the homogeneous case – is justified.

**Assumption 1.** If $\omega_T \subset Q$ fulfills the GCC, then there exists a constant $C > 0$ such that for any $\phi \in H^1(Q)$ we have the following estimates

$$ \|\phi\|_{L^\infty(0,T;L^2(\Omega))} + \|\partial_t \phi\|_{L^2(0,T;H^{-1}(\Omega))} \le C \left( \|\phi\|_{L^2(\omega_T)} + \|\phi\|_{L^2(\Sigma)} + \|\Box_c \phi\|_{H^{-1}(Q)} \right), $$
$$ \|\phi|_{t=0}\|_{L^2(\Omega)} + \|\partial_t\phi|_{t=0}\|_{H^{-1}(\Omega)} \le C (\|\phi\|_{L^2(\omega_T)} + \|\phi\|_{L^2(\Sigma)} + \|\Box_c\phi\|_{H^{-1}(Q)}). $$

In summary, our goal is to understand whether the Lipschitz stability property is preserved in the presence of strong heterogeneity in the medium, provided that the GCC holds. This question is of high practical importance, as discontinuous media appear in many applications, for example, in seismology. If the corresponding data assimilation problems were indeed subject to the poor logarithmic stability property ensured by (1.5), this would pose an enormous challenge for the design of accurate and reliable numerical methods.

**Remark 1.** While Assumption 1 appears natural, it is possible that in certain configurations some additional conditions on the geometry and minimal time T have to be imposed to prove Lipschitz stability. This is due to the complex interactions of rays with the interface $\Gamma$. For further discussion we refer to the reference [12] in which sufficient conditions for boundary controllability of wave equations with a jumping coefficient have been derived and in particular to Remark 1.6 therein.

# 2 Fully discrete discontinuous-in-time discretization

In this section, we present a fully discrete discontinuous-in-time finite element method for discretizing the data assimilation problem (1.1)-(1.2)+(DA), extending the method from [11] to the heterogeneous case. Section 2.1 describes the geometric partition of the space-time domain $Q$ into time slabs, followed by the full discretization in Section 2.2. The chosen discretization allows for an iterative solution of the linear systems which is described in detail in [11, Section 5].

## 2.1 Discretization of the space-time domain

Let $\hat{T}_h$ be a quasi-uniform triangulation with mesh size $h$ such that $\hat{\Omega} := \cup_{K \in \hat{T}_h} K$ contains $\Omega$, i.e. $\Omega \subset \hat{\Omega}$. Setting $\mathcal{T}_h := \{K \cap \Omega | K \in \hat{T}_h\}$ then ensures that $\Omega = \cup_{K \in \mathcal{T}_h} K$ holds true. We recall from [9, Sec. 4.2] that the triangulation $\{\hat{T}_h | h > 0\}$ can be constructed such that for $h$ small enough the following continuous trace inequality
holds: For all $v \in [H^1(K)]^d$ and $K \in \mathcal{T}_h$, we have

$$ \|v\|_{[L^2(\partial K)]^d} \le C (h^{-1/2} \|v\|_{[L^2(K)]^d} + h^{1/2} \|\nabla v\|_{[L^2(K)]^d}). \quad (2.1) $$

We assume additionally that the triangulation $\mathcal{T}_h$ fits the interface $\Gamma$ and the data domain $\omega$. An extension to allow for unfitted interfaces could be possible by utilizing techniques introduced in [13] for a stationary problem, but this is certainly beyond the scope of this article.

For polynomial degree $k \ge 1$, we define the $H^1$-conforming finite element space

$$ V_h^k := \{v \in H^1(\Omega) : v|_K \in \mathcal{P}^k(K) \ \forall K \in \mathcal{T}_h\}, \quad (2.2) $$

where $\mathcal{P}^k(K)$ denotes the space of polynomials of degree at most $k \in \mathbb{N}$ on $K \in \mathcal{T}_h$. Furthermore, we partition the time axis into $N$ subintervals $I_n = (t_n, t_{n+1})$, $n = 0, \dots, N-1$, where $0 = t_0 \le t_1 \le \dots \le t_N = T$. We assume that the intervals are of equal length and denote $\Delta t = |t_{n+1} - t_n|$. Then, we partition $Q$ and $\Sigma$ into time-slabs

$$ \begin{gathered} Q^n := I_n \times \Omega, \quad \Sigma^n := I_n \times \Sigma, \quad n = 0, \dots, N-1, \\ Q = \bigcup_{n=0}^{N-1} Q^n, \quad \Sigma = \bigcup_{n=0}^{N-1} \Sigma^n. \end{gathered} \quad (2.3) $$

In the following, we denote the space-time integrals on the time slabs as

$$ (u, v)_{Q^n} := \int_{I_n} \int_{\Omega} uv \, dxdt, \quad (u, v)_{\Sigma^n} := \int_{I_n} \int_{\Sigma} uv \, dSdt, $$

and define $\|v\|_{Q^n}^2 := (v, v)_{Q^n}$ and $\|v\|_{\Sigma^n}^2 := (v, v)_{\Sigma^n}$. Furthermore, we set

$$ \omega^n := I_n \times \omega, \quad (u, v)_{\omega^n} := \int_{I_n} \int_{\omega} uv \, dxdt, \quad \|v\|_{\omega^T}^2 := \sum_{n=0}^{N-1} (v, v)_{\omega^n}. $$

Finally, we define the time jump operator

$$ v_{\pm}^{n}(x) := \lim_{s \to 0^{+}} v(x, t_{n} \pm s), \quad \llbracket v^{n} \rrbracket := v_{+}^{n} - v_{-}^{n}. $$
## 2.2 Full discretization

We define the discontinuous in time finite element spaces

$$W_{h, \Delta t}^{k,q} := \bigotimes_{n=0}^{N-1} \mathcal{P}^q(I_n) \otimes V_h^k, \quad W_{h, \Delta t}^{k,q} := W_{h, \Delta t}^{k,q} \times W_{h, \Delta t}^{k,q}, \quad q \in \mathbb{N}_0, k \in \mathbb{N}. \quad (2.4)$$

For elements of $W_{h, \Delta t}^{k,q}$, we use the notation $\underline{U}_h = (\underline{u}_1, \underline{u}_2) \in W_{h, \Delta t}^{k,q}$. In the following, we denote

$$a(u, v)_{Q^n} := \int_{I_n} \int_{\Omega} c^2 \nabla u \cdot \nabla v \, dxdt,$$

and introduce a bilinear form $A$ that represents the wave-equation (1.1) in mixed formulation:

$$A[\underline{U}_h, \underline{Y}_h] := \sum_{n=0}^{N-1} \left\{ \begin{aligned} & (\partial_t \underline{u}_2, \underline{y}_1)_{Q^n} + a(\underline{u}_1, \underline{y}_1)_{Q^n} + (\partial_t \underline{u}_1 - \underline{u}_2, \underline{y}_2)_{Q^n} \\ & - (c^2 \nabla \underline{u}_1 \cdot \mathbf{n}, \underline{y}_1)_{\Sigma^n} \end{aligned} \right\}, \quad (2.5)$$

where $\underline{U}_h \in W_{h, \Delta t}^{k,q}$ and $\underline{Y}_h \in W_{h, \Delta t}^{k^*, q^*}$ with $k, k^*, q \in \mathbb{N}$ and $q^* \in \mathbb{N}_0$. To approximate solutions of (1.1)-(1.2)+(DA), we search for stationary points of the Lagrangian

$$\mathcal{L}_h(\underline{U}_h, \underline{Z}_h) := \frac{1}{2} \|\underline{u}_1 - u_\omega\|_{\omega_T}^2 + A[\underline{U}_h, \underline{Z}_h] + \frac{1}{2} S_h(\underline{U}_h, \underline{U}_h) \\ - \frac{1}{2} S_h^*(\underline{Z}_h, \underline{Z}_h) + \frac{1}{2} S_{\Delta t}^{\uparrow\downarrow}(\underline{U}_h, \underline{U}_h),$$

for $\underline{U}_h \in W_{h, \Delta t}^{k,q}$, $\underline{Z}_h \in W_{h, \Delta t}^{k^*, q^*}$, where $S_h, S_h^*$, and $S_{\Delta t}^{\uparrow\downarrow}$ are stabilization terms yet to be defined. Note that the first and the second terms of $\mathcal{L}_h$ incorporate the data and the PDE constraints, respectively. To define the stabilization terms $S_h$ and $S_h^*$, we introduce the following terms:

$$J(\underline{U}_h, \underline{W}_h) := \sum_{n=0}^{N-1} \int_{I_n} \sum_{F \in \mathcal{F}_i} h(\llbracket c^2 \nabla \underline{u}_1 \rrbracket \cdot \mathbf{n}, \llbracket c^2 \nabla \underline{w}_1 \rrbracket \cdot \mathbf{n})_F \, dt,$$

$$G(\underline{U}_h, \underline{W}_h) := \sum_{n=0}^{N-1} \int_{I_n} \sum_{K \in \mathcal{T}_h} h^2 (\partial_t \underline{u}_2 - \operatorname{div}(c^2 \nabla \underline{u}_1), \partial_t \underline{w}_2 - \operatorname{div}(c^2 \nabla \underline{w}_1))_K \, dt,$$

$$I_0(\underline{U}_h, \underline{W}_h) := \sum_{n=0}^{N-1} (\underline{u}_2 - \partial_t \underline{u}_1, \underline{w}_2 - \partial_t \underline{w}_1)_{Q_n}, \quad R(\underline{U}_h, \underline{W}_h) := \sum_{n=0}^{N-1} h^{-1}(\underline{u}_1, \underline{w}_1)_{\Sigma^n},$$
where $[c^2\nabla]\cdot \mathbf{n}$ denotes the jump of the normal derivative over interior facets $\mathcal{F}_i$. $J(\cdot, \cdot)$ is a continuous interior penalty term in space, $G(\cdot, \cdot)$ is a Galerkin least squares term enforcing the PDE locally on each element, $I_0(\cdot, \cdot)$ enforces that $u_2 = \partial_t u_1$, and $R(\cdot, \cdot)$ ensures boundary stability. Then we define the primal stabilizer $S_h$ as

$$S_h(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) := J(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) + I_0(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) + G(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) + R(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h), \quad (2.6)$$

and the dual stabilizer $S_h^*$ through

$$S_h^*(\underline{\mathbf{Y}}_h, \underline{\mathbf{Z}}_h) := \sum_{n=0}^{N-1} \left\{ (\underline{y}_1, \underline{z}_1)_{Q^n} + a(\underline{y}_1, \underline{z}_1)_{Q^n} + (\underline{y}_2, \underline{z}_2)_{Q^n} + h^{-1}(\underline{y}_1, \underline{z}_1)_{\Sigma^n} \right\}. \quad (2.7)$$

The remaining stabilization term $S_{\Delta t}^{\uparrow\downarrow}$ imposes regularity on the discontinuities in time and is defined as

$$S_{\Delta t}^{\uparrow\downarrow}(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) := \underline{I}_1^{\uparrow\downarrow}(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) + \underline{I}_2^{\uparrow\downarrow}(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h), \quad (2.8)$$

where

$$\underline{I}_1^{\uparrow\downarrow}(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) := \sum_{n=1}^{N-1} \left\{ \frac{1}{\Delta t} (\llbracket \underline{u}_1^n \rrbracket, \llbracket \underline{w}_1^n \rrbracket)_{\Omega} + \Delta t (c^2 \llbracket \nabla \underline{u}_1^n \rrbracket, c^2 \llbracket \nabla \underline{w}_1^n \rrbracket)_{\Omega} \right\}, \\ \underline{I}_2^{\uparrow\downarrow}(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) := \sum_{n=1}^{N-1} \frac{1}{\Delta t} (\llbracket \underline{u}_2^n \rrbracket, \llbracket \underline{w}_2^n \rrbracket)_{\Omega}.$$

With the definition of $\mathcal{L}_h$, the first order optimality conditions take the form: Find $(\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) \in \mathcal{W}_{h, \Delta t}^{k,q} \times \mathcal{W}_{h, \Delta t}^{k^*, q^*}$ such that

$$(\underline{u}_1, \underline{w}_1)_{\omega_T} + A[\underline{\mathbf{W}}_h, \underline{\mathbf{Z}}_h] + S_h(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) + S_{\Delta t}^{\uparrow\downarrow}(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) = (u_\omega, \underline{w}_1)_{\omega_T} \quad \forall \underline{\mathbf{W}}_h \in \mathcal{W}_{h, \Delta t}^{k,q}, \\ A[\underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h] - S_h^*(\underline{\mathbf{Y}}_h, \underline{\mathbf{Z}}_h) = 0 \quad \forall \underline{\mathbf{Y}}_h \in \mathcal{W}_{h, \Delta t}^{k^*, q^*}.$$

In a more compact form, we can write these conditions as: Find $(\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) \in \mathcal{W}_{h, \Delta t}^{k,q} \times \mathcal{W}_{h, \Delta t}^{k^*, q^*}$ such that

$$B[(\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h), (\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h)] = (u_\omega, \underline{w}_1)_{\omega_T} \quad \forall (\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h) \in \mathcal{W}_{h, \Delta t}^{k,q} \times \mathcal{W}_{h, \Delta t}^{k^*, q^*}, \quad (2.9)$$
where

$$ B[(\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h), (\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h)] := (\underline{u}_1, \underline{w}_1)_{\omega_T} + A[\underline{\mathbf{W}}_h, \underline{\mathbf{Z}}_h] + S_h(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) \\ + S_{\Delta t}^{\uparrow\downarrow}(\underline{\mathbf{U}}_h, \underline{\mathbf{W}}_h) + A[\underline{\mathbf{U}}_h, \underline{\mathbf{Y}}_h] - S_h^*(\underline{\mathbf{Y}}_h, \underline{\mathbf{Z}}_h). \quad (2.10) $$

# 3 Error analysis

This section outlines how the steps from [11] can be adapted to the heterogeneous setting to analyze the approximation of (1.1)-(1.2)+(DA) through (2.9). In Section 3.1, we establish inf-sup stability and continuity of the discrete problem with respect to suitable discrete norms. In Section 3.2, we show that the error between the first component $\underline{u}_1$ of the discrete solution and the continuous solution $u$ is bounded by the best approximation error in these norms. Finally, we argue that Assumption 1 allows to recover convergence rates in a physically relevant norm.

## 3.1 Inf-sup stability and continuity

In the following, we assume that (1.3) is satisfied. Let $|\cdot|_{S_h}$, $|\cdot|_{\uparrow\downarrow}$, and $\|\cdot\|_{S_h^*}$ be the (semi-)norms induced by the stabilizers $S_h$, $S_{\Delta t}^{\uparrow\downarrow}$, and $S_h^*$:

$$ |\underline{\mathbf{U}}_h|_{S_h}^2 := S_h(\underline{\mathbf{U}}_h, \underline{\mathbf{U}}_h), \quad |\underline{\mathbf{U}}_h|_{\uparrow\downarrow}^2 := S_{\Delta t}^{\uparrow\downarrow}(\underline{\mathbf{U}}_h, \underline{\mathbf{U}}_h), \quad \|\underline{\mathbf{Z}}_h\|_{S_h^*}^2 := S_h^*(\underline{\mathbf{Z}}_h, \underline{\mathbf{Z}}_h). \quad (3.1) $$

Then, we define the discrete norm

$$ \|\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h\|^2 := |\underline{\mathbf{U}}_h|_{S_h}^2 + |\underline{\mathbf{U}}_h|_{\uparrow\downarrow}^2 + \|\underline{u}_1\|_{\omega_T}^2 + \|\underline{\mathbf{Z}}_h\|_{S_h^*}^2, \quad (3.2) $$

and its strengthened version

$$ \|\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h\|_{\square_c}^2 := \|\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h\|^2 + \sum_{n=0}^{N-1} \left\{ \begin{aligned} & \|\partial_t \underline{u}_2\|_{Q^n}^2 + \|c^2 \nabla \underline{u}_1\|_{Q^n}^2 + \|\partial_t \underline{u}_1\|_{Q^n}^2 \\ & + \|\underline{u}_2\|_{Q^n}^2 + \int_{I_n} \sum_{K \in \mathcal{T}_h} h^2 \|c^2 \underline{u}_1\|_{H^2(K)}^2 dt \end{aligned} \right\}. \quad (3.3) $$

**Lemma 3.1.** The expressions $\|\cdot, \cdot\|$ and $\|\cdot, \cdot\|_{\square_c}$ are norms on $W_{h, \Delta t}^{k,q}$ and $W_{h, \Delta t}^{k^*, q^*}$.

*Proof.* It suffices to show that $\|\cdot, \cdot\|$ is a norm. Since $\|\cdot, \cdot\|$ is a semi-norm by definition, we only require that $\|\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h\| = 0$ implies $\underline{\mathbf{U}}_h = \underline{\mathbf{Z}}_h = 0$. Assume that
|||(U_h, Z_h)||| = 0. Then, in particular, |U_h|_↑↓ = 0 and therefore (u_1, u_2) ∈ [H^1(Q)]^2.
Furthermore, by definition of the stabilization terms S_h and S_h^*, we have that Z_h = 0,
u_1|_Σ = 0, u_1|_{ω_T} = 0 and ∂_t u_1 = u_2. Similar to [11, Lemma 2], partial integration yields

$$
\begin{align*}
\|\Box_c \underline{u}_1\|_{H^{-1}(Q)}^2 &:= \sup_{\substack{y \in H_0^1(Q), \\ \|y\|_{H^1(Q)}=1}} \int_Q \{-(\partial_t \underline{u}_1)\partial_t y + c^2 \nabla \underline{u}_1 \nabla y\} \\
&= \sup_{\substack{y \in H_0^1(Q), \\ \|y\|_{H^1(Q)}=1}} \left\{ \sum_{n=0}^{N-1} (\underline{u}_2 - \partial_t \underline{u}_1, \partial_t y)_{Q_n} + \sum_{n=0}^{N-1} \int_{I_n} \sum_{K \in \mathcal{T}_h} (\partial_t \underline{u}_2 - \operatorname{div}(c^2 \nabla \underline{u}_1), y)_K \, dt \right. \\
&\qquad \left. + \sum_{n=0}^{N-1} \int_{I_n} \sum_{F \in \mathcal{F}_i} (\llbracket c^2 \nabla \underline{u}_1 \rrbracket \cdot \mathbf{n}, y)_F \, dt + \sum_{n=0}^{N-1} (\llbracket \underline{u}_1^n \rrbracket, y)_{\Sigma^n} \right\},
\end{align*}
$$

It follows by definition of the stabilizers and |||(U_h, Z_h)||| = 0 that ||\Box_c u_1||_{H^{-1}(Q)} = 0 and that the interface conditions (1.2) are fulfilled. Thus, by taking $\mu \to \infty$ in Theorem 1.2, we conclude that $u_0 = u_1 = 0$. As solutions to (1.1)-(1.2)+(IVP) depend continuously on the data [1], it follows that $\underline{U}_h = 0$. □

From this result and the identity

$$
B[(\underline{U}_h, \underline{Z}_h), (\underline{U}_h, -\underline{Z}_h)] = |||(\underline{U}_h, \underline{Z}_h)|||^2 = |||(\underline{U}_h, \underline{Z}_h)|||||(\underline{U}_h, -\underline{Z}_h)||,
$$

it directly follows that the bilinear form $B$ enjoys inf-sup stability on $W_{h, \Delta t}^{k,q} \times W_{h, \Delta t}^{k^*, q^*}$ with respect to the $||(\cdot, \cdot)||$-norm.

**Corollary 3.2.** *There exists a constant $C_B > 0$ such that*

$$
\sup_{(\underline{W}_h, \underline{Y}_h) \in \mathcal{W}_{h, \Delta t}^{k,q} \times \mathcal{W}_{h, \Delta t}^{k^*, q^*}} \frac{B[(\underline{U}_h, \underline{Z}_h), (\underline{W}_h, \underline{Y}_h)]}{|||(\underline{W}_h, \underline{Y}_h)|||} \geq C_B |||(\underline{U}_h, \underline{Z}_h)||| \quad (3.4)
$$

for all $(\underline{U}_h, \underline{Z}_h) \in \mathcal{W}_{h, \Delta t}^{k,q} \times \mathcal{W}_{h, \Delta t}^{k^*, q^*}$.

The trace inequality (2.1) and the Cauchy-Schwarz inequality immediately yield
the continuity of the bilinear form A.

**Lemma 3.3** (Continuity of A). If $\mathbf{U}|_{Q^n} \in [H^1(Q^n) \cap L^2(0, T; H^2(\mathcal{T}_h))] \times H^1(Q^n)$ for $n = 0, \dots, N-1$ and for all $\mathbf{Y}_h \in \mathcal{W}_{h, \Delta t}^{k^*, q^*}$, we have that

$$
A[\mathbf{U}, \mathbf{Y}_h] \leq C |||(\mathbf{U}, 0)|||_{\Box_c} |||(0, \mathbf{Y}_h)|||.
$$
Proof. First, note that $||(0, \underline{\mathbf{Y}}_h)||| = ||\underline{\mathbf{Y}}_h|||_{S_h^*}$. With the Cauchy-Schwarz inequality, we obtain

$$ \sum_{n=0}^{N-1} \left\{ (\partial_t u_2, \underline{y}_1)_{Q^n} + a(u_1, \underline{y}_1)_{Q^n} \right\} \le C \left( \sum_{n=0}^{N-1} \|\partial_t u_2\|_{Q^n}^2 + \|c^2 \nabla u_1\|_{Q^n}^2 \right)^{1/2} \|\underline{\mathbf{Y}}_h\|_{S_h^*} $$

and

$$ \sum_{n=0}^{N-1} (\partial_t u_1 - u_2, \underline{y}_2)_{Q^n} \le C \left( \sum_{n=0}^{N-1} \|\partial_t u_1\|_{Q^n}^2 + \|u_2\|_{Q^n}^2 \right)^{1/2} \|\underline{\mathbf{Y}}_h\|_{S_h^*} $$

Since the interface $\Gamma$ is fitted by the triangulation, we have that $c^2 \nabla u_1 \in H^1(K)$ for each $K \in \mathcal{T}_h$. Thus, we can apply the trace inequality (2.1) to obtain

$$ \sum_{n=0}^{N-1} (c^2 \nabla u_1 \cdot n, \underline{y}_1)_{\Sigma^n} \\ \le \left( \sum_{n=0}^{N-1} \int_{I_n} \sum_{K \in \mathcal{T}_h} h \|c^2 \nabla v\|_{\partial \Omega \cap \partial K}^2 dt \right)^{1/2} \left( \sum_{n=0}^{N-1} \int_{I_n} \sum_{K \in \mathcal{T}_h} h^{-1} \|\underline{y}_1\|_{\partial \Omega \cap \partial K}^2 dt \right)^{1/2} \\ \le C \left( \sum_{n=0}^{N-1} \left\{ \|c^2 \nabla u_1\|_{Q^n}^2 + \int_{I_n} \sum_{K \in \mathcal{T}_h} h^2 \|c^2 u_1\|_{H^2(K)} dt \right\} \right)^{1/2} \|\underline{\mathbf{Y}}_h\|_{S_h^*} $$

The claim follows from putting all estimates together. □

## 3.2 Approximation results

With Corollary 3.2 and Lemma 3.3, we can show that the approximation error in the $||(\cdot, \cdot)|||$-norm is bounded by the best approximation error in the $||(\cdot, \cdot)|||_{\square_c}$-norm.

**Theorem 3.4.** Let $u$ be a sufficient regular solution of (1.1)-(1.2)+(DA) and $(\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h) \in W_{h, \Delta t}^{k,q} \times W_{h, \Delta t}^{k^*, q^*}$ be the solution to (2.9). Set $\mathbf{U} := (u, \partial_t u)$. Then, there exists a constant $C > 0$ such that

$$ |||(\mathbf{U} - \underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h)||| \le \left(1 + \frac{C}{C_B}\right) \inf_{\underline{\mathbf{Y}}_h \in W_{h, \Delta t}^{k,q}} |||(\mathbf{U} - \underline{\mathbf{V}}_h, 0)|||_{\square_c}. \quad (3.5) $$

Proof. Let $\underline{\mathbf{V}}_h \in W_{h, \Delta t}^{k,q}$ and $(\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h) \in W_{h, \Delta t}^{k,q} \times W_{h, \Delta t}^{k^*, q^*}$ be arbitrary. The triangle inequality yields

$$ |||(\mathbf{U} - \underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h)||| \le |||(\mathbf{U} - \underline{\mathbf{V}}_h, 0)||| + |||(\underline{\mathbf{V}}_h - \underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h)||| $$
For the second term, we consider

$$
\begin{align*}
B[(\underline{\mathbf{U}}_h - \underline{\mathbf{V}}_h, \underline{\mathbf{Z}}_h), (\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h)]
&= B[(\underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h), (\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h)] - (\underline{v}_1, \underline{w}_1)_{\omega_T} - S_h(\underline{\mathbf{V}}_h, \underline{\mathbf{W}}_h) - S_{\Delta t}^{\uparrow\downarrow}(\underline{\mathbf{V}}_h, \underline{\mathbf{W}}_h) - A[\underline{\mathbf{V}}_h, \underline{\mathbf{Y}}_h] \\
&= (u_\omega, \underline{w}_1)_{\omega_T} - (\underline{v}_1, \underline{w}_1)_{\omega_T} - S_h(\underline{\mathbf{V}}_h, \underline{\mathbf{W}}_h) - S_{\Delta t}^{\uparrow\downarrow}(\underline{\mathbf{V}}_h, \underline{\mathbf{W}}_h) - A[\underline{\mathbf{V}}_h, \underline{\mathbf{Y}}_h] \\
&= (u - \underline{v}_1, \underline{w}_1)_{\omega_T} + S_h(\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, \underline{\mathbf{W}}_h) + S_{\Delta t}^{\uparrow\downarrow}(\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, \underline{\mathbf{W}}_h) + A[\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, \underline{\mathbf{Y}}_h],
\end{align*}
$$

where we use the fact that $\mathbf{U}$ is sufficiently smooth and $u = u_\omega$ on $\omega_T$. Then, we can apply Lemma 3.3 to obtain

$$
\begin{align*}
(u - \underline{v}_1, \underline{w}_1)_{\omega_T} + S_h(\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, \underline{\mathbf{W}}_h) + S_{\Delta t}^{\uparrow\downarrow}(\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, \underline{\mathbf{W}}_h) + A[\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, \underline{\mathbf{Y}}_h] \\
&\le C(\|\|(\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, 0)\|\|_{\square_c} \|\|(0, \underline{\mathbf{Y}}_h)\| + \|\|(\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, 0)\|\|\|(\underline{\mathbf{W}}_h, 0)\|) \\
&\le C\|\|(\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, 0)\|\|_{\square_c} \|\|(\underline{\mathbf{W}}_h, \underline{\mathbf{Y}}_h)\|\|.
\end{align*}
$$

Thus, the inf-sup condition on $B$ yields that

$$
\|\|(\underline{\mathbf{U}} - \underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h)\|\| \le \frac{C}{C_B} \|\|(\underline{\mathbf{U}} - \underline{\mathbf{V}}_h, 0)\|\|_{\square_c},
$$

which gives the claim.

<figure><img src="image_3.png" /></figure>

**Corollary 3.5.** Under the assumption of Theorem 3.4 and assuming additionally that $\Delta t = C'h$ for some $C' > 0$ we have that

$$
\|\|(\underline{\mathbf{U}} - \underline{\mathbf{U}}_h, \underline{\mathbf{Z}}_h)\|\| \le C h^s \|u\|_{H^{s+2}(Q)},
$$

where $s := \min\{k, q\}$.

*Proof.* Due to the assumption that $c_i \in C^\infty(\Omega_i)$, $i = 1, 2$, there exists $C > 0$ s.t.

$$
\|\|(\underline{\mathbf{U}}_h - \underline{\mathbf{V}}_h, 0)\|\|_{\square_c} \le C \|\|(\underline{\mathbf{U}}_h - \underline{\mathbf{V}}_h, 0)\|\|_{\square_1},
$$

where $\|\|(\cdot, \cdot)\|_{\square_1}$ is defined by (3.3) with $c_i = 1, i = 1, 2$. Thus, we can apply the interpolation results derived in [11, Lemma 5 & 11] to obtain the claim. $\square$

<figure><img src="image_4.png" /></figure>
If the problem is Lipschitz stable, for instance in the homogeneous case (see Theorem 1.1), we can prove convergence in a physically relevant norm. To deal with the discontinuity in time of the discrete solution, a suitable lifting operator $L_{\Delta t}: W_{h,\Delta t}^{k,q} \to C^0(0,T;V_h^k)$ that ensures that $L_{\Delta t}u_1 \in H^1(Q)$ has been introduced in [11, Sec. 4]. Then, we state the following result.

**Theorem 3.6.** Let Assumption 1 be satisfied and assume that $\Delta t = C'h$. For $u$ being a sufficiently regular solution of (1.1) and $(\underline{U}_h, \underline{Z}_h) \in W_{h,\Delta t}^{k,q} \times W_{h,\Delta t}^{k^*,q^*}$ being the solution to (2.9), there exists a constant $C > 0$ such that

$$ \|u - L_{\Delta t} \underline{u}_1\|_{L^{\infty}(0,T;L^2(\Omega))} + \|\partial_t(u - L_{\Delta t} \underline{u}_1)\|_{L^2(0,T;H^{-1}(\Omega))} \leq C h^s \|u\|_{H^{s+2}(Q)}, \quad (3.6) $$

where $s := \min\{k,q\}$. In particular, we have that

$$ \|u_0 - \underline{u}_1|_{t=0}\|_{L^2(\Omega)} + \|u_1 - (\partial_t \underline{u}_1)|_{t=0}\|_{H^{-1}(\Omega)} \leq C h^s \|u\|_{H^{s+2}(Q)}. \quad (3.7) $$

*Proof.* The proof follows with the same arguments as the proof of [11, Theorem 14]. □

**Remark 2 (Homogenous media).** If $c_i = 1$ for $i = 1,2$, Assumption 1 holds (see Theorem 1.1), and we recover the result from [11, Theorem 14].

**Remark 3 (Logarithmic stability).** If Assumption 1 is not satisfied, we could still hope to utilize the results from Theorem 1.2, provided that $T$ is large enough such that (1.3) holds. To this end, we would have to deal with the following technicalities:

(1) The results from Theorem 1.2 are formulated in terms of $\|\Box_c u\|_{L^2(Q)}$, but we can only control $\|\Box_c \tilde{e}_h\|_{H^{-1}(Q)}$ for the residual $\tilde{e}_h = u - L_{\Delta t} \underline{u}_1$ to the right order. Thus, we would have to rewrite the estimates from Theorem 1.2 with $\|\Box_c u\|_{H^{-1}(Q)}$ instead of $\|\Box_c u\|_{L^2(Q)}$. To this end, one might be able to shift the regularity suitably, cf. [1, Sec. 2.6], or use similar arguments as in [4, Thm. A.1].

(2) We have to ensure that the constants in the estimates from Theorem 1.2 remain bounded when we consider the residual $\tilde{e}_h$. In particular, we have to control $C_{0,1} := \|(\tilde{e}_h|_{t=0}, \partial_t \tilde{e}_h|_{t=0})\|_{H^1 \times L^2}$. We note that we can always ensure that $C_{0,1}$ is bounded by adding a suitable Tikhonov regularization term to (2.6).
Then, we can use the arguments of [11, Theorem 14] to show that

$$ \|(u - \underline{u}_1)|_{t=0}\|_{L^2(\Omega)} + \|\partial_t(u - \underline{u}_1)|_{t=0}\|_{H^{-1}(\Omega)} \le CC_{0,1} \log(1 + C_{0,1}h^{-s}\|u\|_{H^{s+2}(Q)}^{-1})^{-1}. $$

# 4 Numerical experiments

In this section, we present numerical experiments carried out with the proposed method. The examples in Section 4.1-Section 4.2 are implemented using the FEniCSx library [14–17]. The example in Section 4.3 uses the finite element software Netgen/NGSolve [18, 19] and in addition the space-time functionality provided by the add-on ngsxfem [20]. A reference to reproduction material can be found in the section on Code availability.

## 4.1 Examples in one space dimension

We consider the unit interval $\Omega := (0,1) \subset \mathbb{R}$ partitioned into $2^{L+1}$, $L \in \mathbb{N}$, uniform elements of size $h = 1/2^{L+1}$. Unless specified otherwise, we choose the time step as $\Delta t = h$. Further, we set $q = k$ and choose $k^* = 1$, $q^* = 1$ if $k = 1$ and $q^* = 0$ if $k > 1$.

### 4.1.1 A first example

We define two subdomains $\Omega_1 = (0, 0.5)$ and $\Omega_2 = (0.5, 1.0)$ such that $\bar{\Omega} = \bar{\Omega}_1 \cup \bar{\Omega}_2$ and fix the wavespeed $c_2 := c|_{\Omega_2} = 1$ in $\Omega_2$. For variable $c_1 := c|_{\Omega_1}$, we consider the exact solution [21]

$$ u(x,t) := \begin{cases} \cos(w_1 c_1 t) \cos(w_1(x - 0.5)), & x \in \Omega_1, \\ \cos(w_2 c_2 t) \cos(w_2(x - 0.5)), & x \in \Omega_2, \end{cases} \quad (4.1) $$

where we set $w_1 = 3\pi$ and $w_2 = w_1 c_1/c_2$. Furthermore, we define

$$ \omega := [0.0, 0.25] \cup [0.75, 1.0], \quad \text{and} \quad \Omega_R := \bar{\Omega} \setminus \omega = [0.25, 0.75]. \quad (4.2) $$

In one space dimension, the requirement (1.3) that $T > 2 \sup_{x \in \Omega} \text{dist}_c(x, \omega_T)$ already implies the GCC. Thus, we can focus for the one-dimensional experiments on
the condition (1.3). Let $d(x, y) := |x - y|$, $x, y \in \Omega$, be the standard Euclidean metric.
As in [10], we define a metric scaled with the discontinuous wave speed through

$$d_c(x, y) := 1_{\Omega_1} c_1(x)^{-1} d(x, y) + 1_{\Omega_2} d(x, y), \quad x, y \in \Omega. \quad (4.3)$$

Thus, for $\omega$ as above, $c_2 = 1$ fixed and $c_1 > 1.0$, the condition (1.3) translates to $T > 0.25(1 + c_1^{-1})$. First, we consider the contrast $c_1 = 2.5$, which means that the GCC is satisfied if $T > 0.35$.

In Figure 1, we observe that the errors $\|u - L_{\Delta t}u_1\|_{L^{\infty}(0,T;L^2(\Omega_R))}$ and $\|\partial_t(u - L_{\Delta t}u_1)\|_{L^2(0,T;L^2(\Omega_R))}$ converge with optimal order $O(h^k)$ for $T = 0.5 > 0.35$. However, if $T = 0.1 < 0.35$, the order of convergence degrades, while the respective $L^2$-best approximation<fnref n="3" /> errors still converge with optimal order. This suggests that the GCC seems to imply Lipschitz stability even in the heterogeneous case, and therefore Assumption 1 is indeed justified.

To investigate the dependence of the errors on the contrast, we fix the polynomial degree $k = 3$, refinement level $L = 3$, and $c_2 = 1$, and consider increasing values of $c_1 \in \{1.0, 1.5, \dots, 4.5\}$. In Fig. 2, we observe that the approximation errors increase polynomially with $O(c_1^3)$, while the $L^2$-best approximation errors grow quadratically with $O(c_1^2)$. This is an important observation since an exponential growth in the contrast would have rendered the corresponding problems nearly impossible to resolve numerically.

However, this also entails that a strong refinement of the discretization is required to maintain the accuracy of the numerical solution as the contrast increases. We illustrate this problem in Fig. 3. For $c_1 = 2.5$, the exact solution is already well-captured by the numerical solution on refinement level $L = 2$. However, for higher contrast $c_1 = 5.5$ this refinement level is clearly insufficient.

### 4.1.2 Multiple jumps

To increase the complexity of the problem, we modify the exact solution (4.1) to allow for multiple jumps in the wave speed $c$. For two points $p_1, p_2 \in (0, 1)$, we decompose $\Omega$

<footnote n="3">³We compare with the $L^2$-best approximation error, where we solve for the $L^2$-orthogonal projection of the exact solution onto the discrete space, as a measure of the approximation quality of the discrete space.</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 1:** For the approximation of (4.1) with polynomial degree $k \in \{2, 3\}$ and contrast $c_1 = 2.5$ we measure the errors $\|u - L_{\Delta t}u_1\|_{L^{\infty}(0,T;L^2(\Omega_R))}$ (left) and $\|\partial_t(u - L_{\Delta t}u_1)\|_{L^2(0,T;L^2(\Omega_R))}$ (right) and compare to the respective error of the $L^2$-best approximation (dashed). For final time $T = 0.5 > 0.35$ (top), we observe optimal convergence rates, while for $T = 0.1 < 0.35$ (bottom) the order of convergence degrades.</figcaption></figure>

into three subdomains $\Omega_1 = (0, p_1)$, $\Omega_2 = (p_1, p_2)$, and $\Omega_3 = (p_2, 1)$ s.t. $\bar{\Omega} = \cup_{i=1,2,3} \bar{\Omega}_i$.
<figure><img src="image_6.png" /><figcaption>**Fig. 2:** For $k = 3$ and $L = 3$, we consider the errors $\|\partial_t(u - L_{\Delta t}u_1)\|_{L^2(0,T;L^2(\Omega_R))}$ (left) and $\|u - L_{\Delta t}u_1\|_{L^\infty(0,T;L^2(\Omega_R))}$ (right) and the respective $L^2$-best approximation errors (dashed) for increasing contrast $c_1 \in \{1.0, 1.5, \dots, 4.5\}$. We compare the case where $T = 0.5$ is fixed for all choices of $c_1$ and the case where $T$ adapted such that it barely fulfills $T > 0.25(1 + c_1^{-1})$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 3:** At time $t = 0.25$, we compare the exact solution (4.1) and approximations with polynomial degree $k = 3$ on refinement levels $L \in \{1, 2, 3, 4\}$ for contrasts $c_1 = 2.5$ (left) and $c_1 = 5.5$ (right). The scale in the right plot indicates the mesh size $h = 1/2^{L+1}$ for the different refinement levels.</figcaption></figure>
Then, we consider the following ansatz for the exact solution:

$$ u(x, t) := \begin{cases} \cos(w_1 c_1 t) \cos(w_1 (x - p_1)), & x \in \Omega_1, \\ \cos(w_2 c_2 t) \cos(w_2 (x - p_1)), & x \in \Omega_2, \\ \cos(w_3 c_3 t) \cos(w_3 (x - p_2)), & x \in \Omega_3. \end{cases} \qquad (4.4) $$

To ensure the continuity of (4.4), we choose $p_2$ in dependence of $p_1$ and $w_2$. For $n \in \mathbb{Z}$, we set

$$ p_2 := \frac{2\pi n + w_2 p_1}{w_2} \implies \cos(w_2(p_2 - p_1)) = 1. $$

We fix $c_2 = 1.0$ and set $w_1 = 3\pi$, $w_2 = w_1 c_1/c_2, c_3 = c_1$, and $w_3 = w_1$.

As a first example, we consider the data domain $\omega = [0.0, 0.3] \times [0.7, 1.0]$ and choose $p_1 = 0.4$ and $n = 3$ such that $p_2 \approx 0.667$. Thus, the interval $(p_1, p_2)$ is not contained in the data domain.

For contrast $c_1 = c_3 \in \{7.5, 11.5\}$, Fig. 4 shows the exact solution and the approximation with $k = 3$ for fixed time step $\Delta t = 1/32$ on the refinement levels $L \in \{3, 4\}$. In the case of contrast $c_1 = 11.5$, the refinement level $L = 3$ does not seem to be sufficient to resolve the solution accurately, but for $L = 4$ the solution is captured reasonably well, even though we do not prescribe data in the interval $(p_1, p_2)$.

As a second example, we consider the case where the data is only given in one half of the interval, i.e. $\omega = [0, 0.3]$. The requirement (1.3) and thus the GCC is fulfilled if

$$ T > 2 \left[ (p_2 - p_1) + \frac{p_1 - 0.3 + 1 - p_2}{c_1} \right]. \qquad (4.5) $$

We consider $c_1 \in \{2.5, 7.5\}$, for which the requirement $T > \{0.89, 0.65\}$ ensures the GCC. Figure 5 compares the results for polynomial degree $k = 3$ on the refinement level $L = 4$ with a fixed time step $\Delta t = 1/32$, for two final times $T = 0.5$ and $T = 1.0$. We compare the approximation of the exact solution for both final times at time $t = 0.5$. For both contrasts, the approximate solution matches the exact solution well if the GCC is fulfilled ($T = 1.0$), but shows significant deviations when it is not ($T = 0.5$). The deviation is more extreme for contrast $c_1 = 2.5$, where the gap between the considered time $T = 0.5$ and the threshold $T > 0.89$ is larger.
<figure><img src="image_8.png" /><figcaption>**Fig. 4:** Exact solution (4.4) and approximated solution with $k = 3$ of (4.4) at $t = 0.25$ with contrast $c_1 = 7.5$ (left) and $c_1 = 11.5$ (right). We compare the approximations on the refinement levels $L = 3$ and $L = 4$, and set $\Delta t = 1/32$ for both cases.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5:** For fixed time step $\Delta t = 1/32$, refinement level $L = 4$, and polynomial degree $k = 3$, we compare the exact solution (4.4) and approximated solution computed with final time $T = 1.0$ and $T = 0.5$ at time $t = 0.5$. We set $c_2 = 1.0$ and consider the cases $c_1 = c_3 = 2.5$, $n = 1$ (left) and $c_1 = c_3 = 7.5$, $n = 3$ (right).</figcaption></figure>
## 4.2 Example in two space dimensions

We consider another example in two space dimensions. We partition the unit-square $\Omega := [0, 1]^2 \subset \mathbb{R}^2$ into subdomains $\Omega_1 = \{(x, y) \in \Omega : x < 0.5\}$ and $\Omega_2 = \{(x, y) \in \Omega : x > 0.5\}$, and extend the exact solution defined in (4.1) to the two-dimensional case by setting $u(x, y; t) := u(x; t)$.

We define the data domain $\omega := \Omega \setminus [0.25, 0.75]^2$ and note that $\omega$ fulfills the GCC. To ensure that the data domain and the interface $\Gamma$ are meshed exactly, we consider a quadrilateral mesh with $2^{L+1}$, $L \in \{1, 2, 3, 4\}$, elements in each direction. We consider the final time $T = 0.75$ and set the time step size to $\Delta t = (1/2)^{L+1}$. As before, we fix the wave speed $c_2 = 1$ and consider different values for $c_1$ to control the contrast. Further, we choose the parameters $q = k^* = q^* = k$.

The results of the experiment are presented in Fig. 6. With contrast $c_1 = 2.5$, we observe quasi-optimal convergence with order $O(h^k)$ for polynomial degrees $k \in \{2, 3\}$ for the errors $\|\partial_t(u - L_{\Delta t}u_1)\|_{L^2(0,T;L^2(\Omega))}$ and $\|u - L_{\Delta t}u_1\|_{L^\infty(0,T;L^2(\Omega))}$. Similar as already observed in the one-dimensional case in Fig. 2, we observe that these errors are more sensitive to increasing contrast $c_1$ than the corresponding $L^2$-best approximation errors. In particular, the ratio for $\|u - L_{\Delta t}u_1\|_{L^\infty(0,T;L^2(\Omega))}$ with $k=2$ appears to grow as $O(c_1^{3/2})$, which renders it very demanding to maintain a reasonable accuracy as the contrast increases.

## 4.3 Whispering gallery modes in an annulus

As a final example, we would like to consider a two-dimensional setting in which the Lipschitz stability estimates from Assumption 1 are not valid. To this end, we follow reference [22], in which it has been shown that, under certain assumptions on the contrast of the sound speed and the geometry, so-called whispering gallery modes exist. These are solutions of the wave equation that are exponentially decaying away from the interface, showing that the logarithmic stability estimate (1.5) from [10] is in general optimal.

Indeed, in certain domains $\Omega$, we can follow [22] to construct a sequence of eigenfunctions $w_\ell$ of $-\Delta_c := -\operatorname{div}(c^2\nabla\cdot)$ with corresponding eigenvalues $\lambda_\ell \to \infty$, such
<figure><img src="image_10.png" /><figcaption>Fig. 6: On the left: We observe quasi-optimal convergence of $||\partial_t(u - L_{\Delta t}u_1)||_{L^2(0,T;L^2(\Omega))}$ (solid) and $||u - L_{\Delta t}u_1||_{L^\infty(0,T;L^2(\Omega))}$ (dashed) with contrast $c_1 = 2.5$ and polynomial degrees $k \in \{2,3\}$. On the right: We consider the ratio between these errors and the $L^2$-best approximation error for increasing contrast $c_1$.</figcaption></figure>

that for all $\omega \subset \Omega$ with $\text{dist}(\bar{\omega}, \Gamma) > 0$, there exists constants $C, D > 0$ such that

$$ \|w_{\ell}\|_{L^2(\omega)} \le C \exp(-D\sqrt{\lambda_{\ell}}), \quad -\Delta_c w_{\ell} = \lambda_{\ell} w_{\ell}, \quad w_{\ell}|_{\partial\Omega} = 0, \quad \|w_{\ell}\|_{L^2(\Omega)} = 1. $$

Then, we obtain a family of solutions of the wave equation

$$ u_{\ell}(t, x) = \cos(\sqrt{\lambda_{\ell}t})w_{\ell}(x), \quad \|u_{\ell}|_{t=0}\|_{L^2(\Omega)} + \|\partial_t u_{\ell}|_{t=0}\|_{H^{-1}(\Omega)} = 1, \quad (4.6) $$

whose initial condition can not be controlled uniformly in a Lipschitz stable manner by the data. This follows, since

$$ \|u_{\ell}\|_{L^2(\omega_T)} + \|u_{\ell}\|_{L^2(\Sigma)} + \|\Box_c u_{\ell}\|_{H^{-1}(Q)} \le C_T \exp(-D\sqrt{\lambda_{\ell}}) $$

goes to zero as $\ell \to \infty$.

As shown in [22, Section 2], one possible geometrical configuration to accommodate such modes is an annulus, where $\Omega_1 = \{R_0 < r < R_1\}$ and $\Omega_2 = \{R_1 < r < R_2\}$ with $R_0 < R_1 < R_2$. To ensure the existence of whispering Gallery modes the contrast has
to fulfill the condition

$$ \frac{c_1}{R_1} < \frac{c_2}{R_2}. $$

Due to the spherical geometry, we can obtain the mode $w_\ell$ by separation of variables

$$ w_{\ell}(r, \theta) = \exp(i\ell\theta)\psi_{\ell}(r) \qquad (4.7) $$

and by solving the one-dimensional eigenvalue problem

$$ -\frac{c^2}{\ell^2} \left[ \partial_r^2 \psi_\ell + \frac{1}{r} \partial_r \psi_\ell \right] + \frac{c^2}{r^2} \psi_\ell = \frac{\lambda_\ell}{\ell^2} \psi_\ell \qquad (4.8) $$

with Dirichlet boundary conditions $\psi_\ell(r = R_0) = \psi_\ell(r = R_2) = 0$.

While there exists many eigenpairs $(\lambda_\ell, \psi_\ell)$ for fixed $\ell$ that solve (4.8), for our purposes it suffices to pick just one. We fix the parameters

$$ R_0 = 3/5, \quad R_1 = 1, \quad R_2 = 3/2, \quad c_1 = 1, \quad c_2 = \sqrt{20} $$

and find two eigenpairs for $\ell = 1$ and $\ell = 10$ whose (radial) eigenfunctions $\psi_\ell$ are plotted in the upper panel of Fig. 7. The corresponding solutions of the wave equation⁴ $u_\ell(t, x)$ as defined in (4.6) are also displayed at the final time $t = T$ for $T = 1$. The corresponding eigenvalues are $\lambda_1 \approx 190.669$ and $\lambda_{10} \approx 404.347$. Since $\lambda_{10} > \lambda_1$, the corresponding eigenfunction decays much faster towards the right of the interface at $r = 1$ as can be seen in the central panel of Fig. 7. The eigenfunctions shown here do not fall off significantly towards the left of the interface since we consider very low $\ell$ in order to keep the computational costs moderate. For higher $\ell \approx 200$, the corresponding eigenfunctions indeed decay rapidly also towards the left of the interface. However, these functions oscillate heavily in $\theta$-direction and in time making them very challenging to resolve numerically. The two modes we have chosen for $\ell = 1$ and $\ell = 10$ are fortunately sufficient for our purpose.

We want to evaluate the performance of our method in reconstructing the modes $u_\ell$, particularly, in the case where the data domain $\omega$ is far away from the interface so that the mode for $\ell = 10$ is very small in $\omega$. For this purpose, we consider the data domain $\omega_o = \{1.3 < r < R_2\}$ in which the mode is very small according to the central

<footnote n="3">³In the numerical experiments we work only with the real part of $w_\ell$ to obtain real-valued solutions</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 7:** The upper panel displays the modes $\psi(r)$ for $\ell = 1$ (left), respectively $\ell = 10$ (right) obtained by solving the eigenvalue problem (4.8). The insets show the corresponding solutions of the wave equation defined by (4.6) and (4.7) at $t = T$. The absolute value of $\psi(r)$ is shown in the central panel. The lower panel displays the convergence of $\|u_\ell - u_1\|_{L^2(Q)}$. The case where the data is given in $\omega_i$ is shown in magenta (solid lines with crosses) whereas the case where data is given in $\omega_o$ is displayed in cyan (same linestyle). For reference, the circles (marked as '+' i.c.' in the legend) display the errors when the method is run with initial conditions $u_\ell(\cdot, t=0)$ and $\partial_t u_\ell(\cdot, t=0)$ known everywhere. The insets in the lower panel display the numerical reconstructions achieved at the final time corresponding to the refinement level indicated by the arrows.</figcaption></figure>
panel of Fig. 7. For comparison, we also consider a data domain close to the interface
$\omega_i = \{0.85 < r < R_1\}$.

The $L^2$-errors in $Q = (0, T) \times \Omega$ for the unique continuation problem of the mode
$u_\ell$ computed with polynomial degree $q = k^* = q^* = k = 2$ are displayed in the lower
panel of Fig. 7. Here, we have chosen $T = 1$ such that the condition (1.3) for unique
solvability is fulfilled. When the data is given in $\omega_i$ (i.e. close to the interface), then
the method converges with the same order as when the initial condition is given, i.e.
the convergence is as good as for a well-posed problem irrespective of the considered
mode $\ell$. However, when the data is given far away from the interface in $\omega_0$, then the
rate of convergence drastically deteriorates as $\ell$ increases. This is in full agreement
with the discussion in the beginning of this subsection: the Lipschitz stability result
of Assumption 1 cannot hold for $\text{dist}(\bar{\omega}_0, \Gamma) > 0$. The stability is of logarithmic type
as asserted by Theorem 1.2 (Thm. 5.19 of [10]) and this is exemplified in our numer-
ical experiments by the observation that the convergence rate appears to drops to
logarithmic type as $\ell$ increases.

# 5 Conclusion

In this article, we investigated a numerical method for solving the variational data
assimilation problem for the wave equation with discontinuous wave speed. Our numer-
ical experiments lead to two main conclusions. First, obtaining accurate numerical
solutions in the presence of material interfaces involving large contrasts is computa-
tionally demanding, and the resolution required to capture the complex behavior of
the solution is a limiting factor in numerical simulations. This challenge is amplified
by the fact that unique continuation problems are globally coupled in space and time.
Second, even in the case of a discontinuous wave speed, the geometric control condi-
tion seems to ensure Lipschitz stability, allowing the design of accurate and reliable
numerical methods.

**Acknowledgments.** The authors would like to thank Spyridon Filippas and Lauri Oksanen for interesting discussions.
## Declarations

On behalf of all authors, the corresponding author states that there is no conflict of interest.

• **Funding** Funding by EPSRC grant EP/V050400/1 is gratefully acknowledged.

• **Code availability** Reproduction material is available at https://github.com/TimvanBeeck/waveUC_discCoefs.

## References

[1] Stolk, C.: On the modeling and inversion of seismic data. PhD thesis, Utrecht University (2000). https://dspace.library.uu.nl/handle/1874/855

[2] Cindea, N., Münch, A.: Inverse problems for linear hyperbolic equations using mixed formulations. Inverse Problems **31**(7), 075001 (2015) https://doi.org/10.1088/0266-5611/31/7/075001

[3] Burman, E., Feizmohammadi, A., Oksanen, L.: A finite element data assimilation method for the wave equation. Math. Comp. **89**(324), 1681–1709 (2020) https://doi.org/10.1090/mcom/3508

[4] Burman, Erik, Feizmohammadi, Ali, Münch, Arnaud, Oksanen, Lauri: Space time stabilized finite element methods for a unique continuation problem subject to the wave equation. ESAIM Math. Model. Numer. Anal. **55**, 969–991 (2021) https://doi.org/10.1051/m2an/2020062

[5] Mishra, S., Molinaro, R.: Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for PDEs. IMA J. Numer. Anal. **42**(2), 981–1022 (2021) https://doi.org/10.1093/imanum/drab032

[6] Dahmen, Wolfgang, Monsuur, Harald, Stevenson, Rob: Least squares solvers for ill-posed pdes that are conditionally stable. ESAIM Math. Model. Numer. Anal. **57**(4), 2227–2255 (2023) https://doi.org/10.1051/m2an/2023050
[7] Burman, E., Delay, G., Ern, A.: The unique continuation problem for the wave equation discretized with a high-order space-time nonconforming method. Numer. Math. **157**, 1259–1284 (2025) https://doi.org/10.1007/s00211-025-01479-2

[8] Bardos, C., Lebeau, G., Rauch, J.: Sharp sufficient conditions for the observation, control, and stabilization of waves from the boundary. SIAM J. Control Optim. **30**(5), 1024–1065 (1992) https://doi.org/10.1137/0330055

[9] Burman, Erik, Feizmohammadi, Ali, Münch, Arnaud, Oksanen, Lauri: Spacetime finite element methods for control problems subject to the wave equation. ESAIM Control Optim. Calc. Var. **29**, 41 (2023) https://doi.org/10.1051/cocv/2023028

[10] Filippas, S.: Quantitative unique continuation for wave operators with a jump discontinuity across an interface and applications to approximate control. Bull. Soc. Math. France **152**(2), 199–294 (2024) https://doi.org/10.24033/bsmf.2889

[11] Burman, E., Preuss, J.: Unique continuation for the wave equation based on a discontinuous Galerkin time discretization. IMA J. Numer. Anal., 036 (2025) https://doi.org/10.1093/imanum/draf036

[12] Gagnon, L.: Sufficient conditions for the controllability of wave equations with a transmission condition at the interface. ESAIM: COCV **29**, 51 (2023) https://doi.org/10.1051/cocv/2023034

[13] Burman, E., Preuss, J.: Unique continuation for an elliptic interface problem using unfitted isoparametric finite elements. SMAI J. Comput. Math. **11**, 165–202 (2025) https://doi.org/10.5802/smai-jcm.122

[14] Baratta, I.A., Dean, J.P., Dokken, J.S., Habera, M., Hale, J.S., Richardson, C.N., Rognes, M.E., Scroggs, M.W., Sime, N., Wells, G.N.: DOLFINx: the next generation FEniCS problem solving environment. preprint (2023). https://doi.org/10.5281/zenodo.10447666

[15] Scroggs, M.W., Baratta, I.A., Richardson, C.N., Wells, G.N.: Basix: a runtime
finite element basis evaluation library. J. Open Source Softw. **7**(73), 3982 (2022)
https://doi.org/10.21105/joss.03982

[16] Scroggs, M.W., Dokken, J.S., Richardson, C.N., Wells, G.N.: Construction of arbitrary order finite element degree-of-freedom maps on polygonal and polyhedral cell meshes. ACM Trans. Math. Software **48**(2), 18–11823 (2022) https://doi.org/10.1145/3524456

[17] Alnaes, M.S., Logg, A., Ølgaard, K.B., Rognes, M.E., Wells, G.N.: Unified form language: A domain-specific language for weak formulations of partial differential equations. ACM Trans. Math. Software **40** (2014) https://doi.org/10.1145/2566630

[18] Schöberl, J.: NETGEN An advancing front 2D/3D-mesh generator based on abstract rules. Comput. Vis. Sci. **1**(1), 41–52 (1997) https://doi.org/10.1007/s007910050004

[19] Schöberl, J.: C++11 implementation of finite elements in NGSolve. Technical report, ASC-2014-30, Institute for Analysis and Scientific Computing (September 2014). http://hdl.handle.net/20.500.12708/28346

[20] Lehrenfeld, C., Heimann, F., Preuß, J., Wahl, H.: **ngsxfem**: Add-on to NGSolve for geometrically unfitted finite element discretizations. J. Open Source Softw. **6**(64), 3237 (2021) https://doi.org/10.21105/joss.03237

[21] Mattsson, K., Ham, F., Iaccarino, G.: Stable and accurate wave-propagation in discontinuous media. J. Comput. Phys. **227**(19), 8753–8767 (2008) https://doi.org/10.1016/j.jcp.2008.06.023

[22] Filippas, S.: Whispering gallery modes for a transmission problem. ESAIM: COCV **29**, 75 (2023) https://doi.org/10.1051/cocv/2023067