# An ALE Numerical Method for the Waterline Evolution of a Floating Object with Non-linear Shallow-water Equations

Ali Haidar, Fabien Marche, François Vilar

► To cite this version:

Ali Haidar, Fabien Marche, François Vilar. An ALE Numerical Method for the Waterline Evolution of a Floating Object with Non-linear Shallow-water Equations. Journal of Scientific Computing, 2025, 102 (2), pp.53. 10.1007/s10915-024-02778-w. hal-05171576

HAL Id: hal-05171576

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

Submitted on 20 Jul 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An ALE numerical method for the waterline evolution of a floating object with non-linear shallow-water equations

Ali Haidar · Fabien Marche · François Vilar

**Abstract** In this work, a numerical method is introduced for the study of nonlinear interactions between free-surface shallow-water flows and a partly immersed floating object. At the continuous level, the fluid's evolution is modeled with the nonlinear hyperbolic shallow-water equations. The description of the flow beneath the object reduces to an algebraic and nonlinear equation for the free-surface, together with a nonlinear differential equation for the discharge. The object's motion may be either prescribed, or computed as a response to the hydrodynamic forcing. In the later case, with heaving, surging and pitching allowed in the horizontal one-dimensional case, these equations are supplemented with the Newton's second law for the object's motion, involving the force and torque applied by the surrounding fluid, and parts of this external forcing are regarded as an added-mass effect. At the discrete level, we introduce a discontinuous Galerkin approximation, stabilized by a recent *a posteriori* Local Subcell Correction method in the vicinity of the solution's singularities and loss of admissibility. The motion of the fluid-structure contact-points is described with an Arbitrary-Lagrangian-Eulerian strategy, resulting in a global algorithm that ensures the preservation of the water-height positivity at the sub-cell level, preserves the class of motionless steady-states even when the object is allowed to evolve freely, and the Discrete Geometric Conservation Law. Several numerical computations involving wave and floating object interactions are provided, showing the robust computation of the air-water-body contact-points dynamics.

**Keywords** shallow-water equations · discontinuous Galerkin · *a posteriori* limitation · Arbitrary-Lagrangian-Eulerian · floating body · wave-structure interactions

## 1 Introduction

The modeling of a free-surface water domain in the presence of a floating object is an intricate problem, combining the time evolution of an incompressible fluid with the motion of a solid surface-piercing and partly immersed object. Possible applications may include, for instance, the study of offshore structures, near-ship flows, floating breakwaters or floating wave-energy converters. From a mathematical viewpoint, such a problem may fall into the family of *free-boundary* problems, as the immersed part of the object (*the wetted surface*) depends on time and the boundary of the free-surface water domain identifies with the water-object contact-line. In early works [28,29], this issue is overcome by assuming small-amplitude motions, neglecting the time-evolution of the wetted surface and using a linear potential model for the fluid evolution. Among seminal works, let also mention [21,39] where surface-waves generation through prescribed motion of floating or submerged bodies are studied. As of today, dedicated software generally rely on Boundary Element Methods (BEM) applied to linear models in the time-domain, like the Cummins equation [10] with hydrodynamic coefficients coming from linear potential flow, or the (time-dependent)

<footnote>Université côte d'Azur
Nice, France
Tel.: +33788578073
E-mail: ali.haidar@univ-cotedazur.fr

Université de Montpellier
Montpellier, France
E-mail: fabien.marche@umontpellier.fr

Université de Montpellier
Montpellier, France
E-mail: francois.vilar@umontpellier.fr</footnote>
mild-slope equation, see [3]. Let mention, among others, [36] for the study of offshore structures, [34] for wave-energy converters or [41,32] for floating breakwaters.

In some situations however nonlinear effects may become dominant, like for instance in the evolution of energy converters or breakwaters operating inside the resonance domain, or in nearshore areas. Also, non-linear bathymetry effects (generally neglected within linear-BEM methods) should be accurately modeled in shallow-water as they mainly lead the wave's transformations. With this in mind, recent attempts to account for nonlinear effects generally rely on fully-nonlinear potential flow descriptions together with nonlinear-BEM, see for instance [15,38], or on full Reynolds Averaged Navier-Stokes (RANS) simulations, like for instance [43]. Both strategies are generally considered as extremely computationally demanding. Let also mention the recent and general formulation of [33] for the full water-wave equations with a floating body.

Focusing on particular flow regimes, an alternative to the use of the water-wave equations may be to consider simpler asymptotic (but still nonlinear) flow models. Albeit appearing as an interesting compromise between linear modeling and fully-nonlinear 3D modeling, this strategy has received very few attention in the literature. Specifically focusing on the fully-nonlinear and shallow-water regime, a depth-integrated flow model may be used to describe the free-surface flow away from the floating object, like for instance in [31] where floating breakwater are modeled through the use of Boussinesq-type (BT) equations and a Finite-Difference (FD) scheme. Such shallow-water asymptotics may also model the fluid under the floating object: in [26,27,42] a BT model is applied to describe near-ship flows and a BT model is used to compute the heave motion of floating structures with straight-sided vertical boundaries in [4]. Let us also mention [8,14] among the related works involving depth-integrated models. Among shallow-water asymptotics, the Nonlinear Shallow-Water (NSW) equations [11] are certainly one of the most widely used mathematical model for the simulation of long waves hydrodynamic. Considering its hydrostatic and hyperbolic nature, it provides an accurate description of steep-fronted flows, see for instance [35]. Yet, up to our knowledge, very few studies investigate the description of a floating object with NSW equations. We may mention the recent studies [16,17] for congested shallow-water flows with a pseudo-compressible relaxation/projection scheme, or again [4] in which the NSW equations are finally used in the vicinity of the object instead of BT equations. Very recently in [22], Initial Boundary Value Problems (IBVP) for the NSW equations with a partly immersed and floating object are introduced and extensively analyzed, in the simplified horizontal $d = 1$ case. We ground the present study on this last work, providing a firm mathematical basis for the modeling and simulation of floating objects in shallow-water flows. Note also that such a reduction of the wave-structure interaction problem to free-boundary and transmission problems has been extended even more recently for a BT model in [5] in the case of a stationary surface obstacle, and [2] in the case of a floating object with heave motion and vertical sidewalls.

The design of accurate and robust discrete algorithms for such free-boundary problems is also an intricate issue, especially concerning the description of the contact-points dynamics. Among the simplest approaches, one may either neglect the wetted-surface variations, as initially proposed in [28,29], or consider exclusively pure heaving motions for objects with straight-sided vertical walls, like in [31,4]. We can also refer to [23,24,25] in the specific context of floating breakwaters. Going further in the description, a Mixed Lagrangian-Eulerian (MLE) strategy is developed in [30], for the study of wave-induced motions with a floating body, with additional difficulties related to the computation of the time-derivative of the velocity potential in the elliptic equation for the pressure. In practice, the flow interactions with the time-evolving boundaries may result in unsteady phenomena at various scales, and it may be convenient for the discrete formulation to handle moving and deforming spatial domains, in order to avoid interface-tracking methods while still maintaining the required accuracy, conservation and robustness properties. Hence, the Arbitrary Lagrangian-Eulerian (ALE) description appears as an attractive solution and we refer to [13] for a general review, and also to the references provided in [20].

Relying on an ALE description of the contact-points dynamics, we choose to approximate the main flow variables with discontinuous polynomials by means of discontinuous Galerkin (DG) method, see [9] for a general review. The assets of DG methods are numerous, and briefly recalled in [20]. Nonetheless, one major acknowledged drawback of DG methods is their lack of nonlinear stability. In this work, such a limitation is alleviated by using the recent a posteriori Local Subcell Correction (LSC) method introduced in [40] for general hyperbolic equations, in [19] for the NSW equations, and extended very recently to the ALE framework in [20].
In what follows, gathering all these ingredients, we numerically investigate the nonlinear wave-body interactions in shallow-water, in order to adaptively and accurately simulate the time-evolution of a partly immersed object floating in a nonlinear shallow-water flow. To this end, we introduce a novel discrete formulation based on a robust DG-ALE method and directly modeled from the recent IBVP described and analyzed in [22], which is the first well-posed model that fully describes the general floating object problem in the particular case $d = 1$. An explicit mapping between the initial configuration and the current one is introduced, and the NSW equations are recast in the reference domain before being approximated by high-order broken-polynomials. Additionally, the discrete formulation associated with the dynamics of the floating object relies on the introduction of the related added-mass effect, appearing as an added-mass-inertia matrix and acting as an additional stabilization mechanism. We emphasize that any object with a continuously differentiable underside's parameterization may be considered, and this object can move in any of the three degrees of freedom allowed in the $d = 1$ surface dimension case (heaving, surging and pitching). The remaining of this paper is structured as follows: the next section is devoted to the description of the governing model, based on the NSW equations with a varying bathymetry for the fluid description, both for the exterior and interior sub-domains, and the corresponding IBVP is explicitly stated, for the sake of completeness. The next §3 is devoted to the description and study of a discrete formulation associated with this model. We define a general discrete approximation relying on local polynomials of arbitrary order and several interesting properties are highlighted. In particular, assuming that the initial acceleration of the object vanishes, we emphasize that the well-balancing properties is ensured. In §4, we show several numerical assessments of the resulting discrete algorithm, involving the preservation of and convergence towards motionless steady-states, nonlinear interactions between surface-waves and the floating object or the occurrence of dry areas and propagation of wet-dry front.

## 2 The model

In this section, we recall the equations introduced in [22], that will be numerically investigated in the remainder of the paper.

### 2.1 Shallow-water flow with a floating object

We consider a floating non-deformable object, denoted by $O_{bj}$, of mass $m_o$, inertia coefficient $i_o$, center of mass $M_G$ and diameter $d_o$, which is partly immersed in an inviscid, incompressible and irrotational shallow-water flow, with a reference water-depth $H_0$, under the following assumptions:

1. the object is not entirely immersed, and no wave overtopping can occur,

2. there are only two contact-points where the water, the air, and the object meet, see Fig. 1,

3. $d_o$ is much smaller than the water-depth at rest: $d_o \ll H_0$.

We emphasize that the object should be partly immersed, as there are two contact-points where the water, the air, and the object meet. For the sake of simplicity, we assume that the underside's profile is symmetric with respect to the vertical axis passing through $M_G$. This additional assumption is not mandatory but helps to simplify initial computations related to inertia, see D. For any given time value $t \ge 0$, the horizontal spatial coordinate of these contact-points are denoted by $\chi_-(t)$ and $\chi_+(t)$, with $\chi_-(t) < \chi_+(t)$. The horizontal line is split into two time-dependent sub-domains, namely the interior sub-domain denoted by $\mathcal{I}(t)$ and the exterior sub-domain $\mathcal{E}(t)$, $\mathcal{E}(t)$ and $\mathcal{I}(t)$ being the projections on the horizontal line of the areas where the water surface get in touch with the floating structure and the air:

$$ \mathcal{I}(t) := ]\chi_-(t), \chi_+(t)[, \quad \mathcal{E}(t) := \mathcal{E}^-(t) \cup \mathcal{E}^+(t), \quad \mathcal{E}^-(t) := ]-\infty, \chi_-(t)[, \quad \mathcal{E}^+(t) := ]\chi_+(t), +\infty[. \quad (1) $$

and we conveniently gather the contact-points into the set $\partial\mathcal{I}(t) := \{\chi_-(t), \chi_+(t)\}$. The bathymetry variations are parameterized by a regular function denoted by $b: \mathbb{R} \to \mathbb{R}$, $H^i$ and $u^i$ respectively denote the water-height and the water averaged horizontal velocity in $\mathcal{I}(t)$, $H^e$ and $u^e$ the water-height and the velocity in $\mathcal{E}(t)$ and we set $\eta^i := H^i + b$, $\eta^e := H^e + b$, $q^e := H^e u^e$, $q^i := H^i u^i$ the free-surface elevations and the vertically-averaged horizontal discharge respectively in $\mathcal{E}(t)$ and $\mathcal{I}(t)$. We also assume the pressure field to be hydrostatic:

$$ p(x, z, t) := \begin{cases} \text{patm} - \rho g(z - \eta^e(x, t)) & \text{in } \mathcal{E}(t), \\ \underline{p}^i(x, t) - \rho g(z - \eta^i(x, t)) & \text{in } \mathcal{I}(t), \end{cases} \quad (2) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1: shallow-water interacting with a floating object.</figcaption></figure>

where $\rho$ is the density of the water, $p_{atm}$ the atmospheric pressure (at the fluid free-surface) and $\underline{p}^i(x,t)$ is the inner pressure that applies on the underside of the floating object. We consider the following flow model:

$$ \partial_t \mathbf{v}^e + \partial_x \mathbf{F}(\mathbf{v}^e, b) = \mathbf{B}(\mathbf{v}^e, b') \quad \text{in } \mathcal{E}(t) = ]-\infty, \chi_-(t)] \cup (\chi_+(t), +\infty[, \quad (3a) $$

$$ \partial_t \mathbf{v}^i + \partial_x \mathbf{F}(\mathbf{v}^i, b) = \mathbf{B}(\mathbf{v}^i, b') + \mathbf{P}(\mathbf{v}^i, \partial_x \underline{\mathbf{p}}^i) \quad \text{in } \mathcal{I}(t) = ]\chi_-(t), \chi_+(t)], \quad (3b) $$

$$ \mathbf{v}^e = \mathbf{v}^i \quad \text{on} \quad \partial \mathcal{I}(t), \quad (3c) $$

where

$$ \mathbf{v}^e : \mathcal{E}(t) \times \mathbb{R}^+ \ni (x, t) \mapsto \mathbf{v}^e(x, t) \in \Theta := \{(\eta, q) \in \mathbb{R}^2, H = \eta - b \ge 0\}, \quad (4) $$

$$ \mathbf{v}^i : \mathcal{I}(t) \times \mathbb{R}^+ \ni (x, t) \mapsto \mathbf{v}^i(x, t) \in \Theta, \quad (5) $$

respectively gather the flow's main variables in $\mathcal{E}(t)$ and $\mathcal{I}(t)$, $\mathbf{F}: \Theta \times \mathbb{R} \to \mathbb{R}^2$ is the (nonlinear) flux function, $\mathbf{B}: \Theta \times \mathbb{R} \to \mathbb{R}^2$ is a bathymetry source term and $\mathbf{P}(\mathbf{v}^i, \underline{\mathbf{p}}^i): \Theta \times \mathbb{R} \to \mathbb{R}^2$ is a pressure source term, respectively defined as follows:

$$ \mathbf{F}(\mathbf{v}, b) := \left( \begin{array}{c} q \\ uq + \frac{1}{2}g\eta(\eta - 2b) \end{array} \right), \quad \mathbf{B}(\mathbf{v}, b') := \left( \begin{array}{c} 0 \\ -g\eta b' \end{array} \right), \quad \mathbf{P}(\mathbf{v}^i, \partial_x \underline{\mathbf{p}}^i) := \left( \begin{array}{c} 0 \\ -\frac{1}{\rho} H^i \partial_x \underline{\mathbf{p}}^i \end{array} \right). \quad (6) $$

For further use, let introduce the unit vectors $\mathbf{e}_x := (1,0)^T$ and $\mathbf{e}_z := (0,1)^T$ in the plane $(Oxz)$. Let also introduce the following operators, respectively extracting an average and an oscillating part of any regular enough scalar function $v(\cdot, t)$ defined on $\mathcal{I}(t)$, as follows:

$$ \langle\langle v \rangle\rangle_{\mathcal{I}(t)} := \overline{H^i}^{-1} \int_{\mathcal{I}(t)} \frac{v}{H^i} dx, \quad v_{\mathcal{I}(t)}^* := v - \langle\langle v \rangle\rangle_{\mathcal{I}(t)}, \quad (7) $$

with the notation $\overline{H^i}^{-1} = (\int_{\mathcal{I}(t)} \frac{1}{H^i} dx)^{-1}$ and the subscript $\mathcal{I}(t)$ may be forgotten when no confusion is possible.

## 2.2 Object's position and motion

We consider the case of a moving object and face the following alternative: (i) the motion of the object is prescribed and therefore not influenced by the surface waves, (ii) the motion of the object is free and submitted to the influence of the wave's motion while having a reciprocal influence on the flow. In both cases, for any time value, the position of $\mathcal{O}_{bj}$ is completely specified through the knowledge of the spatial coordinates $\mathbf{x}_G(t) = (x_G(t), z_G(t))$ of $\mathcal{M}_G$ ($x_G, z_G$ are respectively the horizontal and vertical coordinates), together with the (signed) value of the rotation (pitch) angle $\theta(t)$ with the vertical direction, see Fig. 2. In a similar way, the motion of $\mathcal{O}_{bj}$ may be entirely defined through the knowledge of the velocity $\mathbf{v}_G(t) = (u_G(t), w_G(t)) = \mathbf{x}_G'(t)$ and the angular velocity $\omega(t) := -\theta'(t)$ (so that $\theta$ is oriented according to the
standard trigonometric convention in the plane (Oxz)) of M<sub>G</sub>. For the sake of convenience, let introduce the position, velocity and acceleration vectors X<sub>G</sub>, v<sub>G</sub> and A<sub>G</sub>, defined as follows:

$$
\mathcal{X}_G := \begin{pmatrix} x_G \\ z_G \\ -\theta \end{pmatrix}, \quad \mathcal{v}_G := \begin{pmatrix} u_G \\ w_G \\ \omega \end{pmatrix} = \frac{d}{dt} \mathcal{X}_G \quad \text{and} \quad \mathcal{A}_G := \frac{d}{dt} \mathcal{v}_G.
$$

<figure><img src="image_3.png" /><figcaption>Fig. 2: Moving object in water</figcaption></figure>

In what follows, either $t \mapsto (\mathcal{X}_G(t), \vartheta_G(t))$ belongs to the provided data, or it has to be explicitly computed
from external forces and torque. We further observe that:

- at $t=0$, the initial location of $\mathcal{M}_G$ is denoted by $(X_G, Z_G) := (x_G(0), z_G(0))$ and the initial pitch angle $\theta(0)$ is arbitrarily set to zero (so that the inertia momentum vanishes for objects with symmetric profiles, see D). We assume that the underside of the floating object is initially parameterized by a smooth function $\eta_{\text{lid}}$ defined on an open interval $\mathcal{I}_{\text{lid}} \subset \mathbb{R}$, with $\eta_{\text{lid}} \in \mathcal{C}^1(\mathcal{I}_{\text{lid}}) \cap W^{s,\infty}(\mathcal{I}_{\text{lid}})$, $s \ge 1$ ($W^{s,\infty}(\mathcal{I})$ being the Sobolev space of functions which are uniformly bounded on $\mathcal{I}$, together with their weak derivatives up to order $s$). We observe that:

$$
\eta^i(\cdot, 0) = \eta_{\text{lid}} \quad \text{on} \quad \mathcal{I}_0 = ]X_0^-, X_0^+[ \subset \mathcal{I}_{\text{lid}},
$$

and for any material point located on the underside of the object, identified by its coordinates (X, Z),
we have Z = η<sub>lid</sub>(X) = η<sup>i</sup>(X, 0). For further use, let define a normal vector on the underside of the
object:

$$
\mathbf{n}_{\text{lid}}(x) := \begin{pmatrix} -\eta'_{\text{lid}}(x) \\ 1 \end{pmatrix}.
$$

- at any time value $t > 0$, we denote by $\mathbf{x} = (x, z)$, with $z = \eta^i(x, t)$, the coordinates of an arbitrary point belonging to the object's underside, and we note by $\mathbf{r}_G := \mathbf{x} - \mathbf{x}_G$ the translated coordinate vector of this point with respect to $\mathcal{M}_G$, and by $\mathbf{n}^i$ a normal vector on the underside of the object:

$$
\mathbf{n}^i(x, t) := \begin{pmatrix} -\partial_x \eta^i(x, t) \\ 1 \end{pmatrix} = -\partial_x \mathbf{r}_G(x, t)^\perp.
$$

When a free-motion is allowed, its response to external force and torque is ruled by the Newton’s second
law for the conservation of linear and angular momentum, which is formulated as follows:

$$
M_0 \mathcal{A}_G = - \begin{pmatrix} m_o g & e_z \\ 0 & 0 \end{pmatrix} - \int_{\mathcal{I}(t)} \left( \underline{p}^i - p_{\text{atm}} \right) \begin{pmatrix} -\mathbf{n}^i \\ \mathbf{r}_G^\perp \cdot \mathbf{n}^i \end{pmatrix}, \quad (8)
$$

where the mass-inertia matrix is defined as:

$$
M_0 := \begin{pmatrix} m_o & \mathrm{Id}_{2 \times 2} & 0 \\ 0 & 0 & i_o \end{pmatrix}.
$$
For further use in (15), we also introduce the vector $\mathcal{T}_G$ defined as follows:

$$ \mathcal{T}_G(x, t) := \begin{pmatrix} -\mathbf{r}_G^{\perp}(x, t) \\ \frac{1}{2} |\mathbf{r}_G(x, t)|^2 \end{pmatrix}, $$

such that the following identities hold:

$$ \partial_x \mathcal{T}_G = \begin{pmatrix} -\mathbf{n}^i \\ \mathbf{r}_G^{\perp} \cdot \mathbf{n}^i \end{pmatrix}, \quad \partial_t \mathcal{T}_G = \mathcal{M}_G \vartheta_G, \qquad (9) $$

with

$$ \mathcal{M}_G := \begin{pmatrix} \mathbf{e}_x \cdot \mathbf{n}_{\text{lid}} & 0 & -\mathbf{r}_G^{\perp} \cdot \mathbf{n}_{\text{lid}} \\ 1 & 0 & 0 \\ -\mathbf{r}_G^{\perp} \cdot \mathbf{n}_{\text{lid}} & 0 & -(\mathbf{e}_z \cdot \mathbf{r}_G)(\mathbf{r}_G^{\perp} \cdot \mathbf{n}_{\text{lid}}) \end{pmatrix}. $$

## 2.3 Interior-flow description

In this sub-section, we reformulate the flow equations in the interior domain. For any point $\mathbf{x}$ of the object's underside, the corresponding initial coordinates $(X, Z) \in \mathcal{I}_0 \times \mathbb{R}$ can be traced back through the following identity:

$$ \mathbf{r}_G(x, t) = \begin{pmatrix} \cos(\theta(t)) & -\sin(\theta(t)) \\ \sin(\theta(t)) & \cos(\theta(t)) \end{pmatrix} \mathbf{r}_G(X, 0), \qquad (10) $$

and as a consequence, we have for any given $t > 0$ and $x \in \mathcal{I}(t)$:

$$ \eta^i(x, t) = z_G(t) + \sin(\theta(t))(X - X_G) + \cos(\theta(t))(\eta_{\text{lid}}(X) - Z_G) =: \tilde{F}(X, t, z_G(t), \theta; X_G, Z_G, \eta_{\text{lid}}), \quad (11) $$

where $X$ satisfies the following nonlinear algebraic equation:

$$ \frac{x - x_G(t) + \sin(\theta(t))(\eta_{\text{lid}}(X) - Z_G)}{\cos(\theta(t))} + X_G - X = 0. \qquad (12) $$

The ~ notation, related to the fact that $\tilde{F}$ depends on $X$ in the initial object's configuration, is explicitly defined in the next section with the introduction of the ALE description. Under the additional assumptions that: (i) $\mathcal{M}_G$ remains close to its initial location, (ii) the pitch angle is small enough, in the following sense:

$$ \forall t \in [0, T_{\text{max}}], |\theta(t)| \le \theta_{\text{max}}, \text{ with } \theta_{\text{max}} \in (0, \pi/2) \text{ such that } \|\eta_{\text{lid}}\|_{\infty} \tan(\theta_{\text{max}}) < 1, $$

then it is possible to show that: (i) there is a unique $X \in \mathcal{I}_0$ satisfying (12), (ii) the discharge can be expressed as:

$$ q^i(x, t) = \vartheta_G(t) \cdot \mathcal{T}_G(x, t) + q^i(t), \qquad (13) $$

where $q^i$ is the solution of the following IVP:

$$ \frac{d}{dt} q^i = - \left( \langle \! \langle f_1 \rangle \! \rangle_{\mathcal{I}(t)} + \langle \! \langle f_2 \rangle \! \rangle_{\mathcal{I}(t)} + \langle \! \langle f_3 \rangle \! \rangle_{\mathcal{I}(t)} \right), \qquad (14a) $$

$$ q^i(0) := q^i_0, \qquad (14b) $$

with the following right-hand sides:

$$ f_1 := \partial_x(u^i q^i) + gH^i \partial_x \eta^i, \quad f_2 := \mathcal{A}_G \cdot \mathcal{T}_G \quad \text{and} \quad f_3 := \vartheta_G \cdot \partial_t \mathcal{T}_G, \qquad (15) $$

where $f_3$ is obtained using (9). The reader is referred to [22] for the detailed derivation of these relations.

**Remark 1** It is known that, for immersed objects, part of the force and torque applied on the object by the surrounding fluid acts as if the mass-inertia matrix in Newton's laws was modified through the addition of a positive matrix, which is the so-called *added-mass* effect, see for instance [6]. Hence, (8) may be reformulated in order to exhibit the corresponding added-mass, as follows:

$$ (\mathcal{M}_0 + \mathcal{M}_a (H^i, \mathcal{T}_G)) \frac{d}{dt} \vartheta_G = \begin{pmatrix} -m_0 g \mathbf{e}_z \\ 0 \end{pmatrix} - \rho \int_{\mathcal{I}(t)} (f_1^* + f_3^*) \frac{\mathcal{T}_G^*}{H^i}, \qquad (16) $$

where the *added-mass-inertia matrix* $\mathcal{M}_a$ is defined as:

$$ \mathcal{M}_a (H^i, \mathcal{T}_G) := \int_{\mathcal{I}(t)} \frac{\mathcal{T}_G^* \otimes \mathcal{T}_G^*}{H^i}. \qquad (17) $$

This reformulation is detailed in C.
Remark 2 We observe that deriving the second equation of the geometric relation (10) with respect to $t$ and $x$ leads to the following identity for the time derivative of the interior free-surface:

$$ \partial_t \eta^i(\cdot, t) = (\mathbf{v}_G(t) - \mathbf{r}_G(\cdot, t)^\perp \omega(t)) \cdot \mathbf{n}^i(\cdot, t) = -\partial_x (\vartheta_G(t) \cdot \mathcal{T}_G(x, t)). \quad (18) $$

Remark 3 The assumption that $\mathcal{M}_G$ remains close to its initial location helps to ensure that some singular behaviors, that may break-down the existence conditions, do not occur. In particular, the motion $t \mapsto (\mathcal{X}_G(t), \vartheta_G(t))$ should ensure that the object is never entirely immersed, and never touches the bottom boundary. As already mentioned, we typically require that the object's diameter $d_o$ is smaller than the water-depth at rest (or $d_o \ll H_b$ with $H_b := H_0 - b$ when the bathymetry is not flat), in order to ensure that, for all time value, $H^i(\cdot, t) > 0$.

## 2.4 IBVP in the case of a free-motion

Accounting for the description and notations introduced above, the coupled problem (3) may be particularized as follows: find $(\mathbf{v}^e, \mathbf{v}^i, \chi_-, \chi_+, \mathcal{X}_G)$ such that:

$$ \partial_t \mathbf{v}^e + \partial_x \mathbf{F}(\mathbf{v}^e, b) = \mathbf{B}(\mathbf{v}^e, b') \quad \text{in } \mathcal{E}(t) = ]-\infty, \chi_-(t)] \cup (\chi_+(t), +\infty[, \quad (19a) $$

$$ \left. \begin{array}{l} \eta^i(x, t) = \tilde{F}(X, t, z_G, X_G, Z_G, \theta, \eta_{\text{lid}}) \quad \text{where } X \text{ solves (12),} \\ q^i(x, t) = \vartheta_G(t) \cdot \mathcal{T}_G(x, t) + \underline{q}^i(t), \\ \frac{d}{dt} \underline{q}^i = - \langle \mathbf{f}_1 \rangle_{\mathcal{I}(t)} - \langle \mathbf{f}_2 \rangle_{\mathcal{I}(t)} - \langle \mathbf{f}_3 \rangle_{\mathcal{I}(t)}, \end{array} \right\} \text{in } \mathcal{I}(t) = ]\chi_-(t), \chi_+(t)[, \quad (19b) $$

$$ \mathbf{v}^e = \mathbf{v}^i \quad \text{on } \partial\mathcal{I}(t), \quad (19c) $$

$$ \left. \begin{array}{l} \frac{d}{dt} \mathcal{X}_G = \vartheta_G, \\ (\mathbb{M}_0 + \mathbb{M}_a (H^i, \mathcal{T}_G)) \frac{d}{dt} \vartheta_G = \begin{pmatrix} -m_0 g e_z \\ 0 \end{pmatrix} - \rho \int_{\mathcal{I}(t)} (f_1^* + f_3^*) \frac{\mathcal{T}_G^*}{H^i}. \end{array} \right\} \quad (19d) $$

Supplementing (19) with some initial-data, for $s \ge 2$:

$$ \mathbf{v}_{|t=0}^e := \mathbf{v}_0^e \in (H^s(\mathcal{E}_0))^2, \quad (20a) $$

$$ (\chi_-, \chi_+)|_{t=0} := (X_0^-, X_0^+), \quad (20b) $$

$$ \mathbf{v}_{|t=0}^i := (\eta_{\text{lid}}, \underline{q}_0^i) \in (C^1(\mathcal{I}_{\text{lid}}) \cap W^{s,\infty}(\mathcal{I}_{\text{lid}})) \times \mathbb{R}, \quad (20c) $$

$$ \mathcal{X}_G(0) := (X_G, Z_G, 0), \quad (20d) $$

$$ \vartheta_G(0) := (u_G^0, w_G^0, \omega_0), \quad (20e) $$

a local well-posedness result is proved in [22]. Specifically under further assumptions on the data, which can be summarized as: (i) there is no dry-state in the vicinity of the floating structure, (ii) the flow is initially sub-critical at $\chi_\pm$, (iii) the first-order spatial derivative of the free-surface is initially discontinuous at contact points:

$$ (\eta_0^e - \eta_0^i)' \neq 0 \quad \text{on} \quad X_0^\pm, \quad (21) $$

then there exists a maximum time $T_{\text{max}} \le T$ and a unique solution of (19)-(20) such that $\mathbf{v}^e \circ \chi \in C^0([0, T_{\text{max}}]; H^s(\mathcal{E}_0)) \cap C^1([0, T_{\text{max}}]; H^{s-1}(\mathcal{E}_0))$, $\underline{q}^i \in H^{s+1}(0, T_{\text{max}})$, $(\chi_-, \chi_+) \in (H^s(0, T_{\text{max}}))^2$ and $\mathcal{X}_G \in (H^{s+2}(0, T_{\text{max}}))^3$.

Remark 4 A well-posedness result is also established in [22] for the case of an object with a prescribed motion. The corresponding model may be straightforwardly deduced from (19)-(20) by: (i) removing the Newton's equation (19d), (ii) removing the initial data on $\vartheta_G$ and $\mathcal{X}_G$ in (20), (iii) supplementing the resulting model with the prescribed evolution law:

$$ \mathcal{X}_G \in (H^{s+2}(0, T))^3, \quad (22a) $$

$$ \mathcal{X}_G(0) := (X_G, Z_G, 0) \in \mathbb{R}^3. \quad (22b) $$
Such prescribed motions are also numerically investigated in §4.

**Remark 5** In the next section, as for the numerical validations of §4, we consider these IBVPs on a bounded computational domain of the form

$$ \Omega_t := \mathcal{E}^{-}(t) \cup \mathcal{I}(t) \cup \mathcal{E}^{+}(t) = ]x_{\text{left}}(t), x_{-}(t)[ \cup ]\mathcal{I}(t) \cup ]x_{+}(t), x_{\text{right}}(t)[, $$

so that the domain's exterior boundary is defined as $\partial\Omega_t := \{x_{\text{left}}(t), x_{\text{right}}(t)\}$.

**Remark 6** As shown in [33], one has conservation of the total energy

$$ \frac{d}{dt} \varepsilon_{\text{tot}} = 0 \quad \text{with} \quad \varepsilon_{\text{tot}} := \varepsilon_{\text{obj}} + \varepsilon_{\text{sw}}, \qquad (23) $$

with, for all $t \ge 0$:

$$ \varepsilon_{\text{sw}}(t) := \int_{\Omega(t)} \frac{H(\cdot, t)u^2(\cdot, t)}{2} dx + g \int_{\Omega(t)} H(\cdot, t)(b(\cdot) + \frac{H(\cdot, t)}{2}) dx, \qquad (24) $$

and

$$ \varepsilon_{\text{obj}}(t) := m_o g Z_G + \frac{1}{2} M_0 \vartheta_G \cdot \vartheta_G. \qquad (25) $$

# 3 Discrete formulations

In this section, we design a discrete formulation associated with the previous model. The reader is also referred to our previous works [19, 40] for additional details.

## 3.1 Discrete setting for DG-ALE on mesh elements and FV-ALE on sub-cells

*Computational domain, sub-domains and mesh*

For any time value $t \in [0, T_{\text{max}}]$, we consider the partition $\mathcal{P}_{\Omega}(t) := \{\mathcal{E}^{-}(t), \mathcal{I}(t), \mathcal{E}^{+}(t)\}$ of $\Omega_t$ into disjoint sub-domains, defined from the knowledge of the location of the contact points $x_{-}(t) < x_{+}(t)$. We introduce a conforming partition $\mathcal{F}_h(t) := \{c_i(t)\}_{1 \le i \le n_{\text{el}}}$ of $\Omega_t$ into $|\mathcal{F}_h(t)| =: n_{\text{el}}$ disjoint segments, such that we have $\overline{\Omega_t} = \bigcup_{c(t) \in \mathcal{F}_h(t)} \overline{c(t)}$. Such a partition is characterized by the (time-dependent) mesh-size

$$ h := \max_{c(t) \in \mathcal{F}_h(t)} h_{c(t)}, $$

where $h_c$ is the length of the element $c$, and where we hide the time dependency of $h$ for the sake of simplicity. We make the following additional assumptions:

1. $n_{\text{el}}$ does not depend on time,

2. $\forall t \in [0, T_{\text{max}}]$, $x_{-}(t) \neq x_{\text{left}}(t)$ and $x_{+}(t) \neq x_{\text{right}}(t)$,

3. $\mathcal{F}_h(t)$ is compatible with $\mathcal{P}_{\Omega}(t)$: each mesh element $c(t) \in \mathcal{F}_h(t)$ is a subset of only one set of the partition $\mathcal{P}_{\Omega}(t)$,

so that we can write:

$$ \mathcal{F}_h(t) = \mathcal{F}_h^e(t) \cup \mathcal{F}_h^i(t), \quad \text{with} \quad \overline{\mathcal{E}(t)} = \bigcup_{c(t) \in \mathcal{F}_h^e(t)} \overline{c(t)} \quad \text{and} \quad \overline{\mathcal{I}(t)} = \bigcup_{c(t) \in \mathcal{F}_h^i(t)} \overline{c(t)}, $$

where $\mathcal{F}_h^e(t)$ and $\mathcal{F}_h^i(t)$ are respective partitions of the sub-domains $\mathcal{E}(t)$ and $\mathcal{I}(t)$, and at any time $t \in [0, T_{\text{max}}]$, the contact points $x_{-}(t), x_{+}(t)$ are uniquely identified with some mesh interfaces. For some specified mesh element $c_i(t) \in \mathcal{F}_h(t)$, we note $c_i(t) := ]x_{i-\frac{1}{2}}(t), x_{i+\frac{1}{2}}(t)[$ (with the convention that $x_{\frac{1}{2}} := x_{\text{left}}$, $x_{n_{\text{el}}+\frac{1}{2}} := x_{\text{right}}$), $x_i(t)$ its barycenter and $\partial c_i(t) := \{x_{i-\frac{1}{2}}(t), x_{i+\frac{1}{2}}(t)\}$ its boundary. Mesh interfaces are collected in the sets $\partial \mathcal{F}_h^e$ and $\partial \mathcal{F}_h^i$, respectively defined as follows:

$$ \partial \mathcal{F}_h^e(t) := \{\partial c(t), c(t) \in \mathcal{F}_h^e(t)\}, \quad \partial \mathcal{F}_h^i(t) := \{\partial c(t), c(t) \in \mathcal{F}_h^i(t)\}, $$

such that we have

$$ \partial \mathcal{I}(t) = \partial \mathcal{F}_h^e(t) \cap \partial \mathcal{F}_h^i(t), \quad \partial \mathcal{F}_h(t) := \partial \mathcal{F}_h^e(t) \cup \partial \mathcal{F}_h^i(t) = \{x_{i+\frac{1}{2}}(t), 0 \le i \le n_{\text{el}}\}. \qquad (26) $$
For any integer $k \ge 0$ and for any $t \in [0, T_{\max}]$, we consider the broken-polynomials space defined on the exterior domain:

$$ \mathbb{P}^k(\mathcal{T}_h^e(t)) := \{v(\cdot, t) \in L^2(\mathcal{E}(t)), \forall c(t) \in \mathcal{T}_h^e(t), v_{|c(t)} \in \mathbb{P}^k(c(t))\}, $$

where $\mathbb{P}^k(c(t))$ denotes the space of polynomials of total degree at most $k$ defined onto $c(t)$, with dim $(\mathbb{P}^k(c(t))) = k + 1$. In what follows, piecewise polynomial functions (and, more generally, any discrete counterpart computed from or acting on piecewise polynomial functions) are denoted with a subscript $h$. Also, for any $c(t) \in \mathcal{T}_h^e(t)$ and $v_h(\cdot, t) \in \mathbb{P}^k(\mathcal{T}_h^e(t))$, we may use the convenient shorthand: $v_c := v_{h|c}$ when no confusion is possible, and we also note $\mathbb{P}^k(\mathcal{T}_h^e(t)) := (\mathbb{P}^k(\mathcal{T}_h^e(t)))^2$.

For any mesh element $c(t) \in \mathcal{T}_h^e(t)$ and any integer $k \ge 0$, we consider a basis for $\mathbb{P}^k(c(t))$ denoted by

$$ \Psi_{c(t)} = \{\psi_m^c(\cdot, t)\}_{m \in [1, k+1]}. $$

We observe that we have:

$$ \forall t \in [0, T_{\max}], \quad \forall c(t) \in \mathcal{T}_h^e(t), \quad \forall m \in [1, k+1], \quad \text{supp}(\psi_m^c(\cdot, t)) \subset \overline{c(t)}. $$

A basis for the global space $\mathbb{P}^k(\mathcal{T}_h^e(t))$ is obtained by gathering the local basis functions:

$$ \Psi_h^e(t) := \bigtimes_{c(t) \in \mathcal{T}_h^e(t)} \Psi_{c(t)} = \{\{\psi_m^c(\cdot, t)\}_{m \in [1, k+1]}\}_{c(t) \in \mathcal{T}_h^e(t)}. $$

Remark 7 Similarly to [20], we choose the set of monomials in the physical space as basis functions, defined as follows:

$$ \forall c_i(t) \in \mathcal{T}_h^e(t), \forall m \in [1, \dots, k+1], \forall x \in c_i(t), \quad \psi_m^{c_i}(x, t) := \left( \frac{x - x_i(t)}{|c_i(t)|} \right)^m. \quad (27) $$

For any given time value, the degrees of freedom are chosen to be the functionals that map a given discrete unknown belonging to $\mathbb{P}^k(\mathcal{T}_h^e(t))$ to the coefficients of its expansion on the chosen basis functions. The degrees of freedom applied to a given function $v_h \in \mathbb{P}^k(\mathcal{T}_h^e)$ return the real numbers

$$ \{\underline{v}_m^c\}_{m \in [1, k+1]}, \quad \text{such that} \quad v_c = \sum_{m=1}^{k+1} \underline{v}_m^c \psi_m^c(\cdot, t), \quad \forall c \in \mathcal{T}_h^e(t). \quad (28) $$

With a little abuse, we refer hereafter to the real numbers (28) as the degrees of freedom associated with $v_h$ and we note $\underline{v}_c \in \mathbb{R}^{k+1}$ the vector that gathers the degrees of freedom associated with $v_c$.

In a similar way, the approximating space $\mathbb{P}^k(\mathcal{T}_h^i(t))$ may be defined for the interior domain, and the global approximation space $\mathbb{P}^k(\mathcal{T}_h^i(t))$ is obtained by gathering the contributions coming from both sub-domains. The product spaces $\mathbb{P}^k(\mathcal{T}_h^i(t))$ and $\mathbb{P}^k(\mathcal{T}_h^i(t))$ are defined accordingly. Reversely, any function $\phi_h(\cdot, t) \in \mathbb{P}^k(\mathcal{T}_h^i(t))$ may be regarded as the gathering of its contributions coming respectively from the exterior and interior domains:

$$ \phi_h(\cdot, t)|_{\mathcal{E}(t)} = \phi_h^e(\cdot, t), \quad \phi_h(\cdot, t)|_{\mathcal{I}(t)} = \phi_h^i(\cdot, t), \quad \phi_h^e(\cdot, t) \in \mathbb{P}^k(\mathcal{T}_h^e(t)), \quad \phi_h^i(\cdot, t) \in \mathbb{P}^k(\mathcal{T}_h^i(t)). $$

### Projection, interpolation, averages and jumps

For $c(t) \in \mathcal{T}_h^e(t)$, we denote by $p_c^k$ the $L^2$-orthogonal projector onto $\mathbb{P}^k(c(t))$ and $p_{\mathcal{T}_h^e}^k$ the $L^2$-orthogonal projector onto $\mathbb{P}^k(\mathcal{T}_h^e(t))$. Similarly, we denote $i_c^k$ the element nodal interpolator into $\mathbb{P}^k(c(t))$, where the nodal distributions in mesh elements are chosen to be the approximate optimal nodes of [7] (leading to better approximation properties than equidistant distributions, and including, for each element, the elements boundaries into the interpolation nodes). The global interpolator into $\mathbb{P}^k(\mathcal{T}_h^e(t))$, denoted by $i_{\mathcal{T}_h^e}^k$, is obtained by gathering the local interpolating polynomials defined on each elements. Similar projector $p_{\mathcal{T}_h^i}^k$ and interpolator $i_{\mathcal{T}_h^i}^k$ may be defined on $\mathcal{I}(t)$, and globally on $\Omega_t$ by gathering the sub-domains contributions.
For any $\phi_h(\cdot, t) \in \mathbb{P}^k(\mathcal{H}_h(t))$ defined on $c_i(t) \cup c_{i+1}(t)$, we introduce the following interface-centered average $\{\cdot\}$ and jump $[\cdot]$ operators defined as follows:

$$ \{\phi_h(\cdot, t)\}_{i+\frac{1}{2}} := \frac{1}{2} \left( \phi_{c_i}(\cdot, t)|_{x_{i+\frac{1}{2}}} + \phi_{c_{i+1}}(\cdot, t)|_{x_{i+\frac{1}{2}}} \right), \quad \{\phi_h(\cdot, t)\}_{i+\frac{1}{2}} := \phi_{c_{i+1}}(\cdot, t)|_{x_{i+\frac{1}{2}}} - \phi_{c_i}(\cdot, t)|_{x_{i+\frac{1}{2}}} $$

and this definition should be supplemented with suitable values for the averages and jumps at exterior boundaries, depending on the chosen type of boundary conditions. For any regular-enough scalar-valued function $v(\cdot, t)$ defined on $c_i(t)$, and extending the convenient notation $v_{c_i}(t)(\cdot) := v(\cdot, t)|_{c_i(t)}$, we also introduce the cell-centered jump value defined as:

$$ \{v(\cdot, t)\}_{\partial c_i(t)} := v_{c_i(t)}|_{x_{i+\frac{1}{2}}} - v_{c_i(t)}|_{x_{i-\frac{1}{2}}}, $$

together with the following shorthands for the exterior scalar-products of functions $v, w \in L^2(\mathcal{H}_h^e)$ and $\mu, \nu \in L^2(\partial \mathcal{H}_h^e)$:

$$ (v, w)_{\mathcal{H}_h^e(t)} := \sum_{c(t) \in \mathcal{H}_h^e(t)} \int_{c(t)} v(x,t)w(x,t)dx, \quad \langle \mu, \nu \rangle_{\partial \mathcal{H}_h^e(t)} := \sum_{c(t) \in \mathcal{H}_h^e(t)} \{\mu\nu\}_{\partial c(t)}, $$

the extension to vector-valued functions being straightforward.

## Discrete derivation and integration

In what follows, we need a consistent and accurate discrete counterpart of the first-order derivative, which may be applied to the broken polynomial functions defined above, while accounting for the domain partition $\mathcal{P}_\Omega(t)$ and the jumps of the functions at interfaces. This may be achieved in the current setting by adapting the liftings and discrete gradient of [1] to the sub-domains partition. Let define the element-by-element first-order derivative of a broken-polynomial belonging to $\mathbb{P}^k(\mathcal{H}_h^e(t))$:

$$ \partial_x^h : \mathbb{P}^k(\mathcal{H}_h^e(t)) \ni \phi_h^e(\cdot, t) \mapsto \partial_x^h \phi_h^e(\cdot, t) \in \mathbb{P}^{k-1}(\mathcal{H}_h^e(t)), $$

such that:

$$ (\partial_x^h \phi_h^e)_{|c(t)} := \partial_x(\phi_{c(t)}^e), \quad \forall c(t) \in \mathcal{H}_h^e(t). $$

Then, for any $\phi_h^e(\cdot, t) \in \mathbb{P}^k(\mathcal{H}_h^e(t))$, we introduce the following global lifting of the jumps on the exterior mesh interfaces $\partial \mathcal{H}_h^e$, defined as follows:

$$ \mathcal{R}_h^k(\{\phi_h^e(\cdot, t)\}) := \sum_{x_{i+\frac{1}{2}} \in \partial \mathcal{H}_h^e} r_{i+\frac{1}{2}}^k (\{\phi_h^e(\cdot, t)\}), $$

where, for all $x_{i+\frac{1}{2}}(t) \in \partial \mathcal{H}_h^e(t)$, the local lifting operator $r_{i+\frac{1}{2}}^k$ applied to the jumps of $\phi_h^e(\cdot, t)$ is defined as the unique solution in $\mathbb{P}^k(\mathcal{H}_h^e(t))$ of the following problem:

$$ (r_{i+\frac{1}{2}}^k(\{\phi_h^e(\cdot, t)\}), \psi_h^e(\cdot, t))_{\mathcal{H}_h^e(t)} = \{\phi_h^e(\cdot, t)\}_{i+\frac{1}{2}} \{\psi_h^e(\cdot, t)\}_{i+\frac{1}{2}}, \quad \forall \psi_h^e(\cdot, t) \in \mathbb{P}^k(\mathcal{H}_h^e(t)). \quad (29) $$

In order to apply the definition (29) to the interfaces corresponding to the contact-points (which are boundaries for $\mathcal{E}(t)$), the definitions of the interface-centered jumps and averages on $\chi_\pm(t)$ have to be provided. Denoting by $\bar{i}$ and $\bar{\iota}$ the respective mesh element labels such that $\chi_-(t) = c_\bar{i}(t) \cap c_{\bar{i}+1}(t)$ and $\chi_+(t) = c_\bar{\iota}(t) \cap c_{\bar{\iota}+1}(t)$, we set:

$$ \{\phi_h^e(\cdot, t)\}_{\chi_-(t)} := \phi_{c_{\bar{i}}}(\cdot, t)|_{\chi_\pm} - \phi_{c_{\bar{i}+1}}^i(\cdot, t)|_{\chi_\pm}, $$

$$ \{\phi_h^e(\cdot, t)\}_{\chi_+(t)} := \phi_{c_{\bar{i}}}(\cdot, t)|_{\chi_\pm} - \phi_{c_{\bar{i}+1}}^i(\cdot, t)|_{\chi_\pm}, $$

$$ \{\phi_h^e(\cdot, t)\}_{\chi_+(t)} := \frac{1}{2} (\phi_{c_{\bar{i}}}(\cdot, t)|_{\chi_\pm} + \phi_{c_{\bar{i}+1}}^i(\cdot, t)|_{\chi_\pm}), $$

$$ \{\phi_h^e(\cdot, t)\}_{\chi_+(t)} := \frac{1}{2} (\phi_{c_{\bar{i}}}(\cdot, t)|_{\chi_\pm} + \phi_{c_{\bar{i}+1}}^i(\cdot, t)|_{\chi_\pm}). $$

Following [12], we define the discrete first-order derivative $\mathcal{G}_h^k : \mathbb{P}^k(\mathcal{H}_h^e(t)) \to \mathbb{P}^k(\mathcal{H}_h^e(t))$ such that, for all $\phi_h^e(\cdot, t) \in \mathbb{P}^k(\mathcal{H}_h^e(t))$,

$$ \mathcal{G}_h^k \phi_h^e(\cdot, t) := \partial_x^h \phi_h^e(\cdot, t) - \mathcal{R}_h^k(\{\phi_h^e(\cdot, t)\}). \quad (30) $$

This operator has better consistency properties than the element-by-element derivative, as it accounts for the jumps of its argument through the second contribution; see [12, Theorem 2.2] for further insight into this
point. In a similar way, a discrete derivative acting on functions of $\mathbb{P}^k(\mathcal{R}_h^i)$ may be defined, and a discrete gradient globally defined on $\mathbb{P}^k(\mathcal{R})$ is obtained by gathering both contributions coming from the exterior and interior domains.

Let also introduce a discrete counterpart for the integration of a regular-enough function $\phi(\cdot, t)$ on $\mathcal{I}(t)$:

$$ \mathcal{G}_{\mathcal{I}(t)}^{h, n_g}[\phi] := \sum_{c(t) \in \mathcal{R}_h^i(t)} \sum_{1 \le r \le n_g} \alpha_r^c(t) \phi(x_r^c(t), t), \quad (31) $$

where $(\alpha_r^c(t))_{1 \le r \le n_g}$ and $(x_r^c(t))_{1 \le r \le n_g}$ respectively refer to some suitable positive Gauss quadrature weights and nodes transported onto the transient mesh element $c(t) \in \mathcal{R}_h^i(t)$, and the degree $n_g$ may be adapted to the polynomial degree of the integrand (or estimated from the regularity of non-polynomial integrands). From (31), we deduce a discrete counterpart of the $\mathcal{I}(t)$-averaging operator (7), as follows:

$$ \langle v \rangle_h := \mathcal{G}_{\mathcal{I}(t)}^{h, n_g} \left[ \frac{1}{H_h^i} \right]^{-1} \mathcal{G}_{\mathcal{I}(t)}^{h, n_g} \left[ \frac{v}{H_h^i} \right], \quad v_h^* := v - \langle v \rangle_h. \quad (32) $$

### FV on sub-cells: sub-partitions, sub-resolution basis and sub-mean values

For any mesh element $c_i(t) \in \mathcal{R}_h^e(t)$, we introduce a sub-partition $\mathcal{R}_i(t)$ into $k+1$ open disjoint sub-cells:

$$ \bar{c}_i(t) = \bigcup_{m=1}^{k+1} \bar{s}_m^{c_i}(t), \quad (33) $$

where the sub-cell $s_m^{c_i}(t) := (\tilde{x}_{m-\frac{1}{2}}^{c_i}(t), \tilde{x}_{m+\frac{1}{2}}^{c_i}(t))$ is of size $|s_m^{c_i}| = |\tilde{x}_{m+\frac{1}{2}}^{c_i} - \tilde{x}_{m-\frac{1}{2}}^{c_i}|$, with the convention $\tilde{x}_{\frac{1}{2}}^{c_i} := x_{i-\frac{1}{2}}$ and $\tilde{x}_{k+\frac{3}{2}}^{c_i} := x_{i+\frac{1}{2}}$, see Fig. 3. When considering a sequence of neighboring mesh cells $c_{i-1}, c_i, c_{i+1}$, the convenient convention $s_0^{c_i} := s_{k+1}^{c_i-1}$ and $s_{k+2}^{c_i} := s_1^{c_i+1}$ may be used. For any regular enough function $v(\cdot, t)$ defined on $s_m^{c_i}$, we introduce the following shorthand for the sub-cell centered jump:

$$ \llbracket v(\cdot, t) \rrbracket_{\partial s_m^{c_i}} := v(\cdot, t)|_{\tilde{x}_{m+\frac{1}{2}}^{c_i}} - v(\cdot, t)|_{\tilde{x}_{m-\frac{1}{2}}^{c_i}}. $$

<figure><img src="image_4.png" /><figcaption>Fig. 3: Partition of a mesh element $c_i$ in $k+1$ sub-cells</figcaption></figure>

Similarly to [40, 20], for $c_i(t) \in \mathcal{R}_h^e(t)$, the sub-cell indicator functions $\{1_{s_m^{c_i}}(\cdot, t), m \in \llbracket 1, k+1 \rrbracket\}$ are defined as follows:

$$ 1_{s_m^{c_i}}(x, t) := \begin{cases} 1 & \text{if } x \in s_m^{c_i}(t) \\ 0 & \text{if } x \notin s_m^{c_i}(t) \end{cases}, \quad \forall j \in \llbracket 1, k+1 \rrbracket, $$

and the sub-resolution basis functions $\{\phi_m^{c_i}(\cdot, t) \in \mathbb{P}^k(c_i(t)), m \in \llbracket 1, k+1 \rrbracket\}$ as follows:

$$ \phi_m^{c_i}(\cdot, t) := \mathbb{P}_{c_i}^k (1_{s_m^{c_i}}(\cdot, t)), \quad \forall m \in \llbracket 1, k+1 \rrbracket. \quad (34) $$

For $v_k \in \mathbb{P}^k(c(t))$, $\{\bar{v}_m^c\}_{m \in \llbracket 1, k+1 \rrbracket}$ are the lowest-order piecewise constant components defined as the mean-values of $v_k$ on the sub-cells belonging to $\mathcal{R}(t)$, called sub-mean values in the following, and which may be conveniently gathered into $\bar{v}_k \in \mathbb{R}^{k+1}$. Whenever a sequence of neighboring mesh elements $c_{i-1}, c_i, c_{i+1}$ and associated neighboring approximations is considered, the following convenient convention may be used: $\bar{v}_0^{c_i} := \bar{v}_{k+1}^{c_i-1}$ and $\bar{v}_{k+2}^{c_i} := \bar{v}_1^{c_i+1}$.
Remark 8 As shown in [40], any polynomial function $v_c \in \mathbb{P}^k(c)$ can be equivalently expressed as a collection of degrees of freedom $\underline{v}_t$ or a collection of sub-means values $\tilde{v}_c$. Defining the set of piecewise constant functions on the sub-grid:

$$ \mathbb{P}^0(\mathcal{F}_i(t)) := \{v(\cdot, t) \in L^2(c_i(t)), v_{|s_m^i}(t) \in \mathbb{P}^0(s_m^i(t)), \forall s_m^i(t) \in \mathcal{F}_i(t)\}, $$

the (one-to-one) projector onto $\mathbb{P}^0(\mathcal{F}_i(t))$ is denoted as follows:

$$ \begin{gathered} \pi_{\mathcal{F}_c}^k : \mathbb{P}^k(c(t)) \to \mathbb{P}^0(\mathcal{F}_c(t)) \tag{35} \\ v_c \mapsto \pi_{\mathcal{F}_c}^k(v_c) := \bar{v}_c. \end{gathered} $$

## 3.2 ALE description

In this section, an ALE description for the problem (19) is introduced. Such a description relies on the definition of a regular (continuously differentiable with respect to time, globally continuous and piecewise continuously differentiable with respect to $X$) mapping which recasts the equations from the initial domain $\Omega_0$ to the current moving one $\Omega_t$:

$$ \Omega_0 \times [0, T_{\max}] \ni (X, t) \mapsto x(X, t) \in \Omega_t, \tag{36} $$

where $X$ refers to the reference coordinate and $x := x(X, t)$ the corresponding physical coordinate. Denoting by $v_g(x, t)$ the frame's velocity at $x := x(X, t)$, we have:

$$ \partial_t x(X, t) = v_g(x(X, t), t). \tag{37} $$

Introducing, for any regular-enough function $v : \Omega_t \times [0, T_{\max}] \to \mathbb{R}$ depending on $(x, t)$, the corresponding function $\tilde{v} : \Omega_0 \times [0, T_{\max}] \to \mathbb{R}$ such that:

$$ \tilde{v}(X, t) := v(x(X, t), t), \tag{38} $$

the central ALE relation between the total time-derivative, the Eulerian time-derivative and the spatial-derivative may be written as follows:

$$ \frac{d}{dt} v(x(X, t), t) := (\partial_t + v_g \partial_x) v(x(X, t), t) := \partial_t \tilde{v}(X, t). \tag{39} $$

### Coordinate transformation and frame's motion

In order to build such a mapping, for any given time value, the velocity of the contact-points may be deduced from the current flow configuration, deriving the free-surface continuity condition (3c) with respect to time, as follows:

$$ \frac{d}{dt} \eta^e = \frac{d}{dt} \eta^i \quad \text{on} \quad \chi_{\pm}, \tag{40} $$

which gives, using the identity $\partial_t \eta^e = -\partial_x q^e$:

$$ v_g|_{\chi_{\pm}} = \left( \frac{\partial_x q^e + \partial_t \eta^i}{\partial_x \eta^e - \partial_x \eta^i} \right) \Big|_{\chi_{\pm}}. \tag{41} $$

Having such contact points velocity at hand, let consider the following smooth diffeomorphism $\chi(\cdot, t) : \mathcal{E}_0 \to \mathcal{E}(t)$, defined as:

$$ \chi(X, t) := \begin{cases} X + \varphi \left( \frac{X - X_0^-}{\epsilon} \right) (\chi_-(t) - X_0^-) & \text{for } X \in \mathcal{E}_0^-, \\ X + \varphi \left( \frac{X - X_0^+}{\epsilon} \right) (\chi_+(t) - X_0^+) & \text{for } X \in \mathcal{E}_0^+, \end{cases} \tag{42} $$

where $\varphi \in C_0^\infty(\mathbb{R})$ is a cut-off function satisfying $\varphi(x) = 1$ for $|x| \le 1$ and $\epsilon := \epsilon_0 \ell$ (the reader is referred to A for the practical definition of $\varphi$, $\epsilon_0$ and Remark 9 for additional considerations regarding the value of $\ell$). Then, for any moving grid's interface $x_{i+\frac{1}{2}}(t) := x(X_{i+\frac{1}{2}}, t)$, we define the corresponding interface's velocity as follows:

$$ v_{g|_{i+\frac{1}{2}}}(t) := \tilde{v}_g(X_{i+\frac{1}{2}}, t), $$
with:

$$
\tilde{v}_g(X_{i+\frac{1}{2}}, t) := \begin{cases}
\partial_t \chi(\cdot, t)|_{X_{i+\frac{1}{2}}} = \begin{cases} \varphi\left(\frac{X_{i+\frac{1}{2}} - X_0^-}{\varepsilon}\right) v_{g|\chi_-} & \text{if } X_{i+\frac{1}{2}} \in \mathcal{E}_0^-, \\ \varphi\left(\frac{X_{i+\frac{1}{2}} - X_0^+}{\varepsilon}\right) v_{g|\chi_+} & \text{if } X_{i+\frac{1}{2}} \in \mathcal{E}_0^+, \end{cases} & (43) \\
\frac{(X_0^+ - X_{i+\frac{1}{2}})}{|\mathcal{I}_0|} v_{g|\chi_-} + \frac{(X_{i+\frac{1}{2}} - X_0^-)}{|\mathcal{I}_0|} v_{g|\chi_+} & \text{if } X_{i+\frac{1}{2}} \in \mathcal{I}_0.
\end{cases}
$$

The updated locations of the mesh interfaces of $\partial \mathcal{F}_h(t)$ are deduced from the frame velocity as the solutions of the following family of IVPs:

$$
\begin{align}
\partial_t x(X_{i+\frac{1}{2}}, t) &= v_{g|_{i+\frac{1}{2}}}(t), \tag{44} \\
x(X_{i+\frac{1}{2}}, 0) &:= X_{i+\frac{1}{2}}. \nonumber
\end{align}
$$

As a result, one has available the following sets of discrete grid's interfaces velocities $(v_{g|_{i+\frac{1}{2}}}(t))_{0\le i \le n_{el}}$ and locations $(x_{i+\frac{1}{2}}(t))_{0\le i \le n_{el}}$.

Remark 9 (42) allows to dispatch the mesh elements in the moving exterior sub-domain $\mathcal{E}(t)$, preventing
from elements collapsing and distorting to occur. It also enables to properly deal with the possible occurrence
of dry areas, provided that such areas are initially far enough from the object, to prevent the water-height
from vanishing at contact points. Indeed, assuming that the distance between $\chi_+(t)$ and the nearest mesh
interface where the water-height vanishes is greater than $\ell$, function (42) then ensures that this mesh
interface's location does not vary over time.

Mapping and geometric parameters

A suitable definition for the mapping (36) can now be provided and we consider

$$
\Omega_0 \times [0, T_{\max}] \ni (X, t) \mapsto x(X, t) \in \Omega_t,
$$

such that, for any $c_i(0) := ]X_{i-\frac{1}{2}}, X_{i+\frac{1}{2}}[ \in \mathcal{T}_h(0)$, $X \in c_i(0)$ and $t \in [0, T_{\max}]$:

$$
x|_{c_i(0)}(X,t) := \frac{(X_{i+\frac{1}{2}} - X)}{|c_i(0)|} x_{i-\frac{1}{2}}(t) + \frac{(X - X_{i-\frac{1}{2}})}{|c_i(0)|} x_{i+\frac{1}{2}}(t). \quad (45)
$$

With this mapping at hand, the velocity of any physical points belonging to the moving frame can be
deduced as follows:

**Proposition 1** The frame's velocity is such that, for all $t \in [0, T_{\max}]$ and all mesh element $c_i(t) =$
$(x_{i-\frac{1}{2}}(t), x_{i+\frac{1}{2}}(t)) \in \mathcal{F}_h(t)$, we have:

$$
\forall x \in c_i(t), v_{g|c_i(t)}(x,t) = \frac{(x_{i+\frac{1}{2}}(t) - x)}{|c_i(t)|} v_{g|_{i-\frac{1}{2}}}(t) + \frac{(x - x_{i-\frac{1}{2}}(t))}{|c_i(t)|} v_{g|_{i+\frac{1}{2}}}(t). \quad (46)
$$

*Proof* Deriving (45) with respect to time gives:

$$
\tilde{v}_{g|\omega_i(0)}(X,t) = \frac{(X_{i+\frac{1}{2}} - X)}{|c_i(0)|} v_{g|_{i-\frac{1}{2}}}(t) + \frac{(X - X_{i-\frac{1}{2}})}{|c_i(0)|} v_{g|_{i+\frac{1}{2}}}(t).
$$

The deformation gradient associated with the grid's motion is obtained as the Jacobian of this mapping. In
particular, the following identities are satisfied:

$$
\partial_X x(X, t)|_{c_i(t)} := \mathcal{J}_{c_i(t)} = \frac{|c_i(t)|}{|c_i(0)|},
$$

$$
\partial_X^k x(X, t)|_{c_i(t)} = 0, \forall k \ge 2,
$$
so that the mapping is invertible and orientation-preserving. Also, for any $(X_a, X_b) \in (c_i(0))^2$, we have:

$$x(X_b, t) = x(X_a, t) + (X_b - X_a)\mathcal{J}_{c_i(t)},$$

and in particular, we deduce (46).

**Remark 10** From (37), we observe that $\mathcal{J} := |\mathcal{J}|$ satisfies the following fundamental relation, generally referred to as Geometric Conservation Law (GCL):

$$\partial_t \mathcal{J}(X, t) = \mathcal{J} \partial_x v_g(x(X, t), t). \quad (48)$$

**Remark 11** In [20], the reader should further find how this particular discretization of the frame motion satisfies the GCL, as well as the Discrete GCL (DGCL) for the fully discrete scheme. It has been also proved that the chosen basis function $\psi_p^e$ as well as the sub-resolution ones $\phi_p^e$ follow the trajectories, as

$$\forall c(t) \in \mathcal{F}_h^e(t), \forall p \in [1, \dots, k+1], \quad \frac{d}{dt} \psi_p^e(x(X, t), t) = \frac{d}{dt} \phi_p^e(x(X, t), t) = 0. \quad (49)$$

This latter property is required to establish the following local weak-formulation satisfied by the flow variables in $\mathcal{E}(t)$:

**Proposition 2** Let $\mathbf{v}_h^e \in \mathbf{P}^k(\mathcal{F}_h^e(t))$ satisfying (3a). Then, for all $\varphi_h(\cdot, t) \in \mathbb{P}^k(\mathcal{F}_h^e(t))$ such that $\frac{d}{dt}\varphi_h(x(X, t), t) = 0$, and all $c(t) \in \mathcal{F}_h^e(t)$, the following identity holds:

$$\frac{d}{dt} \int_{c(t)} \mathbf{v}_h^e \varphi_h dx - \int_{c(t)} \mathbf{G}(\mathbf{v}_h^e, b_h, v_g) \partial_x \varphi_h dx + \left[ \varphi_h \mathbf{G}(\mathbf{v}_h^e, b_h, v_g) \right]_{\partial c(t)} = \int_{c(t)} \varphi_h \mathbf{B}(\mathbf{v}_h^e, b'_h) dx, \quad (50)$$

$$\text{with } \mathbf{G}(\mathbf{v}_h^e, b_h, v_g) := \mathbf{F}(\mathbf{v}_h^e, b_h) - \mathbf{v}_h^e v_g.$$

*Proof* Multiplying (3a) by $\varphi_h(\cdot, t) \in \mathbb{P}^k(\mathcal{F}_h^e(t))$ and integrating over an arbitrary element $c(t) \in \mathcal{F}_h^e(t)$ gives:

$$\int_{c(t)} \varphi_h \partial_t \mathbf{v}_h^e dx + \int_{c(t)} \varphi_h \partial_x \mathbf{F}(\mathbf{v}_h^e, b) dx = \int_{c(t)} \varphi_h \mathbf{B}(\mathbf{v}_h^e, b'_h) dx. \quad (51)$$

Remark that the use of property (49), as well as the GCL (see [20] for further details), leads to

$$\frac{d}{dt} \int_{c(t)} \mathbf{v}_h^e \varphi_h dx = \int_{c(t)} \varphi_h \partial_t \mathbf{v}_h^e dx + \int_{c(t)} \varphi_h \partial_x (\mathbf{v}_h^e v_g) dx,$$

we may reformulate (51) as follows:

$$\frac{d}{dt} \int_{c(t)} \mathbf{v}_h^e \varphi_h dx + \int_{c(t)} \varphi_h \partial_x \mathbf{G}(\mathbf{v}_h^e, b_h, v_g) dx = \int_{c(t)} \varphi_h \mathbf{B}(\mathbf{v}_h^e, b'_h) dx, \quad (52)$$

where we have set $\mathbf{G}(\mathbf{v}_h^e, b_h, v_g) := \mathbf{F}(\mathbf{v}_h^e, b_h) - \mathbf{v}_h^e v_g$. Integrating by parts again gives (50).

### 3.3 A semi-discrete (non-stabilized) DG-ALE formulation

In this sub-section, we introduce a general semi-discrete in space DG formulation associated with the free-boundary problem (19). The case of a prescribed-motion may be deduced by replacing the discrete equations for the object's motion (53e) by some prescribed data (22). The DG-ALE semi-discrete formulation associated with (19) reads as follows:
for all $t \le T_{\max}$, find $(\mathbf{v}_h^e(\cdot, t), \mathbf{v}_h^i(\cdot, t)) \in \mathbf{P}^k(\mathcal{T}_h^e(t)) \times \mathbf{P}^k(\mathcal{T}_h^i(t))$, $(\chi_-(t), \chi_+(t)) \in ]x_{\text{left}}(t), x_{\text{right}}(t)]^2$ and $\mathcal{X}_G(t) \in ]\chi_-(t), \chi_+(t)] \times \mathbb{R}^2$, such that the following system holds:

$$ 
\left\{
\begin{array}{l}
\displaystyle \frac{d}{dt}(\mathbf{v}_h^e, \varphi_h)_{\mathcal{T}_h^e(t)} + (\mathcal{A}_h(\mathbf{v}_h^e), \varphi_h)_{\mathcal{T}_h^e(t)} = 0, \quad \forall \varphi_h(\cdot, t) \in \mathbb{P}^k(\mathcal{T}_h^e(t)) \text{ s.t. } \frac{d}{dt}\varphi_h(x(X,t),t) = 0, \\
\displaystyle \mathbf{v}_h^e(\cdot, 0) \qquad \qquad \qquad \qquad := \mathbf{p}_{\mathcal{T}_h^e, 0}^k(\mathbf{v}_0^e), \\
\displaystyle \mathbf{v}_h^e|_{\chi_\pm} \qquad \qquad \qquad \qquad = \mathbf{v}_h^i|_{\chi_\pm},
\end{array}
\right.
\tag{53a}
$$

$$
\begin{align}
\eta_h^i(\cdot, t) &:= \mathbf{p}_{\mathcal{T}_h^i(t)}^k \circ F_h(\cdot, t, \mathcal{X}_G(t), \mathcal{X}_G(0), \eta_{\text{id}}), \nonumber \\
q_h^i(\cdot, t) &:= \mathbf{p}_{\mathcal{T}_h^i(t)}^k (\vartheta_G(t) \cdot \mathcal{T}_{G,h}(\cdot, t) + q^i(t)), \nonumber \\
\frac{d}{dt} q^i(t) &:= - (\langle \! \langle f_{1,h} \rangle \! \rangle_h + \langle \! \langle f_{2,h} \rangle \! \rangle_h + \langle \! \langle f_{3,h} \rangle \! \rangle_h), \nonumber \\
q^i(0) &:= q_0^i, \tag{53b}
\end{align}
$$

$$
\left\{
\begin{aligned}
v_g|_{\chi_{\pm}} &:= (\mathfrak{G}_h^k q_h^e|_{\chi_{\pm}} + (\mathbf{v}_G - \omega \mathbf{r}_{G,h}^{\perp}) \cdot \mathbf{n}_h^i) (\mathfrak{G}_h^k \eta_{h|\chi_{\pm}}^e - \mathfrak{G}_h^k \eta_{h|\chi_{\pm}}^i)^{-1}, \\
\frac{d}{dt} \chi_{\pm}(t) &= \tilde{v}_g(X_0^{\pm}, t), \\
\chi_{\pm}(0) &= X_0^{\pm},
\end{aligned}
\right.
\tag{53c}
$$

$$ b_h(\cdot, t) := i_{\mathcal{T}_h^k(t)}^k(b), \tag{53d} $$

$$
\left\{
\begin{aligned}
\frac{d}{dt} \mathcal{X}_G &= \vartheta_G, \\
\left( \mathbb{M}_0 + \mathbb{M}_{a,h} \left[ H_h^i, \mathcal{T}_{G,h} \right] \right) \frac{d}{dt} \vartheta_G &= \begin{pmatrix} -m_0 g e_z \\ 0 \end{pmatrix} - \rho \mathfrak{G}_{\mathcal{I}(t)}^{h,n_g} \left[ (f_{1,h}^* + f_{3,h}^*) \frac{\mathcal{T}_{G,h}^*}{H_h^i} \right], \\
\mathcal{X}_G(0) &= (X_0, Z_0, 0), \\
\vartheta_G(0) &= (u_G^0, w_G^0, \omega_0),
\end{aligned}
\right.
\tag{53e}
$$

where:

(i) the discrete nonlinear operator $\mathcal{A}_h$ in (53a) is defined by

$$ (\mathcal{A}_h(\mathbf{v}_h^e), \varphi_h)_{\mathcal{T}_h^e(t)} := -(\mathbf{G}(\mathbf{v}_h^e, b_h, v_g), \partial_x^h \varphi_h)_{\mathcal{T}_h^e(t)} + \langle \mathbf{G}^*, \varphi_h \rangle_{\partial \mathcal{T}_h^e(t)} - (\mathbf{B}(\mathbf{v}_h^e, b'_h), \varphi_h)_{\mathcal{T}_h^e(t)}, \tag{54} $$

where $\mathbf{G}^*$ is an interface numerical flux which aims at approximating $\mathbf{F}(\mathbf{v}, b) - v_g \mathbf{v}$ at an interior element boundary, which is moving with velocity $v_g$. This formulation is built upon the local formulation shown in Proposition 2, by gathering the local contributions of the residuals, and introducing $\mathbf{G}^* := \mathbf{F}^* - v_g \mathbf{v}^*$ as an approximate interface numerical flux, where $\mathbf{F}^*$ and $\mathbf{v}^*$ are also interface fluxes, respectively consistent with $\mathbf{F}$ and $\mathbf{v}$, and computed with the Lax-Friedrichs formulae:

$$ \mathbf{F}^*(\mathbf{v}_R, \mathbf{v}_L, b_R, b_L) := \frac{1}{2} (\mathbf{F}(\mathbf{v}_R, b_R) - \mathbf{F}(\mathbf{v}_L, b_L) - \sigma(\mathbf{v}_R - \mathbf{v}_L)), \tag{55} $$

$$ \mathbf{v}^*(\mathbf{v}_R, \mathbf{v}_L, b_R, b_L) := \frac{1}{2} \left( \mathbf{v}_R + \mathbf{v}_L - \frac{1}{\sigma} (\mathbf{F}(\mathbf{v}_R, b_R) - \mathbf{F}(\mathbf{v}_L, b_L)) \right), \tag{56} $$

with $\sigma := \max_{c \in \mathcal{T}_h^e} \sigma_c$ and

$$ \sigma_c := \max_{m \in [1, k+1]} \left( |\overline{u_m^e}^c - \overline{v_m^e}^c| + \sqrt{g \overline{H_m^e}^c} \right). \tag{57} $$

(ii) the first equation in (53b) offers a way to compute a high-order broken polynomial approximation of the specific part of the object's underside, which projection along the horizontal line at time $t$ identifies to
$\mathcal{I}(t)$. More precisely, $F_h(\cdot, t, \mathcal{X}_G(t), \mathcal{X}_G(0), \eta_{\text{id}})$ is a discrete nonlinear operator such that, for all $x \in \mathcal{I}(t)$:

$$F_h(x, t, \mathcal{X}_G(t), \mathcal{X}_G(0), \eta_{\text{id}}) := \tilde{F}(X_h, t, \mathcal{X}_G(t), \mathcal{X}_G(0), \eta_{\text{id}}), \quad (58)$$

where $X_h$ is an approximation of the following nonlinear equation's unique root, obtained by some Newton's iterations:

$$\frac{x - x_G(t) + \sin(\theta(t))(\eta_{\text{id}}(X_h) - Z_G)}{\cos(\theta(t))} + X_G - X_h = 0. \quad (59)$$

With $\eta_h^i$ in hands, one can compute

$$H_h^i := \eta_h^i - p_{\mathcal{T}_h^i(t)}^k(b), \quad (60)$$

(iii) the second and third equations in (53b) allow to compute an approximation of the discharge in the interior domain, through the evaluation of a purely geometrical term, together with a time-dependent term obtained as the solution of a nonlinear ordinary differential equation. Specifically, we set:

$$\mathbf{r}_{G,h}(\cdot, t) := \begin{pmatrix} \cdot - x_G(t) \\ \eta_h^i(\cdot, t) - z_G(t) \end{pmatrix}, \quad \mathcal{T}_{G,h}(\cdot, t) := \begin{pmatrix} -\mathbf{r}_{G,h}^{\perp}(\cdot, t) \\ \frac{1}{2} |\mathbf{r}_{G,h}(\cdot, t)|^2 \end{pmatrix}, \quad (61)$$

and the discrete versions of the right-hand sides are defined as follows

$$\begin{align}
f_{1,h} &:= \mathfrak{G}_h^k \circ p_{\mathcal{T}_h^i(t)}^k (u_h^i q_h^i) + g H_h^i \partial_x^h \eta_h^i, \notag \\
f_{2,h} &:= \frac{d}{dt} \vartheta_G \cdot \mathcal{T}_{G,h}, \tag{62} \\
f_{3,h} &:= \vartheta_G^T \mathbb{M}_{G,h} \vartheta_G, \notag
\end{align}$$

where we use a discrete version of (9) to evaluate the term $\partial_t \mathcal{T}_{G,h}$ that appears in $f_{3,h}$, with

$$\mathbb{M}_{G,h} := \begin{pmatrix} \mathbf{e}_x \cdot \mathbf{n}_{\text{id}} & 0 & -\mathbf{r}_{G,h}^{\perp} \cdot \mathbf{n}_{\text{id}} \\ 1 & 0 & 0 \\ -\mathbf{r}_{G,h}^{\perp} \cdot \mathbf{n}_{\text{id}} & 0 & -(\mathbf{e}_z \cdot \mathbf{r}_{G,h})(\mathbf{r}_{G,h}^{\perp} \cdot \mathbf{n}_{\text{id}}) \end{pmatrix}.$$

We also recall that the discrete version of the $\mathcal{I}(t)$-averaging operator is provided in (32),

(iv) the discrete contact-points velocity (53c) is obtained from (41), using the expression of the time-derivative (18). Once this velocity is known, the updated computation domain's macro-partition $\Omega_t = \mathcal{E}(t) \cup \mathcal{I}(t)$ may be updated,

(v) the discrete counterpart of the added-mass-inertia matrix, denoted by $\mathbb{M}_{a,h}[H_h^i, \mathcal{T}_{G,h}]$, is defined as follows:

$$\mathbb{M}_{a,h}[H_h^i, \mathcal{T}_{G,h}] := \mathfrak{G}_{\mathcal{I}(t)}^{h,n_g} \left( \frac{\mathcal{T}_{G,h}^* \otimes \mathcal{T}_{G,h}^*}{H_h^i} \right). \quad (63)$$

This matrix may simply be denoted by $\mathbb{M}_{a,h}$ in what follows.

**Remark 12** The boundary conditions on $\partial \Omega_t$ are weakly enforced through the numerical fluxes $\mathbf{G}^*$. In particular, one may enforce any usual boundary condition for the NSW equations, including inflow and outflow conditions relying on local Riemann invariants, periodic conditions or solid-wall conditions. The coupling condition $\mathbf{v}_{h|\chi_{\pm}}^e = \mathbf{v}_{h|\chi_{\pm}}^i$ is also weakly enforced through the definition of $\mathbf{G}^*$ at the contact points $\chi_{\pm}(t)$. It is worth emphasizing here that such a coupling between $\mathbf{v}^e$ and $\mathbf{v}^i$ appears as straightforward and natural, thanks to the discontinuous nature of the chosen approximation spaces and the DG approximations.

**Remark 13** In practice, the positivity of $H_h^i$ (defined in (60)) may be ensured thanks to the sizing introduced in Remark 3, and to the assumptions that $\mathcal{M}_G$ remains close to its initial position and that $\mathcal{O}_{bj}$ can not be entirely immersed. Indeed, considering a flat bottom for the sake of simplicity, assuming that $H_h^{i,0} \gg 0$, together with $d_o \ll H_0$ and considering that the object can not be completely immersed, with an object's motion $\mathcal{X}_G$ such that $\mathcal{M}_G$ remains close to its initial location, then we necessarily have $H_h^i(\cdot, t) > 0$ for $t \in [0, T_{\text{max}}]$. According to (58), this entails that $F_h(\cdot, t, \mathcal{X}_G(t), \mathcal{X}_G(0), \eta_{\text{id}}) > 0$. This can be also verified when the bathymetry is not flat by assuming that $d_o \ll H_b$.
### 3.4 Time-marching algorithms

For a given final computational time $T_{max} > 0$, we consider a partition $(t^n)_{0 \le n \le N}$ of the time interval $[0, T_{max}]$ with $t^0 := 0$, $t^N := T_{max}$ and $t^{n+1} - t^n := \Delta t^n$. For any sufficiently regular function $w$ depending on time, we set $w^n := w(t^n)$ and in what follows, such a "superscript $n$" notation may be used with any time-varying entity, evaluated at discrete time $t^n$. In particular, we note:

$$ \mathcal{E}^n := \mathcal{E}(t^n), \quad \mathcal{I}^n := \mathcal{I}(t^n), \quad \chi_{-}^n := \chi_{-}(t^n), \quad \chi_{+}^n := \chi_{+}(t^n), \quad \mathcal{F}_h^n := \mathcal{F}_h(t^n), \quad \mathcal{F}_h^{e,n} := \mathcal{F}_h^e(t^n), $$

and so on, together with similar notations for the main unknowns of the problem:

$$ \mathbf{v}_h^{e,n} := \mathbf{v}_h^e(\cdot, t^n), \quad \mathbf{v}_h^{i,n} := \mathbf{v}_h^i(\cdot, t^n), \quad \chi_{\pm}^n := \chi_{\pm}(t^n), \quad \mathbf{\chi}_G^n := \mathbf{\chi}_G(t^n), \quad \mathbf{\vartheta}_G^n := \mathbf{\vartheta}_G(t^n). $$

When fully-discrete formulations are considered, the time-stepping is carried out with explicit Strong-Stability-Preserving Runge-Kutta (SSP-RK) schemes [18,37]. Anticipating on the description of our stability-enforcement operator in the next section, which relies on both DG approximations on mesh elements $c^n \in \mathcal{F}_h^{e,n}$ and FV schemes on the subcells $s_m^{c_i,n} \in \mathcal{F}_c^n$, the time step $\Delta t^n$ is computed adaptively using the following CFL condition:

$$ \Delta t^n = \frac{\min_{c^n \in \mathcal{F}_h^{e,n}} \left( \frac{h_c^n}{2k+1}, \min_{s_m^{c_i,n} \in \mathcal{F}_c^n} |s_m^{c_i,n}| \right)}{\sigma}, \quad (64) $$

where $\sigma$ is defined in (57). The same SSP-RK method is used for the time discretization of the IVP involved in (53).

### 3.5 Stabilization of hyperbolic DG through a posteriori LSC method

It is well-known that the solutions of nonlinear hyperbolic equations may develop singularities in finite time and high-order DG methods usually suffer from a lack of robustness, possibly inducing some stability issues. For the problem under study, problems are generally observed in the vicinity of the solution's discontinuities or where the water-height vanishes. To alleviate such a limitation, we equip our DG scheme with the a posteriori Local Subcell Correction (LSC) method introduced in [19], based on the general guidelines provided in [40]. As this stabilization process, adapted to the DG-ALE framework, has been investigated very recently in [20], we only provide general guidelines:

‡ 1 starting from an admissible solution $\mathbf{v}_h^{e,n} \in \mathbf{P}^k(\mathcal{F}_h^{e,n})$ at discrete time $t^n$, compute a candidate solution $\mathbf{v}_h^{e,n+1} \in \mathbf{P}^k(\mathcal{F}_h^{n+1})$ using the uncorrected DG-ALE scheme (53a)-(54), together with a SSP-RK time discretization,

‡ 2 for all cells $c_i^{n+1} \in \mathcal{F}_h^{e,n+1}$, compute the sub-mean values:

$$ \mathbb{P}^0(\mathcal{F}_{c_i}^{e,n+1}) \ni \bar{\mathbf{v}}_{c_i}^{e,n+1} = \pi_{\mathcal{F}_{c_i}^{e,n+1}}(\mathbf{v}_{c_i}^{e,n+1}), $$

and check their admissibility using the specified admissibility criteria,

‡ 3 for a given subcell $s_m^{c_i,n+1}$ in need of stabilization, the related DG reconstructed interface fluxes $\hat{\mathbf{G}}_{m \pm \frac{1}{2}}$ defined in [19] (which are formally used to compute the predictor $\mathbf{v}_h^{e,n+1}$), are replaced by the FV corrected fluxes $\mathcal{G}_{m \pm \frac{1}{2}}^{l/r}$ into the update process. The whole set of substituting rules is not recalled here,

‡ 4 for the cells $c_i(t)$ where some corrections occurred, new updated (and corrected) sub-mean values are now available and used to produce a new high-order polynomial candidate solution, still denoted by $\mathbf{v}_h^{e,n+1}$ for the sake of convenience:

$$ \mathbb{P}^k(c_i^n) \ni \mathbf{v}_{c_i}^{e,n+1} = \pi_{\mathcal{F}_{c_i}^{n+1}}(\bar{\mathbf{v}}_{c_i}^{e,n+1}). $$

The process may go further in time after checking that this new candidate is admissible.

Let us also recall that the definition of an *admissible* solution relies on two criteria: one for the *Physical Admissibility Detection* (PAD) and the other one to address the occurrence of spurious oscillations, called *Subcell Numerical Admissibility Detection* (SubNAD). The whole detection-correction-projection iterative
process may be conveniently summarized through the application of a global reconstruction operator, denoted as follows:

$$ \Lambda_h^{k,n} : \mathbf{P}^k(\mathcal{F}_h^{e,n}) \to \mathbf{P}^k(\mathcal{F}_h^{e,n}), \\ \mathbf{v}_h^{e,n} \mapsto \Lambda_h^{k,n}(\mathbf{v}_h^{e,n}), \qquad (65) $$

where the resulting broken-polynomial $\Lambda_h^{k,n}(\mathbf{v}_h^{e,n})$ satisfies all the requested admissibility criteria. When higher-order time-discretization schemes are used, this stabilization operator may be applied at each sub-step, see [20] for further details. In what follows, the resulting stabilized DG-ALE method is denoted by DG-ALE-LSC method.

## 3.6 RK-DG-ALE-LSC fully-discrete formulation

Let us completely specify the fully-discrete formulation obtained by considering (53) together with a first-order Euler time-marching algorithm for the sake of simplicity. Any higher-order in time RK-DG-ALE-LSC formulation based on §3.4 may be straightforwardly deduced from this lowest-order in time one by adapting accordingly the various time-stages. Assuming that all the needed quantities are available at discrete-time $t^n$, the first-order in time fully-discrete formulation associated with (19) reads as follows:

find $(\mathbf{v}_h^{e,n+1}, \mathbf{v}_h^{i,n+1}) \in \mathbf{P}^k(\mathcal{F}_h^{e,n+1}) \times \mathbf{P}^k(\mathcal{F}_h^{i,n+1}), \mathcal{X}_G^{n+1} \in [\chi_-^{n+1}, \chi_+^{n+1}] \times \mathbb{R}^2$ and $\chi_\pm^{n+1} \in [x_{\text{left}}^{n+1}, x_{\text{right}}^{n+1}]^2$, such that the following system holds:

$$ \chi_{\pm}^{n+1} - \chi_{\pm}^{n} = \Delta t^n \tilde{v}_{\pm}^{n} (X_{0}^{\pm}), \qquad (66a) $$

$$ b_h^{n+1} = i_h^k \mathcal{F}_h^{n+1}(b), \qquad (66b) $$

$$ \begin{gathered} \mathcal{X}_G^{n+1} - \mathcal{X}_G^n = \Delta t^n \vartheta_G^n, \\ \mathcal{A}_G^n = (M_0 + M_{a,h}^n)^{-1} \left\{ \begin{pmatrix} -m_0 g e_z \\ 0 \end{pmatrix} - \rho \mathfrak{G}_{\mathcal{I}^n}^{h,n_g} \left( (f_{1,h}^{\star,n} + f_{3,h}^{\star,n}) \frac{\mathcal{T}_{G,h}^{\star,n}}{H_h^{i,n}} \right) \right\}, \end{gathered} \qquad (66c) $$

$$ \begin{gathered} \vartheta_G^{n+1} - \vartheta_G^n = \Delta t^n \mathcal{A}_G^n, \\ \begin{cases} \eta_h^{i,n+1} = p_{\mathcal{F}_h^{i,n+1}}^k \circ F_h^{n+1}(\mathcal{X}_G^{n+1}; \mathcal{X}_G^0, \eta_{\text{lid}}), \\ q_h^{i,n+1} = p_{\mathcal{F}_h^{i,n+1}}^k (\vartheta_G^{n+1} \cdot \mathcal{T}_{G,h}^{n+1} + \underline{q}^{i,n+1}), \\ \underline{q}^{i,n+1} - \underline{q}^{i,n} = -\Delta t^n (\langle\langle f_{1,h}^n \rangle\rangle_h + \langle\langle f_{2,h}^n \rangle\rangle_h + \langle\langle f_{3,h}^n \rangle\rangle_h), \end{cases} \end{gathered} \qquad (66d) $$

$$ \begin{gathered} f_{1,h}^n := \mathfrak{G}_h^k \circ p_{\mathcal{F}_h^{i,n}}^k (u_h^{i,n} q_h^{i,n}) + g H_h^{i,n} \partial_x \eta_h^{i,n}, \\ f_{2,h}^n := \mathcal{A}_G^n \cdot \mathcal{T}_{G,h}^n, \\ f_{3,h}^n := (\vartheta_G^n)^T M_{G,h}^n \vartheta_G^n, \end{gathered} $$

$$ \begin{gathered} \mathbf{v}_h^{e,n+1} = \Lambda_h^{k,n+1} (\mathbf{v}_h^{e,n} - \Delta t^n \mathcal{A}_h (\mathbf{v}_h^{e,n})), \\ \mathbf{v}_{h|_{\chi_{\pm}}}^{e,n+1} = \mathbf{v}_{h|_{\chi_{\pm}}}^{i,n+1}, \end{gathered} \qquad (66e) $$

$$ \mathbf{v}_{g|\chi_{\pm}}^{n+1} = (\mathfrak{G}_h^k q_{h|\chi_{\pm}}^{e,n+1} + (\mathbf{v}_G^{n+1} - \omega^{n+1} \mathbf{r}_{G,h}^{\perp,n+1}) \cdot \mathbf{n}_h^{i,n+1}) (\mathfrak{G}_h^k \eta_{h|\chi_{\pm}}^{e,n+1} - \mathfrak{G}_h^k \eta_{h|\chi_{\pm}}^{i,n+1})^{-1}, \qquad (66f) $$

where all the discrete-in-time quantities and operators are obtained from (53)-(54)-(62) at $t = t^n$ and with
a first iteration initialized with:

$$ \mathbf{v}_h^{e,0} := \mathrm{p}_{\mathcal{T}_h^{e,0}}^k (\mathbf{v}_0^e), \quad \text{with } \mathbf{v}_0^e \in (H^s(\mathcal{E}_0))^2, \qquad (67a) $$

$$ \mathcal{X}_G^0 := \mathcal{X}_G(0), \qquad (67b) $$

$$ \vartheta_G^0 := \vartheta_G(0), \qquad (67c) $$

$$ \chi_{\pm}^0 := X_0^{\pm}, \qquad (67d) $$

$$ \underline{q}^{i,0} := q_0^i. \qquad (67e) $$

**Remark 14** The DG scheme (66e) is written in the operator form, for the sake of simplicity. In practice, it may be convenient either to express the corresponding scalar-products at the initial time on the reference configuration, or to use the equivalent FV formulation on the sub-cells. It is also implicitly assumed that $\mathcal{G}_h^k \eta_{h|_{\chi_{\pm}}}^{e,n+1} - \mathcal{G}_h^k \eta_{h|_{\chi_{\pm}}}^{i,n+1} \neq 0$, so that (66f) is defined.

**Remark 15** This fully-discrete formulation is only provided for the first-order Euler time-marching algorithm for the sake of simplicity, and this can be easily extended to any higher-order stable time-marching method written as a convex combination of this first-order scheme, like the SSP-RK3 method used in the next section.

## 3.7 Some properties

The global fully-discrete formulation (66) enjoys several properties, which are conveniently gathered in this sub-section.

### Invertibility of the discrete added-mass matrix

**Proposition 3** The matrix $M_{a,h} [H_h^i, \mathcal{T}_{G,h}]$ is symmetric and non-negative.

*Proof* The matrix $\mathcal{T}_{G,h}^* \otimes \mathcal{T}_{G,h}^*$ is obviously symmetric and of rank one, with non-negative eigenvalues. From the implicit additional assumption that $H_h^i > 0$ (see Remarks 3 and 13), and since we are using quadrature rules with positive coefficients, it results that $M_{a,h}$ is a symmetric and non-negative matrix.

### Well-balancing (for motionless steady-state)

A motionless steady-state for the system (19) may be defined as follows, for all $t \ge 0$:

$$ v_{g|_{\chi_{\pm}}} = 0, \quad \mathbf{v}^e(\cdot, t) = \begin{pmatrix} \eta^c \\ 0 \end{pmatrix} = \mathbf{v}^c, \quad \mathbf{v}^i(\cdot, t) = \begin{pmatrix} \eta_{\mathrm{lid}} \\ 0 \end{pmatrix}, \quad \vartheta_G(t) = 0, \qquad (68) $$

$$ \text{and} \quad 0 = -m_o g + \rho \int_{\mathcal{I}(t)} (\underline{\mathbf{p}}^i(x, t) - \mathbf{p}_{\mathrm{atm}}) \, dx. \qquad (69) $$

We emphasize that, in contrast with [20] where the case of a stationary obstacle is considered, the formulation (66) accounts for the motion of the floating object and, as a consequence for objects with a symmetric underside’s profile, the additional identity (69) is required to balance the action of gravity with the Archimedean force and is enough, when combined to (68), to ensure that $\mathcal{A}_G(t) = 0$. The reader is referred to D for additional details concerning (69). At the discrete level, the preservation of (68)-(69) translates into the following result:

**Proposition 4** The fully-discrete formulation (66) preserves the motionless steady-states (68)-(69), provided that the integrals of (66) are exactly computed (at motionless steady-states). Precisely, under this assumption, for all $n \ge 0$,

$$ \left\{ \begin{array}{l} v_{g|_{\chi_{\pm}}}^n = 0, \quad \mathbf{v}_h^{e,n} = \mathbf{v}^c, \quad q_h^{i,n} = 0, \quad \vartheta_G^n = 0, \quad \mathcal{A}_G^n = 0 \\ v_{g|_{\chi_{\pm}}}^{n+1} = 0, \quad \mathbf{v}_h^{e,n+1} = \mathbf{v}^c, \quad q_h^{i,n+1} = 0, \quad \vartheta_G^{n+1} = 0, \quad \mathcal{A}_G^{n+1} = 0 \end{array} \right\} \Rightarrow $$
Proof We first emphasize that, from a fully-discrete viewpoint, the identity $\mathcal{A}_G^n = 0$ ensures the balance between the action of gravity with the Archimedean force at time $t^n$. The issue of how one may enforce such equilibrium, for instance for the initial data, is detailed in D. We also remark that, thanks to the weak coupling condition between $\mathcal{E}(t)$ and $\mathcal{I}(t)$ through the DG interface fluxes, proving that (66) preserves steady-states boils down to prove that: (i) the DG-ALE scheme (66e) is well-balanced on $\mathcal{E}(t) = \mathcal{E}^-(t) \cup \mathcal{E}^+(t)$, (ii) the floating object preserves its initial equilibrium. As such well-balancing property for the DG scheme in ALE description (66e) on $\mathcal{E}(t) = \mathcal{E}^-(t) \cup \mathcal{E}^+(t)$ has been established in [20] for the case of a fixed surface obstacle, we directly conclude here that $v_h^{e,n+1} = v_c$ and we only examine the propagation in discrete time of the object's equilibrium.

From (66a) and $v_{g|\chi_\pm}^n = 0$, we get $\chi_\pm^{n+1} = \chi_\pm^n$ (and $\mathcal{I}^{n+1} = \mathcal{I}^n$), and from (66c) together with $\vartheta_G^n = 0$ and $\mathcal{A}_G^n = 0$, we get that $\mathcal{X}_G^{n+1} = \mathcal{X}_G^n$ and $\vartheta_G^{n+1} = \vartheta_G^n = 0$: the object does not move, then $\eta_h^{i,n+1} = \eta_h^{i,n}$ and as $b_h^{n+1} = i_{\mathcal{F}_h^{i,n}}(b)$, we also have $H_h^{i,n+1} = H_h^{i,n}$. Moreover, we have at time $t^n$:

$$f_{1,h}^n = gH_h^{i,n} \partial_x^h \eta_h^{i,n}, \quad f_{2,h}^n = 0, \quad f_{3,h}^n = 0,$$

so that

$$\langle\langle f_{1,h}^n \rangle\rangle_h = g \overline{H_h^i}^{-1} \int_{\mathcal{I}^n} \partial_x^h \eta_h^{i,n} dx, \quad (70)$$

$$= g \overline{H_h^i}^{-1} \left[ \eta_h^{i,n} \right]_{\mathcal{I}^n} = 0, \quad (71)$$

where the first and second equalities come from the facts that both $\partial_x^h \eta_h^{i,n}$ and $\eta_h^{i,n}$ belong to $\mathbb{P}^k(\mathcal{T}_h^{i,n})$ and that the chosen quadrature rules and discrete gradient are exact for polynomials, and $\eta_h^{i,n}$ is continuous; and the third one comes from $\eta_{h,|\chi_\pm}^{i,n} = \eta^c$ at steady-state. As a consequence, from (66d), we get $\underline{q}^{i,n+1} = \underline{q}^{i,n}$ and then $\underline{q}_h^{i,n+1} = \underline{q}_h^{i,n} = 0$. Straightforward computations also show that

$$\mathcal{S}_{\mathcal{I}^{n+1}}^{h, n_g} \left( (f_{1,h}^{*,n+1} + f_{3,h}^{*,n+1}) \frac{\mathcal{T}_{G,h}^{*,n+1}}{H_h^{1,n+1}} \right) = \mathcal{S}_{\mathcal{I}^n}^{h, n_g} \left( (f_{1,h}^{*,n} + f_{3,h}^{*,n}) \frac{\mathcal{T}_{G,h}^{*,n}}{H_h^{1,n}} \right),$$

which implies that $\mathcal{A}_G^{n+1} = \mathcal{A}_G^n = 0$, and (66f) shows that $v_{g|\chi_\pm}^{n+1} = 0$.

**Remark 16** This well-balanced property can be extended to any higher-order time discretization that can be expressed as a convex combination of first-order Euler schemes.

### Water-height positivity

The a posteriori LSC method offers the ability to compute again, in a more robust way, the time evolution of the flow variables only in the particular mesh elements (and even in the particular sub-cell) for which a straightforward computation (with the uncorrected DG formulation) would generate negative values for the water-height. It would be difficult to state and prove a general robustness result without entering into the details of the substitution and correction rules, which are described in [19]. However, the overall effect of the a posteriori LSC method may be summarized within the following formal result:

Assuming that $\forall c_i \in \mathcal{T}_h^{e,n}, \forall s_m^{c_i} \in \mathcal{T}_{c_i}, \overline{H_m^e}^{c_i,n} \ge 0$, the fully-discrete formulation (66) ensures that

$$\forall c_i \in \mathcal{T}_h^{e,n+1}, \forall s_m^{c_i} \in \mathcal{T}_{c_i}, \overline{H_m^e}^{c_i,n+1} \ge 0,$$

under the CFL-like condition (64).

The proof of such a result would rely on the fact that in the worse scenario, all the interface reconstructed fluxes would be replaced by a lowest-order FV interface flux, which proves to be robust.
## 3.8 About high-order approximations and convergence

While the numerical algorithm (53) is defined (and stabilized) for any arbitrary high-order polynomial approximations, we recall that we aim at approximating the solutions of some coupled problems, relying on various approximation strategies. In particular, the description of the flow beneath the object reduces to a nonlinear ODE and the mesh-related approximations in this domain reduces to polynomial interpolations and projections: the mathematical assumptions on the initial data, together with the well-posedness result of [22] and recalled in §2, ensure that the mesh convergence properties are optimal.

The coupling with the free-surface domains is achieved relying on the boundary conditions approximations, regarded as transmission conditions, modulating the global convergence properties, as the convergence behaviors associated with classical boundary conditions (Dirichlet, Neumann) are already investigated in [40, 19].

In the design of robust-high-order schemes in the numerical modeling of such fully nonlinear coupled problems, and keeping practical applications in mind, we'd rather expect to benefit, in practice, from the ability to use very coarse meshes, while preserving sub-mesh resolution, and robustness facing strong non-linearities and singularities.

# 4 Numerical validations

In this section, we provide several numerical assessments of the formulation (66), using fourth-order local polynomials ($k = 3$) and with the third-order SSP-RK time-marching scheme instead of the first-order one. Unless stated otherwise, we enforce homogeneous Neumann boundary conditions at the right and left exterior boundaries $x_{\text{right}}$ and $x_{\text{left}}$.

## *Plotting options*

In what follows, unless stated otherwise, we choose to display sub-mean values instead of point-wise values of the polynomial approximations, as it allows to precisely illustrate the sub-cell resolution of the scheme. Also, in order to minimize the total number of figures, we sometimes choose to display both the free-surface $\eta$ and the discharge $q$ on the same graphics. As the magnitudes of these two flow quantities are generally not similar, instead of directly plotting $q$, we choose to display the rescaled and translated quantity $\tilde{q} := \frac{q}{H_0} + H_0$.

## *Body's parameterization and initial balance*

We choose to consider the same floating object with an elliptic underside's profile for all the test-cases. Such a choice is of course arbitrary and may be adapted to alternative shapes. In order to avoid unnecessary repetitions, we refer the reader to the Appendix B of the companion paper [20] for explicit formulae. Also, when a floating object is embedded into a motionless steady state-equilibrium (for instance for the validation of the well-balancing property, or for the definition of initial conditions), we rely on the strategy described in D to define the mass $m_o$ of $O_{bj}$ as follows:

$$m_o := \rho \int_{\mathcal{I}(0)} (H_0 - i_{\mathcal{H}}^k \circ \eta_{\text{id}}) dx. \qquad (72)$$

## *Mass and energy balances*

Whenever needed, the total mass of the fluid in $\mathcal{E}(t) \cup \mathcal{I}(t)$ may be computed by approximating the following:

$$m_{\text{sw}}(t) := \int_{\Omega(t)} H(x, t) dx,$$

while the associated energy quantities may be approximately computed from their definitions in (24)-(25). When these quantities are expected to be conserved over time, this may be checked through the computation of the corresponding errors, for instance for the fluid's mass and total energy:

$$E_{m_{\text{sw}}} = \left| \frac{m_{\text{sw}}(t) - m_{\text{sw}}(0)}{m_{\text{sw}}(0)} \right| \quad \text{and} \quad E_{e_{\text{tot}}} = \left| \frac{e_{\text{tot}}(t) - e_{\text{tot}}(0)}{e_{\text{tot}}(0)} \right|. \qquad (73)$$
## 4.1 Well-balancing for motionless steady-state

In this first test-case, we aim at showing that the method is able to preserve some motionless steady states involving a floating object. Such property has already been proved in the particular case of a stationary object in [20], which also encompass the particular case of a prescribed motion with $\vartheta_G = 0$. Hence, we only focus on the case of a freely floating object, which is placed in a computational domain $\Omega = [-50, 200]$, and with a bathymetry defined as follows:

$$b(x) := \begin{cases} A \left( \sin \left( \frac{\pi (x - x_1)}{75} \right) \right)^2, & \text{if } x_1 \le x \le x_2, \\ \frac{1}{\beta} (x - x_3), & \text{if } x \ge x_3, \\ 0, & \text{elsewhere,} \end{cases} \qquad (74)$$

where $A = 1.5 m$, $\beta = 11$, $x_1 = 12.5 m$, $x_2 = 87.5 m$ and $x_3 = 90 m$. The initial data is defined as

$$\eta_0^e(x) := \max(5, b(x)), \quad \eta_0^i := p_{\mathcal{H}^{i,0}}^k(\eta_{\text{id}}), \quad \text{and} \quad q_0^e = q_0^i = 0.$$

We evolve this initial configuration up to $T_{\text{max}} = 50$ s, with $n_{\text{el}}^e = 50$ and $n_{\text{el}}^i = 10$. The numerical results are shown on Fig. 4 for the free-surface, with a zoom in the vicinity of the contact-points and of the wet-dry transition, showing the corrected and uncorrected subcells on Fig. 5.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Test 4.1 - Preservation of a motionless steady-state - Free-surface elevation at $t = 50s$ for $k = 3$ and $n_{\text{el}}^e = 50$, $n_{\text{el}}^i = 10$.</figcaption></figure>

As expected, the steady-state is preserved for any arbitrary number of iterations, up to the machine accuracy. This result is also obtained for alternate values of $k$ and/or $n_{\text{el}}^e, n_{\text{el}}^i$.

## 4.2 Periodic elementary prescribed motions

This second test-case is dedicated to the individual validations of all three elementary prescribed motions along the three geometrical freedom degrees allowed in a two-dimensional space, with a flat bottom. The computational domain is $\Omega = [-100 m, 200 m]$ and the object is initially located at $(X_G, Z_G) = (50 m, H_0 + 2.5 m)$, with $H_0 = 5 m$, in a motionless steady-state configuration. The initial data are provided as follows:

$$\eta_0^e = H_0, \quad \eta_0^i = i_{\mathcal{H}^{i,0}}^k(\eta_{\text{id}}) \quad \text{and} \quad q_0^e = q_0^i = 0.$$
<figure><img src="image_6.png" /><figcaption>Fig. 5: Test 4.1 - Preservation of a motionless steady-state - Free-surface elevation at $t = 50s$, with a zoom near the floating body (left) and near the wet-dry front (right), showing the corrected and the uncorrected subcells, for $k = 3$, $n_{el}^e = 50$ and $n_{el}^i = 10$.</figcaption></figure>

For the following three different configurations, we set $n_{el}^e = 50$ and $n_{el}^i = 10$ and apply one of the elementary prescribed motion and observe the impact on the flow. To this end, for each test, the corresponding free-surface elevation and discharge are shown at various time-step ($t = 3T + \frac{T}{4}$, $3T + \frac{T}{2}$ and $3T + \frac{3T}{4}$, on the left, for the chosen time-period $T$), together a zoom showing the motion of the mesh in the vicinity of the contact-points $\chi_{\pm}$ (on the right).

### Pure heaving

A periodic purely vertical motion may be prescribed as follows:

$$ \mathcal{X}_G(t) := \left( x_G(t) = 50, z_G(t) = H_0 + 3 - \frac{1}{2} \cos(\frac{2\pi t}{T}), \theta(t) = 0 \right). $$

with $T = 15s$ and the corresponding results are shown in Fig. 6.

### Pure surging

A periodic purely horizontal motion may be prescribed as follows:

$$ \mathcal{X}_G(t) := \left( x_G(t) = 48 + 2 \cos(\frac{2\pi t}{T}), z_G(t) = H_0 + 2.5, \theta(t) = 0 \right), $$

with $T = 10s$. The free-surface elevation and discharge are displayed on Fig. 7.

### Pure pitching

A periodic purely rotational motion may be prescribed as follows:

$$ \mathcal{X}_G(t) := \left( x_G(t) = 50, z_G(t) = H_0, \theta(t) = \frac{\pi}{25} \sin(\frac{2\pi t}{T}) \right), $$

with $T = 8s$. The related results are shown on Fig. 8.

### Convergence study

We do not find any analytical solution for the general problems involving floating objects with a $C^1$-shape in shallow-water. Furthermore, trying to provide some error analysis for such a coupled and nonlinear problem is difficult and out of the scope of this work. We may only try to give a hint about the overall global convergence of our algorithm, and we consider a reference solution for the pure heaving problem studied above, obtained by considering a (robust and highly validated) first-order finite-volume scheme on a very fine mesh, consisting of 12,000 elements in $\mathcal{E}(t)$ and 1,000 in $\mathcal{I}(t)$. Then, in order to define a relevant numerical
<figure><img src="image_7.png" /><figcaption>Fig. 6: Test 4.2 - Elementary prescribed motions: pure heaving - Free-surface elevation and discharge computed at $t = 3T + \frac{T}{4}$, $3T + \frac{T}{2}$ and $3T + \frac{3T}{4}$ (left) with a zoom showing the displacement of the mesh nodes in the vicinity of contact-points $\chi_{\pm}$ (right)</figcaption></figure>

error, we compute the discrete $L^2$ norm of the error between the numerical and reference solutions on a sequence of increasingly refined meshes for the surface elevation in the exterior flow domain $\mathcal{E}(t)$ at time $t$. Using high-order polynomial approximations, we observe as expected that the global algorithm does not converge with the same high-order rates usually observed when studying smooth shallow-water flows without free-boundary or surface obstacle: we observe some numerical orders of convergence ranging between 1 and 2, depending of the chosen polynomial approximation. We also observe that considering a fixed number of mesh elements, the numerical errors decrease with respect to the order of polynomial approximation.

## 4.3 Pure heaving with varying bathymetry and wet-dry transition

In this third test-case, we further investigate the ability of our model to handle a body's prescribed motion that generates the flow motion. To make the test-case more challenging, we also add a varying bathymetry mimicking a plane sloping beach and the water-depth at rest is chosen such that a wet-dry transition occurs, modeling a shoreline. Depending on the underside's shape, such a configuration can be seen as an idealized model for the generation of unidirectional (regular or irregular) waves by a wave-paddle and related issues may be the control problems associated with the design of some specific shapes and motions in order to generate some targeted incoming waves in a basin. The computational domain is set to $\Omega = ]0, 140 \, m[$ and
<figure><img src="image_8.png" /><figcaption>Fig. 7: Test 4.2 - Elementary prescribed motions: pure surging - Free-surface elevation and discharge computed at $t = 3T + \frac{T}{4}$, $3T + \frac{T}{2}$ and $3T + \frac{3T}{4}$ (left) with a zoom showing the displacement of the mesh nodes near $\chi_{\pm}$ (right).</figcaption></figure>

the bathymetry is defined as follows:

$$b(x) := \begin{cases} \frac{1}{\beta}(x - x_{\beta}) & \text{if } x < x_{\beta}, \\ 0 & \text{elsewhere,} \end{cases}$$

with $\beta = 11$ and $x_{\beta} = 65m$. The object is placed at $(X_G, Z_G) = (50m, H_0 + 2.5m)$ with $H_0 = 5m$ and the initial condition is defined as follows:

$$\eta_0^e = \max(H_0 - b, 0) + b, \quad \eta_0^i = \mathcal{I}_h^k(\eta_{\text{id}}) \quad \text{and} \quad q_0^e = q_0^i = 0.$$

A periodic purely vertical motion is prescribed as follows:

$$\mathcal{X}_G(t) := \left( x_G(t) = 50, z_G(t) = H_0 + 2.75 - \frac{1}{2} \cos(\frac{2\pi t}{T}), \theta(t) = 0 \right),$$

with $T = 20s$ and we compute the time-evolution of the global system with $n_{\text{el}}^e = 60$ and $n_{\text{el}}^i = 10$. The corresponding solution is shown in Fig. 9, under the form of several snapshots of the free-surface elevation at times $t = 53s$, $60.5s$ and $66s$ respectively (left), with the corrected and uncorrected subcells (respectively plotted with blue squares and green dots), and a zoom on the wet-dry front (right). We observe that some surface waves are generated from the fluid near the moving object and propagates towards the shoreline, inducing some run-up and the back-propagation of a reflected wave. The computation is performed in a very
<figure><img src="image_9.png" /><figcaption>Fig. 8: Test 4.2 - Elementary prescribed motions: pure pitching - Free-surface elevation and discharge computed for different values of time $t = 3T + \frac{T}{4}$, $3T + \frac{T}{2}$ and $3T + \frac{3T}{4}$ for $k = 3$, $n_{el}^e = 50$ and $n_{el}^i = 10$.</figcaption></figure>

stable way, even with this very low number of mesh element and cubic local interpolation. We emphasize that thanks to the chosen smooth coordinate transformation (see the diffeomorphism definition introduced in [20]), the mesh elements move around the floating body, while they remain stationary near the shoreline and in the dry area.

## 4.4 Nonlinear interactions with a single surface wave

In this fourth test-case, we begin to investigate the computation of free-motions, and we model the object’s oscillations generated by a single surface wave over a flat bathymetry. The computational domain is $\Omega = [-150 \text{ m}, 150 \text{ m}]$ and $M_G$ is placed in equilibrium at $(x_G(0), z_G(0)) = (50 \text{ m}, H_0 + 2.5 \text{ m})$, with $H_0 = 5 \text{ m}$. The surface-wave is defined as follows:

$$ \eta_0^e(x) = H_0 + \frac{A_w}{\cosh(\gamma(x - x_0))^2} \quad \text{and} \quad q_0^e = c_{q2} \frac{\sqrt{g}(\eta_0^e - H_0)}{c_{q1}} H_0^e, \qquad (75) $$

with $A_w = 0.92 \text{ m}$, $x_0 = -80 \text{ m}$, $c_{q1} = 0.1$, $c_{q2} = 0.05$ and $\gamma = \frac{c_{q1}}{\sqrt{\frac{4H_0}{3A_w}}}$, see Fig. 10, and the initial free-surface elevation and discharge beneath the floating body are defined as:

$$ \eta_0^i = i_{\mathcal{F}_h^{i,0}}^k (\eta_{\text{lid}}) \quad \text{and} \quad q_0^i = 0. $$

We observe that part of the wave energy is transmitted to the object as a rightward oriented pressure force, inducing a global displacement which mobilizes all three possible types of motion, see Fig. 11. In order to have a better idea of the global behavior of the body, we show the trajectory of $M_G$ on Fig. 12, under the form of time-series of its spatial coordinates $(x_G, z_G)$, together with the time-series of the pitch angle $\theta$. On this last curve, we observe an oscillating behavior associated with the return to equilibrium (which is not entirely displayed) and this motivates the next test-case.
<figure><img src="image_10.png" /><figcaption>Fig. 9: Test 4.3 - Pure heaving with a varying bathymetry and a wet-dry transition - Free-surface elevation computed at $t = 53s$, 60.5s and 66s (left) with corrected and uncorrected subcells respectively plotted with blue squares and green dots, and a zoom on the wet-dry front (right)</figcaption></figure>

## 4.5 Convergence towards equilibrium

In this fifth test case, we focus on the ability of our model to converge towards a motionless steady state involving a floating object. The computational domain is $\Omega = [-50m, 150m]$, the mean water-depth is $H_0 = 8m$ and we consider the following bathymetry profile:

$$b(x) := \begin{cases} A \sin \left( \frac{(x - x_1) \cdot \pi}{75} \right), & \text{if } x_1 \le x \le x_2, \\ 0, & \text{elsewhere,} \end{cases} \qquad (76)$$

with $x_1 = 12.5m$, $x_2 = 87.5m$ and $A = 2.5$. For the chosen mass defined as follows

$$m_o := 2.19 \times \int_{\mathcal{I}(0)} \rho(H_0 - \eta_0^i(x)) dx, \qquad (77)$$

the equilibrium position would be $\mathcal{X}_G(0) = (50m, H_0 + 0.83m, 0)$. However, we place the object at $\mathcal{X}_G(0) = (50m, H_0 + 2.5m, 0)$, such that it is not initially in equilibrium. As a consequence, an oscillating heaving is induced, associated to the return towards equilibrium. The initial-data are defined as follows:

$$\eta_0^e = H_0, \quad \eta_0^i = i_{\mathcal{T}_h^{i,0}}(\eta_{\text{lid}}) \quad \text{and} \quad q_0^e = q_0^i = 0,$$
<figure><img src="image_11.png" /><figcaption>Fig. 10: Test 4.4 - Nonlinear interactions with a single surface wave - Free-surface elevation and discharge computed at initial time for $k = 3$, $n_{el}^e = 50$ and $n_{el}^i = 10$.</figcaption></figure>

and we set $n_{el}^e = 60$ and $n_{el}^i = 10$. The free-surface elevation is shown in Fig. 13, for several times in the range $[0 \text{ s}, 25 \text{ s}]$. We observe that the instantaneous release of the body generates two wave fronts propagating in both directions. These discontinuities are well-captured by the DG scheme and efficiently stabilized by the a posteriori LSC method. In Fig. 14, we display the oscillations of $z_G(t)$ in the time-range $[0 \text{ s}, 25 \text{ s}]$, highlighting the return to equilibrium (on the left), and a zoom on the free-surface configuration at $t = 25 \text{ s}$ (on the right). On Fig. 15, we display the free-surface elevation with the propagating waves at $t = 4.9 \text{ s}$ in order to highlight the resolution of the mesh and the corrected subcells. To complete the picture, we also investigate the conservation of mass and energy over time, by computing the time evolution of their relative errors (using solid-wall boundary conditions for both left and right exterior boundaries to close the domain). Related time-series are displayed on Fig. 16 and 17 for $n_{el} = 200$, $k = 3$ and we observe that mass and total energy are preserved over time, with an error magnitude of 5E-10.

## 4.6 Nonlinear interactions with a single surface-wave: varying bathymetry and wet-dry front

In this last test-case, we study the full interactions between a floating object and a single surface-wave propagating shoreward. This is our more challenging test-case. As in §4.3, we consider a varying bathymetry modeling a sloping beach and choose $H_0$ such that a wet-dry transition occurs. The computational domain is $\Omega = [-300 \text{ m}, 150 \text{ m}]$ and $M_G$ is placed in equilibrium at $(x_G(0), z_G(0)) = (50 \text{ m}, H_0 + 2.5 \text{ m})$, with $H_0 = 5 \text{ m}$. The surface-wave is defined as follows:

$$ \eta_0^e(x) = H_0 + \frac{A_w}{\cosh(\gamma(x - x_0))} \quad \text{and} \quad q_0^e = c_{q2} \frac{\sqrt{g(\eta_0^e - H_0)}}{c_{q1}} H_0^e, \qquad (78) $$

with $A_w = 0.55 \text{ m}$, $x_0 = -150 \text{ m}$, $c_{q1} = 0.1$, $c_{q2} = 0.05$ and $\gamma = \frac{c_{q1}}{\sqrt{\frac{4H_0}{3A_w}}}$, see Fig. 18, and the initial free-surface elevation and discharge beneath the floating body are defined as:

$$ \eta_0^i = i_{\mathcal{H}_h^{i,0}}^k (\eta_{\text{id}}) \quad \text{and} \quad q_0^i = 0. $$

Furthermore, in this test, the object is allowed to move according to the response to the mechanical action of the incoming wave, and not according to a prescribed motion as in §4.3. We set $n_{el}^e = 70$ and $n_{el}^i = 10$. We show on Fig. 19 the free-surface at several values in the time-range $[16.78 \text{ s}, 150 \text{ s}]$. In Fig. 20, we show the trajectory of $M_G$, under the form of time-series of its spatial coordinates $(x_G, z_G)$ and the pitch angle $\theta$. On Fig. 21, the emphasis is put on the free-surface computed at $t = 28.61 \text{ s}$, with the corrected and uncorrected subcells respectively plotted with blue squares and green dots (left) and a zoom on the shoreline (right). We observe that the body is pushed shoreward by the incoming wave, which is almost entirely transmitted to the shore, generating a run-up, before being reflected by the bathymetry and hit the object again, pushing it seaward. The computations are performed in a very robust way, and we observe that as expected, the a posteriori LSC method only operate in the vicinity of the shoreline.
<figure><img src="image_12.png" /><figcaption>Fig. 11: Test 4.4 - Nonlinear interactions with a single surface wave - Free-surface elevation computed for different values of time in the range [0 s, 35 s] for $k = 3$, $n_{el}^e = 50$ and $n_{el}^i = 10$.</figcaption></figure>

## 5 Conclusion

In this paper, we investigate the construction of an efficient and robust coupled numerical strategy for the simulation of surface waves and floating body interactions in shallow-water and we consider in this work the case of freely floating objects. We recall the main lines of the model derivation, starting from the NSW equations and accounting for the dynamics of the object as a response to the action of the surrounding water, or to some external forces and torque. Then, we introduce a discrete setting involving DG discretizations, as well as FV discretization on subcells and the possible motion of the free-boundaries is described through an ALE strategy. The FV discrete setting on subcells is used as a ground for an *a posteriori* LSC method which helps to stabilize the computations and enforce some nonlinear stability. The resulting
<figure><img src="image_13.png" /><figcaption>Fig. 12: Test 4.4 - Nonlinear interactions with a single surface wave - Time-series for $\mathcal{X}_G = (x_G, z_G, \theta)$ (from left to right) in the range of time $[0 \text{ s}, 35 \text{ s}]$, for $k = 3$, $n_{\text{el}}^e = 50$ and $n_{\text{el}}^i = 10$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13: Test 4.5 - Convergence towards equilibrium - Free-surface elevation computed for different values of time in the range $[0 \text{ s}, 25 \text{ s}]$ for $k = 3$, $n_{\text{el}}^e = 60$ and $n_{\text{el}}^i = 10$.</figcaption></figure>

global algorithm gathers all these numerical ingredients and enjoys some welcome properties: well-balancing for motionless steady states, DGCL, positivity of the water-height at the subcell level, invertibility of the discrete added-mass matrix. Additionally, several test-cases assess a good qualitative behavior, as well as the validity of some more quantitative markers, like the global conservation of mass and energy, and an excellent robustness even when using high-order polynomials and low-resolution meshes. The extension towards the two-dimensional horizontal framework with a surface obstacle is currently under investigations, together with the introduction of weakly dispersive effects following [2].
<figure><img src="image_15.png" /><figcaption>Fig. 14: Test 4.5 - Convergence towards equilibrium - showing the variation of $z_G(t)$ in the range of time [0 s, 25 s] (left), showing water-body equilibrium state at $t = 25$ (right), for $k = 3$, $n_{el}^e = 60$ and $n_{el}^i = 10$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15: Test 4.5 - Convergence towards equilibrium - Free-surface elevation computed at $t = 4.9$ s: corrected and uncorrected subcells are respectively plotted with blue squares and green dots, with a zoom on the left wave, for $k = 3$ and $n_{el}^e = 60$ and $n_{el}^i = 10$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16: Test 4.5 - Convergence towards equilibrium - Fluid mass and relative error, for $k = 3$, $n_{el} = 200$ in the range [0 s, 100 s].</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 17: Test 4.5 - Convergence towards equilibrium - Total energy and relative error, for $k = 3$, $n_{el} = 200$ in the range $[0 \text{ s}, 100 \text{ s}]$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18: Test 4.6 - Nonlinear interactions with a single surface-wave: varying bathymetry and wet-dry front - Free-surface elevation at initial time for $n_{el}^e = 70$ and $n_{el}^i = 10$.</figcaption></figure>

## A Cut-off function

The cut-off function $\varphi \in \mathcal{D}(\mathbb{R})$ used in (42) is defined as follows:

$$ \forall x \in \mathbb{R}, \quad \varphi(x) := e \psi_e(\varepsilon_0 x), $$

where

$$ \forall x \in \mathbb{R}, \quad \psi_e(x) := \phi_e(1 - |x|^2), $$

and

$$ \forall t \in \mathbb{R}, \quad \phi_e(t) := \begin{cases} e^{-t} & \text{if } t > 0, \\ 0 & \text{elsewhere.} \end{cases} $$

Note that we have $\text{supp}(\psi_e) \subset \bar{B}(0, 1)$, $\text{supp}(\varphi) \subset [-\frac{1}{\varepsilon_0}, \frac{1}{\varepsilon_0}]$ and $\varepsilon_0$ chosen such that we have $\varphi(x) = 1, \forall |x| \le 1$.

## B Definition of the elliptic obstacle

In this work, we consider a partially immersed obstacle $\mathcal{O}_{bj}$, which center of mass is located at $(x_G, z_G)$ and which boundary is denoted by $\partial \mathcal{O}_{bj}$. Denoted respectively by $a, b$ its major and minor radius, we define $\partial \mathcal{O}_{bj}$ as an ellipse, so that we have:

$$ (x, y) \in \partial \mathcal{O}_{bj} \iff \frac{(x - x_G)^2}{a^2} + \frac{(z - z_G)^2}{b^2} = 1. $$
<figure><img src="image_20.png" /><figcaption>Fig. 19: Test 4.6 - Nonlinear interactions with a single surface-wave: varying bathymetry and wet-dry front - Free-surface elevation for different values of time in the range [16.78 s, 150 s], for $n_{el}^e = 70$ and $n_{el}^i = 10$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20: Test 4.6 - Nonlinear interactions with a single surface wave - Time-series for $\mathcal{X}_G = (x_G, z_G, \theta)$ (from left to right) in the range of time [0 s, 150 s], for $k = 3$, $n_{el}^e = 70$ and $n_{el}^i = 10$.</figcaption></figure>

The underside of the obstacle may be locally parameterized as follows:

$$ \forall x \in \mathcal{I}_{\text{lid}} := [x_G - a, x_G + a], \quad \eta_{\text{lid}}(x) := z_G - b\sqrt{1 - \frac{(x - x_G)^2}{a^2}}. $$
<figure><img src="image_22.png" /><figcaption>Fig. 21: Test 4.6 - Nonlinear interactions with a single surface-wave: varying bathymetry and wet-dry front - Free-surface elevation computed at $t = 28.61$ s for $n_{el}^e = 70$ and $n_{el}^i = 10$: corrected and uncorrected subcells are respectively plotted with blue squares and green dots, with a zoom on the shoreline (right).</figcaption></figure>

Note that denoting $z_G = H_0 + e_0$, we have:

$$X_0^{\pm} := x_G \pm \sqrt{a^2 - \frac{a^2 e_0^2}{b^2}}.$$

### C Newton's second law and added-mass effect

The pressure $\underline{p}^i$ satisfies:

$$\partial_x \underline{p}^i = - \frac{\rho}{H_i} (f_1^* + f_2^* + f_3^*). \qquad (79)$$

Using an integration by parts and the boundary condition $\underline{p}^i = \underline{p}_{atm}$ on $\chi_{\pm}$, we can rewrite (8) as:

$$M_0 \frac{d}{dt} \vartheta_G = - \begin{pmatrix} m_o g e_z \\ 0 \end{pmatrix} + \int_{\mathcal{I}(t)} \left( \partial_x \underline{p}^i \right) \mathcal{T}_G^*, \qquad (80)$$

$$= - \begin{pmatrix} m_o g e_z \\ 0 \end{pmatrix} - \rho \int_{\mathcal{I}(t)} (f_1^* + f_2^* + f_3^*) \frac{\mathcal{T}_G^*}{H_i}, \qquad (81)$$

and by using the definition of $f_2$ in (15), (8) is finally reduced to the following ODE:

$$(M_0 + M_a(H^i, \mathcal{T}_G)) \frac{d}{dt} \vartheta_G = \begin{pmatrix} -m_o g e_z \\ 0 \end{pmatrix} - \rho \int_{\mathcal{I}(t)} (f_1^* + f_3^*) \frac{\mathcal{T}_G^*}{H^i}. \qquad (82)$$

### D Object's mass and inertia

When embedding a partly immersed object in a free-surface flow which is initially in a motionless steady-state, one may desire to place the object in such a way that the whole system "fluid-object" is in equilibrium. In practice, for a given floating structure that comes with its own mechanical properties (center of mass $M_G$, boundary profile $\eta_{lid}$, mass $m_o$ and inertia coefficient $i_o$), this boils down to accurately define the position $\mathcal{X}_G(0)$ which ensures the balancing of Newton's laws. That being said, for our numerical studies, we consider the following "inverse" strategy: for a given object profile $\eta_{lid}$, and a given initial position $\mathcal{X}_G(0)$, we carefully calibrate the mass and inertia $m_o$ and $i_o$ such that the Newton's laws are initially balanced with $\mathcal{A}_G(0) = 0$. Specifically, for an elliptic object with a symmetric profile as those considered in §4, $M_G$ may be placed at $(X_G, Z_G = H_0 + e_0)$ for some chosen $e_0$, with $\theta(0) = 0$. And ensuring the equilibrium boils down to the definition of $m_o$, such that:

$$-m_o g e_z + \int_{\mathcal{I}(0)} (\underline{p}^i - \underline{p}_{atm}) \underline{n}^i = 0. \qquad (83)$$

As $\underline{p}^i$ satisfies the following BVP:

$$\left\{ \partial_x \underline{p}^i = -\rho \left( \frac{1}{H_i} \partial_t q^i + \partial_x \left( \frac{1}{2} \left( \frac{q^i}{H^i} \right)^2 + g \eta^i \right) \right) \quad \text{in } \mathcal{I}(t), \qquad (84a) $$

$$\underline{p}^i = \underline{p}_{atm} \qquad \text{on } \mathcal{E}(t) \cap \mathcal{I}(t), \qquad (84b)$$
we have

$$ \underline{p}^i(x, t) - \mathrm{patm} = -\rho \left\{ \int_{\chi_c(t)}^{x} \frac{\partial_t q^i(x', t)}{H^i(x', t)} dx' + \frac{1}{2} \left( \frac{q^i(x, t)^2}{H^i(x, t)^2} - \frac{q^i(\chi_c(t), t)^2}{H^i(\chi_c(t), t)^2} \right) + g(\eta^i(x, t) - \eta^i(\chi_c(t), t)) \right\}, $$

and therefore at motionless steady-state:

$$ \underline{p}^i(x, 0) - \mathrm{patm} = -\rho g (\eta^i(x, 0) - \eta^c), \qquad (85) $$

so that (83) reduces to the following scalar identity, standing as the definition of the object's mass:

$$ m_o := \rho \int_{\mathcal{I}(0)} (\eta^c - i_{\mathcal{H}}^k(\eta_{\mathrm{lid}})) dx, \qquad (86) $$

and the corresponding inertia coefficient is defined:

$$ i_o := \frac{m_o(a^2 + b^2)}{5}, $$

for the elliptic profiles. At the discrete level, $m_o$ may be equivalently computed from (86) or from the second equation of (66c), for $n=0$ and $\mathcal{A}_G(0)=0$, as soon as the corresponding integrals on $\mathcal{I}(0)$ are exactly computed.

## Declarations

### Data availability

Data sharing not applicable to this article as no dataset was generated or analyzed during the current study.

### Conflicts of interest

The authors declare that they have no conflicts of interest in the present work.

## Acknowledgment

We would like to extend our sincere appreciation to David Lannes for his support and contributions. The second author was supported by the BOURGEONS project, grant ANR-23-CE40-0014-01 of the French National Research Agency (ANR).

## References

1. F. Bassi and S. Rebay. A high-order accurate Discontinuous Finite-Element method for the numerical solution of the compressible Navier-Stokes equations. *J. Comput. Phys.*, 131:267-279, 1997.

2. G. Beck and D. Lannes. Freely floating objects on a fluid governed by the Boussinesq equations. *Annales de l'Institut Henri Poincaré C*, 39(3):575-646, 2022.

3. C. Beels, P. Troch, K. De Visch, J.P. Kofoed, and G De Backer. Application of the time-dependent mild-slope equations for the simulation of wake effects in the lee of a farm of wave dragon wave energy converters. *Renewable Energ*, 35(8):1644-1661, 2010.

4. U. Bosi, A.P. Engsig-Karup, C. Eskilsson, and M. Ricchiuto. A spectral/hp element depth-integrated model for nonlinear wave-body interaction. *Comp. Meth. App. Mech. Eng.*, 348:22-249, 2019.

5. D. Bresch, D. Lannes, and G. Métivier. Waves interacting with a partially immersed obstacle in the Boussinesq regime. *Analysis and PDE*, 14(4):1085-1124, 2021.

6. P. Causin, J.-F. Gerbeau, and F. Nobile. Added-mass effect in the design of partitioned algorithms for fluid-structure problems. *Comput. Methods Appl. Mech. Engrg.*, 194:4506-4527, 2005.

7. Q. Chen and I. Babuska. Approximate optimal points for polynomial interpolation of real functions in an interval and in a triangle. *Comput. Methods Appl. Mech. Engrg*, 128:405-417, 1995.

8. X.N. Chen and S.D. Sharma. A slender ship moving at a near-critical speed in a shallow channel. *J. Fluid Mech.*, 291:263-285, 1995.

9. B. Cockburn and C.-W. Shu. Runge-Kutta Discontinuous Galerkin methods for convection-dominated problems. *J. Sci. Comput.*, 16(3):173-260, 2001.

10. W.E. Cummins. *The impulse response function and ship motions*. Department of the Navy, David Taylor Model Basin, 1962.

11. A.J.-C. de Saint-Venant. Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit. *C.R. Acad. Sci. Paris, Section Mécanique*, 73:147-154, 1871.

12. D. A. Di Pietro and A. Ern. Discrete functional analysis tools for discontinuous Galerkin methods with application to the incompressible Navier-Stokes equations. *Math. Comp.*, 79(271):1303-1330, 2010.

13. J. Donea, A. Huerta, J.-Ph. Ponthot, and A. Rodríguez-Ferran. *Arbitrary Lagrangian-Eulerian Methods, The Encyclopedia of Computational Mechanics*, pages 413-437. Wiley, 2004.
14. R.C. Ertekin, W.C. Webster, and J.V. Wehausen. Waves caused by a moving disturbance in a shallow channel of finite width. *J. Fluid Mech.*, 169:275–292, 1986.

15. A. Feng, Z.M. Chen, and W.G. Price. A rankine source computation for three dimensional wave-body interactions adopting a nonlinear body boundary condition. *Appl. Ocean Res.*, 47(313-321), 2014.

16. E. Godlewski, M. Parisot, J. Sainte-Marie, and F. Wahl. Congested shallow water model: roof modelling in free surface flow. *ESAIM Math. Model. Numer. Anal.*, 52(5):1679 – 1707, 2018.

17. E. Godlewski, M. Parisot, J. Sainte-Marie, and F. Wahl. Congested shallow water model: on floating body. *SMAI Journal of Computational Mathematics*, in press, 2022.

18. S. Gottlieb, C.-W. Shu, and Tadmor E. Strong stability preserving high order time discretization methods. *SIAM Review*, 43:89–112, 2001.

19. A. Haidar, F. Marche, and F. Vilar. A posteriori finite-volume local subcell correction of high-order discontinuous Galerkin schemes for the nonlinear shallow-water equations. *J. Comput. Phys.*, 452:110902, 2022.

20. A. Haidar, F. Marche, and F. Vilar. A robust DG-ALE formulation for nonlinear shallow water interactions with a partially immersed object. preprint - https://hal.archives-ouvertes.fr/hal-03764650v1, 2022.

21. T. H. Havelock. The initial wave resistance of a moving surface pressure. *Proc. R. Soc. A*, 93, 1917.

22. T. Iguchi and D. Lannes. Hyperbolic free boundary problems and applications to wave-structure interactions. *Indiana Univ. Math. J.*, 70:353–464, 2021.

23. M. Isaacson. Fixed and floating axisymmetric structures in waves. *J. Waterw. Port Coastal Ocean Eng.*, 108(2):180–199, 1982.

24. M. Isaacson. Nonlinear effects on fixed and floating bodies. *J. Fluid Mech.*, 120:267–281, 1982.

25. M. Isaacson and S. Bhat. Wave propagation past a pile- restrained floating breakwater. *Int. J. Offshore Polar Eng.*, 8:265–269, 1998.

26. T. Jiang. *Ship Waves in Shallow Water*. Verkehrstechnik, Fahrzeugtechnik. Fortschritt-Berichte VDIReihe, 2001.

27. T. Jiang, R. Henn, and S.D. Sharma. Wash waves generated by ships moving on fairways of varying topography. In *24th Symposium on Naval Hydrodynamics Fukuoka, JAPAN*, 2002.

28. F. John. On the motion of floating bodies I. *Communications on Pure and Applied Mathematics*, 2:13–57, 1949.

29. F. John. On the motion of floating bodies II. simple harmonic motions. *Comm. Pure Appl. Math.*, 3:45–101, 1950.

30. M. Kashiwagi. Non-linear simulations of wave-induced motions of a floating body by means of the mixed Eulerian-Lagrangian method. In *Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science*, volume 214, pages 841–855, 2000.

31. E.V. Koutandos, T.V. Karambas, and C.G. Koutitas. Floating breakwater response to waves action using a Boussinesq model coupled with a 2dv elliptic solver. *J. Waterw. Port Coastal Ocean Eng.*, 130:243–255, 2004.

32. E.E. Kriezi, T. Karambas, P. Prinos, and C. Koutitas. Interaction of floating breakwaters with waves in shallow waters. In *Proc., Int. Conf. IAHR, Beijing*, 2001.

33. D. Lannes. On the dynamics of floating structures. *Annals of PDE*, 2017.

34. C. Lee and J.N. Newman. Computation of wave effects using the panel method. In *Numerical Models in Fluid-Structure Interaction*, 2005.

35. F. Marche, P. Bonneton, P. Fabrie, and N. Seguin. Evaluation of well-balanced bore-capturing schemes for 2d wetting and drying processes. *Internat. J. Numer. Methods Fluids*, 53(5):867–894, 2007.

36. T.F. Ogilvie. Second-order hydrodynamic effects on ocean platforms. In *Proc. Intl. Workshop on Ship and Platform Motions*, ed. R. W. Yeung, University of California, Berkeley, pages 205–265, 1983.

37. C.-W. Shu and S. Osher. Efficient implementation of Essentially Non-Oscillatory shock-capturing schemes. *J. Comput. Phys.*, 77:439–471, 1988.

38. H.G. Sung and S.T. Grilli. Bem computations of 3D fully nonlinear free-surface flows caused by advancing surface disturbances. *Int. J. Offshore Polar Eng.*, 18(292-301), 2008.

39. F. Ursell. Short surface waves due to an oscillating immersed body. *Proc. R. Soc. A*, 220, 1953.

40. F. Vilar. A posteriori correction of high-order discontinuous Galerkin scheme through subcell finite volume formulation and flux reconstruction. *J. Comput. Phys.*, 387:245–279, 2019.

41. A. N. Williams and W. G. McDougal. Flexible floating break-water. *J. Waterw. Port Coastal Ocean Eng.*, 117(5):429–450, 1991.

42. G.Q. Yang, O.M. Faltinsen, and R. Zhao. Wash of ships in finite water depth. In *Proceedings of the FAST 2001, Southampton, UK*, 2001.

43. Y.H. Yu and L. Ye. Reynolds-averaged navier stokes simulation of the heave performance of a two-body floating-point absorber wave energy system. *Computers and Fluids*, 73:104–114, 2013.