# Coupling Dispersive Shallow Water Models by Deriving Asymptotic Interface Operators

José Galaz, Maria Kazolea, Antoine Rousseau

► To cite this version:

José Galaz, Maria Kazolea, Antoine Rousseau. Coupling Dispersive Shallow Water Models by Deriving Asymptotic Interface Operators. 27th International Conference on Domain Decomposition Methods in Science and Engineering - DD27, Jul 2022, Prague, Czech Republic. pp.181-188, 10.1007/978-3-031-50769-4_21. hal-03851031v3

HAL Id: hal-03851031

https://inria.hal.science/hal-03851031v3

Submitted on 21 Feb 2024

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
# Coupling Dispersive Shallow Water Models by Deriving Asymptotic Interface Operators

José Galaz, Maria Kazolea, Antoine Rousseau

## 1 Introduction

We are interested in coupling the linear Green-Naghdi equations (LGNE)

$$ \partial_t \zeta + \partial_x V = 0 \qquad (1) $$

$$ \partial_t V + \partial_x \zeta = \phi \qquad (2) $$

$$ -\frac{\mu}{3} \partial_x^2 \phi + \phi = -\frac{\mu}{3} \partial_x^3 \zeta \qquad (3) $$

for $x < 0$ and $t \in (0, T)$, with the linear shallow water equations (LSWE)

$$ \partial_t \zeta + \partial_x V = 0 \qquad (4) $$

$$ \partial_t V + \partial_x \zeta = 0 \qquad (5) $$

for $x > 0$ and $t \in (0, T)$ with $T > 0$, to represent the 1D propagation of water waves in shallow water. Here $\partial_t, \partial_x$ denote partial derivatives in the time and space variables $t, x$; $V(t, x)$ and $\zeta(t, x)$ stand for the vertically-averaged velocity and the free-surface level over its state at rest; $\phi(t, x)$ is an auxiliary variable for the elliptic part of the problem; and $\mu > 0$ is the asymptotic parameter characterizing the wave dispersion.

José Galaz
INRIA, Team LEMON, Centre Inria d'Université Côte d'Azur, Antenne Montpellier, Bat 5 CC05
017, 860 rue Saint-Priest, 34095 Montpellier Cedex 5 France, e-mail: jose.galaz@inria.fr

Maria Kazolea
INRIA, Team CARDAMOM, INRIA Bordeaux-Sud-Ouest, 200 Avenue de la Vieille Tour, 33405
Talence cedex France e-mail: maria.kazolea@inria.fr

Antoine Rousseau
INRIA, Team LEMON, Centre Inria d'Université Côte d'Azur, Antenne Montpellier, Bat 5 CC05
017, 860 rue Saint-Priest, 34095 Montpellier Cedex 5 France, e-mail: antoine.rousseau@inria.fr
# Coupling Dispersive Shallow Water Models by Deriving Asymptotic Interface Operators

José Galaz, Maria Kazolea, Antoine Rousseau

## 1 Introduction

We are interested in coupling the linear Green-Naghdi equations (LGNE)

$$ \partial_t \zeta + \partial_x V = 0 \qquad (1) $$

$$ \partial_t V + \partial_x \zeta = \phi \qquad (2) $$

$$ -\frac{\mu}{3} \partial_x^2 \phi + \phi = -\frac{\mu}{3} \partial_x^3 \zeta \qquad (3) $$

for $x < 0$ and $t \in (0, T)$, with the linear shallow water equations (LSWE)

$$ \partial_t \zeta + \partial_x V = 0 \qquad (4) $$

$$ \partial_t V + \partial_x \zeta = 0 \qquad (5) $$

for $x > 0$ and $t \in (0, T)$ with $T > 0$, to represent the 1D propagation of water waves in shallow water. Here $\partial_t, \partial_x$ denote partial derivatives in the time and space variables $t, x$; $V(t, x)$ and $\zeta(t, x)$ stand for the vertically-averaged velocity and the free-surface level over its state at rest; $\phi(t, x)$ is an auxiliary variable for the elliptic part of the problem; and $\mu > 0$ is the asymptotic parameter characterizing the wave dispersion.

José Galaz
INRIA, Team LEMON, Centre Inria d'Université Côte d'Azur, Antenne Montpellier, Bat 5 CC05
017, 860 rue Saint-Priest, 34095 Montpellier Cedex 5 France, e-mail: jose.galaz@inria.fr

Maria Kazolea
INRIA, Team CARDAMOM, INRIA Bordeaux-Sud-Ouest, 200 Avenue de la Vieille Tour, 33405
Talence cedex France e-mail: maria.kazolea@inria.fr

Antoine Rousseau
INRIA, Team LEMON, Centre Inria d'Université Côte d'Azur, Antenne Montpellier, Bat 5 CC05
017, 860 rue Saint-Priest, 34095 Montpellier Cedex 5 France, e-mail: antoine.rousseau@inria.fr
In the nonlinear case, Boussinesq-type equations, such as the Green-Naghdi equations (GNE), have been coupled with the nonlinear shallow water equations (NSWE) to take advantage of their physical-modeling features: the dispersive terms in the GNE can be used to accurately represent the phase and amplitude of waves in the shoaling zone, while shock-capturing well-balanced finite volume schemes for the NSWE can mimic the energy dissipation of wave-breaking and provide a robust handling of vanishing water depths without ad-hoc parametrizations. However, this coupled model has been shown to be unstable unless dissipative terms are added [1].

Since this coupling is a "divide and conquer" type of problem, domain decomposition methods (DDM) can help to obtain further insights. Usually, coupling conditions have been derived for each equation on a case-by-case basis (e.g., [2, 3, 4, 5]). Here we explore a different approach, based on the steps of the derivation of the GNE and NSWE [6, ch. 1 and 5]. First, a DDM of the linearized Euler equations in the discrete level is defined, based on the Neumann-Dirichlet method. Then, recalling that both GNE and NSWE derive from the Euler equations, we derive transmission conditions for the children asymptotic equations by taking the vertical average of the original operators and truncating the resulting expression according to an asymptotic expansion of the velocity potential. We examine this approach in the homogeneous case first, when coupling LGNE with LGNE, and then in the heterogeneous case, coupling the LGNE with the LSWE.

## 2 A domain decomposition of the free-surface Euler equations

The linear Euler equations for an incompressible fluid and irrotational flow in one horizontal dimension can be formulated as an elliptic problem for the velocity potential $\Phi(t, x, z)$ and two evolution equations for the free-surface function $\zeta(t, x)$ and the trace of the potential at the surface $\psi(t, x)$ respectively (see [6, ch. 1.1.3]). A finite-difference discretization of the equations for $\zeta_i^n = \zeta(t^n, x_i)$ and $\psi_i^n = \psi(t^n, x_i)$, in a uniform grid $x_i = i\Delta x$ and $t^n = n\Delta t$, with $i \in \mathcal{N} = \{1...N_x - 2\}$ and $n = 1, 2, ...$, is given by

