# MOR-T L: A Novel Model Order Reduction Method for Parametrized Problems with Application to Seismic Wave Propagation

Julien Besset, Hélène Barucq, Rabia Djellouli, Stefano Frambati

► To cite this version:

Julien Besset, Hélène Barucq, Rabia Djellouli, Stefano Frambati. MOR-T L: A Novel Model Order Reduction Method for Parametrized Problems with Application to Seismic Wave Propagation. RR-9583, Inria Bordeaux - Sud Ouest. 2025. hal-05002365v2

HAL Id: hal-05002365

https://inria.hal.science/hal-05002365v2

Submitted on 18 Aug 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
MOR-T<sub>L</sub>: A Novel Model
Order Reduction Method
for Parametrized
Problems with
Application to Seismic
Wave Propagation

Julien Besset, Hélène Barucq, Rabia Djellouli, Stefano Frambati

RESEARCH
REPORT

Nº 9583

March 2025

Project-Team Makutu
<figure><img src="image_2.png" /></figure>
Inria

# MOR-T<sub>L</sub>: A Novel Model Order Reduction Method for Parametrized Problems with Application to Seismic Wave Propagation

Julien Besset*, Hélène Barucq*, Rabia Djellouli†, Stefano Frambati*

Project-Team Makutu

Research Report n° 9583 — version 2‡ — initial version March 2025 — revised version August 2025 — 30 pages

**Abstract:** This paper presents an efficient strategy for constructing Reduced-Order Model (ROM) bases using Taylor polynomial expansions and Fréchet derivatives with respect to model parameters. The proposed approach enables the construction of ROM bases with minimal additional computational cost. By exploiting Fréchet derivatives – solution to the same problem with distinct right-hand sides – the method introduces a streamlined multiple-right-hand-side (RHS) strategy for ROM bases construction. This approach not only reduces overall computational expenses but also improves accuracy during model parameter updates. Numerical experiments on a two-dimensional wave problem demonstrate significant efficiency gains and enhanced performance, highlighting the potential of the proposed method to advance computational cost-effectiveness, particularly in seismic inversion applications.

**Key-words:** Méthode de Réduction d'Ordre (MOR) – Modèle d'Ordre Réduit (ROM) – Dérivées de Fréchet – Polynôme de Taylor – Elements Finis Spectraux (SEM) – Décomposition QR – Recherche Linéaire

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

<footnote>* Makutu, Inria, TotalEnergies, University of Pau et des Pays de l'Adour, 64000 Pau, France.</footnote>

<footnote>† Department of Mathematics and Interdisciplinary Research Institute for the Sciences (IRIS), California State University, Northridge, CA 91330, USA.</footnote>

<footnote>‡ (a) Revision of the MOR-T<sub>L</sub> formulation in the application to the acoustic case. (b) Change for a relative computational time representation in the different tables.</footnote>
# MOR-T<sub>L</sub> : Une nouvelle méthode de réduction d'ordre pour les problèmes paramétrés avec application à la propagation des ondes sismiques.

**Résumé** : Ce rapport présente MOR-T<sub>L</sub>, une nouvelle méthode de réduction d'ordre basée sur le développement de Taylor et les dérivées de Fréchet par rapport aux paramètres du modèle. Cette approche permet de construire des bases réduites avec un coût computationnel minimal à l'aide d'une décomposition QR tout en assurant une robustesse accrue face aux variations des paramètres du modèle.

L'efficacité de la méthode MOR-T<sub>L</sub> est démontrée sur des problèmes de propagation d'ondes sismiques en deux dimensions. Nos expériences numériques montrent que cette méthode permet de réduire significativement le coût de calcul lié à la construction de la base réduite par rapport à une approche classique comme la décomposition orthogonale aux valeurs propres, tout en maintenant une précision élevée, y compris pour des variations du paramètre du modèle importantes.

Ces résultats positionnent MOR-T<sub>L</sub> comme une alternative prometteuse pour les simulations intensives et les problèmes inverses, notamment en inversion sismique, où des mises à jour répétées du modèle sont nécessaires.

**Mots-clés** : Model Order Reduction (MOR) – Reduced Order Model (ROM) – Fréchet derivatives – Taylor polynomial representation – Spectral Element Method (SEM) – QR decomposition – Line Search
# Contents

<table><tr><td>1</td><td>Introduction and Motivation</td><td>4</td></tr><tr><td>2</td><td>Preliminaries: Notations, Assumptions, and Problem Setup</td><td>5</td></tr><tr><td>2.1</td><td>Nomenclature and Assumptions</td><td>5</td></tr><tr><td>2.2</td><td>Mathematical Model</td><td>6</td></tr><tr><td>3</td><td>Introducing MOR-T<sub>L</sub>: A Novel Reduced Order Method for Solving (2.3) Across Parameter Variations</td><td>7</td></tr><tr><td>4</td><td>Application to Two-Dimensional Seismic Wave Problems</td><td>11</td></tr><tr><td>4.1</td><td>The problem prototype</td><td>11</td></tr><tr><td>4.2</td><td>The formulation in a Bounded Domain</td><td>12</td></tr><tr><td>4.3</td><td>The basis functions</td><td>13</td></tr><tr><td>4.4</td><td>The resulting reduced order model</td><td>14</td></tr><tr><td>5</td><td>Assessing cost-effectiveness and performances<sup>1</sup> of MOR-T<sub>L</sub> v.s. SEM in wave problems</td><td>15</td></tr><tr><td>5.1</td><td>Assessing the Cost-Effectiveness of QR and SVD Basis Representations</td><td>16</td></tr><tr><td>5.1.1</td><td>Case of Homogeneous Medium</td><td>16</td></tr><tr><td>5.1.2</td><td>Case of Heterogeneous Medium</td><td>19</td></tr><tr><td>5.2</td><td>Comparative Analysis of MOR-T<sub>L</sub> and SEM in Line Search Minimization</td><td>20</td></tr><tr><td>5.2.1</td><td>Accuracy of the MOR-T<sub>L</sub> in Computing Velocity Perturbations</td><td>21</td></tr><tr><td>5.2.2</td><td>Computational Efficiency and Robustness of the MOR-T<sub>L</sub> During Line Search</td><td>24</td></tr><tr><td>6</td><td>Summary and Conclusion</td><td>28</td></tr></table>

<footnote n="1">¹All these numerical experiments have been performed on a system equipped with an 11th Gen Intel Core i5-1135G7 (4 cores, 8 threads, 2.40GHz), Mesa Intel Xe Graphics (TGL GT2). Python 3.10.12 with Numpy 1.21.5 multi-threaded for BLAS operations</footnote>
# 1 Introduction and Motivation

The development of efficient solution methodologies for mathematical models that describe complex, real-world phenomena is of paramount importance across a wide range of fields, including heat transfer in materials with varying thermal properties, pollutant transfer in fluids, semiconductor applications, seismic and geophysical problems, radar and sonar technologies, and climate dynamics. These phenomena are typically governed by Partial Differential Equations (PDEs) whose discretization results in large, high-dimensional systems that are computationally expensive to solve. Moreover, these models are often highly sensitive to variations in multiple parameters, such as material properties, environmental conditions, or geometrical factors, and there is frequently a need to solve these problems repeatedly for different parameter values. This is especially true in inverse problems, where one seeks to reconstruct a set of parameters (or all of them) from experimental measurements, requiring the repeated solution of direct problems. As the number of parameters increases, the computational burden becomes prohibitive, making this process practically unattainable.

Over the past fifty years, significant efforts have been dedicated to the design of numerical strategies that strike a balance between computational cost and accuracy. Approaches such as adaptive grids [13, 35], where mesh sizes are adjusted based on parameter variations, and variable approximation orders [14, 2], where the order of approximation changes in response to parameter changes, have been employed to reduce the computational load. Additionally, methods like Trefftz [21, 11] and hybridizable discontinuous Galerkin (HDG) methods [9, 8, 12, 22] have been introduced, which use fewer degrees of freedom by exploiting local solutions or static condensation techniques. Although these methods offer computational advantages, their performance remains unsatisfactory in situations where model parameters change frequently, especially in the context of model updates where accuracy must be preserved while also avoiding dramatic increases in computational cost.

A promising recent development in this domain is the use of Model Order Reduction (MOR) techniques, which aim to approximate the high-dimensional solution space of a system using a significantly smaller subspace, thereby reducing the size of the system and improving computational efficiency. Well-known MOR techniques include Proper Orthogonal Decomposition (POD) [10], Dynamic Mode Decomposition (DMD) [32], and neural network-based ROM [27]. However, the major challenge with these approaches is their sensitivity to parameter changes. As parameters vary, the accuracy of the ROM-based solutions can degrade significantly. The simplest ways to address this issue involve reconstructing or enlarging the basis functions, both of which result in increased computational costs and undermine the very efficiency gains that ROM aims to achieve. Alternative methods, such as Reduced Basis (RB) [5] or Greedy algorithms [6], Proper Generalized Decomposition (PGD) [7] and Grassmann Interpolation [26], have been proposed. These so called parametrized-MOR methods offer a powerful and adaptive framework for constructing ROMs that balance accuracy and computational efficiency. However, despite their theoretical appeal and success in many MOR applications, their use in large-scale optimization problems seems to be limited, especially when limited prior information are available. A key challenge lies in defining a representative and sufficiently rich sampling of the parameter space in advance, which is often hindered by the inherent uncertainty and spatial variability of the model parameters. Furthermore, even if a comprehensive sampling strategy could be established, the computational demands of applying these methods to large-scale problems often become prohibitively high.

In response to these challenges, we propose a novel MOR approach based on Taylor polynomial expansions. The key innovation of this method lies in how the ROM basis functions are constructed:
we use Taylor polynomial representations to generate these bases, which leads to a more robust
method that can handle parameter variations more effectively. We have named this new approach
**MOR-T**<sub>*L*</sub> where *T* stands for Taylor polynomials and *L* represents the polynomial degree. The
proposed MOR-T<sub>*L*</sub> strategy provides a unique combination of advantages: it reduces the size of the
system, thereby lowering computational cost, while also being robust to parameter changes. This
dual benefit is critical not only for maintaining accuracy when the model is updated but also for
preventing the computational cost from rising dramatically as parameters vary.

