# Design and analysis of a Schwarz coupling method for 3D Navier-Stokes equations and 2D Shallow Water equations

Manel Tayachi, Céline Acary-Robert, Eric Blayo

► To cite this version:

Manel Tayachi, Céline Acary-Robert, Eric Blayo. Design and analysis of a Schwarz coupling method for 3D Navier-Stokes equations and 2D Shallow Water equations. 2024. hal-04792817

HAL Id: hal-04792817

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

Preprint submitted on 20 Nov 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
Design and analysis of a Schwarz coupling method for 3D
Navier-Stokes equations and 2D Shallow Water equations

Manel Tayachi*, Céline Acary-Robert†and Eric Blayo‡

Abstract

We propose in the present work an iterative coupling method for a dimensionally heterogeneous problem. We consider the 3D linearized hydrostatic Navier-Stokes equations coupled with corresponding 2D linearized shallow water equations. We first show briefly how to derive the 2D linearized shallow water system from the 3D linearized hydrostatic Navier-Stokes system. Then we propose and study a Schwarz-like algorithm to couple the two systems and we prove under some assumptions that the convergence of this Schwarz algorithm is equivalent to the convergence of the classical domain decomposition algorithm of shallow water equations. Finally, we give some theoretical results related to the control of the difference between a global 3D reference solution and the 3D part of the coupled solution. These results are illustrated numerically.

***Keywords***— dimensionally heterogeneous coupling; domain decomposition; multiscale analysis; hydrostatic Navier-Stokes equations, shallow water equations

# 1 Introduction

Modeling complex phenomena, such as some hydrodynamical ones, may require the use of several mathematical and numerical models rather than using a single system of equations. One can for example replace locally the most general (often complex) model, such as the 3D Navier-Stokes system, with simpler models when physics allows it. The simpler models could in particular work in lower dimensions than the dimension of the most general model. For example, the 2D shallow water equations can locally replace the 3D Navier-Stokes system, as they are derived from this set of equations by vertical integration under the assumption of a small domain aspect ratio. One has thus to deal with a dimensionally heterogeneous coupling problem. Such a coupling between a *m*D model and a *n*D model where *n* > *m* may limit heavy computations and lead to efficient results. Such coupled problems were studied in many works and used in several situations. For example, to model blood flows in compliant vessels, Formaggia et al. proposed in [8] to couple the 1D and the 3D Navier-Stokes equations. Miglio et al. [17], Marin and Monnier [15], Finaud-Guyot et al. [7] and Malleron et al. [14] have coupled the 1D and 2D shallow water equations in the context of river flows. Leiva et al. [13] also used dimensionally heterogeneous models in the context of fluid mechanics. In [22], we investigated such a dimensionally heterogeneous coupling problem, and proposed and analysed an efficient Schwarz-like algorithm to couple a 2D Laplace equation with non-symmetric boundary conditions with a corresponding 1D Laplace equation obtained by vertical integration of the 2D model. Similar test cases were addressed by Blanco et al. in [3] and by Panasenko [18] but with different coupling methodologies (variational approach in [3] and asymptotic partial decomposition of domain in [18]). The main difference between the approaches presented in [3] and [18] and those of [22] is the choice of the coupling

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, Grenoble INP, LJK, 38000 Grenoble, France, manel.tayachi@univ-grenoble-alpes.fr</footnote>

<footnote n="2">²Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France, celine.acary-robert@univ-grenoble-alpes.fr</footnote>

<footnote n="3">³Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France, eric.blayo@univ-grenoble-alpes.fr</footnote>
method. In [22] the use of a Schwarz-like method was motivated by the need of a non-intrusive computational approach. Indeed the Schwarz algorithms, developed initially in the context of domain decomposition, do not require changes in the models, but only some exchange of information through boundary conditions. Different strategies can be used for space discretization. Besides, the choice of different time steps is possible by using the so-called Schwarz global-in-time or waveform relaxation algorithms [10, 9, 11]. However, due to their iterative nature, these algorithms may generate huge computations. It is therefore useful to optimize the interface conditions in order to accelerate the convergence. From a theoretical point of view, the so-called perfectly transparent or perfectly absorbing boundary conditions allow an exact convergence in only two iterations. Nevertheless, they are generally non-local in space and/or time, which represents a major obstacle to their direct use, and leads to the quest for local approximations. In the context of dimensionally heterogeneous coupled problems, the situation is further complicated by the addition of extension and reduction operators.

The objective of the present work is to design and analyse an efficient Schwarz-like algorithm to couple the 2D linearized shallow water equations with the 3D linearized hydrostatic Navier-Stokes equations in the context of river flow modeling. Indeed one can consider 2D (or even 1D) shallow water equations in many parts of rivers, and 3D Navier-Stokes equations where 3D effects cannot be neglected. In this regard, this work is in the continuity of [22] and [5]. These two works were the first steps to study and design an efficient Schwarz algorithm for coupling 1D or 2D shallow water models with 3D Navier-Stokes models. In [22], the study of a Schwarz algorithm to couple a 2D Laplace equation with a 1D Laplace equation allowed to understand the main questions that we have to face when dealing with the coupling of dimensionally heterogeneous models and when using a Schwarz-like algorithm in this context. In [21], we proved that under some assumptions the convergence of the coupling algorithm developed in the present work is equivalent to the convergence of the Schwarz waveform relaxation algorithm for linear viscous shallow water equations designed and studied in [5].

This paper is organized as follows: in Section 2 we introduce the 3D linearized hydrostatic Navier-Stokes equations and we briefly show how to derive the 2D linearized shallow water equations. We then define the coupling problem. In Section 3 we define the Schwarz waveform relaxation algorithm and we prove that the convergence of this algorithm is equivalent to the convergence of the classic domain decomposition algorithm for the shallow water system if the bottom friction is neglected. Then we briefly study the well-posedness of the algorithm with Robin-like interface conditions. In Section 4, we prove a theoretical result related to the control of the difference between the 3D part of the coupled solution and the global 3D solution. Finally, Section 5 presents numerical illustrations of the convergence of the classic domain decomposition method for the linearized shallow water system with viscosity, and of the convergence of the coupling algorithm. It also illustrates the influence of the interface position on the coupled solution.

# 2 Coupled models

In order to derive an efficient coupling algorithm between the 3-D hydrostatic Navier-Stokes equations and the 2-D shallow water system, we first write the linearized approximation of the 3-D hydrostatic Navier-Stokes equations and derive the corresponding 2-D linearized shallow water system.

## 2.1 Hydrostatic Navier-Stokes equations

We consider the hydrostatic Navier-Stokes equations:

$$ 
\left\{ 
\begin{array}{l}
\partial_t \mathbf{U}_h + \mathbf{U} \cdot \nabla_h \mathbf{U}_h - \mu \Delta \mathbf{U}_h + \frac{1}{\rho} \nabla_h p = 0 \\
\mathrm{div}_h \mathbf{U}_h + \partial_z W = 0 \\
\partial_z p = -\rho g
\end{array}
\right.
\qquad (1)
$$

in the domain $(x, y, z, t) \in \Omega_t \times (0, T) = \omega \times [-H, \zeta(x, y, t)] \times (0, T)$, where $\omega$ is an open domain of $\mathbb{R}^2$, $T$ denotes the length of the considered time period $(0 < T \le +\infty)$, $\zeta(x, y, t)$ denotes the free surface height, and the depth
$H$ is constant (flat bottom). The unknowns are the velocity $\mathbf{U} = (\mathbf{U}_h, W) = (U, V, W)$ and the pressure $p$. The density $\rho$ and the kinematic viscosity $\mu$ are assumed constant.

This set of equations is derived from the incompressible Navier-Stokes equations by introducing the small aspect ratio $\epsilon = \frac{H}{L_c}$, where $L_c$ is the horizontal characteristic length, and considering the hydrostatic approximation which consists in making $\epsilon = 0$ in the equation for $W$ in the nondimensional system.

These equations are supplemented with initial and boundary conditions. At the bottom of the domain $(\Gamma_B = \{z = -H\})$, we impose a non-penetration condition and a frictionless condition:

$$W(x, y, -H, t) = 0 \qquad (2)$$

and

$$\mu \partial_z \mathbf{U}_{h|z=-H} = 0 \qquad (3)$$

At the free surface $\Gamma_T = \{z = \zeta(x, y, t)\}$, we impose a kinematic boundary condition and the balance of the stresses:

$$\partial_t \zeta + \mathbf{U}_{h|z=\zeta} \cdot \nabla_h \zeta - W_{|z=\zeta} = 0 \qquad (4)$$

and

$$\sigma \cdot \mathbf{n} = 0 \qquad (5)$$

where $\mathbf{n}$ is the outward normal vector to the free surface and $\sigma = -pI + \mu (\nabla \mathbf{U} + \nabla \mathbf{U}^t)$ is the constraint tensor. We neglect here the atmospheric pressure.

At the lateral boundaries, we impose Dirichlet boundary conditions:

$$\mathbf{U}_h = \mathbf{U}^d \quad \text{on} \quad \partial \Omega_t \setminus (\Gamma_B \cup \Gamma_T) \times (0, T) \quad \text{and} \quad \zeta = \zeta^d \quad \text{on} \quad \partial \omega \times (0, T)$$

and in the case of an unbounded domain $\omega$, we impose homogeneous Dirichlet boundary conditions when $\|(x, y)\| \to \infty$.

Finally, initial conditions are provided:

$$\mathbf{U}_h(\cdot, 0) = \mathbf{U}_h^{\text{ini}} \quad \text{in} \quad \Omega_0 \quad \text{and} \quad \zeta(\cdot, 0) = \zeta^{\text{ini}} \quad \text{in} \quad \omega$$

**Remark 1.** In the case of the ocean, the set of hydrostatic equations is supplemented with equations for temperature, salinity and density, and a Coriolis term is also added, leading to the so-called “primitive equations”.

## 2.2 Linearized hydrostatic Navier-Stokes equations

The hydrostatic Navier-Stokes system can be transformed into an equivalent form, as in [2] or in [6]. Using the continuity equation and the non-penetration condition, the vertical velocity $W$ can be written as a function of the horizontal velocity $\mathbf{U}_h$ and, as shown in [6], one can prove by using the Leibniz rule that the kinematic condition (4) is equivalent to a free surface equation. One can then write the pressure $p$ as a function of the free surface $\zeta$ using the hydrostatic condition, and the set of equations becomes:

$$\left\{
\begin{array}{ll}
\partial_t \mathbf{U}_h + \mathbf{U} \cdot \nabla \mathbf{U}_h - \mu \Delta \mathbf{U}_h + g \nabla_h \zeta = 0 & \text{in } \Omega_t \times (0, T) \\
\text{div}_h \mathbf{U}_h + \partial_z W = 0 & \text{in } \Omega_t \times (0, T) \\
\partial_t \zeta + \text{div}_h \left( \int_{-H}^{\zeta} \mathbf{U}_h dz \right) = 0 & \text{in } \omega \times (0, T)
\end{array}
\right.
\qquad (6)$$
In order to set up a Schwarz algorithm and to prove its convergence, we now linearize the problem around a constant velocity $(\mathbf{U}_0, w_0) = (u_0, v_0, 0)$ and the reference value $\zeta_0 = 0$ of the free surface. The linearized system reads (see [2] for more details):

$$
\left\{
\begin{array}{ll}
\partial_t \mathbf{U}_h^G + \mathbf{U}_0 \cdot \nabla_h \mathbf{U}_h^G - \mu \Delta \mathbf{U}_h^G + g \nabla_h \zeta = 0 & \text{in } \Omega \times (0, T) = \omega \times [-H, 0] \times (0, T) \\
\partial_t \zeta + H \operatorname{div}_h (\bar{\mathbf{U}}_h^G) + \mathbf{U}_0 \cdot \nabla_h \zeta = 0 & \text{in } \omega \times (0, T) \\
\mu \partial_z \mathbf{U}_h^G = 0 & \text{at } z=0 \\
\mu \partial_z \mathbf{U}_h^G = 0 & \text{at } z=-H \\
\mathbf{U}_h^G = \mathbf{U}^d & \text{on } \partial \Omega \setminus (\Gamma_B \cup \Gamma_T) \times (0, T) \\
\zeta = \zeta^d & \text{on } \partial \omega \times \mathbb{R}^+ \times (0, T) \\
\mathbf{U}_h^G(., 0) = \mathbf{U}_h^{\text{ini}} & \text{in } \Omega, \text{ and } \zeta(., 0) = \zeta^{\text{ini}} & \text{in } \omega
\end{array}
\right.
\tag{7}
$$

where the superscript G refers to what we will call the *global reference solution*, i.e. the solution of the linearized hydrostatic Navier-Stokes equations throughout the whole domain $\Omega$. The overbar denotes the averaging operator in the vertical direction, defined by:

$$
\bar{f} = \frac{1}{H} \int_{-H}^{0} f(z) dz
$$

The vertical velocity $W^G$ can be obtained from:

$$
\left\{
\begin{array}{ll}
\operatorname{div}_h \mathbf{U}_h^G + \partial_z W^G = 0 & \text{in } \Omega \times (0, T) \\
W^G(x, y, -H, t) = 0 &
\end{array}
\right.
\qquad (8)
$$

In the sequel, we assume that the domain $\Omega$ is divided into two subdomains: a subdomain $\Omega_{2D}$ where the 3-D effects can be neglected, and a subdomain $\Omega_{3D} = \Omega \setminus \Omega_{2D}$ where they cannot. As indicated in [22], the definition of $\Omega_{2D}$ depends on several features such as the domain aspect ratio or forcing terms. We assume that there exists $L_1$ such that:

$$
\Omega_{2D} = \Omega \cap \{x < L_1\} \quad \text{and} \quad \Omega_{3D} = \Omega \cap \{x > L_1\}
$$

Thus, the interface $\Gamma$ between the two subdomains $\Omega_{2D}$ and $\Omega_{3D}$ is located at $x = L_1$. For the sake of simplicity, we suppose that $L_1$ is independent of time. We will now replace the 3D linearized hydrostatic Navier-Stokes equations in $\Omega_{2D}$ by 2D linearized shallow water equations, and study the resulting coupled 2D/3D problem.

## 2.3 Linearized Shallow Water equations

The linearized shallow water system is obtained by averaging system (7). Due to the frictionless condition at the bottom, we do not need to make any other approximation as in [12]. With $U_c$ a characteristic scale for the horizontal velocity, we can introduce the following dimensionless variables and numbers:

$$
(x, y) = L_c(\tilde{x}, \tilde{y}), \quad z = H\tilde{z}, \quad t = \frac{L_c}{U_c}\tilde{t}, \quad \mathbf{U}_h = U_c\tilde{\mathbf{U}}_h, \quad \zeta = H\tilde{\zeta}
$$

and

$$
\nu = \frac{1}{Re} = \frac{\mu}{L_c U_c}, \quad Fr = \frac{U_c}{\sqrt{gH_c}}
$$

Therefore we have the following result:
**Lemme 1.** The linearized shallow water system

$$ egin{cases} 
\partial_t \zeta + H \operatorname{div}_h(\mathbf{u}) + \mathbf{U}_0 \cdot \nabla_h \zeta = 0 \\
\partial_t \mathbf{u} + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{u} - \mu \Delta_h \mathbf{u} + g \nabla_h \zeta = 0
\end{cases} 
\qquad (9) $$

results from an approximation in $O(\varepsilon^2)$ of the linearized hydrostatic Navier-Stokes equations in $\Omega_{2D}$ (where $\varepsilon = \frac{H}{L_c}$ is the previously defined aspect ratio and assuming $L_c = L_1$).

Moreover $\mathbf{U}_h^G(x, y, z, t) = \mathbf{u}(x, y, t) + O(\varepsilon^2)$ for all $(x, y, z, t) \in \Omega_{2D} \times (0, T)$

**Proof** The proof being quite similar to the one given in [12], we only give here the outline and refer to [12] for more details. First, by rescaling the system (7), we obtain (we omit the ~ for the sake of clarity):

$$ egin{cases} 
\partial_t \mathbf{U}_h^G + \mathbf{U}_0 \cdot \nabla_h \mathbf{U}_h^G - \nu \Delta_h \mathbf{U}_h^G - \frac{\nu}{\varepsilon^2} \partial_z^2 \mathbf{U}_h^G + \frac{1}{Fr^2} \nabla_h \zeta = 0 \quad \text{in } \omega_{2D} \times [-1, 0] \times (0, T) \\
\partial_t \zeta + \operatorname{div}_h(\bar{\mathbf{U}}_h) + \mathbf{U}_0 \cdot \nabla_h \zeta = 0 \quad \text{in } \omega_{2D} \times (0, T) \\
\nu \partial_z \mathbf{U}_h^G = 0 \quad \text{at } z=0 \\
\nu \partial_z \mathbf{U}_h^G = 0 \quad \text{at } z=-1 \\
\mathbf{U}_h^G(\cdot, 0) = \mathbf{U}^{\text{ini}} \quad \text{in } \omega_{2D} \times [-1, 0] \text{, and } \zeta(\cdot, 0) = \zeta^{\text{ini}} \quad \text{in } \omega_{2D}
\end{cases} 
\qquad (10) $$

Let suppose as in [12] that $\nu = \varepsilon\nu_0$. Thus, we deduce from the momentum equation and from the boundary condition at the bottom that $\partial_z \mathbf{U}_h^G = O(\varepsilon)$. Then integrating between $-1$ and $z$, the following first order relation holds:

$$ \mathbf{U}_h^G(x, y, z, t) = \mathbf{U}_h^G(x, y, -1, t) + O(\varepsilon) \qquad (11) $$

$$ = \bar{\mathbf{U}}_h^G(x, y, t) + O(\varepsilon) \qquad (12) $$

Now, averaging the momentum equation in (10) between $-1$ and $0$ yields:

$$ \partial_t \bar{\mathbf{U}}_h^G + (\mathbf{U}_0 \cdot \nabla_h) \bar{\mathbf{U}}_h^G + \frac{1}{Fr^2} \nabla_h \zeta = \nu \Delta_h \bar{\mathbf{U}}_h^G \qquad (13) $$

By keeping just the first order terms on $\varepsilon$ in (13), we obtain:

$$ \partial_t \bar{\mathbf{U}}_h^G + (\mathbf{U}_0 \cdot \nabla_h) \bar{\mathbf{U}}_h^G + \frac{1}{Fr^2} \nabla_h \zeta = O(\varepsilon) \qquad (14) $$

The first equations of (10), (11) and (14) give:

$$ egin{aligned} 
\frac{\nu_0}{\varepsilon} \partial_z^2 \mathbf{U}_h^G &= \partial_t \mathbf{U}_h^G + \mathbf{U}_0 \cdot \nabla \mathbf{U}_h^G - \nu \Delta_h \mathbf{U}_h^G + \frac{1}{Fr^2} \nabla_h \zeta \\ 
&= \partial_t \bar{\mathbf{U}}_h^G + \mathbf{U}_0 \cdot \nabla \bar{\mathbf{U}}_h^G + \frac{1}{Fr^2} \nabla_h \zeta + O(\varepsilon) \\ 
&= O(\varepsilon) 
\end{aligned} $$

By integrating between $z$ and $0$ we obtain:

$$ \partial_z \mathbf{U}_h^G = O(\varepsilon^2) $$

Then, by integrating this time between $-1$ and $z$, we have:

$$ \mathbf{U}_h^G = \mathbf{U}_h^G|_{z=-1} + O(\varepsilon^2) \qquad (15) $$

Thus we deduce the second order relation:

$$ \bar{\mathbf{U}}_h^G = \mathbf{U}_h^G|_{z=-1} + O(\varepsilon^2) \qquad (16) $$

Hence the final result, after a back transformation into dimensional variables. $\square$
## 2.4 The coupling problem

Once a criterion has been defined to quantify 3D effects, we have seen that $\Omega_{2D}$ and $\Omega_{3D}$ are supposed to be separated by an interface $x = L_1$. However, in practice, we do not know its exact location. On the other hand, we generally know a value $L_2$ such that we are sure that the 3D effects are significant for $x > L_2$. So, in practice, we will choose a value $L_0 < L_2$ and place our interface at $x = L_0$ (but with no full guarantee that $L_0$ will indeed be smaller than the optimal value $L_1$).

For the sake of simplicity we suppose here that $L_0 = 0$. Thus we consider the linearized shallow water system in $\Omega^- = \omega^- \times [-H, 0]$, where $\omega^- = \omega \cap \{x < 0\}$, and we keep the 3D linearized hydrostatic Navier-Stokes equations in $\Omega^+ = \omega^+ \times [-H, 0]$, where $\omega^+ = \omega \cap \{x > 0\}$. We denote by $\Gamma$ the common 2D interface between these two non-overlapping subdomains $\Omega^-$ and $\Omega^+$: $\Gamma = \{x = 0\}$. We also denote by $\gamma$ the 1D interface between $\omega^-$ and $\omega^+$.

We then consider the two following systems:

$$ \partial_t \mathbf{u} + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{u} + g \nabla_h \zeta - \mu \Delta_h \mathbf{u} = 0 \quad \text{in } \omega^- \times (0, T) \quad (17a) $$

$$ \partial_t \zeta + H \operatorname{div}_h(\mathbf{u}) + \mathbf{U}_0 \cdot \nabla_h \zeta = 0 \quad \text{in } \omega^- \times (0, T) \quad (17b) $$

$$ (\mathbf{u}, \zeta) = (\bar{\mathbf{U}}_{-}^{d}, \zeta_{-}^{d}) \quad \text{on} \quad \partial\omega_{out}^{-} \times (0, T) \quad (17c) $$

$$ \mathbf{u}(\cdot, 0) = \bar{\mathbf{U}}_{-}^{\text{ini}} \quad \text{in } \omega^{-} \quad \text{and} \quad \zeta_{-}(\cdot, 0) = \zeta_{-}^{\text{ini}} \quad \text{in } \omega^{-} \quad (17d) $$

and

$$ \partial_t \mathbf{U}_h + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{U}_h - \mu \Delta \mathbf{U}_h + g \nabla_h \zeta_+ = 0 \quad \text{in } \Omega^+ \times (0, T) \quad (18a) $$

$$ \partial_t \zeta_+ + H \operatorname{div}_h(\bar{\mathbf{U}}_h) + \mathbf{U}_0 \cdot \nabla_h \zeta_+ = 0 \quad \text{in } \omega^+ \times (0, T) \quad (18b) $$

$$ \mu \partial_z \mathbf{U}_h = 0 \quad \text{at} \quad z = 0 \quad (18c) $$

$$ \mu \partial_z \mathbf{U}_h = 0 \quad \text{at} \quad z = -H \quad (18d) $$

$$ \mathbf{U}_h = \mathbf{U}_+^d \quad \text{in} \quad (\partial \Omega_{out}^+ \setminus \Gamma_B \cup \Gamma_T) \times (0, T) \quad (18e) $$

$$ \zeta = \zeta_+^d \quad \text{in} \quad \partial \omega_{out}^+ \times (0, T) \quad (18f) $$

$$ \mathbf{U}_h(\cdot, 0) = \mathbf{U}_+^{\text{ini}} \quad \text{in } \omega^+ \quad \text{and} \quad \zeta(\cdot, 0) = \zeta_2^{\text{ini}} \quad \text{in } \omega^+ \quad (18g) $$

where the boundaries $\partial\omega_{out}^-$ and $\partial\Omega_{out}^+$ are respectively $\partial\omega^- \setminus \gamma$ and $\partial\Omega^+ \setminus \Gamma$. These two systems are to be coupled through the interfaces $\gamma$ and $\Gamma = \gamma \times [-H, 0]$ by defining interface conditions under the form:

$$ \mathcal{B}_- (\mathbf{u}, \zeta_-) = \mathcal{B}_- (\mathcal{R}(\mathbf{U}_h, \zeta_+)) \quad \text{on} \quad \gamma \times (0, T) \quad (19) $$

$$ \mathcal{B}_+ (\mathbf{U}_h, \zeta_+) = \mathcal{B}_+ (\mathcal{E}(\mathbf{u}, \zeta_-)) \quad \text{on} \quad \Gamma \times (0, T) \quad (20) $$

Operators $\mathcal{R}$ and $\mathcal{E}$ are respectively a restriction (or projection) on $\gamma$ and an extension on $\Gamma$ to be defined, and the boundary operators $\mathcal{B}_-$ and $\mathcal{B}_+$ are to be determined.

In order to set up an efficient Schwarz coupling algorithm in the next section, we first need to define the coupling
notion by itself, that is, defining the quantities or values to be exchanged between the two models through the coupling interfaces. Let us first define the natural transmission conditions of system (7):

**Lemme 2.** The transmission condition of system (7) through the interface $\Gamma$ is the continuity of:

$$ \left( \mu \partial_x \mathbf{U}_h^G - u_0 \mathbf{U}_h^G - g \begin{pmatrix} \zeta \\ 0 \end{pmatrix}, u_0 \zeta + h \bar{\mathbf{U}}_h^G \right) \qquad (21) $$

The proof is quite similar to the one given in [5], so we refer to this work for more details. As a consequence in our case, and from a physical point of view, one may propose the following definition of the coupled problem:

**Definition 1.** The coupled problem is defined by systems (17) and (18), with the following interface conditions on $\gamma$:

$$ \mu \partial_x \mathbf{u} - u_0 \mathbf{u} - g \begin{pmatrix} \zeta_0 \end{pmatrix} = \mu \partial_x \bar{\mathbf{U}}_h - u_0 \bar{\mathbf{U}}_h - g \begin{pmatrix} \zeta_+ \\ 0 \end{pmatrix} \qquad (22) $$

$$ u_0 \zeta_0 + H u = u_0 \zeta_+ + H \bar{\mathbf{U}}_h \qquad (23) $$

# 3 Schwarz multidimensional coupling algorithm

Several approaches can be used to couple different models: variational, algebraic or domain decomposition methods. We choose here a Schwarz domain decomposition approach, which has several practical advantages (simplicity for set-up, few changes in the numerical codes). The Schwarz domain decomposition methods have been generalized to the case of heterogeneous spatial dimension coupling in [22].

## 3.1 Schwarz coupling algorithm

After having defined the coupling notion, we now set up a Schwarz coupling algorithm and solve it. We propose the following algorithm:

**Initialization** : for $\mathbf{U}_h^0$ and $\zeta_+^0$ given

At each iteration $k$ ($k \ge 0$), solve:

$$ \left\{ \begin{array}{ll} \mathcal{L}_{LSW} (\mathbf{u}^{k+1}, \zeta_-^{k+1}) & = \quad 0 \quad \text{in} \quad (\omega^- \times (0, T))^2 \\ \mathcal{B}_-^{\text{out}} (\mathbf{u}^{k+1}, \zeta_-^{k+1}) & = \quad G_-^{\text{out}} \quad \text{on} \quad (\partial \omega_-^{\text{out}} \times (0, T))^2 \\ \mathcal{B}_- (\mathbf{u}^{k+1}, \zeta_-^{k+1}) & = \quad \mathcal{B}_- (\mathcal{R}(\mathbf{U}_h^k, \zeta_+^k)) \quad \text{on} \quad (\gamma \times (0, T))^2 \\ \mathbf{u}^{k+1}(., 0) = \bar{\mathbf{U}}_-^{\text{ini}} \quad \text{in} \quad \omega^- \quad \text{and} \quad \zeta_-^{k+1}(., 0) = \zeta_-^{\text{ini}} \quad \text{in} \quad \omega^- \end{array} \right. $$

then solve:

$$ \left\{ \begin{array}{ll} \mathcal{L}_{LHNS} (\mathbf{U}_h^{k+1}, \zeta_+^{k+1}) & = \quad 0 \quad \text{in} \quad (\Omega^+ \times (0, T)) \times (\omega^+ \times (0, T)) \\ \mathcal{B}_+^{\text{out}} (\mathbf{U}_h^{k+1}, \zeta_+^{k+1}) & = \quad G_+^{\text{out}} \quad \text{on} \quad (\partial \Omega_+^{\text{out}} \times (0, T)) \times (\partial \omega_+^{\text{out}} \times (0, T)) \\ \mathcal{B}_+ (\mathbf{U}_h^{k+1}, \zeta_+^{k+1}) & = \quad \mathcal{B}_+ (\mathcal{E}(\mathbf{u}^{k+1}, \zeta_-^{k+1})) \quad \text{on} \quad (\Gamma \times (0, T)) \times (\gamma \times (0, T)) \\ \mathbf{U}_h^{k+1}(., 0) = \mathbf{U}_+^{\text{ini}} \quad \text{in} \quad \Omega^+ \quad \text{and} \quad \zeta_+^{k+1}(., 0) = \zeta_+^{\text{ini}} \quad \text{in} \quad \omega^+ \end{array} \right. $$

<caption>Table 1: Schwarz multidimensional algorithm</caption>
The operator $\mathcal{L}_{LSW}$ denotes the set of equations (17a) and (17b) and the operator $\mathcal{L}_{LHNS}$ denotes the set of equations (18a) and (18b). The exterior boundaries are $\partial\omega_{out}^- = \partial\omega^- \setminus \gamma$, $\partial\Omega_{out}^+ = \partial\Omega^+ \setminus \Gamma$ and the operators $\mathcal{B}_{-}^{out}$ and $\mathcal{B}_{+}^{out}$ denote the exterior boundary conditions. The operators $G_{-}^{out}$ and $G_{+}^{out}$ take the values 0 or $(\bar{U}_{-}^{d}, \zeta_{-}^{d})$ and $(\bar{U}_{+}^{d}, \zeta_{+}^{d})$ depending on if we have Neumann or Dirichlet conditions.

The operators $\mathcal{B}_{-}$ and $\mathcal{B}_{+}$ will be determined such that the algorithm converges and the physical transmission conditions (22) and (23) are satisfied. Their choice is crucial in order to set up an efficient coupling algorithm with optimal convergence.

As in [3] and [22] we can define the operators $\mathcal{R}$ and $\mathcal{E}$ by:

$$ \begin{array}{l} \mathcal{R}: \Lambda_{3D} \times (0, T) \rightarrow \Lambda_{2D} \times (0, T) \quad \text{and} \quad \mathcal{E}: \Lambda_{2D} \times (0, T) \rightarrow \Lambda_{3D} \times (0, T) \\ (\mathbf{U})|_{\Gamma \times (0,T)} \mapsto (\mathcal{R}\mathbf{U})|_{\gamma \times (0,T)} \quad \mathbf{u}|_{\gamma \times (0,T)} \mapsto (\mathcal{E}\mathbf{u})|_{\Gamma \times (0,T)} \end{array} $$

where $\Lambda_{2D}$ and $\Lambda_{3D}$ denote the trace spaces on the interface $\gamma$ for the 2D spatial functions and on $\Gamma$ for the 3D spatial functions. In view of the derivation of the linearized shallow water system from linearized Navier-Stokes equations, we define the operator $\mathcal{R}$ as the vertical average:

$$ \begin{array}{l} \mathcal{R}: \Lambda_{3D} \times (0, T) \rightarrow \Lambda_{2D} \times (0, T) \\ (\mathbf{U}_h)|_{\Gamma \times (0,T)} \mapsto \frac{1}{H} \int_{-H}^{0} (\mathbf{U}_h)|_{\Gamma \times (0,T)} dz \end{array} $$

that is $\mathcal{R}(\mathbf{U}_h)(0, y, t) = \frac{1}{H} \int_{-H}^{0} \mathbf{U}_h(0, y, z, t) dz$.

Now, in view of the second order relations (15) and (16), we may define the operator $\mathcal{E}$ by:

$$ \begin{array}{l} \mathcal{E}: \Lambda_{2D} \times (0, T) \rightarrow \Lambda_{3D} \times (0, T) \\ (\mathbf{u})|_{\gamma \times (0,T)} \mapsto (\mathbf{u})|_{\Gamma \times (0,T)} \end{array} $$

that is $\mathcal{E}(\mathbf{u})(0, y, z, t) = \mathbf{u}(0, y, t)$. In other words, as we assumed a frictionless condition at the bottom, the quantities coming from the 2D model are extended uniformly on the vertical through the interface $\Gamma$.

## 3.2 Rewriting the Schwarz algorithm

In this section we study the convergence of the Schwarz algorithm 1. But in order to simplify the theoretical study of the coupling algorithm, we first rewrite it by decomposing the 3D velocity into vertical modes, which implies to find the spectrum of the operator $-\partial_z^2$ in $[-H, 0]$ with homogeneous Neumann boundary conditions (see [4], [19] and [20]).

Then we look for $\mathbf{U}_h$ under the form:

$$ \mathbf{U}_h(x, y, z, t) = \sum_{n=0}^{\infty} \mathbf{U}_h^n(x, y, t)w_n(z) = \bar{\mathbf{U}}_h(x, y, t) + \sum_{n=1}^{\infty} \mathbf{U}_h^n(x, y, t)w_n(z) $$

where $w_n(z) = \alpha_n \cos(\frac{n\pi z}{H})$ with $\alpha_0 = 1$ and $\alpha_n = \sqrt{2} (n > 0)$, and $\bar{\mathbf{U}}_h(x, y, t) = \frac{1}{H} \int_{-H}^{0} \mathbf{U}_h(x, y, z, t) dz$ (the exponent $n$ denotes here the rank of the vertical mode, not to be confused with Schwarz iteration index $k$).

By injecting this decomposition into (18), then multiplying by $w_n$ and integrating between $-H$ and 0, we obtain (using the fact that $(w_n)_{n \ge 0}$ is an orthonormal basis):

* For $n = 0$, the first vertical mode (also called barotropic mode) $\mathbf{U}_h^0 = \bar{\mathbf{U}}_h$ coupled with the free surface is
solution of the 2D linearized shallow water system in $\omega^+$:

$$ egin{cases} 
\partial_t \bar{\mathbf{U}}_h + (\mathbf{U}_0 \cdot \nabla_h) \bar{\mathbf{U}}_h + g \nabla_h \zeta - \mu \Delta_h \bar{\mathbf{U}}_h = 0 & \text{in } \omega^+ \times (0, T) \\[1.5em] 
\partial_t \zeta + H \operatorname{div}_h (\bar{\mathbf{U}}_h) + \mathbf{U}_0 \cdot \nabla_h \zeta = 0 & \text{in } \omega^+ \times (0, T) \\[1.5em] 
(\bar{\mathbf{U}}_h, \zeta) = (\bar{\mathbf{U}}_+^d, \zeta_+^d) & \text{on } (\partial \omega_{out}^+ \times (0, T))^2 \\[1.5em] 
(\bar{\mathbf{U}}_h(\cdot, 0), \zeta(\cdot, 0)) = (\bar{\mathbf{U}}_+^{ini}, \zeta_+^{ini}) & \text{in } (\omega^+)^2 
\end{cases} 
\tag{24} $$

• For $n \ge 1$ (baroclinic modes):

$$ \partial_t \mathbf{U}_h^n + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{U}_h^n - \mu \Delta_h \mathbf{U}_h^n + \frac{\mu(n\pi)^2}{H^2} \mathbf{U}_h^n = 0 \quad \text{in } \Omega^+ \times (0, T) \tag{25} $$

The sum of the baroclinic modes $\mathbf{U}_b = \mathbf{U}_h - \bar{\mathbf{U}}_h$ is solution of the system:

$$ egin{cases} 
\frac{\partial \mathbf{U}_b}{\partial t} + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{U}_b - \mu \Delta \mathbf{U}_b = 0 & \text{in } \Omega^+ \times (0, T) \\[1.5em] 
\mathcal{B}_+^{\text{out}}(\mathbf{U}_b) = G_+^{\text{ext}} - \bar{G}_+^{\text{ext}} & \text{on } \partial \Omega_{\text{out}}^+ \times (0, T) \\[1.5em] 
\mathbf{U}_b(\cdot, 0) = \mathbf{U}_+^{\text{ini}} - \bar{\mathbf{U}}_+^{\text{ini}} & \text{in } \Omega^+ 
\end{cases} 
\tag{26} $$

We can then rewrite the Schwarz algorithm 1 as follows:

**Initialization**: $(\mathbf{U}_h)^0$ and $\zeta_+^0$ given

At each iteration $k$ ($k \ge 0$), solve:

$$ egin{cases} 
\mathcal{L}_{LSW} (\mathbf{u}^{k+1}, \zeta_{-}^{k+1}) = 0 & \text{in } (\omega^{-} \times (0, T))^2 \\[1.5em] 
\mathcal{B}_{-}^{\text{out}} (\mathbf{u}^{k+1}, \zeta_{-}^{k+1}) = G_{-}^{\text{out}} & \text{on } (\partial \omega_{\text{out}}^{-} \times (0, T))^2 \\[1.5em] 
\mathcal{B}_{-} (\mathbf{u}^{k+1}, \zeta_{-}^{k+1}) = \mathcal{B}_{-} (\bar{\mathbf{U}}_{h}^{k}, \zeta_{+}^{k}) & \text{on } (\gamma \times (0, T))^2 \\[1.5em] 
\mathbf{u}^{k+1}(\cdot, 0) = \bar{\mathbf{U}}_{-}^{\text{ini}} & \text{in } \omega^{-} \text{ and } \zeta_{-}^{k+1}(\cdot, 0) = \zeta_{-}^{\text{ini}} \text{ in } \omega^{-} 
\end{cases} $$

then solve

$$ egin{cases} 
\mathcal{L}_{LSW} (\bar{\mathbf{U}}_h^{k+1}, \zeta_+^{k+1}) = 0 & \text{in } (\omega^+ \times (0, T))^2 \\[1.5em] 
\mathcal{B}_+^{\text{out}} (\bar{\mathbf{U}}_h^{k+1}, \zeta_+^{k+1}) = \bar{G}_+^{\text{out}} & \text{on } (\partial \omega_{\text{out}}^+ \times (0, T))^2 \\[1.5em] 
\mathcal{B}_+^0 (\bar{\mathbf{U}}_h^{k+1}, \zeta_+^{k+1}) = \mathcal{B}_+^0 (\mathbf{u}^{k+1}, \zeta_-^{k+1}) & \text{on } (\gamma \times (0, T))^2 \\[1.5em] 
\bar{\mathbf{U}}_h^{k+1}(\cdot, 0) = \bar{\mathbf{U}}_+^{\text{ini}} & \text{in } \omega^+ \text{ and } \zeta_+^{k+1}(\cdot, 0) = \zeta_+^{\text{ini}} \text{ in } \omega^+ 
\end{cases} $$

and

$$ egin{cases} 
\mathcal{L}_{CD} ((\mathbf{U}_b)^{k+1}) = 0 & \text{in } \Omega^+ \times (0, T) \\[1.5em] 
\mathcal{B}_+^{\text{out}} ((\mathbf{U}_b)^{k+1}) = G_+^{\text{out}} - \bar{G}_+^{\text{out}} & \text{on } \partial \Omega_{\text{out}}^+ \times (0, T) \\[1.5em] 
\mathcal{B}_+^0 ((\mathbf{U}_b)^{k+1}) = 0 & \text{on } \Gamma \times (0, T) \\[1.5em] 
(\mathbf{U}_b)^{k+1}(\cdot, 0) = \mathbf{U}_+^{\text{ini}} - \bar{\mathbf{U}}_+^{\text{ini}} & \text{in } \Omega^+ 
\end{cases} $$

<caption>Table 2: Schwarz multidimensional coupling algorithm - Version 2</caption>
where $\mathcal{B}_+^0$ denotes the restriction of the operator $\mathcal{B}_+$ for 2D spatial functions, $\mathcal{B}_+'$ is obtained from $\mathcal{B}_+$ by vanishing all terms containing $\zeta_+$, and $\mathcal{L}_{CD}$ denotes the convection-diffusion operator $\partial_t + (\mathbf{U}_0 \cdot \nabla_h) - \mu\Delta$. Due to this new form of the algorithm, we have then the following convergence result as a direct consequence of the algorithm 2:

**Lemme 3.** The coupling algorithm 1 converges if and only if the classic domain decomposition algorithm of the linearized shallow water system defined by:

$\bar{\mathbf{u}}_+^0$ and $\zeta_+^0$ given

$$ \left\{ \begin{array}{lclcl} \mathcal{L}_{LSW} (\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) & = & 0 & \text{in} & (\omega^- \times (0, T))^2 \\ \mathcal{B}_-^{\text{out}} (\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) & = & \mathcal{C}_-^{\text{out}} & \text{on} & (\partial\omega_{\text{out}}^- \times (0, T))^2 \\ \mathcal{B}_- (\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) & = & \mathcal{B}_- (\mathbf{u}_+^k, \zeta_+^k) & \text{on} & (\gamma \times (0, T))^2 \\ (\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) (., 0) & = & (\bar{\mathbf{U}}_-^{\text{ini}}, \zeta_-^{\text{ini}}) & \text{in} & \omega^- \end{array} \right. \qquad (27) $$

then

$$ \left\{ \begin{array}{lclcl} \mathcal{L}_{LSW} (\bar{\mathbf{U}}_h^{k+1}, \zeta_+^{k+1}) & = & 0 & \text{in} & (\omega^+ \times (0, T))^2 \\ \mathcal{B}_+^{\text{out}} (\bar{\mathbf{U}}_h^{k+1}, \zeta_+^{k+1}) & = & \bar{\mathcal{G}}_+^{\text{out}} & \text{on} & \partial(\partial\omega_{\text{out}}^+ \times (0, T))^2 \\ \mathcal{B}_+^0 (\bar{\mathbf{U}}_h^{k+1}, \zeta_+^{k+1}) & = & \mathcal{B}_+^0 (\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) & \text{on} & (\gamma \times (0, T))^2 \\ (\bar{\mathbf{U}}_h^{k+1}, \zeta_+^{k+1}) (., 0) & = & (\bar{\mathbf{U}}_+^{\text{ini}}, \zeta_+^{\text{ini}}) & \text{in} & \omega^+ \end{array} \right. \qquad (28) $$

converges, and the baroclinic velocity $\mathbf{U}_b = \mathbf{U}_h - \bar{\mathbf{U}}_h$ is solution of the system:

$$ \left\{ \begin{array}{lclcl} \frac{\partial \mathbf{U}_b}{\partial t} + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{U}_b - \mu \Delta \mathbf{U}_b & = & 0 & \text{in} & \Omega^+ \times (0, T) \\ \mathcal{B}_+^{\text{out}} (\mathbf{U}_b) & = & \mathcal{G}_+^{\text{out}} - \bar{\mathcal{G}}_+^{\text{out}} & \text{in} & \partial \Omega_{\text{out}}^+ \times (0, T) \\ \mathcal{B}_+^\prime (\mathbf{U}_b) & = & 0 & \text{on} & \Gamma \times (0, T) \\ \mathbf{U}_b (., 0) & = & \mathbf{U}_+^{\text{ini}} - \bar{\mathbf{U}}_+^{\text{ini}} & \text{in} & \Omega^+ \end{array} \right. $$

**Remark 2.** Note that this result is no longer true if there is an additional bottom friction term, since the vertical modes of the operator $-\partial_z^2$ cannot be decoupled in this case, see [2].

### 3.3 Study of the coupling algorithm with “Robin boundary” conditions

We investigated in [5] the convergence of the domain decomposition algorithm for the linearized shallow water system. We first studied the approximation of the optimal transmission conditions by assuming a large Reynolds number and a small ratio aspect $\epsilon$, as in section 2.3. Unfortunately we were not able to find a useful approximation of these operators. Therefore we studied the domain decomposition algorithm with “Robin” boundary conditions. More precisely we studied the convergence of the algorithm with the following boundary conditions (for $u_0 > 0$):

$$ \mathcal{B}_- (\mathbf{u}, \zeta) = \begin{pmatrix} \mu \frac{\partial u}{\partial x} - g\zeta + \frac{(\lambda - u_0)}{2} u \\ \mu \frac{\partial v}{\partial x} + \frac{(\lambda - u_0)}{2} v \end{pmatrix} \qquad (29) $$
and

$$ \mathcal{B}_+^0(\mathbf{u}, \zeta) = \begin{pmatrix} -\mu \frac{\partial u}{\partial x} + g\zeta + \frac{(\lambda + u_0)}{2} u \\ -\mu \frac{\partial v}{\partial x} + \frac{(\lambda + u_0)}{2} v \\ u_0 \zeta \end{pmatrix} \qquad (30) $$

where $\lambda$ is a positive constant.

We can extend the operator $\mathcal{B}_+$ defined by (30) in the case of 3D spatial functions:

$$ \mathcal{B}_+(\mathbf{U}_h, \zeta) = \begin{pmatrix} -\mu \frac{\partial U}{\partial x} + g\zeta + \frac{(\lambda + u_0)}{2} U \\ -\mu \frac{\partial V}{\partial x} + \frac{(\lambda + u_0)}{2} V \\ u_0 \zeta \end{pmatrix} \qquad (31) $$

and redefine the Schwarz algorithm 2 with the operators (29) and (31).

**Proposition 1.** The coupling algorithms 1 and 2 defined with the operators (29) and (30) are well-posed. The sequences $(\mathbf{u}^{k+1}, \zeta_{-}^{k+1})$ and $(\mathbf{U}_h^{k+1}, \zeta_{+}^{k+1})$ defined by the algorithm 2 converge to $(\mathbf{u}|_{\omega-}, \zeta|_{\omega-})$ and $(\mathbf{U}_h^{\lambda}, \zeta_{+}^{\lambda})$ respectively in $(C(0, T; L^2(\omega^{-}, \mathbb{R}^2) \cap L^2(0, T; H^1(\omega^{-}, \mathbb{R}^2))) \times (L^2(\omega^{-} \times (0, T)) \cap C(0, T; L^2(\omega^{-})))$ and $(C(0, T; L^2(\Omega^{+}, \mathbb{R}^2) \cap L^2(0, T; H^1(\Omega^{+}, \mathbb{R}^2))) \times (L^2(\omega^{+} \times (0, T)) \cap C(0, T; L^2(\omega^{+})))$, where $\mathbf{u}|_{\omega-}$ et $\zeta|_{\omega-}$ denote the restriction in $\omega^{-}$ of $(\mathbf{u}, \zeta)$ solution of the linearized shallow water system throughout the domain $\omega$. The sequence of barotropic velocities $(\overline{\mathbf{U}}_h^{k+1})_{k \ge 0}$ and the sequence $(\zeta_{+}^{k+1})_{k \ge 0}$ converge respectively to $\mathbf{u}|_{\omega+}$ and $\zeta|_{\omega+}$. At convergence, the physical constraints (22) and (23) are satisfied.

**Proof** The convergence of the coupling algorithm results from Lemma 3.

The proof of the well-posedness of the two systems in algorithms (1) and (2) is quite similar to the one given in [2] et [5]. We will only give here its outline.

For $(\mathbf{U}_h)^0$ and $\zeta_+^0$ given and for all $k \ge 0$, the study of the coupling algorithm 1 is equivalent to the study of the following systems:

* In the domain $\omega^{-}$, we solve the parabolic system:

$$ \partial_t \mathbf{u}^{k+1} + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{u}^{k+1} - \mu \Delta_h \mathbf{u}^{k+1} = -g \nabla_h \zeta_{-}^{k+1} \quad \text{in } \omega^{-} \times (0, T) \qquad (32a) $$

$$ \mathcal{B}_-(\mathbf{u}^{k+1}, \zeta_{-}^{k+1}) = \mathcal{B}_-(\overline{\mathbf{U}}_h^k, \zeta_{+}^k) \quad \text{on } \gamma \times (0, T) \qquad (32b) $$

$$ \mathbf{u}^{k+1}(., 0) = \overline{\mathbf{U}}_{-}^{\text{ini}} \quad \text{in } \omega^{-} \qquad (32c) $$

and the transport equation:

$$ \partial_t \zeta_{-}^{k+1} + \mathbf{U}_0 \cdot \nabla_h \zeta_{-}^{k+1} = -H \operatorname{div}_h(\mathbf{u}^{k+1}) \quad \text{in } \omega^{-} \times (0, T) \qquad (33a) $$

$$ \zeta_{-}^{k+1}(., 0) = \zeta_{-}^{\text{ini}} \quad \text{in } \omega^{-} \qquad (33b) $$

We do not consider here a boundary condition at $x = L_0 = 0$ for the transport equation as we assumed $u_0 > 0$.
• In the domain $\Omega^+$, we solve the parabolic system:

$$ \left\{ \partial_t \mathbf{U}_h^{k+1} + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{U}_h^{k+1} - \mu \Delta_h \mathbf{U}_h^{k+1} = -g \nabla_h \zeta_+^{k+1} \quad \text{in } \Omega^+ \times (0, T) \quad (34a) \right. $$

$$ \mu \partial_z \mathbf{U}_h = 0 \quad \text{at} \quad z=0 \qquad (34b) $$

$$ \mu \partial_z \mathbf{U}_h = 0 \quad \text{at} \quad z = -H \qquad (34c) $$

$$ \mathcal{B}_+^u(\mathbf{U}_h^{k+1}, \zeta_+^{k+1}) = \mathcal{B}_+^u(\mathbf{u}^{k+1}, \zeta_-^{k+1}) \quad \text{on} \quad \Gamma \times (0, T) \qquad (34d) $$

$$ \mathbf{U}_h^{k+1}(., 0) = \mathbf{U}_+^{\text{ini}} \quad \text{in} \quad \Omega^+ \qquad (34e) $$

and the transport equation:

$$ \left\{ \partial_t \zeta_+^{k+1} + \mathbf{U}_0 \cdot \nabla_h \zeta_+^{k+1} = -H \operatorname{div}_h (\bar{\mathbf{U}}_h^{k+1}) \quad \text{in} \quad \omega^+ \times (0, T) \quad (35a) \right. $$

$$ \mathcal{B}_+^\zeta (\mathbf{U}_h^{k+1}, \zeta_+^{k+1}) = \mathcal{B}_+^\zeta (\mathbf{u}^{k+1}, \zeta_-^{k+1}) \quad \text{on} \quad \gamma \times (0, T) \qquad (35b) $$

$$ \zeta_+^{k+1}(., 0) = \zeta_+^{\text{ini}} \quad \text{in} \quad \omega^+ \qquad (35c) $$

Then we study in each subdomain the parabolic system with a prescribed water height and the transport equation with a prescribed velocity. Finally one can use the fixed point theorem to conclude that the Schwarz algorithm is well-posed (See [2] and [5] for details). □

## 4 Modeling error

Unlike the usual case of domain decomposition, at convergence of the Schwarz algorithm, we do not have $\mathbf{U}_h^\lambda$ the limit of $(\mathbf{U}_h^k)$ equal to $\mathbf{U}_{h|\Omega^+}^G$, where $\mathbf{U}_h^G$ is the solution of system (7).

Due to the way the coupling algorithm is rewritten, the error is obviously contained in the baroclinic mode. We will thus investigate, as in [22], the amplitude of the modeling error between the baroclinic mode $\mathbf{U}_b^\lambda$ of the coupled solution and the baroclinic mode $\mathbf{U}_b$ of the global reference solution. As in [22], the choice of Robin-like operators allows its control. Actually we have the following result:

**Theorem 1.** For every $\lambda > 0$, let $\mathbf{U}_b^\lambda$ denote the baroclinic mode of the coupled solution. If $L_0 < L_1$, then there exists $M = M(\lambda)$ such that

$$ \|\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G - \mathbf{U}_b^\lambda\|_{L^2((0,T);H^1(\Omega^+)\times\Omega^+)} + \|\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G - \mathbf{U}_b^\lambda\|_{C([0,T];L^2(\Omega^+)\times\Omega^+)} \leq M(\lambda)\varepsilon\sqrt{1+\delta^2} \quad (36) $$

where $\delta = \frac{L_1}{L_1-L_0}$

**Proof** Let introduce in the sequel the function $\mathbf{E}_+^\lambda$ defined by $\mathbf{E}_+^\lambda = \mathbf{U}_h^G - \bar{\mathbf{U}}_h^G - \mathbf{U}_b^\lambda$. Therefore $\mathbf{E}_+^\lambda$ is solution of the system:

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{E}_+^\lambda}{\partial t} + (\mathbf{U}_0 \cdot \nabla_h) \mathbf{E}_+^\lambda - \mu \Delta \mathbf{E}_+^\lambda = 0 \quad \text{in } \Omega^+ \times (0, T) \\[2ex] \displaystyle \mathcal{B}_+^{\text{out}} (\mathbf{E}_+^\lambda) = 0 \quad \text{in } \partial \Omega_{\text{out}}^+ \times (0, T) \\[2ex] \displaystyle \mathcal{B}_+^\prime (\mathbf{E}_+^\lambda) = \mathcal{B}_+^\prime (\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G) \quad \text{on } \Gamma \times (0, T) \\[2ex] \displaystyle (\mathbf{E})_+^\lambda(., 0) = 0 \quad \text{in } \Omega^+ \end{array} \right. \qquad (37) $$
where $U_0 = (u_0, v_0)$ with $u_0 > 0$ and the interface operator $B'_+$ is defined by:

$$B'_+(U) = \begin{pmatrix} -\mu \frac{\partial U}{\partial x} + \frac{(\lambda + u_0)}{2} U \\ -\mu \frac{\partial V}{\partial x} + \frac{(\lambda + u_0)}{2} V \end{pmatrix}$$

where $U = (U, V)$. Multiplying (37) by $E_+^\lambda$ and integrating over $\Omega^+$ leads to:

$$\int_{\Omega^+} \frac{\partial \mathbf{E}_+^\lambda}{\partial t} \cdot \mathbf{E}_+^\lambda + \int_{\Omega^+} (\mathbf{U}_0 \cdot \nabla_h) \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda - \mu \int_{\Omega^+} \Delta \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda = 0 \quad (38)$$

Integrating by parts and using the following relation:

$$-\mu \int_{\Omega^+} \Delta \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda = \mu \int_{\Omega^+} \nabla_h \mathbf{E}_+^\lambda : \nabla_h \mathbf{E}_+^\lambda + \mu \int_{\Omega^+} \frac{\partial \mathbf{E}_+^\lambda}{\partial z} \cdot \frac{\partial \mathbf{E}_+^\lambda}{\partial z} - \mu \int_{\Gamma} \frac{\partial \mathbf{E}_+^\lambda}{\partial \mathbf{n}^+} \cdot \mathbf{E}_+^\lambda$$

where $\mathbf{n}^+ = (n_1^+, n_2^+, n_3^+)^T$ denotes the unit outward vector normal to $\Omega^+$, equation (38) then becomes:

$$\frac{1}{2} \frac{d}{dt} \| \mathbf{E}_+^\lambda \|_{\Omega^+}^2 + \int_{\Omega^+} (\mathbf{U}_0 \cdot \nabla_h) \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda + \mu \| \nabla_h \mathbf{E}_+^\lambda \|_{\Omega^+}^2 + \mu \left\| \frac{\partial \mathbf{E}_+^\lambda}{\partial z} \right\|_{\Omega^+}^2 - \mu \int_{\Gamma} \frac{\partial \mathbf{E}_+^\lambda}{\partial \mathbf{n}^+} \cdot \mathbf{E}_+^\lambda = 0$$

Now, one has:

$$\int_{\Omega^+} (\mathbf{U}_0 \cdot \nabla_h) \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda = - \int_{\Omega^+} (\mathbf{U}_0 \cdot \nabla_h) \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda + \int_{\Gamma} u_0 \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda n_1^+ + \int_{\Gamma} v_0 \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda n_2^+$$

and as $\mathbf{n}^+ = (-1, 0, 0)$, this implies:

$$\int_{\Omega^+} (\mathbf{U}_0 \cdot \nabla_h) \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda = -\frac{1}{2} \int_{\Gamma} u_0 \mathbf{E}_+^\lambda \cdot \mathbf{E}_+^\lambda n_1^+$$

Equation (38) leads to:

$$\frac{1}{2} \frac{d}{dt} \| \mathbf{E}_+^\lambda \|_{\Omega^+}^2 + \mu \| \nabla_h \mathbf{E}_+^\lambda \|_{\Omega^+}^2 + \mu \left\| \frac{\partial \mathbf{E}_+^\lambda}{\partial z} \right\|_{\Omega^+}^2 = \int_{\Gamma} \left( \frac{1}{2} u_0 \mathbf{E}_+^\lambda - \mu \frac{\partial \mathbf{E}_+^\lambda}{\partial x} \right) \cdot \mathbf{E}_+^\lambda$$

Due to the boundary condition on $\Gamma$, one has:

$$\frac{1}{2} \frac{d}{dt} \| \mathbf{E}_+^\lambda \|_{\Omega^+}^2 + \mu \| \nabla_h \mathbf{E}_+^\lambda \|_{\Omega^+}^2 + \mu \left\| \frac{\partial \mathbf{E}_+^\lambda}{\partial z} \right\|_{\Omega^+}^2 = \int_{\Gamma} \left( B'_+(\mathbf{E}_+^\lambda) - \frac{\lambda}{2} \mathbf{E}_+^\lambda \right) \cdot \mathbf{E}_+^\lambda$$

and therefore:

$$\frac{1}{2} \frac{d}{dt} \| \mathbf{E}_+^\lambda \|_{\Omega^+}^2 + \mu \| \nabla_h \mathbf{E}_+^\lambda \|_{\Omega^+}^2 + \mu \left\| \frac{\partial \mathbf{E}_+^\lambda}{\partial z} \right\|_{\Omega^+}^2 + \frac{\lambda}{2} \| \mathbf{E}_+^\lambda \|_{\Gamma}^2 = \int_{\Gamma} B'_+(\mathbf{E}_+^\lambda) \cdot \mathbf{E}_+^\lambda \quad (39)$$

The right-hand side of equation (39) reads as follows:

$$\int_{\Gamma} B'_+(\mathbf{E}_+^\lambda) \cdot \mathbf{E}_+^\lambda = \int_{\Gamma} B'_+(\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G) \cdot (\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G - \mathbf{U}_b^\lambda) \quad (40)$$

$$= \int_{\Gamma} B'_+(\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G) \cdot (\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G) - \int_{\Gamma} B'_+(\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G) \cdot \bar{\mathbf{U}}_b^\lambda \quad (41)$$

Due to the relations (15) and (16), one has:

$$\mathbf{U}_h^G - \bar{\mathbf{U}}_h^G = O(\varepsilon^2)$$
Then, due again to the relations (15) and (16), it is reasonable to assume that $\mathcal{B}'_+(\mathbf{U}_h^G - \overline{\mathbf{U}}_h^G) = O(1)$, so that one has:

$$ \int_{\Gamma} \mathcal{B}'_+(\mathbf{U}_h^G - \overline{\mathbf{U}}_h^G) \cdot (\mathbf{U}_h^G - \overline{\mathbf{U}}_h^G) = O(\varepsilon^2) $$

Now we will focus on the term $\int_{\Gamma} \mathcal{B}'_+(\mathbf{U}_h^G - \overline{\mathbf{U}}_h^G) \cdot \overline{\mathbf{U}}_h^\lambda$. In the same way and as in [22], if we assume that $L_0 < L_1$, so that the 3D effects are insignificant in $\Omega^+ \cap \{L_0 \le x \le L_1\}$, and applying a similar asymptotic analysis as in the first section to the 3D model, we can deduce:

$$ \mathbf{U}_h^\lambda = \overline{\mathbf{U}}_h^\lambda + O(\varepsilon'^2) $$

where $\varepsilon' = \frac{H}{L_1 - L_0} = \frac{L_1}{L_1 - L_0} \varepsilon$. We denote in the sequel $\delta = \frac{L_1}{L_1 - L_0}$, therefore:

$$ \mathbf{U}_b^\lambda = \mathbf{U}_h^\lambda - \overline{\mathbf{U}}_h^\lambda = O(\delta^2 \varepsilon^2) $$

So that, there exists a positive constant $C_1$ depending continuously on $\lambda$ and $t$ such that:

$$ \frac{1}{2} \frac{d}{dt} \|\mathbf{E}_+^\lambda\|_{\Omega^+}^2 + \mu \|\nabla_h \mathbf{E}_+^\lambda\|_{\Omega^+}^2 + \mu \left\| \frac{\partial \mathbf{E}_+^\lambda}{\partial z} \right\|_{\Omega^+}^2 + \frac{\lambda}{2} \|\mathbf{E}_+^\lambda\|_{\Gamma}^2 \le C_1 (1 + \delta^2) \varepsilon^2 \quad (42) $$

Integrating between 0 and $t$ for $t \in [0; T]$ and using the initial conditions, (42) leads to:

$$ \frac{1}{2} \|\mathbf{E}_+^\lambda\|_{\Omega^+}^2 + \mu \int_0^t \|\nabla_h \mathbf{E}_+^\lambda\|_{\Omega^+}^2 + \mu \int_0^t \left\| \frac{\partial \mathbf{E}_+^\lambda}{\partial z} \right\|_{\Omega^+}^2 + \frac{\lambda}{2} \int_0^t \|\mathbf{E}_+^\lambda\|_{\Gamma}^2 \le C_2(\lambda, t)(1 + \delta^2)\varepsilon^2 \quad (43) $$

and then one has:

$$ \|\mathbf{E}_+^\lambda\|_{\Omega^+}^2 \le 2C_2(\lambda, t)(1 + \delta^2)\varepsilon^2 \quad (44) $$

and we can also deduce that:

$$ \int_0^T \|\mathbf{E}_+^\lambda\|_{\Omega^+}^2 + \int_0^T \|\nabla \mathbf{E}_+^\lambda\|_{\Omega^+}^2 \le C_3(\lambda)(1 + \delta^2)\varepsilon^2 \quad (45) $$

Finally, we can establish the error majoration (36). □

# 5 Numerical schemes and tests

The aim of this section is to illustrate the previous theoretical results with some numerical experiments. In a first step and to illustrate the theoretical convergence result of the Schwarz coupling algorithm, we study the numerical convergence of the domain decomposition method applied to the shallow water equations. Then, in a second step, experiments with the dimensionally heterogeneous coupling method are conducted.

## 5.1 Domain decomposition for the 2D shallow water equations

In this subsection, we numerically study the convergence of the domain decomposition method for 2D linearized shallow water equations. We first present the monodomain solution, i.e. the numerical solution to the system (17) throughout the whole domain $\omega$, and then we numerically set up the Schwarz algorithm and illustrate its convergence.

### 5.1.1 Reference monodomain solution

In this paragraph, we detail the parameters of the numerical monodomain configuration and the resulting simulated solution. This reference solution will allow to estimate the convergence of the domain decomposition algorithm.
**Domain and boundary conditions** The 2D linearized shallow water equations are solved in $\omega = [-L; L] \times [0; L_y]$. Although the mathematical analysis has been developed with a null right-hand side, we consider for the numerical tests a non-zero source term $\tau = (\tau_x, 0)^T$ (wind stress). We impose also homogeneous Dirichlet conditions for the velocity $u_h$ and for the water height $\zeta$ on the boundary $\partial\omega$.

**Numerical schemes** Model equations are discretized using a finite difference scheme on a $N_x \times N_y$ cartesian staggered Arakawa C-grid [1] (velocities are computed on the edges of each cell and water height in the center - see Figure 1). The space steps are $\Delta x$ and $\Delta y$ and the time step is $\Delta t = \frac{T}{N_t}$. The cell $C_{i,j}$ is defined by $C_{i,j} = [x_i; x_{i+1}] \times [y_i; y_{i+1}]$, where $x_i = -L + i\Delta x$ and $y_j = j\Delta y$ for $0 \le i \le N_x$ and $0 \le j \le N_y$. We define also $t_n = n\Delta t$ for $0 \le n \le N_t$. The discrete unknowns are: $u_{i,j}^n \approx u(t_n, x_i, y_j)$, $v_{i,j}^n \approx v(t_n, x_{i+\frac{1}{2}}, y_{j+\frac{1}{2}})$ and $\zeta_{i,j}^n \approx \zeta(t_n, x_{i+\frac{1}{2}}, y_j)$.

<figure><img src="image_2.png" /><figcaption>Figure 1: Staggered grid in space</figcaption></figure>

* Discrete $u$-equation: the equation for $u$

$$ \frac{\partial u}{\partial t} + u_0 \frac{\partial u}{\partial x} + v_0 \frac{\partial u}{\partial y} + g \frac{\partial \zeta}{\partial x} = \mu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) + \tau_x \quad (46) $$

is discretized using standard explicit schemes as follows:

$$ \frac{1}{\Delta t} (u_{i,j}^{n+1} - u_{i,j}^n) + \frac{u_0}{\Delta x} (u_{i,j}^n - u_{i-1,j}^n) + \frac{v_0}{\Delta y} (u_{i,j}^n - u_{i,j-1}^n) + \frac{g}{\Delta x} (\zeta_{i,j}^n - \zeta_{i-1,j}^n) \\ = \frac{\mu}{\Delta x} \left(\frac{\partial u}{\partial x}\right)_{i,j}^n - \frac{\mu}{\Delta x} \left(\frac{\partial u}{\partial x}\right)_{i-1,j}^n + \mu \frac{u_{i,j+1}^n - 2u_{i,j}^n + u_{i,j-1}^n}{\Delta y^2} + (\tau_x)_{i,j} \quad (47) $$

Note that the discretization of $\partial^2 u/\partial x^2$ can also be written as

$$ \frac{u_{i+1,j}^n - 2u_{i,j}^n + u_{i-1,j}^n}{\Delta x^2} = \frac{1}{\Delta x} \left( \left( \frac{\partial u}{\partial x} \right)_{i,j}^n - \left( \frac{\partial u}{\partial x} \right)_{i-1,j}^n \right) \quad (48) $$

where $(\partial u/\partial x)_{i,j}^n = \frac{u_{i,j}^n - u_{i-1,j}^n}{\Delta x} \approx \frac{\partial u}{\partial x}(t_n, x_i, y_j)$. This will be used later for the discretization of the Schwarz algorithm.
* Similarly the equation for $v$

$$ \frac{\partial v}{\partial t} + u_0 \frac{\partial v}{\partial x} + v_0 \frac{\partial v}{\partial y} + g \frac{\partial \zeta}{\partial y} = \mu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} \right) \quad (49) $$

is discretized as

$$ \begin{aligned} & \frac{1}{\delta t} (v_{i,j}^{n+1} - v_{i,j}^n) + \frac{u_0}{\Delta x} (v_{i,j}^n - v_{i-1,j}^n) + \frac{v_0}{2\Delta y} (v_{i,j+1}^n - v_{i,j-1}^n) + \frac{g}{2\Delta y} (\zeta_{i,j+1}^n - \zeta_{i,j-1}^n) \\ & = \frac{\mu}{\Delta x} \left(\frac{\partial v}{\partial x}\right)_{i,j}^n - \frac{\mu}{\Delta x} \left(\frac{\partial v}{\partial x}\right)_{i,j}^n + \mu \frac{v_{i,j+1}^n - 2v_{i,j}^n + v_{i,j-1}^n}{\Delta y^2} \end{aligned} \quad (50) $$

* Finally, the discrete equation for $\zeta$ is:

$$ \frac{1}{\delta t} (\zeta_{i,j}^{n+1} - \zeta_{i,j}^n) + \frac{H}{\Delta x} (u_{i,j}^n - u_{i-1,j}^n) + \frac{H}{\Delta y} (v_{i,j}^n - v_{i-1,j}^n) + \frac{u_0}{\Delta x} (\zeta_{i,j}^n - \zeta_{i-1,j}^n) + \frac{v_0}{\Delta y} (\zeta_{i,j}^n - \zeta_{i,j-1}^n) = 0 \quad (51) $$

**Monodomain reference solution.** For this reference simulation, we choose a 600 km × 120 km rectangular basin. The values of the physical parameters are $\mu = 1100 m^2 s^{-1}$, $\tau_x(x, y) = -\tau_0 \cos(\frac{2\pi}{L_y}y)$ with $\tau_0 = 10^{-5} ms^{-2}$, $u_0 = v_0 = 0.0075 ms^{-1}$ and $T = 3.5$ days. The numerical scheme described above is used with the space steps $\Delta x = \Delta y = 1.2 km$ and the time space $\Delta t = 1min$. The simulation is started at rest $u_0 = v_0 = 0$. The simulation lasts for 3.5 days. The zonal velocity $u$ of this reference solution at $t = 3.5$ days is displayed in Figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2: Zonal velocity $u$ of the reference monodomain solution at $t = 3.5$ days</figcaption></figure>

### 5.1.2 Numerical scheme for the domain decomposition algorithm

We now numerically study the domain decomposition algorithm (27) and (28) with Robin-like boundary conditions at the interface. Let split $\omega = [-L; L] \times [0; L_y]$ into two non-overlapping subdomains $\omega^- = [-L; 0] \times [0; L_y]$ and $\omega^+ = [0; L] \times [0; L_y]$. Boundary conditions for $(u, v, \zeta)$ are exchanged at the interface $\gamma$ located at $x = 0$. The grids of $\omega^-$ and $\omega^+$ have $N_x^-$ and $N_y^-$ and $N_x^+$ and $N_y^+$ points respectively. As the grids are staggered, the physical boundary on the $y$-axis is located on $u$-points, while it is located on $v$-points on the $x$-axis, resulting in the need for ghost points on other axes. The boundaries are respectively called north, south, west and east for each
subdomain. Since the two problems are similar, we describe only the discretization in $\omega^-$. For the domain $\omega^-$, and at each step of the algorithm, one has to solve the following problem:

$$ 
\left\{ 
\begin{array}{lcl}
\mathcal{L}_{LSW}(\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) & = & \tau \quad \text{in} \quad (\omega^- \times (0, T))^2 \\
\mathcal{B}_-^{\text{out}}(\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) & = & G_-^{\text{out}} \quad \text{on} \quad (\omega_-^{\text{out}} \times (0, T))^2 \\
\mathcal{B}_- (\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) & = & \mathcal{B}_- (\mathbf{u}_+^k, \zeta_+^k) \quad \text{on} \quad (\gamma \times (0, T))^2 \\
(\mathbf{u}_-^{k+1}, \zeta_-^{k+1}) (\cdot, 0) & = & (\bar{\mathbf{U}}_-^{\text{ini}}, \zeta_-^{\text{ini}}) \quad \text{in} \quad \omega^-
\end{array}
\right. 
$$

The index $k$ refers to the Schwarz iterations. As in [16], we focus on two key points: the discretization of the boundary condition $\mathcal{B}_- (\mathbf{u}_-, \zeta_-) = \mathbf{G}$ for a given $\mathbf{G}$, and the extraction of the quantity $\mathcal{B}_- (\mathbf{u}_+, \zeta_+^k)$ from $\omega^+$. Recall that the transmission condition on $\gamma$ reads:

$$ 
\left\{ 
\begin{aligned}
& \mu \frac{\partial u_-^{k+1}}{\partial x} - g \zeta_-^{k+1} + \frac{1}{2}(\lambda - u_0) u_-^{k+1} = \mu \frac{\partial u_+^k}{\partial x} - g \zeta_+^k + \frac{1}{2}(\lambda - u_0) u_+^k \\
& \mu \frac{\partial v_-^{k+1}}{\partial x} + \frac{1}{2}(\lambda - u_0) v_-^{k+1} = \mu \frac{\partial v_+^k}{\partial x} + \frac{1}{2}(\lambda - u_0) v_+^k
\end{aligned}
\right. 
\quad (52) 
$$

The discrete value of $u$ and $v$ at the interface $\gamma$ should satisfy simultaneously this transmission condition and the discrete interior equations for $u$ and $v$. Let consider the discretization of the first equation of (52):

$$ \mu \left( \frac{\partial u_-}{\partial x} \right)_{N_x, j}^{k+1, n} - g \zeta_{-, N_x, j}^{k+1, n} + \frac{1}{2}(\lambda - u_0) u_{-, N_x, j}^{k+1, n} = \mu \left( \frac{\partial u_+}{\partial x} \right)_{0, j}^{k, n} - g \zeta_{+, 0, j}^{k, n} + \frac{1}{2}(\lambda - u_0) u_{+, 0, j}^{k, n} $$

We note that the two first left terms are present in the discrete equation for $u$ (47). We isolate then these two terms:

$$ \mu \left( \frac{\partial u_-}{\partial x} \right)_{N_x, j}^{k+1, n} - g \zeta_{-, N_x, j}^{k+1, n} = A_+^k - \frac{1}{2}(\lambda - u_0) u_{-, N_x, j}^{k+1, n} \quad (53) $$

where $A_+^k$ denotes the right-hand side of the first equation of (52), i.e. all the terms depending on $u_+^k$. The main objective in the sequel is to identify $A_-^k$ in the discrete equation for $u$ in $\omega^-$ and then to extract it from the model defined in $\omega^+$.

From the discrete interior equation for $u$ in $\omega^-$, we deduce:

$$ 
\begin{aligned}
& \frac{1}{\Delta t} (u_{-, N_x, j}^{k+1, n+1} - u_{-, N_x, j}^{k+1, n}) + \frac{u_0}{\Delta x} (u_{-, N_x, j}^{k+1, n} - u_{-, N_x-1, j}^{k+1, n}) + \frac{v_0}{2\Delta y} (u_{-, N_x, j+1}^{k+1, n} - u_{-, N_x, j-1}^{k+1, n}) \\
& - \frac{g}{\Delta x} \zeta_{-, N_x-1, j}^{k+1, n} + \frac{\mu}{\Delta x} \left( \frac{\partial u_-}{\partial x} \right)_{N_x-1, j}^{k+1, n} - \mu \frac{u_{-, N_x, j+1}^{k+1, n} - 2u_{-, N_x, j}^{k+1, n} + u_{-, N_x, j-1}^{k+1, n}}{\Delta y^2} - (\tau_x)_{N_x, j} \\
& = \frac{\mu}{\Delta x} \left( \frac{\partial u_-}{\partial x} \right)_{N_x, j}^{k+1, n} - \frac{g}{\Delta x} \zeta_{-, N_x, j}^{k+1, n}
\end{aligned}
\quad (54) 
$$

Therefore, we replace the right-hand side by (53):

$$ 
\begin{aligned}
& \frac{1}{\Delta t} (u_{-, N_x, j}^{k+1, n+1} - u_{-, N_x, j}^{k+1, n}) + \frac{u_0}{\Delta x} (u_{-, N_x, j}^{k+1, n} - u_{-, N_x-1, j}^{k+1, n}) + \frac{v_0}{2\Delta y} (u_{-, N_x, j+1}^{k+1, n} - u_{-, N_x, j-1}^{k+1, n}) \\
& - \frac{g}{\Delta x} \zeta_{-, N_x-1, j}^{k+1, n} + \frac{\mu}{\Delta x} \left( \frac{\partial u_-}{\partial x} \right)_{N_x-1, j}^{k+1, n} - \mu \frac{u_{-, N_x, j+1}^{k+1, n} - 2u_{-, N_x, j}^{k+1, n} + u_{-, N_x, j-1}^{k+1, n}}{\Delta y^2} - (\tau_x)_{N_x, j} \\
& = \frac{A_+^k}{\Delta x} - \frac{1}{2\Delta x} (\lambda - u_0) u_{-, N_x, j}^{k+1, n}
\end{aligned}
\quad (55) 
$$

The next step consists in extracting $A_+^k$ from the model defined in $\omega^+$. We have:

$$ \frac{A_+^k}{\Delta x} = \frac{\mu}{\Delta x} \left( \frac{\partial u_+}{\partial x} \right)_{0, j}^{k, n} - \frac{g}{\Delta x} \zeta_{+, 0, j}^{k, n} + \frac{1}{2\Delta x} (\lambda - u_0) u_{+, 0, j}^{k, n} \quad (56) $$
The equation for $u$ in $\omega^+$ at iteration $k$ of the Schwarz algorithm and at the interface $\gamma$ reads:

$$ \frac{1}{\Delta t}(u_{+,0,j}^{k,n+1} - u_{+,0,j}^{k,n}) + \frac{u_0}{\Delta x}(u_{+,1,j}^{k,n} - u_{+,0,j}^{k,n}) + \frac{v_0}{2\Delta y}(u_{+,0,j+1}^{k,n} - u_{+,0,j-1}^{k,n}) + \frac{g}{\Delta x}(\zeta_{+,1,j}^{k,n} - \zeta_{+,0,j}^{k,n}) - (\tau_x)_{0,j} - \mu \frac{(u_{+,0,j+1}^{k,n} - 2u_{+,0,j}^{k,n} + u_{+,0,j-1}^{k,n})}{\Delta y^2} = \frac{\mu}{\Delta x} \left(\frac{\partial u^k}{\partial x}\right)_{+,1,j}^n - \frac{\mu}{\Delta x} \left(\frac{\partial u^k}{\partial x}\right)_{+,0,j}^n \quad (57) $$

We deduce then that:

$$ \frac{\mu}{\Delta x} \left(\frac{\partial u^k}{\partial x}\right)_{+,0,j}^n - \frac{g}{\Delta x} \zeta_{+,0,j}^{k,n} = -\frac{1}{\Delta t}(u_{+,0,j}^{k,n+1} - u_{+,0,j}^{k,n}) - \frac{u_0}{\Delta x}(u_{+,1,j}^{k,n} - u_{+,0,j}^{k,n}) + \mu \frac{(u_{+,0,j+1}^{k,n} - 2u_{+,0,j}^{k,n} + u_{+,0,j-1}^{k,n})}{\Delta y^2} - \frac{v_0}{2\Delta y}(u_{+,0,j+1}^{k,n} - u_{+,0,j-1}^{k,n}) - \frac{g}{\Delta x} \zeta_{+,1,j}^{k,n} + \frac{\mu}{\Delta x} \left(\frac{\partial u^k}{\partial x}\right)_{+,1,j}^n + (\tau_x)_{0,j} \quad (58) $$

and then we have:

$$ \frac{A_{+}^{k}}{\Delta x} = \frac{1}{2\Delta x}(\lambda - u_0)u_{+,0,j}^{k,n} - \frac{1}{\Delta t}(u_{+,0,j}^{k,n+1} - u_{+,0,j}^{k,n}) - \frac{u_0}{\Delta x}(u_{+,1,j}^{k,n} - u_{+,0,j}^{k,n}) - \frac{v_0}{2\Delta y}(u_{+,0,j+1}^{k,n} - u_{+,0,j-1}^{k,n}) + \mu \frac{(u_{+,0,j+1}^{k,n} - 2u_{+,0,j}^{k,n} + u_{+,0,j-1}^{k,n})}{\Delta y^2} - \frac{g}{\Delta x} \zeta_{+,1,j}^{k,n} + \frac{\mu}{\Delta x} \left(\frac{\partial u^k}{\partial x}\right)_{+,1,j}^n + (\tau_x)_{0,j} \quad (59) $$

Substituting $A_{+}^{k}$ in (55), we obtain the final discrete equation for $u$ in $\omega^{-}$:

$$ \frac{1}{\Delta t}(u_{-,N_x,j}^{k+1,n+1} - u_{-,N_x,j}^{k+1,n}) + \frac{u_0}{\Delta x}(u_{-,N_x,j}^{k+1,n} - u_{-,N_x,j-1}^{k+1,n}) + \frac{v_0}{2\Delta y}(u_{-,N_x,j+1}^{k+1,n} - u_{-,N_x,j-1}^{k+1,n}) - \frac{g}{\Delta x}\zeta_{-,N_x,j}^{k+1,n} + \frac{\mu}{\Delta x}\left(\frac{\partial u_-}{\partial x}\right)_{N_x-1,j}^{k+1,n} - \mu\frac{u_{-,N_x,j+1}^{k+1,n} - 2u_{-,N_x,j}^{k+1,n} + u_{-,N_x,j-1}^{k+1,n}}{\Delta y^2} - (\tau_x)_{N_x,j} = -\frac{1}{2\Delta x}(\lambda - u_0)u_{-,N_x,j}^{k+1,n} + \frac{1}{2\Delta x}(\lambda - u_0)u_{+,0,j}^{k,n} - \frac{1}{\Delta t}(u_{+,0,j}^{k,n+1} - u_{+,0,j}^{k,n}) - \frac{u_0}{\Delta x}(u_{+,1,j}^{k,n} - u_{+,0,j}^{k,n}) - \frac{v_0}{2\Delta y}(u_{+,0,j+1}^{k,n} - u_{+,0,j-1}^{k,n}) + \mu\frac{(u_{+,0,j+1}^{k,n} - 2u_{+,0,j}^{k,n} + u_{+,0,j-1}^{k,n})}{\Delta y^2} \quad (60) $$

Finally, we note that the discretization of the equations in $\omega^+$ and of the boundary condition $B^+$ is the same at each iteration of the algorithm. In fact, the third component of $B^+$ implies a Dirichlet-Dirichlet transmission condition. No special treatment is then needed.

### 5.1.3 Numerical results

We split the channel into two subdomains of similar dimensions, see figure 3, and implement the discrete domain decomposition algorithm described in the previous paragraph. In order to avoid heavy computations by solving the algorithm over a long time interval and because Schwarz algorithm is more efficient on small windows, we decompose the time interval $(0, T)$, where $T = 3.5$ days, into $P$ windows of equal lengths, as in [2] and [16]. Since the aim of this numerical test is to validate the convergence of the decomposition domain algorithm for any $\lambda > 0$, Figure 4a displays the the relative $L^2$ norm, integrated in time over the period, of the difference between the solution of the eastern model and the corresponding reference solution as a function of the Schwarz iteration $k$, for different values of $\lambda$. We can see that the algorithm converges in relatively few iterations. The errors are of the order of $10^{-2}$, which is similar to those obtained in [16] in the case of a non-overlapping algorithm. As
<figure><img src="image_4.png" /><figcaption>Figure 3: Configuration of the splitted domain</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(a) Relative $L^2$ norm of the error as a function of (b) $L^2$ norm of the error at iteration $k = 3$ for the iteration index, for several values of $\lambda$ different values of $\lambda$</figcaption></figure>

<caption>Figure 4</caption>

mentioned in this work, an overlapping of one grid point could accelerate the convergence for a low extra cost. We do not test this possibility in the present work. We also notice that some “extremal" values of $\lambda$ imply a very slow convergence ($\lambda = 0.5$) or even a non-convergence ($\lambda = 8.15$) of the algorithm. This could be explained by the fact that these values make the algorithm behaving approximately as a non-overlapping domain decomposition algorithm with “Dirichlet-Dirichlet" interface conditions or with “Neumann-Neumann" interface conditions. We also show in Figure 4b the $L^2$ norm of the error at the third iteration of the algorithm as a function of $\lambda$. The numerical results suggest that some values of $\lambda$ give a better convergence of the algorithm. One can then numerically optimize the rate of the convergence.

## 5.2 Coupling hydrostatic linearized Navier-Stokes system with corresponding linearized shallow water system

First, let us note that for the sake of simplicity, and since the boundary conditions (29) and (31) do not depend on $y$, we will numerically illustrate the coupling of the one dimensional $x$ version of (17) with the two dimensional $(x, z)$ version of (18). For obvious reasons, we consider a finite domain. We set then in this section $\Omega = [-L; L] \times [-H; 0]$ and $\omega = [-L; L]$.
### 5.2.1 Monodomain reference solution

The numerical reference solution is obtained by solving the x-z version of the system (7) on the whole domain $\Omega$. We do not detail the discretization of the models and of the interface boundary conditions, since it is quite similar to the discretization of the linearized shallow water equations. The physical dimensions of the whole domain are $L = 160m$, and $H = 2m$. Consequently the aspect ratio for the 2D model is equal to 0.0125. The other physical parameters of the simulation are detailed in the following table:

<table><thead><tr><th>Δx</th><th>Δz</th><th>Δt</th><th>μ</th><th>τ₀</th><th>u₀</th></tr></thead><tbody><tr><td>0.1 m</td><td>0.1 m</td><td>0.001 s</td><td>0.8 m²/s⁻¹</td><td>10⁻⁹ N/m²</td><td>0.01 m/s</td></tr></tbody></table>

Due to the explicit Euler time scheme discretization, the time-step is very small. We first perform a preliminary 100 000 time-step simulation, starting with an initial constant horizontal velocity ($u_0 = 0.01m/s$) and with a local surface forcing $\partial_n u(x_i) = -\tau_0 \exp(x_i/L_1)$, in order to create some vertical motion in the right part of the domain ($L_1$ is the length of the zone in which we apply the boundary condition). The final state of this preliminary simulation, displayed in Figure 5, will be used as the initial condition for the following numerical experiments: solution of the discrete version of (7) on the whole domain $\Omega$ (reference 2D solution), and solution of the Schwarz coupling algorithm between the 1D model and the 2D model. In these experiments, the surface friction is reset to zero, in order to fulfill the boundary conditions of the 2D model. As we can see on the figure, the vertical effects of the preliminary friction on a part of the surface are still present in the horizontal velocity field.

### 5.2.2 Numerical results of the coupling algorithm

We split now the domain $\Omega$ into two subdomains of equal size, and in the left subdomain, we replace the x-z version of the linearized hydrostatic Navier-Stokes system (designed in the sequel by the 2D model) by a x version of the linearized shallow water system (designed by the 1D model) — see Figure 6. The discretization and the boundary conditions for both 1D and 2D models are calculated in the same way as in the domain decomposition method.

**Convergence** Figure 7 presents the $L^2$ and the relative $L^2$ norms of the difference between the successive iterates of the Schwarz algorithm. We can see that the algorithm converges for the different tested values of $\lambda$. However the convergence is slow.

**Sensitivity to the parameter $\lambda$** As mentioned and studied in section 2.4, unlike the classical domain decomposition methods, the converged coupled solution is not equal to the restriction of the reference solution on each subdomain. In theorem 1 we highlight a dependence of this modeling error with respect to $\lambda$, to $\varepsilon$ and to the interface position. Figure 8 illustrates this dependence w.r.t. $\lambda$. We can see that, for the tested values of $\lambda$, this error remains of the order of $10^{-2}$ all along the coupling windows. And we observe that the $L^2$ amplitude of these variations actually does not depend on $\lambda$, similarly to the numerical results obtained in [21] in the ellipitic case.

**Sensitivity to the interface position** We also performed experiments to illustrate the dependence of the modeling error on the interface position. We fixed the value of $\lambda$ to 0.2, already used in the previous experiments (Figures 7 and 8). We can see on Figure 9 that the error regularly increases with the shift of the coupling interface. The error begins to increase significantly for an interface located around 75% of the entire 2D domain length, which means that the coupled solution is not able to mimic a real 2D behaviour when the 2D subdomain becomes too small. This numerical result is similar to the one obtained in [22] in the case of a 1D Laplace equation coupled with a 2D Laplace equation.
<figure><img src="image_6.png" /><figcaption>Figure 5: Horizontal velocity at initial and final time of the 2D reference simulation.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Grids of the 1D/2D domains</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Left panel: $L^2$ norm of the difference between successive iterates $\|U^{k+1} - U^k\|$; Right panel: relative $L^2$ norm of the difference between successive iterates $\frac{\|U^{k+1} - U^k\|}{\|U^k\|}$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Evolution of the $L^2$ norm of the modeling error, as a function of time, for some relevant values of $\lambda$</figcaption></figure>

# 6 Conclusion

We presented in this work a Schwarz-like algorithm to couple the 3D linearized hydrostatic Navier-Stokes system with corresponding 2D linearized shallow water system obtained from the 3D equations under a small aspect ratio hypothesis. After introducing the iterative coupling method, we prove that, if we assume a frictionless condition at the bottom, the convergence of the coupling algorithm is equivalent to the convergence of classical domain decomposition method applied to the shallow water system. The main contribution of this work is the study of the role of the coupling interface location. This work can be extended in several directions: numerical optimization of the convergence rate of the domain decomposition algorithm for the linearized shallow water equation, study of the sensitivity of the modeling error to the aspect ratio, or set-up of a method to calculate the optimal interface location. Moreover, if we consider a non-zero friction on the bottom, which is often the case in real applications, there is no more equivalence between the convergence of the coupling algorithm and the convergence of the classical domain decomposition of the shallow water system. Therefore the convergence of the multi-dimensional coupling algorithm with Robin-type conditions, as well as the control of the modeling error, could be more complicated to obtain. This is essentially due to the expression of the extension operator:
<figure><img src="image_10.png" /><figcaption>Figure 9: $L^2$ norm, integrated in time, of the modeling error as a function of the interface position (expressed as a percentage of the total length of the domain)</figcaption></figure>

the quantities coming from the 2D model are no more uniformly extended on the vertical through the interface, but rather follow a parabolic distribution. The study of this case is therefore an interesting alternative. This work is also a first step for more realistic coupling problem of 2D shallow water system, 3D hydrostatic Navier-Stokes system and 3D nonhydrostatic Navier-Stokes system. The most challenging perspective of the present work remains the design of a Schwarz algorithm to couple dimensionally heterogeneous non linear systems.

## Acknowledgements

This work was partially supported by the research department of the French national electricity company, EDF R&D.

## References

[1] A. Arakawa and V. Lamb. Computational design of the basic dynamical processes of the ucla general circulation model. *Methods in Computational Physics*, 17, 1977.

[2] E. Audusse, P. Dreyfuss, and B. Merlet. Schwarz wave form relaxation for Primitive equations of the ocean. *SIAM Journal on Scientific Computing*, 32(5): 2908–2936, 2010.

[3] P.J. Blanco, M. Discacciati, and A. Quarteroni. Modeling dimensionally-heterogeneous problems: analysis, approximation and applications. *Numerische Mathematik*, 119(2): 299–335, 2011.

[4] E. Blayo and L. Debreu. Revisiting open boundary conditions from the point of view of characteristic variables. *Elsevier Ocean Modelling*, 9: 231–252, 2005.

[5] E. Blayo, A. Rousseau, and M. Tayachi-Pigeonnat. Boundary conditions and schwarz waveform relaxation method for linear viscous shallow water equations in hydrodynamics. *SMAI Journal of Computational Mathematics*, 3:117–137, 2017.

[6] A. Decoene. Modèle hydrostatique pour les écoulements à surface libre tridimensionnels et schémas numériques. PhD thesis, Université Paris 6, 2006.

[7] P. Finaud-Guyot, C. Delenne, V. Guinot, and L. Cecile. 1d-2d coupling for river flow modeling. *Comptes Rendus Mecanique - C R MEC*, 339:226-234, 2011.
[8] L. Formaggia, J.-F. Gerbeau, F. Nobile, and A. Quarteroni. On the coupling of 3d and 1d navier-stokes equations for flow problems in compliant vessels. *Computer Methods in Applied Mechanics and Engineering*, 191:561–582, 2001.

[9] M.J. Gander. Optimized Schwarz methods. *SIAM Journal on Numerical Analysis*, 44(2):699–731, 2006.

[10] M.J. Gander. Schwarz methods over the course of time. *Electron. Trans. Numer. Anal.*, 31:228–255, 2008.

[11] M.J. Gander, L. Halpern, and F. Nataf. Optimal convergence for overlapping and non-overlapping schwarz waveform relaxation. In *ddm.org, editor, Eleventh International Conference on Domain Decomposition Methods*, pages 27–36, Augsburg, 1999.

[12] J.-F. Gerbeau and B. Perthame. Derivation of Viscous Saint-Venant System for Laminar Shallow Water; Numerical Validation. *Discrete and Continuous Dynamical Systems*, 1(1): 89–102, 2001.

[13] J. Leiva, P. Blanco, and G. Buscaglia. Partitioned analysis for dimensionally-heterogeneous hydraulic networks. *Multiscale Modeling and Simulation*, 9:872–903, 2011.

[14] N. Malleron, F. Zaoui, N. Goutal, and T. Morel. On the use of a high-performance framework for efficient model coupling in hydroinformatics. *Environmental Modelling and Software - ENVSOFT*, 26, 2011.

[15] J. Marin and J. Monnier. Superposition of local zoom models and simultaneous calibration for 1d–2d shallow water flows. *Mathematics and Computers in Simulation*, 80:547–560, 2009.

[16] V. Martin. Schwarz waveform relaxation algorithms for the linear viscous equatorial shallow water equations. *SIAM J. Scientific Computing*, 31:3595–3625, 2009.

[17] E. Miglio, S. Perotto, and F. Saleri. Model coupling techniques for free-surface flow problems: Part i. *Nonlinear Analysis-theory Methods and Applications - NONLINEAR ANAL-THEOR METH APP*, 63, 2005.

[18] G. Panasenko. Method of Asymptotic Partial Decomposition of Domain. *Mathematical Models and Methods in Applied Sciences*, 8(1): 139–156, 1998.

[19] A. Rousseau. *Etudes théoriques et numériques des équations primitives de l'océan sans viscosité*. PhD thesis, Université Paris XI Orsay, 2005.

[20] A. Rousseau, R. Temam, and J. Tribbia. Boundary value problems for the inviscid primitive equations in limited domain. In Philippe G. Ciarlet, Roger Temam, and Joe Tribbia, editors, *Computational Methods for the Atmosphere and the Oceans*, volume 14 of *Handbook of Numerical Analysis*, pages 481–576. Elsevier, November 2008.

[21] M. Tayachi. *Couplage de modèles de dimensions hétérogènes et application en hydrodynamique*. PhD thesis, University of Grenoble, 2013.

[22] M. Tayachi, Pigeonnat, A. Rousseau, E. Blayo, N. Goutal, and V. Martin. Design and analysis of a schwarz coupling method for a dimensionally heterogeneous problem. *International Journal for Numerical Methods in Fluids*, 75(6):446–465, 2014.