$$ \frac{\zeta_i^{n+1} - \zeta_i^n}{\Delta t} + \frac{V_i^n - V_{i-1}^n}{\Delta x} = 0 \quad \text{for } i \in \mathcal{N} \qquad (6) $$

$$ \frac{\psi_i^{n+1} - \psi_i^n}{\Delta t} + \zeta_i^n = 0 \quad \text{for } i \in \mathcal{N} \qquad (7) $$

where $V_i^n$ is the discrete vertically-averaged velocity cf. [6, ch. 3.31] given by

$$ V_i^n = \sum_{j=1}^{N_z-1} \frac{\Phi_{i+1,j}^n - \Phi_{i,j}^n}{\Delta x} \Delta z \qquad (8) $$

Boundary conditions at nodes $i = 0$ and $i = N_x - 1$ can be $V_i^n = 0$ and $\zeta_0^n = \zeta_1^n$, $\zeta_{N_x-1}^n = \zeta_{N_x-2}^n$ [7, eqs. (28), (37)]. In equation (8), $\Phi_{i,j}^n = \Phi(t^n, x_i, z_j)$ is the discrete
In the nonlinear case, Boussinesq-type equations, such as the Green-Naghdi equations (GNE), have been coupled with the nonlinear shallow water equations (NSWE) to take advantage of their physical-modeling features: the dispersive terms in the GNE can be used to accurately represent the phase and amplitude of waves in the shoaling zone, while shock-capturing well-balanced finite volume schemes for the NSWE can mimic the energy dissipation of wave-breaking and provide a robust handling of vanishing water depths without ad-hoc parametrizations. However, this coupled model has been shown to be unstable unless dissipative terms are added [1].

Since this coupling is a "divide and conquer" type of problem, domain decomposition methods (DDM) can help to obtain further insights. Usually, coupling conditions have been derived for each equation on a case-by-case basis (e.g., [2, 3, 4, 5]). Here we explore a different approach, based on the steps of the derivation of the GNE and NSWE [6, ch. 1 and 5]. First, a DDM of the linearized Euler equations in the discrete level is defined, based on the Neumann-Dirichlet method. Then, recalling that both GNE and NSWE derive from the Euler equations, we derive transmission conditions for the children asymptotic equations by taking the vertical average of the original operators and truncating the resulting expression according to an asymptotic expansion of the velocity potential. We examine this approach in the homogeneous case first, when coupling LGNE with LGNE, and then in the heterogeneous case, coupling the LGNE with the LSWE.

## 2 A domain decomposition of the free-surface Euler equations

The linear Euler equations for an incompressible fluid and irrotational flow in one horizontal dimension can be formulated as an elliptic problem for the velocity potential $\Phi(t, x, z)$ and two evolution equations for the free-surface function $\zeta(t, x)$ and the trace of the potential at the surface $\psi(t, x)$ respectively (see [6, ch. 1.1.3]). A finite-difference discretization of the equations for $\zeta_i^n = \zeta(t^n, x_i)$ and $\psi_i^n = \psi(t^n, x_i)$, in a uniform grid $x_i = i\Delta x$ and $t^n = n\Delta t$, with $i \in \mathcal{N} = \{1...N_x - 2\}$ and $n = 1, 2, ...$, is given by

$$ \frac{\zeta_i^{n+1} - \zeta_i^n}{\Delta t} + \frac{V_i^n - V_{i-1}^n}{\Delta x} = 0 \quad \text{for } i \in \mathcal{N} \qquad (6) $$

$$ \frac{\psi_i^{n+1} - \psi_i^n}{\Delta t} + \zeta_i^n = 0 \quad \text{for } i \in \mathcal{N} \qquad (7) $$

where $V_i^n$ is the discrete vertically-averaged velocity cf. [6, ch. 3.31] given by

$$ V_i^n = \sum_{j=1}^{N_z-1} \frac{\Phi_{i+1,j}^n - \Phi_{i,j}^n}{\Delta x} \Delta z \qquad (8) $$

Boundary conditions at nodes $i = 0$ and $i = N_x - 1$ can be $V_i^n = 0$ and $\zeta_0^n = \zeta_1^n$, $\zeta_{N_x-1}^n = \zeta_{N_x-2}^n$ [7, eqs. (28), (37)]. In equation (8), $\Phi_{i,j}^n = \Phi(t^n, x_i, z_j)$ is the discrete
velocity potential computed from

$$ \mu \frac{\Phi_{i+1,j}^n + \Phi_{i-1,j}^n - 2\Phi_{i,j}^n}{\Delta x^2} + \frac{\Phi_{i,j+1}^n + \Phi_{i,j-1}^n - 2\Phi_{i,j}^n}{\Delta z^2} = 0 \quad \text{for } i, j \in \mathcal{N} \times \{1, \dots, N_z - 1\} \quad (9) $$

on a grid $(x_i, z_j)$ such that $z_j = j\Delta z - 1$ and $z_{N_z-1} = 0$ and its boundary conditions are $\Phi_{i,N_z-1}^n = \psi_i^n, \Phi_{i,0} = \Phi_{i,1}, \Phi_{0,j} = \Phi_{1,j}$ and $\Phi_{N_x-1,j} = \Phi_{N_x-2,j}$

We decompose the domain in two components with a vertical interface located at $i = l \in \mathcal{N}$. To do this let $\zeta_{i,s}^{n+1}, \psi_{i,s}^{n+1}, \Phi_{i,j,s}^{n+1}, V_{i,s}^n$ be the values of the unknowns on each subdomain $s = 1, 2$. These variables are computed from equations (6), (7), (8) but with $i$ in $\mathcal{N}_s$ instead of $\mathcal{N}$, with $\mathcal{N}_1 = \{1, \dots, l-1\}$ and $\mathcal{N}_2 = \{l, \dots, N_x - 2\}$; equation (9) is solved with $i \in \mathcal{N}_1$ for $s = 1$ and $i \in \mathcal{N}_2 \setminus \{l\}$ for $s = 2$. To obtain the same solution as the monodomain problem, equations (6) and (7) are complemented with Dirichlet transmission conditions

$$ \begin{aligned} \zeta_{l,1}^{n+1} &= \zeta_{l,2}^{n+1} \quad V_{l,1}^{n+1} = V_{l,2}^{n+1} \\ \zeta_{l-1,2}^{n+1} &= \zeta_{l-1,1}^{n+1} \quad V_{l-1,2}^{n+1} = V_{l-1,1}^{n+1} \end{aligned} \quad (10) $$

