# Barotropic-Baroclinic Splitting for Multilayer Shallow Water Models with Exchanges

Nina Aguillon, Sophie Hörnschemeyer

► To cite this version:

Nina Aguillon, Sophie Hörnschemeyer. Barotropic-Baroclinic Splitting for Multilayer Shallow Water Models with Exchanges. 2026. hal-05483700

HAL Id: hal-05483700

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

Preprint submitted on 29 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Barotropic-Baroclinic Splitting for Multilayer Shallow Water Models with Exchanges

Nina Aguillon¹, Sophie Hörnschemeyer², Jacques Sainte-Marie¹,³

¹ LJLL, Sorbonne Université, Paris, France

² IGPM, RWTH Aachen University, Germany

³ INRIA, Paris, France

January 15, 2026

## Contents

<table><tr><td>1</td><td>Barotropic-baroclinic splitting at the continuous level</td><td>4</td></tr><tr><td>1.1</td><td>Barotropic-baroclinic splitting for the multilayer shallow water model with mass exchanges</td><td>6</td></tr><tr><td>1.2</td><td>Two hyperbolicity results</td><td>8</td></tr><tr><td>2</td><td>Numerical scheme</td><td>10</td></tr><tr><td>2.1</td><td>Overview of the scheme</td><td>10</td></tr><tr><td>2.1.1</td><td>Baroclinic step</td><td>11</td></tr><tr><td>2.1.2</td><td>Barotropic step</td><td>13</td></tr><tr><td>2.2</td><td>CFL conditions and properties of the scheme</td><td>14</td></tr><tr><td>2.2.1</td><td>Properties of the baroclinic step</td><td>14</td></tr><tr><td>2.2.2</td><td>Properties of the barotropic step</td><td>15</td></tr><tr><td>2.2.3</td><td>Properties of the entire scheme</td><td>16</td></tr><tr><td>3</td><td>Numerical experiments</td><td>19</td></tr><tr><td>3.1</td><td>Convergence study: analytical solution for the 2d hydrostatic Euler equations</td><td>19</td></tr><tr><td>3.2</td><td>Computational cost</td><td>21</td></tr><tr><td>3.3</td><td>Wind driven cavity</td><td>21</td></tr><tr><td>3.4</td><td>Well-balancing for the lake at rest</td><td>23</td></tr><tr><td>3.5</td><td>Western boundary current and "linear well-balancing" for the geostrophic equilibrium</td><td>25</td></tr><tr><td>4</td><td>Proofs of energy (in)equalities</td><td>26</td></tr><tr><td>4.1</td><td>Continuous energy equalities for smooth solutions</td><td>26</td></tr><tr><td>4.2</td><td>Discrete entropy inequality</td><td>28</td></tr></table>

## Abstract

This work presents the numerical analysis of a barotropic-baroclinic splitting in a nonlinear multilayer framework with exchanges between the layers in terrain-following coordinates. The splitting is formulated as an exact operator splitting. The barotropic step handles free surface evolution and depth-averaged velocity via a well-balanced one-layer model, while the baroclinic step manages vertical exchanges between layers and adjusts velocities to their mean values. We show that the barotropic-baroclinic splitting preserves total energy conservation and meets both a
discrete maximum principle and a discrete entropy inequality. Several numerical experiments are presented showing the gain in computational cost, particularly in low Froude simulations, and the benefits of using well-balancing strategies for the barotropic step.

**Keywords:** multilayer models, Saint-Venant equations, splitting, free surface flows, well-balancing, density-stratified flows

Introduction

Ocean global circulation models (OGCM) are numerical cores able to simulate the ocean with atmospheric forcing on the whole planet for very long periods of time. They are used in the CMIP simulations (Coupled Model Intercomparison Project) that inform IPPC reports (Inter-governmental Panel on Climate Change). In such simulations the primary focus is on the evolution of stratification and ocean currents through time, and the space and time scales are enormous. The underlying equations are the free surface Navier Stokes equations with hydrostatic pressure and the Boussinesq assumption. Without any specific adjustment the time step is constrained by the fastest wave speed, namely of the surface gravity waves with order of magnitude of 200 ms⁻¹. The fastest current velocity is of about 3 ms⁻¹.

Such a severe constraint is unacceptable and from the very beginning of global ocean circulation simulation efforts have been made to alleviate this constraint. The old rigid lid approximation [Bry69] quickly showed its limits and was replaced either by the implicit surface height method [DSM93] or by the barotropic-baroclinic decomposition [KWSP91] under study in the present work. The main idea is that the surface gravity waves arise from a vertically integrated (barotropic) two-dimensional system resembling to the shallow water model. The other phenomenon reflect the deviation from this barotropic state and correspond to a fully three-dimensional system with much slower behavior (baroclinic). We also refer the reader to [Hig99] and [SM05] for different vertical coordinates and to [DDM⁺¹⁹], [LJW⁺²²] and references therein for more recent developments.

In this paper we revisit the barotropic-baroclinic splitting in the setting of an inviscid terrain- following vertical coordinate (σ-model) and in the constant density case for simplicity. We propose a nonlinear framework where it reads as an operator splitting. We study the hyperbolic properties of both parts of the splitting and illustrate the ability of the splitting in considerably reducing the computational cost without losing accuracy. This is achieved in the framework of multilayer (or layer- averaged) shallow water models ([ABPSM11b], [ABPSM11a], [ABF⁺¹⁹] and references therein). Such models are also used for regional applications for which the barotropic-baroclinic splitting can also have an advantage as soon as the physical phenomenon under scrutiny is not about the free surface waves but rather on a slower material time scale.

Let us briefly recall how multilayer shallow water models are obtained. The starting point are the inviscid Euler equations with a passive tracer T

$$ \left\{ \begin{aligned} & \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0, \\ & \rho_0 \left( \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z} \right) + \frac{\partial p}{\partial x} = 0, \\ & \rho_0 \left( \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z} \right) + \frac{\partial p}{\partial y} = 0, \\ & \rho_0 \left( \frac{\partial w}{\partial t} + u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z} \right) + \frac{\partial p}{\partial z} = -\rho_0 g, \\ & \frac{\partial T}{\partial t} + u \frac{\partial T}{\partial x} + v \frac{\partial T}{\partial y} + w \frac{\partial T}{\partial z} = 0 \end{aligned} \right. \quad (1) \quad \text{eq:Euler} $$

where (u, v, w) is the 3D velocity field of the fluid of constant density ρ₀ and g is the acceleration of gravity. We suppose that at any point (x, y) in the horizontal plane the water fills the vertical interval [z₉, η] where z₉ is the bottom topography and η is the free surface. We now perform a semi discretization of (1) along the vertical. The column of water is divided into N layers with interfaces

$$ z_b = z_{\frac{1}{2}} < z_{\frac{3}{2}} < \cdots < z_{\alpha - \frac{1}{2}} < z_{\alpha + \frac{1}{2}} < z_{N + \frac{1}{2}} = \eta. $$
The α-th layer has total height $h_α = z_{α+1/2} - z_{α+1/2}$, mean horizontal velocity

$$ (u_α, v_α) = \left( \frac{1}{h_α} \int_{z_{α-1/2}}^{z_{α+1/2}} u(t, x, y, z) dz, \frac{1}{h_α} \int_{z_{α-1/2}}^{z_{α+1/2}} v(t, x, y, z) dz \right), $$

and we also denote the mean tracer in the layer $T_α = \frac{1}{h_α} \int_{z_{α-1/2}}^{z_{α+1/2}} T(t, x, y, z) dz$.

The total water height is $h = \sum_{α=1}^N h_α$ and the mean horizontal velocities are given by $\bar{u} = \frac{1}{h} \sum_{α=1}^N h_α u_α$ and $\bar{v} = \frac{1}{h} \sum_{α=1}^N h_α v_α$. Integrating equation (1) for $z \in [z_{α-1/2}, z_{α+1/2}]$ and making the approximations of quadratic terms

$$ \frac{1}{h_α} \int_{z_{α-1/2}}^{z_{α+1/2}} r s dz \approx r_α s_α \quad \text{for all } r, s \in \{u, v, T\} $$

we obtain the system

