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

HAL
open science

# High order ImEx method for the shallow water model

Maria Kazolea, Martin Parisot, Ralph Lteif

► To cite this version:

Maria Kazolea, Martin Parisot, Ralph Lteif. High order ImEx method for the shallow water model. XXVIII CEDYA / XVIII CMA 2024 Congress, Jun 2024, Bilbao, Spain. hal-04884462

HAL Id: hal-04884462

https://inria.hal.science/hal-04884462v1

Submitted on 13 Jan 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

XXVIII
CONGRESS OF DIFFERENTIAL
EQUATIONS AND APPLICATIONS

XVIII
CONGRESS OF
APPLIED MATHEMATICS

BILBAO ✠ JUNE 24-28, 2024

# High order ImEx method for the shallow water model

M. Kazolea, R. Lteif, M. Parisot

INRIA, Univ. Bordeaux, UMR 5251, 200 Avenue de la Vieille Tour, 33405 Talence cedex, France

CEDYA 2024

Bilbao - June 24-28 2024

Inria

anr®
Climate change is expected to increase coastal flooding hazard in years to come.

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

La Faute-sur-Mer, western France
after Xynthia storm, February 28
2010

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

Sources : carte de simulation des niveaux d'eau établis par le BRGM (years 2010, BRGM-RP-58261-FR, p26), niveau d'eau extrême mesuré dans le bourg de Agglution selon rapport CSEDD (16 sep 2010, Maison n°00736-01, p21)

Water level simulation map established by BRGM

Recent storm events (Xynthia [2010], Klaus [2009], Lothar [1999]) have shown the very destructive effects that flooding can have in terms of human victims and economic losses.
**Objective:** Design a multi-scale tool capable of simulating different features of oceanic flows, from large scale linear waves offshore to small scale non-linear flows in coastal areas.

<figure><img src="image_6.png" /><figcaption>Zone division</figcaption></figure>

Offshore, tidal waves and atmospheric surges are very linear processes under the influence of Earth rotation.
A model is essential for understanding complex flow dynamics.