while equation (9) uses Neumann and Dirichlet transmission conditions

$$ \begin{aligned} & \frac{\Phi_{l,j,1}^n - \Phi_{l-1,j,1}^n}{\Delta x} - \frac{\Delta x}{2\mu\Delta z^2} (\Phi_{l,j+1,1}^n + \Phi_{l,j-1,1}^n - 2\Phi_{l,j,1}^n) && (11) \\ = & \frac{\Phi_{l+1,j,2}^n - \Phi_{l,j,2}^n}{\Delta x} + \frac{\Delta x}{2\mu\Delta z^2} (\Phi_{l,j+1,2}^n + \Phi_{l,j-1,2}^n - 2\Phi_{l,j,2}^n) \\ & \Phi_{l,j,2}^n = \Phi_{l,j,2}^1 && (12) \end{aligned} $$

which include an $O(\Delta x)$ term necessary in finite-difference schemes to preserve the monodomain solution [8]. This scheme satisfies

* $\Phi_{i,j,*}^n = \Phi_{i,j}^n$ with $\Phi_{i,j,*}^n = \Phi_{i,j,1}^n$ if $i \le l$ and $\Phi_{i,j,2}^n = \Phi_{i,j}^n$ if $i > l$. This means that the solution $\Phi_{i,j,*}^n$ formed by both subdomains will be equal to the monodomain solution $\Phi_{i,j}^n$

* A parallel or alternating method to solve (9) with (12) will be convergent if $L_1 < L_2$, with $L_1 = l\Delta x$ and $L_2 = (N_x - l)\Delta x$ (see [9, eq. (2.5)] with $\theta = 1$] for example).

### 3 Asymptotic domain-decomposition method

In the first part of this section we drop the time superscript $n$ and introduce the asymptotic-degree superscript $(k) = (1), (2), \dots$. We need an asymptotic expansion $\Phi_{i,j} = \Phi_{i,j}^{(0)} + \mu\Phi_{i,j}^{(1)} + \dots + \mu^k\Phi_{i,j}^{(k)}$ of the solution to the discrete Laplace equation (9).
velocity potential computed from

$$ \mu \frac{\Phi_{i+1,j}^n + \Phi_{i-1,j}^n - 2\Phi_{i,j}^n}{\Delta x^2} + \frac{\Phi_{i,j+1}^n + \Phi_{i,j-1}^n - 2\Phi_{i,j}^n}{\Delta z^2} = 0 \quad \text{for } i, j \in \mathcal{N} \times \{1, \dots, N_z - 1\} \quad (9) $$

on a grid $(x_i, z_j)$ such that $z_j = j\Delta z - 1$ and $z_{N_z-1} = 0$ and its boundary conditions are $\Phi_{i,N_z-1}^n = \psi_i^n, \Phi_{i,0} = \Phi_{i,1}, \Phi_{0,j} = \Phi_{1,j}$ and $\Phi_{N_x-1,j} = \Phi_{N_x-2,j}$

We decompose the domain in two components with a vertical interface located at $i = l \in \mathcal{N}$. To do this let $\zeta_{i,s}^{n+1}, \psi_{i,s}^{n+1}, \Phi_{i,j,s}^{n+1}, V_{i,s}^n$ be the values of the unknowns on each subdomain $s = 1, 2$. These variables are computed from equations (6), (7), (8) but with $i$ in $\mathcal{N}_s$ instead of $\mathcal{N}$, with $\mathcal{N}_1 = \{1, \dots, l-1\}$ and $\mathcal{N}_2 = \{l, \dots, N_x - 2\}$; equation (9) is solved with $i \in \mathcal{N}_1$ for $s = 1$ and $i \in \mathcal{N}_2 \setminus \{l\}$ for $s = 2$. To obtain the same solution as the monodomain problem, equations (6) and (7) are complemented with Dirichlet transmission conditions

$$ \begin{aligned} \zeta_{l,1}^{n+1} &= \zeta_{l,2}^{n+1} \quad V_{l,1}^{n+1} = V_{l,2}^{n+1} \\ \zeta_{l-1,2}^{n+1} &= \zeta_{l-1,1}^{n+1} \quad V_{l-1,2}^{n+1} = V_{l-1,1}^{n+1} \end{aligned} \quad (10) $$

while equation (9) uses Neumann and Dirichlet transmission conditions

$$ \begin{aligned} & \frac{\Phi_{l,j,1}^n - \Phi_{l-1,j,1}^n}{\Delta x} - \frac{\Delta x}{2\mu\Delta z^2} (\Phi_{l,j+1,1}^n + \Phi_{l,j-1,1}^n - 2\Phi_{l,j,1}^n) && (11) \\ = & \frac{\Phi_{l+1,j,2}^n - \Phi_{l,j,2}^n}{\Delta x} + \frac{\Delta x}{2\mu\Delta z^2} (\Phi_{l,j+1,2}^n + \Phi_{l,j-1,2}^n - 2\Phi_{l,j,2}^n) \\ & \Phi_{l,j,2}^n = \Phi_{l,j,2}^1 && (12) \end{aligned} $$

which include an $O(\Delta x)$ term necessary in finite-difference schemes to preserve the monodomain solution [8]. This scheme satisfies

* $\Phi_{i,j,*}^n = \Phi_{i,j}^n$ with $\Phi_{i,j,*}^n = \Phi_{i,j,1}^n$ if $i \le l$ and $\Phi_{i,j,1}^n = \Phi_{i,j,2}^n$ if $i > l$. This means that the solution $\Phi_{i,j,*}^n$ formed by both subdomains will be equal to the monodomain solution $\Phi_{i,j}^n$

* A parallel or alternating method to solve (9) with (12) will be convergent if $L_1 < L_2$, with $L_1 = l\Delta x$ and $L_2 = (N_x - l)\Delta x$ (see [9, eq. (2.5)] with $\theta = 1$] for example).

### 3 Asymptotic domain-decomposition method

In the first part of this section we drop the time superscript $n$ and introduce the asymptotic-degree superscript $(k) = (1), (2), \dots$. We need an asymptotic expansion $\Phi_{i,j} = \Phi_{i,j}^{(0)} + \mu\Phi_{i,j}^{(1)} + \dots + \mu^k\Phi_{i,j}^{(k)}$ of the solution to the discrete Laplace equation (9).
At first order $\Phi_{i,j} = \Phi_{i,j}^{(0)} + O(\mu)$ which substituted on equation (9) and discarding $O(\mu)$ terms leads to $\Phi_{i,j+1}^{(0)} + \Phi_{i,j-1}^{(0)} - 2\Phi_{i,j}^{(0)} = 0$, whose solution is $\Phi_{i,j}^{(0)} = \psi_i$ so the first order expansion is

