# Rotational shallow water equations with viscous damping and boundary control: structure-preserving spatial discretization

Flávio Luiz Cardoso-Ribeiro, Ghislain Haine, Laurent Lefèvre, Denis Matignon

► To cite this version:

Flávio Luiz Cardoso-Ribeiro, Ghislain Haine, Laurent Lefèvre, Denis Matignon. Rotational shallow water equations with viscous damping and boundary control: structure-preserving spatial discretization. Mathematics of Control, Signals, and Systems, 2024, 10.1007/s00498-024-00404-6. hal-04835390

HAL Id: hal-04835390

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

Submitted on 13 Dec 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Rotational shallow water equations with
viscous damping and boundary control:
structure-preserving spatial discretization

Flávio Luiz Cardoso-Ribeiro¹*†, Ghislain Haine²†, Laurent Lefèvre³† and Denis Matignon²†

¹*Divisão de Engenharia Aeronáutica e Aeroespacial, Instituto Tecnológico de Aeronáutica, São José dos Campos, Brazil.

²Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, Toulouse, France.

³Laboratoire de Conception et d'Intégration des Systèmes (LCIS), Univ. Grenoble Alpes, Grenoble INP, Valence, France.

*Corresponding author(s). E-mail(s): flaviocr@ita.br;

Contributing authors: ghislain.haine@isae-supaero.fr;

laurent.lefevre@lcis.grenoble-inp.fr;

denis.matignon@isae-supaero.fr;

†These authors contributed equally to this work.

**Abstract**

This paper is dedicated to structure-preserving spatial discretiza-
tion of shallow water dynamics. First, a port-Hamiltonian for-
mulation is provided for the two-dimensional rotational shallow
water equations with viscous damping. Both tangential and nor-
mal boundary port variables are introduced. Then the corre-
sponding weak form is derived and a partitioned finite ele-
ment method is applied to obtain a finite-dimensional continuous-
time port-Hamiltonian approximation. Four simulation scenarios are
investigated to illustrate the approach and show its effectiveness.

**Keywords:** Shallow Water Equations (SWE), Port-Hamiltonian Systems (pHs), Viscous Damping, Partitioned Finite Element Method (PFEM)

MSC Classification: 76D55, 35Q35, 76M10
2 *Shallow water equations with viscous damping*

# 1 Introduction

The shallow water equations (SWE) are a set of partial differential equations that describe the motion of water in shallow areas like lakes, rivers, and coastlines. These equations are useful for understanding the dynamics of water flow, as well as predicting flooding and other related phenomena. They are also used to simulate the effects of wind and waves on the ocean surface. The study of SWE can be beneficial for a variety of purposes, including stability and control of moving tanks [1], coastal structure design, tsunami simulation [2, 3], flood forecasting and mitigation [4], and water resource management.

The two-dimensional rotational SWE (2D-SWE) are a useful model to study magnetohydrodynamic flows [5], geophysical fluids such as on large scale oceanic flows [6], or even the flow induced by the breaking of a large dam on a dry bed [7]. We are motivated by control applications including the boundary control of sloshing water in a rotating tank or non-radial boundary control of large circular water basin [8] such as in the FloWave Ocean Energy Research Facility [9], in Nantes, France.

A common approach to include dissipation in the shallow water model is to consider a fluid friction with the water bed, represented through some empirical formulae, which relates the local tension to the fluid velocity, and results in an additional nonlinear algebraic term in the SWE (see for instance [10, 11] or [12], for a list of such formulas for various physical configurations). However, these empirical formulas will not lead to the possibility of controlling the boundary tangential component of the water flow and will be useless in the applications motivating this work. Fortunately, the derivation of shallow water models from the Navier-Stokes equations for free surface flows with viscous damping has been proposed either in [13], for the one-dimensional case, or in [14], for the general three-dimensional rotational case with irregular bed topography. From the asymptotic analysis of the Navier-Stokes free-surface problem, a 2D-SWE is then obtained in which a second-order Laplacian-like differential operator accounts for the viscous dissipation (see also Section 2.2, hereafter). This latter model, which encompasses viscous damping and rotational fluid dynamics, incorporates all the physical phenomena which are needed to address the control problems that motivate this work. Note that [15] addresses the global existence for the Cauchy problem for this equation in the whole plane, for small initial data.

For analysis or control design, it may be crucial to obtain information on the underlying geometric structure for the considered dynamical systems, together with an account for the stored and dissipated energy. In our case we wish to make use of the port-Hamiltonian (pH) formulation of the dynamics [16, 17], combined with passivity-based control designs (namely impedance matching and control by interconnection) to achieve boundary control of the shallow water flows through boundary actuation. This pH approach has been used for the modelling, simulation and control of distributed parameter systems for more than twenty years [18]. The SWE were previously studied in this pH framework (see for instance [19, 20]). In [21], the pH 1D model of an
open channel irrigation system is spatially discretized and energy-based control designs are used to regulate water levels and flows throughout the reaches via sluice gates actuation. In [22, 23] the 1D and 2D sloshing problems are considered: in both these references, the pH formulation is being used for the modelling, simulation and control of the fluid-structure coupled system. In [8] surface waves in a circular water tank are modelled using the pH approach and boundary passive feedback control is designed. The pH model and energy-based decomposition have also been proposed for the general Navier-Stokes equations [24].

A classical approach to embed viscous friction forces in the SWE pH model is the approximation of the viscous term that arises in the Navier-Stokes equations with a Laplacian [16]. In a previous work [25], we used this approach where the viscous terms are given as a function of the Laplace diffusion operator only. First results of simulations in 1D were promising. However, as reported in [26], the use of the Laplacian formulation in these equations is misleading and may lead to non-physical solutions (which violates the fundamental objectivity principle of continuum mechanics [27]) and oversimplified boundary conditions. Therefore, in the present paper, we start from the asymptotic approach and the resulting SWE with viscous friction, as presented in [13] (1D case) and [14] (2D case), to derive their corresponding pH formulation.

This paper is dedicated to the pH formulation of the 2D SWE with viscous damping and its structure-preserving spatial discretization. Preserving the underlying Dirac interconnection structure in the pH model results (among others) in energy conservation and associated dynamical properties such as stability and controllability. Mixed finite element methods were introduced a long time ago to perform such structure-preserving spatial discretization for the Maxwell field equations [28]. In [29] this approach is applied to the weak formulation of pH models to obtain a quite general class of structure-preserving spatial discretization methods. In [30] we have shown that a partitioned version of these mixed finite element methods directly leads to a finite-dimensional Dirac interconnection structure and no further projection is required to obtain finite-dimensional pH equations with a flexible choice for causality (i.e. for the boundary port variables). This partitioned finite element method (PFEM) method has been applied to the discretization of various 2D and 3D pH models with non-autonomous boundary conditions [31]. Accurate convergence results in the sense of numerical analysis have been obtained [32] which suggest a heuristic for the optimal choice of finite element conforming spaces.

In this paper, the PFEM is applied to the rotational SWE with viscous damping. It is structured as follows. We propose first, in Section 2, a pH model for the rotational SWE with the appropriate viscous damping operator (see Section 2.2.1), derived from [13, 14]. We derive, in Section 2.2.2, an extended Stokes-Dirac interconnection structure, associated with a weighted inner product (designed specifically to simplify the SWE dynamics, see Section 2.1.2), and the appropriate factorization for the dissipation operator. In particular, it is shown that the correct normal and tangential boundary port-variables
4 Shallow water equations with viscous damping

both emerge from the power balance associated with this Stokes-Dirac inter-
connection structure. Finally, we compute the weak formulation (Section
3.1) and apply the PFEM spatial discretization to the obtained pH model
(Section 3.2). The finite-dimensional power balance is computed in Subsec-
tion 3.3. In Section 4, we validate the approach on four different simulation
scenarios:

1. a 2D simulation for the SWE in a rectangular closed tank, without dissipation and with a zero normal boundary control;

2. the same scenario with viscous dissipation, which produces tangential components for the velocities and a boundary limit layer;

3. a scenario where the water tank is emptied by forcing a normal component of the velocity (i.e., normal velocity boundary control);

4. a 2D simulation for a 2D rotating circular tank with uniform initial water level profile (i.e., tangential velocity boundary control).

To the best of our knowledge, the main contributions in this work are:

* the pH formulation of a physically coherent model for the 2D rotational SWE with viscous damping;

* the derivation of an extended 2D Stokes-Dirac interconnection structure with dissipation ports such that normal and tangential boundary-port variables may be identified in the power balance equation;

* the structure-preserving spatial discretization (with a PFEM-Galerkin approach) of the previous model to produce a finite-dimensional pH model with discrete power balance;

* the numerical results obtained for the 2D rotational SWE with viscous damping and normal/tangential boundary controls.

2 Dissipative SWE as a pH system

This section reviews the Shallow Water Equations written in the pH frame-
work. Firstly, Section 2.1 presents the frictionless model, without any source of
dissipation. This form leads to a power-preserving pH system, where the only
energy exchange occurs through the boundary ports. Secondly, the viscous
dissipative model is presented in Section 2.2.

## 2.1 Frictionless model

We consider a bounded domain in Ω ⊂ R², with a C² boundary ∂Ω, see e.g. [33, Chap. III]. For any finite time horizon T > 0, we assume strong solutions belonging to the functional space C¹([0, T]; C¹(Ω)). Moreover, at this stage, we formulate the hypothesis that, provided the initial water height h₀ is sufficiently high, the height h(t, x) will remain bounded from below by some h_min > 0, which might depend on T.
### 2.1.1 Classical setting

The frictionless irrotational SWE is usually written as a system of two conservation laws as:

$$ \left( \begin{matrix} \partial_t h \\ \partial_t \boldsymbol{p} \end{matrix} \right) = \left[ \begin{matrix} 0 & -\operatorname{div} \\ -\operatorname{grad} & 0 \end{matrix} \right] \left( \begin{matrix} e_h \\ e_p^0 \end{matrix} \right), \qquad (1) $$

where $h$ is the height of the fluid, $v$ is the velocity, $\rho$ is the fluid density (supposed constant), $\boldsymbol{p} := \rho v$ is the linear momentum, $e_h = \frac{1}{2}\rho \|v\|^2 + \rho gh$ is the total pressure and $e_p^0 = hv$ is the volumetric flow of the fluid. Thus, the first line of the matrix equation represents the conservation of the mass (or volume, since the fluid is assumed to be incompressible) and the second represents the conservation of linear momentum.

Furthermore, one can define the system Hamiltonian (or total energy) as a functional of $h$ and $\boldsymbol{p}$, which are thus called energy variables:

$$ \mathcal{H}(h, \boldsymbol{p}) := \int_{\Omega} \frac{1}{2} \rho h \|\boldsymbol{v}\|^2 + \frac{1}{2} \rho g h^2 \, d\Omega = \int_{\Omega} \frac{1}{2} \rho h \|\boldsymbol{p}\|^2 + \frac{1}{2} \rho g h^2 \, d\Omega. \quad (2) $$

The co-energy variables can be computed from the variational derivative of the Hamiltonian such that:

$$ \begin{aligned} e_h &:= \delta_h \mathcal{H} = \frac{1}{2} \rho \|\boldsymbol{v}\|^2 + \rho gh, \\ e_p^0 &:= \delta_p \mathcal{H} = h \boldsymbol{v}. \end{aligned} \qquad (3) $$

The time-derivative of the Hamiltonian can then be obtained as follows:

$$ \frac{d}{dt} \mathcal{H} = \int_{\Omega} \left( \partial_t h e_h + \partial_t \boldsymbol{p} \cdot \boldsymbol{e}_p^0 \right) d\Omega. \qquad (4) $$

From (1), we get that the time-derivative of the Hamiltonian depends only on the boundary variables:

$$ \frac{d}{dt} \mathcal{H} = - \int_{\partial \Omega} e_h \boldsymbol{e}_p^0 \cdot \boldsymbol{n} \, ds, \qquad (5) $$

which enables to define collocated control and observation distributed ports along the boundary $\partial \Omega$. For example:

$$ \begin{aligned} u_{\partial} &= -e_p^0 \cdot \boldsymbol{n}, & \text{volumetric flow,} \\ y_{\partial} &= e_h, & \text{total pressure,} \end{aligned} \qquad (6) $$

and the power-balance is given by a product between input and output boundary ports:

$$ \frac{d}{dt} \mathcal{H} = \int_{\partial \Omega} u_{\partial} y_{\partial} \, ds. \qquad (7) $$
6 *Shallow water equations with viscous damping*

For the rotational case, a slightly modified version of (1) can be defined,
that takes into account the vorticity of the fluid, see e.g., [20, 23]:

$$
\begin{pmatrix} \partial_t h \\ \partial_t \boldsymbol{p} \end{pmatrix} = \begin{bmatrix} 0 & -\operatorname{div} \\ -\operatorname{grad} & \frac{1}{h} G(\omega) \end{bmatrix} \begin{pmatrix} e_h \\ \boldsymbol{e}_p^0 \end{pmatrix}, \qquad (8)
$$

with (scalar) vorticity $\omega(t, \boldsymbol{x})$ and gyroscopic term $G(\omega)$ defined by:

$$
\omega := \mathrm{curl}_{2D} \boldsymbol{v} = \partial_x v_2 - \partial_y v_1, \quad \text{and} \quad G(\omega) := \rho \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \omega. \qquad (9)
$$

Since the matrix $G$ is skew-symmetric, it will play no role in the power balance
(and it computes exactly as (7)).

2.1.2 A new choice of scalar product

As mentioned in [34], following e.g., [35], it can be interesting to adapt the
chosen scalar product to the physical problem: indeed, since the shallow water
model performs an average in one dimension, i.e., on the height of the water
column, it is natural to introduce for the velocities the scalar product in L²h(Ω):

$$
(\mathbf{v}_1, \mathbf{v}_2)_h := \int_{\Omega} \mathbf{v}_1 \cdot \mathbf{v}_2 \, h \, d\Omega . \qquad (10)
$$

Here the hypothesis $\forall x \in \Omega, \forall t \in [0,T], h(t, x) \ge h_{\min} > 0$ is of utmost importance, otherwise (10) would not define a scalar product. In this case, the computation of the co-energy variables must be adapted, since the definition of the variational derivative, see e.g. [36, 37], requires the use of Riesz representation theorem, which is dependent on the choice of the scalar product. Thus, scalar fields belong to $L^2(\Omega)$, while vector fields belong to $L_h^2(\Omega)$:

$$
e_h := \delta_h \mathcal{H} = \frac{1}{2} \rho \|\mathbf{v}\|^2 + \rho g h, \qquad (11)
$$

$$
e_p := \delta_p^h \mathcal{H} = \mathbf{v}.
$$

Indeed the characterization $(\delta_p^h \mathcal{H}, \varepsilon \mathbf{q})_h = \mathcal{H}(h, \mathbf{p} + \varepsilon \mathbf{q}) - \mathcal{H}(h, \mathbf{p}) + O(\varepsilon^2)$,
$\forall \mathbf{q} \in L_h^2(\Omega)$, leads to the new definition $e_p := \delta_p^h \mathcal{H} = \mathbf{v}$, which will be used
in the sequel, and provides some advantages at the numerical level.

In this case, the differential operators have to be adapted, and we can find:

$$
\begin{pmatrix} \partial_t h \\ \partial_t \mathbf{p} \end{pmatrix} = \begin{bmatrix} 0 & -\operatorname{div}(h \cdot) \\ -\operatorname{grad} & G(\omega) \end{bmatrix} \begin{pmatrix} e_h \\ \mathbf{e}_p^0 \end{pmatrix}. \qquad (12)
$$

In what follows, the *formal* adjoint is denoted by *, while the formal adjoint
with respect to the new scalar product is denoted by *$^h$*. Both should not be
confused with the usual adjoint. By formal adjoint, it is meant in the sense of distributions, i.e., for all compactly-supported functions in $C^\infty(\Omega)$.

An operator that depends on the energy variables is said to be *modulated*, in accordance with the terminology of *modulated Dirac structure*, see e.g., [16, Section 2.2.2] or [17].

Let us denote by $H^1(\Omega) \subset L^2(\Omega)$ the Sobolev space of scalar fields with first order weak partial derivative in $L^2(\Omega)$, and:

$$H_h^{\text{div}}(\Omega) := \{\psi \in L_h^2(\Omega), \text{div}(h \psi) \in L^2(\Omega)\} \subset L_h^2(\Omega).$$

**Proposition 1** Assume that there exist regular energy variables, $h \in C^1([0, T]; C^1(\bar{\Omega}))$, and $p \in C^1([0, T]; (C^1(\bar{\Omega}))^2)$, solutions of the SWE (12)-(11). The unbounded matrix-valued differential operator $\mathcal{J}(h, p) : \mathcal{D}(\mathcal{J}(h, p)) := H^1(\Omega) \times H_h^{\text{div}}(\Omega) \to L^2(\Omega) \times L_h^2(\Omega)$:

$$\mathcal{J}(h, p) := \begin{bmatrix} 0 & -\text{div}(h \cdot) \\ -\text{grad} & G(\omega) \end{bmatrix}, \qquad (13)$$

is modulated by the energy variables $h$ and $p$, and is formally skew-symmetric.

**Proof**

To prove the formal skew-symmetry of $\mathcal{J}(h, p)$, let us first recall that the bounded term $J_{22} := G(\omega)$ is a real-valued $2 \times 2$ skew-symmetric matrix. Note also that $\rho\omega = \text{curl}_{2D}(p) \in C^1([0, T]; C^0(\bar{\Omega}))$; thus, $\forall \psi \in L_h^2(\Omega), G(\omega)\psi \in L_h^2(\Omega)$.

Now, for the unbounded terms $J_{12}$ and $J_{21}$, taking care of the new scalar product (10) for vector-valued variables, one has $\forall \phi \in H^1(\Omega), \forall \psi \in H_h^{\text{div}}(\Omega)$:

$$(-\text{grad}\,\phi, \psi)_h = (-\text{grad}\,\phi, h\psi), \\ = +(\phi, \text{div}(h\psi)) - \int_{\partial\Omega} \phi\psi \cdot n \, h \, ds. \qquad (14)$$

Thus, for fields $\phi$, $\psi$ vanishing at the boundary, the formal adjoint of $-\text{grad}$ now is $\text{div}(h \cdot)$. And in the general case, the boundary term reads $-\int_{\partial\Omega} \phi \psi \cdot n \, h \, ds$. Finally, this also shows that $J_{12} \in \mathcal{L}(H_h^{\text{div}}(\Omega), L^2(\Omega))$, and that $J_{21} \in \mathcal{L}(H^1(\Omega), L_h^2(\Omega))$. Altogether, $\mathcal{J}(h, p) \in \mathcal{L}(H^1(\Omega) \times H_h^{\text{div}}(\Omega), L^2(\Omega) \times L_h^2(\Omega))$. $\square$

The possible choice (6) of collocated input – output ports has to be adapted accordingly:

$$\begin{align}
u_{\partial} &= -e_{p} \cdot n = -v \cdot n, \quad \text{normal velocity,} \tag{15} \\
y_{\partial} &= e_{h}, \quad \text{total pressure,}
\end{align}$$

and the power-balance (7) now reads:

$$\frac{d}{dt} \mathcal{H} = \int_{\partial \Omega} u_{\partial} y_{\partial} h \, ds := (u_{\partial}, y_{\partial})_{h, \partial \Omega}. \qquad (16)$$
8 *Shallow water equations with viscous damping*

## 2.2 Linear viscous (differential) dissipation

Now let us investigate the linear dissipation model induced by the Navier-Stokes equation, averaged on a slice of fluid: the full model will be recalled first in Section 2.2.1, and recast as a pH system (an explicit distributed pH system, giving rise to a pH-ODE once spatially discretized) under the form $J - R$, helping prove dissipativity in an easy way. Then in Section 2.2.2, enlightening a factorization of $R$ as $GSG^*$ gives rise to the definition of physically meaningful dissipation ports and a pH system with an extended structure matrix $J_e$ (an implicit distributed pH system, leading to a pH-DAE once spatially discretized, see [38] for a review), allowing a more straightforward computation of the boundary terms in the energy balance: here the tangential component of the velocity does play a role as an additional boundary control port.

### 2.2.1 Modelling as a dissipative dynamical system

The viscous effects on the SWE can be obtained from the Navier-Stokes equations, which involves a diffusion term with water viscosity $\mu > 0$ in the momentum equation:

$$ \partial_t p = - \mathrm{grad}(e_h) + G(\omega)e_p - h^{-1}\tau. \quad (17) $$

However, the characterization of the viscous stress $\tau$ proves intricate. Already in the 1D case, following e.g. [39], it was suggested in [25] to use $\tau := -\mu \partial_{xx} v$; however, the careful derivation of the viscous model from the Navier-Stokes equations averaged on a slice proved in [13], gives instead:

$$ \tau := -4\mu \partial_x (h \partial_x v). $$

In the 2D case, things become even more intricate: it was first suggested in [25] to use $\tau := -\mu \Delta v$, involving the vector Laplacian $\Delta$, which appears classically in the Navier-Stokes equations; however the careful derivation of the model performed in [14] leads to a much more complex definition of the viscous stress, namely:

$$ \tau := - \mathrm{grad}(2\mu h \mathrm{div}(v)) - \mathrm{Div}(2\mu h \mathrm{Grad}(v)), \quad (18) $$

where $\mathrm{Grad}(v)$ is the symmetric tensor $\frac{1}{2}(\nabla v + \nabla v^\top)$, and $d = \mathrm{Div}(\Sigma)$ is the vector of the divergences of the columns of the tensor; i.e. $d_i = \mathrm{div}(\Sigma(\cdot, i))$.

**Proposition 2** *The unbounded matrix-valued differential operator: $\mathcal{R}(h) : \mathcal{D}(\mathcal{R}(h)) := L^2(\Omega) \times \mathcal{D}(R_{22}(h)) \to L^2(\Omega) \times L_h^2(\Omega)$:*

$$ \mathcal{R}(h) := \begin{bmatrix} 0 & 0 \\ 0 & -\frac{1}{h} \mathrm{grad}(2\mu h \mathrm{div}(\cdot)) - \frac{1}{h} \mathrm{Div}(2\mu h \mathrm{Grad}(\cdot)) \end{bmatrix}, \quad (19) $$

$\mathcal{D}(R_{22}(h)) := \{w \in L_h^2(\Omega) \mid \mathrm{grad}(2\mu h \mathrm{div}(w)) + \mathrm{Div}(2\mu h \mathrm{Grad}(w)) \in L^2(\Omega)\}$, is modulated by the energy variable $h$, and formally symmetric and non-negative.
**Proof** Since $\mathcal{R}(h)$ acts on the second variable only through $R_{22}(h)$, let us compute for $C^\infty$ fields $v$, $w$ vanishing at the boundary $\partial\Omega$:

$$ 
\begin{aligned}
(R_{22}(h)v, w)_h &= \left(-\frac{1}{h}\mathbf{grad}(2\mu h \operatorname{div}(v)) - \frac{1}{h}\operatorname{Div}(2\mu h \operatorname{Grad}(v)), w\right)_h \\
&= \left(-\mathbf{grad}(2\mu h \operatorname{div}(v)) - \operatorname{Div}(2\mu h \operatorname{Grad}(v)), w\right) \\
&= (2\mu h \operatorname{div}(v), \operatorname{div}(w)) + (2\mu h \operatorname{Grad}(v), \operatorname{Grad}(w)) \\
&= 2\mu (\operatorname{div}(v), \operatorname{div}(w))_h + 2\mu (\operatorname{Grad}(v), \operatorname{Grad}(w))_h,
\end{aligned}
$$

from this latter line, it is obvious that $R_{22}(h)$ is formally symmetric; moreover, for the special case $w = v$ we obtain:

$$ 
\begin{aligned}
(R_{22}(h)v, v)_h &= 2\mu (\operatorname{div}(v), \operatorname{div}(v))_h + 2\mu (\operatorname{Grad}(v), \operatorname{Grad}(v))_h \\
&= \int_{\Omega} 2\mu (\operatorname{div}(v))^2 + 2\mu \operatorname{Grad}(v) : \operatorname{Grad}(v) h \, d\Omega \geq 0,
\end{aligned}
$$

which shows that the operator is formally non-negative. □

As a consequence, for fields vanishing at the boundary $\partial\Omega$, the dynamical system proves dissipative:

$$ \frac{d}{dt} \mathcal{H} = - \int_{\Omega} 2\mu (\operatorname{div}(v))^2 + 2\mu \operatorname{Grad}(v) : \operatorname{Grad}(v) h \, d\Omega \leq 0; \quad (20) $$

this kind of mechanical energy balance is in perfect accordance with classical results from the literature, see e.g. [33, 40].

Since we are most interested in boundary control, we must go further and make the boundary terms explicit, either by one integration by parts from this first compound formulation, or by extending the formulation with dissipative ports, which will make the computation a bit easier.

### 2.2.2 A port-Hamiltonian system with dissipation ports

The 2D rotational SWE with viscous damping boils down to:

$$ 
\left\{ 
\begin{array}{l}
\partial_t h = -\operatorname{div}(h e_p), \\
h \partial_t \boldsymbol{p} = -h \mathbf{grad}(e_h) + h \operatorname{curl}_{2D}(\boldsymbol{p}) \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} e_p \\
\qquad + 2\mu \operatorname{Div}(h \operatorname{Grad}(e_p)) + 2\mu \mathbf{grad}(h \operatorname{div}(e_p)),
\end{array}
\right. 
\quad (21) $$

together with the constitutive relations (11) given by:

$$ 
\left\{ 
\begin{array}{l}
e_h = \rho gh + \frac{\|\boldsymbol{p}\|^2}{2\rho}, \\
e_p = \frac{\boldsymbol{p}}{\rho}.
\end{array}
\right. 
\quad (22) $$
10 Shallow water equations with viscous damping

The general factorization of dissipation operators proposed in [41] under the form $R = GSG^*$ can be illustrated here, and provides a useful step towards the pH modelling of dissipative systems. Indeed, it is known that $\mathrm{div}^* = - \mathbf{grad}$, and it has been proved in [42] that $\mathrm{Div}^* = - \mathrm{Grad}$ for symmetric tensors. Note here that for $\mathbf{p}$-type variables, the $h$-dependent scalar product must be used, and the formally adjoint operators recomputed adequately.

We are now in a position to define dissipation flows, as follows: (they both apply to $e_p = v$),

* $F_d := \mathrm{Grad}(v) \in L^2_{\mathrm{sym}}(\Omega)$ a symmetric tensor, the so-called strain rate tensor,

* $f_d := \mathrm{div}(v) \in L^2(\Omega)$, a scalar field, the velocity divergence,

which are both physically meaningful. In the above, $L^2_{\mathrm{sym}}(\Omega)$ is the space of second order symmetric tensor fields with components in $L^2(\Omega)$. The $h$-dependent version of these spaces will be useful in the sequel; $L^2_{\mathrm{sym},h}(\Omega)$ with scalar product $(E, F)_h := \int_{\Omega} E : F h \, d\Omega$, and $L^2_h(\Omega)$ with scalar product $(d_1, d_2)_h := \int_{\Omega} d_1 d_2 h \, d\Omega$.

Let us denote by $H^1(\Omega) \subset L^2(\Omega)$, $H^1(\Omega) \subset L^2(\Omega)$ and $H^1_{\mathrm{sym}}(\Omega) \subset L^2_{\mathrm{sym}}(\Omega)$ the Sobolev spaces of fields with first order partial derivative of each component in $L^2(\Omega)$, and their $h$-dependent version $H^1_{\mathrm{sym},h}(\Omega) \subset L^2_{\mathrm{sym},h}(\Omega)$ and $H^1_h(\Omega) \subset L^2_h(\Omega)$.

**Proposition 3** The 2D-SWE with viscous damping (21) can be recast as the following constrained dynamical system:

$$ \begin{pmatrix} \partial_t h \\ \partial_t \mathbf{p} \\ F_d \\ f_d \end{pmatrix} = \mathcal{J}_e(h, \mathbf{p}) \begin{pmatrix} e_h \\ e_p \\ E_d \\ e_d \end{pmatrix}, \qquad (23) $$

together with the constitutive relations (22) and the additional closure relations:

$$ e_d = 2 \mu f_d, \quad \text{and} \quad E_d = 2 \mu F_d. \qquad (24) $$

The extended operator $\mathcal{J}_e(h, \mathbf{p}) : \mathcal{D}(\mathcal{J}_e(h, \mathbf{p})) \to \mathcal{X}$, defined by:

$$ \begin{aligned} \mathcal{D}(\mathcal{J}_e(h, \mathbf{p})) &:= H^1(\Omega) \times H^1(\Omega) \times H^1_{\mathrm{sym},h}(\Omega) \times H^1_h(\Omega), \\ \mathcal{X} &:= L^2(\Omega) \times L^2_h(\Omega) \times L^2_{\mathrm{sym},h}(\Omega) \times L^2_h(\Omega), \end{aligned} $$

and:

$$ \mathcal{J}_e(h, \mathbf{p}) := \begin{bmatrix} 0 & -\mathrm{div}(h \cdot) & 0 & 0 \\ -\mathbf{grad} & G(\omega) & \frac{1}{h} \mathrm{Div}(h \cdot) & \frac{1}{h} \mathbf{grad}(h \cdot) \\ 0 & \mathrm{Grad} & 0 & 0 \\ 0 & \mathrm{div} & 0 & 0 \end{bmatrix}, \qquad (25) $$

is formally skew-symmetric.

**Proof.** Firstly, let us verify that we can recast the dynamic equations as (23). Starting from (21), and substituting the dissipation flows $F_d = \mathrm{Grad}(e_p)$ and $f_d = \mathrm{div}(e_p)$ leads to:
$$ \begin{pmatrix} \partial_t h \\ \partial_t \mathbf{p} \end{pmatrix} = \begin{bmatrix} 0 & -\operatorname{div}(h \cdot) \\ -\operatorname{grad} & G(\omega) \end{bmatrix} \begin{pmatrix} e_h \\ \mathbf{e}_p \end{pmatrix} + \begin{pmatrix} 0 \\ h^{-1} \operatorname{grad}(2\mu h f_d) \end{pmatrix} + \begin{pmatrix} 0 \\ h^{-1} \operatorname{Div}(2\mu h F_d) \end{pmatrix}. $$

With the closure relations $E_d = 2\mu F_d$ and $e_d = 2\mu f_d$, this rearranges as:

$$ \begin{pmatrix} \partial_t h \\ \partial_t \mathbf{p} \end{pmatrix} = \begin{bmatrix} 0 & -\operatorname{div}(h \cdot) & 0 & 0 \\ -\operatorname{grad} & G(\omega) & \frac{1}{h} \operatorname{Div}(h \cdot) & \frac{1}{h} \operatorname{grad}(h \cdot) \end{bmatrix} \begin{pmatrix} e_h \\ \mathbf{e}_p \\ E_d \\ e_d \end{pmatrix}, $$

which yields (23) by definition of $F_d$ and $f_d$.

Secondly, to verify that $\mathcal{J}_e(h, \mathbf{p})$ is formally skew-symmetric, it proves necessary to compute its formal adjoint, where the $h$-dependent scalar product (10) is used for vector fields, but also for tensor fields on $L_{\text{sym},h}^2(\Omega)$ and scalar fields on $L_h^2(\Omega)$.

This proceeds in 3 steps: computation of the formal adjoint of $\frac{1}{h} \operatorname{grad}(h \cdot)$, of $\frac{1}{h} \operatorname{Div}(h \cdot)$, and factorization of the damping $\mathcal{R}(h)$ leading to formal skew-symmetry of $\mathcal{J}_e(h, \mathbf{p})$.

**Step 1** for any scalar field $\varphi \in H_h^1(\Omega)$ and vector field $\psi \in H^1(\Omega)$,

$$ \begin{aligned} \left(-\frac{1}{h} \operatorname{grad}(h \varphi), \psi\right)_h &= \left(-\operatorname{grad}(h \varphi), \psi\right) \\ &= (h \varphi, \operatorname{div}(\psi)) - \int_{\partial \Omega} h \varphi \psi \cdot \mathbf{n} \, ds \\ &= (\varphi, \operatorname{div}(\psi))_h - \int_{\partial \Omega} \varphi \psi \cdot \mathbf{n} \, h \, ds. \end{aligned} \quad (26) $$

Thus, the formal adjoint of $-\frac{1}{h} \operatorname{grad}(h \cdot)$ on $L_h^2(\Omega)$ is $\operatorname{div}$; and $-\int_{\partial \Omega} \varphi \psi \cdot \mathbf{n} \, h \, ds$ is the boundary term.

**Step 2** for any vector field $\mathbf{w} \in H^1(\Omega)$ and symmetric tensor field $E \in \mathbb{H}_{\text{sym},h}^1(\Omega)$,

$$ \begin{aligned} \left(-\frac{1}{h} \operatorname{Div}(h E), \mathbf{w}\right)_h &= \left(-\operatorname{Div}(h E), \mathbf{w}\right) \\ &= (h E, \operatorname{Grad}(\mathbf{w})) - \int_{\partial \Omega} h \operatorname{En} \cdot \mathbf{w} \, ds \\ &= (E, \operatorname{Grad}(\mathbf{w}))_h - \int_{\partial \Omega} \mathbf{w} \cdot \operatorname{En} \, h \, ds. \end{aligned} \quad (27) $$

Thus, the formal adjoint of $-\frac{1}{h} \operatorname{Div}(h \cdot)$ on $L_{\text{sym},h}^2(\Omega)$ is $\operatorname{Grad}$; and $-\int_{\partial \Omega} \mathbf{w} \cdot \operatorname{En} \, h \, ds$ is the boundary term.
12 Shallow water equations with viscous damping

**Step 3:** Denoting $\mathcal{G} := \begin{bmatrix} 0 & 0 \\ \frac{1}{h} \mathrm{Div}(h \cdot) & \frac{1}{h} \mathrm{grad}(h \cdot) \end{bmatrix}$, from steps 1 and 2, we can conclude that $\mathcal{G}^{*h} = \begin{bmatrix} 0 & -\mathrm{Grad} \\ 0 & -\mathrm{div} \end{bmatrix}$. Then, we can rewrite (25) using the classical decomposition $\mathcal{J}_e(h, \mathbf{p}) = \begin{bmatrix} \mathcal{J}(h, \mathbf{p}) & \mathcal{G} \\ -\mathcal{G}^{*h} & 0 \end{bmatrix}$, where $\mathcal{J}(h, \mathbf{p})$ is defined in (13). We conclude to the formal skew-symmetry of $\mathcal{J}_e(h, \mathbf{p})$ with the help of Proposition 1. $\square$

**Remark 1** Note that if we define $S := \mathrm{diag}(2\mu I, 2\mu)$ and apply the closure relation $\begin{pmatrix} E_d \\ e_d \end{pmatrix} = S \begin{pmatrix} F_d \\ f_d \end{pmatrix}$ in (23), we get:

$$ \left( \frac{\partial_t h}{\partial t \mathbf{p}} \right) = (\mathcal{J}(h, \mathbf{p}) - \mathcal{R}) \begin{pmatrix} e_h \\ e_p \end{pmatrix}, \qquad (28) $$

where $\mathcal{R} = \mathcal{G} \mathcal{S} \mathcal{G}^{*h}$, which is equivalent to (19), and (28) is equivalent to (21).

Thanks to this strongly structured formulation of the original system, we are now in a position to state the following:

**Theorem 4** The evolution of the Hamiltonian (2) along the trajectories of the dynamical system (23) with the constitutive relations (11) and the closure relations (24) is given by:

$$ \frac{d}{dt} \mathcal{H} = - \int_{\Omega} \left[ 2\mu (\mathrm{div}(\mathbf{v}))^2 + 2\mu \mathrm{Grad}(\mathbf{v}) : \mathrm{Grad}(\mathbf{v}) \right] h \, d\Omega \qquad (29) $$

$$ + \int_{\partial\Omega} \left[ (-e_h + 2\mu \mathrm{div}(\mathbf{v}))\mathbf{v} \cdot \mathbf{n} + 2\mu \mathbf{v} \cdot \mathrm{Grad}(\mathbf{v})\mathbf{n} \right] h \, ds. \qquad (30) $$

The boundary term (30) alone can be further decomposed into:

$$ \int_{\partial\Omega} \left[ (-e_h + 2\mu \mathrm{div}(\mathbf{v}) + 2\mu \mathbb{N}(s)) : \mathrm{Grad}(\mathbf{v})\mathbf{v} \cdot \mathbf{n} + (2\mu \mathbb{R}(s) : \mathrm{Grad}(\mathbf{v}))\mathbf{v} \cdot \mathbf{t} \right] h \, ds, \qquad (31) $$

where $\mathbf{t}$ is the tangential vector on the boundary (+90 degrees rotation of $\mathbf{n}$, the outward normal vector), and the tensors at the boundary are defined locally by $\mathbb{N}(s) = \mathbf{n}\mathbf{n}^{\top}$, and $\mathbb{R}(s) = \mathbf{t}\mathbf{n}^{\top}$.

**Proof.** The proof proceeds in 3 steps: 1. computation of the scalar product $(e, \mathcal{J}_e e)$ (where $e$ is the vector with all co-energy variables) and identification of the derivative of the Hamiltonian, together with the in-domain terms, 2. careful computation of the right-hand side involving boundary terms only, 3. equivalent expression of the boundary integral after projection into local coordinates on the normal and tangential components.
Step 1

$$
\left( \begin{pmatrix} e_h \\ \mathbf{e}_p \\ E_d \\ e_d \end{pmatrix} \right), \mathcal{J}_e \left( \begin{pmatrix} e_h \\ \mathbf{e}_p \\ E_d \\ e_d \end{pmatrix} \right) = \int_{\Omega} (e_h \partial_t h + \mathbf{e}_p \cdot \partial_t \mathbf{p} + F_d : E_d + f_d e_d) \, h \, d\Omega \\
= \frac{d}{dt} \mathcal{H} + \int_{\Omega} [2\mu (\mathrm{div}(\mathbf{v}))^2 + 2\mu \mathrm{Grad}(\mathbf{v}) : \mathrm{Grad}(\mathbf{v})] \, h \, d\Omega .
$$

**Step 2** Computing $(e, \mathcal{J}_e e)$ carefully, one can group terms in three pairs and make use of the previous results on the formal adjoints (14), (26) and (27); indeed the three following contributions must be added:

$$
(e_h, -\mathrm{div}(h \mathbf{e}_p)) + (\mathbf{e}_p, -\mathrm{grad}(e_h))_h = -\int_{\partial\Omega} e_h \mathbf{v} \cdot \mathbf{n} \, h \, ds,
$$

$$
(e_p, \frac{1}{h} \mathrm{Div}(h E_d))_h + (E_d, \mathrm{Grad}(e_p))_h = \int_{\partial \Omega} 2\mu \mathbf{v} \cdot \mathrm{Grad}(\mathbf{v}) \mathbf{n} \, h \, ds,
$$

$$
\text{and } (e_p, \frac{1}{h} \mathbf{grad}(h e_d))_h + (e_d, \text{div}(e_p))_h = \int_{\partial \Omega} 2\mu \text{div}(\mathbf{v}) \mathbf{v} \cdot \mathbf{n} \, h \, ds .
$$

**Step 3** Using the identity $\mathbf{v} \cdot \mathrm{Grad}(\mathbf{v})\mathbf{n} = \mathrm{Grad}(\mathbf{v}) : (\mathbf{v}\mathbf{n}^{\top})$, and plugging the decomposition of the velocity into a local basis $\mathbf{v} = (\mathbf{v} \cdot \mathbf{n})\mathbf{n} + (\mathbf{v} \cdot \mathbf{t})\mathbf{t}$, leads to (31).

The conclusion of Theorem 4 follows.

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

Remark 2 When comparing (20) with (29)–(30), we can see that the closed dynamical system (i.e., without boundary control) is dissipative, while the open dynamical system is passive (up to the definition of appropriate collocated inputs and outputs). The latter decomposition (31) is physically meaningful: both the normal and the tangential components of the velocity vector at the boundary play a role. However, this result differs from the original computation made in [25] from the vectorial Laplacian formulation, and also from the energy balance for a compressible isentropic fluid given in [43].

A possible choice for the boundary ports could be the boundary velocity
(with both normal and tangential components) for control:

$$
\mathbf{u}(t, \mathbf{s}) = \mathbf{e}_p(t, \mathbf{s}), \quad \forall t \ge 0, \mathbf{s} \in \partial\Omega. \tag{32}
$$

Such a choice of boundary control implies, according to (30), that the colocated
boundary observation is given by:

$$
\mathbf{y}(t, \mathbf{s}) = (-e_h(t, \mathbf{s}) + 2\mu \mathrm{div}(\mathbf{e}_p)) \mathbf{n} + 2\mu \mathrm{Grad}(\mathbf{e}_p) \cdot \mathbf{n}, \quad \forall t \ge 0, \mathbf{s} \in \partial\Omega. \quad (33)
$$
14 Shallow water equations with viscous damping

Then, the power balance related with the boundary control (30) becomes:

$$ \frac{d}{dt} \mathcal{H} \leq \int_{\partial \Omega} \mathbf{y} \cdot \mathbf{u} \, h \, ds. \qquad (34) $$

*Remark 3* In the power balance (34), one can further decompose the density $\mathbf{y} \cdot \mathbf{u}$ into a non viscous term $y^0 u^0$ when $\mu = 0$, and a purely viscous contribution $y^\mu \cdot \mathbf{u}$ when $\mu > 0$. This motivates splitting the output variable (33) as follows:

$$ \mathbf{y} = y^0 \mathbf{n} + \mathbf{y}^\mu, \quad \text{with } y^0 := -e_h, \quad \text{and } \mathbf{y}^\mu := 2\mu (\text{Grad}(e_p) \cdot \mathbf{n} + \text{div}(e_p)\mathbf{n}). \quad (35) $$

The notation $u^0 := \mathbf{u} \cdot \mathbf{n} = u_n$ has also been introduced, since this latter quantity is the only control variable which is available when $\mu = 0$.

*Remark 4* In this paper, we assume that the fluid height is always positive. This hypothesis is necessary considering the definition of the inner product of Section 2.1.2. Proving that the solutions fulfill $h > 0$ is beyond the scope of this paper. For this reason, special care must be taken to ensure that the simulations meet this condition. The interested reader is referred to [15, 44, 45] and references therein for several results on the global existence of solutions of the viscous SWE. In particular, [15] presents a proof that the height is always positive, at least when assuming small initial perturbations around the equilibrium and small forcing. In this section, the strong regularity $C^1([0, T]; C^1(\bar{\Omega}))$ is assumed both on the height $h$, and on the linear momentum $\mathbf{p}$ (hence, on the velocity $\mathbf{v}$). In light of the results obtained in $\Omega = \mathbb{R}^2$, proven in [15], this strong assumption seems reasonable.

## 3 Structure-preserving discretization

We are now in a position to propose a structure-preserving discrete approximation of the dissipative SWE. At least two approaches may be used to guarantee a discrete power balance mimicking (29)–(30) on the infinite-dimensional level. The first one would rely on the flow/effort formulation of the problem, with dissipative ports together with linear constitutive relations to take the dissipation into account, see e.g. [46], the second one relies on the $J-R$ state formulation, as in Section 2.2.1, where $R$ is symmetric non-negative and represents the dissipation as one single block. While the former is very useful to enlighten the underlying Stokes-Dirac structure, it suffers from an increased number of unknowns to be computed. The latter seems more appropriate for numerical purpose, since it strongly reduces the number of degrees of freedom in the discrete system, the only price to pay is that a Lagrange multiplier will be needed to access the desired tangent control of the velocity.
## 3.1 Weak formulation

Let us consider the SWE (21)–(22) with the boundary control (32):

$$e_p(t, s) = u(t, s) = u_n(t, s)n(s) + u_t(t, s)t(s), \quad \forall t \ge 0, s \in \partial\Omega,$$

where $u_n$ and $u_t$ are the normal and tangential velocity respectively, and the colocated boundary observation (33):

$$y(t, s) = (-e_h(t, s) + 2\mu \operatorname{div}(e_p)) n + 2\mu \operatorname{Grad}(e_p) \cdot n, \quad \forall t \ge 0, s \in \partial\Omega.$$

Let us consider sufficiently regular test functions $\varphi$, $\Phi$ on $\Omega$. The weak form of (21) reads:

$$\left\{
\begin{array}{l}
(\partial_t h, \varphi)_{L^2} = -(\operatorname{div}(he_p), \varphi)_{L^2}, \\
(h \partial_t \mathbf{p}, \mathbf{\Phi})_{L^2} = -(h \operatorname{grad}(e_h), \mathbf{\Phi})_{L^2} + \left(h \operatorname{curl}_{2D}(\mathbf{p}) \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} e_p, \mathbf{\Phi} \right)_{L^2} \\
\qquad + 2\mu (\operatorname{Div}(h \operatorname{Grad}(e_p)), \mathbf{\Phi})_{L^2} \\
\qquad + 2\mu (\operatorname{grad}(h \operatorname{div}(e_p)), \mathbf{\Phi})_{L^2}.
\end{array}
\right.$$

Applying integration by parts, once on the first line, and on each dissipative terms, leads to:

$$\left\{
\begin{array}{l}
(\partial_t h, \varphi)_{L^2} = + (he_p, \operatorname{grad}(\varphi))_{L^2} - (he_p \cdot n, \varphi)_{\partial\Omega}, \\
(h \partial_t \mathbf{p}, \mathbf{\Phi})_{L^2} = - (h \operatorname{grad}(e_h), \mathbf{\Phi})_{L^2} + \left( h \operatorname{curl}_{2D}(\mathbf{p}) \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} e_p, \mathbf{\Phi} \right)_{L^2} \\
\qquad - 2\mu (h \operatorname{Grad}(e_p), \operatorname{Grad}(\mathbf{\Phi}))_{\mathbb{L}^2} \\
\qquad + 2\mu (h \operatorname{Grad}(e_p) \cdot n, \mathbf{\Phi})_{\partial\Omega} \\
\qquad - 2\mu (h \operatorname{div}(e_p), \operatorname{div}(\mathbf{\Phi}))_{L^2} \\
\qquad + 2\mu (\mathbf{\Phi} \cdot n, h \operatorname{div}(e_p))_{\partial\Omega},
\end{array}
\right.$$

where we recall that both solutions and test functions are assumed regular enough to ensure that the boundary duality brackets reduce to boundary $L^2$-scalar products.

Gathering the boundary terms of the second equation gives:

$$\left\{
\begin{array}{l}
(\partial_t h, \varphi)_{L^2} = + (he_p, \operatorname{grad}(\varphi))_{L^2} - (he_p \cdot n, \varphi)_{\partial\Omega}, \\
(h \partial_t \mathbf{p}, \mathbf{\Phi})_{L^2} = - (h \operatorname{grad}(e_h), \mathbf{\Phi})_{L^2} + \left( h \operatorname{curl}_{2D}(\mathbf{p}) \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} e_p, \mathbf{\Phi} \right)_{L^2} \\
\qquad - 2\mu (h \operatorname{Grad}(e_p), \operatorname{Grad}(\mathbf{\Phi}))_{\mathbb{L}^2} \\
\qquad - 2\mu (h \operatorname{div}(e_p), \operatorname{div}(\mathbf{\Phi}))_{L^2} \\
\qquad + 2\mu (h (\operatorname{Grad}(e_p) \cdot n + \operatorname{div}(e_p)n), \mathbf{\Phi})_{\partial\Omega}.
\end{array}
\right.$$

As a result, one can observe that parts of the boundary input/output, defined in (32) and (33), appear partitioned across the two lines of the system. In the first line, the boundary normal velocity $u^0 = e_p \cdot n$ appears, whereas in the second line, $y^\mu$ appears, as the viscous part of the output $y$, see (35).
16 Shallow water equations with viscous damping

Finally, including the boundary terms, the weak formulation of the pH system reads:

$$ egin{equation} 
\left\{ 
\begin{aligned}
(\partial_t h, \varphi)_{L^2} &= + (h e_p, \operatorname{grad}(\varphi))_{L^2} - (h u^0, \varphi)_{\partial\Omega}, \\
(h \partial_t p, \mathbf{\Phi})_{L^2} &= - (h \operatorname{grad}(e_h), \mathbf{\Phi})_{L^2} + \left( h \operatorname{curl}_{2D}(p) \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} e_p, \mathbf{\Phi} \right)_{L^2} \\
&\quad -2\mu (h \operatorname{Grad}(e_p), \operatorname{Grad}(\mathbf{\Phi}))_{\partial L^2} \\
&\quad -2\mu (h \operatorname{div}(e_p), \operatorname{div}(\mathbf{\Phi}))_{L^2} \\
&\quad + (h y^\mu, \mathbf{\Phi})_{\partial\Omega}.
\end{aligned}
\right. 
\tag{36}
\end{equation} $$

Remark 5 In order to obtain a finite-dimensional pH system, the integration by parts of the dissipative terms is mandatory (otherwise the induced matrices would not be symmetric and non-negative). On the contrary, the integration of $-(\operatorname{div}(he_p), \varphi)_{L^2}$ could have been replaced by the integration, on the second line, of the term $-(h \operatorname{grad}(e_h), \mathbf{\Phi})_{L^2}$. Numerical tests have shown a better behavior (e.g. on stability) for the former choice. Indeed, consider for instance the inviscid case, i.e., $\mu = 0$, then the formulation with integration by parts performed on the second line implies that the normal velocity control is imposed by an algebraic constraint, the output being the associated Lagrange multiplier. However, the latter integration by parts leads to some simplification as $y^\mu$ is then replaced by the whole output vector $y$. Further investigations are required to discriminate the two possibilities.

The weak form of the constitutive relations (22) may be written, making use of the appropriate weighted inner product introduced in Section 2.1.2:

$$ 
\left\{ 
\begin{aligned}
(e_h, \varphi)_{L^2} &= (\rho g h, \varphi)_{L^2} + \left( \frac{\|p\|^2}{2\rho}, \varphi \right)_{L^2}, \\
(h e_p, \mathbf{\Phi})_{L^2} &= \left( h \frac{p}{\rho}, \mathbf{\Phi} \right)_{L^2}.
\end{aligned}
\right. 
\tag{37} $$

To conclude to a pH system, it remains to take into account the power-conjugated boundary quantities, namely those related to $u^0$ and $y^\mu$. Let $\psi = \begin{pmatrix} \psi_n \\ \psi_t \end{pmatrix}$ be regular test functions at the boundary $\partial\Omega$, which are considered in the local basis formed by $(n, t)$.

The definition of the scalar observation $y^0$, given in (35), together with an $h$-dependent $L^2$ inner product, leads to:

$$ (h\psi_n, y^0)_{\partial\Omega} = -(h\psi_n, e_h)_{\partial\Omega}, $$

which turns out to be, as expected, the counterpart of the control term $(h u^0, \varphi)_{\partial\Omega}$ of (36).

On the other line of (36), the observation $y^\mu$ appears to be the Lagrange multiplier of the constraint prescribing the control in (32), which one writes
again with h embedded in:

$$ (h\psi, u)_{\partial\Omega} = (h\psi, e_p)_{\partial\Omega}. $$

Hence, the weak formulation at the boundary reads:

$$ \left\{ \begin{array}{l} (h\psi_n, y^0)_{\partial\Omega} = -(h\psi_n, e_h)_{\partial\Omega}, \\ (h\psi, u)_{\partial\Omega} = (h\psi, e_p)_{\partial\Omega}. \end{array} \right. \qquad (38) $$

Then, the full pH system in weak form is given by the dynamical equations (36), the constitutive relations (37) and the boundary ports (38).

Remark 6 Finally, the compatibility condition $u^0 = u \cdot n$ must be fulfilled, for the above system of equations to be well-defined.

## 3.2 Galerkin approximation

Consider two finite elements families ($\varphi^i$)$_{1\le i\le N_h}$ and ($\Phi^k$)$_{1\le k\le N_p}$ defined on $\Omega$. Let us approximate the energy and co-energy variables in these families:

$$ h(t, \boldsymbol{x}) \approx h^d(t, \boldsymbol{x}) := \sum_{i=1}^{N_h} h^i(t) \varphi^i(\boldsymbol{x}) = \varphi_h(\boldsymbol{x})^\top \underline{h}(t), $$

$$ e_h(t, \boldsymbol{x}) \approx e_h^d(t, \boldsymbol{x}) := \sum_{i=1}^{N_h} e_h^i(t) \varphi^i(\boldsymbol{x}) = \varphi_h(\boldsymbol{x})^\top \underline{e}_h(t), $$

$$ p(t, \boldsymbol{x}) \approx p^d(t, \boldsymbol{x}) := \sum_{k=1}^{N_p} p^k(t) \Phi^k(\boldsymbol{x}) = \Phi_p(\boldsymbol{x})^\top \underline{p}(t), $$

$$ e_p(t, \boldsymbol{x}) \approx e_p^d(t, \boldsymbol{x}) := \sum_{k=1}^{N_p} e_p^k(t) \Phi^k(\boldsymbol{x}) = \Phi_p(\boldsymbol{x})^\top \underline{e}_p(t). $$

Then, consider two finite elements families ($\psi_n^m$)$_{1\le m\le N_\partial}$ and ($\psi_t^m$)$_{1\le m\le N_\partial}$ defined at the boundary $\partial\Omega$, and let us approximate the boundary variables in these families as follows:

$$ \begin{aligned} \boldsymbol{u}(t, \boldsymbol{s}) &\approx \boldsymbol{u}^d(t, \boldsymbol{s}) := \sum_{m=1}^{N_\partial} u_n^m(t) \begin{pmatrix} \psi_n^m(\boldsymbol{s}) \\ 0 \end{pmatrix} + \sum_{m=1}^{N_\partial} u_t^m(t) \begin{pmatrix} 0 \\ \psi_t^m(\boldsymbol{s}) \end{pmatrix} \\ &= \psi_n(\boldsymbol{s})^\top \underline{u}_n(t) + \psi_t(\boldsymbol{s})^\top \underline{u}_t(t) = \boldsymbol{\Psi}(\boldsymbol{s})^\top \underline{u}(t), \end{aligned} $$

and its colocated output:

$$ \begin{aligned} \boldsymbol{y}^\mu(t, \boldsymbol{s}) &\approx \boldsymbol{y}^{\mu,d}(t, \boldsymbol{s}) := \sum_{m=1}^{N_\partial} y_n^{\mu,m}(t) \begin{pmatrix} \psi_n^m(\boldsymbol{s}) \\ 0 \end{pmatrix} + \sum_{m=1}^{N_\partial} y_t^{\mu,m}(t) \begin{pmatrix} 0 \\ \psi_t^m(\boldsymbol{s}) \end{pmatrix} \\ &= \psi_n(\boldsymbol{s})^\top \underline{y}_n^\mu(t) + \psi_t(\boldsymbol{s})^\top \underline{y}_t^\mu(t) = \boldsymbol{\Psi}(\boldsymbol{s})^\top \underline{y}^\mu(t). \end{aligned} $$
18 Shallow water equations with viscous damping

Furthermore, the other part of the boundary observation given in (35) has to be discretized in the normal family:

$$y^0(t, \mathbf{s}) \approx y^{0,d}(t, \mathbf{s}) := \sum_{m=1}^{N_{\partial}} y_n^{0,m}(t) \psi_n^m(\mathbf{s}) = \psi_n(\mathbf{s})^{\top} \underline{y}_n^0(t),$$

together with its colocated control, which will be denoted with the same exponent 0:

$$u^0(t, \mathbf{s}) \approx u^{0,d}(t, \mathbf{s}) := \sum_{m=1}^{N_{\partial}} u_n^{0,m}(t) \psi_n^m(\mathbf{s}) = \psi_n(\mathbf{s})^{\top} \underline{u}_n^0(t).$$

We may now focus back to Remark 6 and confirm that the extra control $u^0$ and the normal component $u_n$ must be set to the same value for compatibility. Therefore, numerically, the number of equations will be equal to the number of degrees of freedom, but physically, no additional control will be available.

Substituting the discrete counterpart of the variables and using the finite elements families as test functions in (36) gives rise to the following finite-dimensional system:

$$\left\{
\begin{array}{l}
\displaystyle M_h \frac{d}{dt} \underline{h} = D[\underline{h}] e_p(t) - B_n[\underline{h}] \underline{u}^0, \\[1.5ex]
\displaystyle M_p[\underline{h}] \frac{d}{dt} \underline{p} = -D[\underline{h}]^{\top} e_h + G[\underline{h}, \underline{p}] e_p \\
\displaystyle \qquad -2\mu \mathbf{R}_{\text{Grad}}[\underline{h}] e_p - 2\mu \mathbf{R}_{\text{div}}[\underline{h}] e_p + B_{\mu}[\underline{h}] \underline{y}^{\mu},
\end{array}
\right.
\qquad (39)$$

where the mass matrices are given by:

$$(M_h)_{i,j} := (\varphi^j, \varphi^i)_{L^2}, \quad (M_p[\underline{h}])_{k,l} := (\Phi^l, h^d \Phi^k)_{L^2},$$

the $D$ and $G$ matrices are defined by:

$$(D[\underline{h}])_{i,l} := (\Phi^l, h^d \mathrm{grad}(\varphi^i))_{L^2},$$

$$(G[\underline{h}, \underline{p}])_{k,l} := \left( \mathrm{curl}_{2D}(\underline{p}^d) \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} \underline{\Phi}^l, h^d \underline{\Phi}^k \right)_{L^2},$$

the symmetric non-negative dissipative matrices are obtained with:

$$(R_{\text{Grad}}[\underline{h}])_{k,l} := (\text{Grad}(\underline{\Phi}^l), h^d \text{Grad}(\underline{\Phi}^k))_{L^2},$$

$$(R_{\text{div}}[\underline{h}])_{k,l} := (\text{div}(\underline{\Phi}^l), h^d \text{div}(\underline{\Phi}^k))_{L^2},$$

and finally the control matrices are defined by:

$$(B_n[\underline{h}])_{i,n} := (\psi_n^n, h^d \varphi^i)_{\partial\Omega}, \quad (B_\mu[\underline{h}])_{k,m} := (\Psi^m, h^d \Phi^k)_{\partial\Omega}.$$
Regarding the constitutive relation, the system (37) becomes:

$$ \left\{ \begin{array}{l} M_h e_h = Q_h \underline{h} + N_h [\underline{p}] \underline{p}, \\ M_p [\underline{h}] e_p = Q_p [\underline{h}] \underline{p}, \end{array} \right. \qquad (40) $$

where the matrices on the right-hand side are given by:

$$ (Q_h)_{i,j} := (\varphi^j, \rho g \varphi^i)_{L^2}, \quad (Q_p[\underline{h}])_{k,l} := \left( \frac{\Phi^l}{\rho}, h^d \Phi^k \right)_{L^2}, $$

both square and symmetric, and:

$$ (N_h[\underline{p}])_{i,l} := \left( \frac{\underline{p}^d \cdot \Phi^l}{2\rho}, \varphi^i \right)_{L^2}, \quad \text{a rectangular matrix.} $$

Finally, the finite-dimensional counterpart of the colocated control and observation system (38) reads:

$$ \left\{ \begin{array}{l} M_n[\underline{h}] \underline{y}^0 = -B_n[\underline{h}]^\top e_h, \\ M_\partial[\underline{h}] \underline{u} = B_\mu[\underline{h}]^\top e_p, \end{array} \right. \qquad (41) $$

where the boundary mass matrices are:

$$ (M_n[\underline{h}])_{m,n} := (\psi_n^n, h^d \psi_n^m)_{\partial\Omega}, \quad (M_\partial[\underline{h}])_{m,n} := (\psi^m, h^d \psi^m)_{\partial\Omega}. $$

In compact form, the dynamics (39) together with the boundary variables (41) rewrite:

$$ \mathrm{Diag} \begin{bmatrix} M_h \\ M_p[\underline{h}] \\ M_n[\underline{h}] \\ M_\partial[\underline{h}] \end{bmatrix} \begin{pmatrix} \frac{d}{dt} \underline{h} \\ \frac{d}{dt} \underline{p} \\ -\underline{y}^0 \\ \underline{u} \end{pmatrix} = (\mathbf{J}[\underline{h}, \underline{p}] - \mathbf{R}[\underline{h}]) \begin{pmatrix} \frac{e_h}{\underline{u}} \\ \frac{e_p}{\underline{u}} \\ 0 \\ 0 \end{pmatrix}, \qquad (42) $$

where the extended structure matrix J is defined as follows:

$$ \mathbf{J}[\underline{h}, \underline{p}] := \begin{bmatrix} 0 & \mathbf{D}[\underline{h}] & -\mathbf{B}_n[\underline{h}] & 0 \\ -\mathbf{D}[\underline{h}]^\top & \mathbf{G}[\underline{h}, \underline{p}] & 0 & -\mathbf{B}_\mu[\underline{h}] \\ \mathbf{B}_n[\underline{h}]^\top & 0 & 0 & 0 \\ 0 & \mathbf{B}_\mu[\underline{h}]^\top & 0 & 0 \end{bmatrix}, $$

and the symmetric non-negative resistive matrix R is given by:

$$ \mathbf{R}[\underline{h}] := \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & \mathbf{R}_\mu[\underline{h}] & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}, \quad \mathbf{R}_\mu[\underline{h}] := 2\mu (\mathbf{R}_{\mathrm{Grad}}[\underline{h}] + \mathbf{R}_{\mathrm{div}}[\underline{h}]). $$
20 Shallow water equations with viscous damping

Remark 7 It is important to remember that in (42), $\underline{u}$ and $\underline{u}^0$ require the compatibility condition $\underline{u} \cdot \underline{n} = \underline{u}^0$.

The constitutive relations (40) may be written in compact form as:

$$ \begin{bmatrix} M_h & 0 \\ 0 & M_p[\underline{h}] \end{bmatrix} \begin{pmatrix} \underline{e}_h \\ \underline{e}_p \end{pmatrix} = \begin{bmatrix} Q_h & N_h[\underline{p}] \\ 0 & Q_p[\underline{h}] \end{bmatrix} \begin{pmatrix} \underline{h} \\ \underline{p} \end{pmatrix} \quad (43) $$

The final discrete pH system is composed of (42), preserving the underlying Dirac structure, and (43).

Remark 8 The discrete system (42)–(43) is a pH Differential Algebraic Equation (pH-DAE), see [47]. In the linear-quadratic case, making use of the so-called co-energy formulation is proven to be very efficient [32]; however, in our case, the non-linearities forbid the substitution of the constitutive relations (43) into the dynamics (42); furthermore, algebraic equations appear due to the arrangement of the boundary ports. For these reasons, the difficulties raised by DAEs should be carefully considered. Among them, compatibility conditions at the initial time are critical and will be carefully checked in Section 4. Another source of difficulty would be to use the formulation involving extra dissipation ports, some of them being symmetric tensors, see Section 2.2.2; thus, it has been preferred to work with the $J-R$ formulation, which is of much smaller dimension. The wide class of linear pH-DAEs is already challenging and constitutes an active research topic: in particular, the index of such DAEs can be at most 2 [48, Corollary 51]; the interested reader may refer to [38, 49, 50] and the many references therein.

## 3.3 Discrete power balance

One of the main features of a pH system being the power balance satisfied by the Hamiltonian, a discrete counterpart of Theorem 4 is expected. This preservation property carries over for linear systems with a quadratic Hamiltonian, while in the modulated case or nonlinear case, only a few attempts have already been made, see e.g., [51, 52], where the polynomial structure of the equations makes a difference in the proof of the energy preservation.

Let us define the discrete Hamiltonian as follows:

$$ \mathcal{H}^d(t) := \mathcal{H}(h^d(t, \boldsymbol{x}), \boldsymbol{p}^d(t, \boldsymbol{x})) = \frac{1}{2} \int_{\Omega} \left[ h^d(t, \boldsymbol{x}) \frac{\|\boldsymbol{p}^d(t, \boldsymbol{x})\|^2}{\rho} + \rho g(h^d(t, \boldsymbol{x}))^2 \right] d\Omega, $$

which may be rewritten as:

$$ \mathcal{H}^d(t) = \frac{1}{2} \underline{p}(t)^\top \underline{Q}_p[\underline{h}(t)] \underline{p}(t) + \frac{1}{2} \underline{h}(t)^\top \underline{Q}_h \underline{h}(t). \quad (44) $$

Before stating the power balance satisfied by the discrete Hamiltonian, let us enlighten a trivial but crucial equality.
**Lemma 1** The following holds true:

$$ \frac{1}{2} \underline{p}(t)^{\top} \frac{d}{dt} (\underline{Q}_{\underline{p}}[\underline{h}(t)]) \underline{p}(t) = \underline{p}(t)^{\top} \underline{N}_{h}[\underline{p}(t)]^{\top} \frac{d}{dt} \underline{h}(t). $$

*Proof* Indeed, one has:

$$ \begin{aligned} \frac{1}{2} \underline{p}(t)^{\top} \frac{d}{dt} (\underline{Q}_{\underline{p}}[\underline{h}(t)]) \underline{p}(t) &= \frac{1}{2} \left( \frac{\underline{p}^d}{\rho}, \partial_t h^d \underline{p}^d \right)_{L^2}, \\ &= \left( \frac{\|\underline{p}^d\|^2}{2\rho}, \partial_t h^d \right)_{L^2}, \\ &= \underline{p}(t)^{\top} \underline{N}_h[\underline{p}(t)]^{\top} \frac{d}{dt} \underline{h}(t). \end{aligned} $$

□

**Theorem 5** The discrete Hamiltonian satisfies the following discrete power balance:

$$ \begin{aligned} \frac{d}{dt} \mathcal{H}^d(t) = &-2\mu e_p(t)^{\top} (\mathbf{R}_{\text{Grad}}[\underline{h}(t)] + \mathbf{R}_{\text{div}}[\underline{h}(t)]) e_p(t) \\ &+ \underline{u}^0(t)^{\top} \mathbf{M}_n[\underline{h}(t)] \underline{y}^0(t) + \underline{y}^{\mu}(t)^{\top} \mathbf{M}_{\partial}[\underline{h}(t)] \underline{u}(t), \\ &\le \underline{u}^0(t)^{\top} \mathbf{M}_n[\underline{h}(t)] \underline{y}^0(t) + \underline{y}^{\mu}(t)^{\top} \mathbf{M}_{\partial}[\underline{h}(t)] \underline{u}(t). \end{aligned} $$

*Proof* On the one hand, one easily computes from the time-derivative of the discrete Hamiltonian (44):

$$ \frac{d}{dt} \mathcal{H}^d(t) = \underline{p}(t)^{\top} \underline{Q}_{\underline{p}}[\underline{h}(t)] \frac{d}{dt} \underline{p}(t) + \underline{h}(t)^{\top} \underline{Q}_{\underline{h}} \frac{d}{dt} \underline{h}(t) + \frac{1}{2} \underline{p}(t)^{\top} \frac{d}{dt} (\underline{Q}_{\underline{p}}[\underline{h}(t)]) \underline{p}(t), \quad (45) $$

On the other hand, multiplying (42) by $\begin{pmatrix} \frac{e_h(t)}{e_p(t)} \\ \frac{\underline{u}^0(t)}{\underline{u}^{\mu}(t)} \\ -\underline{y}^{\mu}(t) \end{pmatrix}$ by the left side leads to (since $\underline{J}[\underline{h}(t), \underline{p}(t)]$ is skew-symmetric):

$$ \begin{aligned} \underline{e}_h(t)^{\top} \mathbf{M}_h \frac{d}{dt} \underline{h}(t) + \underline{e}_p(t)^{\top} \mathbf{M}_p[\underline{h}(t)] \frac{d}{dt} \underline{p}(t) &- \underline{u}^0(t)^{\top} \mathbf{M}_n[\underline{h}(t)] \underline{y}^0(t) - \underline{y}^{\mu}(t)^{\top} \mathbf{M}_{\partial}[\underline{h}(t)] \underline{u}(t) = - \underline{e}_p(t)^{\top} \mathbf{R}_{\mu}[\underline{h}(t)] \underline{e}_p(t). \end{aligned} $$

By symmetry of the mass matrices, (43) implies:

$$ \begin{aligned} \left( \underline{h}(t)^{\top} \underline{Q}_{\underline{h}}^{\top} + \underline{p}(t)^{\top} \underline{N}_{h}[\underline{p}(t)]^{\top} \right) \frac{d}{dt} \underline{h}(t) + \underline{p}(t)^{\top} \underline{Q}_{\underline{p}}[\underline{h}(t)]^{\top} \frac{d}{dt} \underline{p}(t) &- \underline{u}^0(t)^{\top} \mathbf{M}_n[\underline{h}(t)] \underline{y}^0(t) - \underline{y}^{\mu}(t)^{\top} \mathbf{M}_{\partial}[\underline{h}(t)] \underline{u}(t) = - \underline{e}_p(t)^{\top} \mathbf{R}_{\mu}[\underline{h}(t)] \underline{e}_p(t). \end{aligned} $$

Rearranging the terms and using Lemma 1, one gets:

$$ \begin{aligned} \underline{p}(t)^{\top} \underline{Q}_{\underline{p}}[\underline{h}(t)]^{\top} \frac{d}{dt} \underline{p}(t) + \underline{h}(t)^{\top} \underline{Q}_{\underline{h}}^{\top} \frac{d}{dt} \underline{h}(t) + \frac{1}{2} \underline{p}(t)^{\top} \frac{d}{dt} (\underline{Q}_{\underline{p}}[\underline{h}(t)]) \underline{p}(t) &- \underline{e}_p(t)^{\top} \mathbf{R}_{\mu}[\underline{h}(t)] \underline{e}_p(t) + \underline{u}^0(t)^{\top} \mathbf{M}_n[\underline{h}(t)] \underline{y}^0(t) + \underline{y}^{\mu}(t)^{\top} \mathbf{M}_{\partial}[\underline{h}(t)] \underline{u}(t) = 0. \end{aligned} $$

Substituting (45) and the definition of the matrix $\mathbf{R}_{\mu}[\underline{h}(t)]$ in the latter gives the result.

□
# 4 Numerical results

This section is devoted to simulation tests in order to show the feasibility of our approach. Note that the difficult question of time integration of the resulting non-linear DAE is not treated in this work, where the Backward Differentiation Formula (BDF) of order 2 provided by the PETSc library [53] has been used, with a maximal time step $\Delta t = 0.005$. In particular, the fully discrete system is not a discrete pH system *a priori*, implying that the discrete power balance is only satisfied up to the time scheme error.

Videos of the simulations will be available on the website of the journal as supplementary materials. The codes are available as *published examples* on the github https://github.com/g-haine/scrimp of the SCRIMP environment used for the simulations [31].

For the following tests, continuous Lagrange finite elements have been chosen to approximate all quantities (both scalar-fields and vector-fields). Indeed, we assume a smooth solution in the previous sections, and the use of discontinuous finite elements would require further analysis. Moreover, the definition of the matrix $D$ involves the gradient of $\varphi^i$, while $R$ involves first order derivatives of $\Phi$ (except for the inviscid case where $R \equiv 0$). Therefore, an efficient approach is to take the same order $k$ for both $h$-type (scalar) fields and $p$-type (vectorial) fields, when dissipation occurs, and an order $k+1$ for the $h$-type variables otherwise. Regarding the boundary fields, to ensure compatibility with the traces of in-domain fields, they are approximated with Lagrange elements of order $k$. In what follows, we present the results obtained for $k = 2$.

## 4.1 Closed system: the sloshing problem

Let us consider a rectangular tank of size $(0, 2) \times (0, 0.5)$ filled with a liquid of mass density $\rho = 1 \text{ kg m}^{-3}$, lying in a gravity field$^1 g = 0.01 \text{ m s}^{-2}$.

### 4.1.1 The inviscid case

In this first example, we aim to validate our approach by recovering known behaviors in 1D.

Since there is no viscosity, the boundary control that is directly available in the model is the normal velocity $u_n$, as shown in Section 3.1. In particular, the vector control $u$ from (32) and the Lagrange multiplier $y^\mu$ from (35) are removed from the model in this example. Taking $u_n = 0$ and an adequate initial data, the solution will only vary in the $x$ direction. More precisely, we choose a continuous piecewise linear function that is constant equal to $h_0 = 55 \text{ m}$ before some $x_0$ and constant equal to $h_1 = 50 \text{ m}$ after some $x_1 > x_0$, and a null initial velocity. With these values at hand, one may compute at which velocity the front wave will propagate, see e.g., [54, Chapter 10]: $\sqrt{gh_0} \approx \sqrt{0.55} \approx$

<footnote>$^1$These values have been chosen in order to be able to use a relatively large time step w.r.t. the mesh size. Indeed, with $\rho = 1000 \text{ kg m}^{-3}$ and $g = 10 \text{ m s}^{-2}$, the characteristic velocity would be about $100 \text{ m s}^{-1}$, which would require a very small time step to be approximated in a tank of length 2 m. Nevertheless, we choose to keep a ratio of 100 between these two parameters.</footnote>
0.74 m s⁻¹, which is recovered numerically, as seen on Fig. 1. Comparing the second and third plots of Fig. 1, corresponding to time t = 0.95 s and t = 1.85 s, one may appreciate how we recover this characteristic speed.

### 4.1.2 Homogeneous boundary velocity control

We now add a dissipation, represented by a viscosity $\mu = 0.001 \text{ kg m}^{-1} \text{ s}^{-1}$, as introduced in Section 2.2. It is then possible to access the boundary term $\mathbf{u}$, and in particular the tangent component of the velocity $u_t$. In addition to $u_n = 0$, let us impose furthermore $u_t = 0$ in this test case.

It has been seen in Section 3.1 that $\mathbf{u}$ can be prescribed in an implicit manner: the boundary term appearing after integrations by parts is $y^\mu$, which turns out to be the Lagrange multiplier of the algebraic constraint (41) used to prescribe $\mathbf{u}$. It has to be noted that the use of an extended system with dissipation ports, as presented in Section 2.2.2, instead of the $J-R$ approach would have made possible to prescribe directly the velocity at the boundary, but to the price of a larger finite-dimensional system.

While the inviscid case ($\mu = 0$) should give a constant Hamiltonian, it can be observed on the top of Fig. 2 that it varies only at time scheme precision. The viscous case should show a loss of energy, that is indeed observable on the bottom of Fig. 2, a decay which shows proportional to the value of $\mu$ in our tests.

The tangent velocity is correctly prescribed, as can be seen on the first plot of Fig. 3, and the vorticity indeed develops near the boundaries as expected (see the second plot of Fig. 3).

## 4.2 Emptying a tank: normal control

In this test, we applied an outgoing flow at the right side of the rectangle, to simulate the emptying of the tank by opening of a sluice gate at the right end. The initial height is constant at 50m, with a null velocity field. The control is a parabolic profile that slowly increases as $t$ goes to 1, and constant after this instant. More precisely, the normal control is $u_n(t, y) = 0.1y(y - 0.5)\min(1, t)$ at the right-end side of the rectangle, and null elsewhere. The tangent control is still imposed at zero: $u_t = 0$. The parabolic profile of $u_n$ at $x = 2$ guarantees the compatibility with the null tangent velocity prescribed on both the top and the bottom boundaries of the tank.

One may appreciate the expected behavior of the height of fluid on Fig. 4. In particular, the loss of potential energy due to the outgoing boundary flow is clearly the dominant effect in the Hamiltonian decay.

## 4.3 Rotating circular tank: tangential control

In this last numerical experiment, we consider a circular tank that is rotating, as an example of tangential velocity control. The velocity field is initialized with a rotating field, compatible with the boundary control, while the height of the fluid starts at a constant value, see the first plot of Fig. 5.
24 Shallow water equations with viscous damping

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Evolution of the height and velocity of an inviscid fluid contained in a tank at different time steps (t=0 s, t=1.02 s, t=2.09 s, t=3.09 s).</figcaption></figure>
At first, by the centrifugal effect, one may observe the fluid is pushed against the boundary, as shown on the middle plot of Fig. 5 at $t = 1.01$ s. Second, it bounces at the boundary and goes back to generate a spike at the center of the tank, as seen on the last plot of Fig. 5 at $t = 1.77$ s.

Starting at $t = 100$ s, the tangential boundary control begins to reduce, reaching zero near $t = 250$ s, as depicted in the first plot of Fig. 6.

Regarding the time evolution of the energies in this last experiment, one may observe on Fig. 6 that, despite small oscillations of the total energy in the beginning of the simulation, they quickly stabilize and the total energy behaves as expected, even under continuous external energy supply and dissipation. A long-time simulation (with final time $t_f = 500$ s) has been performed, showing excellent global behavior (see the lower part of Fig. 6). One might notice that the spurious oscillations that occur in the beginning of the simulation coincide with high-amplitude height oscillations, leading to a strong nonlinear behavior of the Shallow Water Equations (see (42)), which would require some dedicated time-stepping schemes.

Finally, on Fig. 7, one can see the effect of the decay of the boundary control at $t = 250$ s (top); this decay then propagates inside the circular domain at $t = 300$ s (middle); and eventually, stabilization is achieved around $t = 375$ s (bottom): the kinetic energy (dashed blue curve on Fig. 6) indeed decays to zero, while the potential energy (solid blue curve on Fig. 6) stabilizes at its minimum value. One may appreciate the constant value of the supplied energy $\mathbf{S}u$ (dashed purple curve on Fig. 6) when the power supplied at the boundary through $\mathbf{u}$ reaches zero. Regarding the dissipation, it mainly occurs via $\mathcal{D}_{\text{Grad}}$ (solid green curve on Fig. 6), which is consistent with the chosen geometry, initial data and applied control.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Variation of the different energies in play (homogeneous boundary velocity case): in the inviscid case (top) and the viscous case (bottom). The notation $Su^0$ represents the supplied energy provided by the normal velocity control $u^0$, $Su$ the supplied energy provided by the vectorial control $u$, $D_{Grad}$ the dissipated energy via the matrix $R_{Grad}[\bar{h}]$ and $D_{div}$ via the matrix $R_{div}[\bar{h}]$, see (39).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 3** (Top) Height and velocity of a viscous fluid in a tank at time $t = 1.02$ s, with homogeneous velocity boundary control. (Bottom) Vorticity at time $t = 1.02$ s.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4** (Top) Variation of the different energies in play (emptying case). (Bottom) Height and velocity profiles at time $t = 2.09$ s.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Height of fluid and velocity field of a rotating circular tank at different times *t*: on top *t* = 0 s, in the middle *t* = 1.01 s, at bottom *t* = 1.77 s.</figcaption></figure>
30

Shallow water equations with viscous damping

<figure><img src="image_9.png" /><figcaption>**Fig. 6** Time evolution of the tangential control (top) and variation of the different energies (bottom) for a long-time simulation.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Height of fluid and velocity field of a rotating circular tank at different times t: on top t = 250 s, in the middle t = 300 s, at bottom t = 375 s.</figcaption></figure>
# 5 Conclusions and further work

Firstly, this paper presented a pH formulation for the 2D rotational SWE with viscous damping. In previous work of SWE using pH formulations, only the normal velocity is taken into account as boundary control. With the inclusion of viscous effects, this framework allows to explicitly take into account the tangential velocity as boundary control port.

Secondly, using a partitioned finite element method, we achieve a finite-
dimensional continuous-time pH approximation that preserves the system's
structure. Our method ensures the incorporation of normal/tangential bound-
ary control ports for the 2D rotational SWE with viscous damping after
discretization.

Finally, four simulation scenarios were investigated to illustrate the approach and proved to be encouraging. However, a careful attention to the time-marching numerical scheme should be addressed in further work, for instance, through the study of energy-preserving time-steppers. This issue is especially important to tackle severe nonlinearities when the Froude number becomes greater than one, i.e. for supercritical flow, see e.g. [54]. Furthermore, the relation of compatibility condition between the finite element families should be explored, e.g., adjusting the order of the elements.

Further work should also address the use of the models developed here for control design, especially by taking advantage of the availability of the tangential boundary control port.

# 6 Declarations

Ethical Approval

Not applicable.

Competing interests

The authors declare that they have no competing interests as defined by
Springer, or other interests that might be perceived to influence the results
and/or discussion reported in this paper.

Authors’ contributions

The authors contributed equally to this work.

Funding

This study was funded in part by Agence Nationale de la Recherche (France), project ANR-21-CE48-0018: IMPACTS, and by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - CAPES (Brazil), project 88887.717002/2022-00.
Acknowledgments

We sincerely thank both reviewers for their thorough and thoughtful reviews. Their constructive feedback and insightful suggestions, both on the manuscript and the numerical simulations, have significantly improved the quality of this paper. We deeply appreciate their dedication and effort.

Availability of data and materials

The simulation codes are available at: https://github.com/g-haine/scrimp/.

References

[1] Robles, R.R., Serrano, J.P.: Sloshing mechanical model for stability and handling qualities evaluation of the C295 aircraft with the OSD system. In: Proceedings of 29th Congress of the International Council of the Aeronautical Sciences, St. Petersburg, Russia (2014)

[2] Cho, Y.-S., Sohn, D.-H., Lee, S.O.: Practical modified scheme of linear shallow-water equations for distant propagation of tsunamis. Ocean Engineering **34**(11-12), 1769–1777 (2007)

[3] Clifford, M., Horton, C., Schmitz, J.: Swafs: Shallow water analysis and forecast system. In: Proceedings of OCEANS’94, vol. 3, p. 82 (1994). IEEE

[4] Barati, R., Rahimi, S., Akbari, G.H.: Analysis of dynamic wave model for flood routing in natural rivers. Water Science and Engineering **5**(3), 243–258 (2012)

[5] Gilman, P.A.: Magnetohydrodynamic “shallow water” equations for the solar tachocline. The Astrophysical Journal **544**(1), 79 (2000)

[6] Zeitlin, V.: Geophysical Fluid Dynamics: Understanding (almost) Everything with Rotating Shallow Water Models. Oxford University Press, Oxford (2018)

[7] Fraccarollo, L., Toro, E.F.: Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems. Journal of hydraulic research **33**(6), 843–864 (1995)

[8] Cardoso-Ribeiro, F.L., Brugnoli, A., Matignon, D., Lefèvre, L.: Port-Hamiltonian modeling, discretization and feedback control of a circular water tank. In: 2019 IEEE 58th Conference on Decision and Control (CDC), pp. 6881–6886 (2019). https://doi.org/10.1109/CDC40024.2019.9030007
[9] Noble, D.R., Davey, T., Smith, H.C., Kaklis, P., Robinson, A., Bruce, T.: Spatial variation in currents generated in the flowave ocean energy research facility. In: Proceedings of the 11th European Wave and Tidal Energy Conference (EWTEC2015), pp. 1–8 (2015). Nantes France

[10] Kirstetter, G., Hu, J., Delestre, O., Darboux, F., Lagrée, P.-Y., Popinet, S., Fullana, J.M., Josserand, C.: Modeling rain-driven overland flow: Empirical versus analytical friction terms in the shallow water approximation. Journal of Hydrology **536**, 1–9 (2016). https://doi.org/10.1016/j.jhydrol.2016.02.022

[11] James, F., Lagrée, P.-Y., Le, M.H., Legrand, M.: Towards a new friction model for shallow water equations through an interactive viscous layer. ESAIM: M2AN **53**(1), 269–299 (2019). https://doi.org/10.1051/m2an/2018076

[12] Lagrée, P.-Y.: Equations de Saint Venant et application aux mouvements de fonds érodables, Paris, France. Class notes (2021). http://www.lmm.jussieu.fr/%7Elagree/COURS/MFEnv/MFEnv.pdf

[13] Gerbeau, J.-F., Perthame, B.: Derivation of viscous Saint-Venant system for laminar shallow water; numerical validation. Discrete and Continuous Dynamical Systems - B **1**(1), 89–102 (2001). https://doi.org/10.3934/dcdsb.2001.1.89

[14] Marche, F.: Derivation of a new two-dimensional viscous shallow water model with varying topography, bottom friction and capillary effects. European Journal of Mechanics B/Fluids **26**, 49–63 (2007). https://doi.org/10.1016/j.euromechflu.2006.04.007

[15] Sundbye, L.: Global existence for the cauchy problem for the viscous shallow water equations. The Rocky Mountain Journal of Mathematics **28**(3), 1135–1152 (1998)

[16] Duindam, V., Macchelli, A., Stramigioli, S., Bruyninckx, H.: Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach, p. 442. Springer, Berlin, Heidelberg (2009). https://doi.org/10.1007/978-3-642-03196-0

[17] van der Schaft, A.J., Jeltssema, D.: Port-Hamiltonian Systems Theory: An Introductory Overview. Foundations and Trends® in Systems and Control **1**(2), 173–378 (2014). https://doi.org/10.1561/2600000002

[18] Rashad, R., Califano, F., van der Schaft, A.J., Stramigioli, S.: Twenty years of distributed port-Hamiltonian systems: a literature review. IMA Journal of Mathematical Control and Information **37**(4), 1400–1422 (2020)
[19] Hamroun, B., Lefèvre, L., Mendes, E.: Port-based modelling and geometric reduction for open channel irrigation systems. In: Proceedings of the 46th IEEE Conference on Decision and Control, pp. 1578–1583. IEEE, New Orleans, LA, USA (2007). https://doi.org/10.1109/CDC.2007.4434237

[20] Pasumarthy, R., Ambati, V.R., van der Schaft, A.J.: Port-Hamiltonian discretization for open channel flows. Systems & Control Letters 61(9), 950–958 (2012)

[21] Hamroun, B., Dimofte, A., Lefèvre, L., Mendes, E.: Control by Interconnection and Energy-Shaping Methods of Port Hamiltonian Models. Application to the Shallow Water Equations. European Journal of Control 16(5), 545–563 (2010). https://doi.org/10.3166/ejc.16.545-563

[22] Cardoso-Ribeiro, F.L., Matignon, D., Pommier-Budinger, V.: A port-Hamiltonian model of liquid sloshing in moving containers and application to a fluid-structure system. Journal of Fluids and Structures 69, 402–427 (2017). https://doi.org/10.1016/j.jfluidstructs.2016.12.007

[23] Cardoso-Ribeiro, F.L., Matignon, D., Pommier-Budinger, V.: Port-Hamiltonian model of two-dimensional shallow water equations in moving containers. IMA Journal of Mathematical Control and Information 37(4), 1348–1366 (2020). https://doi.org/10.1093/imamci/dnaa016

[24] Califano, F., Rashad, R., Schuller, F.P., Stramigioli, S.: Geometric and energy-aware decomposition of the Navier-Stokes equations: A port-Hamiltonian approach. Physics of Fluids 33(4), 047114 (2021). https://doi.org/10.1063/5.0048359

[25] Cardoso-Ribeiro, F.L., Matignon, D., Lefèvre, L.: Dissipative Shallow Water equations: a port-Hamiltonian formulation. IFAC-PapersOnLine 54(19), 167–172 (2021). https://doi.org/10.1016/j.ifacol.2021.11.073. 7th IFAC Workshop on Lagrangian and Hamiltonian Methods for Nonlinear Control LHMNC 2021

[26] Limache, A., Sánchez, P.J., Dalcín, L., Idelsohn, S.: Objectivity tests for Navier–Stokes simulations: The revealing of non-physical solutions produced by Laplace formulations. Comput. Methods Appl. Mech. Engrg. 197, 4180–4192 (2008). https://doi.org/10.1016/j.cma.2008.04.020

[27] Leborgne, G.: Objectivity in continuum mechanics, an introduction. motions, Eulerian and Lagrangian variables and functions, deformation gradient, Lie derivatives, velocity-addition formula, Coriolis. arXiv preprint arXiv:2301.01056 (2023)
[28] Bossavit, A.: Computational Electromagnetism: Variational Formula-
tions, Complementarity, Edge Elements. Academic Press, San Diego
(1998). https://doi.org/10.1016/B978-0-12-118710-1.X5000-4

[29] Kotyczka, P., Maschke, B., Lefèvre, L.: Weak form of Stokes–Dirac struc-
tures and geometric discretization of port-Hamiltonian systems. Journal
of Computational Physics **361**, 442–476 (2018)

[30] Cardoso-Ribeiro, F.L., Matignon, D., Lefèvre, L.: A Partitioned Finite-
Element Method for power-preserving discretization of open systems of
conservation laws. IMA J. Mathematical Control and Information **38**(2),
493–533 (2021). https://doi.org/10.1093/imamci/dnaa038

[31] Ferraro, G., Fournié, M., Haine, G.: Simulation and control of interac-
tions in multi-physics, a Python package for port-Hamiltonian systems.
IFAC-PapersOnLine **58**(6), 119–124 (2024). https://doi.org/10.1016/j.
ifacol.2024.08.267. 8th IFAC Workshop on Lagrangian and Hamiltonian
Methods for Nonlinear Control LHMNC 2024

[32] Haine, G., Matignon, D., Serhani, A.: Numerical analysis of a structure-
preserving space-discretization for an anisotropic and heterogeneous
boundary controlled $N$-dimensional wave equation as a port-Hamiltonian
system. Int. J. Numer. Anal. Mod. **20**(1), 92–133 (2023). https://doi.org/
10.4208/ijnam2023-1005

[33] Boyer, F., Fabrie, P.: Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models. Applied Mathematical Sciences, vol. 183, p. 525. Springer, New York (2013)

[34] Cotter, C.J.: Compatible finite element methods for geophysical fluid dynamics. Acta Numerica **32**, 291–393 (2023). https://doi.org/10.1017/S0962492923000028

[35] Arnold, D.N., Falk, R.S., Winther, R.: Finite element exterior calcu-
lus, homological techniques, and applications. Acta Numerica **15**, 1–155
(2006). https://doi.org/10.1017/S0962492906210018

[36] Olver, P.J.: Applications of Lie Groups to Differential Equations, 2nd edn.
Graduate Texts in Mathematics, vol. 107. Springer, New York (1993)

[37] van der Schaft, A.J., Maschke, B.M.: Hamiltonian formulation of distributed-parameter systems with boundary energy flow. Journal of Geometry and Physics **42**, 166–194 (2002)

[38] Mehrmann, V., Unger, B.: Control of port-Hamiltonian differential-algebraic systems and applications. Acta Numerica **32**, 395–515 (2023). https://doi.org/10.1017/S0962492922000083
[39] Altmann, R., Schulze, P.: A port-Hamiltonian formulation of the Navier-Stokes equations for reactive flows. Systems & Control Letters 100, 51-55 (2017). https://doi.org/10.1016/j.sysconle.2016.12.005

[40] Chorin, A.J., Marsden, J.E.: A Mathematical Introduction to Fluid Mechanics. Texts in Applied Mathematics, vol. 4. Springer, New York (1992)

[41] Moses Badlyan, A., Maschke, B., Beattie, C.A., Mehrmann, V.: Open Physical Systems: From GENERIC to Port-Hamiltonian Systems. In: Proc. of the International Symposium on Mathematical Theory of Networks and Systems, Hong Kong, China, pp. 204-211 (2018)

[42] Brugnoli, A., Alazard, D., Pommier-Budinger, V., Matignon, D.: Port-Hamiltonian formulation and symplectic discretization of plate models Part I: Mindlin model for thick plates. Applied Mathematical Modelling 75, 940-960 (2019). https://doi.org/10.1016/j.apm.2019.04.035

[43] Mora, L.A., Le Gorrec, Y., Matignon, D., Ramirez, H., Yuz, J.I.: On port-Hamiltonian formulations of 3-dimensional compressible Newtonian fluids. Physics of Fluids 33(11), 117117 (2021). https://doi.org/10.1063/5.0067784

[44] Sundbye, L.: Global existence for the Dirichlet problem for the viscous shallow water equations. Journal of Mathematical Analysis and Applications 202(1), 236-258 (1996). https://doi.org/10.1006/jmaa.1996.0315

[45] Kloeden, P.E.: Global existence of classical solutions in the dissipative shallow water equations. SIAM Journal on Mathematical Analysis 16(2), 301-315 (1985). https://doi.org/10.1137/0516022

[46] Serhani, A., Matignon, D., Haine, G.: A partitioned finite element method for the structure-preserving discretization of damped infinite-dimensional port-Hamiltonian systems with boundary control. In: Geometric Science of Information, pp. 549-558. Springer, Cham (2019). https://doi.org/10.1007/978-3-030-26980-7_57

[47] van der Schaft, A.J.: Port-Hamiltonian Differential-Algebraic Systems, pp. 173-226. Springer, Berlin, Heidelberg (2013). https://doi.org/10.1007/978-3-642-34928-7_5

[48] Mehrmann, V., van der Schaft, A.: Differential-algebraic systems with dissipative Hamiltonian structure. Math. Control Signals Syst. 35, 541-584 (2023). https://doi.org/10.1007/s00498-023-00349-2

[49] Beattie, C., Mehrmann, V., Xu, H., Zwart, H.: Linear port-Hamiltonian
descriptor systems. Math. Control Signals Syst. **30**, 17 (2018). https://doi.org/10.1007/s00498-018-0223-3

[50] van der Schaft, A., Mehrmann, V.: Linear port-Hamiltonian DAE systems revisited. Systems & Control Letters **177**, 105564 (2023). https://doi.org/10.1016/j.sysconle.2023.105564

[51] Bendimerad-Hohl, A., Haine, G., Matignon, D., Maschke, B.: Structure-preserving discretization of a coupled Allen-Cahn and heat equation system. IFAC-PapersOnLine **55**(18), 99-104 (2022). https://doi.org/10.1016/j.ifacol.2022.08.037

[52] Bendimerad-Hohl, A., Haine, G., Matignon, D.: Structure-preserving discretization of the Cahn-Hilliard equations recast as a port-Hamiltonian system. In: Geometric Science of Information, pp. 192–201. Springer, Cham (2023). https://doi.org/10.1007/978-3-031-38299-4_21

[53] Balay, S., Abhyankar, S., Adams, M.F., Benson, S., Brown, J., Brune, P., Buschelman, K., Constantinescu, E., Dalcin, L., Dener, A., Eijkhout, V., Faibussowitsch, J., Gropp, W.D., Hapla, V., Isaac, T., Jolivet, P., Karpeev, D., Kaushik, D., Knepley, M.G., Kong, F., Kruger, S., May, D.A., McInnes, L.C., Mills, R.T., Mitchell, L., Munson, T., Roman, J.E., Rupp, K., Sanan, P., Sarich, J., Smith, B.F., Zampini, S., Zhang, H., Zhang, H., Zhang, J.: PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.20, Argonne National Laboratory (2023). https://doi.org/10.2172/1968587

[54] Stoker, J.J.: Water Waves. John Wiley & Sons, Ltd, New York (1992). https://doi.org/10.1002/9781118033159