The theoretical foundation and numerical implementation of MOR-T<sub>L</sub> are presented in this paper. We demonstrate its efficiency in reducing computational time while maintaining high accuracy, even as parameters change. By leveraging the Fréchet derivatives of the solution with respect to the problem parameters - solutions to the same problem with distinct right-hand sides - we enable a streamlined multiple-right-hand-side (RHS) strategy for constructing the ROM bases. This method allows the adaptive construction of reduced bases with minimal additional computational cost. As we show in Section 5, MOR-T<sub>L</sub> is highly effective in addressing parameter changes, making it a promising candidate for solving large-scale problems in a variety of applications.

The remainder of this paper is organized as follows. In Section 2, we introduce the adopted nota-
tions, assumptions, and the class of parametrized partial differential equations under consideration.
Section 3 presents the theoretical foundation of the MOR-T<sub>L</sub> method and discusses its computa-
tional complexities and cost. In Section 4, we apply MOR-T<sub>L</sub> to a parametrized two-dimensional
wave problem arising in seismic applications. Finally, Section 5 provides a numerical performance
comparison of MOR-T<sub>L</sub> with the standard POD approach and the Spectral Element Method (SEM),
highlighting the significant computational savings and robustness of MOR-T<sub>L</sub>, which make it a
strong candidate for large-scale seismic applications.

**2 Preliminaries: Notations, Assumptions, and Problem Setup**

In this section, we present the adopted notations, assumptions, and the class of boundary value
problems under consideration.

**2.1 Nomenclature and Assumptions**

Throughout this paper, we adopt the following notations and hypothesis:

• x represents the spatial variable. x = (x₁, x₂) ∈ R².

• *t* denotes a time variable. *t* ≥ 0.

• $N_x$ is a natural number representing the number of degrees of freedom associated with the discretization scheme.

• $N_t$ is a natural number representing the number of time steps.

• $\Delta t$ is a positive number representing the time step.

• $\mathbf{x}_j, j = 1, \dots, N_x$.

• $t_j = j\Delta t, j = 1, \dots, N_t$.

• T is the timelength of the computations. T = N<sub>t</sub>Δt.
* $\Omega$ is a bounded domain in $\mathbb{R}^2$ whose boundary $\partial\Omega$ is Lipschitz continuous [17].

* $X, Y$, and $Z$ are respectively Banach spaces, corresponding to the space of solutions, source terms and model parameters, respectively.

* $\theta$ represents the model parameter. $\theta \in Z$.

* $f$ is a source term. $f \in Y$.

* $u_\theta$ is the solution of the considered boundary value problem with a parameter $\theta$. $u_\theta \in X$.

* $\delta\theta$ represents the model parameter perturbation direction. $\delta\theta \in Z$.

* The derivative of $u_\theta$ at $\theta = \theta_0$ in the direction of $\delta\theta_0$ is given by:

$$Du_{\theta_0} \delta\theta_0 = \lim_{\epsilon \to 0} \frac{1}{\epsilon} (u_{\theta_0 + \epsilon \delta\theta_0} - u_{\theta_0}) \quad (2.1)$$

Hence, we denote the Fréchet derivative of order $\ell$ of $u_\theta$ at $\theta = \theta_0$ in the direction of $\delta\theta_0$ by $v_{\theta_0}^\ell = D^\ell u_{\theta_0} \delta\theta_0^\ell$, and $v_{\theta_0}^0 = u_{\theta_0}$.

* The Taylor polynomial of degree $L$ of $u_\theta$ centered at $\theta_0$ with a perturbation $\delta\theta_0$ is given by:

$$P_{\theta_0}(u_\theta) = \sum_{\ell=0}^{L} \frac{1}{\ell!} v_{\theta_0}^{\ell} \quad (2.2)$$

* $\mathcal{L}_\theta$ is a linear operator from $X$ to $Y$. $\mathcal{L}_\theta$ depends on the parameter $\theta$.

* The derivative of order $\ell$ of $\mathcal{L}_\theta$ at $\theta = \theta_0$ in the direction of $\delta\theta$ is denoted by $D^\ell \mathcal{L}_{\theta_0} \delta\theta_0^\ell$, and $D^0 \mathcal{L}_{\theta_0} \delta\theta_0^0 = \mathcal{L}_{\theta_0}$.

* $\delta_{ij}$ defines the Kronecker symbol, i.e. $\delta_{ij} = 1$ if $i = j$, and 0 otherwise.

* $\delta(\mathbf{x} - \mathbf{x}_0)$ is the delta Dirac function at $\mathbf{x} = \mathbf{x}_0$.

## 2.2 Mathematical Model

Given $f \in Y$ and for $\theta \in Z$, we consider the following class of linear PDE systems:

$$\left\{ \begin{array}{l} \text{Find } u_\theta \in X \text{ such that} \\ \mathcal{L}_\theta u_\theta = f \end{array} \right. \quad (2.3)$$

The goal is to develop an efficient computational procedure to solve (2.3) for a range of values of the parameter $\theta$ in a bounded domain. This is particularly important for applications like inverse problems, including inverse seismic scenarios, where solving (2.3) for varying $\theta$ is often required. It is essential to note that classical representations of the solution $u_{\theta_0}$ using POD for instance, only applies for the fixed value of $\theta = \theta_0$. This limitation arises because the basis functions depend on $\theta_0$, and thus this representation cannot be used to approximate accurately solution to (2.3) for $\theta \neq \theta_0$. As noted in [3], even a small perturbation in $\theta_0$ can lead to significant inaccuracies. Therefore, standard MOR approaches at this stage are inadequate for solving (2.3) for multiple values of $\theta$. To overcome this challenge, based on Taylor polynomial expansion, we propose to construct a basis for the Fréchet derivatives of $u_\theta$ evaluated at $\theta = \theta_0$ in the direction of perturbation $\delta\theta = \delta\theta_0$, then reconstruct the solution $u_\theta$ evaluated in $\theta_\alpha = \theta_0 + \alpha\delta\theta_0$ with $\alpha \in \mathbb{R}$.
## 3 Introducing MOR-T$_L$: A Novel Reduced Order Method for Solving (2.3) Across Parameter Variations

The method we propose for solving (2.3) for a set of admissible parameters $\theta$ falls within the category of MOR methods. What sets our approach apart is the novel way in which we construct the basis functions. Specifically, we use Taylor polynomial representations to build this basis. For this reason, we have named the proposed method MOR-T$_L$: T stands for Taylor polynomials, and L represents the polynomial degree. The strength of the MOR-T$_L$ strategy lies in its ability to solve (2.3) for different values of $\theta$ without sacrificing accuracy. In contrast to existing ROM methods, where the constructed basis for a specific $\theta$ is not suitable for solving (2.3) when $\theta_\alpha = \theta_0 + \alpha\delta\theta_0$ (even for small perturbations $\alpha\delta\theta_0$) our approach retains its accuracy. In traditional MOR methods, such perturbations often result in a significant loss of accuracy [3]. To correct this, the simplest options are to either reconstruct the basis or enlarge it, both of which significantly increase computational costs and undermine the efficiency gains that MOR methods aim to provide.

MOR-T$_L$'s approach is to construct the basis in two distinct stages, outlined below:

### Stage 1: Construction of the Bases for the Fréchet Derivatives Solutions

The first step is to reformulate (2.3), evaluated in $\theta_\alpha = \theta_0 + \alpha\delta\theta_0$, as a series of $L + 1$ PDEs using the Taylor polynomial of order $L$ centered at $\theta_0$ given by (2.2) for $\alpha \in \mathbb{R}$. In the following, we consider the classical Taylor expansion of order $L$ as described by the polynomial defined in (2.2).

$$ \left\{ \begin{array}{l} \text{Find } v_{\theta_0}^{\ell} \in X, \ell = 0, \dots, L, \text{ such that} \\ \mathcal{L}_{\theta_0}(v_{\theta_0}^{\ell}) = f^{(\ell)} \end{array} \right. \qquad (3.1) $$

where $v_{\theta_0}^{\ell} = D^{\ell}u_{\theta_0}\delta\theta_0^{\ell}$, and $f^{(\ell)} = -\sum_{k=0}^{\ell-1} \binom{\ell}{k} (D^{\ell-k}\mathcal{L}_{\theta_0}\delta\theta_0^{\ell-k})(v_{\theta_0}^k)$ with $\binom{\ell}{k} = \frac{\ell!}{k!(\ell-k)!}$; $0 \le k \le \ell-1$.

Note that $\alpha$ drops in (3.1) by linearity of $D$ and $\mathcal{L}_{\theta_0}$:

$$ (D^{\ell-k}\mathcal{L}_{\theta_0}(\alpha\delta\theta_0)^{\ell-k})(\alpha^k v_{\theta_0}^k) = \alpha^\ell (D^{\ell-k}\mathcal{L}_{\theta_0}\delta\theta^{\ell-k})(v_{\theta_0}^k) \qquad (3.2) $$

This property is of great practical interest since it allows the reconstruction of the solution of (2.3) for $\theta = \theta_\alpha$ with different values of $\alpha$ using (2.2) up to multiplication factors.

Also note that $v_{\theta_0}^{\ell}$ are solutions to the PDE (2.3), but with different right-hand sides, resulting in significant computational cost savings.

As with any domain-based discretization scheme, if (3.1) is set in an infinite domain, the problem must first be reformulated in a bounded domain $\Omega$. This is achieved by introducing an artificial external boundary to which a boundary condition is applied. The boundary condition is characterized by the boundary operator $B_\theta$ (see, e.g., [33, 18] and the references therein). Therefore, the boundary value problems considered in the bounded domain are given by:
$$ (\mathbf{BVP}) \quad \left\{ \begin{array}{ll} \text{Find } v_{\theta_0}^{\ell} \in X, \ell = 0, \dots, L & \text{, such that} \\ \mathcal{L}_{\theta_0} (v_{\theta_0}^{\ell}) = f^{(\ell)} & \text{in } \Omega \times (0, T) \\ \mathcal{B}_{\theta_0} (v_{\theta_0}^{\ell}) = 0 & \text{on } \partial\Omega \times (0, T) \end{array} \right. \qquad (3.3) $$

Then, the construction of the reduced bases follows a two-steps procedure.

### (a) Initialization Step.

For $\theta = \theta_0$ fixed, we solve (3.3) using a full order high-fidelity space discretization scheme, e.g., finite element approximation (FEM), that we can represent by the following Galerkin decomposition:

$$ v_{\theta_0}^{\ell}(\mathbf{x}, t) = \sum_{i=1}^{N_x} v_{\theta_0, i}^{\ell}(t) \psi_i(\mathbf{x}); \quad \ell = 0, \dots, L \qquad (3.4) $$

where $\{\psi_i\}_{i=1}^{N_x}$ represents a nodal finite element basis associated with the considered discretization scheme. To simplify the notation, we omit the explicit reference to the dependence of the approximate solution on the domain discretization scheme. Thus, the notation $v_{\theta_0}^{\ell}$ should be understood as the solution to the discretized PDE. Let $\mathbf{v}_{\theta_0}^{\ell}$ be the approximate solution vector. Hence:

$$ \mathbf{v}_{\theta_0}^{\ell}(t) = [v_{\theta_0,1}^{\ell}(t), \dots, v_{\theta_0,N_x}^{\ell}(t)]^T; \quad \ell = 0, \dots, L. \qquad (3.5) $$

Then, for a sequence of times $t_i$, $i = 1, \dots, N_t$, we define the snapshot vector $\mathbf{s}_i^{\ell}$ as:

$$ \mathbf{s}_i^{\ell} = \mathbf{v}_{\theta_0}^{\ell}(t_i), \quad i = 1, 2, \dots, N_t; \quad \ell = 0, \dots, L \qquad (3.6) $$

which represents the discrete version of the solution field to (3.3) at time $t_i$. We can then define the snapshot matrix $S^{\ell}$ whose columns vectors are $\mathbf{s}_i^{\ell}$, viz.,

$$ S^{\ell} = \left[ \begin{array}{c|c} \mathbf{s}_1^{\ell} & \mathbf{s}_2^{\ell} \cdots \mathbf{s}_{N_t}^{\ell} \\ | & | \\ | & | \end{array} \right]; \quad \ell = 0, \dots, L \qquad (3.7) $$

Note that $S^{\ell}$ is an $N_x \times N_t$ real matrix. Then, the solution vector $\mathbf{v}_{\theta_0}^{\ell}(t)$ at any given time is represented as a finite sum, which is expressed as:

$$ \mathbf{v}_{\theta_0}^{\ell}(t) = \sum_{i=1}^{N_t} \tilde{a}_i^{\ell}(t) \mathbf{s}_i^{\ell}; \quad \ell = 0, \dots, L \qquad (3.8) $$

Here, the time-dependent coefficients $\{\tilde{a}_i^{\ell}\}_{i=1}^{N_t}$ are solutions to the system of ordinary differential equations, which emerge when substituting (3.8) into the problem (3.3). It is important to note that the field representation in equation (3.8) assumes that $\mathbf{v}_{\theta_0}^{\ell}$ lies in the span of the vectors $\{\mathbf{s}_1^{\ell}, \mathbf{s}_2^{\ell}, \dots, \mathbf{s}_{N_t}^{\ell}\}$. Hence, the latter need not be linearly independent, so it is crucial to address the potential redundancy within the set and reduce it to retain only the independent vectors. This reduction is necessary to minimize both storage and computational costs.

To emphasize the need for this reduction, consider that in seismic applications, the number of time iterations can be very large ($\sim 10^4$) and the value of $N_x$ can reach up to $10^9$. Storing all computed snapshots during the simulation becomes computationally prohibitive, if not outright impossible.
(b) Reduction Step.

There are several well-established methods for constructing an orthonormal basis from a given set of vectors, such as Singular Value Decomposition (SVD) [19, 30, 31, 15] and QR decomposition methods [19, 30, 31, 25]. Applying any of these techniques results in the construction of an orthonormal reduced set of modified snapshot vectors which form the reduced basis

$$[\phi_1^\ell, \phi_2^\ell, \dots, \phi_{M_\ell}^\ell]; \quad M_\ell \le N_t; \quad \ell = 0, \dots, L \qquad (3.9)$$

This reduction also leads to the full-rank modified snapshots matrix $\tilde{S}^\ell$ given by:

$$\tilde{S}^{\ell} = \left[ \begin{array}{cccc} \phi_{1}^{\ell} & & & \\ & \phi_{2}^{\ell} & & \\ & & \ddots & \\ & & & \phi_{M_{\ell}}^{\ell} \end{array} \right]; \quad \ell = 0, \dots, L \qquad (3.10)$$

with $M_\ell$ being always less than $N_t$. Hence, we represent $\mathbf{v}_{\theta_0}^\ell$ as follows:

$$\mathbf{v}_{\theta_0}^{\ell}(t) \approx \sum_{i=1}^{M_{\ell}} \hat{a}_i^{\ell}(t) \phi_i^{\ell}; \quad \ell = 0, \dots, L \qquad (3.11)$$

**Remark 3.1** Note that MOR-T$_0$ equipped with SVD is the classical POD.

Stage 2: Construction of the reduced model

This stage consists in selecting the final basis for MOR-T$_L$, following the reduction process and incorporating the Fréchet derivatives. It follows from **Stage 1**, the obtention of an augmented modified snapshot matrix given by

$$S^{\text{aug}} = [[\tilde{S}^0] \ [\tilde{S}^1] \ \dots \ [\tilde{S}^L]] \qquad (3.12)$$

where the block matrix $\tilde{S}^0, \dots, \tilde{S}^L$ are given by (3.10).

Note that the matrix $S^{\text{aug}}$ is therefore a $N_x \times \sum_{\ell=0}^{L} M_\ell$ matrix. In addition, each block matrix $\tilde{S}^\ell$ is a full-rank matrix since the column vectors form an orthonormal basis. However, this does not ensure that $S^{\text{aug}}$ is a full-rank matrix. Hence, modulo a reordering of the columns of (3.12), we apply again a reduction order procedure and obtain a new orthonormal set of vectors:

$$[\phi_1, \phi_2, \dots, \phi_N]; \quad N \le \sum_{\ell=0}^{L} M_\ell. \qquad (3.13)$$

The number $N$ is the result of applying the reduction method to the final augmented snapshot matrix. Therefore, for $\theta_\alpha = \theta_0 + \alpha \delta \theta_0$, $\mathbf{u}_{\theta_\alpha}$ can be approximated, using Taylor polynomial expansions of $u_\theta$ centered at $\theta_0$, as follows:

$$\mathbf{u}_{\theta_\alpha}(t) = \sum_{i=1}^{N} a_i(t) \phi_i \qquad (3.14)$$
The computation of ${a_i}_{i=1}^N$ requires solving the ODE system resulting from substituting (3.14) into the space-discretized problem (2.3) for the new value $\theta_\alpha$. Consequently, the proposed representation (3.14) is used to compute $u_\theta$ for arbitrary values of $\theta$ along the direction of $\delta\theta$, in a substantially faster way by solving a reduced system.

**Remark 3.2** The following points are worth nothing:

a. The right-hand sides of the boundary value problems in (3.3) do not depend on the source term $f$ because $f$ is independent of the parameter $\theta$. If $f$ depends on $\theta$, the method still applies, but with updated right-hand sides that include the successive Fréchet derivatives of $f$ with respect to $\theta$ up to order $L-1$.

b. MOR-$T_L$ can also handle perturbations of the parameter in multiple directions of the form:
$$ \theta_{\alpha} = \theta_0 + \sum_{k=1}^{K} \alpha_k \delta\theta_k \text{ with } \langle \delta\theta_i, \delta\theta_j \rangle = 0, \forall i \neq j. \text{ However, this introduces the need to compute the Fréchet derivatives for the cross-directions, thereby increasing the number of problems to solve combinatorially.} $$

## Final remark on the Computational costs and complexity of the MOR-$T_L$

In this section, we provide an overview of the computational costs associated with the MOR-$T_L$ method, which can be divided into two categories: (i) the Offline computational cost, which covers the time required to construct the reduced basis at a fixed parameter value $\theta_0$ as well as the reduced order system itself, and (ii) the Online computational cost, which accounts for the computation of the reduced-order system solution and the reconstruction of the numerical solution for an updated parameter $\theta_\alpha$.

### Offline Computational Cost

The offline computational cost primarily involves solving the discrete system with multiple right-hand sides and construct the reduced-order basis. These tasks are considered pre-processing steps since they are performed once for a fixed parameter $\theta_0$. Specifically, the components of the offline cost are as follows:

* **Construction of $L+1$ snapshots matrices $S^\ell$ (3.7):** This task requires solving a full boundary value problem with multiple right hand sides. The computational cost depends on the chosen numerical scheme for solving (3.3) and the solver (direct of iterative) employed for the resulting system.

* **Construction of $L+1$ modified (orthonormal) snapshots matrices $\tilde{S}^\ell$ (3.10):** This step involves a reduction technique. When using randomized SVD [19, 30, 31], the computational cost is of order $\sum_{\ell=0}^L \mathcal{O}(N_x N_t (M_\ell + p) + (N_x + N_t)(M_\ell + p)^2)$ with $p$ an oversampling, whereas using QR decomposition with the Gram-Schmidt orthonormalization procedure [19, 30, 31, 25] results in a cost of order $\sum_{\ell=0}^L \mathcal{O}(N_x M_\ell)$.
* **Construction of the final basis (3.13):** This requires applying a reduction technique (SVD or QR) to the augmented modified snapshot matrix (3.14). The computational cost of this task is similar to the previous step.

* **Construction of the reduced order system: (3.14):** This involves matrix products arising when substituting (3.14) into (2.3) for $\theta = \theta_{\alpha}$.

## Online Computational Cost

The online phase includes the computation of the reduced-order system solution and the reconstruction of the numerical solution. Specifically:

* **Computation of the coefficients {$a_i$}$_{i=1}^N$ in the representation (3.14):** This step requires solving the ROM, which is an ordinary differential system derived from substituting the field representation (3.14) into (2.3) for the updated parameter values $\theta_{\alpha} = \theta_0 + \alpha \delta \theta$. The computational cost for this task depends on the selected methodology for solving the ODE.

* **Reconstruction of the solution $u_{\theta_{\alpha}}$ (3.14):** This task involves computing scalar vector products to reconstruct the solution.

# 4 Application to Two-Dimensional Seismic Wave Problems

In this section, we demonstrate how to apply MOR-T_L to a specific problem, selecting a class of 2D seismic wave problems as our test case. The choice of seismic wave stems from our interest in monitoring CO₂ injection in subsurface environments, a process that involves solving an inverse problem where the goal is to reconstruct the velocity profile of the medium. Efficiently solving the corresponding direct problem is a crucial step in achieving this objective. In many inversion processes, such as Full Waveform Inversion (FWI) (see, e.g., [34, 28, 23]), it is necessary to solve the same direct problem multiple times, but with different velocity models. These models typically represent perturbations of the original medium, which are determined by solving an adjoint problem. As such, each iteration of the inversion requires solving multiple direct wave problems, often with fine discretization meshes to ensure high accuracy – especially important when solving ill-posed problems. This results in large-scale, computationally expensive problems, which can render the inversion process computationally infeasible.

We propose applying MOR-T_L to this context, demonstrating that a carefully designed ROM can significantly reduce computational costs while maintaining accuracy. This, in turn, can restore efficiency to inversion procedures like FWI.

## 4.1 The problem prototype

Here, we define the boundary value problem for the acoustic wave problem, followed by the set of problems that characterize the Fréchet derivatives of the solutions. The initial value problem we consider is given by:

$$ (\text{IVP}) \quad \begin{cases} \theta \frac{\partial^2 u_{\theta}}{\partial t^2}(\mathbf{x}, t) - \Delta u_{\theta}(\mathbf{x}, t) = f(\mathbf{x}, t) & \text{in } \mathbb{R}^2 \times (0, T) \\ u_{\theta}(\mathbf{x}, 0) = u_0(\mathbf{x}), \quad \frac{\partial u_{\theta}}{\partial t}(\mathbf{x}, 0) = u_1(\mathbf{x}) & \text{in } \mathbb{R}^2 \end{cases} \quad (4.1) $$
The parameter $\theta := \theta(\mathbf{x})$, $\mathbf{x} \in \mathbb{R}^2$, represents the propagation velocity as it writes $\theta(\mathbf{x}) = \frac{1}{c(\mathbf{x})^2}$ with $c(\mathbf{x}) > 0$ the functional representing the acoustic wave speed in the domain. The unknown $u_\theta$ is the pressure wave field, $f^{(0)} = f$ is the source term independent of the parameter $\theta$. Then, $v_\theta^\ell$, the Fréchet derivative of order $\ell$ in the direction of $\delta\theta$, satisfies the following collection of initial value problems

$$ \left\{ \begin{array}{ll} \theta \frac{\partial^2 v_\theta^\ell}{\partial t^2} - \Delta v_\theta^\ell = f^{(\ell)} & \text{in } \mathbb{R}^2 \times (0, T) \\ v_\theta^\ell(\mathbf{x}, 0) = 0, \quad \frac{\partial v_\theta^\ell}{\partial t}(\mathbf{x}, 0) = 0 & \text{in } \mathbb{R}^2 \end{array} \right. \qquad (4.2) $$

with $f^{(0)} = f$ and $f^{(\ell)} = -\ell\delta\theta \frac{\partial^2 v_\theta^{\ell-1}}{\partial t^2}$, $1 \le \ell \le L$.

**Remark 4.1** Due to the separable structure of the acoustic wave equation, with the wave operator splitting into space and time components, each right-hand side of the successive Fréchet derivatives depends on the solution from the previous order. As a result, all systems can be solved sequentially at each time step.

## 4.2 The formulation in a Bounded Domain

To solve the initial value problems (4.1) and (4.2) we use the SEM [29, 24]. SEM has the advantage of giving rise to a diagonal mass matrix which in turns can easily accommodate the use of explicit time discretization scheme. This is a clear advantage for large scale wave applications, hence helping in controlling an intensive use of memory.

As with any domain-based discretization scheme, the problem must first be reformulated within a bounded domain $\Omega$. This is achived by introducing an artificial external boundary, to which a boundary condition is applied. Therefore, the boundary value problems considered in the bounded domain are given by:

$$ \left\{ \begin{array}{ll} \theta \frac{\partial^2 v_\theta^\ell}{\partial t^2} - \Delta v_\theta^\ell = -\ell \frac{\partial^2 v_\theta^{\ell-1}}{\partial t^2} & \text{in } \Omega \times (0, T) \\ B_\theta v_\theta^\ell = 0 & \text{on } \partial\Omega \times (0, T) \\ v_\theta^\ell(\mathbf{x}, 0) = 0, \quad \frac{\partial v_\theta^\ell}{\partial t}(\mathbf{x}, 0) = 0 & \text{in } \Omega \\ v_\theta^0 = u_\theta & \end{array} \right. \qquad (4.3) $$

We have chosen two boundary conditions for this formulation. The first is the Dirichlet boundary condition, where $B_\theta$ is the identity operator, chosen for its simplicity in implementation. The second condition is a first order Absorbing Boundary Condition (ABC) [18, 16] where the operator $B_\theta$ depends on the problem's parameter $\theta$. This choice enables us to assess the efficacy and robustness of the reduced basis with respect to changes in the boundary condition resulting from parameter variations. We aim to verify that MOR-T$_L$ continues to provide a solution with an accuracy level comparable to that obtained by solving the full problem using SEM when both methods are equipped with the same exterior boundary condition. It is important to emphasize that when solving practical or more realistic problems, more accurate boundary conditions would be necessary.
## 4.3 The basis functions

In this section, we describe the practical construction of basis functions used in MOR-T_L. It is common in reduced-order methods for the basis functions to be computed using a truncated SVD of the snapshot matrices $S^l$, given by equation (3.7). This approach is based on POD [1, 36, 4, 20]. However, there are several limitations to this method, with the following three being particularly noteworthy: (a) the lack of a rigorous criterion - or at least practical guidelines - for selecting the snapshots in $S^l$, (b) the potential for $S^l$ to become too large to store in memory for large-scale systems, and (c) even when can be stored, the decomposition process itself can become computationally excessive and ultimately prohibitive.

Hence, we propose a strategy where the SVD is replaced with a truncated QR decomposition based on a Gram-Schmidt (GS) orthonormalization process [25] allowing us to compute the basis as we solve (4.1) without storing the full matrices $S^l$. This is the purpose of the following.

Let $s_1^l$ be the snapshot of the (nonzero) solution at a time $t_1^l$ given by (3.6). Let $\phi_1^l$ be the first basis vector defined as:

$$ \phi_1^l = \frac{s_1^l}{\|s_1^l\|}; \quad l = 0, \dots, L $$

where $\|\cdot\|$ should be understood as the discrete norm associated with the discrete chosen inner product $\langle \cdot, \cdot \rangle$.

We compute the $i_\ell^{th}$ basis functions $i_\ell > 1$ in terms of the solution $s_j^l, j = 1, \dots, N_t$, and the previous basis vectors $\phi_k^l, k = 1, \dots, i_\ell - 1$ as:

$$ \tilde{\phi}_{i_\ell}^l = s_j^l - \sum_{k=1}^{i_\ell-1} \langle s_j^l, \phi_k^l \rangle \phi_k^l; \quad \ell = 0, \dots, L. \quad (4.4) $$

If the following condition is satisfied

$$ \|\tilde{\phi}_{i_\ell}^l\|^2 > \epsilon \|s_j^l\|^2; \quad \ell = 0, \dots, L \quad (4.5) $$

for a prescribed value for the positive number $\epsilon$, we proceed to the normalization of the vector.

$$ \phi_{i_\ell}^l = \frac{\tilde{\phi}_{i_\ell}^l}{\|\tilde{\phi}_{i_\ell}^l\|}; \quad \ell = 0, \dots, L. \quad (4.6) $$

If not, the basis function is not selected. We re-iterate the process ((4.4) and (4.5)) and compute $\tilde{\phi}_{i_\ell}^l$ with $s_{j+1}^l$.

The stopping criterion (4.5) gives us a quantitative way to select the snapshots. However, GS process is subject to instability error for a large amount of data and may bias this criterion. To reduce this effect, re-orthogonalization must be applied every few selected snapshots. At the end of the process/simulation, we have constructed a basis $\{\phi_i^l\}_{i=1}^{M_\ell}$, with $M_\ell \le N_t$. In practice, $M_\ell$ depends on the source considered $f$ as well as on the values of the parameter $\theta_0$.

To create the final basis, we gather all the modified snapshot matrices $\tilde{S}^l$ in the augmented snapshot matrix $S^{aug}$ given by (3.12). Then, before applying the reduction order technique on $S^{aug}$, we proceed to a re-ordering of the column vectors by increasing time, that is, we start from the vector $\phi_1^0$ obtained at time $t_1^0$, the next vectors are all the vectors $\phi_1^l$ obtained at time $t_1^l$, if
any, Then we consider the vectors obtained at $t_2^\ell$, if any, starting from $\ell = 0$ up to $\ell = L$, etc. We then obtain after the completion of Gram-Schmidt process a modified full-rank augmented snapshot matrix $\tilde{S}^{\text{aug}}$ given by (3.14) with our final reduced basis being $\{\phi_i\}_{i=1}^N$, $N \le \sum_{\ell=0}^L M_\ell$.

**Remark 4.2** The following points are noteworthy:

a. In addition to being memory-efficient compared to the SVD approach, the QR strategy is easy to parallelize, as each scalar product can be computed independently over subdomains. This parallelization accelerates the process and helps reduce numerical instabilities in the scalar products.

b. The method's efficiency depends on the choice of $\epsilon$. If $\epsilon$ is too small, it may lead to the selection of an excessively large number of basis functions, increasing computation time and potentially causing numerical instabilities. On the other hand, if $\epsilon$ is too large, the number of basis functions may be insufficient, resulting in an overly limited approximation subspace. Therefore, selecting $\epsilon$ requires a balance. In practice, we find that setting $\epsilon = 0.01$ offers a good compromise, yielding excellent results.

## 4.4 The resulting reduced order model

We now introduce the reduced-order model, represented by the boundary value problem

$$ 
\left\{
\begin{array}{ll}
\displaystyle \sum_{i=1}^{N} \theta \frac{\partial^2 a_i(t)}{\partial t^2} \phi_i(\mathbf{x}) - a_i(t) \Delta \phi_i(\mathbf{x}) = f(\mathbf{x}, t) & \text{in } \Omega \times ]0, T[ \\
\\
B_{\theta} \sum_{i=1}^{N} a_i(t) \phi_i(\mathbf{x}) = 0 & \text{on } \partial \Omega \times ]0, T[ \\
\\
a_i(0) = \langle u_0(\mathbf{x}), \phi_i(\mathbf{x}) \rangle, \quad \frac{\partial a_i}{\partial t}(0) = \langle u_1(\mathbf{x}), \phi_i(\mathbf{x}) \rangle & \text{in } \Omega
\end{array}
\right.
\quad (4.7) $$

resulting from substituting (3.14) into (4.1).

Since $u_{\theta_0}(\mathbf{x}, t)$ is supposed to be obtained through a SEM formulation, then the $\{\phi_i(\mathbf{x})\}_{i=1}^N$ can be written in the SEM basis as $\phi_i(\mathbf{x}) = \sum_{k=1}^{N_x} \phi_{k,i} \psi_k(\mathbf{x})$. Considering the case where $B_{\theta} = \sqrt{\theta} \frac{\partial}{\partial t} + \frac{\partial}{\partial \mathbf{n}}$,

and gathering all the $\phi_k, i$ in a matrix $\phi$, the $f_k$ in a vector $\mathbf{f}$, and the $a_i(t)$ in a vector $\mathbf{a}(t)$, $\forall k = 1, \dots, N_x, \forall i = 1, \dots, N$, one can write after introducing a weak formulation of (4.7):

$$ \phi^T \mathbf{M}_\theta \phi \frac{\partial^2 \mathbf{a}(t)}{\partial t^2} + \phi^T \mathbf{D}_{\sqrt{\theta}} \phi \frac{\partial \mathbf{a}(t)}{\partial t} + \phi^T \mathbf{K} \phi \mathbf{a}(t) = \phi^T \mathbf{f}(t) \quad (4.8) $$

where $\mathbf{M}_\theta, \mathbf{D}_{\sqrt{\theta}}$, and $\mathbf{K}$ are the mass, the damping, and the stiffness matrices whose respective entries are:

$$ \int_{\Omega} \theta \psi_k(\mathbf{x}) \psi_l(\mathbf{x}) d\mathbf{x}; \quad \int_{\partial \Omega} \sqrt{\theta} \psi_k(\mathbf{x}) \psi_l(\mathbf{x}) d\mathbf{S}; \quad \int_{\Omega} \nabla \psi_k(\mathbf{x}) \cdot \nabla \psi_l(\mathbf{x}) d\mathbf{x} \quad (4.9) $$

For the time discretization we choose a leap-frog scheme for simplicity, but any other suitable scheme may be used. We set $\mathbf{M}_\phi = \phi^T \mathbf{M}_\theta \phi$, $\mathbf{D}_\phi = \phi^T \mathbf{D}_{\sqrt{\theta}} \phi$, $\mathbf{K}_\phi = \phi^T \mathbf{K} \phi$ and $\mathbf{f}_\phi = \phi^T \mathbf{f}$. Therefore, for $n = 1, \dots, N_t - 1$,
$$ \mathbf{a}^{n+1} = \left( \mathbf{M}_{\phi} + \frac{\Delta t}{2} \mathbf{D}_{\phi} \right)^{-1} \left( \left( 2\mathbf{M}_{\phi} - \Delta t^2 \mathbf{K}_{\phi} \right) \mathbf{a}^n - \left( \mathbf{M}_{\phi} - \frac{\Delta t}{2} \mathbf{D}_{\phi} \right) \mathbf{a}^{n-1} + \Delta t^2 \mathbf{f}_{\phi} \right). \quad (4.10) $$

**Remark 4.3** In the case of a homogeneous Dirichlet boundary condition, that is, $\sum_{i=1}^{N} a_i(t)\phi_i(\mathbf{x}) = 0, \forall \mathbf{x} \in \partial\Omega, \mathbf{D}_{\sqrt{\bar{\theta}}} = 0$. Then, modulo a change in the stiffness matrix and the source term, we obtain:

$$ \mathbf{a}^{n+1} = 2\mathbf{a}^n - \mathbf{a}^{n-1} - \Delta t^2 \mathbf{M}_{\phi}^{-1} (\mathbf{K}_{\phi}\mathbf{a}^n - \mathbf{f}_{\phi}); \quad n = 1, \dots, N_t - 1 \quad (4.11) $$

**Remark 4.4** A judicious choice of inner product during the QR decomposition, as described in Section 4.3, can enhance the computational efficiency of the MOR-T_L strategy. For the application under consideration, we select the inner product associated with the stiffness matrix $\mathbf{K}$, defined as $\langle u, v \rangle_{\mathbf{K}} = \mathbf{u}^T \mathbf{K} \mathbf{v}$. This choice offers two key computational advantages: (a) due to the leap-frog time-stepping scheme, only half of the inner product is inherently computed at each time step when solving the resulting SEM linear system, and (b) it eliminates the need to compute the projections $\phi^T \mathbf{K} \phi$, as the basis is orthonormal with respect to $\mathbf{K}$ by construction, thereby resulting in significant computational cost savings. The only projections that need to be computed, as indicated in (4.7), are those associated with the mass $\mathbf{M}_{\theta}$ and damping $\mathbf{D}_{\sqrt{\bar{\theta}}}$ matrices. These projections are much less computationally demanding, since both the mass and damping matrices are diagonal, a property arising from the SEM formulation.

# 5 Assessing cost-effectiveness and performances²of MOR-T_L v.s. SEM in wave problems

In this section, we present a numerical investigation aimed at evaluating the cost-effectiveness of the proposed MOR strategy compared to the standard SEM approach for solving the class of two-dimensional wave problems (4.1). Our analysis is divided into two parts. First, we examine the cost-effectiveness of using QR versus SVD basis representations in the context of MOR. Next, we compare the computational cost of MOR-T_L with that of the SEM method, assessing performance for a range of velocity values and a prescribed accuracy level. We conduct these experiments for both homogeneous and heterogeneous media to capture a broad set of scenarios. Additionally, we investigate the impact of the exterior artificial boundary choice on computational cost, presenting results with both the simplest homogeneous Dirichlet boundary condition and a more realistic absorbing boundary condition.

We consider $u(\mathbf{x}, 0) = 0, \frac{\partial u}{\partial t}(\mathbf{x}, 0) = 0$, but a Ricker source of the form:

$$ f(\mathbf{x}, t) = \left(1 - \frac{1}{2}\omega^2(t - t_0)^2\right) e^{-\frac{1}{4}\omega^2(t - t_0)^2} \delta(\mathbf{x} - \mathbf{x}_s) \quad (5.1) $$

where $\omega = 2\pi f_s$ the pulsation with $f_s$ the central source frequency and $t_0$ a delay.

<footnote n="2">²All these numerical experiments have been performed on a system equipped with an 11th Gen Intel Core i5-1135G7 (4 cores, 8 threads, 2.40GHz), Mesa Intel Xe Graphics (TGL GT2). Python 3.10.12 with Numpy 1.21.5 multi-threaded for BLAS operations</footnote>
## 5.1 Assessing the Cost-Effectiveness of QR and SVD Basis Representations

In this subsection, we present numerical results comparing the performance and cost-effectiveness of the QR and SVD basis representations. We consider the initial boundary value problem (4.1), defined over a 50m × 50m square computational domain Ω, with a source function f represented by the Ricker wavelet (5.1) with a central frequency of $f_s = 2.5$Hz and homogeneous initial conditions, i.e., $u_0 = u_1 = 0$. The source is located at $\mathbf{x}_s = (25, 5)$ as depicted in Figure 15.

The reference solutions are computed using the SEM method with a 201 × 201 Q1 finite element mesh [24]. The experiments are run over 4000 time steps with $\Delta t = 0.001$, and the model parameter is $\theta = \frac{1}{c^2}$, where c is the velocity in $\Omega$. These experiments are conducted on a fixed, unperturbed velocity, with both homogeneous and heterogeneous media considered. For clarity, results are presented for $\epsilon = 0.01$, though results for other values of $\epsilon$ are provided in [3].

### 5.1.1 Case of Homogeneous Medium

We begin with a homogeneous medium where the velocity is constant, $c(x) = 15$ m/s in $\Omega$. Here, we assess the cost-effectiveness of QR and SVD basis representations in two boundary condition scenarios: a simple homogeneous boundary condition and a more realistic absorbing boundary condition.

**Case of Homogeneous Dirichlet boundary condition.**

The results for the homogeneous Dirichlet boundary condition are shown in Figures 1 - 3 and Table 1. Notable observations include:

* Figure 1 illustrates sample basis functions for both the SVD and QR decompositions. With $\epsilon = 0.01$, the QR method yields $N = 81$ basis elements. The figure highlights that the QR basis functions capture the finite speed propagation of the wave support, while the SVD functions preserve the frequency content. The QR approach maintains time causality, as each basis function is computed from linear combinations of the solution $u_{\theta_0}(\mathbf{x}, t)$ at the current and previous time steps (see (4.4)).

* Figure 2 and 3 show comparisons of the accuracy of the ROM solutions using QR and SVD bases. Time traces (Figure 2) at a selected location $\mathbf{x}_r = (39.05, 15.42)$ demonstrate excellent agreement with the reference SEM solution for both QR and SVD methods. Similar results were obtained at other source and receiver locations. Figure 3 plots the relative error in $L^2$ norm over the computational domain $\Omega$ at every $100^{th}$ time step, confirming that both representations achieve comparable accuracy (below $10^{-2}$%). Initially, QR significantly outperforms SVD, likely due to the QR method's first basis function being the exact normalized solution $u_{\theta_0}(\mathbf{x}, t)$ at $t = t_1$. Additionally, the results in [3] confirm that varying the value of $\epsilon$ affects the solution's accuracy, as expected. However, it is also demonstrated in [3] that decreasing $\epsilon$ below 0.01 has little to no impact on accuracy. In terms of online computational, the ROM simulation has a speed up of approximately $7.3 \times 10^2$ relatively to the SEM one.

* Table 1 demonstrates the superior computational efficiency of the MOR-T$_0$ formulation. The MOR-T$_0$ method requires approximately 2 less computational time to achieve comparable,
or even better, accuracy than the POD approach with randomized SVD (rSVD). It is also shown in [3] that, as the number of degrees of freedom increases, the QR method outperforms the rSVD in terms of both computational speed and memory efficiency.

* Final remark: when the source function ceases emission, increasing $\epsilon$ causes the QR algorithm to separate consecutive basis functions by an amount corresponding to approximately half a wavelength. This behavior allows significant computational cost savings, as fewer snapshots need to be evaluated, reducing the overall computational load.

<figure><img src="image_3.png" /><figcaption>Figure 1: Sample of basis functions: SVD (left) v.s. QR (right).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2: Pressure Field Comparison at $x_r = (39.05, 15.42)$ for a Homogeneous Medium.

Pressure values over $[0, 410^3]$ are compared between SEM (reference), MOR-T$_0$, and POD for a homogeneous medium with boundary condition $B_\theta = I$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Relative Errors Between MOR-T$_0$, POD, and SEM.

Relative errors over the time interval $[0, 4]$ are shown for MOR-T$_0$, and POD compared to SEM. The error is computed as the standard relative $L^2$-norm error across the computational domain $\Omega$ at each time step, for a homogeneous medium $\Omega$ with boundary condition $B_\theta = I$.</figcaption></figure>

<table><caption>Table 1: Comparison of MOR-T$_0$ and POD: Total Computational Time and Relative Errors.</caption><thead><tr><th>Method</th><th>Offline Computational Time</th><th>1/N<sub>t</sub> &sum;<sub>i=1</sub><sup>N<sub>t</sub></sup> ||u<sub>SEM</sub>(x, t<sub>i</sub>) - u<sub>ROM</sub>(x, t<sub>i</sub>)||<sub>2</sub> / ||u<sub>SEM</sub>(x, t<sub>i</sub>)||<sub>2</sub></th></tr></thead><tbody><tr><td>rSVD</td><td>110.0s</td><td>6.54 &times; 10<sup>-4</sup></td></tr><tr><td>QR</td><td>57.1s</td><td>5.19 &times; 10<sup>-5</sup></td></tr></tbody></table>

This table presents the total computational time (offline time for constructing the basis and online time for solving the system) and the average relative $L^2$-norm error for MOR-T$_0$ and POD. The errors are computed across the computational domain $\Omega$ and averaged over time steps. These results are obtained in the case of a homogeneous medium $\Omega$ and a homogeneous boundary condition $B_\theta = I$.

In conclusion, this experiment clearly establishes the MOR-T$_0$ formulation's superiority over POD in terms of computational efficiency. The QR-based approach also offers strong parallelizability, which can further reduce computation time on parallel platforms. Furthermore, QR decomposition requires less storage since only one snapshot of the solution needs to be stored at any given time, while the SVD approach requires storing the entire matrix $S$. Thus, MOR-T$_0$ outperforms POD both in terms of computational cost and storage requirements – a crucial consideration for high-frequency wave simulations.

### Case of Absorbing Boundary Condition.

We now turn to the case where the boundary condition is replaced by an absorbing boundary condition, characterized by the operator $B_\theta = \sqrt{\theta} \frac{\partial}{\partial t} + \frac{\partial}{\partial \mathbf{n}}$. The results, shown in Figures 4 and 5, as well as Table 2, reveal similar behavior to the Dirichlet case:

* Figure 4 shows a perfect match between the time trace solutions from both MOR-T$_0$ and POD methods and the reference SEM solution.
* Figure 5, demonstrates that both methods achieve a comparable accuracy (within $10^{-3}\%$), although the POD method exhibits a noticeable error increase at the initial and final time steps.

* Table 2 again confirms the superior computational efficiency of MOR-T<sub>0</sub>. Specifically, MOR-T<sub>0</sub> reduces computational time by approximately 4 times, while providing an accuracy three orders of magnitude better than POD.

In summary, the use of a more sophisticated absorbing boundary condition does not diminish the performance of MOR-T<sub>0</sub>, maintaining its computational superiority over the POD method.

<figure><img src="image_6.png" /><figcaption>Figure 4: Pressure Field Comparison at $x_r = (39.05, 15.42)$ for a Heterogeneous Medium.</figcaption></figure>

Pressure values over the time interval [0, 4] are compared between SEM (reference), MOR-T<sub>0</sub>, and POD for a homogeneous medium with boundary condition $B_\theta = \sqrt{\theta} \frac{\partial}{\partial t} + \frac{\partial}{\partial \mathbf{n}}$.

<figure><img src="image_7.png" /><figcaption>Figure 5: Relative Errors Between MOR-T<sub>0</sub>, POD, and SEM.</figcaption></figure>

Relative errors over the time interval [0, 4] are shown for MOR-T<sub>0</sub> and POD compared to SEM. The error is computed as the standard relative $L^2$-norm error across the computational domain $\Omega$ at each time step, for a homogeneous medium $\Omega$ with boundary condition $B_\theta = \sqrt{\theta} \frac{\partial}{\partial t} + \frac{\partial}{\partial \mathbf{n}}$.

<table><caption>Table 2: Comparison of MOR-T<sub>0</sub> and POD: Total Computational Time and Relative Errors.</caption><thead><tr><th>Method</th><th>Total Computational Time</th><th>1 / N<sub>t</sub> &sum;<sub>i=1</sub><sup>N<sub>t</sub></sup> ||u<sub>SEM</sub>(<b>x</b>, t<sub>i</sub>) - u<sub>ROM</sub>(<b>x</b>, t<sub>i</sub>)||<sub>2</sub> / ||u<sub>SEM</sub>(<b>x</b>, t<sub>i</sub>)||<sub>2</sub></th></tr></thead><tbody><tr><td>SVD</td><td>132.2s</td><td>4.04 &times; 10<sup>-2</sup></td></tr><tr><td>QR</td><td>32.4s</td><td>5.56 &times; 10<sup>-5</sup></td></tr></tbody></table>

This table presents total computational time (offline time for constructing the basis and online time for solving the system) and the average relative $L^2$-norm error for MOR-T<sub>0</sub> and POD. The errors are computed across the computational domain $\Omega$ and averaged over time steps. These results are obtained in the case of a homogeneous medium $\Omega$ and a boundary condition $B_\theta = \sqrt{\theta} \frac{\partial}{\partial t} + \frac{\partial}{\partial \mathbf{n}}$.

### 5.1.2 Case of Heterogeneous Medium

For the heterogeneous medium, we consider a smooth velocity field $c(\mathbf{x})$, as shown in Figure 6. Results are presented for the absorbing boundary condition, with further details available in [3]. The time trace solutions and accuracy levels follow the same pattern as the homogeneous case, with
perfect agreement and impressive accuracy (below 10⁻²%), as depicted in Figures 7 and 8. This further supports the robustness of the MOR-T₀ method across different medium configurations.

<figure><img src="image_8.png" /><figcaption>Figure 6: Velocity Distribution in a Heterogeneous Medium.</figcaption></figure>

The figure shows the velocity $c_0 = \frac{1}{\sqrt{\theta_0}}$ in the heterogeneous medium $\Omega$, with depth represented on the y-axis and the x-coordinate on the x-axis. The color map indicates the velocity values.

<figure><img src="image_9.png" /><figcaption>Figure 7: Pressure Field Comparison at $x_r = (39.05, 15.42)$ for a Heterogeneous Medium.</figcaption></figure>

Pressure values over the time interval [0, 4] are compared between SEM (reference), MOR-T₀, and POD for a heterogeneous medium with boundary condition $B_\theta = \sqrt{\theta} \frac{\partial}{\partial t} + \frac{\partial}{\partial n}$.

<figure><img src="image_10.png" /><figcaption>Figure 8: Relative Errors Between MOR-T₀, POD, and SEM.</figcaption></figure>

Relative errors over the time interval [0, 4] are shown for MOR-T₀ and POD compared to SEM. The error is computed as the standard relative L²-norm error across the computational domain $\Omega$ at each time step, for a heterogeneous medium $\Omega$ with boundary condition $B_\theta = \sqrt{\theta} \frac{\partial}{\partial t} + \frac{\partial}{\partial n}$.

## 5.2 Comparative Analysis of MOR-T_L and SEM in Line Search Minimization

In this subsection, we evaluate both the accuracy and computational efficiency of the proposed MOR-T_L strategy within the context of line search minimization, a commonly employed iterative
process in FWI. To clarify, a line search is an iterative optimization method that minimizes a
functional *J* along a given perturbation direction, starting from an initial model θ₀. In the following
analysis, we consider the perturbation velocity *δθ*₀, which is obtained through an adjoint state
method. The line search process is then applied to determine the optimal step size α*, such that
θ₀α* *δθ*₀ represents an admissible updated velocity.

<figure><img src="image_11.png" /><figcaption>Figure 9: Normalized Perturbation Distribution.
The figure shows the normalized velocity perturbation $\delta\theta_0$
in the heterogeneous medium $\Omega$, with depth represented
on the y-axis and the x-coordinate on the x-axis. The
color map indicates the normalized velocity perturbation
values.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Updated Velocity Distribution in a Heterogeneous Medium.
The figure shows the updated velocity perturbation $c_1 = \frac{1}{\sqrt{\theta_0 + 2.12 \times 10^{-3} \delta\theta_0}}$ in the heterogeneous medium $\Omega$, with depth represented on the y-axis and the x-coordinate on the x-axis. The color map indicates the velocity values.</figcaption></figure>

**5.2.1 Accuracy of the MOR-T<sub>L</sub> in Computing Velocity Perturbations**

To perform a line search minimization, we consider the perturbation $\delta\theta_0$ as shown in Figure 9. The updated velocity takes the form $\theta_1 = \theta_0 + \alpha\delta\theta_0$, with $\alpha = 2.12 \times 10^{-3}$ resulting in a perturbation size $\alpha \frac{\|\delta\theta_0\|_2}{\|\theta_0\|_2}$ of approximately 10%. This perturbation represents a substantial change in the initial velocity (see Figure 10 for the updated $\theta_1$, compared with the original velocity $\theta_0$ in Figure 6). We then apply the proposed MOR-T$_L$ strategy to compute the pressure field $u_{\theta_1}$ and compare the results with those from the standard SEM approach, which serves as the reference solution. In particular, we investigate the sensitivity of MOR-T$_0$'s accuracy to the degree $L$ of the Taylor polynomial (see Equation (2.2)). The results are summarized in Figures 11-14 and Table 3. Key observations include:

* Figure 11 compares the pressure field at the point $\mathbf{x}_r = (39.05, 15.42)$ over the time interval $[0, 4]$ for Taylor polynomial of different degrees. When $L=0$, corresponding to the standard MOR method, a phase-shift error is observed relative to the reference solution. This discrepancy is concerning, as it could significantly impact the convergence of the line search procedure. As expected, increasing the degree $L$ improves accuracy. Remarkably, the MOR-T$_L$ method achieves excellent accuracy with as few as four terms in the Taylor expansion
$(L = 3).$

* Figure 12 shows the relative error in the L²-norm over the entire spatial domain Ω and the time step interval. Unlike Figure 11, which focuses on local errors, Figure 12 presents global error estimates. The results confirm that, for $L = 0$, the standard ROM method produces unacceptable errors, while the MOR-T_L method achieves high accuracy even for small L values, particularly for $L = 3$.

* Figure 13 and Table 3 provide an initial assessment of the computational cost of the proposed MOR-T_L method. The number of basis functions increases linearly with L, and for $L = 3$, the offline computational time is about 4.6 times the SEM simulation, with a 6% accuracy level.

* Figure 14 indicates that the MOR-T_L method maintains acceptable accuracy even when the perturbation $\alpha\delta\theta_0$ reaches 15%, a relatively high value for the perturbation. This suggests that the proposed method is robust to large perturbations.

<figure><img src="image_13.png" /><figcaption>Figure 11: Accuracy of MOR-T_L and Sensitivity to Taylor Polynomial Degree.</figcaption></figure>

This figure consists of six subfigures illustrating the accuracy of the proposed MOR-T_L method and its sensitivity to the degree L of the Taylor polynomial expansion (from $L = 0$ to $L = 5$). Each subfigure shows a comparison of the pressure field at $x_r = (39.05, 15.42)$ over the time $[0, 4]$ for the updated velocity $c_1$ with results from SEM (reference) and MOR-T_L.
<figure><img src="image_14.png" /><figcaption>Figure 12: Relative L2-norm Errors for MOR-T<sub>L</sub> with Varying Taylor Polynomial Degree.</figcaption></figure>

This figure presents six curves, each representing the relative L²-norm error for MOR-T<sub>L</sub> with different degrees L of the Taylor polynomial expansion (from L = 0 to L = 5). The errors are computed across the computational domain Ω at each time step in the interval [0, 4 × 10³] for the updated velocity c₁.

<figure><img src="image_15.png" /><figcaption>Figure 13: Dependence of MOR-T<sub>L</sub> Basis Functions on Taylor Polynomial Degree L.</figcaption></figure>

This figure illustrates the linear relationship between the number of basis functions required by the MOR-T<sub>L</sub> method and the degree L of the Taylor polynomial expansion. As L increases, the number of basis functions grows linearly.

<table><caption>Table 3: MOR-T<sub>L</sub>: Offline Computational Time and Average Relative Errors for Varying Taylor Polynomial Degree L relatively to the SEM.</caption><thead><tr><th>L</th><th>Offline Computational Time<br>SEM Computational Time</th><th>1 / N<sub>t</sub> &sum;<sub>i=1</sub><sup>N<sub>t</sub></sup> ||u<sub>SEM</sub>(x, t<sub>i</sub>) - u<sub>ROM</sub>(x, t<sub>i</sub>)||<sub>2</sub> / ||u<sub>SEM</sub>(x, t<sub>i</sub>)||<sub>2</sub></th></tr></thead><tbody><tr><td>0</td><td>1.08</td><td>0.535</td></tr><tr><td>1</td><td>2.21</td><td>0.198</td></tr><tr><td>2</td><td>3.27</td><td>0.113</td></tr><tr><td>3</td><td>4.59</td><td>0.061</td></tr><tr><td>4</td><td>5.85</td><td>0.055</td></tr><tr><td>5</td><td>7.24</td><td>0.053</td></tr></tbody></table>

This table compares the offline computational time and the average relative L²-norm error for MOR-T<sub>L</sub> across different values of L, the Taylor polynomial degree. The errors are computed over the computational domain Ω and averaged over time steps for the updated velocity c₁.
<figure><img src="image_16.png" /><figcaption>Figure 14: Sensitivity of MOR-T<sub>L</sub> Accuracy to Perturbation Size for Varying Taylor Polynomial Degree L.</figcaption></figure>

This figure shows the sensitivity of MOR-T<sub>L</sub> accuracy to the size of the perturbation αδθ. The x-axis represents the perturbation size in percentage, while the y-axis shows the average relative errors. Six curves are plotted, each corresponding to a different Taylor polynomial degree L (from L=0 to L=5). The results are computed for the updated velocity c<sub>1</sub>.

In summary, these results emphasize the inadequacy of the standard ROM approach for accurately computing pressure fields at updated velocities. The proposed MOR-T<sub>L</sub> method, however, delivers high accuracy even with a cubic Taylor expansion (L = 3), and its computational cost remains manageable due to the linear growth of the basis functions with L.

5.2.2 Computational Efficiency and Robustness of the MOR-T<sub>L</sub> During Line Search

In this section, we evaluate the computational efficiency and robustness of the MOR-T_L method in the line search, particularly in the presence of noisy measurements. Specifically, the measurements of the pressure field vector $u_{obs}$ is obtained by computing the pressure $u_{\theta^*}$ corresponding to the updated model parameter $\theta_0 + \alpha^* \delta\theta_0$ with $\alpha^* = 2.12 \times 10^{-3}$ over the time interval [0, 4] on a line of 202 equidistant receivers from $\mathbf{x}_{r1} = (0, 15.42)$ to $\mathbf{x}_{r202} = (50, 15.42)$.
<figure><img src="image_17.png" /><figcaption>Figure 15: 2D computational domain Ω showing its depth and the horizontal axis.
The figure illustrates the location of the source f and the positions of the 202 receivers, which are aligned along a
horizontal line to measure the pressure field.</figcaption></figure>

These synthetic measurements are obtained using the SEM method to avoid the so-called inverse crime. We then compute the pressure fields $u_{\theta_{\alpha}}$ for a set of velocities $\theta_{\alpha} = \theta_0 + \alpha \delta \theta_0$ using the MOR-T$_L$ method. The goal of the line search is to minimize the cost function $J(\alpha)$ defined as:

$$J(\alpha) = \frac{1}{2} \|u_{\theta_{\alpha}} - u_{\text{obs}}\|_{2}^{2} \qquad (5.2)$$

by evaluating it at several values of $\alpha$. We also introduce white noise into the vector measurements $u_{\text{obs}}$ (see Figure 16) to assess the method's robustness under noisy conditions. The results are presented in Figures 16 and Tables 3-5, with the following key findings:

* Preliminary Remark: The line search is conducted by evaluating $J(\alpha)$ at 200 points within the interval $[0, 4 \times 10^{-3}]$ which is uniformly subdivided. This interval is chosen to account for perturbations $\alpha \delta \theta_0$ of up to 20% on the relative errors, representing relatively large perturbations. Negative values of $\alpha$ are excluded to ensure that the perturbation $\alpha \delta \theta_0$ aligns with the direction of $-\nabla_{\theta} J$, as required by the steepest descent method.

* Figure 16(a) shows results for noise-free measurements. The standard MOR method ($L=0$) fails to identify the optimal $\alpha^*$, as expected given the sensitivity of the ROM basis to model parameters. In contrast, the MOR-T$_L$ method accurately determines $\alpha^*$ even when using a linear Taylor expansion ($L=1$). For all non-constant Taylor polynomials, $J(\alpha)$ exhibits a distinct minimum around the correct value of $\alpha^*$, facilitating an unambiguous solution.

* Figure 16(b)-(c) show results for noise levels of 1% and 5%, respectively. Despite the noise, the MOR-T$_L$ method remains robust and continues to determine $\alpha^*$ with impressive accuracy, achieving a relative error of less than 1% in both cases, as reported in Table 4.

* Table 5 compares the computational costs of the MOR-T$_L$ and SEM methods when evaluating $J(\alpha)$ at 10 parameter values, consistent with the use of Wolfe conditions in FWI to limit the number of evaluations. We observe that MOR-T$_3$ reduces the computational time by about
50%. While evaluating $J(\alpha)$ at a small number of points may compromise the accuracy of $\alpha^*$ and slow down convergence, the MOR-T$_L$ method still outperforms SEM in terms of computational cost and accuracy.

* Final Remark: The use of SEM in the line search stage can result in prohibitive computational costs when a large number of parameter points are needed. For example, evaluating 200 points with SEM requires over three hours of computation, while MOR-T$_L$ with a Taylor polynomial degree $L$ smaller than 5 completes the task in under 10 minutes. This efficiency allows for the use of many points during the line search, which is crucial for ensuring fast convergence and accuracy in the FWI process. In contrast, relying on only a few points based on the Wolfe conditions – especially if the initial guess is outside the pre-asymptotic region – can significantly hinder convergence, or even prevent it entirely. Consequently, any computational savings from using fewer points may be offset by the need for more iterations, ultimately negating the intended efficiency gains.

<figure><img src="image_18.png" /><figcaption>(a) Measurement (left) and misfit function $J(\alpha)$ (right) for noise level: 0%</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Measurement (left) and misfit function $J(\alpha)$ (right) for noise level: 1%</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 16: Pressure Field Measurements and Cost Function with Varying Noise Levels.</figcaption></figure>

This figure consists of six subfigures. The left three show pressure field measurements at $x_r = (39.05, 15.42)$ over the time interval $[0, 4]$ for the velocity $c^*$ obtained with SEM and corrupted with 0%, 1%, and 5% white noise (top to bottom). The right three subfigures display the cost function $J(\alpha)$ evaluated at 200 parameter values, with six curves corresponding to Taylor polynomial degrees $L = 0$ to $L = 5$ for each noise level.

<table><caption>Table 4: Retrieved Parameter $\alpha^*$ by MOR-T<sub>L</sub> for Varying Noise Levels and Taylor Polynomial Degree L.</caption><thead><tr><th rowspan="2">L</th><th colspan="3">argmin J(α) for noise level</th></tr><tr><th>0%</th><th>1%</th><th>5%</th></tr></thead><tbody><tr><td>0</td><td>3.7 × 10<sup>-3</sup></td><td>3.6 × 10<sup>-3</sup></td><td>3.6 × 10<sup>-3</sup></td></tr><tr><td>1</td><td>2.12 × 10<sup>-3</sup></td><td>2.04 × 10<sup>-3</sup></td><td>2.04 × 10<sup>-3</sup></td></tr><tr><td>2</td><td>2.08 × 10<sup>-3</sup></td><td>2.2 × 10<sup>-3</sup></td><td>2.2 × 10<sup>-3</sup></td></tr><tr><td>3</td><td>2.12 × 10<sup>-3</sup></td><td>2.14 × 10<sup>-3</sup></td><td>2.14 × 10<sup>-3</sup></td></tr><tr><td>4</td><td>2.12 × 10<sup>-3</sup></td><td>2.14 × 10<sup>-3</sup></td><td>2.14 × 10<sup>-3</sup></td></tr><tr><td>5</td><td>2.12 × 10<sup>-3</sup></td><td>2.14 × 10<sup>-3</sup></td><td>2.14 × 10<sup>-3</sup></td></tr></tbody></table>

This table shows the computed values of $\alpha^*$ for different noise levels (0%, 1%, and 5%) and Taylor polynomial degrees L. It illustrates the effect of noise and L on the accuracy of $\alpha^*$ and the robustness of the MOR-T<sub>L</sub> strategy.

<table><caption>Table 5: Computational Time Ratio between MOR-T<sub>L</sub> and SEM for Line Search with 10 Parameter Values.</caption><thead><tr><th rowspan="2">L</th><th>MOR-T<sub>L</sub></th><th>Total Computational Time</th></tr><tr><th>SEM Computational Time</th></tr></thead><tbody><tr><td>0</td><td>0.12</td><td></td></tr><tr><td>1</td><td>0.27</td><td></td></tr><tr><td>2</td><td>0.39</td><td></td></tr><tr><td>3</td><td>0.55</td><td></td></tr><tr><td>4</td><td>0.74</td><td></td></tr><tr><td>5</td><td>1.16</td><td></td></tr></tbody></table>

This table shows the time ratio between MOR-T<sub>L</sub> and SEM during a line search with 10 parameter values $\alpha$, highlighting the computational cost efficiency of MOR-T<sub>L</sub> for different Taylor polynomial degrees L.

In conclusion, the results confirm that the standard MOR method is inadequate for line search due
to its poor performance in determining the optimal $\alpha^*$. The MOR-T$_L$ method, however, excels in both accuracy and computational efficiency, making it a strong candidate for practical FWI applications. Furthermore, the MOR-T$_L$ method remains robust even in the presence of noise, offering significant advantages over SEM, particularly in terms of computational feasibility.

# 6 Summary and Conclusion

In this work, we introduced a novel MOR approach, called MOR-T$_L$, based on Taylor polynomial expansions. The key innovation of MOR-T$_L$ lies in the construction of the ROM basis functions, which are generated using Taylor polynomial representations. This approach provides a more robust method for handling parameter variations effectively. Numerical experiments on a two-dimensional wave prototype problem, relevant to seismic applications, demonstrated the robustness of the MOR-T$_L$ basis to parameter changes and noise. Specifically, MOR-T$_L$ excels in both accuracy and computational efficiency, particularly in a line search minimization context, and can accommodate parameter perturbations of up to 20%, which in the presence of noise, offers substantial advantages over SEM in terms of computational feasibility. These results position MOR-T$_L$ as a promising candidate for practical applications in FWI, demonstrating its potential to improve computational cost-effectiveness and robustness.

# References

[1] M. BERGMANN AND L. CORDIER, *Optimal control of the cylinder wake in the laminar regime by trust-region methods and POD reduced-order models*, Journal of Computational Physics, 227 (2008), pp. 7813–7840.

[2] H. BÉRIOT AND G. GABARD, *Anisotropic adaptivity of the p-fem for time-harmonic acoustic wave propagation*, Journal of Computational Physics, 378 (2019), pp. 234–256.

[3] J. BESSET, *A Model Order Reduction Strategy for Parametrized PDEs: A New Paradigm for Efficient Subsurface Imaging.*, PhD thesis, Université de Pau et des Pays de l'Adour, June 2025.

[4] T. BUI-THANH, K. WILLCOX, AND O. GHATTAS, *Model Reduction for Large-Scale Systems with High-Dimensional Parametric Input Space*, SIAM Journal on Scientific Computing, 30 (2008), pp. 3270–3288. Publisher: Society for Industrial and Applied Mathematics.

[5] Prud'Homme, C., Rovas, D. V., Veroy, K., Machiels, L., Maday, Y., Patera, A. T., and Turinici, G. (2002). Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods. *J. Fluids Eng.*, 124(1):70–80.

[6] Buffa, A., Maday, Y., Patera, A. T., Prud'Homme, C., and Turinici, G. (2012). A priori convergence of the greedy algorithm for the parametrized reduced basis method. *ESAIM: Mathematical modelling and numerical analysis*, 46(3):595–603.

[7] F. CHINESTA, A. AMMAR, AND E. CUETO, *Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models*, Archives of Computational Methods in Engineering, 17 (2010), pp. 327–350.
[8] B. COCKBURN, G. E. KARNIADAKIS, AND C.-W. SHU, *The development of discontinuous galerkin methods*, in *Discontinuous Galerkin methods: theory, computation and applications*, Springer, 2000, pp. 3–50.

[9] —, *Discontinuous Galerkin methods: theory, computation and applications*, vol. 11, Springer Science & Business Media, 2012.

[10] L. CORDIER AND M. BERGMANN, *Réduction de dynamique par décomposition orthogonale aux valeurs propres (POD)*, 2006.

[11] B. DESPRÉS, M. EL GHAOUI, AND T. SAYAH, *A trefftz method with reconstruction of the normal derivative applied to elliptic equations*, Mathematics of Computation, 91 (2022), pp. 2645–2679.

[12] D. A. DI PIETRO AND A. ERN, *Mathematical aspects of discontinuous Galerkin methods*, vol. 69, Springer Science & Business Media, 2011.

[13] P. DÍEZ AND A. HUERTA, *A unified approach to remeshing strategies for finite element h-adaptivity*, Computer methods in applied mechanics and engineering, 176 (1999), pp. 215–229.

[14] M. DUMBSER AND M. KÄSER, *An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes*—ii. the three-dimensional isotropic case, Geophysical Journal International, 167 (2006), pp. 319–336.

[15] C. ECKART AND G. YOUNG, *The Approximation of One Matrix by Another of Lower Rank*, Psychometrika, 1 (1936), pp. 211–218.

[16] B. ENGQUIST AND A. MAJDA, *Absorbing boundary conditions for numerical simulation of waves*, Proceedings of the National Academy of Sciences, 74 (1977), pp. 1765–1766. Publisher: Proceedings of the National Academy of Sciences.

[17] L. C. EVANS AND R. GARIEPY, *Measure theory and fine properties of functions*, 1992.

[18] D. GIVOLI, *High-order local non-reflecting boundary conditions: a review*, Wave Motion, 39 (2004), pp. 319–326.

[19] G. H. GOLUB AND C. F. V. LOAN, *Matrix Computations*, Johns Hopkins University Press, Feb. 2013.

[20] R. HAWKINS, M. H. KHALID, M. SCHLOTTBOM, AND K. SMETANA, *Model order reduction for seismic applications*, June 2024. arXiv:2406.07207.

[21] I. HERRERA, *Trefftz method: a general theory*, Numerical Methods for Partial Differential Equations: An International Journal, 16 (2000), pp. 561–580.

[22] J. S. HESTHAVEN AND T. WARBURTON, *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications*, Springer Science & Business Media, 2007.

[23] Q. HU, D. GRANA, AND K. A. INNANEN, *Feasibility of seismic time-lapse monitoring of CO2 with rock physics parametrized full waveform inversion*, Geophysical Journal International, 233 (2022), pp. 402–419.
[24] D. KOMATITSCH AND J.-P. VILOTTE, *The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures*, Bulletin of the Seismological Society of America, 88 (1998), pp. 368-392.

[25] S. J. LEON, A. BJÖRCK, AND W. GANDER, *Gram-Schmidt orthogonalization: 100 years and more*, Numerical Linear Algebra with Applications, 20 (2013), pp. 492-532.

[26] R. M. MEZA, *Interpolation sur les variétés grassmanniennes et applications à la réduction de modèles en mécanique*, PhD thesis, Université de La Rochelle, 2018.

[27] A. T. MOHAN AND D. V. GAITONDE, *A Deep Learning based Approach to Reduced Order Modelling for Turbulent Flow Control using LSTM Neural Networks*, Apr. 2018. arXiv:1804.09269.

[28] L. MÉTIVIER, R. BROSSIER, J. VIRIEUX, AND S. OPERTO, *Full Waveform Inversion and the Truncated Newton Method*, SIAM Journal on Scientific Computing, 35 (2013), pp. B401-B437.

[29] A. T. PATERA, *A spectral element method for fluid dynamics: Laminar flow in a channel expansion*, Journal of Computational Physics, 54 (1984), pp. 468-488.

[30] G. STRANG, *Introduction to Linear Algebra*, Wellesley, Aug. 2016.

[31] L. N. TREFETHEN AND D. B. III, *Numerical Linear Algebra*, Society for Industrial and Applied Mathematics, Philadelphia, Pa, 1997.

[32] J. H. TU, C. W. ROWLEY, D. M. LUCHTENBURG, S. L. BRUNTON, AND J. N. KUTZ, *On dynamic mode decomposition: Theory and applications*, Journal of Computational Dynamics, 1 (2014), pp. 391-421.

[33] E. TURKEL, *Boundary conditions and iterative schemes for the Helmholtz equation in unbounded regions*, Computational Methods for Acoustics Problems, (2008), pp. 127-158.

[34] J. VIRIEUX AND S. OPERTO, *An overview of full-waveform inversion in exploration geophysics*, GEOPHYSICS, 74 (2009), pp. WCC1-WCC26. Publisher: Society of Exploration Geophysicists.

[35] O. ZIENKIEWICZ AND J. ZHU, *Adaptivity and mesh generation*, International Journal for Numerical Methods in Engineering, 32 (1991), pp. 783-810.

[36] R. ŞTEFĂNESCU AND I. M. NAVON, *POD/DEIM nonlinear model order reduction of an ADI implicit shallow water equations model*, Journal of Computational Physics, 237 (2013), pp. 95-114.
Inria

RESEARCH CENTRE
BORDEAUX - SUD-OUEST

200 avenue de la Vieille Tour
33405 Talence Cedex

Publisher
Inria
Domaine de Voluceau - Rocquencourt
BP 105 - 78153 Le Chesnay Cedex
inria.fr

ISSN 0249-6399