$$ \Phi_{i,j} = \psi_i + O(\mu) \quad (13) $$

Similarly, at second order, replacing $\Phi_{i,j} = \Phi_{i,j}^{(0)} + \mu\Phi_{i,j}^{(1)} + O(\mu^2)$ into (9)

$$ \frac{\Phi_{i,j+1}^{(1)} + \Phi_{i,j-1}^{(1)} - 2\Phi_{i,j}^{(1)}}{\Delta x^2} = - \frac{\psi_{i+1} + \psi_{i-1} - 2\psi_i}{\Delta x^2} \quad (14) $$

whose solution gives us the second order expansion

$$ \Phi_{i,j} = \psi_i - \mu \frac{\Delta z^2}{2\Delta x^2} (\psi_{i+1} + \psi_{i-1} - 2\psi_i) (j^2 - j - (N_z - 1)(N_z - 2)) + O(\mu^2) \quad (15) $$

Substituting (13) into (8) and using that $(N_z - 1)\Delta z = 1$ one obtains that at first order

$$ \frac{\psi_{i+1} - \psi_i}{\Delta x} = V_i + O(\mu) \quad (16) $$

From (15) we can proceed similarly to obtain the second order expansion for $V_i$

$$ V_i = \frac{\psi_{i+1} - \psi_i}{\Delta x} - \nu T \left( \frac{\psi_{i+1} - \psi_i}{\Delta x} \right) + O(\mu^2) \quad (17) $$

where $TV_i = -(V_{i+1} + V_{i-1} - 2V_i)/(3\Delta x^2)$ and $\nu = \mu(1 - \Delta z/2)(1 - \Delta z)$. We can now substitute the first order approximation (16) into (17) and isolate $(\psi_{i+1} - \psi_i)/\Delta x$ to obtain

$$ \frac{\psi_{i+1} - \psi_i}{\Delta x} = V_i + \nu T V_i + O(\mu^2) \quad (18) $$

To substitute (18) into (7) let us apply a forward finite-difference in $x$ to equation (7). Introducing the notation $D_x^+ f_i = (f_{i+1} - f_i)/\Delta x$, $D_z^+ f_j = (f_{j+1} - f_j)/\Delta z$ and $D_z^- f_j = (f_j - f_{j-1})/\Delta z$, and the time superscript $n$, (8) becomes $V_i^n = \sum_{j=1}^{N_z-1} D_x^+ \Phi_{i,j}^n \Delta z$. Using the asymptotic expansion (18), discarding terms of size $O(\mu^2)$ and rearranging one finally obtains the discrete momentum equation of the LGNE

$$ \frac{V_i^{n+1} - V_i^n}{\Delta t} + \frac{\zeta_{i+1}^n - \zeta_i^n}{\Delta x} = \phi_i^n \\ (1 + \nu T)\phi_i^n = \nu T(D_x^+ \zeta_i^n) \quad (19) $$

with $\phi_i^n$ an auxiliary variable for the new elliptic problem. And discarding all $O(\mu)$ terms the discrete momentum equation of the LSWE reads

$$ \frac{V_i^{n+1} - V_i^n}{\Delta t} + \frac{\zeta_{i+1}^n - \zeta_i^n}{\Delta x} = 0 \quad (20) $$
At first order $\Phi_{i,j} = \Phi_{i,j}^{(0)} + O(\mu)$ which substituted on equation (9) and discarding $O(\mu)$ terms leads to $\Phi_{i,j+1}^{(0)} + \Phi_{i,j-1}^{(0)} - 2\Phi_{i,j}^{(0)} = 0$, whose solution is $\Phi_{i,j}^{(0)} = \psi_i$ so the first order expansion is

$$ \Phi_{i,j} = \psi_i + O(\mu) \qquad (13) $$

Similarly, at second order, replacing $\Phi_{i,j} = \Phi_{i,j}^{(0)} + \mu\Phi_{i,j}^{(1)} + O(\mu^2)$ into (9)

$$ \frac{\Phi_{i,j+1}^{(1)} + \Phi_{i,j-1}^{(1)} - 2\Phi_{i,j}^{(1)}}{\Delta x^2} = -\frac{\psi_{i+1} + \psi_{i-1} - 2\psi_i}{\Delta x^2} \qquad (14) $$

whose solution gives us the second order expansion

$$ \Phi_{i,j} = \psi_i - \mu \frac{\Delta z^2}{2\Delta x^2} (\psi_{i+1} + \psi_{i-1} - 2\psi_i) (j^2 - j - (N_z - 1)(N_z - 2)) + O(\mu^2) \quad (15) $$

Substituting (13) into (8) and using that $(N_z - 1)\Delta z = 1$ one obtains that at first order

$$ \frac{\psi_{i+1} - \psi_i}{\Delta x} = V_i + O(\mu) \qquad (16) $$

From (15) we can proceed similarly to obtain the second order expansion for $V_i$

$$ V_i = \frac{\psi_{i+1} - \psi_i}{\Delta x} - \nu T \left( \frac{\psi_{i+1} - \psi_i}{\Delta x} \right) + O(\mu^2) \qquad (17) $$

where $TV_i = -(V_{i+1} + V_{i-1} - 2V_i)/(3\Delta x^2)$ and $\nu = \mu(1 - \Delta z/2)(1 - \Delta z)$. We can now substitute the first order approximation (16) into (17) and isolate $(\psi_{i+1} - \psi_i)/\Delta x$ to obtain

$$ \frac{\psi_{i+1} - \psi_i}{\Delta x} = V_i + \nu T V_i + O(\mu^2) \qquad (18) $$

To substitute (18) into (7) let us apply a forward finite-difference in $x$ to equation (7). Introducing the notation $D_x^+ f_i = (f_{i+1} - f_i)/\Delta x$, $D_z^+ f_j = (f_{j+1} - f_j)/\Delta z$ and $D_z^- f_j = (f_j - f_{j-1})/\Delta z$, and the time superscript $n$, (8) becomes $V_i^n = \sum_{j=1}^{N_z-1} D_x^+ \Phi_{i,j}^n \Delta z$. Using the asymptotic expansion (18), discarding terms of size $O(\mu^2)$ and rearranging one finally obtains the discrete momentum equation of the LGNE

$$ \frac{V_i^{n+1} - V_i^n}{\Delta t} + \frac{\zeta_{i+1}^n - \zeta_i^n}{\Delta x} = \phi_i^n \\ (1 + \nu T)\phi_i^n = \nu T(D_x^+ \zeta_i^n) \qquad (19) $$

