# Modélisation avancée des tsunamis : une approche spectrale pour intégrer la dynamique sismique

Thomas Melkior, Harsha S. Bhat, Faisal Amlani

► To cite this version:

Thomas Melkior, Harsha S. Bhat, Faisal Amlani. Modélisation avancée des tsunamis: une approche spectrale pour intégrer la dynamique sismique. 26e Congrès Français de Mécanique (CFM 2025), Association Française de Mécanique (AFM), Aug 2025, Metz, France. hal-05229546

HAL Id: hal-05229546

https://hal.science/hal-05229546v1

Submitted on 29 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
# Modélisation avancée des tsunamis : une approche spectrale pour intégrer la dynamique sismique

T. MELKIORª, H. S. BHATᵇ, F. AMLANIᶜ

a. Laboratoire de Mécanique Paris-Saclay; Ecole Normale Supérieure Paris-Saclay; Centrale Supélec; CNRS (thomas.melkior@ens-paris-saclay.fr)

b. Laboratoire de géologie de l'ENS; Ecole Normale Supérieure de Paris; CNRS (harshasbhat@gmail.com)

c. Laboratoire de Mécanique Paris-Saclay; Ecole Normale Supérieure Paris-Saclay; Centrale Supélec; CNRS (faisal.amlani@ens-paris-saclay.fr)

## Résumé

Malgré les progrès en modélisation des tsunamis, certains événements historiques restent mal compris [1]. Notamment, la dynamique du mouvement sismique est souvent négligée, les perturbations étant modélisées par un déplacement instantané du fond marin. Ce travail explore l'impact de cette dynamique sur la génération de tsunamis, en tenant compte des échelles distinctes entre les ondes sismiques (haute fréquence) et les vagues de tsunami (basse fréquence). Pour cela, un nouveau solveur basé sur la méthode spectrale de Continuation de Fourier (FC) [2] est proposé. Cette méthode d'ordre élevé élimine les erreurs de dispersion numérique, offrant une précision accrue pour la propagation à longue distance par rapport aux méthodes classiques (différences finies, éléments finis, etc.), qui nécessitent un raffinement coûteux. Le modèle repose sur les équations non linéaires des eaux peu profondes, modifiées pour intégrer des fonds marins en mouvement. Deux validations sont proposées : un cas analytical visant à illustrer l'absence de dispersion numérique et un cas expérimental inspiré des travaux de Hammack, où un piston déforme le fond marin [3]. Les résultats montrent une excellente précision, validant l'approche pour la modélisation réaliste des tsunamis d'origine sismique.

## Abstract :

Despite significant progress in tsunami modeling, certain historical events remain poorly understood [1]. In particular, the dynamic nature of seismic motion is often neglected, with the perturbation typically represented as an instantaneous vertical displacement of the seafloor. This work investigates the impact of seismic dynamics on tsunami generation, taking into account the contrasting scales between seismic waves (high-frequency) and tsunami waves (low-frequency). To address this, a new solver based on the spectral Fourier Continuation (FC) method [2] is proposed. This high-order method eliminates numerical dispersion errors, offering enhanced accuracy for long-distance wave propagation when compared to conventional methods (finite differences, finite elements, etc.), which require substantial mesh refinement to mitigate dispersion. The governing model is based on the nonlinear shallow water equations, modified to include time-dependent seafloor motion. Two validation cases are presented: an analytical test case demonstrating the absence of numerical dispersion, and an experimental scenario inspired by the work of Hammack, in which a piston induces vertical seabed motion [3]. The results
exhibit excellent accuracy, supporting the suitability of the proposed approach for realistic modeling of earthquake-generated tsunamis.

# Mots clefs : Tsunami, Simulation numérique, Dynamique sismique, Benchmark, Méthode spectrale.

## 1 Introduction

Afin de mieux comprendre les mécanismes de génération et de propagation des tsunamis pour une évaluation rigoureuse des risques, la simulation numérique est un outil essentiel. Des événements comme le tsunami atypique de Palu (2018) [1] soulignent l'importance de tels efforts de modélisation.