$$ \left\{ \begin{aligned} & \frac{\partial h_\alpha}{\partial t} + \frac{\partial}{\partial x}(h_\alpha u_\alpha) + \frac{\partial}{\partial y}(h_\alpha v_\alpha) = G_{\alpha+1/2} - G_{\alpha-1/2}, \\ & \frac{\partial}{\partial t}(h_\alpha u_\alpha) + \frac{\partial}{\partial x}(h_\alpha u_\alpha^2) + gh_\alpha \frac{\partial h}{\partial x} + \frac{\partial}{\partial y}(h_\alpha u_\alpha v_\alpha) = -gh_\alpha \frac{\partial z_\alpha}{\partial x} + u_{\alpha+1/2} G_{\alpha+1/2} - u_{\alpha-1/2} G_{\alpha-1/2}, \\ & \frac{\partial}{\partial t}(h_\alpha v_\alpha) + \frac{\partial}{\partial x}(h_\alpha u_\alpha v_\alpha) + \frac{\partial}{\partial x}(h_\alpha v_\alpha^2) + gh_\alpha \frac{\partial h}{\partial y} = -gh_\alpha \frac{\partial z_\alpha}{\partial y} + v_{\alpha+1/2} G_{\alpha+1/2} - v_{\alpha-1/2} G_{\alpha-1/2}, \\ & \frac{\partial}{\partial t}(h_\alpha T_\alpha) + \frac{\partial}{\partial x}(h_\alpha u_\alpha T_\alpha) + \frac{\partial}{\partial y}(h_\alpha v_\alpha T_\alpha) = T_{\alpha+1/2} G_{\alpha+1/2} - T_{\alpha-1/2} G_{\alpha-1/2} \end{aligned} \right. \quad (2) \quad \text{eq:MSW} $$

with mass exchange terms $G_{\alpha \pm 1/2}$, interface velocities $u_{\alpha \pm 1/2}$, $v_{\alpha \pm 1/2}$ and interface tracer concentrations $T_{\alpha \pm 1/2}$ discussed later on. Note that here we considered the simplest possible case by omitting forcings, diffusive terms, and by taking a constant density. A no flux condition at the surface and bottom $G_{1/2} = G_{N+1/2} = 0$ is imposed. The notation are gathered on Figure 1.

The barotropic-baroclinic splitting found in the ocean community corresponds to the separation of two different types of waves: the fast surface gravity waves that arise in a depth-independent shallow water model and the adjustment with respect to this barotropic evolution in a fully baroclinic model with much slower wave speed. Here, we adopt a presentation very close to [LJW+22] in the framework of multilayers models (2), adapting it to the case of a constant density fluid. Firstly the summation over the layers $\alpha = 1, \dots, N$ yields the barotropic system with baroclinic forcing

$$ \left\{ \begin{aligned} & \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h\bar{u}) + \frac{\partial}{\partial y}(h\bar{v}) = 0, \\ & \frac{\partial}{\partial t}(h\bar{u}) + \frac{\partial}{\partial x}(h\bar{u}^2 + \frac{g}{2}h^2) + \frac{\partial}{\partial y}(h\bar{u}\bar{v}) = -gh\frac{\partial z_\alpha}{\partial x} + \frac{\partial}{\partial x}(h\bar{u}^2 - h\bar{u}^2) + \frac{\partial}{\partial y}(h\bar{u}\bar{v} - h\bar{u}\bar{v}), \\ & \frac{\partial}{\partial t}(h\bar{v}) + \frac{\partial}{\partial x}(h\bar{u}\bar{v}) + \frac{\partial}{\partial y}(h\bar{v}^2 + \frac{g}{2}h^2) = -gh\frac{\partial z_\alpha}{\partial y} + \frac{\partial}{\partial x}(h\bar{u}\bar{v} - h\bar{u}\bar{v}) + \frac{\partial}{\partial y}(h\bar{v}^2 - h\bar{v}^2) \end{aligned} \right. \quad (3) \quad \text{eq:oceanoBT} $$

where we denote $\overline{u^2} := \frac{1}{h} \sum_{\alpha=1}^N h_\alpha u_\alpha^2$, $\overline{v^2} := \frac{1}{h} \sum_{\alpha=1}^N h_\alpha v_\alpha^2$ and $\overline{u\bar{v}} := \frac{1}{h} \sum_{\alpha=1}^N h_\alpha u_\alpha v_\alpha$. Secondly some straightforward manipulations yield the baroclinic equation

$$ \left\{ \begin{aligned} & \frac{\partial}{\partial t}(u_\alpha - \bar{u}) + u_\alpha \frac{\partial u_\alpha}{\partial x} + v_\alpha \frac{\partial u_\alpha}{\partial y} + \frac{\bar{u}}{h} \left( \frac{\partial}{\partial x}(h\bar{u}) + \frac{\partial}{\partial y}(h\bar{v}) \right) - \frac{1}{h} \left( \frac{\partial}{\partial x}(h\overline{u^2}) + \frac{\partial}{\partial y}(h\overline{u\bar{v}}) \right) \\ & = \frac{1}{h_\alpha} \left( (u_{\alpha+1/2} - u_\alpha)G_{\alpha+1/2} - (u_{\alpha-1/2} - u_\alpha)G_{\alpha-1/2} \right), \\ & \frac{\partial}{\partial t}(v_\alpha - \bar{v}) + u_\alpha \frac{\partial v_\alpha}{\partial x} + v_\alpha \frac{\partial v_\alpha}{\partial y} + \frac{\bar{v}}{h} \left( \frac{\partial}{\partial x}(h\bar{u}) + \frac{\partial}{\partial y}(h\bar{v}) \right) - \frac{1}{h} \left( \frac{\partial}{\partial x}(h\overline{u\bar{v}}) + \frac{\partial}{\partial y}(h\overline{v^2}) \right) \\ & = \frac{1}{h_\alpha} \left( (v_{\alpha+1/2} - v_\alpha)G_{\alpha+1/2} - (v_{\alpha-1/2} - v_\alpha)G_{\alpha-1/2} \right). \end{aligned} \right. \quad (4) \quad \text{eq:oceanoBC} $$

This latter system describes the evolution of the deviation of the individual velocities from the mean vertical averaged velocities. It is a fully 3D system but no pressure or bathymetry terms appear, thus it can be solved with a large time step denoted by $\Delta t$. In parallel, the barotropic system (3) is used to predict the evolution of $h$ and $(\bar{u}, \bar{v})$ by performing many small time steps $\delta t$ constrained by the velocity of the gravity waves until the large baroclinic time step $\Delta t$ is reached. The left-hand side is frozen at its initial value and (3) reads as the 2D shallow water equations with source term. Thus it remains computationally affordable despite the many time step.
A major issue is that at time $\Delta t$ one has to correct the update of the individual velocities ($u_\alpha, v_\alpha$) obtained in the baroclinic step (4) so that their mean value ($\bar{u}, \bar{v}$) coincides with the one obtained in the barotropic step (3). In order to stabilize the splitting OGCMs also incorporate some time filtering in the barotropic step. The mean velocity ($\bar{u}, \bar{v}$) predicted at time $\Delta t$ is not taken as the final substep, but as an average of many substeps. Two particular filters are described in [DDM+19]: one uses all the substeps in the interval $[0, 2\Delta t]$ and the other the ones between $[\Delta t/2, 3\Delta t/2]$. See also [HB96] for isopycnal models.

In this paper we reformulate the barotropic-baroclinic splitting as an operator splitting in Section 1.1. The evolution in time in (2) is decomposed into two parts:

* a barotropic part gathering the transport at the mean speed ($\bar{u}, \bar{v}$) and the pressure terms;

* a baroclinic part gathering the deviation from the mean velocity and the vertical exchange terms.

The sum of the two contributions exactly gives the source term and space derivatives of (2). Similarly to the original splitting (3)-(4) the barotropic part is essentially depth-independent and contains the surface gravity waves, while the baroclinic part is fully 3D but has much slower dynamics. The main differences are that first, the barotropic step does not have any baroclinic forcing but encompasses a redistribution of the evolution of the mean velocities over the individual velocities and second the mean velocity also evolves in the baroclinic step.

The numerical treatment presented in Section 2 is classically based on a Lie splitting (we only present first order schemes in this article) with a subcycling procedure in the fast barotropic step, see [DL16] and [AB03] for a presentation on such techniques. The numerical treatment of the exchange terms in the baroclinic step is similar to [ABSM18]. We exhibit the time step restrictions of all the parts of the scheme and prove that the total energy decreases at each time step. In other words, the scheme fulfills a discrete entropy inequality, an important property that ensures the convergence towards physical solutions of the scheme (provided that the sequence of approximation does have a limit, Lax theorem).

In Section 3 we present a variety of testcases to illustrate the good properties of the scheme. As expected, the computational cost are reduced drastically for low Froude number simulations. Finally, we take full advantage of proximity between the barotropic system and the simple one-layer shallow water equations to implement state of the art well-balanced schemes for this part of the system. As a result the overall scheme preserves the lake-at-rest equilibrium and the positivity of the water height with almost no additional work (see [ABB+04], [CN17]). We also show that recent work to better maintain the geostrophic equilibrium such as [GTCDGR24] and [ADD+21] provide more accurate results for coarse-mesh long-time integration in the multilayer environment in the presence of the Coriolis force.

# 1 Barotropic-baroclinic splitting at the continuous level

This section is devoted to the presentation and the analysis of a splitting operator of the layer-averaged Euler equations (2) of [ABPSM11b] that shares the same properties as the barotropic-baroclinic splitting of the ocean model (3)-(4) described in the introduction. We present the procedure on the 1D ($x, z$)-model to alleviate the notation, see Figure 1.

$$ \left\{ \begin{aligned} \frac{\partial h_\alpha}{\partial t} + \frac{\partial}{\partial x}(h_\alpha u_\alpha) &= G_{\alpha+\frac{1}{2}} - G_{\alpha-\frac{1}{2}}, \\ \frac{\partial}{\partial t}(h_\alpha u_\alpha) + \frac{\partial}{\partial x}(h_\alpha u_\alpha^2) + gh_\alpha \frac{\partial h}{\partial x} &= -gh_\alpha \frac{\partial z_b}{\partial x} + u_{\alpha+\frac{1}{2}} G_{\alpha+\frac{1}{2}} - u_{\alpha-\frac{1}{2}} G_{\alpha-\frac{1}{2}}, \\ \frac{\partial}{\partial t}(h_\alpha T_\alpha) + \frac{\partial}{\partial x}(h_\alpha u_\alpha T_\alpha) &= T_{\alpha+\frac{1}{2}} G_{\alpha+\frac{1}{2}} - T_{\alpha-\frac{1}{2}} G_{\alpha-\frac{1}{2}}. \end{aligned} \right. \quad (5) \quad \text{eq:MSW1D} $$

The addition of a second horizontal direction $y$ brings no further difficulties.

Before presenting the splitting let us give some more details on model (5). Written in this form the exchange terms are unknown and the system is unclosed with $4N-1$ unknowns ($h_\alpha, u_\alpha, T_\alpha$)$_{1\le\alpha\le N}$ and ($G_{\alpha+\frac{1}{2}}$)$_{1\le\alpha\le N-1}$.

There are two families of multilayer models. One can consider that there are no exchange terms. In that case we obtain a family of shallow water equations stacked on each other and coupled through
the pressure term. This model is ill-posed in the present constant density case, but is useful when the stratification is strong with respect to the velocity shear and corresponds to the so called “isopycnal models” ([Val17, Chapter 3] or [CDV17]). A second possibility is to enforce a constraint on the water heights $h_\alpha$. For example the inner separations between the layers $z_{\alpha+\frac{1}{2}}$, $1 \le \alpha \le N-1$, are horizontal (z-coordinates in ocean models), or all the layers are a given fraction of the total water height (“terrain- following models” or $\sigma$-coordinates).

<figure><img src="image_2.png" /><figcaption>Figure 1: Notations for the multilayer approach.</figcaption></figure>

We focus on the latter family and place our footsteps in the family of multilayer models developed in [ABPSM11b], [ABPSM11a], [ABF$^+$19] and references therein. The notation is gathered in Figure 1. The weights ($l_\alpha$)$_{\alpha \in \{1, \dots, N\}}$ with $\sum_{\alpha=1}^N l_\alpha = 1$ and $l_\alpha > 0$ are fixed once and for all. In such models with mass exchanges, the unknowns in system (5) are the total water height $h$ and the individual velocities ($u_\alpha, v_\alpha$)$_{\alpha \in \{1, \dots, N\}}$. The individual water heights $h_\alpha$ verify the constraint $h_\alpha = l_\alpha h$. The evolution of the total water height is given by the summation over all the layers of the equation on $h_\alpha$ in (5), namely

$$ \frac{\partial h}{\partial t} + \frac{\partial(h\bar{u})}{\partial x} = 0. \qquad (6) \quad \text{eq:h} $$

The exchange terms $G_{\alpha \pm \frac{1}{2}}$ are Lagrange multipliers enforcing the constraints $h_\alpha = l_\alpha h$ at all time. In other words the individual equation on $h_\alpha$ should be understood as a definition of the exchange terms $G_{\alpha \pm \frac{1}{2}}$. Multiplying equation (6) by $l_\alpha$ and subtracting it from the equation on $h_\alpha$ in (5) to cancel the time derivative we obtain

$$ G_{\alpha+\frac{1}{2}} - G_{\alpha-\frac{1}{2}} = l_{\alpha} \frac{\partial}{\partial x} (h(u_{\alpha} - \bar{u})). \qquad (7) \quad \text{eq:G} $$

The final closed system reads

$$ \left\{ \begin{array}{l} \displaystyle \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h\bar{u}) = 0, \\[1.5ex] \displaystyle \frac{\partial}{\partial t}(hu_\alpha) + \frac{\partial}{\partial x}(hu_\alpha^2 + \frac{g}{2}h^2) = -gh\frac{\partial z_b}{\partial x} + \frac{1}{l_\alpha}(u_{\alpha+\frac{1}{2}}G_{\alpha+\frac{1}{2}} - u_{\alpha-\frac{1}{2}}G_{\alpha-\frac{1}{2}}), \\[1.5ex] \displaystyle \frac{\partial}{\partial t}(hT_\alpha) + \frac{\partial}{\partial x}(hu_\alpha T_\alpha) = \frac{1}{l_\alpha}(T_{\alpha+\frac{1}{2}}G_{\alpha+\frac{1}{2}} - T_{\alpha-\frac{1}{2}}G_{\alpha-\frac{1}{2}}), \\[1.5ex] \displaystyle G_{\alpha+\frac{1}{2}} = \sum_{j=1}^\alpha l_j \frac{\partial}{\partial x}(h(u_j - \bar{u})). \end{array} \right. \qquad (8) \quad \text{eq:MSW1D-closed} $$

The interface velocities $u_{\alpha+\frac{1}{2}}$ and interface tracers $T_{\alpha+\frac{1}{2}}$ are defined using an upwinding strategy

$$ u_{\alpha+\frac{1}{2}} = \begin{cases} u_{\alpha}, & \text{if } G_{\alpha+\frac{1}{2}} \le 0 \\ u_{\alpha+1}, & \text{if } G_{\alpha+\frac{1}{2}} > 0 \end{cases} \quad \text{and} \quad T_{\alpha+\frac{1}{2}} = \begin{cases} T_{\alpha}, & \text{if } G_{\alpha+\frac{1}{2}} \le 0 \\ T_{\alpha+1}, & \text{if } G_{\alpha+\frac{1}{2}} > 0 \end{cases}. \qquad (9) \quad \text{eq:interu} $$

We now propose a barotropic-baroclinic splitting for the set of equations (8).
## 1.1 Barotropic-baroclinic splitting for the multilayer shallow water model with mass exchanges

### (S:multilayer) Definition 1.1. The barotropic step corresponds to the advection of h, $u_\alpha$ and $T_\alpha$ at the mean velocity $\bar{u}$ together with the pressure force and the topography source term. It writes

$$ \left\{ \begin{aligned} & \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h\bar{u}) = 0, \\ & \frac{\partial}{\partial t}(hu_\alpha) + \frac{\partial}{\partial x}(hu_\alpha\bar{u} + \frac{g}{2}h^2) = -gh\frac{\partial z_b}{\partial x}, \\ & \frac{\partial}{\partial t}(hT_\alpha) + \frac{\partial}{\partial x}(h\bar{u}T_\alpha) = 0. \end{aligned} \right. \qquad (10) \text{eq:BT-closed} $$

The baroclinic step gathers the nonlinear adjustment around the mean vertical velocity $\bar{u}$ and the exchange terms to maintain the constraint $h_\alpha = l_\alpha h$. It writes

$$ \left\{ \begin{aligned} & \frac{\partial h}{\partial t} = 0, \\ & \frac{\partial}{\partial t}(hu_\alpha) + \frac{\partial}{\partial x}(hu_\alpha(u_\alpha - \bar{u})) = \frac{1}{l_\alpha}(u_{\alpha+\frac{1}{2}}G_{\alpha+\frac{1}{2}} - u_{\alpha-\frac{1}{2}}G_{\alpha-\frac{1}{2}}), \\ & \frac{\partial}{\partial t}(hT_\alpha) + \frac{\partial}{\partial x}(hT_\alpha(u_\alpha - \bar{u})) = \frac{1}{l_\alpha}(T_{\alpha+\frac{1}{2}}G_{\alpha+\frac{1}{2}} - T_{\alpha-\frac{1}{2}}G_{\alpha-\frac{1}{2}}), \\ & G_{\alpha+\frac{1}{2}} = \sum_{j=1}^\alpha l_j \frac{\partial}{\partial x}(h(u_j - \bar{u})) \end{aligned} \right. \qquad (11) \text{eq:BC} $$

where the definition of the interface velocities is given by (9).

An important feature of this splitting is that it is an operator splitting.

**Proposition 1.2.** Write the multilayer operator (8) as

$$ \partial_t \begin{pmatrix} h \\ hu_\alpha \\ hT_\alpha \end{pmatrix} = \text{ML} \begin{pmatrix} h \\ hu_\alpha \\ hT_\alpha \end{pmatrix} $$

and similarly the barotropic operator (10) and the baroclinic operator (11) respectively as

$$ \partial_t \begin{pmatrix} h \\ hu_\alpha \\ hT_\alpha \end{pmatrix} = \text{BT} \begin{pmatrix} h \\ hu_\alpha \\ hT_\alpha \end{pmatrix} \quad \text{and} \quad \partial_t \begin{pmatrix} h \\ hu_\alpha \\ hT_\alpha \end{pmatrix} = \text{BC} \begin{pmatrix} h \\ hu_\alpha \\ hT_\alpha \end{pmatrix}. $$

Then $\text{ML} = \text{BT} + \text{BC}$.

*Proof.* In the original model (8) the first term $u_\alpha$ in the equation on $hu_\alpha$ is written as $u_\alpha = \bar{u} + (u_\alpha - \bar{u})$. Then the mean term in $\bar{u}$ yields the barotropic part (10) while the deviation term $u_\alpha - \bar{u}$ yields the baroclinic part (11). $\square$

We now emphasize the proximity of the barotropic step (10) with the one layer shallow water model.

**Proposition 1.3.** System (10) in variables $(h, hu_\alpha, hT_\alpha)_{\alpha \in \{1, \dots, N\}}$ is equivalent to the system in variables $(h, h\bar{u}, (h\sigma_\alpha)_{\alpha \in \{1, \dots, N-1\}}, (hT_\alpha)_{\alpha \in \{1, \dots, N\}})$ with $\sigma_\alpha = u_\alpha - \bar{u}$:

$$ \left\{ \begin{aligned} & \frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h\bar{u}) = 0, \\ & \frac{\partial}{\partial t}(h\bar{u}) + \frac{\partial}{\partial x}(h\bar{u}^2 + \frac{g}{2}h^2) = -gh\frac{\partial z_b}{\partial x}, \\ & \frac{\partial}{\partial t}(h\sigma_\alpha) + \frac{\partial}{\partial x}(h\bar{u}\sigma_\alpha) = 0, \\ & \frac{\partial}{\partial t}(hT_\alpha) + \frac{\partial}{\partial x}(h\bar{u}T_\alpha) = 0. \end{aligned} \right. \qquad (12) \text{eq:BT} $$

The evolution of $h\sigma_N$ is similar as in the other layers, but this quantity is deduced from the others as $\sum_{\alpha=1}^N l_\alpha h\sigma_\alpha = 0$.

*Proof.* The equation on $h\bar{u}$ is obtained by summing the equations on $hu_\alpha$ weighted by $l_\alpha$. The subtraction of the individual equation on $hu_\alpha$ with the new equation on $h\bar{u}$ gives the equation on $h\sigma_\alpha$. $\square$
**Remark 1.4.** Another possibility is to treat the tracer in a similar fashion, separating the evolution of its mean value $\bar{T} = \sum_{\alpha=1}^{N} l_{\alpha} T_{\alpha}$ and of the deviations $\theta_{\alpha} = T_{\alpha} - \bar{T}$:

$$ 
\left\{
\begin{aligned}
\frac{\partial h}{\partial t} + \frac{\partial}{\partial x}(h\bar{u}) &= 0, \\
\frac{\partial}{\partial t}(h\bar{u}) + \frac{\partial}{\partial x}(h\bar{u}^2 + \frac{g}{2}h^2) &= -gh\frac{\partial z_b}{\partial x}, \\
\frac{\partial}{\partial t}(h\bar{T}) + \frac{\partial}{\partial x}(h\bar{u}\bar{T}) &= 0, \\
\frac{\partial}{\partial t}(h\sigma_{\alpha}) + \frac{\partial}{\partial x}(h\bar{u}\sigma_{\alpha}) &= 0, \\
\frac{\partial}{\partial t}(h\theta_{\alpha}) + \frac{\partial}{\partial x}(h\bar{u}\theta_{\alpha}) &= 0.
\end{aligned}
\right.
$$

This form could be useful if $T$ was an active tracer, meaning that the pressure would depend on the values of $(T_{\alpha})_{\alpha \in \{1, \dots, N\}}$. The extension to this case is left for future work.

Let us take a closer look on the individual equation on $h_{\alpha}$.

**Remark 1.5.** For the barotropic system (10) the individual equation on $h_{\alpha}$ is given by

$$ \frac{\partial h_{\alpha}}{\partial t} + \frac{\partial}{\partial x}(h_{\alpha}\bar{u}) = 0. $$

Together with the equation on it holds

$$ \frac{\partial}{\partial t}(h_{\alpha} - l_{\alpha}h) + \frac{\partial}{\partial x}((h_{\alpha} - l_{\alpha}h)\bar{u}) = 0 $$

and it is clear that if the constraint $h_{\alpha} = l_{\alpha}h$ initially holds this property is preserved through time. This explains why there are no exchange terms in this step.

For the baroclinic step the definition of the exchange terms yields $G_{\alpha+\frac{1}{2}} - G_{\alpha-\frac{1}{2}} = l_{\alpha} \frac{\partial}{\partial x}(h(u_{\alpha} - \bar{u}))$ and

$$ \frac{\partial h_{\alpha}}{\partial t} + \frac{\partial}{\partial x}(h_{\alpha}(u_{\alpha} - \bar{u})) = G_{\alpha+\frac{1}{2}} - G_{\alpha-\frac{1}{2}}. $$

Once again the exchange terms can be seen as Lagrange multipliers enforcing the constraints $h_{\alpha} = l_{\alpha}h$.

We now turn to the classical question of energy conservation or dissipation for both parts of the splitting. In the unsplit multilayer system (8)-(9) the total energy $E = \frac{gh^2}{2} + ghz_b + \sum_{\alpha=1}^{N} \frac{h_{\alpha}u_{\alpha}^2}{2}$ verifies the energy inequality [ABPSM11b]

$$ \frac{\partial}{\partial t}E + \frac{\partial}{\partial x} \left( gh^2 \bar{u} + ghz_b \bar{u} + \sum_{\alpha=1}^{N} \frac{h_{\alpha} u_{\alpha}^3}{2} \right) \leq 0. \qquad (13) \text{ [eq:entropyMSW]} $$

**Prop:energyequality** Proposition 1.6. The barotropic-baroclinic splitting (10)-(11) is associated with a splitting of (13): the smooth solutions of the barotropic step (10) verify the energy equality

$$ \frac{\partial}{\partial t}E + \frac{\partial}{\partial x} \left( gh^2 \bar{u} + ghz_b \bar{u} + \sum_{\alpha=1}^{N} \frac{h_{\alpha} u_{\alpha}^2 \bar{u}}{2} \right) = 0, \qquad (14) \text{ [eq:EIBT]} $$

while smooth solutions of the baroclinic step (11) with interface velocities (9) verify the energy inequality

$$ \frac{\partial}{\partial t}E + \frac{\partial}{\partial x} \left( \sum_{\alpha=1}^{N} \frac{h_{\alpha} u_{\alpha}^2 (u_{\alpha} - \bar{u})}{2} \right) \leq 0. \qquad (15) \text{ [eq:EIBC]} $$

**? (prop:entropy)**

Proof. The proof can be found in Section 4.1. □

<figure><img src="image_3.png" /></figure>
## 1.2 Two hyperbolicity results

In this section we investigate the hyperbolicity of the systems (10) and (11). We ignore the passive tracer T and focus on the nonlinear part.

**Proposition 1.7.** The barotropic system (10) with flat topography $z_b = 0$ is hyperbolic with eigenvalues $\bar{u} - \sqrt{gh}$, $\bar{u}$ with $N-1$ multiplicity and $\bar{u} + \sqrt{gh}$.

*Proof.* We study the hyperbolicity in the variables $(h, \bar{u}, \sigma_1, \dots, \sigma_{N-1})$, see (12). The quasilinear form of this system is

$$ \partial_t \begin{pmatrix} h \\ \bar{u} \\ \sigma_1 \\ \sigma_2 \\ \vdots \\ \sigma_{N-1} \end{pmatrix} + \begin{pmatrix} \bar{u} & h & 0 & \dots & \dots & 0 \\ g & \bar{u} & 0 & \dots & \dots & 0 \\ 0 & 0 & \bar{u} & 0 & \dots & 0 \\ \vdots & \vdots & 0 & \ddots & \ddots & \vdots \\ \vdots & \vdots & \vdots & \ddots & \ddots & 0 \\ 0 & 0 & 0 & \dots & 0 & \bar{u} \end{pmatrix} \partial_x \begin{pmatrix} h \\ \bar{u} \\ \sigma_1 \\ \sigma_2 \\ \vdots \\ \sigma_{N-1} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ \vdots \\ 0 \end{pmatrix} $$

and the result immediately follows.

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

The study of the hyperbolicity of the baroclinic system (11) shares many similarities with the one of the original system (5). It is known from the seminal paper [ABPSM11b] that the two-layer shallow water system is hyperbolic for any choice of interface velocity $u_1$ lying between $u_1$ and $u_2$. For the general N-layer system (5) with upwind interface velocities (9) numerical investigations indicate that there are always two extremal eigenvalues resembling $\bar{u} \pm \sqrt{gh}$ (with a perturbation due to the shear). However "in the middle" complex eigenvalues arise in general. The same behavior arises when considering the hyperbolicity of the free-surface Euler equations (or Benney model) [DMEHG+25].

The following proposition concerns the hyperbolicity of the baroclinic part without mass exchange. At first glance this system seems somehow artificial but:

* it is very close to the multilayer semi-Lagrangian model (5.2) in [DMEHG+25],

* this system is solved in the numerical approximation of the baroclinic step.

**Proposition 1.8.** Consider the following system with the 2N unknowns $(h_\alpha, u_\alpha)_{\alpha \in \{1, \dots, N\}}$

$$ \left\{ \begin{aligned} & \frac{\partial h_\alpha}{\partial t} + \frac{\partial}{\partial x}(h_\alpha(u_\alpha - \bar{u})) = 0, \\ & \frac{\partial}{\partial t}(h_\alpha u_\alpha) + \frac{\partial}{\partial x}(h_\alpha u_\alpha (u_\alpha - \bar{u})) = 0. \end{aligned} \right. \qquad (16) \quad \text{eq:BCnoME} $$

The total water height $h$ and the mean velocity $\bar{u}$ are defined using the other unknowns with $h = \sum_{\alpha=1}^N h_\alpha$ and $\bar{u} = \frac{1}{h} \sum_{\alpha=1}^N h_\alpha u_\alpha$.

Suppose that all the velocity deviations $u_\alpha - \bar{u}$ are different and that none of them vanishes. Suppose moreover that $\sum_{\alpha=1}^N \frac{h_\alpha}{u_\alpha - \bar{u}} \neq 0$. Then system (16) is strictly hyperbolic with the 2N distinct eigenvalues

$$ \{u_1 - \bar{u}, \dots, u_N - \bar{u}\} \cup \left\{ \lambda \in \mathbb{R}, h = \sum_{\alpha=1}^{N} \frac{(u_{\alpha} - \bar{u})h_{\alpha}}{u_{\alpha} - \bar{u} - \lambda} \right\}. $$

Moreover, all the eigenvalues belong to the interval $[\min_\alpha u_\alpha - \bar{u}, \max_\alpha u_\alpha - \bar{u}]$.

*Proof.* Combining the two equations of (16), for regular solutions we obtain the transport equation on $u_\alpha$, namely

$$ \frac{\partial u_{\alpha}}{\partial t} + (u_{\alpha} - \bar{u}) \frac{\partial u_{\alpha}}{\partial x} = 0. \qquad (17) \{?\} $$
We now work on the equation on $h_\alpha$. With the chain rule and $\bar{u} = \frac{\sum_{k=1}^N h_k u_k}{\sum_{k=1}^N h_k}$ we get

$$ 
\begin{aligned}
\frac{\partial}{\partial x}(h_\alpha(u_\alpha - \bar{u})) &= (u_\alpha - \bar{u})\frac{\partial h_\alpha}{\partial x} + h_\alpha \frac{\partial u_\alpha}{\partial x} - h_\alpha \frac{\partial}{\partial x} \left( \frac{\sum_{k=1}^N h_k u_k}{\sum_{k=1}^N h_k} \right) \\
&= (u_\alpha - \bar{u})\frac{\partial h_\alpha}{\partial x} + h_\alpha \frac{\partial u_\alpha}{\partial x} - \frac{h_\alpha}{h} \left( \sum_{k=1}^N (u_k - \bar{u})\frac{\partial h_k}{\partial x} \right) - \frac{h_\alpha}{h} \left( \sum_{k=1}^N h_k \frac{\partial u_k}{\partial x} \right).
\end{aligned}
$$

We define $H = (h_1, \dots, h_N)^T$, $U = (u_1, \dots, u_N)^T$, $\sigma_\alpha = u_\alpha - \bar{u}$, $\Sigma = (\sigma_1, \dots, \sigma_N)^T$, $D_\Sigma = \text{diag}(\Sigma)$ and $L = (\ell_1, \dots, \ell_N)^T = \frac{1}{h}(h_1, \dots, h_N)^T$. Then the quasilinear linear form of (16) is

$$ \frac{\partial}{\partial t} \left( \frac{H}{U} \right) + \left( \frac{D_{\Sigma} - L\Sigma^T}{0} \middle| \frac{hLL^T}{D_{\Sigma}} \right) \frac{\partial}{\partial x} \left( \frac{H}{U} \right) = \frac{\partial}{\partial x} \left( \frac{0}{0} \right). $$

All the $\sigma_\alpha$ are eigenvalues of the matrix. We look for the values of $\lambda \in \mathbb{R}$ such that $D_\sigma - L\Sigma^T - \lambda I$ is non invertible. As long as $\lambda \neq \sigma_k$ the matrix determinant lemma yields

$$ \det(D_{\sigma} - \lambda I - L\Sigma^{T}) = (1 - \Sigma^{T}(D_{\sigma} - \lambda I)^{-1}L) \det(D_{\sigma} - \lambda I). $$

We now prove that the following function has $N$ roots, denoted by $\lambda_1 < \dots < \lambda_N$

$$ \varphi(\lambda) = 1 - \Sigma^T (D_\sigma - \lambda I)^{-1} = 1 - \sum_{\alpha=1}^{N} \frac{\sigma_\alpha \ell_\alpha}{\sigma_\alpha - \lambda}. $$

The fact that $\sum_{\alpha=1}^N h_\alpha \sigma_\alpha = 0$ implies the existence of at least one non-positive and one non-positive $\sigma_\alpha$. We order the velocity deviations: for some $k \in \{1, \dots, N-1\}$ and some permutation $(p_1, \dots, p_N)$ of $\{1, \dots, N\}$,

$$ \sigma_{p_1} < \dots < \sigma_{p_k} < 0 < \sigma_{p_{k+1}} < \dots < \sigma_{p_N}. $$

The limits

$$ \lim_{\lambda \uparrow \sigma_{p_i}} \varphi(\lambda) = \begin{cases} +\infty & \text{if } i \le k \\ -\infty & \text{if } i \ge k+1 \end{cases}, \quad \lim_{\lambda \downarrow \sigma_{p_i}} \varphi(\lambda) = \begin{cases} -\infty & \text{if } i \le k \\ +\infty & \text{if } i \ge k+1 \end{cases} \quad \text{and} \quad \lim_{\lambda \to \pm\infty} \varphi(\lambda) = 1 $$

imply the existence of a root in between each $N-2$ consecutive pairs $(\sigma_{p_1}, \sigma_{p_2}), \dots, (\sigma_{p_{k-1}}, \sigma_{p_k})$, $(\sigma_{p_{k+1}}, \sigma_{p_{k+2}}), \dots, (\sigma_{p_{N-1}}, \sigma_{p_N})$ sharing the same sign. On the interval $[\sigma_{p_k}, \sigma_{p_{k+1}}]$ the function $\varphi$ has two roots. Indeed its limits on the boundaries are $-\infty$ while at $\lambda = 0$ the function $\varphi$ vanishes with derivative $\varphi'(0) = \sum_{\alpha=1}^N \frac{\ell_\alpha}{\sigma_\alpha} = \frac{1}{h} \sum_{\alpha=1}^N \frac{h_\alpha}{u_\alpha - \bar{u}} \neq 0$. Alltogether the system has $2N$ distinct eigenvalues

$$ \sigma_{p_1} < \lambda_1 < \sigma_{p_2} < \dots < \sigma_{p_k} < \lambda_k < \lambda_{k+1} < \sigma_{p_{k+1}} < \lambda_{k+2} < \sigma_{p_{k+2}} < \dots < \lambda_N < \sigma_{p_N}. $$

□

**Remark 1.9.** If $u_\alpha - \bar{u} = 0$ for some $\alpha$ or if $\sum_{\alpha=1}^N \frac{h_\alpha}{u_\alpha - \bar{u}} = 0$ the system is no longer strictly hyperbolic in general. The most obvious examples are the case $u_\alpha - \bar{u} = 0$ for all $\alpha$ with the apparition of the Jordan block $hLL^T$. In the case $N=2$, $u_1 = -u$, $u_2 = u$, $h_1 = h_2$ the matrix reduces to

$$ \begin{pmatrix} -\frac{u}{2} & -\frac{u}{2} & \frac{h}{4} & \frac{h}{4} \\ \frac{u}{2} & \frac{u}{2} & \frac{h}{4} & \frac{h}{4} \\ 0 & 0 & -u & u \\ 0 & 0 & 0 & u \end{pmatrix} $$

and 0 is a double eigenvalue with only one eigenvector.
## 2 Numerical scheme

### (S:scheme)

We propose a numerical discretization of the multilayer system (5) based on the barotropic-baroclinic splitting (11)-(12) presented above. We use a finite volume framework and discretize the computational domain by a grid with equidistant distributed nodes $x_j$. We denote by $C_j = [x_{j-\frac{1}{2}}, x_{j+\frac{1}{2}}]$ the cell of length $\Delta x_j = x_{j+\frac{1}{2}} - x_{j-\frac{1}{2}}$ with $x_{j+\frac{1}{2}} = \frac{x_j+x_{j+1}}{2}$. At the beginning of the time step the quantities ($h_j^n, u_{\alpha,j}^n, T_{\alpha,j}^n$) at time $t^n$ are known. They are updated at time $t^{n+1} = t^n + \Delta t^n$. The baroclinic (11) and the barotropic (12) systems are discretized one after the other with a Lie splitting. More precisely,

1. **Baroclinic step:** System (11) is solved on a large time step $\Delta t^n$. We denote the variables obtained at the end of this step with the superscript $n + \frac{1}{2}$.

2. **Barotropic step:** System (12) is then solved on the time interval $[0, \Delta t^n]$ with $K$ small time substeps $\delta t^k$ using $h_j^{n+\frac{1}{2}}, u_{\alpha,j}^{n+\frac{1}{2}}$ and $T_{\alpha,j}^{n+\frac{1}{2}}$ as initial data. The last sub timestep gives the values $h_j^{n+1}, u_{\alpha,j}^{n+1}$ and $T_{\alpha,j}^{n+1}$. It serves as the initial data for step 1 (the baroclinic step) at the next time step.

<figure><img src="image_5.png" /><figcaption>Figure 2: Overview of the scheme. The large baroclinic step with $\Delta t^n$ that consists of a prediction step where system (16) is solved and a correction step where the mass exchange terms are defined and applied. The computed data of the baroclinic step ($h_\alpha, u_\alpha, T_\alpha$)$^{n+\frac{1}{2}}$ is used as initial data for the barotropic step. We perform many small barotropic time steps $\delta t^k$ until we reach the final time $\Delta t^n$. During this loop we only update the vertically averaged quantities $h$ and $\bar{u}$. The transported variables $\sigma_\alpha$ and $T_\alpha$ are updated only once with a large time step $\Delta t^n$.</figcaption></figure>

We detail these two steps in Section 2.1. The timestep restrictions on the timesteps are discussed in Section 2.2. We prove that both steps ensure a discrete maximum principle on the velocity and on the tracer and a discrete entropy inequality. Special choices of the numerical fluxes are discussed in Section 2.2.3.

### 2.1 Overview of the scheme

### (S:overview)

An overview of the scheme is sketched in Figure 2 and in Algorithm 1. The details of the baroclinic and the barotropic steps are given in the following subsections.
**Algorithm 1** Barotropic-baroclinic splitting. Lines 4 to 6 correspond to the baroclinic part (11). Lines 8 to 13 correspond to the barotropic part (12)

**Input:** Initial data $(h^0, u_\alpha^0, T_\alpha^0)$, final time $t^{\text{final}}$

**Output:** $(h^{\text{final}}, u_\alpha^{\text{final}}, T_\alpha^{\text{final}})$

1: $t \leftarrow 0$

2: $(h, u_\alpha, T_\alpha) \leftarrow (h^0, u_\alpha^0, T_\alpha^0)$

3: **while** $t < t^{\text{final}}$ **do**

4: $\Delta t \leftarrow \min(\text{computeDT\_baroclinic}(h, u_\alpha), t^{\text{final}} - t)$ ▷ see CFL conditions (19), (31), (32)

5: $(h, u_\alpha, T_\alpha) \leftarrow \text{baroclinic}(h, u_\alpha, T_\alpha, \Delta t)$ ▷ see Section 2.1.1

6: $\bar{u} \leftarrow \sum_{\alpha=1}^{N} l_\alpha u_\alpha$

7: $t^{\text{BT}} \leftarrow 0$

8: **while** $t^{\text{BT}} < \Delta t$ **do**

9: $\delta t \leftarrow \min(\text{computeDT\_barotropic}(h, \bar{u}), \Delta t - t^{\text{BT}})$ ▷ see CFL condition (27)

10: $(h, \bar{u}) \leftarrow \text{barotropic\_shallow\_water}(h, \bar{u}, \delta t)$ ▷ see Section 2.1.2

11: $t^{\text{BT}} \leftarrow t^{\text{BT}} + \delta t$

12: **end while**

13: $(\sigma_\alpha, T_\alpha) \leftarrow \text{barotropic\_adjustment}(h, \bar{u}, \sigma_\alpha, T_\alpha, \Delta t)$ ▷ see update (30)

14: $u_\alpha \leftarrow \bar{u} + \sigma_\alpha$

15: $t \leftarrow t + \Delta t$

16: **end while**

17: return $(h, u_\alpha, T_\alpha)$

The notation and the most important steps of the scheme are summarised in Figure 2. Compared to the original scheme without splitting, in which all variables are updated with a small time step $\delta t^n$, there are two advantages in terms of computational cost. Firstly, the exchange terms are only calculated in each large time step $\Delta t^n$ in the baroclinic step. Secondly, only the water height $h$ and the mean velocity $\bar{u}$ are calculated with small time steps. The transported variables $\sigma_\alpha$ and $T_\alpha$, on the other hand, are updated with only one large time step $\Delta t^n$.

### 2.1.1 Baroclinic step

This step corresponds to (11). At the beginning of the timestep the values of $h_j^n, u_{\alpha,j}^n$ and $T_{\alpha,j}^n$ are known. We follow the prediction/correction strategy of [ABSM18]. An overview of the baroclinic step is given in Algorithm 2.

**Algorithm 2 Baroclinic step**

**Input:** $h^n, u_\alpha^n, T_\alpha^n, \Delta t^n$

**Output:** $h^{n+\frac{1}{2}}, u_\alpha^{n+\frac{1}{2}}, T_\alpha^{n+\frac{1}{2}}$

1: $(h^{n+\frac{1}{2}}, h_\alpha^*, u_\alpha^*, T_\alpha^*) \leftarrow \text{prediction}(h^n, u_\alpha^n, T_\alpha^n, \Delta t^n)$ ▷ see (18), (21)

2: $G_{\alpha \pm \frac{1}{2}}^* \leftarrow \text{exchanges}(h^{n+\frac{1}{2}}, h_\alpha^*)$ ▷ see (22)

3: $(u_\alpha^{n+\frac{1}{2}}, T_\alpha^{n+\frac{1}{2}}) \leftarrow \text{correction}(h^*, u_\alpha^*, T_\alpha^*, G_{\alpha \pm \frac{1}{2}}^*, \Delta t^n)$ ▷ see (24)

4: return $h^{n+\frac{1}{2}}, u_\alpha^{n+\frac{1}{2}}, T_\alpha^{n+\frac{1}{2}}$

### Prediction of the total water height

The prediction step updates the total water height $h_j^{n+\frac{1}{2}}$ and defines the exchange terms in a way that the constraint $h_\alpha = l_\alpha h$ is preserved through time. It consists in solving system (16) without mass exchanges with $2N$ unknowns studied in Proposition 1.8. We recall that in the absence of exchange terms each individual water height $h_\alpha$ has its own evolution. The numerical fluxes are denoted with
the letter $F$. At the end of this step we obtain

$$ egin{cases} h_{\alpha,j}^* = h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} (F_{j+\frac{1}{2}}^{h_\alpha} - F_{j-\frac{1}{2}}^{h_\alpha}), \\ (h_\alpha u_\alpha)_j^* = (h_\alpha u_\alpha)_j^n - \frac{\Delta t^n}{\Delta x} (F_{j+\frac{1}{2}}^{h_\alpha u_\alpha} - F_{j-\frac{1}{2}}^{h_\alpha u_\alpha}), \\ (h_\alpha T_\alpha)_j^* = (h_\alpha T_\alpha)_j^n - \frac{\Delta t^n}{\Delta x} (F_{j+\frac{1}{2}}^{h_\alpha T_\alpha} - F_{j-\frac{1}{2}}^{h_\alpha T_\alpha}) \end{cases} \quad (18) \quad \text{eq:predictionstep} $$

with $\Delta t^n$ being constrained by the CFL condition motivated by Proposition 1.8 about the eigenvalues of the baroclinic prediction step

$$ \frac{\Delta t^n}{\Delta x} \max_{(\alpha,j)} |u_{\alpha,j}^n - \bar{u}_j^n| \le \text{CFL} \le \frac{1}{2} \quad (19) \quad \text{eq:CFL

for a constant CFL-number.

Two specific choices of numerical mass flux $F_{j+\frac{1}{2}}^{h_\alpha}$ are discussed later on in Propositions 2.10 and 2.11 to ensure several desirable properties of the overall scheme. On the other hand, we remark that $u_\alpha$ and $T_\alpha$ are transported at speed $\sigma_\alpha = u_\alpha - \bar{u}$ and choose the numerical fluxes [AB03]

$$ F_{j+\frac{1}{2}}^{h_\alpha \varphi_\alpha} = \varphi_{\alpha,j}^n (F_{j+\frac{1}{2}}^{h_\alpha})^+ - \varphi_{\alpha,j+1}^n (F_{j+\frac{1}{2}}^{h_\alpha})^- \quad \text{for } \varphi_\alpha \in \{u_\alpha, T_\alpha\}, \quad (20) \quad \text{eq:fluxdec} $$

where we denote for any real number $a$, its positive part by $a^+ := \max(a, 0)$ and its negative part by $a^- := -\min(a, 0)$, thus $a = a^- - a^+$. The total water height is defined as

$$ h_j^{n+\frac{1}{2}} = \sum_{\alpha=1}^{N} h_{\alpha,j}^* = h_j^n - \frac{\Delta t^n}{\Delta x} (F_{j+\frac{1}{2}}^h - F_{j-\frac{1}{2}}^h) \quad \text{where} \quad F_{j\pm\frac{1}{2}}^h = \sum_{\alpha=1}^{N} F_{j\pm\frac{1}{2}}^{h_\alpha}. \quad (21) \quad \text{eq:hjd} $$

## Definition of the exchange terms

The individual water heights $h_{\alpha,j}^*$ obtained in the prediction step (18) are not proportional to each other: $h_{\alpha,j}^* \ne l_\alpha h_j^{n+\frac{1}{2}}$. The exchange terms are defined to restore this constraint.

### Proposition 2.1. Define the exchange terms as $G_{\frac{1}{2},j}^* := 0$ and for all $\alpha \in \{1, \dots, N-1\}$

$$ G_{\alpha+\frac{1}{2},j}^* := \frac{1}{\Delta t^n} \left( l_\alpha h_j^{n+\frac{1}{2}} - h_{\alpha,j}^* \right) + G_{\alpha-\frac{1}{2},j}^* \quad (22) \quad \text{eq:def_exchange_t} $$

Then $G_{N+\frac{1}{2},j}^* := 0$ and for any $\alpha \in \{1, \dots, N\}$ it holds

$$ h_{\alpha,j}^{n+\frac{1}{2}} = h_{\alpha,j}^* + \Delta t^n (G_{\alpha+\frac{1}{2},j}^* - G_{\alpha-\frac{1}{2},j}^*) = l_\alpha h_j^{n+\frac{1}{2}}. \quad (23) \quad \text{eq:hajd} $$

*Proof.* Equality (23) is straightforwardly deduced from Definition (22). Definition (21) yields

$$ \Delta t^n G_{N+\frac{1}{2},j}^* = \sum_{\alpha=1}^{N} \left( l_\alpha h_j^{n+\frac{1}{2}} - h_{\alpha,j}^* \right) + \Delta t^n G_{\frac{1}{2},j}^* = h_j^{n+\frac{1}{2}} - \sum_{\alpha=1}^{N} h_{\alpha,j}^* = 0. \quad \square $$

## Correction step

The exchange terms defined in Proposition 2.1 are applied to correct the velocities and tracer values:

$$ (h_\alpha \varphi_\alpha)_j^{n+\frac{1}{2}} = (h_\alpha \varphi_\alpha)_j^* + \Delta t^n \left( \varphi_{\alpha+\frac{1}{2},j}^* G_{\alpha+\frac{1}{2},j}^* - \varphi_{\alpha-\frac{1}{2},j}^* G_{\alpha-\frac{1}{2},j}^* \right) \quad \varphi \in \{u, T\}, \quad (24) \quad \text{eq:correctionstep} $$

with either $\# = \star$ (explicit scheme) or $\# = n + \frac{1}{2}$ (implicit scheme). The interface velocities and tracers are defined in accordance with (9) by

$$ \varphi_{\alpha+\frac{1}{2},j}^{\#} = \begin{cases} \varphi_{\alpha,j}^{\#} & \text{if } G_{\alpha+\frac{1}{2},j}^{*} \le 0, \\ \varphi_{\alpha+1,j}^{\#} & \text{if } G_{\alpha+\frac{1}{2},j}^{*} > 0, \end{cases} \quad \varphi \in \{u, T\}. \quad (25) \quad \text{eq:decu} $$
### 2.1.2 Barotropic step

The barotropic step corresponds to system (12). Its first two lines correspond to the shallow water equations on $(h, h\bar{u})$ and have a severe CFL restriction on the timestep. Thus the timestep $\Delta t^n$ set in the baroclinic step is decomposed in $K$ small subtimesteps. We initialize the barotropic step with

$$h_j^{n+\frac{1}{2},0} := h_j^{n+\frac{1}{2}}, \quad \bar{u}_j^{n+\frac{1}{2},0} := \sum_{\alpha=1}^{N} l_{\alpha} u_{\alpha,j}^{n+\frac{1}{2}} \quad \text{and} \quad \sigma_{\alpha,j}^{n+\frac{1}{2}} := (u_{\alpha} - \bar{u})_j^{n+\frac{1}{2}}.$$

For $0 \le k \le K-1$, the $k$-th time subtimestep updates the values $(h_j^{n+\frac{1}{2},k}, \bar{u}_j^{n+\frac{1}{2},k})$ at time $\tau^k$ to the values $(h_j^{n+\frac{1}{2},k+1}, \bar{u}_j^{n+\frac{1}{2},k+1})$ at time $\tau^{k+1} = \tau^k + \delta t^k$. At the last subtimestep $K$ we have $\sum_{k=0}^{K-1} \delta t^k = \tau^K = \Delta t^n$ and we define $h_j^{n+1} := h_j^{n+\frac{1}{2},K}, \bar{u}_{\alpha,j}^{n+1} := \bar{u}_{\alpha,j}^{n+\frac{1}{2},K}$. The deviation $\sigma_\alpha$ and the tracer $T_\alpha$ are updated directly from $(\sigma_{\alpha,j}^{n+\frac{1}{2}}, T_{\alpha,j}^{n+\frac{1}{2}})$ to $(\sigma_{\alpha,j}^{n+1}, T_{\alpha,j}^{n+1})$.

The shallow water numerical fluxes for the density $h$ and the mean momentum $h\bar{u}$ at interfaces $j \pm \frac{1}{2}$ are denoted respectively with $f_{j \pm \frac{1}{2}}^{h,k}$ and $f_{j \pm \frac{1}{2}}^{h\bar{u},k}$ and the discretization of the topography source term with $S_j^k$. The scheme reads as

$$\left\{ \begin{aligned} h_j^{n+\frac{1}{2},k+1} &= h_j^{n+\frac{1}{2},k} - \frac{\delta t^k}{\Delta x} (f_{j+\frac{1}{2}}^{h,k} - f_{j-\frac{1}{2}}^{h,k}), \\ (h\bar{u})_j^{n+\frac{1}{2},k+1} &= (h\bar{u})_j^{n+\frac{1}{2},k} - \frac{\delta t^k}{\Delta x} (f_{j+\frac{1}{2}}^{h\bar{u},k} - f_{j-\frac{1}{2}}^{h\bar{u},k}) + \delta t^k S_j^k, \end{aligned} \right. \quad (26) \quad \text{eq:FV-SW}$$

with timestep $\delta t^k$ being constrained by a classical CFL condition of the form

$$\delta t^k \le \frac{\text{CFL}\Delta x}{\max_j |\bar{u}_j^{n+\frac{1}{2},k}| + \sqrt{gh_j^{n+\frac{1}{2},k}}} \quad (27) \quad \text{eq:CFL\_SWE}$$

for some constant Courant-Friedrichs-Lewy number CFL depending on the scheme.

The simplest method would be to update the remaining unknowns $T_\alpha$ and $\sigma_\alpha$ as frequently as $h$ and $\bar{u}$ with

$$\left\{ \begin{aligned} (h\sigma_\alpha)_j^{n+\frac{1}{2},k+1} &= (h\sigma_\alpha)_j^{n+\frac{1}{2},k} - \frac{\delta t^k}{\Delta x} (f_{j+\frac{1}{2}}^{h\sigma_\alpha,k} - f_{j-\frac{1}{2}}^{h\sigma_\alpha,k}), \\ (hT_\alpha)_j^{n+\frac{1}{2},k+1} &= (hT_\alpha)_j^{n+\frac{1}{2},k} - \frac{\delta t^k}{\Delta x} (f_{j+\frac{1}{2}}^{hT_\alpha,k} - f_{j-\frac{1}{2}}^{hT_\alpha,k}). \end{aligned} \right.$$

However it is preferable both from the computational cost point of view and to reduce the numerical diffusion to take advantage of the fact that equations on $hT_\alpha$ and $h\sigma_\alpha$ are decoupled from the ones on $(h, \bar{u})$ in (12). The mass fluxes of all the subtimesteps are gathered in

$$\Delta t^n \mathcal{F}_{j+\frac{1}{2}}^h := \sum_{k=0}^{K-1} \delta t^k f_{j+\frac{1}{2}}^{h,k} \quad \text{so that} \quad h_j^{n+1} = h_j^{n+\frac{1}{2}} - \frac{\Delta t^n}{\Delta x} \left( \mathcal{F}_{j+\frac{1}{2}}^h - \mathcal{F}_{j-\frac{1}{2}}^h \right). \quad (28) \quad \text{eq:summed_massfluc}$$

We decenter the numerical fluxes as before with

$$\mathcal{F}_{j+\frac{1}{2}}^{h\varphi_\alpha} := \varphi_{\alpha,j}^{n+\frac{1}{2}} \left( \mathcal{F}_{j+\frac{1}{2}}^{h} \right)^{+} - \varphi_{\alpha,j+1}^{n+\frac{1}{2}} \left( \mathcal{F}_{j+\frac{1}{2}}^{h} \right)^{-} \quad \text{for } \varphi_\alpha \in \{\sigma_\alpha, T_\alpha\} \quad (29) \quad \text{eq:fluxdecBT}$$

and update $\sigma_\alpha$ and $T_\alpha$ with only one large timestep $\Delta t^n$

$$\left\{ \begin{aligned} (h\sigma_\alpha)_j^{n+1} &= (h\sigma_\alpha)_j^{n+\frac{1}{2}} - \frac{\Delta t^n}{\Delta x} (\mathcal{F}_{j+\frac{1}{2}}^{h\sigma_\alpha} - \mathcal{F}_{j-\frac{1}{2}}^{h\sigma_\alpha}), \\ (hT_\alpha)_j^{n+1} &= (hT_\alpha)_j^{n+\frac{1}{2}} - \frac{\Delta t^n}{\Delta x} (\mathcal{F}_{j+\frac{1}{2}}^{hT_\alpha} - \mathcal{F}_{j-\frac{1}{2}}^{hT_\alpha}). \end{aligned} \right. \quad (30) \quad \text{eq:subcycling}$$

The strategy is detailed in Algorithm 3. A condition on $\Delta t^n$ ensuring the maximum principle on $\sigma_\alpha$ and $T_\alpha$ is given in Proposition 2.6 below. We observe in practice that this condition is always fulfilled for the baroclinic timestep $\Delta t^n$, see Section 3. For the sake of completeness, we however describe a strategy in Section 2.2.2 handling all the cases.
**Algorithm 3** Barotropic loop with subcycling strategy

**Input:** Initial data for barotropic loop $(h, u_\alpha, T_\alpha)^{n+\frac{1}{2}}$, final time for barotropic loop $\Delta t^n$

**Output:** $(h, u_\alpha, T_\alpha)^{n+1}$

1: $t \leftarrow 0$

2: $\Delta t \mathcal{F}^h \leftarrow 0$

3: $(h, \bar{u}, \sigma_\alpha, T_\alpha) \leftarrow (h^{n+\frac{1}{2}}, \sum_{\alpha=1}^N l_\alpha u_\alpha^{n+\frac{1}{2}}, u_\alpha^{n+\frac{1}{2}} - \bar{u}^{n+\frac{1}{2}}, T_\alpha^{n+\frac{1}{2}})$

4: **while** $t < \Delta t^n$ **do**

5: $\delta t \leftarrow \min(\text{computeDT\_barotropic}(h, \bar{u}), \Delta t^n - t)$ ▷ see (27)

6: $(h, \bar{u}, f^h) \leftarrow \text{SWE\_solver}(h, \bar{u}, \delta t)$ ▷ see (26)

7: $\Delta t \mathcal{F}^h \leftarrow \Delta t \mathcal{F}^h + \delta t f^h$ ▷ see (28)

8: $t \leftarrow t + \delta t$

9: **end while**

10: $(\sigma_\alpha, T_\alpha) \leftarrow \text{update\_transported\_quantities}(\sigma_\alpha, T_\alpha, \Delta t \mathcal{F}^h)$ ▷ see (30)

11: return $(h, u_\alpha = \sigma_\alpha + \bar{u}, T_\alpha)$

## 2.2 CFL conditions and properties of the scheme

In this section, we examine the additional time step restrictions and properties of the barotropic-baroclinic splitting. We focus on the positivity of the water heights, the discrete maximum principle for the tracer and the obtention of a discrete entropy inequality.

### 2.2.1 Properties of the baroclinic step

We exhibit two timestep restrictions that guarantee the non-negativity of the water heights and a maximum principle on $u_\alpha$ and $T_\alpha$ in the baroclinic step.

**Proposition 2.2.** Consider the prediction step (18) with the decentered choice (20) for the momentum flux. Suppose that the time step $\Delta t^n$ is small enough so that for all $\alpha$ and for all $j$

$$h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left( \left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^+ + \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^- \right) \geq 0. \quad (31) \quad \text{eq:positivity_pre}$$

Then the water heights $h_{\alpha,j}^*$ are non-negative and a discrete maximum principle on $u_\alpha$ and $T_\alpha$ holds:

$$\min(\varphi_{\alpha,j-1}^n, \varphi_{\alpha,j}^n, \varphi_{\alpha,j+1}^n) \leq \varphi_{\alpha,j}^* \leq \max(\varphi_{\alpha,j-1}^n, \varphi_{\alpha,j}^n, \varphi_{\alpha,j+1}^n) \quad \forall \varphi \in \{u, T\}.$$

*Proof.* We reproduce the proof of [AB03, Theorem 4.3]. The positivity of the water height is easy:

$$\begin{align*}
h_{\alpha,j}^* &= h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left(F_{j+\frac{1}{2}}^{h_\alpha} - F_{j-\frac{1}{2}}^{h_\alpha}\right) \\
&= h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left(\left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^+ - \left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^- - \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^+ + \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^-\right) \\
&\geq h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left(\left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^+ + \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^-\right)
\end{align*}$$

which is non-negative under hypothesis (31). We now prove the lower bound for $\varphi \in \{u, T\}$ updated with (20)

$$\begin{align*}
(h\varphi)_{\alpha,j}^* &= (h\varphi)_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left[ \varphi_{\alpha,j}^n \left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^+ - \varphi_{\alpha,j+1}^n \left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^- - \varphi_{\alpha,j-1}^n \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^+ + \varphi_{\alpha,j}^n \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^- \right] \\
&= \varphi_{\alpha,j}^n \left[ h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left( \left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^+ - \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^- \right) \right] + \varphi_{\alpha,j+1}^n \frac{\Delta t^n}{\Delta x} \left( \left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^- + \varphi_{\alpha,j-1}^n \frac{\Delta t^n}{\Delta x} \left( \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^+ \right) \right) \\
&\geq \min(\varphi_{\alpha,j-1}^n, \varphi_{\alpha,j}^n, \varphi_{\alpha,j+1}^n) \left( h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left[ \left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^+ - \left(F_{j+\frac{1}{2}}^{h_\alpha}\right)^- - \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^+ + \left(F_{j-\frac{1}{2}}^{h_\alpha}\right)^- \right] \right) \\
&= \min(\varphi_{\alpha,j-1}^n, \varphi_{\alpha,j}^n, \varphi_{\alpha,j+1}^n) h_{\alpha,j}^*.
\end{align*}$$
The inequality holds because the coefficient in front of $\varphi_{\alpha,j}^n$ is non-negative as soon as (31) holds. The upper bound is proved in a similar fashion. □

**Remark 2.3.** If condition (31) should be violated, it can be enforced by modifying the mass fluxes using the draining time technique [BNLM11]. This sets the mass flux to zero when the water height reaches zero within a time step. Note that it does not reduce the overall time step, but is an improved quadrature rule for the time integral of the numerical flux. For the terrain-following coordinates, this will only happen at wet-dry fronts.

#### (P:sumsa) Proposition 2.4. The correction on the velocities and the tracers (24) with the decentered choice (25) for the interface velocities verifies the maximum principle

$$ \min_{\beta}(\varphi_{\beta,j}^*) \leq \varphi_{\alpha,j}^{n+\frac{1}{2}} \leq \max_{\beta}(\varphi_{\beta,j}^*), \quad \varphi_{\alpha} \in \{u, T\} $$

for any time step $\Delta t^n$ in the implicit case $\sharp = n + \frac{1}{2}$. In the explicit case $\sharp = \star$ it verifies the maximum principle

$$ \min(\varphi_{\alpha-1,j}^*, \varphi_{\alpha,j}^*, \varphi_{\alpha+1,j}^*) \leq \varphi_{\alpha,j}^{n+\frac{1}{2}} \leq \max(\varphi_{\alpha-1,j}^*, \varphi_{\alpha,j}^*, \varphi_{\alpha+1,j}^*) $$

if the time step satisfies the bound

$$ \Delta t^n (G_{\alpha+\frac{1}{2},j}^*)^{-} + \Delta t^n (G_{\alpha-\frac{1}{2},j}^*)^{+} \leq h_{\alpha,j}^{n+\frac{1}{2}}. \quad (32) \quad \text{eq:CFLcor} $$

*Proof.* We refer the reader to [ABSM18, Lemma 3.11] for the implicit case. For the explicit case the proof is the same as in Proposition 2.2, with $\alpha$ playing the role of $j$ and $G^*$ the role of $F^{h_\alpha}$.

□

### 2.2.2 Properties of the barotropic step

#### (S:propBT)
We first check that the weighted sum of the deviation $\sum_{\alpha=1}^N l_\alpha \sigma_\alpha = \sum_{\alpha=1}^N l_\alpha (u_\alpha - \bar{u})$ remains null during the entire barotropic step.

#### (P:sumsigma)?
**Proposition 2.5.** The decentered choice (29) for the numerical flux ensures that $\sum_{\alpha=1}^N l_\alpha \sigma_{\alpha,j}^{n+1} = 0$ holds for all $j$ and for all $n$.

*Proof.* At the beginning of the barotropic step we have $\sum_{\alpha=1}^N l_\alpha \sigma_{\alpha,j}^{n+\frac{1}{2}} = 0$. Then

$$ \sum_{\alpha=1}^{N} l_{\alpha} \mathcal{F}_{j+\frac{1}{2}}^{h \sigma_{\alpha}} = (\mathcal{F}_{j+\frac{1}{2}}^{h})^{+} + \sum_{\alpha=1}^{N} l_{\alpha} \sigma_{\alpha,j}^{n+\frac{1}{2}} - (\mathcal{F}_{j+\frac{1}{2}}^{h})^{-} - \sum_{\alpha=1}^{N} l_{\alpha} \sigma_{\alpha,j+1}^{n+\frac{1}{2}} = 0 $$

and the result follows by summing (30) over all the layers. □

□

Second, we proof the non-negativity of the water height $h_j^{n+1}$ and a maximum principle for the velocity and the tracer.

#### (P:maxBT) Proposition 2.6. Suppose that

$$ h_j^{n+\frac{1}{2}} - \frac{\Delta t^n}{\Delta x} \left( \left( \mathcal{F}_{j+\frac{1}{2}}^h \right)^+ - \left( \mathcal{F}_{j-\frac{1}{2}}^h \right)^- \right) \geq 0 \quad (33) \quad \text{eq:non-negativity} $$

holds. Then $h_j^{n+1}$ is non-negative and the following maximum principle on $\varphi \in \{T_\alpha, \sigma_\alpha\}$ obtained with (30) holds:

$$ \min \left( \varphi_{j-1}^{n+\frac{1}{2}}, \varphi_j^{n+\frac{1}{2}}, \varphi_{j+1}^{n+\frac{1}{2}} \right) \leq \varphi_j^{n+1} \leq \max \left( \varphi_{j-1}^{n+\frac{1}{2}}, \varphi_j^{n+\frac{1}{2}}, \varphi_{j+1}^{n+\frac{1}{2}} \right). $$

*Proof.* The proof is the same as for Proposition 2.2. □

□
<caption>**Algorithm 4** Barotropic loop (with subcycling and non-negativity check)</caption>

**Input:** Initial data for barotropic loop $(h, u_\alpha, T_\alpha)^{n+\frac{1}{2}}$, final time for barotropic loop $\Delta t^n$
**Output:** $(h, \bar{u}, \sigma_\alpha, T_\alpha)^{n+1}$
1: $t \leftarrow 0$
2: $\Delta t \mathcal{F}^h \leftarrow 0$
3: $(h, \bar{u}, \sigma_\alpha, T_\alpha) \leftarrow (h^{n+\frac{1}{2}}, \sum_{\alpha=1}^N l_\alpha u_\alpha^{n+\frac{1}{2}}, u_\alpha^{n+\frac{1}{2}} - \bar{u}^{n+\frac{1}{2}}, T_\alpha^{n+\frac{1}{2}})$
4: $h_{\text{old}} \leftarrow h$
5: **while** $t < \Delta t^n$ **do**
6: &nbsp;&nbsp;&nbsp;&nbsp; $\delta t \leftarrow \min(\text{computeDT\_barotropic}(h, \bar{u}), \Delta t^n - t)$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ see (27)
7: &nbsp;&nbsp;&nbsp;&nbsp; $(f^h, f^{h\bar{u}}) \leftarrow \text{compute\_numerical\_SWE\_fluxes}(h, h\bar{u})$
8: &nbsp;&nbsp;&nbsp;&nbsp; **if** $h_{\text{old},j} - \frac{1}{\Delta x} \left( (\Delta t \mathcal{F}_{j+\frac{1}{2}}^h + \delta t f_{j+\frac{1}{2}}^h)^+ + (\Delta t \mathcal{F}_{j-\frac{1}{2}}^h + \delta t f_{j-\frac{1}{2}}^h)^- \right) \ge 0$ for all $j$ **then** &nbsp;&nbsp;&nbsp;&nbsp; ▷ see (33)
9: &nbsp;&nbsp;&nbsp;&nbsp; $(h, \bar{u}) \leftarrow \text{SWE\_solver}(h, \bar{u}, f^h, f^{h\bar{u}}, \delta t)$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ see (26)
10: &nbsp;&nbsp;&nbsp;&nbsp; $\Delta t \mathcal{F}^h \leftarrow \Delta t \mathcal{F}^h + \delta t f^h$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ see (28)
11: &nbsp;&nbsp;&nbsp;&nbsp; $t \leftarrow t + \delta t$
12: **else**
13: &nbsp;&nbsp;&nbsp;&nbsp; $(\sigma_\alpha, T_\alpha) \leftarrow \text{update\_transported\_quatities}(\sigma_\alpha, T_\alpha, \Delta t \mathcal{F}^h)$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ see (30)
14: &nbsp;&nbsp;&nbsp;&nbsp; $h_{\text{old}} \leftarrow h$
15: &nbsp;&nbsp;&nbsp;&nbsp; $\Delta t \mathcal{F}^h \leftarrow 0$
16: **end if**
17: **end while**
18: $(\sigma_\alpha, T_\alpha) \leftarrow \text{update\_transported\_quatities}(\sigma_\alpha, T_\alpha, \Delta t \mathcal{F}^h)$ &nbsp;&nbsp;&nbsp;&nbsp; ▷ see (30)
19: **return** $(h, u_\alpha = \sigma_\alpha + \bar{u}, T_\alpha)$

**Remark 2.7.** In practice, there is no guarantee that condition (33) applies. We therefore check condition (33) using the summed flux $\sum_{k=0}^m \delta t^k f^{h,k}$ instead of $\Delta t^n \mathcal{F}^h = \sum_{k=0}^{K-1} \delta t^k f^{h,k}$ for each $m$, $1 \le m \le K-1$. If condition (33) is not fulfilled for one $m$, we update the transported quantities $(h\sigma_\alpha)^{n+\frac{1}{2},m}$ and $(hT_\alpha)^{n+\frac{1}{2},m}$ in between according to (30), but using the mass fluxes summed up to $m$ and restart the summation of the mass fluxes from here. The strategy is the same as in [AB03] and is illustrated in more detail in Algorithm 4. The main argument of the proof is once again the one of Proposition 2.2.

### 2.2.3 Properties of the entire scheme

In (S:choices) We now consider the global scheme based on the barotropic-baroclinic splitting, going from $(h, u_\alpha, T_\alpha)^n$ to $(h, u_\alpha, T_\alpha)^{n+1}$ with the timestep $\Delta t^n$, consisting of

* the baroclinic part with

- the prediction step (18)-(21) with choice (20) for the momentum flux;

- the correction step (24)-(25) with Definition (22) of the exchange terms;

* the baroclinic part with the update of the shallow water system (26) and the update (30) of the deviation $\sigma_\alpha$ and the tracer $T_\alpha$ with the choice of flux (29).

We suppose that the CFL conditions (31), (32) and (33) are fulfilled. The entropy flux in the entropy inequality (13) for the continuous solutions is denoted by $F^E = gh^2\bar{u} + ghz_b\bar{u} + \sum_{\alpha=1}^N \frac{l_\alpha h u_\alpha^3}{2}$. We now show that this inequality also holds at the discrete level.

**(thm:DEI)** **Theorem 2.8.** Suppose that the scheme verifies the following additional assumptions:
* In the prediction step, there exists a numerical potential energy flux $F_{j \pm \frac{1}{2}}^{E^p}$ consistent with 0 such that

$$ (E^p)_j^* \le (E^p)_j^n - \frac{\Delta t^n}{\Delta x} (F_{j+\frac{1}{2}}^{E^p} - F_{j-\frac{1}{2}}^{E^p}) \quad \text{where} \quad E^p = g\frac{h^2}{2} + ghzb_j. \quad (34) \quad \text{eq:Ep} $$

* The shallow water scheme (26) satisfies a discrete entropy inequality (see Lemma 4.3 for a precise definition) under the subtimestep restriction (27) on $\delta t^k$.

Then there exists a numerical entropy flux $F_{j \pm \frac{1}{2}}^{E^p}$ consistent with the continuous entropy flux $F^{E^p}$ such that, under the same restrictions on $\Delta t^n$ and $\delta t^k$, the following discrete entropy inequality holds:

$$ E_j^{n+1} \le E_j^n + \frac{\Delta t^n}{\Delta x} (F_{j+\frac{1}{2}}^E - F_{j-\frac{1}{2}}^E). $$

**(rmk:sumnull)** **Remark 2.9.** Let us recall that in the baroclinic step (11) the total water height and thus the potential energy $E^p = g\frac{h^2}{2} + ghzb_j$ remain constant. At the discrete level the mass flux (18) may introduce some numerical diffusion and a decay of the potential energy. If the total mass flux $\sum_{\alpha=1}^N F_{j \pm \frac{1}{2}}^{h_\alpha}$ vanishes, we automatically get $h_j^{n+\frac{1}{2}} = h_j^n$ and an exact numerical counterpart of $\partial_t h = 0$ and $\partial_t E^p = 0$. In other words in this case (34) is an equality with $F_{j \pm \frac{1}{2}}^{E^p} = 0$.

**Proof of Theorem 2.8.** The proof is detailed in Section 4.2 and heavily relies on the splitting of the total entropy inequality “(13)=(14)+(15)” obtained in Proposition 1.6. We obtain a discrete counterpart of (15) for the baroclinic part in Lemma 4.1 and Lemma 4.2 for the prediction and correction steps respectively. Lemma 4.3 mimics (14) through the barotropic step with the subcycling strategy. The three inequalities are then put together to obtain the decay of the entropy over the entire timestep. $\square$

**(P:choice)** **Proposition 2.10.** Suppose that the shallow water scheme (26) preserves both the lake at rest equilibrium with null velocity $\bar{u}$ and constant free surface $h + z_b$ and the non-negativity of the water height under a CFL condition of the form (27). Consider the baroclinic mass flux for (18) given by

$$ F_{j+\frac{1}{2}}^{h_\alpha} = (h_\alpha \sigma_\alpha)_j^n := \begin{cases} (h_\alpha \sigma_\alpha)_j^n, & \text{if } h_{\alpha,j}^n < h_{\alpha,j+1}^n, \\ (h_\alpha \sigma_\alpha)_j^n, & \text{if } h_{\alpha,j}^n \ge h_{\alpha,j+1}^n. \end{cases} \quad (35) \quad \text{eq:choice} $$

Then the overall scheme

* preserves the lake at rest equilibrium;

* preserves the shallow water equation solutions:

$$ (\forall \alpha, \forall j, u_{\alpha,j}^n = \bar{u}_j^n) \implies (\forall \alpha, \forall j, u_{\alpha,j}^{n+1} = \bar{u}_j^{n+1}); $$

* verifies the first hypothesis of Theorem 2.8;

* preserves the non-negativity of the water height under the CFL condition $\Delta t^n \le \frac{0.5\Delta x}{\max_{\alpha,j} |\sigma_{\alpha,j}|}$.

**Proof.** We recall that $\sum_{\alpha=1}^N h_{\alpha,j}^n \sigma_{\alpha,j}^n = 0$ which yields $\sum_{\alpha=1}^N F_{j+\frac{1}{2}}^{h_\alpha} = 0$ and $h_j^{n+\frac{1}{2}} = h_j^n$, see Remark 2.9. In particular (34) holds with $F_{j \pm \frac{1}{2}}^{E^p} = 0$.

Consider a lake at rest initial data: there exists $C$ such that for all $j$, $h_j^n + (z_b)_j = C$ and for all $\alpha$ and for all $j$, $u_{\alpha,j}^n = 0$. Then all the deviations $\sigma_{\alpha,j}^n$ are equal to 0. The decentered choice (20) for the fluxes yields $F_{j+\frac{1}{2}}^{h_\alpha u_\alpha} = 0$ and thus $u_{\alpha,j}^* = 0$. Moreover, with the mass flux (35) we not only have $h_j^{n+\frac{1}{2}} = h_j^n$ but also $h_{\alpha,j}^* = h_{\alpha,j}^n = l_\alpha h_j^{n+\frac{1}{2}}$. Thus all the exchange terms (22) vanish and we also have $u_{\alpha,j}^{n+\frac{1}{2}} = 0$ after the correction step (24). The lake at rest equilibrium is thus preserved through the baroclinic
step. It is also preserved in the shallow water part by hypothesis. The redistribution (29)-(30) yields $\sigma_{\alpha,j}^{n+1} = 0$ and all the individual velocities are still null after the barotropic step.

We now turn to the case of a purely barotropic setting, where the velocity is uniform along every vertical column of fluid (but may vary along the horizontal): for all $\alpha$ and for all $j$, $u_{\alpha,j}^n = \bar{u}_j^n$ or in other words $\sigma_{\alpha,j}^n = 0$ identically. Then once again in the baroclinic step $F_{j+\frac{1}{2}}^{h_\alpha} = 0$, yielding $F_{j+\frac{1}{2}}^{h_\alpha u_\alpha} = 0$, $h_{\alpha,j}^* = h_{\alpha,j}^n$ and $u_{\alpha,j}^* = u_{\alpha,j}^n$. The exchange terms vanish and the correction step does nothing. In the barotropic step, $h$ and $\bar{u}$ vary in the resolution of the shallow water system (26). The redistribution over the layers (30) gives $\sigma_{\alpha,j}^{n+1} = 0$ and thus $u_{\alpha,j}^{n+1} = \bar{u}_j^{n+1}$.

We now turn to the non-negativity condition (31). We treat the worst case for the non-negativity, the case $F_{j+\frac{1}{2}}^{h_\alpha} > 0$, $F_{j-\frac{1}{2}}^{h_\alpha} < 0$. Then

$$h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left( \left( F_{j+\frac{1}{2}}^{h_\alpha} \right)^+ + \left( F_{j-\frac{1}{2}}^{h_\alpha} \right)^- \right) = h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left( h_{\alpha,j+\frac{1}{2}}^n |\sigma_{\alpha,j+\frac{1}{2}}^n| + h_{\alpha,j-\frac{1}{2}}^n |\sigma_{\alpha,j-\frac{1}{2}}^n| \right) \\ \geq h_{\alpha,j}^n \left( 1 - \frac{\Delta t^n}{\Delta x} (|\sigma_{\alpha,j+\frac{1}{2}}^n| + |\sigma_{\alpha,j-\frac{1}{2}}^n|) \right)$$

which is non-negative as soon as

$$\Delta t^n \leq \frac{0.5\Delta x}{\max_{\alpha,j} |\sigma_{\alpha,j}^n|}.$$

□

〈P:Rusanov〉

**Proposition 2.11.** Proposition 2.10 is also valid for the Rusanov mass flux

$$F_{j+\frac{1}{2}}^{h_\alpha} = \frac{h_{\alpha,j}^n \sigma_{\alpha,j}^n + h_{\alpha,j+1}^n \sigma_{\alpha,j+1}^n}{2} - A_{j+\frac{1}{2}} \frac{h_{\alpha,j+1}^n - h_{\alpha,j}^n}{2} \quad \text{with} \quad A_{j+\frac{1}{2}} = \max_{\substack{1 \le \alpha \le N, \\ 0 \le k \le 1}} |\sigma_{\alpha,j+k}^n| \quad (36) \text{eq:Rusanov}$$

with the difference that $\Delta t^n \leq \frac{\Delta x}{\max_{\alpha,j} |\sigma_{\alpha,j}^n|}$ is enough to ensure the non-negativity of the water height.

*Proof.* The water height does not vary during the baroclinic step on the lake at rest or shallow water solution because $\sigma_{\alpha,j}^n = 0$ and $A_{j \pm \frac{1}{2}} = 0$ in this case. As in the proof before, $F_{j \pm \frac{1}{2}}^{h_\alpha} = 0$ and thus $F_{j \pm \frac{1}{2}}^{h_\alpha u_\alpha} = 0$ and the exchange terms vanish. Thus also the velocity does not change during the baroclinic step. The barotropic step preserves the lake at rest solution and the shallow water solution by hypothesis. In the case of the shallow water solution $h$ and $\bar{u}$ vary but since the redistribution gives $\sigma_{\alpha,j}^{n+1} = 0$ it holds $u_{\alpha,j}^{n+1} = u_{\alpha,j}^n = 0$ in the case of the lake at rest and $u_{\alpha,j}^{n+1} = \bar{u}_j^{n+1}$ in the case of the shallow water solution.

Now we investigate the non-negativity condition (31). The worst case is when $F_{j+\frac{1}{2}}^{h_\alpha} > 0$ and $F_{j-\frac{1}{2}}^{h_\alpha} < 0$. Then

$$h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left( \left( F_{j+\frac{1}{2}}^{h_\alpha} \right)^+ + \left( F_{j-\frac{1}{2}}^{h_\alpha} \right)^- \right) = h_{\alpha,j}^n - \frac{\Delta t^n}{\Delta x} \left( h_{\alpha,j}^n \frac{\sigma_{\alpha,j}^n + A_{j+\frac{1}{2}}}{2} + h_{\alpha,j+1}^n \frac{\sigma_{\alpha,j+1}^n - A_{j+\frac{1}{2}}}{2} - h_{\alpha,j-1}^n \frac{\sigma_{\alpha,j-1}^n + A_{j-\frac{1}{2}}}{2} - h_{\alpha,j}^n \frac{\sigma_{\alpha,j}^n - A_{j-\frac{1}{2}}}{2} \right) \\ = h_{\alpha,j}^n \left( 1 - \frac{\Delta t^n}{\Delta x} \frac{A_{j+\frac{1}{2}} + A_{j-\frac{1}{2}}}{2} \right) + h_{\alpha,j+1}^n \frac{\Delta t^n}{\Delta x} \frac{A_{j+\frac{1}{2}} - \sigma_{\alpha,j+1}^n}{2} + h_{\alpha,j-1}^n \frac{\Delta t^n}{\Delta x} \frac{\sigma_{\alpha,j-1}^n + A_{j-\frac{1}{2}}}{2}$$

With $A_{j \pm 1/2} = \max(|\sigma_{\alpha,j}^n|, |\sigma_{\alpha,j \pm 1}^n|)$ the coefficients in front of $h_{\alpha,j+1}^n$ and $h_{\alpha,j-1}^n$ are non-negative. The coefficient in front of $h_{\alpha,j}^n$ is non-negative as soon as

$$\Delta t^n \leq \frac{\Delta x}{\max_{\alpha,j} |\sigma_{\alpha,j}^n|} \leq \frac{2\Delta x}{A_{j+\frac{1}{2}} + A_{j-\frac{1}{2}}} \quad (37) \text{eq:dt_rusanov}$$
We now turn to the variation of potential energy (34) in the prediction step. Under the hypothesis (37)

$$h_j^* = h_j^n \left( 1 - \frac{\Delta t^n A_{j+\frac{1}{2}} + A_{j-\frac{1}{2}}}{\Delta x} \right) + h_{j+1}^n \frac{\Delta t^n A_{j+\frac{1}{2}}}{\Delta x} + h_{j-1}^n \frac{\Delta t^n A_{j-\frac{1}{2}}}{\Delta x}$$

is a convex combination of $h_j^n, h_{j+1}^n$ and $h_{j-1}^n$. Hence we get with Jensen's inequality

$$(h_j^*)^2 \le (h_j^n)^2 \left( 1 - \frac{\Delta t^n A_{j+\frac{1}{2}} + A_{j-\frac{1}{2}}}{\Delta x} \right) + (h_{j+1}^n)^2 \frac{\Delta t^n A_{j+\frac{1}{2}}}{\Delta x} + (h_{j-1}^n)^2 \frac{\Delta t^n A_{j-\frac{1}{2}}}{\Delta x} \\ = (h_j^n)^2 + \frac{\Delta t^n}{\Delta x} \left[ A_{j+\frac{1}{2}} \frac{(h_{j+1}^n)^2 - (h_j^n)^2}{2} - A_{j-\frac{1}{2}} \frac{(h_j^n)^2 - (h_{j-1}^n)^2}{2} \right]$$

and we obtain (34) with $F_{j+\frac{1}{2}}^{E,p} = A_{j+\frac{1}{2}} \frac{E_{j+1}^{p,n} - E_j^{p,n}}{2}$. □

# 3 Numerical experiments

In this section we perform several numerical test cases. We investigate the convergence behaviour of the splitting towards an analytical solution of the hydrostatic Euler equations, the reduction of computational cost and the well-balancing property of the scheme for the lake at rest and for the geostrophic equilibrium.

The results presented here are all obtained using the Rusanov mass flux (36). In addition, we ran the same test cases again with the mass flux (35), which gave us very similar results.

## 3.1 Convergence study: analytical solution for the 2d hydrostatic Euler equations

According to [BBSM13] stationary solutions of the 2d $(x, z)$ hydrostatic and incompressible Euler equations are given by

$$z_b(x) = \bar{z}_b - h(x) - \frac{\alpha^2 \beta^2}{2g \sin^2(\beta h(x))}, \\ u(x, z) = \frac{\alpha \beta}{\sin(\beta h(x))} \cos(\beta(z - z_b(x))), \\ w(x, z) = \alpha \beta \left( \frac{\partial z_b}{\partial x} \frac{\cos(\beta(z - z_b(x)))}{\sin(\beta h(x))} + \frac{\partial h}{\partial x} \frac{\sin(\beta(z - z_b(x))) \cos(\beta h(x))}{\sin^2(\beta h(x))} \right)$$

for $\alpha, \beta \in \mathbb{R}$ and $h \in C^1(\Omega)$ non-negative such that $\sin(\beta h(x)) \ne 0$ for all $x \in \Omega$. For our testcase we choose the domain $\Omega = [-5, 5]$, the parameters $\alpha = 0.1, \beta = 1, \bar{z}_b = 2$ and $h(x) = 2 - e^{-x^2}$ and use this as initial data. In this setting the Froude number is $Fr = 0.04$. We compute until final time $t = 60$ s using the barotropic-baroclinic splitting with implicit treatment of the exchange terms and subcycling strategy and investigate the behaviour of the numerical solution when refining the cells in $x$-direction and refining the layers. Figure 3 illustrates the velocity field and the vertical velocity profile of $u$ at various positions $x$. The upper section of the figure displays the velocity field, while the red lines indicate the locations where the convergence behavior of the vertical velocity profile of $u$ is investigated. This analysis is presented in the second row of Figure 3. The $L^1$-errors and experimental orders of convergence are given in Table 3.1.

If we modify this testcase such that it becomes a non stationary testcase by adding a smooth perturbation of the water height, the scheme still converges with first order.

To examine numerical diffusion, we introduce a column of tracer on the left side of the bump in the bottom topography in the initial condition and monitor the $L^2$-norm of the tracer over time. The results are presented in Figure 4a. This investigation reveals that, in this testcase, the original scheme without splitting exhibits lower diffusion compared to the barotropic-baroclinic splitting.
<figure><img src="image_6.png" /><figcaption>Figure 3: Upper row: The velocity field for testcase 3.1, along with the highlighted positions for a detailed examination of the vertical velocity profile of $u$. Lower row: The convergence behavior of the vertical velocity profiles of $u$ at various positions $x$.</figcaption></figure>

<table><caption>Table 1: $L^1$-errors and experimental orders of convergence for testcase 3.1.</caption><thead><tr><th>cells × layers</th><th>h error</th><th>EOC</th><th>u error</th><th>EOC</th></tr></thead><tbody><tr><td>50 × 5</td><td>6.10e-04</td><td>—</td><td>8.37e-02</td><td>—</td></tr><tr><td>100 × 10</td><td>1.98e-04</td><td>1.62</td><td>6.05e-02</td><td>0.47</td></tr><tr><td>200 × 20</td><td>1.02e-04</td><td>0.96</td><td>3.74e-02</td><td>0.70</td></tr><tr><td>400 × 40</td><td>5.50e-05</td><td>0.88</td><td>2.09e-02</td><td>0.84</td></tr><tr><td>800 × 80</td><td>2.88e-05</td><td>0.94</td><td>1.11e-02</td><td>0.92</td></tr><tr><td>1600 × 160</td><td>1.58e-05</td><td>0.86</td><td>5.70e-03</td><td>0.96</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 4: L²-norm of the tracer to examine numerical diffusion. The solid line represents the splitting method, while the dashed line indicates the scheme without splitting.</figcaption></figure>

## 3.2 Computational cost

For the comparison of the different schemes (barotropic-baroclinic splitting (with or without subcycling strategy) or unsplit scheme, implicit or explicit treatment of exchange terms) we perform several runs of the respective scheme and measure the averaged runtime and the total error. As a testcase we use the analytical solution for the 2d hydrostatic Euler equations given in section 3.1 until final time $t = 1$ s. We perform the simulation three times with different Froude numbers given by different values for $\alpha$. The averaged (over 10 runs) measured runtime for the different schemes is plotted against the total error $\|(h_{err}, hu_{err})^T\|_2$ in Figure 5. Note that the gain in computational cost is especially large in the low Froude setting. In addition, in low Froude situations the subcycling strategy has a bigger impact on the computational cost than the explicit or implicit treatment of the exchange terms.

Furthermore we investigated the ratio of the computational cost of the barotropic-baroclinic splitting with subcycling and the scheme without splitting for different Froude numbers. The results are plotted in Figure 6 and show that the gain in computational cost is proportional to the Froude number.

## 3.3 Wind driven cavity

In this section we perform the test case of a two-layer fluid flow, initially at rest, subject to wind stress as it was performed in [ABPSM11a], but with constant density. We consider a fluid in a rectangular basin: the computational domain is given by the interval $[0, 3]$ with reflecting boundary conditions. Initially the fluid is at rest with a total water height of $h = 1$ m. The initial temperature (tracer) distribution is given by

$$T_0(x, z) = \begin{cases} 25, & \text{if } z - b \ge \frac{h}{2} \\ 8, & \text{otherwise.} \end{cases}$$

We impose a constant uniform wind stress in $x$-direction. As done in [ABPSM11a] we use a vertical viscosity of $\nu = 0.003 \text{ m}^2\text{s}^{-1}$, the friction coefficient $\kappa = 0.1 \text{ m}^{-1}$ and a wind velocity of $6 \text{ m}^{-1}$. We let the scheme run until final time $t = 600$ s.

The numerical results for the tracers are shown in Figure 7. The barotropic-baroclinic splitting converges much faster than the scheme without splitting. Note that the splitting seems to be less diffusive than the scheme without splitting. To investigate this further we again tracked the L²-norm of the tracer over time. The results are plotted in Figure 4b. This investigation reveals that, in contrast to testcase 3.1 in this testcase, the barotropic-baroclinic splitting is less diffusive than the original scheme without splitting.
<figure><img src="image_8.png" /><figcaption>Figure 5: Comparison of the computational cost and the total error $||(h_{\text{err}}, hu_{\text{err}})^{\top}||_2$ for testcase 3.1 with different Froude numbers. Top left: $Fr = 1.87$ (using $\alpha = 5$), top right: $Fr = 0.38$ (using $\alpha = 1$), bottom left: $Fr = 0.04$ (using $\alpha = 0.1$).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6: Comparison of the computational cost ratio (cost of barotropic-baroclinic splitting with subcycling over cost of the scheme without splitting) for various Froude numbers. The solid line represents the implicit version of the scheme, while the dashed line indicates the explicit version.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 7: Tracer distribution of the wind driven cavity (testcase 3.3) at time t = 600 s.</figcaption></figure>

So far, horizontal viscosity was neglected in the systems. We add the horizontal viscosity ($\nu_{hor} = \nu = 0.003 \text{ m}^2\text{s}^{-1}$) both to the splitting and the scheme without splitting. The results are plotted in Figure 8.

Since the horizontal viscosity is described by second order terms it gives a severe time step restriction depending on $\Delta x^2$. Thus the grids in the simulations are not as refined as before. But again, the splitting converges faster than the scheme without splitting.

## 3.4 Well-balancing for the lake at rest

The barotropic-baroclinic splitting is well-balanced for the lake at if the shallow water solver (26) used in the barotropic step is well-balanced. Since we used the hydrostatic reconstruction of Audusse et al. [ABB+04, CN17], which is well-balanced for the lake at rest, our barotropic-baroclinic splitting is also well-balanced. We verify this property with the following numerical test case taken from [CN22]: On the computational domain $[0, 4] \times [0, 2]$ we consider a volcano shaped bottom topography with a lake in the crater of the volcano and water surrounding the volcano. The initial data is given by

$$ 
\begin{aligned}
r(x, y) &= 2(x - 2)^2 + 4(y - 1)^2, \\
z_b(x, y) &= \begin{cases} 1 - 0.8 \exp(-r(x, y)), & r(x, y) \le \ln \frac{8}{5} \\ 0.8 \exp(-r(x, y)), & \text{else} \end{cases} \\
h(x, y) &= \begin{cases} \max(0.45 - z_b(x, y), 0), & r(x, y) \le \ln \frac{8}{5} \\ \max(0.3 - z_b(x, y), 0), & \text{else} \end{cases} \\
u \equiv 0, \quad v \equiv 0
\end{aligned}
$$

and is plotted in Figure 9. We perform the simulation on a grid of 200 × 100 cells with 10 layers until final time $t = 2$. The results are shown in Figure 10. Since the deviation of the numerical solution from the lake at rest is in order of the machine accuracy the scheme is well-balanced.
<figure><img src="image_11.png" /></figure>

(a) unsplit scheme, 200 cells, 50 layers

(b) unsplit scheme, 800 cells, 50 layers

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

(c) splitting, 200 cells, 50 layers

(d) splitting, 800 cells, 50 layers

Figure 8: Tracer distribution of the wind driven cavity (testcase 3.3) at time $t = 600$ s with horizontal viscosity $\nu_{hor} = \nu = 0.003 \ \text{m}^2\text{s}^{-1}$.

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

(a) Bottom $z_b$ and water surface $\eta$.

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

(b) Cut at position $y = 1$.

Figure 9: Initial data of lake at rest (testcase 3.4).

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

(a) Error of water surface $\eta$.

(b) $\sum_{\alpha=1}^{N} |u_{\alpha}|$

(c) $\sum_{\alpha=1}^{N} |v_{\alpha}|$

Figure 10: Deviation from the lake at rest steady state in testcase 3.4 at time $t = 2$.
## 3.5 Western boundary current and "linear well-balancing" for the geostrophic equilibrium

In this section we investigate the effect of the "linearly well-balancing" strategy of [GTCDGR24] on the preservation of geostrophic equilibrium in the context of the barotropic-baroclinic splitting.

The goal of this strategy is to preserve the geostrophic equilibrium $g\nabla(h + b) = f(v, -u)^{\top}$ with a linear velocity field. To obtain this we first perform a linear bottom reconstruction $b_{ij}^R(x, y) = b_{ij} + \mathcal{M}_{ij}(x - x_i) + \mathcal{N}_{ij}(y - y_j)$ to have access to $\nabla b_{ij}^R(x, y)$. Using the cell values $u_{ij}$ and $v_{ij}$ we construct a piecewise linear divergence free velocity field

$$u_{ij}^R(x, y) = u_{ij} + \mathcal{A}_{ij}(x - x_i) + \mathcal{B}_{ij}(y - y_j),$$

$$v_{ij}^R(x, y) = v_{ij} + \mathcal{C}_{ij}(x - x_i) - \mathcal{A}_{ij}(y - y_j)$$

solving a least squares minimization problem using the neighbouring cells. Note that with this choice of the coefficients the velocity field $(u^R, v^R)$ is already divergence free. Then we construct

$$h_{ij}^R(x, y) = h_{ij} + \alpha_{ij}(x - x_i) + \beta_{ij}(y - y_j) + \gamma_{ij}(x - x_i)^2 + \delta_{ij}(y - y_j)^2 + \lambda_{ij}(x - x_i)(y - y_j)$$

such that

$$\nabla h_{ij}^R = \frac{f}{g}(v_{ij}^R, -u_{ij}^R)^{\top} - \nabla b_{ij}.$$

This condition is used to determine the coefficients

$$\alpha_{ij} = \frac{f}{g}v_{ij} - \mathcal{M}_{ij}, \quad \beta_{ij} = -\frac{f}{g}u_{ij} - \mathcal{N}_{ij}, \quad \gamma_{ij} = \frac{f}{2g}\mathcal{C}_{ij}, \quad \delta_{ij} = -\frac{f}{2g}\mathcal{B}_{ij}, \quad \lambda_{ij} = -\frac{f}{g}\mathcal{A}_{ij}.$$

After having reconstructed $h_{ij}^R$, $u_{ij}^R$ and $v_{ij}^R$, we compute the numerical fluxes using the evaluations of the reconstructions at the midpoint of the cell interfaces. Finally we compute the source term via

$$S_{ij} = \int_{C_{ij}} f h_{ij}^R (v_{ij}^R, -u_{ij}^R)^{\top} = \int_{C_{ij}} g h_{ij}^R (\nabla h_{ij}^R + \nabla b_{ij}) \\ = \int_{C_{ij}} \frac{1}{2} g (\nabla h_{ij}^R)^2 + \int_{C_{ij}} g h_{ij}^R \nabla b_{ij} = \int_{\Gamma_{ij}} \frac{1}{2} g (h_{ij}^R)^2 \cdot \vec{n} + \int_{C_{ij}} g h_{ij}^R \nabla b_{ij}.$$

To test this strategy, we perform a large scale simulation and perform the testcase Henry Stommel investigated in his paper [Sto48]. Stommel considered a rectangular ocean $\Omega = [0, \lambda] \times [0, b]$ with $\lambda = 10^7$ m and $b = 2\pi \cdot 10^6$ m. The ocean is considered as a homogeneous layer of constant depth of 200 m. At initial time the ocean is at rest. The imposed wind stress is given by $\tau^x = -F \cos(\frac{\pi y}{b})$ and $\tau^y = 0$ with $F = 0.1$ Nm$^{-2}$. The bottom friction coefficient is given by $\kappa = 0.02$. As boundary conditions we use walls. Stommel has shown that using the beta-plane approximation for the Coriolis parameter $f = f_0 + \beta_0 y$ results in a stationary western boundary current that is in geostrophic equilibrium. In contrast to Stommel we choose a midlatitude setting for our simulation, so we set $f_0 = 2.5 \cdot 10^{-5}$ and $\beta_0 = 10^{-11}$. The exact linearized solution of Stommel is plotted in Figure 11. The maximum speed of the western boundary current is $2$ ms$^{-1}$. The width of the western boundary current is less than $100$ km.

For our simulation we use the fastest version of the barotropic-baroclinic splitting (with explicit treatment of the exchange terms and the subcycling strategy). We run the scheme on a grid of $200 \times 120$ cells and 50 layers until final time $t = 10$ years and compare the numerical results obtained with and without the use of the well-balancing strategy detailed above. Since the western boundary current is only created when the beta-plane approximation $f = f_0 + \beta_0 y$ is done, we use a variable Coriolis coefficient $f = f_j$ in each cell. The western boundary current appears more narrow and with higher speed (approximately $1.39$ ms$^{-1}$ compared to $0.51$ ms$^{-1}$) using the well-balancing strategy. The results of the simulations are plotted in Figure 12. One can observe that the modification of the shallow water solver by using the linear well-balancing strategy has a big impact on the solution, even on with a first order scheme on a coarse mesh.
<figure><img src="image_16.png" /><figcaption>Figure 11: Exact linearized solution for a western boundary current of Stommel.</figcaption></figure>

It is quite difficult to compare the theoretical results of Stommel’s illustrated on Figure 11 with our numerical simulations on Figure 12 for at least two reasons:

• Stommel model is a one layer linearized model in primitive variables, while we solve a multilayer nonlinear model with conservative variables. The vertical diffusion coefficient has been chosen in accordance with the choice in OGCMs.

• The boundary layer is strongly influenced by the dissipating mechanisms taken in consideration in the equation and the presence of nonlinear effects [Bry63], [Bö86]. In our simulation the boundary current width and the grid size approximately are comparable. As a result the numerical dissipation is quite strong near the western boundary, resulting in a non-controlled dissipating mechanism not explicitly appearing in the equation but of the same order of magnitude, if not higher, than the frictionnal term.

Our interpretation of the numerical results is the following. Since we are using a Rusanov scheme we can argue that the numerical diffusion writes as a (space and time varying) laplacian term and that the solution of the model is closer to the solution of (nonlinear) Munk model. Then the shape of the western boundary layer depends on the Reynolds number $Re = \frac{L_sL}{A}$ where $L_s$ is the interior velocity flow and $A$ is the diffusion coefficient. With the original Rusanov scheme the numerical diffusion is quite high and the Reynolds number is low. On the other hand the well-balancing strategy results in a much lower numerical diffusion and a higher Reynolds number. In those circumstances the elongation of the western boundary current in $h$ and the eddy in velocity is documented [Bry63, Figure 3].

# 4 Proofs of energy (in)equalities

In the following we are presenting the proofs concerning the continuous and discrete energy (in)equalities presented in the sections before.

## 4.1 Continuous energy equalities for smooth solutions

We start with the continuous energy equalities.

*Proof of Proposition 1.6.* The energy equality (13) for the original multilayer shallow water model can be found in [ABPSM11b, Proposition 4.2].

We now turn to the energy equality (14) fulfilled by the smooth solutions of the barotropic sys-
<figure><img src="image_17.png" /><figcaption>(a) Deviation from mean height.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(b) Depth averaged zonal velocity $\bar{u}$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(c) Depth averaged meridional velocity $\bar{v}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(d) Deviation from mean height.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(e) Depth averaged zonal velocity $\bar{u}$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(f) Depth averaged meridional velocity $\bar{v}$.</figcaption></figure>

Figure 12: Numerical solution on a grid with 400 × 240 cells ($\Delta x = 25$ km without well-balancing strategy for the geostrophic equilibrium (top) and with well-balancing strategy (bottom)).

<figure><img src="image_23.png" /><figcaption>(a) No well-balancing.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) With well-balancing.</figcaption></figure>

Figure 13: Streamplot of the depth averaged velocities obtained without (a) and with (b) linear well-balancing strategy for the geostrophic equilibrium.
tem (10). We have

$$
\begin{align*}
-\frac{\partial}{\partial t} \left( \frac{gh^2}{2} + ghzb + \sum_{\alpha=1}^{N} \frac{h_\alpha u_\alpha^2}{2} \right) &= - \left( gh + gz_b + \sum_{\alpha=1}^{N} \frac{l_\alpha u_\alpha^2}{2} \right) \frac{\partial h}{\partial t} - \sum_{\alpha=1}^{N} h_\alpha u_\alpha \frac{\partial u_\alpha}{\partial t} \\
&= \left( gh + gz_b + \sum_{\alpha=1}^{N} \frac{l_\alpha u_\alpha^2}{2} \right) \frac{\partial}{\partial x} (h\bar{u}) + \sum_{\alpha=1}^{N} h_\alpha u_\alpha \left( \bar{u} \frac{\partial u_\alpha}{\partial x} + g \frac{\partial h}{\partial x} + g \frac{\partial z_b}{\partial x} \right) \\
&= \left( 2gh\bar{u} + gz_b\bar{u} + \sum_{\alpha=1}^{N} \frac{l_\alpha \bar{u} u_\alpha^2}{2} \right) \frac{\partial h}{\partial x} + \left( gh^2 + ghzb + \sum_{\alpha=1}^{N} \frac{h_\alpha u_\alpha^2}{2} \right) \frac{\partial \bar{u}}{\partial x} + \sum_{\alpha=1}^{N} h_\alpha u_\alpha \bar{u} \frac{\partial u_\alpha}{\partial x} + gh\bar{u} \frac{\partial z_b}{\partial x} \\
&= \frac{\partial}{\partial x} \left( gh^2 \bar{u} + ghzb \bar{u} + \sum_{\alpha=1}^{N} \frac{h_\alpha u_\alpha^2 \bar{u}}{2} \right).
\end{align*}
$$

We now turn to the baroclinic System (11) and Inequality (15). With $h_\alpha = l_\alpha h$ and using $\frac{\partial}{\partial t}h_\alpha = 0$ and Identity (7) we obtain

$$
\frac{\partial}{\partial t} \left( \frac{h_\alpha u_\alpha^2}{2} \right) + \frac{\partial}{\partial x} \left( \frac{h_\alpha u_\alpha^2 (u_\alpha - \bar{u})}{2} \right) = \left( u_\alpha u_{\alpha+\frac{1}{2}} - \frac{u_\alpha^2}{2} \right) G_{\alpha+\frac{1}{2}} - \left( u_\alpha u_{\alpha-\frac{1}{2}} - \frac{u_\alpha^2}{2} \right) G_{\alpha-\frac{1}{2}}.
$$

The summation over the layers and the upwinding choice (9) yield the result

$$
\begin{align*}
\frac{\partial}{\partial t} \left( \sum_{\alpha=1}^{N} \frac{h_\alpha u_\alpha^2}{2} \right) + \frac{\partial}{\partial x} \left( \sum_{\alpha=1}^{N} \frac{h_\alpha u_\alpha^2 (u_\alpha - \bar{u})}{2} \right) &= \sum_{\alpha=1}^{N} (u_\alpha - u_{\alpha+1}) \left( u_{\alpha+\frac{1}{2}} - \frac{u_{\alpha+1} + u_\alpha}{2} \right) G_{\alpha+\frac{1}{2}} \\
&= -\frac{1}{2} \sum_{\alpha=1}^{N} (u_\alpha - u_{\alpha+1})^2 |G_{\alpha+\frac{1}{2}}| \le 0.
\end{align*}
$$

□

4.2 Discrete entropy inequality

(S:DEIproofS) We now prove Theorem 2.8.

Prediction step of the baroclinic step

We start with the discrete counterpart of (15) through the baroclinic part. With $E_\alpha^c := \frac{h_\alpha u_\alpha^2}{2}$ the kinetic energy in layer $\alpha$ and $E_\alpha^c \sigma_\alpha := \frac{h_\alpha u_\alpha^2 (u_\alpha - \bar{u})}{2}$ the baroclinic energy flux this inequality simply writes

$$
\frac{\partial}{\partial t} \left( E^p + \sum_{\alpha=1}^{N} E_{\alpha}^{c} \right) + \frac{\partial}{\partial x} \left( \sum_{\alpha=1}^{N} E_{\alpha}^{c} \sigma_{\alpha} \right) \leq 0.
$$

(L:pred) **Lemma 4.1.** Suppose that (34) holds, i.e. that the potential energy $E^p$ decreases through the prediction step (18). Define

$$
F_{j+\frac{1}{2}}^{E_{\alpha}^{c} \sigma_{\alpha}} = \frac{1}{2} (u_{\alpha,j})^{n} (F_{j+\frac{1}{2}}^{h_{\alpha}})^{+} - \frac{1}{2} (u_{\alpha,j+1})^{n} (F_{j+\frac{1}{2}}^{h_{\alpha}})^{-}
$$

Then $F_{j+\frac{1}{2}}^{E_\alpha^c \sigma_\alpha} + F_{j+\frac{1}{2}}^{E^p}$ is consistent with $E_\alpha^c \sigma_\alpha$ and the following discrete entropy inequality holds

$$
E_j^* \le E_j^n - \frac{\Delta t^n}{\Delta x} \left( F_{j+\frac{1}{2}}^{E^p} \sum_{\alpha=1}^N F_{j+\frac{1}{2}}^{E_\alpha^c \sigma_\alpha} - F_{j-\frac{1}{2}}^{E^p} \sum_{\alpha=1}^N F_{j-\frac{1}{2}}^{E_\alpha^c \sigma_\alpha} \right).
$$
Proof. We prove that

$$ (E_{\alpha}^{c})^{*}_{j} \leq (E_{\alpha}^{c})^{n}_{j} - \frac{\Delta t^{n}}{\Delta x} (F_{j+\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}} - F_{j-\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}}) $$

and the result follows by summation over the layers and with (34). With $\lambda = \frac{\Delta t^n}{\Delta x}$ we have

$$ \begin{aligned} h_{\alpha,j}^{*} \Delta := & h_{\alpha,j}^{*} \left[ \left( \frac{h_{\alpha} u_{\alpha}^{2}}{2} \right)_{j}^{*} - \left( \frac{h_{\alpha} u_{\alpha}^{2}}{2} \right)_{j}^{n} + \lambda \left( F_{j+\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}} - F_{j-\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}} \right) \right] \\ = & -\lambda h_{\alpha,j}^{n} u_{\alpha,j}^{n} \left( F_{j+\frac{1}{2}}^{h_{\alpha} u_{\alpha}} - F_{j-\frac{1}{2}}^{h_{\alpha} u_{\alpha}} \right) + \frac{\lambda^2}{2} \left( F_{j+\frac{1}{2}}^{h_{\alpha} u_{\alpha}} - F_{j-\frac{1}{2}}^{h_{\alpha} u_{\alpha}} \right)^2 + \frac{\lambda}{2} h_{\alpha,j}^{n} (u_{\alpha,j})^2 \left( F_{j+\frac{1}{2}}^{h_{\alpha}} - F_{j-\frac{1}{2}}^{h_{\alpha}} \right) \\ & + \lambda h_{\alpha,j}^{n} \left( F_{j+\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}} - F_{j-\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}} \right) - \lambda^2 \left( F_{j+\frac{1}{2}}^{h_{\alpha}} - F_{j-\frac{1}{2}}^{h_{\alpha}} \right) \left( F_{j+\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}} - F_{j-\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}} \right) \end{aligned} $$

We write $F_{j+\frac{1}{2}}^{h_{\alpha} u_{\alpha}} = u_{\alpha,j+\frac{1}{2}}^{n} F_{j+\frac{1}{2}}^{h_{\alpha}}$ and $F_{j+\frac{1}{2}}^{E_{\alpha}^{c}\sigma_{\alpha}} = \frac{(u_{\alpha,j+\frac{1}{2}}^{n})^2}{2} F_{j+\frac{1}{2}}^{h_{\alpha}}$ where $u_{\alpha,j+\frac{1}{2}}^{n} = \begin{cases} u_{\alpha,j}^{n}, & \text{if } F_{j+\frac{1}{2}}^{h_{\alpha}} \geq 0, \\ u_{\alpha,j+1}^{n}, & \text{otherwise.} \end{cases}$

Straightforward computations yield

$$ h_{\alpha,j}^{*} \Delta = \frac{\lambda}{2} h_{\alpha,j}^{n} \left[ (u_{\alpha,j}^{n} - u_{\alpha,j+\frac{1}{2}}^{n})^{2} F_{j+\frac{1}{2}}^{h_{\alpha}} - (u_{\alpha,j}^{n} - u_{\alpha,j-\frac{1}{2}}^{n})^{2} F_{j-\frac{1}{2}}^{h_{\alpha}} \right] + \frac{\lambda^2}{2} F_{j+\frac{1}{2}}^{h_{\alpha}} F_{j-\frac{1}{2}}^{h_{\alpha}} (u_{\alpha,j+\frac{1}{2}}^{n} - u_{\alpha,j-\frac{1}{2}}^{n})^{2}. $$

When $F_{j+\frac{1}{2}}^{h_{\alpha}} F_{j-\frac{1}{2}}^{h_{\alpha}} \leq 0$ the righthand side is clearly nonpositive. Otherwise

$$ \left\{ \begin{array}{ll} \Delta = -\frac{\lambda}{2h_{\alpha,j}^{*}}(u_{\alpha,j}^{n} - u_{\alpha,j-1}^{n})^{2}F_{j-\frac{1}{2}}^{h_{\alpha}} \left(h_{\alpha,j}^{n} - \lambda(F_{j+\frac{1}{2}}^{h_{\alpha}})^{+}\right) & \text{if } F_{j+\frac{1}{2}}^{h_{\alpha}} > 0 \text{ and } F_{j-\frac{1}{2}}^{h_{\alpha}} > 0 \\ \\ \Delta = \frac{\lambda}{2h_{\alpha,j}^{*}}(u_{\alpha,j}^{n} - u_{\alpha,j+1}^{n})^{2}F_{j+\frac{1}{2}}^{h_{\alpha}} \left(h_{\alpha,j}^{n} + \lambda(F_{j-\frac{1}{2}}^{h_{\alpha}})^{-}\right) & \text{if } F_{j+\frac{1}{2}}^{h_{\alpha}} < 0 \text{ and } F_{j-\frac{1}{2}}^{h_{\alpha}} < 0 \end{array} \right. $$

which in nonpositive under the CFL condition (31). $\square$

### Correction step of the baroclinic step

#### (L:cor) Lemma 4.2. Consider Definition (22) for the exchange terms and Definition (9) for the interface velocities. The total energy decreases in the correction step (24) with implicit interface velocities $\sharp = n + \frac{1}{2}$. It is also true in the explicit case $\sharp = \star$ under the additional restriction on the timestep (32).

Proof. The total water height $h$ does not vary during the correction step, so we only have to prove that the total kinetic energy decreases. Definitions (23) and (24) yield in the explicit case

$$ \begin{aligned} \frac{2h_{\alpha,j}^{n+\frac{1}{2}}}{\Delta t^n} \left( (E_{\alpha}^c)_j^{n+\frac{1}{2}} - (E_{\alpha}^c)_j^* \right) = & -h_{\alpha,j}^* G_{\alpha+\frac{1}{2},j}^* \left( u_{\alpha+\frac{1}{2},j}^* - u_{\alpha,j}^* \right)^2 + h_{\alpha,j}^* G_{\alpha+\frac{1}{2},j}^* \left( u_{\alpha+\frac{1}{2},j}^* \right)^2 \\ & + h_{\alpha,j}^* G_{\alpha-\frac{1}{2},j}^* \left( u_{\alpha-\frac{1}{2},j}^* - u_{\alpha,j}^* \right)^2 - h_{\alpha,j}^* G_{\alpha-\frac{1}{2},j}^* \left( u_{\alpha-\frac{1}{2},j}^* \right)^2 \\ & + \Delta t^n \left( u_{\alpha+\frac{1}{2},j}^* G_{\alpha+\frac{1}{2},j}^* - u_{\alpha-\frac{1}{2},j}^* G_{\alpha-\frac{1}{2},j}^* \right)^2 \end{aligned} $$

We replace $h_{\alpha,j}^*$ by $h_{\alpha,j}^{n+\frac{1}{2}} - \Delta t^n (G_{\alpha+\frac{1}{2},j}^* - G_{\alpha-\frac{1}{2},j}^*)$ in the second and fourth terms and we arrive at

$$ \begin{aligned} \frac{h_{\alpha,j}^{n+\frac{1}{2}} (u_{\alpha,j}^{n+\frac{1}{2}})^2 - h_{\alpha,j}^*(u_{\alpha,j}^*)^2}{\Delta t^n} = & -\frac{h_{\alpha,j}^*}{h_{\alpha,j}^{n+\frac{1}{2}}} G_{\alpha+\frac{1}{2},j}^* \left( u_{\alpha+\frac{1}{2},j}^* - u_{\alpha,j}^* \right)^2 + \frac{h_{\alpha,j}^*}{h_{\alpha,j}^{n+\frac{1}{2}}} G_{\alpha-\frac{1}{2},j}^* \left( u_{\alpha-\frac{1}{2},j}^* - u_{\alpha,j}^* \right)^2 \\ & + G_{\alpha+\frac{1}{2},j}^* (u_{\alpha+\frac{1}{2},j}^*)^2 - G_{\alpha-\frac{1}{2},j}^* (u_{\alpha-\frac{1}{2},j}^*)^2 \\ & + \frac{\Delta t^n}{h_{\alpha,j}^{n+\frac{1}{2}}} G_{\alpha+\frac{1}{2},j}^* G_{\alpha-\frac{1}{2},j}^* (u_{\alpha+\frac{1}{2},j}^* - u_{\alpha-\frac{1}{2},j}^*)^2 \end{aligned} $$
We now write $G_{\alpha+\frac{1}{2},j}^*G_{\alpha-\frac{1}{2},j}^* = ((G_{\alpha+\frac{1}{2},j}^*)^+ - (G_{\alpha+\frac{1}{2},j}^*)^-)((G_{\alpha-\frac{1}{2},j}^*)^+ - (G_{\alpha-\frac{1}{2},j}^*)^-)$ and use the de-centered choice (9) to get

$$2 \frac{(E^c)_j^{n+\frac{1}{2}} - (E^c)_j^*}{\Delta t^n} \le -\sum_{\alpha=1}^N (G_{\alpha+\frac{1}{2},j}^* + (u_{\alpha+1,j}^* - u_{\alpha,j}^*)^2) \left[ \frac{h_{\alpha,j}^* - \Delta t^n (G_{\alpha-\frac{1}{2},j}^*)^+}{h_{\alpha,j}^{n+\frac{1}{2}}} \right. \\ \left. - \sum_{\alpha=1}^N (G_{\alpha-\frac{1}{2},j}^* - (u_{\alpha-1,j}^* - u_{\alpha,j}^*)^2) \left[ \frac{h_{\alpha,j}^* - \Delta t^n (G_{\alpha+\frac{1}{2},j}^*)^-}{h_{\alpha,j}^{n+\frac{1}{2}}} \right] \right].$$

The bracket terms are non-negative under the CFL condition (32).

The implicit case $\# = n + \frac{1}{2}$ which is simpler. We write

$$\frac{h_{\alpha,j}^*(u_{\alpha,j}^{n+\frac{1}{2}})^2 - h_{\alpha,j}^*(u_{\alpha,j}^*)^2}{\Delta t^n} = -G_{\alpha+\frac{1}{2},j}^* (u_{\alpha+\frac{1}{2},j}^{n+\frac{1}{2}} - u_{\alpha,j}^{n+\frac{1}{2}})^2 + G_{\alpha-\frac{1}{2},j}^* (u_{\alpha-\frac{1}{2},j}^{n+\frac{1}{2}} - u_{\alpha,j}^{n+\frac{1}{2}})^2 \\ + G_{\alpha+\frac{1}{2},j}^* (u_{\alpha+\frac{1}{2},j}^{n+\frac{1}{2}})^2 - G_{\alpha-\frac{1}{2},j}^* (u_{\alpha-\frac{1}{2},j}^{n+\frac{1}{2}})^2 \\ - G_{\alpha+\frac{1}{2},j}^* (u_{\alpha+\frac{1}{2},j}^{n+\frac{1}{2}})^2 + G_{\alpha-\frac{1}{2},j}^* (u_{\alpha-\frac{1}{2},j}^{n+\frac{1}{2}})^2$$

Choice (9) yields the result

$$2 \frac{(E^c)_j^{n+\frac{1}{2}} - (E^c)_j^*}{\Delta t^n} \le -\sum_{\alpha=1}^N \left[ \left( G_{\alpha+\frac{1}{2},j}^* \right)^+ \left( u_{\alpha,j+1}^{n+\frac{1}{2}} - u_{\alpha,j}^{n+\frac{1}{2}} \right)^2 + \left( G_{\alpha-\frac{1}{2},j}^* \right)^- \left( u_{\alpha,j-1}^{n+\frac{1}{2}} - u_{\alpha,j}^{n+\frac{1}{2}} \right)^2 \right]. \quad \square$$

## Barotropic step

We now turn to the discrete entropy inequality in the barotropic loop. Let us recall that the physically relevant solutions of the shallow water equation (first two equations of (12)) verify the entropy inequality

$$\partial_t E^{sw} + \partial_x f^{sw} \le 0$$

where denote by $E^{sw} = \frac{gh^2}{2} + ghzb + \frac{h\bar{u}^2}{2}$ and $f^{sw} = gh^2\bar{u} + ghzb\bar{u} + \frac{h\bar{u}^3}{2}$

**(L:EIBT) Lemma 4.3.** Suppose that the shallow water scheme (26) is entropy satisfying under Condition (27): there exists some discrete numerical entropy fluxes $f_{j \pm \frac{1}{2}}^{sw,k}$ consistent with the continuous entropy flux $f^{sw}$ such that

$$E_j^{sw,n+\frac{1}{2},k+1} \le E_j^{sw,n+\frac{1}{2},k} + \frac{\delta t^k}{\Delta x} (f_{j+\frac{1}{2}}^{sw,k} - f_{j-\frac{1}{2}}^{sw,n,k}). \quad (38) \quad \text{eq:DEI-SW}$$

If Condition (33) holds, there exists some numerical entropy fluxes $\mathcal{F}_{j \pm \frac{1}{2}}^E$ consistent with the entropy flux of (14) such that the following discrete entropy inequality holds

$$E_j^{n+1} \le E_j^{n+\frac{1}{2}} + \frac{\Delta t^n}{\Delta x} (\mathcal{F}_{j+\frac{1}{2}}^E - \mathcal{F}_{j-\frac{1}{2}}^E).$$

*Proof.* Proposition 2.4 implies that the identity $E = E^{sw} + \sum_{\alpha=1}^N E_{\alpha}^c$ holds at the discrete level. The summation of the inequalities (38) on all the sub timestep gives

$$E_j^{sw,n+1} := E_j^{sw,n+1,K} \le E_j^{sw,n+\frac{1}{2}} + \frac{\Delta t^n}{\Delta x} (\mathcal{F}_{j+\frac{1}{2}}^{sw} - \mathcal{F}_{j-\frac{1}{2}}^{sw}) \quad \text{where} \quad \Delta t^n \mathcal{F}_{j+\frac{1}{2}}^{sw} = \sum_{k=0}^{K-1} \delta t^k f_{j+\frac{1}{2}}^{sw,k}.$$
We now prove that

$$ (h\sigma_{\alpha}^2)_j^{n+1} \le (h\sigma_{\alpha}^2)_j^{n+\frac{1}{2}} - \frac{\Delta t^n}{\Delta x} \left( (\sigma_{\alpha,j+\frac{1}{2}}^n)^2 \mathcal{F}_{j+\frac{1}{2}}^h - (\sigma_{\alpha,j-\frac{1}{2}}^n)^2 \mathcal{F}_{j-\frac{1}{2}}^h \right) $$

with the familiar decentered choice

$$ \sigma_{\alpha,j+\frac{1}{2}}^n = \begin{cases} \sigma_{\alpha,j}^n & \text{if } \mathcal{F}_{j+\frac{1}{2}}^h \ge 0 \\ \sigma_{\alpha,j+1}^n & \text{if } \mathcal{F}_{j+\frac{1}{2}}^h < 0 \end{cases} $$

Multiplying by $h_j^{n+1}$ which is nonnegative if Condition (33) holds, the inequality is equivalent to non-negativity of

$$ -\mathcal{F}_{j+\frac{1}{2}}^h h_j^n \left[ \sigma_{\alpha,j+\frac{1}{2}}^{n+\frac{1}{2}} - \sigma_{\alpha,j}^{n+\frac{1}{2}} \right]^2 + \mathcal{F}_{j-\frac{1}{2}}^h h_j^n \left[ \sigma_{\alpha,j-\frac{1}{2}}^{n+\frac{1}{2}} - \sigma_{\alpha,j}^{n+\frac{1}{2}} \right]^2 - \frac{\Delta t^n}{\Delta x} \mathcal{F}_{j+\frac{1}{2}}^h \mathcal{F}_{j-\frac{1}{2}}^h \left[ \sigma_{\alpha,j-\frac{1}{2}}^{n+\frac{1}{2}} - \sigma_{\alpha,j+\frac{1}{2}}^{n+\frac{1}{2}} \right]^2 $$

This quantity is equal to

$$ \begin{cases} \mathcal{F}_{j-\frac{1}{2}}^h \left[ \sigma_{\alpha,j-1}^{n+\frac{1}{2}} - \sigma_{\alpha,j}^{n+\frac{1}{2}} \right]^2 \left[ h_j^n - \frac{\Delta t^n}{\Delta x} \left( \mathcal{F}_{j+\frac{1}{2}}^h \right)^+ \right]^+ & \text{if } \mathcal{F}_{j+\frac{1}{2}}^h \ge 0 \text{ and } \mathcal{F}_{j-\frac{1}{2}}^h \ge 0, \\ -\mathcal{F}_{j+\frac{1}{2}}^h \left[ \sigma_{\alpha,j+1}^{n+\frac{1}{2}} - \sigma_{\alpha,j}^{n+\frac{1}{2}} \right]^2 \left[ h_j^n + \frac{\Delta t^n}{\Delta x} \mathcal{F}_{j-\frac{1}{2}}^h \right]^+ & \text{if } \mathcal{F}_{j+\frac{1}{2}}^h \le 0 \text{ and } \mathcal{F}_{j-\frac{1}{2}}^h \le 0, \\ \text{clearly nonnegative otherwise.} \end{cases} $$

Condition (33) yields the nonnegativity in the first two cases.

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

**Acknowledgements** This work was supported in part by Deutsche Forschungsgemeinschaft under DFG grant 320021702 / GRK2326. The authors would like to thank Emmanuel Audusse and Carlos Parés for fruitful discussions. I addition SH would like thank Sebastian Noelle for a careful reading and discussions of the manuscript.

## References

[AB03] E. Audusse and M.-O. Bristeau. Transport of pollutant in shallow water: a two time steps kinetic method. M2AN Math. Model. Numer. Anal., 37(2):389-416, 2003.

[ABBKPO4] E. Audusse, F. Bouchut, M.-O. Bristeau, R. Klein, and B. Perthame. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM J. Sci. Comput., 25(6):2050-2065, 2004.

[ABF19] S. Allgeyer, M.-O. Bristeau, D. Froger, R. Hamouda, V. Jauzein, A. Mangeney, J. Sainte-Marie, F. Souillé, and M. Vallée. Numerical approximation of the 3D hydrostatic Navier-Stokes system with free surface. ESAIM Math. Model. Numer. Anal., 53(6):1981-2024, 2019.

[ABPSM11b] E. Audusse, M.-O. Bristeau, M. Pelanti, and J. Sainte-Marie. Approximation of the hydrostatic Navier-Stokes system for density stratified flows by a multilayer model: kinetic interpretation and numerical solution. J. Comput. Phys., 230(9):3453-3478, 2011.

[ABPSM11a] E. Audusse, M.-O. Bristeau, B. Perthame, and J. Sainte-Marie. A multilayer Saint-Venant system with mass exchanges for shallow water flows. Derivation and numerical validation. ESAIM Math. Model. Numer. Anal., 45(1):169-200, 2011.

[ABS18] E. Audusse, M.-O. Bristeau, and J. Sainte-Marie. Kinetic entropy for the layer-averaged hydrostatic navier-stokes equations. https://arxiv.org/abs/1710.04054v1, 2018.
[ADDGNP21] [ADD+21] E. Audusse, V. Dubos, A. Duran, N. Gaveau, Y. Nasseri, and Y. Penel. Numerical approximation of the shallow water equations with Coriolis source term. In *CEMRACS 2019—geophysical fluids, gravity flows*, volume 70 of *ESAIM Proc. Surveys*, pages 31–44. EDP Sci., Les Ulis, 2021.

[BBSM2013] [BBSM13] A. C. Boulanger, M.-O. Bristeau, and J. Sainte-Marie. Analytical solutions for the free surface hydrostatic Euler equations. *Commun. Math. Sci.*, 11(4):993–1010, 2013.

[ukacovaNoelle2011] [BNLM11] A. Bollermann, S. Noelle, and M. Lukáčová-Medvid'ová. Finite volume evolution Galerkin methods for the shallow water equations with dry beds. *Commun. Comput. Phys.*, 10(2):371–404, 2011.

[B63] [Bry63] Kirk Bryan. A numerical investigation of a nonlinear model of a wind-driven ocean. *Journal of Atmospheric Sciences*, 20(6):594–606, 1963.

[B69] [Bry69] K. Bryan. A numerical method for the study of the circulation of the world ocean. *Journal of Computational Physics*, 4(3):347–376, 1969.

[B86] [Bö86] Claus W. Böning. On the influence of frictional parameterization in wind-driven ocean circulation models. *Dynamics of Atmospheres and Oceans*, 10(1):63–92, 1986.

[CDV17] [CDV17] F. Couderc, A. Duran, and J.-P. Vila. An explicit asymptotic preserving low Froude scheme for the multilayer shallow water model with density stratification. *J. Comput. Phys.*, 343:235–270, 2017.

[CN17] [CN17] G. Chen and S. Noelle. A new hydrostatic reconstruction scheme based on subcell reconstructions. *SIAM J. Numer. Anal.*, 55(2):758–784, 2017.

[CN22] [CN22] G. Chen and S. Noelle. A unified surface-gradient and hydrostatic reconstruction scheme for the shallow water equations. *J. Comput. Phys.*, 467:Paper No. 111463, 16, 2022.

[DDMLBE19] [DDM+19] J. Demange, L. Debreu, P. Marchesiello, F. Lemarié, E. Blayo, and C. Eldred. Stability analysis of split-explicit free surface ocean models: implication of the depth-independent barotropic mode approximation. *J. Comput. Phys.*, 398:108875, 26, 2019.

[DL11] [DL16] G. Dujardin and P. Lafitte. Asymptotic behaviour of splitting schemes involving time-subcycling techniques. *IMA J. Numer. Anal.*, 36(4):1804–1841, 2016.

[DMEHGGSM25] [DMEHG+25] B. Di Martino, C. El Hassanieh, E. Godlewski, J. Guillod, and J. Sainte-Marie. Hyperbolicity of a semi-Lagrangian formulation of the hydrostatic free-surface Euler system. *Nonlinearity*, 38(1):Paper No. 015018, 48, 2025.

[DR94] [DSM93] J. K. Dukowicz, R. D. Smith, and R. C. Malone. A reformulation and implementation of the bryan-cox-semtner ocean model on the connection machine. *Journal of Atmospheric and Oceanic Technology*, 10(2):195–208, 1993.

[GCGP24] [GTCDGR24] V. González-Tabernero, M. J. Castro Díaz, and J. A. García-Rodríguez. High-order well-balanced finite volume schemes for 1d and 2d shallow-water equations with variable bathymetry and Coriolis forces. In *Hyperbolic problems: theory, numerics, applications. Vol. II*, volume 35 of *SEMA SIMAI Springer Ser.*, pages 97–109. Springer, Cham, [2024] ©2024.

[H96] [HB96] Robert L. Higdon and Andrew F. Bennett. Stability analysis of operator splitting for large-scale ocean modeling. *Journal of Computational Physics*, 123(2):311–329, 1996.

[H99] [Hig99] R. L. Higdon. Implementation of a barotropic-baroclinic time splitting for isopycnic coordinate ocean modeling. *J. Comput. Phys.*, 148(2):579–604, 1999.
[KSP91] [KWSP91] P. D. Killworth, D. J. Webb, D. Stainforth, and S. M. Paterson. The development of a free-surface bryan-cox-semtner ocean model. *Journal of Physical Oceanography*, 21(9):1333–1348, 1991.

[LJWG22] [LJW$^+$22] R. Lan, L. Ju, Z. Wang, M. Gunzburger, and P. Jones. High-order multirate explicit time-stepping schemes for the baroclinic-barotropic split dynamics in primitive equations. *J. Comput. Phys.*, 457:Paper No. 111050, 25, 2022.

[SM05] [SM05] A. F. Shchepetkin and J. C. McWilliams. The regional oceanic modeling system (roms): a split-explicit, free-surface, topography-following-coordinate oceanic model. *Ocean Modelling*, 9(4):347–404, 2005.

[Stommel48] [Sto48] H. Stommel. The westward intensification of wind-driven ocean currents. *Eos, Transactions American Geophysical Union*, 29(2):202–206, 1948.

[Vallis [Val17]] G. K. Vallis. *Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation*. Cambridge University Press, 2 edition, 2017.