with $\phi_i^n$ an auxiliary variable for the new elliptic problem. And discarding all $O(\mu)$ terms the discrete momentum equation of the LSWE reads

$$ \frac{V_i^{n+1} - V_i^n}{\Delta t} + \frac{\zeta_{i+1}^n - \zeta_i^n}{\Delta x} = 0 \qquad (20) $$
We can proceed in a similar fashion to derive asymptotic versions of the Neumann boundary condition. To do this let us multiply equation (12) by $\Delta z$ and sum it up from $j = 1$ to $j = N_z - 2$. Using the formula for $V_i$ and simplifying, one obtains that (12) can be written as

$$ 
\begin{aligned}
& V_{l-1,1}^n - D_x^+ \Phi_{l-1,N_z-1,1}^n \Delta z - \frac{\Delta x}{2\mu} D_z^- \Phi_{l,N_z-1,1}^n \\
& = V_{l,2}^n - D_x^+ \Phi_{l,N_z-1,2}^n \Delta z + \frac{\Delta x}{2\mu} D_z^- \Phi_{l,N_z-1,2}^n
\end{aligned} 
\quad (21) $$

To further simplify the remaining $\Phi$ terms we will use that $\partial_z \Phi_{z=0}/\mu = \partial_x V + O(\epsilon)$ (from Ref. [6, eqs. (1.29) and Proposition 3.35]), and substitute the discrete derivatives to finally obtain $1/\mu D_z^- \Phi_{i,N_z-1} = -D_x^p V_i + O(\epsilon, \Delta z/\mu, \Delta x^p)$, where $D_x^p V_i = \partial_x V|_{x=x_i} + O(\Delta x^p)$ is a finite-difference operator to be defined. If we replace this back into the equation, use equation (13) written as $D_x \Phi_{i,N_z-1} = V_i + O(\mu)$, and replace $D_x^p$ with the backward finite difference $D_x^-$ in the left-hand-side of the equation and a forward finite difference $D_x^+$ in the right-hand-side, we can write the original Neumann boundary condition as

$$ \frac{1}{2}(V_{l,1}^n + V_{l-1,1}^n) - V_{l-1,1}^n \Delta z = \frac{1}{2}(3V_{l,2}^n - V_{l+1,2}^n) - V_{l,2}^n \Delta z \quad (22) $$

Substituting the auxiliary variable for the LGNE, $\phi_{i,s}^n = D_t^+ V_{i,s}^n + D_x^+ \zeta_{i,s}^n$.

$$ \frac{1}{2}(\phi_{l,1}^n + \phi_{l-1,1}^n) - \phi_{l-1,1}^n \Delta z = \frac{1}{2}(3\phi_{l,2}^n - \phi_{l+1,2}^n) - \phi_{l,2}^n \Delta z \\ - \left( \frac{3\Delta x^2}{2} T(D_x^+ \zeta_l^n) - (D_x^+ \zeta_l^n - D_x^+ \zeta_{l-1}^n) \Delta z \right) \quad (23) $$

Summarizing, for the homogeneous case, the domain decomposition of the LGNE reads, for each subdomain $s = 1, 2$,

$$ 
\begin{aligned}
& \frac{\zeta_i^{n+1} - \zeta_i^n}{\Delta t} + \frac{V_i^n - V_{i-1}^n}{\Delta x} = 0 \quad \text{for } i \in \mathcal{N}_s \\
& \frac{V_i^{n+1} - V_i^n}{\Delta t} + \frac{\zeta_{i+1}^n - \zeta_i^n}{\Delta x} = \phi_i \quad \text{for } i \in \mathcal{N}_s \quad (24) \\
& \phi_{i,s}^n + \nu T \phi_{i,s}^n = \nu T (D_x^+ \zeta_i) \quad \text{for } i \in \mathcal{N}_s \\
& \phi_{0,1}^n = \phi_{Left} \\
& \phi_{N_x-1,2}^n = \phi_{Right}
\end{aligned}
$$

and at the interface $i = l$ equation (10) holds for $V_{l,s}^n$ and $\zeta_{l,s}^n$, while for $\phi_{l,s}^n$
We can proceed in a similar fashion to derive asymptotic versions of the Neumann boundary condition. To do this let us multiply equation (12) by $\Delta z$ and sum it up from $j = 1$ to $j = N_z - 2$. Using the formula for $V_i$ and simplifying, one obtains that (12) can be written as

$$ 
\begin{aligned}
V_{l-1,1}^n - D_x^+ \Phi_{l-1,N_z-1,1}^n \Delta z - \frac{\Delta x}{2\mu} D_z^- \Phi_{l,N_z-1,1}^n \\
= V_{l,2}^n - D_x^+ \Phi_{l,N_z-1,2}^n \Delta z + \frac{\Delta x}{2\mu} D_z^- \Phi_{l,N_z-1,2}^n
\end{aligned} 
\quad (21) $$

To further simplify the remaining $\Phi$ terms we will use that $\partial_z \Phi_{z=0}/\mu = \partial_x V + O(\epsilon)$ (from Ref. [6, eqs. (1.29) and Proposition 3.35]), and substitute the discrete derivatives to finally obtain $1/\mu D_z^- \Phi_{i,N_z-1} = -D_x^p V_i + O(\epsilon, \Delta z/\mu, \Delta x^p)$, where $D_x^p V_i = \partial_x V|_{x=x_i} + O(\Delta x^p)$ is a finite-difference operator to be defined. If we replace this back into the equation, use equation (13) written as $D_x \Phi_{i,N_z-1} = V_i + O(\mu)$, and replace $D_x^p$ with the backward finite difference $D_x^-$ in the left-hand-side of the equation and a forward finite difference $D_x^+$ in the right-hand-side, we can write the original Neumann boundary condition as

$$ \frac{1}{2}(V_{l,1}^n + V_{l-1,1}^n) - V_{l-1,1}^n \Delta z = \frac{1}{2}(3V_{l,2}^n - V_{l+1,2}^n) - V_{l,2}^n \Delta z \quad (22) $$

Substituting the auxiliary variable for the LGNE, $\phi_{i,s}^n = D_t^+ V_{i,s}^n + D_x^+ \zeta_{i,s}^n$.

$$ \frac{1}{2}(\phi_{l,1}^n + \phi_{l-1,1}^n) - \phi_{l-1,1}^n \Delta z = \frac{1}{2}(3\phi_{l,2}^n - \phi_{l+1,2}^n) - \phi_{l,2}^n \Delta z \\ - \left( \frac{3\Delta x^2}{2} T(D_x^+ \zeta_l^n) - (D_x^+ \zeta_l^n - D_x^+ \zeta_{l-1}^n) \Delta z \right) \quad (23) $$