La propagation des tsunamis est traditionnellement étudiée à deux échelles : les zones proches des côtes (near-field), dominées par de fortes non-linéarités, et les zones lointaines (far-field), où certains phénomènes non-linéaires deviennent négligeables grâce à la profondeur accrue [4]. Les solveurs complets reposent souvent sur des maillages emboîtés, avec des méthodes numériques adaptées à chaque zone : les volumes finis [5], sont couramment utilisés pour la zone proche, et les différences finies (FD) [4] pour les zones lointaines. Cependant, ces dernières approches souffrent d'erreurs liées à la dispersion numérique qui peuvent s'accumuler significativement lors de la propagation des ondes dans l'espace et le temps, nécessitant un raffinement coûteux pouvant devenir prohibitif [6].

Nous proposons un nouveau solveur 2D non linéaire des équations des eaux peu profondes, basé sur la méthode Fourier Continuation (FC) [2, 7]. Cette méthode spectrale permet une résolution précise des équations aux dérivées partielles (EDP), tout en rendant quasi-inexistante la dispersion numérique, ce qui en fait un choix prometteur pour la modélisation de tsunamis sismogéniques. De plus, sa précision d'ordre élevé facilite la capture des multiples échelles spatiales et temporelles en jeu. Ce travail introduit également une extension des équations classiques permettant d'inclure le mouvement dynamique du fond marin [8]. Alors que les modèles classiques supposent un déplacement instantané, des études récentes ont montré que la dynamique du fond marin peut être déterminante lorsque la vitesse de propagation du tsunami n'est pas négligeable devant la vitesse de rupture sismique [1, 9]. Cette dernière introduit des défis numériques liés au caractère multi-échelle du problème.

Ce travail prolonge une preuve de concept 1D [1] appliquée au tsunami de Palu, en développant un solveur 2D d'ordre élevé, rapide (comparable à la FFT) et quasi exempt de dispersion. Après la présentation du système d'équations (Section 2) et de la méthode numérique (Section 3), nous illustrons les performances du solveur par un cas test mettant en lumière l'absence de dispersion ainsi qu'un benchmark basé sur des données expérimentales de fond en mouvement (Section 4).

## 2 Modélisation de tsunami avec déplacement du fond dans le mod-èle des eaux peu profondes

Le nouveau solveur d'ordre élevé présenté dans cette étude repose sur les équations des eaux peu profondes de Saint-Venant [10], une version moyennée en profondeur des équations d'Euler. La validité de ces équations pour modéliser les écoulements fluides repose sur plusieurs hypothèses fondamentales [11] : fluide incompressible, pression hydrostatique ainsi qu'une profondeur $h(x, y, t)$ faible comparée à
la longueur d'onde des vagues. Les composantes horizontales du champ de vitesses $u(x, y, t)$ (selon $x$) et $v(x, y, t)$ (selon $y$) sont supposées constantes selon la verticale $z$. Dans cette contribution, le système est formulé dans un repère cartésien — approximation valable tant que les distances de propagation restent raisonnables par rapport à la courbure terrestre.

La configuration est illustrée Figure 1, où la hauteur totale d'eau s'écrit $h(x, y, t) = h_0(x, y) + \eta(x, y, t)$ — $\xi(x, y, t)$. Ici, $h_0(x, y)$ désigne la profondeur au repos (bathymétrie), $\eta(x, y, t)$ est la surélévation de la surface libre, et $\xi(x, y, t)$ représente le déplacement vertical du fond marin (perturbation sismique). L'objectif de ce travail est de développer un solveur capable de prendre en compte l'évolution temporelle de la source sismique afin de modéliser les effets dynamiques du séisme sur la génération de tsunamis, représentés par le déplacement $\xi(x, y, t)$ et la vitesse $\xi_t(x, y, t)^1$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Définition des grandeurs modélisant les eaux peu profondes.</figcaption></figure>

Dans la version unidimensionnelle (1D) de la configuration illustrée à la Figure 1, la dimension transverse $y$ est supposée invariante : les quantités ne dépendent pas de $y$, et la composante de vitesse correspondante est donc nulle, i.e. $v = 0$. Le système 1D utilisé dans cette étude est donné par [1]

$$ \begin{cases} h_t + (hu)_x = 0, \\ (hu)_t + (hu^2)_x + gh\eta_x = 0, \end{cases} \qquad (1) $$