Departing from the incompressible Navier-Stokes equations with gravity and
following [Gerbeau, Perthame'00], we get the 2D shallow water system under the
following assumptions:
A model is essential for understanding complex flow dynamics.

Departing from the incompressible Navier-Stokes equations with gravity and
following [Gerbeau, Perthame'00], we get the 2D shallow water system under the
following assumptions:

-> shallow water assumption (characteristic depth << wavelength)
A model is essential for understanding complex flow dynamics.

Departing from the incompressible Navier-Stokes equations with gravity and
following [Gerbeau, Perthame'00], we get the 2D shallow water system under the
following assumptions:

→ shallow water assumption (characteristic depth ≪ wavelength)

→ vertically averaged horizontal velocity (velocity profile varies slowly vertically)
A model is essential for understanding complex flow dynamics.

Departing from the incompressible Navier-Stokes equations with gravity and
following [Gerbeau, Perthame'00], we get the 2D shallow water system under the
following assumptions:

→ shallow water assumption (characteristic depth ≪ wavelength)

-> vertically averaged horizontal velocity (velocity profile varies slowly vertically)

-> hydrostatic pressure ($P_{bot} = P_{atm} + g \times$ water column weight)
A model is essential for understanding complex flow dynamics.

Departing from the incompressible Navier-Stokes equations with gravity and following [Gerbeau, Perthame'00], we get the 2D shallow water system under the following assumptions:

→ shallow water assumption (characteristic depth ≪ wavelength)

→ vertically averaged horizontal velocity (velocity profile varies slowly vertically)

→ hydrostatic pressure ($P_{bot} = P_{atm} + g \times$ water column weight)

$$ \left\{ \begin{array}{l} \partial_t h + \nabla \cdot (hU) = 0, \\ \\ \partial_t (hU) + \nabla \cdot (hU \otimes U) + \nabla \left( \frac{g}{2} h^2 \right) = -gh \nabla z. \end{array} \right. \qquad (\text{SW}) $$

<figure><img src="image_7.png" /></figure>
Consider the following rescaling:

$$ \tilde{x} = \frac{\mathbf{x}}{l^*}, \quad \tilde{h} = \frac{h}{h^*}, \quad \tilde{U} = \frac{U}{U^*}, \quad \tilde{t} = \frac{t}{t^*}, \quad \tilde{z} = \frac{z}{h^*} $$

System (SW) becomes (after dropping the tildes):

$$ \left\{ \begin{array}{l} \partial_t h + \nabla \cdot (hU) = 0, \\ \\ \partial_t(hU) + \nabla \cdot (hU \otimes U) + \frac{1}{\text{Fr}^2} \nabla \left(\frac{h^2}{2}\right) = -\frac{h}{\text{Fr}^2} \nabla z, \end{array} \right. \quad (\text{SW}_{\text{Fr}}) $$

with the characteristic Froude number

$$ \text{Fr} := \frac{U^*}{\sqrt{gh^*}} = \frac{\text{flow velocity}}{\text{velocity of gravity waves}}. $$
Consider the following rescaling:

$$ \bar{x} = \frac{\mathbf{x}}{l^*}, \quad \bar{h} = \frac{h}{h^*}, \quad \bar{U} = \frac{U}{U^*}, \quad \bar{t} = \frac{t}{t^*}, \quad \bar{z} = \frac{z}{h^*} $$

System (SW) becomes (after dropping the tildes):

$$ \left\{ \begin{array}{l} \partial_t h + \nabla \cdot (hU) = 0, \\ \partial_t(hU) + \nabla \cdot (hU \otimes U) + \frac{1}{\text{Fr}^2} \nabla \left( \frac{h^2}{2} \right) = - \frac{h}{\text{Fr}^2} \nabla z, \end{array} \right. \qquad (\text{SW}_{\text{Fr}}) $$

with the characteristic Froude number

$$ \text{Fr} := \frac{U^*}{\sqrt{gh^*}} = \frac{\text{flow velocity}}{\text{velocity of gravity waves}}. $$

* We focus in particular on the variation of the Froude number Fr that goes from 1 at the coastline to two or three orders less offshore.

$$ \left. \begin{array}{l} \text{mean depth } h^* > 3000m \\ \text{current speed } U^* \approx 1 \text{ m/s} \end{array} \right\} \Rightarrow \text{Fr} \approx 10^{-3} $$
Consider the following rescaling:

$$ \bar{x} = \frac{\mathbf{x}}{l^*}, \quad \bar{h} = \frac{h}{h^*}, \quad \bar{U} = \frac{U}{U^*}, \quad \bar{t} = \frac{t}{t^*}, \quad \bar{z} = \frac{z}{h^*} $$

System (SW) becomes (after dropping the tildes):

$$ \left\{ \begin{array}{l} \partial_t h + \nabla \cdot (hU) = 0, \\ \\ \partial_t(hU) + \nabla \cdot (hU \otimes U) + \frac{1}{\text{Fr}^2} \nabla \left( \frac{h^2}{2} \right) = - \frac{h}{\text{Fr}^2} \nabla z, \end{array} \right. \qquad (\text{SW}_{\text{Fr}}) $$

with the characteristic Froude number

$$ \text{Fr} := \frac{U^*}{\sqrt{gh^*}} = \frac{\text{flow velocity}}{\text{velocity of gravity waves}}. $$

* We focus in particular on the variation of the Froude number Fr that goes from 1 at the coastline to two or three orders less offshore.

$$ \left. \begin{array}{l} \text{mean depth } h^* > 3000m \\ \text{current speed } U^* \approx 1 \text{ m/s} \end{array} \right\} \Rightarrow \text{Fr} \approx 10^{-3} $$

* Limiting system in the limit $Fr \to 0$ satisfy the “lake equations”.
The system ($SW_{Fr}$) is a hyperbolic system of conservation laws:

$$ \frac{\partial U}{\partial t} + \nabla \cdot F(U) = S(U, z) $$

$$ U = \begin{pmatrix} h \\ hU \end{pmatrix}, F(U) = \begin{pmatrix} hU^T \\ hU \otimes U + h^2 \mathbf{l}_2/(2Fr^2) \end{pmatrix}, S(U, z) = \begin{pmatrix} 0 \\ -\frac{h}{Fr^2} \nabla z \end{pmatrix} $$

Its eigenvalues in the direction $n$ are $\lambda_j = (U \cdot n) + j \frac{c}{Fr}, j \in \{-1, 0, 1\}$ with $c = \sqrt{h}$ being the scaled speed of sound.
The system (SW<sub>Fr</sub>) is a hyperbolic system of conservation laws:

$$ \frac{\partial U}{\partial t} + \nabla \cdot F(U) = S(U, z) $$

$$ U = \begin{pmatrix} h \\ hU \end{pmatrix}, F(U) = \begin{pmatrix} hU^T \\ hU \otimes U + h^2 \mathbf{I}_2 / (2\mathrm{Fr}^2) \end{pmatrix}, S(U, z) = \begin{pmatrix} 0 \\ -\frac{h}{\mathrm{Fr}^2} \nabla z \end{pmatrix} $$

Its eigenvalues in the direction **n** are $\lambda_j = (U \cdot \mathbf{n}) + j \frac{c}{\mathrm{Fr}}$, $j \in \{-1, 0, 1\}$ with $c = \sqrt{h}$ being the scaled speed of sound.

Conventional explicit numerical methods are inefficient and often impractical.

**Drawback:** The time step constraint of an explicit time discretization satisfies

$$ \Delta t = \mathrm{CFL} \min \left( \frac{\Delta x}{\max(|u| + c/\mathrm{Fr})}, \frac{\Delta y}{\max(|v| + c/\mathrm{Fr})} \right) \sim \mathrm{Fr} \Delta \min $$
The system ($SW_{Fr}$) is a hyperbolic system of conservation laws:

$$ \frac{\partial U}{\partial t} + \nabla \cdot F(U) = S(U, z) $$

$$ U = \left( \frac{h}{hU} \right), F(U) = \left( hU \otimes U + h^2 \mathbf{l}_2 / (2\mathrm{Fr}^2) \right), S(U, z) = \left( \begin{array}{c} 0 \\ -\frac{h}{\mathrm{Fr}^2} \nabla z \end{array} \right) $$

Its eigenvalues in the direction $n$ are $\lambda_j = (U \cdot n) + j \frac{c}{\mathrm{Fr}}$, $j \in \{-1, 0, 1\}$ with $c = \sqrt{h}$ being the scaled speed of sound.

Conventional explicit numerical methods are inefficient and often impractical.

**Drawback:** The time step constraint of an explicit time discretization satisfies

$$ \Delta t = \mathrm{CFL} \min \left( \frac{\Delta x}{\max(|u| + c/\mathrm{Fr})}, \frac{\Delta y}{\max(|v| + c/\mathrm{Fr})} \right) \sim \mathrm{Fr} \Delta \min $$

Issues related to explicit methods in low Froude number regime:

D1 excessive numerical diffusion on the coarser elements ⇒ damp the water level $\phi$.
The system ($SW_{Fr}$) is a hyperbolic system of conservation laws:

$$ \frac{\partial U}{\partial t} + \nabla \cdot F(U) = S(U, z) $$

$$ U = \left( \frac{h}{hU} \right), F(U) = \left( hU \otimes U + h^2 \mathbf{l}_2 / (2Fr^2) \right), S(U, z) = \left( \begin{array}{c} 0 \\ -\frac{h}{Fr^2} \nabla z \end{array} \right) $$

Its eigenvalues in the direction $n$ are $\lambda_j = (U \cdot n) + j \frac{c}{Fr}, j \in \{-1, 0, 1\}$ with $c = \sqrt{h}$ being the scaled speed of sound.

Conventional explicit numerical methods are inefficient and often impractical.

**Drawback:** The time step constraint of an explicit time discretization satisfies

$$ \Delta t = \text{CFL min} \left( \frac{\Delta x}{\max(|u| + c/Fr)}, \frac{\Delta y}{\max(|v| + c/Fr)} \right) \sim \text{Fr}\Delta\text{min} $$

Issues related to explicit methods in low Froude number regime:

D1 excessive numerical diffusion on the coarser elements ⇒ damp the water level $\Phi$.

D2 not asymptotically consistent. Divergence between continuous system and its discretization ⇒ predominant large numerical dissipation [Guillard, Murrone'04].
The system ($SW_{Fr}$) is a hyperbolic system of conservation laws:

$$ \frac{\partial U}{\partial t} + \nabla \cdot F(U) = S(U, z) $$

$$ U = \begin{pmatrix} h \\ hU \end{pmatrix}, F(U) = \begin{pmatrix} hU^T \\ hU \otimes U + h^2 \mathbf{l}_2/(2Fr^2) \end{pmatrix}, S(U, z) = \begin{pmatrix} 0 \\ -\frac{h}{Fr^2} \nabla z \end{pmatrix} $$

Its eigenvalues in the direction $n$ are $\lambda_j = (U \cdot n) + j \frac{c}{Fr}, j \in \{-1, 0, 1\}$ with $c = \sqrt{h}$ being the scaled speed of sound.

Conventional explicit numerical methods are inefficient and often impractical.

**Drawback:** The time step constraint of an explicit time discretization satisfies

$$ \Delta t = \text{CFL min} \left( \frac{\Delta x}{\max(|u| + c/Fr)}, \frac{\Delta y}{\max(|v| + c/Fr)} \right) \sim \text{Fr}\Delta\text{min} $$

Issues related to explicit methods in low Froude number regime:

D1 excessive numerical diffusion on the coarser elements ⇒ damp the water level $\Phi$.

D2 not asymptotically consistent. Divergence between continuous system and its discretization ⇒ predominant large numerical dissipation [Guillard, Murrone'04]. $\Phi$

D3 high gravity wave speed require small time steps ⇒ does not allow long time simulations (the simulation of tide-surge), even with high performance computing. $\Phi$
Specialized numerical schemes for low-Froude number flows found in literature:

* The first family of schemes respond to $D2 \Rightarrow$ asymptotic consistency when $Fr$ tends to zero . [Bruel et al.'18, Dellacherie et al.'16]
Specialized numerical schemes for low-Froude number flows found in literature:

* The first family of schemes respond to $D2 \Rightarrow$ asymptotic consistency when Fr tends to zero $\approx$. [Bruel et al.'18, Dellacherie et al.'16]

* The second family of schemes respond to $D3 \Rightarrow$ long time simulation, CFL independent of Fr $\approx$. [Degond, Tang'11, Haack et al.'12, Herbin et al.'14, Parisot, Vila'16]
Specialized numerical schemes for low-Froude number flows found in literature:

* The first family of schemes respond to $D2 \Rightarrow$ asymptotic consistency when Fr tends to zero $\approx$. [Bruel et al.'18, Dellacherie et al.'16]

* The second family of schemes respond to $D3 \Rightarrow$ long time simulation, CFL independent of Fr $\approx$. [Degond, Tang'11, Haack et al.'12, Herbin et al.'14, Parisot, Vila'16]

Stability requires implicit treatment (gravity waves) $\Rightarrow$ schemes based on Implicit Explicit (ImEx) or fully implicit.
Specialized numerical schemes for low-Froude number flows found in literature:

* The first family of schemes respond to $D2 \Rightarrow$ asymptotic consistency when Fr tends to zero $\approx$. [Bruel et al.'18, Dellacherie et al.'16]

* The second family of schemes respond to $D3 \Rightarrow$ long time simulation, CFL independent of Fr $\approx$. [Degond, Tang'11, Haack et al.'12, Herbin et al.'14, Parisot, Vila'16]

Stability requires implicit treatment (gravity waves) $\Rightarrow$ schemes based on Implicit Explicit (ImEx) or fully implicit.

**Fully implicit time integration** is too costly since one needs to solve nonlinear systems.
Specialized numerical schemes for low-Froude number flows found in literature:

* The first family of schemes respond to $D2 \Rightarrow$ asymptotic consistency when Fr tends to zero ◊ [Bruel et al.'18, Dellacherie et al.'16]

* The second family of schemes respond to $D3 \Rightarrow$ long time simulation, CFL independent of Fr ◊ [Degond, Tang'11, Haack et al.'12, Herbin et al.'14, Parisot, Vila'16]

Stability requires implicit treatment (gravity waves) ⇒ schemes based on Implicit Explicit (ImEx) or fully implicit.

**Fully implicit time integration** is too costly since one needs to solve nonlinear systems.

The CPR (centered-potential regularization) ImEx scheme [Parisot, Vila'16] is a good first order candidate.

* A fully diagonal segregated method.

* Avoid resolution of large linear systems and limits the number of linear systems to be solved.
**Numerical strategy:** Cell-centered FV method on uniform Cartesian grid

in cell:

$$ \psi_k = \frac{1}{|V_k|} \int_{V_k} \psi dx $$

at the face:

$$ 2(\psi)_f = \psi_k + \psi_{k_f} \quad \text{and} $$

$$ 2[\psi]_k^{k_f} = \psi_{k_f} - \psi_k $$

parameters:

$$ l_k = \frac{|V_k|}{|\partial V_k|} \quad \text{and} $$

$$ \mu_f^k = \frac{|f|}{|\partial V_k|} $$

<figure><img src="image_8.png" /></figure>
**Numerical strategy:** Cell-centered FV method on uniform Cartesian grid

$$
\begin{array}{l@{\quad}l@{\quad}l}
\text{in cell:} & \psi_k = \frac{1}{|V_k|} \int_{V_k} \psi \, dx & \\
\text{at the face:} & 2(\psi)_f = \psi_k + \psi_{k_f} \quad \text{and} \quad 2[\psi]_k^k = \psi_{k_f} - \psi_k & \\
\text{parameters:} & l_k = \frac{|V_k|}{|\partial V_k|} \quad \text{and} \quad \mu_t^k = \frac{|f|}{|\partial V_k|} &
\end{array}
$$

**Step 1:**

→ Use an **AUSM** (Advection Upstream Splitting Method) scheme [Liou, Steffen'93]

Splitting of the equation into the advection part and the potential forces

$$
\partial_t \begin{pmatrix} h_k \\ h_k U_k \end{pmatrix} + \frac{1}{l_k} \sum_{f \in \mathbb{F}_k} \begin{pmatrix} \mathcal{F}_f^h \\ \mathcal{F}_f^{hu} \end{pmatrix} \cdot N_f^k \mu_t^k = \begin{pmatrix} 0 \\ \mathcal{D}_k \end{pmatrix}
$$

with $\mathcal{F}_f^h = \int_f hU \, d\sigma$, $\mathcal{F}_f^{hU} = \int_f hU \otimes U \, d\sigma$ and $\mathcal{D}_k = -\frac{1}{|V_k|} \int_{V_k} h \frac{\nabla \phi}{Fr^2} dx$, where
$\phi = h + z$.
**Numerical strategy:** Cell-centered FV method on uniform Cartesian grid

$$
\begin{array}{lll}
\text{in cell:} & \psi_k = \frac{1}{|V_k|} \int_{V_k} \psi dx & \\
\text{at the face:} & 2(\psi)_f = \psi_k + \psi_{k_f} \quad \text{and} \quad 2[\psi]_k^{k_f} = \psi_{k_f} - \psi_k & \\
\text{parameters:} & l_k = \frac{|V_k|}{|\partial V_k|} \quad \text{and} \quad \mu_t^k = \frac{|f|}{|\partial V_k|} &
\end{array}
$$

**Step 1:**

→ Use an **AUSM** (Advection Upstream Splitting Method) scheme [Liou, Steffen'93]

Splitting of the equation into the advection part and the potential forces

$$
\partial_t \begin{pmatrix} h_k \\ h_k U_k \end{pmatrix} + \frac{1}{l_k} \sum_{f \in \mathbb{F}_k} \begin{pmatrix} \mathcal{F}_f^h \\ \mathcal{F}_f^{hu} \end{pmatrix} \cdot N_f^k \mu_t^k = \begin{pmatrix} 0 \\ \mathcal{D}_k \end{pmatrix}
$$

with $\mathcal{F}_f^h = \int_f hU \, d\sigma$, $\mathcal{F}_f^{hU} = \int_f hU \otimes U \, d\sigma$ and $\mathcal{D}_k = -\frac{1}{|V_k|} \int_{V_k} h \nabla^2 \phi \, dx$, where
$\phi = h + z$.

→ Advect the velocity with an **up-wind** scheme

$$
\partial_t(h_k U_k) + \frac{1}{l_k} \sum_{f \in \mathbb{F}_k} \left( U_k (\mathcal{F}_f^h \cdot N_f^k)_+ - U_{k_f} (\mathcal{F}_f^h \cdot N_f^k)_- \right) \mu_f^k = \mathcal{D}_k
$$
**Numerical strategy:** Cell-centered FV method on uniform Cartesian grid

$$
\begin{array}{l@{\quad}l@{\quad}l}
\text{in cell:} & \psi_k = \frac{1}{|V_k|} \int_{V_k} \psi \, dx & \\
\text{at the face:} & 2(\psi)_f = \psi_k + \psi_{k_f} \quad \text{and} \quad 2[\psi]_k^{k_f} = \psi_{k_f} - \psi_k & \\
\text{parameters:} & l_k = \frac{|V_k|}{|\partial V_k|} \quad \text{and} \quad \mu_t^k = \frac{|f|}{|\partial V_k|} &
\end{array}
$$

**Step 1:**

→ Use an **AUSM** (Advection Upstream Splitting Method) scheme [Liou, Steffen'93]

Splitting of the equation into the advection part and the potential forces

$$
\partial_t \begin{pmatrix} h_k \\ h_k U_k \end{pmatrix} + \frac{1}{l_k} \sum_{f \in \mathbb{F}_k} \begin{pmatrix} \mathcal{F}_f^h \\ \mathcal{F}_f^{hu} \end{pmatrix} \cdot N_f^k \mu_t^k = \begin{pmatrix} 0 \\ \mathcal{D}_k \end{pmatrix}
$$

with $\mathcal{F}_f^h = \int_f hU \, d\sigma$, $\mathcal{F}_f^{hU} = \int_f hU \otimes U \, d\sigma$ and $\mathcal{D}_k = -\frac{1}{|V_k|} \int_{V_k} h \nabla^2 \phi \, dx$, where
$\phi = h + z$.

→ Advect the velocity with an **up-wind** scheme

$$
\partial_t(h_k U_k) + \frac{1}{l_k} \sum_{f \in \mathcal{F}_k} \left( U_k (\mathcal{F}_f^h \cdot N_f^k)_+ - U_{k_f} (\mathcal{F}_f^h \cdot N_f^k)_- \right) \mu_f^k = \mathcal{D}_k
$$

● Ensure the dissipation of the discrete kinetic energy.
Step 2:

→ Use a **centered discretization** of the potential for any $Fr > 0$. [Dellacherie'10]

$$D_k \approx -\frac{h_k}{l_k} \sum_{f \in \mathbb{F}_k} \frac{[\phi]_k^{k_f}}{Fr^2} N_f^k |f|$$
Step 2:

→ Use a **centered discretization** of the potential for any $Fr > 0$. [Dellacherie'10]

$$D_k \approx -\frac{h_k}{l_k} \sum_{f \in \mathbb{F}_k} \frac{[\phi]_k^{k_f}}{Fr^2} N_f^k |f|$$

* Leads to a **consistent** numerical scheme when $Fr$ goes to 0.
### Step 2:

→ Use a **centered discretization** of the potential for any $Fr > 0$. [Dellacherie'10]

$$ \mathcal{D}_k \approx -\frac{h_k}{l_k} \sum_{f \in \mathbb{F}_k} \frac{[\phi]_k^{k_f}}{Fr^2} N_f^k |f| $$

* Leads to a **consistent** numerical scheme when $Fr$ goes to 0.

### Step 3:

→ **Regularization** using the potential jump. [Grenier et al.'13, Parisot and Vila'15]

$$ \mathcal{F}_f^h = (hU)_f - \gamma \Delta t \left(\frac{h}{l}\right)_f \frac{[\phi]_k^{k_f}}{Fr^2} N_f^k, \quad (\Delta t : \text{time scale}, \quad \gamma : \text{regu. param.}) $$
### Step 2:

→ Use a **centered discretization** of the potential for any $Fr > 0$. [Dellacherie'10]

$$ \mathcal{D}_k \approx -\frac{h_k}{l_k} \sum_{f \in \mathbb{F}_k} \frac{[\phi]_k^{k_f}}{Fr^2} N_f^k |f| $$

* Leads to a **consistent** numerical scheme when $Fr$ goes to 0.

### Step 3:

→ **Regularization** using the potential jump. [Grenier et al.'13, Parisot and Vila'15]

$$ \mathcal{F}_f^h = (hU)_f - \gamma \Delta t \left(\frac{h}{l}\right)_f \frac{[\phi]_k^{k_f}}{Fr^2} N_f^k, \quad (\Delta t : \text{time scale}, \quad \gamma : \text{regu. param.}) $$

Ensure:

* The stability of the **steady state at rest** ($\phi = Cst$ and $u = 0$).
### Step 2:

→ Use a **centered discretization** of the potential for any $Fr > 0$. [Dellacherie'10]

$$D_k \approx -\frac{h_k}{l_k} \sum_{f \in \mathbb{F}_k} \frac{[\phi]_k^{k_f}}{Fr^2} N_f^k |f|$$

* Leads to a **consistent** numerical scheme when $Fr$ goes to 0.

### Step 3:

→ **Regularization** using the potential jump. [Grenier et al.'13, Parisot and Vila'15]

$$\mathcal{F}_f^h = (hU)_f - \gamma \Delta t \left(\frac{h}{l}\right)_f \frac{[\phi]_k^{k_f}}{Fr^2} N_f^k, \quad (\Delta t : \text{time scale}, \quad \gamma : \text{regu. param.})$$

Ensure:

* The stability of the **steady state at rest** ($\phi = Cst$ and $u = 0$).

* The **dissipation of the discrete potential energy**.
Q) How to get a stable numerical scheme independently of Fr?
Q) How to get a stable numerical scheme independently of $Fr$?

A) Using an ImEx scheme: **Implicit** for the water level $h$ / **Explicit** for the velocity $u$
Q) How to get a stable numerical scheme independently of Fr?

A) Using an ImEx scheme: **Implicit** for the water level *h* / **Explicit** for the velocity *u*

**Step 1:**

→ **water level: implicit** scheme of type **non-linear** advection-diffusion

$$h_k^{n+1} - h_k^n + \frac{\Delta t}{l_k} \sum_{f \in \mathbb{F}_k} \left( \left( (h^{n+1} U^n)_f \cdot N_f^k - \gamma \Delta t \left( \frac{h^{n+1}}{l} \right)_f \cdot \frac{[\phi^{n+1}]_k^{k_f}}{Fr^2} \right) \mu_f^k = 0 \right. \\ \left. \qquad \mathcal{F}_f^{n+1} \cdot N_f^k \right)$$

Solving a **non-linear** system by using the Newton-Raphson method, transforming the nonlinear system into a series of linear systems.
Q) How to get a stable numerical scheme independently of Fr?

A) Using an ImEx scheme: **Implicit** for the water level *h* / **Explicit** for the velocity *u*

**Step 1:**

→ water level: **implicit** scheme of type **non-linear** advection-diffusion

$$h_k^{n+1} - h_k^n + \frac{\Delta t}{l_k} \sum_{f \in \mathbb{F}_k} \underbrace{\left( (h^{n+1} U^n)_f \cdot N_f^k - \gamma \Delta t \left( \frac{h^{n+1}}{l} \right)_f \cdot \frac{[\phi^{n+1}]_k^{k_f}}{\text{Fr}^2} \right)}_{\mathcal{F}_f^{n+1}} \mu_f^k = 0$$

Solving a **non-linear** system by using the Newton-Raphson method, transforming the nonlinear system into a series of linear systems.

**Step 2:**

→ velocity: **explicit** up-wind scheme with source term.

$$h_k^{n+1} U_k^{n+1} - h_k^n U_k^n + \frac{\Delta t}{l_k} \sum_{f \in \mathbb{F}_k} \left( U_k^n \left( \bar{\mathcal{F}}_f^{n+1} \cdot N_f^k \right)_+ - U_{k_f}^n \left( \bar{\mathcal{F}}_f^{n+1} \cdot N_f^k \right)_- \right) \mu_f^k \\ = - \frac{\Delta t}{l_k} h_k^{n+1} \sum_{f \in \mathbb{F}_k} \frac{[\phi^{n+1}]_k^{k_f}}{\text{Fr}^2} N_f^k \mu_f^k$$
Entropy dissipation [Parisot and Vila '16]

Let $\gamma \ge 1$ and assume the following CFL-like condition is satisfied

$$ \left( \left| (h^{n+1} U^n)_t \cdot N_r^k \right| + (h_k^{n+1})_r \sqrt{\frac{\gamma \left| \left[ \phi^{n+1} \right)_k^{k_f} \right|}{Fr^2}} \right) \Delta t \le \frac{\min(l_k, l_{k_f}) \min(h_k^{n+1}, h_{k_f}^{n+1})}{2} $$

then the discrete mechanic energy is **decreasing**.
The efficient high order time integration is more challenging in the context of ImEx schemes.

Several Runge-Kutta schemes can be found in the literature:

* [Dimarco et al.'17, Boscheri et al.'20] second order schemes based on several multi-stage or step ImEx formulations.

* [Dimarco et al.'18, Michel-Dansac, Thomann'21] higher order TVD ImEx schemes.

Large number of linear systems to solve.
Q) How to extend the CPR ImEx scheme to second order in time while limiting the number of linear systems to solve?
Q) How to extend the CPR ImEx scheme to second order in time while limiting the number of linear systems to solve?

A) We focus on **Crank Nicolson** schemes.
Q) How to extend the CPR ImEx scheme to second order in time while limiting the number of linear systems to solve?

A) We focus on **Crank Nicolson** schemes.

**Step 1:**

→ water level: **implicit** + **extrapolation**

$$h_k^{n+1} - h_k^n + \frac{\Delta t}{l_k} \sum_{t \in \Pi_k} \left( (h^n U^n)_t \cdot N_t^k - \frac{\Delta t}{2} \left( \frac{h^{n+\frac{1}{2}}}{2l} \right)_t \frac{[\phi^{n+\frac{1}{2}}]_k}{Fr^2} - \frac{\Delta t}{2} \left( \nabla(\bar{h} | \bar{U} \cdot N_t^k|^2)^{n+\frac{1}{2}} \cdot N_t^k \right)_t - \frac{\Delta t}{2} \left( \nabla^{\perp}(\tilde{h} \tilde{u} v)^{n+\frac{1}{2}} \cdot N_t^k \right)_t \right) \mu_t^k = 0$$
Q) How to extend the CPR ImEx scheme to second order in time while limiting the number of linear systems to solve?

A) We focus on **Crank Nicolson** schemes.

Step 1:

→ water level: **implicit** + extrapolation

$$h_k^{n+1} - h_k^n + \frac{\Delta t}{\ell_k} \sum_{t \in \Pi_k} \left( (h^n U^n)_t \cdot N_f^k - \frac{\Delta t}{2} \left( \frac{h^{n+\frac{1}{2}}}{2\ell} \right)_f \frac{|\phi^{n+\frac{1}{2}}|}{Fr^2} N_f^k - \frac{\Delta t}{2} \left( \nabla (\bar{h} |\bar{U} \cdot N_f^k|^2)^{n+\frac{1}{2}} \cdot N_f^k \right)_f - \frac{\Delta t}{2} \left( \nabla^{\perp} (\tilde{h} \tilde{u} v)^{n+\frac{1}{2}} \cdot N_f^k \right)_f \right) \mu_f^k = 0$$

$$ \text{with } (\nabla \phi \cdot N_f^k)_f = \frac{\phi_{f_k} - \phi_k}{\Delta \bar{\phi} \cdot N_f^k} \quad \text{and} \quad (\nabla^\perp \phi \cdot N_f^k)_f = \frac{\phi_{v_1} - \phi_{v_2}}{|f|} $$

$$ \psi^{n+1/2} := \frac{\psi^n + \psi^{n+1}}{2}, \quad \tilde{\psi}^{n+1/2} := \frac{\psi^n + \tilde{\psi}^{n+1}}{2} \quad \text{and} \quad \tilde{\psi}^{n+1} := 2\psi^n - \psi^{n-1}. $$
Step 2:

→ velocity: explicit + extrapolation, up-wind scheme with source term.

$$h_k^{n+1} U_k^{n+1} - h_k^n U_k^n + \frac{\Delta t}{l_k} \sum_{f \in \mathbb{F}_k} \left( \tilde{U}_f^{n+\frac{1}{2},+} \cdot \left( \mathcal{F}_f^{n+\frac{1}{2}} \cdot N_f^k \right)_+ + \tilde{U}_f^{n+\frac{1}{2},-} \cdot \left( \mathcal{F}_f^{n+\frac{1}{2}} \cdot N_f^k \right)_- \right) \mu_f^k \\ = - \frac{\Delta t}{l_k} h_k^{n+\frac{1}{2}} \sum_{f \in \mathbb{F}_k} \frac{[\phi^{n+\frac{1}{2}}]_k^{k_f}}{F_f^2} N_f^k \mu_f^k$$
**Step 2:**

→ velocity: explicit + extrapolation, up-wind scheme with source term.

$$h_k^{n+1} U_k^{n+1} - h_k^n U_k^n + \frac{\Delta t}{l_k} \sum_{f \in \mathbb{F}_k} \left( \tilde{U}_f^{n+\frac{1}{2}} + \left( \mathcal{F}_f^{n+\frac{1}{2}} \cdot N_f^k \right)_+ - \tilde{U}_f^{n+\frac{1}{2}} - \left( \mathcal{F}_f^{n+\frac{1}{2}} \cdot N_f^k \right)_- \right) \mu_f^k \\ = - \frac{\Delta t}{l_k} h_k^{n+\frac{1}{2}} \sum_{f \in \mathbb{F}_k} \frac{[\phi^{n+\frac{1}{2}}]_k^{k_f}}{\mathrm{Fr}^2} N_f^k \mu_f^k$$

**Order 2 in space: MUSCL reconstruction of the velocity (with an unlimited slope).**

* Edge *f* with neighboring cells *k* and *k*<sub>*f*</sub>

$$\nabla_{ij} U = \left( \frac{U_{i+1,j} - U_{i,j}}{\Delta x}, \frac{U_{i,j+1} - U_{i,j}}{\Delta y} \right), \quad \delta_f U = \frac{U_{k_f} - U_k}{\mathrm{dist}(k, k_f)}$$

$$U_f^+ = U_k + \mathrm{dist}(k, f) \times \frac{1}{2} \left( \nabla_k U \cdot N_f^k + \delta_f' U \right)$$

$$U_f^- = U_{k_f} - \mathrm{dist}(f, k_f) \times \frac{1}{2} \left( \nabla_{k_f} U \cdot N_f^k + \delta_f' U \right)$$

$$\begin{array}{c|c|c} & f' & f \quad \begin{array}{r} N_f^k \\ \rightarrow \end{array} \\ & \multicolumn{2}{c}{k \quad U_f^+ \quad U_f^-} \end{array}$$
Transcritical flow

Initial data: $h(x, 0) + z(x) = 0.33$, $h(x, 0)u(x, 0) = 0$.

Boundary conditions: $h(0, t)u(0, t) = 0.18$, $h(25, t) = 0.33$.

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

* Good agreement with exact solution. ImEx 2 provides a better solution compared to ImEx CPR. This is clear in the flow discharge figure.
Subcritical flow over a bump

Initial data: $h(x, 0) + z(x) = 2$, $h(x, 0)u(x, 0) = 0$.

Boundary conditions: $h(0, t)u(0, t) = 4.42$, $h(25, t) = 2$.

<figure><img src="image_10.png" /><figcaption>a) Subcritical flow</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>b) $L^1$ errors for $h$ and $hu$</figcaption></figure>

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

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

Initial condition for the traveling vortex with $Fr = 1$, computed on the 80 × 80 grid.

* Exact solution is available. [Ricchiuto, Bollermann'09]
Traveling vortex

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

Initial condition for the traveling vortex with $Fr = 1$, computed on the 80 × 80 grid.

* Exact solution is available. [Ricchiuto, Bollermann'09]

* $\Omega = [0, 1] \times [0, 1]$. Periodic BC in x, absorbing BC in y.
Traveling vortex

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

Initial condition for the traveling vortex with $Fr = 1$, computed on the 80 × 80 grid.

* Exact solution is available. **[Ricchiuto, Bollermann'09]**

* $\Omega = [0, 1] \times [0, 1]$. Periodic BC in x, absorbing BC in y.

* A rotating vortex initially located at the domain center (0.5, 0.5) is transported repeatedly from left to right with a period $T = 1/6$ by an x-directional uniform flow with the advection velocity $U_{ref} = (6, 0)$.
Order of accuracy

Table: Traveling vortex test: EOC for the ImEx 2 scheme, CFL = 2, T = 1/6

Fr = 1

<table><thead><tr><th>N</th><th>L<sup>1</sup>-error in h</th><th>EOC</th><th>L<sup>1</sup>-error in hu</th><th>EOC</th><th>L<sup>1</sup>-error in hv</th><th>EOC</th></tr></thead><tbody><tr><td>40</td><td>1.262e-02</td><td></td><td>2.497e-01</td><td></td><td>4.027e-01</td><td></td></tr><tr><td>80</td><td>3.0520e-03</td><td><strong>2.0474</strong></td><td>7.070e-02</td><td><strong>1.8206</strong></td><td>1.259e-01</td><td><strong>1.6771</strong></td></tr><tr><td>160</td><td>7.002e-04</td><td><strong>2.1239</strong></td><td>1.827e-02</td><td><strong>1.9515</strong></td><td>3.507e-02</td><td><strong>1.8442</strong></td></tr><tr><td>320</td><td>1.677e-04</td><td><strong>2.0614</strong></td><td>4.830e-03</td><td><strong>1.9201</strong></td><td>9.491e-03</td><td><strong>1.8857</strong></td></tr></tbody></table>

Fr = 0.1

<table><thead><tr><th>N</th><th>L<sup>1</sup>-error in h</th><th>EOC</th><th>L<sup>1</sup>-error in hu</th><th>EOC</th><th>L<sup>1</sup>-error in hv</th><th>EOC</th></tr></thead><tbody><tr><td>40</td><td>1.106e-04</td><td></td><td>2.614e-01</td><td></td><td>4.276e-01</td><td></td></tr><tr><td>80</td><td>2.988e-05</td><td><strong>1.8889</strong></td><td>7.251e-02</td><td><strong>1.8501</strong></td><td>1.324e-01</td><td><strong>1.6906</strong></td></tr><tr><td>160</td><td>7.283e-06</td><td><strong>2.0367</strong></td><td>1.882e-02</td><td><strong>1.9457</strong></td><td>3.658e-02</td><td><strong>1.8565</strong></td></tr><tr><td>320</td><td>1.808e-06</td><td><strong>2.0098</strong></td><td>4.964e-03</td><td><strong>1.9227</strong></td><td>9.859e-03</td><td><strong>1.8918</strong></td></tr></tbody></table>
Order of accuracy:

<table><caption>Table: Traveling vortex test: EOC for the ImEx 2 scheme, CFL = 2, T = 1/6</caption><thead><tr><th>N</th><th>L<sup>1</sup>-error in h</th><th>EOC</th><th>L<sup>1</sup>-error in h<sub>u</sub></th><th>EOC</th><th>L<sup>1</sup>-error in h<sub>v</sub></th><th>EOC</th></tr></thead><tbody><tr><td>40</td><td>1.445e-06</td><td></td><td>2.617e-01</td><td></td><td>4.274e-01</td><td></td></tr><tr><td>80</td><td>3.311e-07</td><td><strong>2.1258</strong></td><td>7.255e-02</td><td><strong>1.8509</strong></td><td>1.324e-01</td><td><strong>1.6901</strong></td></tr><tr><td>160</td><td>7.295e-08</td><td><strong>2.1825</strong></td><td>1.884e-02</td><td><strong>1.9450</strong></td><td>3.659e-02</td><td><strong>1.8560</strong></td></tr><tr><td>320</td><td>1.756e-08</td><td><strong>2.0547</strong></td><td>4.971e-03</td><td><strong>1.9221</strong></td><td>9.864e-03</td><td><strong>1.8912</strong></td></tr></tbody></table>

<table><caption>Table: Traveling vortex test: EOC for the ImEx 2 scheme, CFL = 0.5, T = 1/6</caption><thead><tr><th>N</th><th>L<sup>1</sup>-error in h</th><th>EOC</th><th>L<sup>1</sup>-error in h<sub>u</sub></th><th>EOC</th><th>L<sup>1</sup>-error in h<sub>v</sub></th><th>EOC</th></tr></thead><tbody><tr><td>40</td><td>1.282e-08</td><td></td><td>2.564e-01</td><td></td><td>4.182e-01</td><td></td></tr><tr><td>80</td><td>2.850e-09</td><td><strong>2.1697</strong></td><td>7.029e-02</td><td><strong>1.8670</strong></td><td>1.283e-01</td><td><strong>1.7047</strong></td></tr><tr><td>160</td><td>7.656e-10</td><td><strong>1.8962</strong></td><td>1.817e-02</td><td><strong>1.9514</strong></td><td>3.527e-02</td><td><strong>1.8629</strong></td></tr><tr><td>320</td><td>1.701e-10</td><td><strong>2.1698</strong></td><td>4.784e-03</td><td><strong>1.9255</strong></td><td>9.490e-03</td><td><strong>1.8941</strong></td></tr></tbody></table>
<table><caption>Table: CPU times in seconds consumed by the ImEx 2 and HLL-MUSCL-RK2 explicit schemes on different grids for different values of Fr. The final time is $T = 1/6$.</caption><thead><tr><th rowspan="2">Grid</th><th colspan="2">Fr = 1</th><th colspan="2">Fr = 0.1</th><th colspan="2">Fr = 0.01</th><th colspan="2">Fr = 0.001</th></tr><tr><th>ImEx2</th><th>Explicit</th><th>ImEx2</th><th>Explicit</th><th>ImEx2</th><th>Explicit</th><th>ImEx2</th><th>Explicit</th></tr></thead><tbody><tr><td>40 × 40</td><td>8.65</td><td>0.62</td><td>5.02</td><td>2.29</td><td>4.83</td><td>13.89</td><td>5.14</td><td>138.44</td></tr><tr><td>80 × 80</td><td>38.36</td><td>7.54</td><td>39.67</td><td>31.79</td><td>41.59</td><td>278.30</td><td>41.89</td><td>2,052.39</td></tr><tr><td>200 × 200</td><td>788.57</td><td>41.30</td><td>714.20</td><td>220.91</td><td>715.94</td><td>2,133.62</td><td>746.72</td><td>16,718.47</td></tr></tbody></table>
**Conclusion:**

* Extended the CPR scheme to second order in time while limiting the number of linear systems to solve.
**Conclusion:**

* Extended the CPR scheme to second order in time while limiting the number of linear systems to solve.

* We developed a stable and well-balanced scheme with weak numerical dissipation.
**Conclusion:**

* Extended the CPR scheme to second order in time while limiting the number of linear systems to solve.

* We developed a stable and well-balanced scheme with weak numerical dissipation.

* We demonstrated the theoretical second-order rate of convergence of ImEx 2 schemes in one and two dimensional settings and in low Froude number regimes.
**Conclusion:**

* Extended the CPR scheme to second order in time while limiting the number of linear systems to solve.

* We developed a stable and well-balanced scheme with weak numerical dissipation.

* We demonstrated the theoretical second-order rate of convergence of ImEx 2 schemes in one and two dimensional settings and in low Froude number regimes.

**Ongoing work:**

* Investigate the stability of the scheme by a Von Neumann analysis, for the entropy stability.
**Conclusion:**

* Extended the CPR scheme to second order in time while limiting the number of linear systems to solve.

* We developed a stable and well-balanced scheme with weak numerical dissipation.

* We demonstrated the theoretical second-order rate of convergence of ImEx 2 schemes in one and two dimensional settings and in low Froude number regimes.

**Ongoing work:**

* Investigate the stability of the scheme by a Von Neumann analysis, for the entropy stability.

**Perspectives:**

* Adaptation of the scheme to implement in AEROSOL/UHAINA. with V. Perrier and L. Cirrottola
**Conclusion:**

* Extended the CPR scheme to second order in time while limiting the number of linear systems to solve.

* We developed a stable and well-balanced scheme with weak numerical dissipation.

* We demonstrated the theoretical second-order rate of convergence of ImEx 2 schemes in one and two dimensional settings and in low Froude number regimes.

**Ongoing work:**

* Investigate the stability of the scheme by a Von Neumann analysis, for the entropy stability.

**Perspectives:**

* Adaptation of the scheme to implement in AEROSOL/UHAINA. *with V. Perrier and L. Cirrottola*

* Adaptation of the scheme to the complete multi-physic model accounting for atmospheric and astronomical forces... *with V. Pilorget*
**Conclusion:**

* Extended the CPR scheme to second order in time while limiting the number of linear systems to solve.

* We developed a stable and well-balanced scheme with weak numerical dissipation.

* We demonstrated the theoretical second-order rate of convergence of ImEx 2 schemes in one and two dimensional settings and in low Froude number regimes.

**Ongoing work:**

* Investigate the stability of the scheme by a Von Neumann analysis, for the entropy stability.

**Perspectives:**

* Adaptation of the scheme to implement in AEROSOL/UHAINA. *with V. Perrier and L. Cirrottola*

* Adaptation of the scheme to the complete multi-physic model accounting for atmospheric and astronomical forces... *with V. Pilorget*

**Thank you for your attention!**