Summarizing, for the homogeneous case, the domain decomposition of the LGNE reads, for each subdomain $s = 1, 2$,

$$ 
\begin{aligned}
\frac{\zeta_i^{n+1} - \zeta_i^n}{\Delta t} + \frac{V_i^n - V_{i-1}^n}{\Delta x} &= 0 \quad \text{for } i \in \mathcal{N}_s \\
\frac{V_i^{n+1} - V_i^n}{\Delta t} + \frac{\zeta_{i+1}^n - \zeta_i^n}{\Delta x} &= \phi_i \quad \text{for } i \in \mathcal{N}_s \quad (24) \\
\phi_{i,s}^n + \nu T \phi_{i,s}^n &= \nu T (D_x^+ \zeta_i) \quad \text{for } i \in \mathcal{N}_s \\
\phi_{0,1}^n &= \phi_{Left} \\
\phi_{N_x-1,2}^n &= \phi_{Right}
\end{aligned}
$$

and at the interface $i = l$ equation (10) holds for $V_{l,s}^n$ and $\zeta_{l,s}^n$, while for $\phi_{l,s}^n$
$$ \frac{1}{2}(\phi_{l,1}^n + \phi_{l-1,1}^n) - \phi_{l-1,1}^n \Delta z = \frac{1}{2}(3\phi_{l,2}^n - \phi_{l+1,2}^n) - \phi_{l,2}^n \Delta z \\ - \left(\frac{3\Delta x^2}{2}T(D_x^+ \zeta_l^n) - (D_x^+ \zeta_l^n - D_x^+ \zeta_{l-1}^n)\Delta z\right) \\ \phi_{l,2}^n = \phi_{l,1}^n \qquad (25) $$

To test this asymptotic domain decomposition method with an additive iterative scheme we compare the monodomain solution with the DDM solution for $x \in (0, 1)$, $\phi_{Left} = 0$, $\phi_{Right} = -0.5$, $\partial_x^3 \zeta^n = -1$, $\mu = 3$, $\Delta x = 0.05$, $\Delta z = \Delta x^2$. These parameters are convenient to avoid overflow in intermediate calculations in the formula of the analytical solution. Figure 1 (left) shows the $L^2$ distance between the discrete monodomain solution and each subdomain solution at each iteration for an interface located at $x = 0.4$ where we can see that the DDM diverges.

<figure><img src="image_1.png" /><figcaption>Fig. 1: L2 distance to monodomain of the solution on each subdomain as a function of the iteration number when using (25) and (26) respectively.</figcaption></figure>

To fix this situation we can subtract $\phi_{l,1}^n$ from each side of the first equation of (25), use that $\phi_{l,1}^n = \phi_{l,2}^n$, and multiply the equation by $-\frac{2}{\Delta x}$. The coupling conditions become

$$ \frac{\phi_{l,1}^n - \phi_{l-1,1}^n}{\Delta x} + 2 \frac{\Delta z}{\Delta x} \phi_{l-1,1}^n = \frac{\phi_{l+1,2}^n - \phi_{l,2}^n}{\Delta x} + 2 \frac{\Delta z}{\Delta x} \phi_{l,2}^n \\ + 3 \Delta x T(D_x^+ \zeta_l^n) - 2 \frac{(D_x^+ \zeta_l^n - D_x^+ \zeta_{l-1}^n)}{\Delta x} \Delta z \qquad (26) \\ \phi_{l,2}^n = \phi_{l,1}^n $$

As before, Figure 1 (right) shows the L2 distance to the monodomain solution at each iteration of the DDM. In contrast with the previous case now the algorithm converges. This happens because (26) is also consistent at $O(\Delta x)$ with a Neumann boundary condition. To see this notice that $3\Delta x T(D_x^+ \zeta_l^n) = O(\Delta x)$ and $(D_x^+ \zeta_l^n - D_x^+ \zeta_{l-1}^n)\Delta z/\Delta x = O(\Delta z)$, so if $\Delta z = O(\Delta x^2)$, when taking the limit $\Delta x \to 0$, equation (26) will satisfy $\partial_x \phi_1^n = \partial_x \phi_2^n + O(\Delta x)$ at $x = x_l$, with $\phi_s^n$ the limit of $\phi_{i,s}^n$ when
$$ \frac{1}{2}(\phi_{l,1}^n + \phi_{l-1,1}^n) - \phi_{l-1,1}^n \Delta z = \frac{1}{2}(3\phi_{l,2}^n - \phi_{l+1,2}^n) - \phi_{l,2}^n \Delta z \\ - \left( \frac{3\Delta x^2}{2} T(D_x^+ \zeta_l^n) - (D_x^+ \zeta_l^n - D_x^+ \zeta_{l-1}^n) \Delta z \right) \\ \phi_{l,2}^n = \phi_{l,1}^n \quad (25) $$

To test this asymptotic domain decomposition method with an additive iterative scheme we compare the monodomain solution with the DDM solution for $x \in (0, 1)$, $\phi_{Left} = 0$, $\phi_{Right} = -0.5$, $\partial_x^3 \zeta^n = -1$, $\mu = 3$, $\Delta x = 0.05$, $\Delta z = \Delta x^2$. These parameters are convenient to avoid overflow in intermediate calculations in the formula of the analytical solution. Figure 1 (left) shows the $L^2$ distance between the discrete monodomain solution and each subdomain solution at each iteration for an interface located at $x = 0.4$ where we can see that the DDM diverges.

<figure><img src="image_2.png" /><figcaption>Fig. 1: L2 distance to monodomain of the solution on each subdomain as a function of the iteration number when using (25) and (26) respectively.</figcaption></figure>

To fix this situation we can subtract $\phi_{l,1}^n$ from each side of the first equation of (25), use that $\phi_{l,1}^n = \phi_{l,2}^n$, and multiply the equation by $-\frac{2}{\Delta x}$. The coupling conditions become

$$ \frac{\phi_{l,1}^n - \phi_{l-1,1}^n}{\Delta x} + 2 \frac{\Delta z}{\Delta x} \phi_{l-1,1}^n = \frac{\phi_{l+1,2}^n - \phi_{l,2}^n}{\Delta x} + 2 \frac{\Delta z}{\Delta x} \phi_{l,2}^n \\ + 3 \Delta x T(D_x^+ \zeta_l^n) - 2 \frac{(D_x^+ \zeta_l^n - D_x^+ \zeta_{l-1}^n)}{\Delta x} \Delta z \quad (26) \\ \phi_{l,2}^n = \phi_{l,1}^n $$