où $g$ désigne l'accélération gravitationnelle (prise égale à 9,81 m/s²), $h = h(x, t) = h_0(x) + \eta(x, t)$ — $\xi(x, t)$, et où les inconnues sont la vitesse $u = u(x, t)$ et la surface libre $\eta = \eta(x, t)$. En développant, réorganisant, et en supposant $h > 0$ (hypothèse usuelle [12]), on obtient

$$ \begin{cases} \eta_t - \xi_t + (hu)_x = 0, \\ u_t + uu_x + g\eta_x = 0, \\ h(x, t) = h_0(x) + \eta(x, t) - \xi(x, t). \end{cases} \qquad (2) $$

<footnote n="1">¹Par souci de concision, les dérivées partielles selon $x$, $y$ ou $t$ seront notées en indice dans l'ensemble de ce manuscrit, par exemple : $\xi_t := \partial\xi/\partial t$.</footnote>
Sous forme vectorielle, en posant $\mathbf{U} = [\eta, u]^T$, on peut réécrire le système comme

$$ \mathbf{U}_t + \mathbf{F}(\mathbf{U}, \mathbf{U}_x) = \mathbf{S}, \qquad (3) $$

où

$$ \mathbf{F}(\mathbf{U}, \mathbf{U}_x) = \begin{bmatrix} h_x u + h u_x \\ u u_x + g \eta_x \end{bmatrix} = \begin{bmatrix} (\eta_x + h_{0x} - \xi_x)u + (\eta + h_{0} - \xi)u_x \\ u u_x + g \eta_x \end{bmatrix} \quad \text{et} \quad \mathbf{S} = \begin{bmatrix} \xi_t \\ 0 \end{bmatrix}. $$

Le caractère bien posé de ce système est complété par des conditions aux limites.

Le système bidimensionnel (2D) correspondant est donné par [13, 14, 15]