As before, Figure 1 (right) shows the L2 distance to the monodomain solution at each iteration of the DDM. In contrast with the previous case now the algorithm converges. This happens because (26) is also consistent at $O(\Delta x)$ with a Neumann boundary condition. To see this notice that $3\Delta x T(D_x^+ \zeta_l^n) = O(\Delta x)$ and $(D_x^+ \zeta_l^n - D_x^+ \zeta_{l-1}^n) \Delta z / \Delta x = O(\Delta z)$, so if $\Delta z = O(\Delta x^2)$, when taking the limit $\Delta x \to 0$, equation (26) will satisfy $\partial_x \phi_1^n = \partial_x \phi_2^n + O(\Delta x)$ at $x = x_l$, with $\phi_s^n$ the limit of $\phi_{i,s}^n$ when
$\Delta x \to 0$. However, solutions of the subdomains are different from the monodomain discrete solution, since the $O(\Delta x)$ terms in equation (26) lead to a linear system that is not equivalent to the monodomain's.

Defining $\phi_i^*$ as $\phi_i^* = \phi_i^1$ if $i \le l$ and $\phi_i^* = \phi_i^2$ if $i > l$, we conclude that in the homogeneous case the asymptotic domain decomposition method:

* Similarly to its parent DDM, it also corresponds to the Neumann-Dirichlet method, so the additive scheme will be convergent when $L_1 < L_2$, with $L_1 = l\Delta x$ and $L_2 = (N_x - l)\Delta x$.

* The monodomain solution $\phi_i$ is different than $\phi_i^*$, the solution formed by each subdomain solution. This is because the asymptotic boundary condition includes $O(\Delta x)$ terms that induce a different linear system than the monodomain. Additional $O(\Delta x)$ terms must be added to fix this situation.

**The heterogeneous case** Now we want to use (26) to couple the LGNE with the LSWE. To do this we can add the constraint $\phi_i^2 = 0$ for $i \ge l+1$ to impose the LSWE on the right side of the domain, and $\phi_i^1$ satisfying the third equation of system (24), to impose the LGNE on the left side of the domain. If we write the Neumann boundary condition (26) as $(\phi_{l,1}^n - \phi_{l-1,1}^n)/\Delta x = (\phi_{l+1,2}^n - \phi_{l,2}^n)/\Delta x + O(\Delta x)$, use the LSWE $\phi_{l+1,2}^n = 0$ and the Dirichlet boundary condition $\phi_{l,2}^n = \phi_{l,1}^n$, we arrive to the condition $\phi_{l,1}^n = \frac{1}{2}\phi_{l-1,1}^n + O(\Delta x^2)$. An interpretation of this formula is that the Neumann-Dirichlet condition becomes a linear interpolation between $\phi_{l-1,1}^n$ and $\phi_{l+1,2}^n = 0$ plus the $O(\Delta x)$ term.

To test this heterogeneous DDM we manufacture the solitary wave of the GNE into the LGNE and LSWE $\zeta(t,x) = \operatorname{sech}^2(x-t+3)$, $V(t,x) = \zeta(t,x)/(1+\varepsilon\zeta(t,x))$ with $\varepsilon = 0.2$, for $(x,t) \in (-11, 11) \times (-3, 3)$, and an interface located at $x=0$. By definition any change on its shape must be due to the influence of the asymptotic transmission conditions. Figure 2 shows the results for $\Delta x = 0.05, 0.02, 0.01$, when half of the solitary wave has crossed the interface at $t=0$ and later at $t=3$. We see that the interface boundary conditions have introduced oscillations and a discontinuity at the interface whose amplitudes grow as $\Delta x$ decreases. This is similar to the results reported by [1].

<figure><img src="image_2.png" /><figcaption>Fig. 2: Comparison of the free surface of a solitary wave calculated with the asymptotic heterogeneous DDM at two time steps and different grids.</figcaption></figure>
$\Delta x \to 0$. However, solutions of the subdomains are different from the monodomain discrete solution, since the $O(\Delta x)$ terms in equation (26) lead to a linear system that is not equivalent to the monodomain's.

Defining $\phi_i^*$ as $\phi_i^* = \phi_i^1$ if $i \le l$ and $\phi_i^* = \phi_i^2$ if $i > l$, we conclude that in the homogeneous case the asymptotic domain decomposition method:

* Similarly to its parent DDM, it also corresponds to the Neumann-Dirichlet method, so the additive scheme will be convergent when $L_1 < L_2$, with $L_1 = l\Delta x$ and $L_2 = (N_x - l)\Delta x$.

* The monodomain solution $\phi_i$ is different than $\phi_i^*$, the solution formed by each subdomain solution. This is because the asymptotic boundary condition includes $O(\Delta x)$ terms that induce a different linear system than the monodomain. Additional $O(\Delta x)$ terms must be added to fix this situation.

**The heterogeneous case** Now we want to use (26) to couple the LGNE with the LSWE. To do this we can add the constraint $\phi_i^2 = 0$ for $i \ge l+1$ to impose the LSWE on the right side of the domain, and $\phi_i^1$ satisfying the third equation of system (24), to impose the LGNE on the left side of the domain. If we write the Neumann boundary condition (26) as $(\phi_{l,1}^n - \phi_{l-1,1}^n)/\Delta x = (\phi_{l+1,2}^n - \phi_{l,2}^n)/\Delta x + O(\Delta x)$, use the LSWE $\phi_{l+1,2}^n = 0$ and the Dirichlet boundary condition $\phi_{l,2}^n = \phi_{l,1}^n$, we arrive to the condition $\phi_{l,1}^n = \frac{1}{2}\phi_{l-1,1}^n + O(\Delta x^2)$. An interpretation of this formula is that the Neumann-Dirichlet condition becomes a linear interpolation between $\phi_{l-1,1}^n$ and $\phi_{l+1,2}^n = 0$ plus the $O(\Delta x)$ term.

To test this heterogeneous DDM we manufacture the solitary wave of the GNE into the LGNE and LSWE $\zeta(t,x) = \operatorname{sech}^2(x-t+3)$, $V(t,x) = \zeta(t,x)/(1+\varepsilon\zeta(t,x))$ with $\varepsilon = 0.2$, for $(x,t) \in (-11, 11) \times (-3, 3)$, and an interface located at $x=0$. By definition any change on its shape must be due to the influence of the asymptotic transmission conditions. Figure 2 shows the results for $\Delta x = 0.05, 0.02, 0.01$, when half of the solitary wave has crossed the interface at $t=0$ and later at $t=3$. We see that the interface boundary conditions have introduced oscillations and a discontinuity at the interface whose amplitudes grow as $\Delta x$ decreases. This is similar to the results reported by [1].

<figure><img src="image_3.png" /><figcaption>Fig. 2: Comparison of the free surface of a solitary wave calculated with the asymptotic heterogeneous DDM at two time steps and different grids.</figcaption></figure>
## 4 Conclusions

Transmission conditions have been derived for the LGNE and LSWE by taking a vertical average and truncating an asymptotic expansion of the transmission conditions of a DDM of the discrete linear free-surface Euler equations. This DDM uses the Neumann-Dirichlet method on the elliptic problem for the velocity potential. In the homogeneous case, when the LGNE are solved on both sides of the interface, we recover the Neumann-Dirichlet method of the elliptic part of the LGNE, plus an $O(\Delta x)$ term. The method has the same convergence property as its parent method but the $O(\Delta x)$ terms make the limit of the subdomain iterations different from the monodomain solution, even though this was imposed on the parent DDM. Also, using more than 2 subdomains could be handled with a relaxation parameter as in [9, section 4.]. In the heterogeneous case the Neumann-Dirichlet method corresponds to a linear interpolation of the elliptic variable between its last value in the LGNE domain and 0, the condition that defines the LSWE, plus an $O(\Delta x^2)$ term. Numerical results show that this induces an unstable scheme, due to oscillations and discontinuities in the interface that grow in amplitude as $\Delta x$ decreases. The next steps could be the introduction of a free parameter in the boundary conditions to optimize the convergence of the method, for example through a Robin boundary condition, and the analysis of this approach in the continuous case.

## References

1. Maria Kazolea and Mario Ricchiuto. On wave breaking for boussinesq-type models. *Ocean Modelling*, 123:16-39, 2018.

2. Maria Kazakova and Pascal Noble. Discrete transparent boundary conditions for the linearized green-naghdi system of equations. *SIAM Journal on Numerical Analysis*, 58(1):657-683, 2020.

3. Christophe Besse, Jean-François Coulombel, and Pascal Noble. Discrete transparent boundary conditions for the two-dimensional leap-frog scheme: approximation and fast implementation. *ESAIM: Mathematical Modelling and Numerical Analysis*, 55:S535-S571, 2021.

4. Chunxiong Zheng, Xin Wen, and Houde Han. Numerical solution to a linearized kdv equation on unbounded domain. *Numerical Methods for Partial Differential Equations: An International Journal*, 24(2):383-399, 2008.

5. Joao Guilherme Caldas Steinstraesser, Rodrigo Alberto Cienfuegos Carrasco, Jose Daniel Galaz Mora, and Antoine Rousseau. A schwarz-based domain decomposition method for the dispersion equation. *Journal of Applied Analysis & Computations*, 8(3):859, 2018.

6. David Lannes. *The water waves problem: mathematical analysis and asymptotics*, volume 188. American Mathematical Soc., 2013.

7. Ge Wei and James T Kirby. Time-dependent numerical code for extended boussinesq equations. *Journal of waterway, port, coastal, and ocean engineering*, 121(5):251-261, 1995.

8. Martin J Gander, Laurence Halpern, and Frédéric Nataf. Optimal schwarz waveform relaxation for the one dimensional wave equation. *SIAM Journal on Numerical Analysis*, 41(5):1643-1681, 2003.

9. Daniele Funaro, Alfio Quarteroni, and Paola Zanolli. An iterative procedure with interface relaxation for domain decomposition methods. *SIAM Journal on Numerical Analysis*, 25(6):1213-1236, 1988.
## 4 Conclusions

Transmission conditions have been derived for the LGNE and LSWE by taking a vertical average and truncating an asymptotic expansion of the transmission conditions of a DDM of the discrete linear free-surface Euler equations. This DDM uses the Neumann-Dirichlet method on the elliptic problem for the velocity potential. In the homogeneous case, when the LGNE are solved on both sides of the interface, we recover the Neumann-Dirichlet method of the elliptic part of the LGNE, plus an $O(\Delta x)$ term. The method has the same convergence property as its parent method but the $O(\Delta x)$ terms make the limit of the subdomain iterations different from the monodomain solution, even though this was imposed on the parent DDM. Also, using more than 2 subdomains could be handled with a relaxation parameter as in [9, section 4.]. In the heterogeneous case the Neumann-Dirichlet method corresponds to a linear interpolation of the elliptic variable between its last value in the LGNE domain and 0, the condition that defines the LSWE, plus an $O(\Delta x^2)$ term. Numerical results show that this induces an unstable scheme, due to oscillations and discontinuities in the interface that grow in amplitude as $\Delta x$ decreases. The next steps could be the introduction of a free parameter in the boundary conditions to optimize the convergence of the method, for example through a Robin boundary condition, and the analysis of this approach in the continuous case.

## References

1. Maria Kazolea and Mario Ricchiuto. On wave breaking for boussinesq-type models. *Ocean Modelling*, 123:16-39, 2018.

2. Maria Kazakova and Pascal Noble. Discrete transparent boundary conditions for the linearized green-naghdi system of equations. *SIAM Journal on Numerical Analysis*, 58(1):657-683, 2020.

3. Christophe Besse, Jean-François Coulombel, and Pascal Noble. Discrete transparent boundary conditions for the two-dimensional leap-frog scheme: approximation and fast implementation. *ESAIM: Mathematical Modelling and Numerical Analysis*, 55:S535-S571, 2021.

4. Chunxiong Zheng, Xin Wen, and Houde Han. Numerical solution to a linearized kdv equation on unbounded domain. *Numerical Methods for Partial Differential Equations: An International Journal*, 24(2):383-399, 2008.

5. Joao Guilherme Caldas Steinstraesser, Rodrigo Alberto Cienfuegos Carrasco, Jose Daniel Galaz Mora, and Antoine Rousseau. A schwarz-based domain decomposition method for the dispersion equation. *Journal of Applied Analysis & Computations*, 8(3):859, 2018.

6. David Lannes. *The water waves problem: mathematical analysis and asymptotics*, volume 188. American Mathematical Soc., 2013.

7. Ge Wei and James T Kirby. Time-dependent numerical code for extended boussinesq equations. *Journal of waterway, port, coastal, and ocean engineering*, 121(5):251-261, 1995.

8. Martin J Gander, Laurence Halpern, and Frédéric Nataf. Optimal schwarz waveform relaxation for the one dimensional wave equation. *SIAM Journal on Numerical Analysis*, 41(5):1643-1681, 2003.

9. Daniele Funaro, Alfio Quarteroni, and Paola Zanolli. An iterative procedure with interface relaxation for domain decomposition methods. *SIAM Journal on Numerical Analysis*, 25(6):1213-1236, 1988.