$$ \left\{ \begin{array}{l} h_t + (hu)_x + (hv)_y = 0, \\ u_t + uu_x + vu_y + g\eta_x = 0, \\ v_t + uv_x + vv_y + g\eta_y = 0, \\ h = h_0 + \eta - \xi, \end{array} \right. \quad \text{ou,} \quad \left\{ \begin{array}{l} \eta_t - \xi_t + (hu)_x + (hv)_y = 0 \\ u_t + uu_x + vu_y + g\eta_x = 0 \\ v_t + uv_x + vv_y + g\eta_y = 0 \\ h = h_0 + \eta - \xi \end{array} \right. \qquad (4) $$

avec $u = u(x, y, t)$, $h_0 = h_0(x, y)$, $\eta = \eta(x, y, t)$ et $\xi = \xi(x, y, t)$.

Sous forme vectorielle, en posant $\mathbf{U} = [\eta, u, v]^T$, on obtient le système suivant :

$$ \mathbf{U}_t + \mathbf{F}(\mathbf{U}, \nabla \mathbf{U}) = \mathbf{S} \qquad (5) $$

où

$$ \mathbf{F}(\mathbf{U}, \nabla \mathbf{U}) = \begin{bmatrix} (\eta_x + h_{0x} - \xi_x)u + (\eta + h_0 - \xi)u_x + (\eta_y + h_{0y} - \xi_y)v + (\eta + h_0 - \xi)v_y \\ uu_x + vu_y + g\eta_x \\ uv_x + vv_y + g\eta_y \end{bmatrix} $$

et

$$ \mathbf{S} = \begin{bmatrix} \xi_t \\ 0 \\ 0 \end{bmatrix}. $$

Comme pour le système 1D, le caractère bien posé du système 2D est assuré par des conditions aux limites. Dans le cadre de cette étude, la profondeur au repos $h_0(x, y)$ ainsi que la perturbation du fond marin $\xi(x, y, t)$ sont supposées connues a priori (cette dernière pouvant provenir, par exemple, de simulations de rupture sismique ou de données géophysiques).

# 3 Méthodologie

Le solveur présenté repose sur Fourier Continuation pour le calcul des dérivées spatiales, ainsi qu'un schéma temporel explicite. Fourier Continuation est présenté dans la section qui suit avant de décrire le solveur complet.

## 3.1 Fourier Continuation

Soit $N \in \mathbb{N}^*$ le nombre de points utilisés pour discrétiser une fonction régulière (éventuellement non périodique) $f: [0, 1] \to \mathbb{R}$. On considère sans perte de généralité un intervalle unitaire (les résultats
pouvant être étendus à tout intervalle par une simple transformation affine). En notant $x_j$ les points de grille correspondant à une discrétisation uniforme de l'intervalle $[0, 1]$, i.e. $x_j = j\Delta x$, $0 \le j \le N - 1$ avec $\Delta x = 1/(N - 1)$, la méthode Fourier Continuation (FC) consiste à prolonger $f$ en une fonction $f^c$ qui est périodique sur un intervalle étendu $[0, b]$, avec $b > 1$, laquelle s'écrit sous forme d'une série de Fourier à bande limitée :

$$f^c: \begin{cases} [0, b] \to \mathbb{C} \\ x \mapsto \sum_{k=-M}^{M} a_k \exp \left( \frac{2\pi ikx}{b} \right), \end{cases} \qquad (6)$$

où $M = (N + C)/2$, avec $C$ le nombre de points d'extension qui définissent l'intervalle $b$, c'est-à-dire $b = (N + C - 1)\Delta x$. Les coefficients $a_k$, $k = -M, \dots, M$ sont déterminés en minimisant l'erreur discrète entre $f^c$ et la fonction originale $f$ aux points numériques $\mathbf{x} = [x_0, \dots, x_{N-1}]^T$, c'est-à-dire :

$$(a_k)_{k \in \{-M, \dots, M\}} = \arg \min \|f^c(\mathbf{x}) - f(\mathbf{x})\|_2. \qquad (7)$$

Autrement dit, $f^c$ approxime une interpolation trigonométrique de $f$ sur ses points de discrétisation. La Figure 2 illustre l'extension correspondante obtenue à partir d'une fonction discrétisée initiale $f = [f(x_0), \dots, f(x_{N-1})]^T$, où $f_{\text{ext}}^c$ représente les $C$ valeurs interpolées dans l'intervalle $[1, b]$. La représentation discrète complète de la fonction de continuation $f^c$ est alors donnée par

$$\mathbf{f}^c = \begin{bmatrix} \mathbf{f} \\ \mathbf{f}_{\text{ext}}^c \end{bmatrix} = \begin{bmatrix} f(x_0) \\ \vdots \\ f(x_{N-1}) \\ f_{\text{ext}}^c(x_N) \\ \vdots \\ f_{\text{ext}}^c(x_{N+C-1}) \end{bmatrix}.$$

Une fois les coefficients de la série de Fourier $f^c$ (Équation (6)) déterminés, les dérivées spatiales correspondantes peuvent être obtenues analytiquement par dérivation terme à terme, c'est-à-dire

$$(f^c)'(x) = \sum_{k=-M}^{M} \frac{2\pi ik}{b} a_k \exp\left(\frac{2\pi ikx}{b}\right). \qquad (8)$$

Ainsi, on obtient une approximation de la dérivée de $f$ en restreignant $f^c$ à l'intervalle original, soit

$$f'(x) \approx (f^c)'(x) \quad \forall x \in [0, 1]. \qquad (9)$$

Une propriété intéressante des méthodes spectrales est que si $f$ est infiniment différentiable ($C^\infty$), l'erreur de troncature d'une telle approximation converge vers zéro plus rapidement que n'importe quelle puissance finie de $\Delta x$ [6].

L'Équation (7) peut être résolue par exemple à l'aide d'une décomposition en valeurs singulières (SVD). Cependant, la complexité computationnelle associée est de l'ordre $O(N^3)$, ce qui est trop coûteux pour des problèmes dynamiques nécessitant l'application d'une telle opération à chaque pas de temps. Pour résoudre cette difficulté, une méthode accélérée appelée FC(Gram) décrite dans [2, 16] est utilisée dans
<figure><img src="image_3.png" /><figcaption>Figure 2: Schéma de la continuation de Fourier (extension périodique d'une fonction).</figcaption></figure>

la présente étude. Cette version accélérée, dorénavant désignée par FC, permet de réduire la complexité à $O(N \log(N))$. Pour la suite de l'étude nous considérons la version accélérée FC d'ordre 5.

Cette méthode peut se généraliser à deux dimensions en considérant le domaine 2D comme un ensemble de tranches 1D, lesquelles peuvent être successivement prolongées de manière périodique par la procédure décrite afin de calculer les dérivées partielles.

## 3.2 Le solveur

Le solveur développé résout le système des eaux peu profondes à l'aide d'un schéma temporel explicite, au sein duquel les dérivées spatiales du terme d'évolution sont calculées par FC sous sa version accélérée. Afin de rester cohérent avec la haute précision et l'absence de dispersion de la méthode FC — nécessaires pour modéliser les tsunamis d'origine sismique — un schéma temporel d'ordre élevé est utilisé dans le solveur. Bien que le schéma de Runge-Kutta d'ordre quatre (RK4) soit un choix courant, un choix plus avantageux en termes de rapidité de calcul est le schéma d'Adams-Bashforth d'ordre quatre (AB4), dont les régions de stabilité sont comparables [17]. Ce schéma explicite d'ordre 4 permet une bonne précision temporelle tout en maintenant une charge de calcul réduite (une seule évaluation du terme d'évolution de l'EDP), ce qui est particulièrement avantageux pour des simulations longues ou en grande dimension. Étant donné que l'utilisation du schéma AB4 nécessite la connaissance des trois pas de temps précédant le pas courant, le solveur des équations d'eau peu profonde proposé est initialisé en calculant les trois premiers pas de temps à l'aide du schéma RK4.

Le terme d'évolution de l'EDP pour le système des équations d'eau peu profonde tel que défini dans l'Équation (5), est donné par :

$$ \mathbf{rhs} = -\mathbf{F}(\mathbf{U}, \nabla \mathbf{U}) + \mathbf{S}. \qquad (10) $$

Le schéma temporel consiste alors à calculer un pas de temps $U^{n+1}$ en utilisant la relation $U_t = \mathbf{rhs}$ et les pas de temps $U^n, U^{n-1}, U^{n-2}$ et $U^{n-3}$. Le calcul du vecteur rhs (right handside) implique le calcul des dérivées spatiales, ce qui est réalisé à l'aide de la méthode FC.
# 4 Résultats

Nous allons maintenant présenter une étude de la dispersion du solveur, pour montrer que celui-ci hérite des propriétés intéressantes de FC, ainsi qu'une évaluation des performances du solveur sur un benchmark classique mettant en jeu un mouvement dynamique du fond.

## 4.1 Traitement de la dispersion

On considère ici le système 1D, les conclusions de cette section pouvant être étendues au cas 2D. Les erreurs de dispersion numérique (ou erreurs de pollution) du solveur basé sur la méthode FC peuvent être comparées à celles d'autres méthodes en considérant une solution forcée (solution manufacturée, voir [18]) sous la forme d'une onde sinusoidale en propagation, donnée par :

$$ \left\{ \begin{array}{l} \eta_{\text{ref}}(x, t) = \sin(2\pi(x - t)), \\ u_{\text{ref}}(x, t) = \cos(2\pi(x - t)), \end{array} \right. \qquad (11) $$

sur les intervalles espace-temps correspondant à $x \in X_n = [0, n - 0.1]$ et $t \in T_n = [0, n]$ pour un nombre $n \in \mathbb{N}^*$ de longueurs d'onde dans le domaine spatial, chaque longueur d'onde étant discrétisée avec un nombre $N_\lambda$ de points par longueur d'onde. Le décalage de 0.1 dans l'intervalle spatial est intentionnel afin de produire une solution non périodique sur le domaine, rendant ainsi la procédure FC pleinement pertinente. La profondeur d'eau au repos est prise constante et égale à $h_0(x) = 5$ sur tout le domaine, et aucun mouvement du fond marin n'est considéré (c'est-à-dire que $\xi(x, t) = 0$ à tout instant).

Lorsque $n$ augmente, la solution se propage sur une distance plus longue. Pour les tests qui suivent, la résolution par longueur d'onde est fixe : $N_\lambda = 30$ points, ce qui correspond à une discrétisation spatiale de $\Delta x = 1/N_\lambda$. Le pas de temps est alors ajusté pour chaque taille de problème conformément à la condition de CFL du solveur. La Figure 3 présente les erreurs en norme $L^\infty$ pour $n$ allant jusqu'à $n = 40$ longueurs d'onde, en comparant le solveur basé sur la méthode FC et le même solveur auquel on substitue FC par des différences finies d'ordre 4 et 6. Ces erreurs sont calculées sur l'ensemble de l'espace et du temps selon l'expression suivante (ici pour $\eta$) :

$$ \text{erreur } L^\infty = \frac{\max_{x,t} (|\eta(x,t) - \eta_{\text{ref}}(x,t)|)}{\max_{x,t} (|\eta_{\text{ref}}(x,t)|)} \qquad (12) $$

Il est clair que l'erreur associée à la méthode par différences finies (FD) augmente avec la taille de $n$, et ce même en augmentant l'ordre du schéma FD de quatrième à sixième ordre (un ordre plus élevé que celui utilisé dans la méthode FC employée, à savoir cinq, et bien plus élevé que les ordres couramment utilisés dans les solveurs tsunami). En effet, les méthodes FD et FE sont bien connues pour être sujettes à des erreurs de dispersion numérique (pollution) [6, 19], nécessitant un raffinement substantiel des maillages associés afin de capturer la solution avec une précision donnée sur des durées et distances de propagation importantes.

En revanche, l'erreur issue du solveur basé sur la méthode FC demeure constante, ce qui indique une absence totale de dispersion numérique. Cela signifie également que le niveau d'erreur est uniquement déterminé par le nombre de points par longueur d'onde, ce qui constitue un avantage considérable pour la résolution de problèmes de propagation d'ondes.
<figure><img src="image_4.png" /><figcaption>Figure 3: Courbes de dispersion des méthodes FD4, FD6 et FC pour l'équation 1D des eaux peu profondes avec la solution forcée $\eta_{\text{ref}}(x, t) = \sin(2\pi(x - t))$ ; $u_{\text{ref}}(x, t) = \cos(2\pi(x - t))$ sur les intervalles (variables) $X_n = [0, n - 0.1]$ et $T_n = [0, n]$. Le nombre de points par longueur d'onde est $N_{\lambda} = 30$.</figcaption></figure>

## 4.2 Génération de vague par brusque mouvement vertical du fond

Dans cette section, le solveur est confronté à la prise en compte de mouvements dynamiques du fond marin. Une étude expérimentale réalisée par Hammack et al. [3] a mesuré en laboratoire l'élévation de la surface libre de l'eau générée par un piston déplaçant verticalement une partie du fond marin. Un schéma de la situation est donné dans la Figure 4.

La bathymétrie est supposée plane, la perturbation dynamique étant donnée par la loi suivante :

$$ \xi(x, t) = f(t)\mathcal{H}(x^2 - x_1^2) \quad (13) $$

où $\mathcal{H}$ désigne la fonction de Heaviside:

$$ \mathcal{H}(x) = \begin{cases} 1 & \text{si } x \ge 0, \\ 0 & \text{si } x < 0. \end{cases} \quad (14) $$

Étant donné que le solveur utilise une méthode spectrale, ce rectangle a été lissé en approximant la fonction de Heaviside $\mathcal{H}$ par des tangentes hyperboliques. En réalité, l'utilisation de la fonction de Heaviside n'empêche pas le solveur de donner une solution précise, mais elle pourrait poser un problème pour la convergence. L'évolution temporelle est définie par

$$ f(t) = \begin{cases} 0 & \text{si } t < 0, \\ \xi_m (1 - \cos(\pi t / t_r)) / 2 & \text{si } 0 \le t \le t_r, \\ \xi_m & \text{si } t > t_r. \end{cases} \quad (15) $$

Les équations considérées dans le schéma sont celles du système 1D, mais les calculs ont également été réalisés en utilisant le modèle 2D afin de tester le comportement du solveur 2D dans ce cas purement 1D (sans dépendance en $y$). Dans l'étude précédente [3], l'auteur a déterminé la solution du problème
<figure><img src="image_5.png" /><figcaption>Figure 4: Schéma du problème inspiré par Hammack [3].</figcaption></figure>

<table><caption>Table 1: Paramètres du benchmark.</caption><thead><tr><th>h<sub>0</sub><br>[m]</th><th>g<br>[m.s<sup>-2</sup>]</th><th>x<sub>1</sub>/h<br>[-]</th><th>τ<br>[-]</th><th>ξ<sub>m</sub>/h<br>[-]</th><th>x<sub>min</sub>/x<sub>1</sub><br>[-]</th><th>x<sub>max</sub>/x<sub>1</sub><br>[-]</th><th>y<sub>min</sub>/x<sub>1</sub><br>[-]</th><th>y<sub>max</sub>/x<sub>1</sub><br>[-]</th><th>t<sub>max</sub>√g/h<br>[-]</th></tr></thead><tbody><tr><td>200</td><td>9.81</td><td>12.2</td><td>0.11</td><td>0.4</td><td>-10</td><td>10</td><td>-0.03</td><td>0.03</td><td>70</td></tr></tbody></table>

linéarisé sous la forme suivante :

$$ \eta(x, t) = \frac{\xi_m}{\pi} \int_0^\infty \frac{\kappa^2 \cos(kx) \sin(kb)}{k \cosh(kh) (\kappa^2 - \omega^2)} [\cos(\omega t) - \cos(\kappa t) + \mathcal{H}(t - t_r)(\cos(\omega(t - t_r))) + \cos(\kappa t)] dk \quad (16) $$

$$ \text{Où } t_c = x_1/\sqrt{gh}, t_r/t_c = \tau, \kappa = \pi/t_r \text{ et } \omega = \sqrt{gk \tanh(kh)}. $$

Les paramètres de calcul sont résumés dans le Tableau 1 et ont été inspirés par une étude de Derakhti et al. qui a réussi à reproduire les résultats expérimentaux [20]. Ces paramètres sont de plusieurs ordres de grandeur différents par rapport à l'expérience réelle (la profondeur ici est de 200 m, alors que dans l'installation expérimentale, la profondeur variait entre 5 cm et 50 cm), par conséquent, les résultats seront tracés en considérant des valeurs normalisées.

Les données de laboratoire disponibles étant localisées au centre du piston et à l'extrémité du piston, les résultats du solveur ont été tracés à ces emplacements. Les courbes correspondant au centre du piston et à l'extrémité du piston sont présentées dans la Figure 5 (respectivement dans la Fig. 5a et la Fig. 5b). Les résultats du solveur s'accordent bien avec les points de données sur les deux courbes, avec une erreur $L^\infty$ de 11 % pour la courbe au centre du piston et de 14 % pour la courbe à l'extrémité du piston. La figure montre que la théorie linéaire, qui évolue comme un système oscillant amorti, ne modélise pas correctement la réponse dynamique. L'étude de Derakhti et al. [20] explique que cette théorie linéaire aurait correspondu aux résultats expérimentaux pour $\xi_m \ll h$, ce qui n'est pas le cas ici étant donné que $\xi_m/h = 0.4$. Le solveur permet donc de capturer les non-linéarités de ce mouvement.

# 5 Conclusion

Dans ce travail, nous développons et validons un solveur spectral d'ordre élevé basé sur la méthode Fourier Continuation (FC) pour les équations non linéaires d'eaux peu profondes avec fond marin dynamique. L'absence des erreurs de dispersion numérique en fait un outil adapté à la modélisation à longue distance des tsunamis. Comparée aux méthodes classiques par différences ou éléments finis,
<figure><img src="image_6.png" /><figcaption>Figure 5: Évolution de la surface libre normalisée η au centre et au bord du piston. Les résultats expérimentaux et la solution analytique de Hammack (1973) [3] sont comparés aux simulations obtenues avec le solveur pour deux tailles de maillage, montrant une convergence vers les données expérimentales.</figcaption></figure>

l'approche FC conserve fidèlement les formes d'ondes sur de grandes échelles spatio-temporelles sans surcoût numérique majeur — un atout pour l'alerte précoce et l'évaluation des risques. Le solveur a été validé sur des données expérimentales, confirmant sa capacité à capturer les effets physiques liés à la génération et à la propagation des tsunamis. Une hybridation entre solveurs spectraux et techniques d'apprentissage automatique (de type PINNs) représente une piste prometteuse pour concilier vitesse et fidélité. Les perspectives incluent l'intégration de modèles de rupture sismique dynamiques pour une simulation couplée séisme-tsunami à grande échelle, ainsi que l'utilisation de raffinement de maillage adaptatif (AMR) et d'architectures de calcul haute performance (HPC), notamment via l'accélération GPU.

## References

[1] F. Amlani, H. S. Bhat, W. J. F. Simons, A. Schubnel, C. Vigny, A. J. Rosakis, J. Efendi, A. E. Elbanna, P. Dubernet, H. Z. Abidin, Supershear shock front contribution to the tsunami from the 2018 Mw 7.5 Palu, Indonesia earthquake, Geophysical Journal International 230 (3) (2022) 2089–2097.

[2] F. Amlani, O. P. Bruno, An FC-based spectral solver for elastodynamic problems in general three-dimensional domains, Journal of Computational Physics 307 (2016) 333–354.

[3] J. L. Hammack, A note on tsunamis: their generation and propagation in an ocean of uniform depth, Journal of Fluid Mechanics 60 (4) (1973) 769–799.

[4] F. Imamura, A. C. Yalciner, G. Ozyurt, Tsunami modelling manual (TUNAMI model) (1995).

[5] R. J. LeVeque, D. L. George, M. J. Berger, Tsunami modelling with adaptively refined finite volume methods, Acta Numerica 20 (2011) 211–289.

[6] D. Durran, Numerical Methods for Wave Equations in Geophysical Fluid Dynamics, Texts in Applied Mathematics, ISSN 0939-2475, Springer New York, NY, 1999.
[7] O. P. Bruno, M. Lyon, High-order unconditionally stable fc-ad solvers for general smooth domains i. basic elements, Journal of Computational Physics 229 (6) (2010) 2009–2033.

[8] D. Dutykh, D. Clamond, Modified shallow water equations for significantly varying seabeds, Applied mathematical modelling 40 (23-24) (2016) 9767–9787.

[9] L. S. Abrahams, L. Krenz, E. M. Dunham, A.-A. Gabriel, T. Saito, Comparison of methods for coupled earthquake and tsunami modelling, Geophysical Journal International 234 (1) (2023) 404–426.

[10] B. De St Venant, Théorie du mouvement non permanent des eaux, avec application aux crues des rivières et à l'introduction de marées dans leurs lits, Comptes rendus hebdomadaires des séances de l'Académie des sciences 73 (99) (1871) 148–154.

[11] A. Kundu, Tsunami and nonlinear waves, Vol. 316, Springer, 2007.

[12] D. Dutykh, D. Mitsotakis, X. Gardeil, F. Dias, On the use of the finite fault solution for tsunami generation problems, Theoretical and Computational Fluid Dynamics 27 (2013) 177–199.

[13] A. Delis, T. Katsaounis, Numerical solution of the two-dimensional shallow water equations by the application of relaxation methods, Applied Mathematical Modelling 29 (8) (2005) 754–783.

[14] Y. Wei, X.-Z. Mao, K. F. Cheung, Well-balanced finite-volume model for long-wave runup, Journal of Waterway, Port, Coastal, and Ocean Engineering 132 (2) (2006) 114–124.

[15] F. Dias, D. Dutykh, Dynamics of tsunami waves, in: Extreme man-made and natural hazards in dynamics of structures, Springer, 2007, pp. 201–224.

[16] F. Amlani, H. Wei, N. M. Pahlevan, A fourier-based methodology without numerical diffusion for conducting dye simulations and particle residence time calculations, Journal of Computational Physics 493 (2023) 112472.

[17] J.-P. Demailly, Analyse numérique et équations différentielles., Collection Grenoble sciences, ISSN 0767-371X, Les Ulis - EDP Sciences, 2016.

[18] P. J. Roache, Code Verification by the Method of Manufactured Solutions, Journal of Fluids Engineering 124 (1) (2001) 4–10.

[19] A. Deraemaeker, I. Babuška, P. Bouillard, Dispersion and pollution of the fem solution for the helmholtz equation in one, two and three dimensions, International journal for numerical methods in engineering 46 (4) (1999) 471–499.

[20] M. Derakhti, R. A. Dalrymple, E. A. Okal, C. E. Synolakis, Temporal and topographic source effects on tsunami generation, Journal of Geophysical Research: Oceans 124 (2019) 5270–5288.