# Modelling beaches morphodynamic by Hadamard sensitivity analysis

Ronan Dupont, Frédéric Bouchette, Bijan Mohammadi

► To cite this version:

Ronan Dupont, Frédéric Bouchette, Bijan Mohammadi. Modelling beaches morphodynamic by Hadamard sensitivity analysis. Ocean Modelling, 2024, 189, pp.102370. 10.1016/j.ocemod.2024.102370. hal-04559204

HAL Id: hal-04559204

https://hal.science/hal-04559204

Submitted on 25 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
MODELLING BEACHES MORPHODYNAMIC BY
HADAMARD SENSITIVITY ANALYSIS

R. Dupont<sup>1,2,3,*</sup>, F. Bouchette<sup>1,3</sup>, and B. Mohammadi<sup>2,3</sup>

<sup>1</sup>* GEOSCIENCES-M, Univ Montpellier, CNRS, Montpellier, France.

<sup>2</sup> IMAG, Univ Montpellier, CNRS, Montpellier, France.

<sup>3</sup> GLADYS, Univ Montpellier, Montpellier, France.

*Corresponding author(s). E-mail(s): ronan.dupont@umontpellier.fr;
Contributing authors: frederic.bouchette@umontpellier.fr;
bijan.mohammadi@umontpellier.fr.

**Abstract**

The paper presents a morphodynamic model which can be coupled
with any wave model capable of producing time/spectral averaged wave
quantities. This model based on a wave energy minimization principle
highlights the morphodynamic phenomenology, such as the sandbar cre-
ation. Such a model can be used in solving engineering optimization
problems. It is also developed to illustrate the idea that beach sand
transport can be thought as a non-local phenomenon. We used wave
calculations from SWAN and XBeach in our model, and we compared
the morphodynamic results to LIP and SANDS hydro-morphodynamic
benchmark as well as open-sea simulations. Using supplementary math-
ematical development, we improved the minimization method using the
Hadamard derivative.

**Keywords.** Hydro-morphodynamics, Optimization, Model validation, Coastal, Variational approach, Energy minimization, Optimal transport, Waves.
MODELLING BEACHES MORPHODYNAMIC BY
HADAMARD SENSITIVITY ANALYSIS

R. Dupont¹,²,³*, F. Bouchette¹,³, and B. Mohammadi²,³

¹\* GEOSCIENCES-M, Univ Montpellier, CNRS, Montpellier, France.

² IMAG, Univ Montpellier, CNRS, Montpellier, France.

³ GLADYS, Univ Montpellier, Montpellier, France.

\*Corresponding author(s). E-mail(s): ronan.dupont@umontpellier.fr;
Contributing authors: frederic.bouchette@umontpellier.fr;
bijan.mohammadi@umontpellier.fr.

Abstract

The paper presents a morphodynamic model which can be coupled
with any wave model capable of producing time/spectral averaged wave
quantities. This model based on a wave energy minimization principle
highlights the morphodynamic phenomenology, such as the sandbar cre-
ation. Such a model can be used in solving engineering optimization
problems. It is also developed to illustrate the idea that beach sand
transport can be thought as a non-local phenomenon. We used wave
calculations from SWAN and XBeach in our model, and we compared
the morphodynamic results to LIP and SANDS hydro-morphodynamic
benchmark as well as open-sea simulations. Using supplementary math-
ematical development, we improved the minimization method using the
Hadamard derivative.

**Keywords.** Hydro-morphodynamics, Optimization, Model validation, Coastal, Variational approach, Energy minimization, Optimal transport, Waves.
# 1 Introduction

Morphodynamic models are generally very complex and highly parameterized. They separately solve the physical equations of hydrodynamics and morphodynamics at a very small scale of the order of second in time and of the wave length in space. The OptiMorph model that we presented in Cook (2021) and Dupont et al. (2023) proposes a more global approach based on an optimization principle.

The optimization theory is the study of the evolution of a system while searching systematically for the minimum of a function derived from some of its physical properties. Using a certain number of mathematical optimization developments devoted to coastal sciences (Isèbe et al. 2014; Isèbe et al. 2008b; Isèbe et al. 2008a; Bouharguane et al. 2010; Mohammadi et al. 2014; Mohammadi et al. 2011; Cook et al. 2021; Mohammadi 2017; Bouharguane et al. 2012; Dupont et al. 2023), we have designed a model that describes the evolution of the sea bottom elevation while taking into account the coupling between morphodynamic and wave processes. This study is based on the assumption that the sea bottom adapts in time to minimize a certain wave-related function. The choice of this function determines the driving force behind the morphological evolution of the seabed. This optimization problem is subjected to a limited number of constraints, allowing for a more accurate description of the morphodynamic evolution.

The purpose of this study is to use the Hadamard (1914) derivative in order to calculate the gradient of any cost function $\mathcal{J}$ with respect to the shape $\psi$, which allows us to solve the optimization problem at the core of the model. This strategy aims to create a generic morphodynamic model that can be used with any wave model.

The paper starts with a description of the OptiMorph model. Then we introduce Hadamard's strategy by presenting the different ways to compute the gradient with respect to the shape $\psi$. Hadamard's strategy is verified with analytical cases. Finally, applications are performed with OptiMorph model using Hadamard strategy. We show that we can therefore use complex wave models such as XBeach (D. J. Roelvink et al. 2009) and SWAN (Booij et al. 1996). Part of the simulations are linked to the LIP and SANDS flume experiments (JA Roelvink et al. 1995; Eichentopf et al. 2018). Another part concerns simulations in open-sea configurations.

# 2 Presentation of a hydro-morphodynamic model by minimization principle

In this section, we introduce the model presented in Cook (2021) and Dupont et al. (2023), which uses the notations in Figure 1.

We consider a coordinate system composed of a horizontal axis $x$ and a vertical axis $z$. We denote $\Omega := [0, x_{\max}]$ the domain of the cross-shore profile of the active coastal zone, where $x = 0$ is a fixed point in deep water where
# 1 Introduction

Morphodynamic models are generally very complex and highly parameterized. They separately solve the physical equations of hydrodynamics and morphodynamics at a very small scale of the order of second in time and of the wave length in space. The OptiMorph model that we presented in Cook (2021) and Dupont et al. (2023) proposes a more global approach based on an optimization principle.

The optimization theory is the study of the evolution of a system while searching systematically for the minimum of a function derived from some of its physical properties. Using a certain number of mathematical optimization developments devoted to coastal sciences (Isèbe et al. 2014; Isèbe et al. 2008b; Isèbe et al. 2008a; Bouharguane et al. 2010; Mohammadi et al. 2014; Mohammadi et al. 2011; Cook et al. 2021; Mohammadi 2017; Bouharguane et al. 2012; Dupont et al. 2023), we have designed a model that describes the evolution of the sea bottom elevation while taking into account the coupling between morphodynamic and wave processes. This study is based on the assumption that the sea bottom adapts in time to minimize a certain wave-related function. The choice of this function determines the driving force behind the morphological evolution of the seabed. This optimization problem is subjected to a limited number of constraints, allowing for a more accurate description of the morphodynamic evolution.

The purpose of this study is to use the Hadamard (1914) derivative in order to calculate the gradient of any cost function $\mathcal{J}$ with respect to the shape $\psi$, which allows us to solve the optimization problem at the core of the model. This strategy aims to create a generic morphodynamic model that can be used with any wave model.

The paper starts with a description of the OptiMorph model. Then we introduce Hadamard's strategy by presenting the different ways to compute the gradient with respect to the shape $\psi$. Hadamard's strategy is verified with analytical cases. Finally, applications are performed with OptiMorph model using Hadamard strategy. We show that we can therefore use complex wave models such as XBeach (D. J. Roelvink et al. 2009) and SWAN (Booij et al. 1996). Part of the simulations are linked to the LIP and SANDS flume experiments (JA Roelvink et al. 1995; Eichentopf et al. 2018). Another part concerns simulations in open-sea configurations.

# 2 Presentation of a hydro-morphodynamic model by minimization principle

In this section, we introduce the model presented in Cook (2021) and Dupont et al. (2023), which uses the notations in Figure 1.

We consider a coordinate system composed of a horizontal axis $x$ and a vertical axis $z$. We denote $\Omega := [0, x_{\max}]$ the domain of the cross-shore profile of the active coastal zone, where $x = 0$ is a fixed point in deep water where
<figure><img src="image_1.png" /><figcaption>**Figure 1:** Diagram of a cross-shore profile in the case of an experimental flume.</figcaption></figure>

no significant change in bottom elevation can occur, and $x_{max}$ is an arbitrary point at the shore beyond the shoreline, as shown by Figure 1. The elevation of the sea bottom is a one-dimensional positive function, defined by: $\psi : \Omega \times [0, T_f] \times \Psi \rightarrow \mathbb{R}^+$ where $[0, T_f]$ is the duration of the simulation (s) and $\Psi$ is the set of physical parameters describing the characteristics of the beach profile. In order to model the evolution over time of $\psi$ and given the assumption that $\psi$ changes over time in response to the wave energy, a description of the surface waves is needed.

## 2.1 Hypotheses

The model is based on the principle that nature seeks to minimize the energy it expends. A cost function $J$ governs the evolution of the seafloor and has been developed according to $E_H$, the total energy of the waves. This hypothesis is inspired by minimal surfaces in nature, as for instance in soap bubbles surfaces (Taylor 1976).

To illustrate this, 1) we run a simulation on XBeach with a LIP 11C profile (JA Roelvink et al. 1995) with the parameters presented in the table A1: an offshore significant water height $H_s = 0.6$ m, a wave period $T_0 = 8$ s and a simulation time of 13h. This forcing is constant over the duration of the experiment. A more detailed description of the experiment can be found below in the 5.1 section. 2) We use LIDAR DynaRev data (Blenkinsopp et al. 2021; Schimmels et al. 2020; Martins et al. 2020) which requires no simulation, since all data (including hydrodynamics) come from measurements. By calculating the wave energy $E_H = \frac{1}{16}\rho_w g H^2$ distribution at the beginning and at the end of the experiment, we obtained the figure 2.
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Diagram of a cross-shore profile in the case of an experimental flume.</figcaption></figure>

no significant change in bottom elevation can occur, and $x_{max}$ is an arbitrary point at the shore beyond the shoreline, as shown by Figure 1. The elevation of the sea bottom is a one-dimensional positive function, defined by: $\psi: \Omega \times [0, T_f] \times \Psi \rightarrow \mathbb{R}^+$ where $[0, T_f]$ is the duration of the simulation (s) and $\Psi$ is the set of physical parameters describing the characteristics of the beach profile. In order to model the evolution over time of $\psi$ and given the assumption that $\psi$ changes over time in response to the wave energy, a description of the surface waves is needed.

## 2.1 Hypotheses

The model is based on the principle that nature seeks to minimize the energy it expends. A cost function $J$ governs the evolution of the seafloor and has been developed according to $E_H$, the total energy of the waves. This hypothesis is inspired by minimal surfaces in nature, as for instance in soap bubbles surfaces (Taylor 1976).

To illustrate this, 1) we run a simulation on XBeach with a LIP 11C profile (JA Roelvink et al. 1995) with the parameters presented in the table A1: an offshore significant water height $H_s = 0.6$ m, a wave period $T_0 = 8$ s and a simulation time of 13h. This forcing is constant over the duration of the experiment. A more detailed description of the experiment can be found below in the 5.1 section. 2) We use LIDAR DynaRev data (Blenkinsopp et al. 2021; Schimmels et al. 2020; Martins et al. 2020) which requires no simulation, since all data (including hydrodynamics) come from measurements. By calculating the wave energy $E_H = \frac{1}{16}\rho_w g H^2$ distribution at the beginning and at the end of the experiment, we obtained the figure 2.
Comparison of Wave Energy $E_H$ between
the beginning and the end of the experience

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

**Figure 2:** 1) LIP 1C experiment with *H* generated by XBeach. 2) DynaRev experiment with *H* measured by LIDAR. a) Bottom profile and averaged water height at the beginning of the experiment (grey), Bottom profile and averaged water height at the end of the experiment (brown). b) Wave energies associated with water heights. The energy is calculated on the black rectangle $\Omega_A$.

The quantity $I = \int_{\Omega_A} \mathcal{E}_H dx$ corresponds to the integral of $\mathcal{E}_H$ on the subdomain $\Omega_A$. Unsurprisingly, in this part of the domain $\Omega_A$ (mainly the breaking zone), the energy of the final profile is lower than the energy of the initial profile (≈ 10%).

Other assumptions assess the behaviour of the sea bottom and originate
from general observations. In particular, the bed-load sediment transport is
controlled by the orbital displacement of water particles (Soulsby 1987); thus
a greater sediment mobility has to be considered in shallower waters. Another
natural observation concerns the slope of the seabed, which cannot be overly
steep without an avalanching process occurring (Reineck et al. 1973). Last, in
an experimental wave flume, the quantity of sand must remain constant over
time, with no inflow or outflow of sand to alter the sand stock.

## 2.2 Wave model

The time evolution of the sea bottom elevation is based on the assumption that the bottom evolves to minimize a certain wave quantity. Thus, a wave model providing a description of the surface wave state is essential. In this study, we use the wave models XBeach (D. J. Roelvink et al. 2009), SWAN (Booij et al. 1996) and an extended shoaling model presented in appendix D.

2.2.1 XBeach

The XBeach model is a process-based model developed by the Delft University
of Technology. It is a two-dimensional, depth-integrated numerical model that
Comparison of Wave Energy $E_H$ between
the beginning and the end of the experience

<figure><img src="image_3.png" /><figcaption>**Figure 2:** 1) LIP 1C experiment with *H* generated by XBeach. 2) DynaRev experiment with *H* measured by LIDAR. a) Bottom profile and averaged water height at the beginning of the experiment (grey), Bottom profile and averaged water height at the end of the experiment (brown). b) Wave energies associated with water heights. The energy is calculated on the black rectangle $\Omega_A$.</figcaption></figure>

The quantity $I = \int_{\Omega_A} \mathcal{E}_H dx$ corresponds to the integral of $\mathcal{E}_H$ on the subdomain $\Omega_A$. Unsurprisingly, in this part of the domain $\Omega_A$ (mainly the breaking zone), the energy of the final profile is lower than the energy of the initial profile (≈ 10%).

Other assumptions assess the behaviour of the sea bottom and originate
from general observations. In particular, the bed-load sediment transport is
controlled by the orbital displacement of water particles (Soulsby 1987); thus
a greater sediment mobility has to be considered in shallower waters. Another
natural observation concerns the slope of the seabed, which cannot be overly
steep without an avalanching process occurring (Reineck et al. 1973). Last, in
an experimental wave flume, the quantity of sand must remain constant over
time, with no inflow or outflow of sand to alter the sand stock.

## 2.2 Wave model

The time evolution of the sea bottom elevation is based on the assumption that the bottom evolves to minimize a certain wave quantity. Thus, a wave model providing a description of the surface wave state is essential. In this study, we use the wave models XBeach (D. J. Roelvink et al. 2009), SWAN (Booij et al. 1996) and an extended shoaling model presented in appendix D.

2.2.1 XBeach

The XBeach model is a process-based model developed by the Delft University
of Technology. It is a two-dimensional, depth-integrated numerical model that
simulates the hydrodynamics, sediment transport, and morphological changes of coastal systems. XBeach is a flexible model that can be used to simulate a variety of coastal processes, including wave breaking, bedload transport, and nearshore morphological changes. The model is based on the principles of conservation of mass, momentum, and energy and uses a finite-difference numerical scheme to solve the governing equations. XBeach has been widely used in coastal studies due to its flexibility and accuracy, and it has been applied to a wide range of coastal systems, including estuaries, beaches, and coastal wetlands. The model can be used as a profile model in 1D (Pender et al. 2013), or as an area model in 2D (McCall et al. 2010), and today, there are three modes in which the hydrodynamics can be resolved in XBeach, being:

* **Stationary** – All wave group variations, and thereby all infragravity motions, are neglected, and only the mean motions are included. This type can be applied for modeling morphological changes under moderate wave conditions;

* **Surfbeat** – This in-stationary, hydrostatic mode, is wave group resolving, and is hence also applicable when one is interested in the swash zone processes;

* **Non-hydrostatic** – The non-linear Shallow-Water equations are solved, and hence individual short wave propagation and transformation is resolved.

In our case, we use the Stationary mode.

### 2.2.2 SWAN

The SWAN model, also developed by the Delft University of Technology, is a spectral numerical model designed to simulate waves evolving in coastal regions, lakes, and estuaries under defined wind, bathymetry, and current conditions. It is based on the Energy Density Balance equation linking the advection term to the source and sink terms. Therefore, the wave energy evolves in both geographic and spectral space and changes its aspect due to the presence of wind at the surface, friction with the bottom, or during the breaking of the waves. The SWAN model is a stable model based on unconditionally stable numerical schemes (implicit schemes). SWAN, in its third version, is in stationary mode (optionally non-stationary) and is formulated in Cartesian or spherical coordinates. The unconditional numerical stability of the SWAN model makes its application more effective in shallow water.

### 2.2.3 Extended Shoaling

This model is presented in appendix D. It was inspired by the former wave model which was based on the linear wave theory (Dean et al. 2004). This simple model determines the significant wave height, noted *H* along the cross-shore profile. Let $\Omega = [0, x_{max}]$ be the domain of the cross-shore profile, where
simulates the hydrodynamics, sediment transport, and morphological changes of coastal systems. XBeach is a flexible model that can be used to simulate a variety of coastal processes, including wave breaking, bedload transport, and nearshore morphological changes. The model is based on the principles of conservation of mass, momentum, and energy and uses a finite-difference numerical scheme to solve the governing equations. XBeach has been widely used in coastal studies due to its flexibility and accuracy, and it has been applied to a wide range of coastal systems, including estuaries, beaches, and coastal wetlands. The model can be used as a profile model in 1D (Pender et al. 2013), or as an area model in 2D (McCall et al. 2010), and today, there are three modes in which the hydrodynamics can be resolved in XBeach, being:

* **Stationary** – All wave group variations, and thereby all infragravity motions, are neglected, and only the mean motions are included. This type can be applied for modeling morphological changes under moderate wave conditions;

* **Surfbeat** – This in-stationary, hydrostatic mode, is wave group resolving, and is hence also applicable when one is interested in the swash zone processes;

* **Non-hydrostatic** – The non-linear Shallow-Water equations are solved, and hence individual short wave propagation and transformation is resolved.

In our case, we use the Stationary mode.

### 2.2.2 SWAN

The SWAN model, also developed by the Delft University of Technology, is a spectral numerical model designed to simulate waves evolving in coastal regions, lakes, and estuaries under defined wind, bathymetry, and current conditions. It is based on the Energy Density Balance equation linking the advection term to the source and sink terms. Therefore, the wave energy evolves in both geographic and spectral space and changes its aspect due to the presence of wind at the surface, friction with the bottom, or during the breaking of the waves. The SWAN model is a stable model based on unconditionally stable numerical schemes (implicit schemes). SWAN, in its third version, is in stationary mode (optionally non-stationary) and is formulated in Cartesian or spherical coordinates. The unconditional numerical stability of the SWAN model makes its application more effective in shallow water.

### 2.2.3 Extended Shoaling

This model is presented in appendix D. It was inspired by the former wave model which was based on the linear wave theory (Dean et al. 2004). This simple model determines the significant wave height, noted *H* along the cross-shore profile. Let $\Omega = [0, x_{max}]$ be the domain of the cross-shore profile, where
$x = 0$ is an arbitrary point in deep water, and $x_{max}$ is an arbitrary point beyond the coastline. The $\Omega$ domain is divided into two disjoint subsets: the $\Omega_S$ shoaling zone and the $\Omega_B$ breaking zone. The wave height $H$ on $\Omega_S$ is based on the shoaling equation (1), where $H_0$ is the deep water wave height and $K_S$ is a shoaling coefficient. This simple model is described as follows:

$$H(x,t) = \begin{cases} H_0(t) K s(x,h) & \text{for } x \in \Omega_S \\ \gamma h(x,t) & \text{for } x \in \Omega_B \end{cases}. \quad (1)$$

For the breaking in the $\Omega_B$ area, the model is essentially based on Munk's breaking criterion $\gamma$ (Munk 1949).

## 2.3 Morphodynamic model based on Wave Energy Optimization

The evolution of the sea bottom is assumed to be driven by the minimization of a cost function $J$ (J.s.m$^{-1}$). Recalling the hypotheses made in section 2.1, the shape of the beach profile is determined by the minimization of the potential energy of shoaling waves, for all $t \in [0, T_f]$:

$$J(\psi, t) = \frac{1}{16} \int_{t-T_{coupl}}^{t} \int_{\Omega} \rho_w g H^2(\psi, x, \tau) dx d\tau, \quad (2)$$

where $H$ denotes the height of the waves over the cross-shore profile (m), $\rho_w$ is water density (kg.m$^{-3}$), and $g$ is the gravitational acceleration (m.s$^{-2}$). $T_{coupl}$ (s) defines the coupling time interval between wave and morphodynamic models so that we have $T_f/T_{coupl}$ iterations between the two models involved. Unlike the previous paper (Dupont et al. 2023), this time the functional is calculated over the entire domain. Indeed, wave energy is transferred to other contributions: bottom friction, current, turbulence, etc. These contributions are shown in (Sous et al. 2020). In other words, minimizing $J$ over $\Omega$ with such models means that bed changes occur to maximize dissipation and transfer of wave energy to current.

In order to describe the evolution of the beach profile, whose initial state is given by $\psi_0$, we assume that the sea bottom elevation $\psi$ defined as a function of $x$. In its effort to minimize $J$, the dynamic of $\psi$ is described by:

$$\begin{cases} \psi_t = \Upsilon \Lambda d \\ \psi(t=0) = \psi_0 \end{cases}, \quad (3)$$

where $\psi_t$ is the time derivative of $\psi$, that is the evolution of the bottom elevation over time (m.s$^{-1}$). $\Upsilon$ is a measure of the sand mobility expressed in m.s.kg$^{-1}$. This parameter is defined on the basis of an Exner-type flow model, and its definition can be found in the appendix B. It has the same functionality as XBeach's morphological factor (J.A. Roelvink 2006) where it is possible to divide simulation times by 18 as performed in (Shafiei et al. 2023; Marchesiello et al. 2022) on the LIP-1B experiment. $\Lambda$ is a local function which represents
$x = 0$ is an arbitrary point in deep water, and $x_{max}$ is an arbitrary point beyond the coastline. The $\Omega$ domain is divided into two disjoint subsets: the $\Omega_S$ shoaling zone and the $\Omega_B$ breaking zone. The wave height $H$ on $\Omega_S$ is based on the shoaling equation (1), where $H_0$ is the deep water wave height and $K_S$ is a shoaling coefficient. This simple model is described as follows:

$$H(x,t) = \begin{cases} H_0(t) K s(x,h) & \text{for } x \in \Omega_S \\ \gamma h(x,t) & \text{for } x \in \Omega_B \end{cases}. \quad (1)$$

For the breaking in the $\Omega_B$ area, the model is essentially based on Munk's breaking criterion $\gamma$ (Munk 1949).

## 2.3 Morphodynamic model based on Wave Energy Optimization

The evolution of the sea bottom is assumed to be driven by the minimization of a cost function $J$ (J.s.m$^{-1}$). Recalling the hypotheses made in section 2.1, the shape of the beach profile is determined by the minimization of the potential energy of shoaling waves, for all $t \in [0, T_f]$:

$$J(\psi, t) = \frac{1}{16} \int_{t-T_{coupl}}^{t} \int_{\Omega} \rho_w g H^2(\psi, x, \tau) dx d\tau, \quad (2)$$

where $H$ denotes the height of the waves over the cross-shore profile (m), $\rho_w$ is water density (kg.m$^{-3}$), and $g$ is the gravitational acceleration (m.s$^{-2}$). $T_{coupl}$ (s) defines the coupling time interval between wave and morphodynamic models so that we have $T_f/T_{coupl}$ iterations between the two models involved. Unlike the previous paper (Dupont et al. 2023), this time the functional is calculated over the entire domain. Indeed, wave energy is transferred to other contributions: bottom friction, current, turbulence, etc. These contributions are shown in (Sous et al. 2020). In other words, minimizing $J$ over $\Omega$ with such models means that bed changes occur to maximize dissipation and transfer of wave energy to current.

In order to describe the evolution of the beach profile, whose initial state is given by $\psi_0$, we assume that the sea bottom elevation $\psi$ defined as a function of $x$. In its effort to minimize $J$, the dynamic of $\psi$ is described by:

$$\begin{cases} \psi_t = \Upsilon \Lambda d \\ \psi(t=0) = \psi_0 \end{cases}, \quad (3)$$

where $\psi_t$ is the time derivative of $\psi$, that is the evolution of the bottom elevation over time (m.s$^{-1}$). $\Upsilon$ is a measure of the sand mobility expressed in m.s.kg$^{-1}$. This parameter is defined on the basis of an Exner-type flow model, and its definition can be found in the appendix B. It has the same functionality as XBeach's morphological factor (J.A. Roelvink 2006) where it is possible to divide simulation times by 18 as performed in (Shafiei et al. 2023; Marchesiello et al. 2022) on the LIP-1B experiment. $\Lambda$ is a local function which represents
the influence of the relative water depth $kh$ on the beach profile dynamics and is defined after the term describing the vertical attenuation of the velocity potential according to linear wave theory (Soulsby 1987):

$$ \Lambda: \begin{array}{rcl} \Omega \times [0, h_0] & \to & \mathbb{R}^+ \\ (x, z) & \mapsto & \frac{\cosh(k(x)(h(x) - (h_0 - z)))}{\cosh(k(x)h(x))}. \end{array} \qquad (4) $$

And $d$ is the direction of the descent (J.s.m$^{-2}$), which indicates the manner in which the sea bottom changes. Under unconstrained circumstances, we have

$$ d = -\nabla_{\psi} \mathcal{J}. \qquad (5) $$

It is therefore important to know how to calculate this quantity; this will be the subject of most of the developments in the sections 3 and 4.

**Remark:** This dynamic described (eq (3)), only modifies the bottom elevation and does not account for lateral displacements. It permits, for instance, the apparition of sedimentary bars but cannot predict their lateral displacements. This will be discussed more thoroughly in the discussion section 6.2.

Constraints are added to the model to incorporate minimal physics required to deliver realistic results. The first constraint concerns the local slope of the bottom. Depending on the composition of the sediment, the bottom slope is bounded by a grain-dependent threshold $M_{\text{slope}}$ (Dean et al. 2004). This is conveyed by the following equation involving the local bottom slope:

$$ \left| \frac{\partial \psi}{\partial x} \right| \le M_{\text{slope}}. \qquad (6) $$

The dimensionless parameter $M_{\text{slope}}$ represents the critical angle of repose of the sediment. This angle is based on observed angles in natural beach environments, which are often between 0.01 and 0.2 (Bascom 1951; Vos et al. 2020; Short 1996). We have considered the observed critical angle of 0.2.

A second example concerns the sand stock in the case of an experimental flume. In a flume, the quantity of sand must be constant over time, as given by (7), contrarily to an open-sea configuration where sand can be transported between the nearshore zone and a domain beyond the closure water depth where sediment is definitely lost for beach morphodynamics (Hattori et al. 1980; Quick 1991). This constraint can be written as :

$$ \int_{\Omega} \psi(t, x)dx = \int_{\Omega} \psi_0(x)dx \quad \forall t \in [0, T_f] \qquad (7) $$

This constraint is necessary for verifying and validating the numerical model with the wave flume experimental data.
the influence of the relative water depth $kh$ on the beach profile dynamics and is defined after the term describing the vertical attenuation of the velocity potential according to linear wave theory (Soulsby 1987):

$$ \Lambda: \begin{array}{rcl} \Omega \times [0, h_0] & \to & \mathbb{R}^+ \\ (x, z) & \mapsto & \frac{\cosh(k(x)(h(x) - (h_0 - z)))}{\cosh(k(x)h(x))}. \end{array} \qquad (4) $$

And $d$ is the direction of the descent (J.s.m$^{-2}$), which indicates the manner in which the sea bottom changes. Under unconstrained circumstances, we have

$$ d = -\nabla_{\psi} \mathcal{J}. \qquad (5) $$

It is therefore important to know how to calculate this quantity; this will be the subject of most of the developments in the sections 3 and 4.

**Remark:** This dynamic described (eq (3)), only modifies the bottom elevation and does not account for lateral displacements. It permits, for instance, the apparition of sedimentary bars but cannot predict their lateral displacements. This will be discussed more thoroughly in the discussion section 6.2.

Constraints are added to the model to incorporate minimal physics required to deliver realistic results. The first constraint concerns the local slope of the bottom. Depending on the composition of the sediment, the bottom slope is bounded by a grain-dependent threshold $M_{\text{slope}}$ (Dean et al. 2004). This is conveyed by the following equation involving the local bottom slope:

$$ \left| \frac{\partial \psi}{\partial x} \right| \le M_{\text{slope}}. \qquad (6) $$

The dimensionless parameter $M_{\text{slope}}$ represents the critical angle of repose of the sediment. This angle is based on observed angles in natural beach environments, which are often between 0.01 and 0.2 (Bascom 1951; Vos et al. 2020; Short 1996). We have considered the observed critical angle of 0.2.

A second example concerns the sand stock in the case of an experimental flume. In a flume, the quantity of sand must be constant over time, as given by (7), contrarily to an open-sea configuration where sand can be transported between the nearshore zone and a domain beyond the closure water depth where sediment is definitely lost for beach morphodynamics (Hattori et al. 1980; Quick 1991). This constraint can be written as :

$$ \int_{\Omega} \psi(t, x)dx = \int_{\Omega} \psi_0(x)dx \quad \forall t \in [0, T_f] \qquad (7) $$

This constraint is necessary for verifying and validating the numerical model with the wave flume experimental data.
### 3 Gradient calculation with respect to the shape
$$ \nabla_{\psi} \mathcal{J} $$

Calculation of $\nabla_{\psi} \mathcal{J}$ is necessary to do shape optimization with descent method equation (5). This quantity is not easy to compute since we do not differentiate on an axis but on a shape $\psi$. $\mathcal{J}$ depends on wave height $H$; it is thus advisable to have a very simple wave model in order to differentiate it easily. We assume at first that $\mathcal{J}$ is of the form $\mathcal{J}(H(\psi(x)))$ involving dependencies with respect to wave quantities $H$. This sensitivity is given by:

$$ 
\begin{aligned}
\nabla_{\psi} \mathcal{J} &= \nabla_H \mathcal{J} \nabla_{\psi} H, \\
&= \nabla_H \left( \frac{1}{16} \rho g H^2 \right) \nabla_{\psi} H, \\
&= \frac{1}{8} \rho g H \nabla_{\psi} H.
\end{aligned}
\quad (8) $$

Calculating $\nabla_{\psi} \mathcal{J}$ reduces to that of $\nabla_{\psi} H$. It can be done analytically using the simple shoaling model described in equation (1) as described later in section 3.1. One can also use a heavy formalism like automatic differentiation (Hascoet et al. 2004; Mohammadi et al. 2011). These strategies are described below and in section 4, where we show how to obtain $\nabla_{\psi} H$ whatever may be our functions $H$ and $\psi$. For example, we can directly calculate $\nabla_{\psi} \mathcal{J}$ only from $\mathcal{J}$ and $\psi$.

#### 3.1 Analytical calculation of $\nabla_{\psi} H$

The analytical method is the most precise (because it gives the exact value) and the fastest in calculation time. To illustrate the purpose, we take equations (1) of $H$ and we differentiate them in the following way:

$$ \nabla_{\psi} H = \begin{cases} H_0(t) \nabla_{\psi} K_s(x, h) & \text{for } x \in \Omega_S \\ \gamma \nabla_{\psi} h(x, t) & \text{for } x \in \Omega_B \end{cases}. \quad (9) $$

The problem is reduced to the calculation of $\nabla_{\psi} K_S(x, t)$ and $\nabla_{\psi} h(x, t)$. The relation $h = h_0 - \psi$ ensures that $\nabla_{\psi} h(x, t) = -1$. Moreover, we have:

$$ K_S = \left[ \tanh(kh) \left( 1 + \frac{2kh}{\sinh(2kh)} \right) \right]^{-1/2} . \quad (10) $$

Let $U(X) = \tanh(X) \left(1 + \frac{2X}{\sinh(2X)}\right)$ and $X = kh$. Introducing $U$ in equation (10) and derivating $\psi$ results in:

$$ \nabla_{\psi} K_S = -\frac{1}{2} U^{-3/2} \nabla_{\psi} U. \quad (11) $$

By trigonometric transformation, we can demonstrate that:
### 3 Gradient calculation with respect to the shape
$$ \nabla_{\psi} \mathcal{J} $$

Calculation of $\nabla_{\psi} \mathcal{J}$ is necessary to do shape optimization with descent method equation (5). This quantity is not easy to compute since we do not differentiate on an axis but on a shape $\psi$. $\mathcal{J}$ depends on wave height $H$; it is thus advisable to have a very simple wave model in order to differentiate it easily. We assume at first that $\mathcal{J}$ is of the form $\mathcal{J}(H(\psi(x)))$ involving dependencies with respect to wave quantities $H$. This sensitivity is given by:

$$ 
\begin{aligned}
\nabla_{\psi} \mathcal{J} &= \nabla_H \mathcal{J} \nabla_{\psi} H, \\
&= \nabla_H \left( \frac{1}{16} \rho g H^2 \right) \nabla_{\psi} H, \\
&= \frac{1}{8} \rho g H \nabla_{\psi} H.
\end{aligned}
\quad (8) $$

Calculating $\nabla_{\psi} \mathcal{J}$ reduces to that of $\nabla_{\psi} H$. It can be done analytically using the simple shoaling model described in equation (1) as described later in section 3.1. One can also use a heavy formalism like automatic differentiation (Hascoet et al. 2004; Mohammadi et al. 2011). These strategies are described below and in section 4, where we show how to obtain $\nabla_{\psi} H$ whatever may be our functions $H$ and $\psi$. For example, we can directly calculate $\nabla_{\psi} \mathcal{J}$ only from $\mathcal{J}$ and $\psi$.

#### 3.1 Analytical calculation of $\nabla_{\psi} H$

The analytical method is the most precise (because it gives the exact value) and the fastest in calculation time. To illustrate the purpose, we take equations (1) of $H$ and we differentiate them in the following way:

$$ \nabla_{\psi} H = \begin{cases} H_0(t) \nabla_{\psi} K_s(x, h) & \text{for } x \in \Omega_S \\ \gamma \nabla_{\psi} h(x, t) & \text{for } x \in \Omega_B \end{cases}. \quad (9) $$

The problem is reduced to the calculation of $\nabla_{\psi} K_S(x, t)$ and $\nabla_{\psi} h(x, t)$. The relation $h = h_0 - \psi$ ensures that $\nabla_{\psi} h(x, t) = -1$. Moreover, we have:

$$ K_S = \left[ \tanh(kh) \left( 1 + \frac{2kh}{\sinh(2kh)} \right) \right]^{-1/2} . \quad (10) $$

Let $U(X) = \tanh(X) \left(1 + \frac{2X}{\sinh(2X)}\right)$ and $X = kh$. Introducing $U$ in equation (10) and derivating $\psi$ results in:

$$ \nabla_{\psi} K_S = -\frac{1}{2} U^{-3/2} \nabla_{\psi} U. \quad (11) $$

By trigonometric transformation, we can demonstrate that:
$$ \nabla_{\psi} U = \nabla_{\psi} X \frac{2 \cosh^2(X) - X \sinh(2X)}{\cosh^4(X)}, \qquad (12) $$

we also have:

$$ \nabla_{\psi} X = h \nabla_{\psi} k + k \nabla_{\psi} h = h \nabla_{\psi} k - k. \qquad (13) $$

Moreover, differentiating both sides of the dispersion equation $\sigma^2 = gk \tanh(kh)$ by $\psi$ gives

$$ \nabla_{\psi} k = \frac{k^2}{\cosh(kh) \sinh(kh) + kh}. \qquad (14) $$

Combining (11), (12), and (14), we obtain $\nabla_{\psi} K_S$, and therefore $\nabla_{\psi} H$, on $\Omega_S$.

This method is the most accurate and robust because it gives the analytical solution directly. However, it is applicable to a very limited number of wave models. Indeed, they must be very simple to be differentiated by hand. Our ambition is to have a strategy that may allow to differentiate the mathematical representation of any hydrodynamic.

## 3.2 Finite difference calculation of $\nabla_{\psi} H$

Finite difference methods are based on the idea of approximating the derivative of a function at a point by taking the difference between the values of the function at two adjacent points. Considering the directional gradient formula at $\psi$ along direction $l$:

$$ \nabla_{\psi} H(\psi) = \lim_{\epsilon \to 0} \frac{H(\psi + \epsilon l) - H(\psi)}{\epsilon}. \qquad (15) $$

We could define a first order finite difference approximation of the gradient at $\psi \in \mathbb{R}^N$ taking for $l$ the vectors $e_{i=1,...,N}$ of the canonical basis of $\mathbb{R}^N$. The $i^{th}$ evaluation provides the corresponding component of the gradient vector. This method requires $N + 1$ evaluation of the wave model which makes the method computationally expensive, as it can be classically of the order of several thousand runs in practice.

## 3.3 Automatic differentiation (AD) method to calculate $\nabla_{\psi} H$

Automatic differentiation (AD) of programs is an important tool for numerical optimization and scientific computing. It is a technique for computing derivatives of a given program by successive derivation of the lines of the code. AD can be used to compute derivatives of functions with respect to both scalar and vector variables (Griewank et al. 2008).
$$ \nabla_{\psi} U = \nabla_{\psi} X \frac{2 \cosh^2(X) - X \sinh(2X)}{\cosh^4(X)}, \qquad (12) $$

we also have:

$$ \nabla_{\psi} X = h \nabla_{\psi} k + k \nabla_{\psi} h = h \nabla_{\psi} k - k. \qquad (13) $$

Moreover, differentiating both sides of the dispersion equation $\sigma^2 = gk \tanh(kh)$ by $\psi$ gives

$$ \nabla_{\psi} k = \frac{k^2}{\cosh(kh) \sinh(kh) + kh}. \qquad (14) $$

Combining (11), (12), and (14), we obtain $\nabla_{\psi} K_S$, and therefore $\nabla_{\psi} H$, on $\Omega_S$.

This method is the most accurate and robust because it gives the analytical solution directly. However, it is applicable to a very limited number of wave models. Indeed, they must be very simple to be differentiated by hand. Our ambition is to have a strategy that may allow to differentiate the mathematical representation of any hydrodynamic.

## 3.2 Finite difference calculation of $\nabla_{\psi} H$

Finite difference methods are based on the idea of approximating the derivative of a function at a point by taking the difference between the values of the function at two adjacent points. Considering the directional gradient formula at $\psi$ along direction $l$:

$$ \nabla_{\psi} H(\psi) = \lim_{\epsilon \to 0} \frac{H(\psi + \epsilon l) - H(\psi)}{\epsilon}. \qquad (15) $$

We could define a first order finite difference approximation of the gradient at $\psi \in \mathbb{R}^N$ taking for $l$ the vectors $e_{i=1,...,N}$ of the canonical basis of $\mathbb{R}^N$. The $i^{th}$ evaluation provides the corresponding component of the gradient vector. This method requires $N + 1$ evaluation of the wave model which makes the method computationally expensive, as it can be classically of the order of several thousand runs in practice.

## 3.3 Automatic differentiation (AD) method to calculate $\nabla_{\psi} H$

Automatic differentiation (AD) of programs is an important tool for numerical optimization and scientific computing. It is a technique for computing derivatives of a given program by successive derivation of the lines of the code. AD can be used to compute derivatives of functions with respect to both scalar and vector variables (Griewank et al. 2008).
### 3.3.1 Direct and reverse modes of Automatic Differentiation

Direct AD uses the chain rule to compute derivatives of a program with respect to the input parameters of the code. The direct AD method can be used to compute derivatives of functions of any order, including higher-order derivatives. This method is relatively simple to implement, and is often used when the number of input variables is small. On the other hand, when the size of input variables is large, the reverse mode of AD is used. The computation cost is independent of the size of the inputs. A typical AD tool is the TAPENADE program (Hascoet et al. 2004) which provides Fortran or C codes for the derivatives of programs in direct and reverse modes. This means that we need to provide the source code. As a consequence, the main limitation of this approach is that it cannot be applied to a commercial code when the source code is not provided. Even when the code is provided (open source), it is written in a modular way, which makes it very difficult to isolate the variables to differentiate.

## 4 Using Hadamard for the calculation of $\nabla_{\psi}\mathcal{J}$

In this section, we focus on the calculation of $\nabla_{\psi}H$ in order to obtain $\nabla_{\psi}\mathcal{J}$ (as illustrated in section 3). This method can be applied to whatever may be the variables: we can directly calculate $\nabla_{\psi}\mathcal{J}$. However, in this case, the approximation would be less good because the analytical derivative of $\mathcal{J}$ (equation (8)) is always more accurate.

### 4.1 Principle

We use the approximation described in (Hadamard 1914; Mohammadi 2007; Mohammadi 2010). We consider $\nabla_{\psi}H$ in the sense of Hadamard following the definition:

$$ \nabla_{\psi}H = \lim_{\epsilon \to 0} \frac{H(\psi + \epsilon n) - H(\psi)}{\epsilon}, \qquad (16) $$

where $n$ is the normal to the shape $\psi$. This can be seen as applying a Gâteaux (1913) derivation in the direction normal to the shape. The principle is illustrated in figure 3.

Using the Taylor-Young formula at order 1, we consider the following approximation:

$$ \begin{aligned} \nabla_{\psi}H &= \lim_{\epsilon \to 0} \frac{H(\psi) + \epsilon \nabla_{X}H.n - H(\psi)}{\epsilon}, \qquad (17) \\ &\approx (\nabla_{X}H).n, \end{aligned} $$

with $X = (x, z)^{\mathrm{T}}$. This approximation is illustrated in the appendix C on simple analytical examples; and also on the simple shoaling model in section 4.3.
### 3.3.1 Direct and reverse modes of Automatic Differentiation

Direct AD uses the chain rule to compute derivatives of a program with respect to the input parameters of the code. The direct AD method can be used to compute derivatives of functions of any order, including higher-order derivatives. This method is relatively simple to implement, and is often used when the number of input variables is small. On the other hand, when the size of input variables is large, the reverse mode of AD is used. The computation cost is independent of the size of the inputs. A typical AD tool is the TAPENADE program (Hascoet et al. 2004) which provides Fortran or C codes for the derivatives of programs in direct and reverse modes. This means that we need to provide the source code. As a consequence, the main limitation of this approach is that it cannot be applied to a commercial code when the source code is not provided. Even when the code is provided (open source), it is written in a modular way, which makes it very difficult to isolate the variables to differentiate.

## 4 Using Hadamard for the calculation of $\nabla_{\psi}\mathcal{J}$

In this section, we focus on the calculation of $\nabla_{\psi}H$ in order to obtain $\nabla_{\psi}\mathcal{J}$ (as illustrated in section 3). This method can be applied to whatever may be the variables: we can directly calculate $\nabla_{\psi}\mathcal{J}$. However, in this case, the approximation would be less good because the analytical derivative of $\mathcal{J}$ (equation (8)) is always more accurate.

### 4.1 Principle

We use the approximation described in (Hadamard 1914; Mohammadi 2007; Mohammadi 2010). We consider $\nabla_{\psi}H$ in the sense of Hadamard following the definition:

$$ \nabla_{\psi}H = \lim_{\epsilon \to 0} \frac{H(\psi + \epsilon n) - H(\psi)}{\epsilon}, \qquad (16) $$

where $n$ is the normal to the shape $\psi$. This can be seen as applying a Gâteaux (1913) derivation in the direction normal to the shape. The principle is illustrated in figure 3.

Using the Taylor-Young formula at order 1, we consider the following approximation:

$$ \begin{aligned} \nabla_{\psi}H &= \lim_{\epsilon \to 0} \frac{H(\psi) + \epsilon \nabla_{X}H.n - H(\psi)}{\epsilon}, \qquad (17) \\ &\approx (\nabla_{X}H).n, \end{aligned} $$

with $X = (x, z)^{\mathrm{T}}$. This approximation is illustrated in the appendix C on simple analytical examples; and also on the simple shoaling model in section 4.3.
<figure><img src="image_3.png" /><figcaption>**Figure 3:** Representation of two sea bottom profiles $\psi$ and $\psi + \varepsilon n$. To calculate the gradient, we need to calculate at all points the associated normal vector $n$.</figcaption></figure>

## 4.2 Numerical validation

The approximation (17) can be verified by calculating numerically the solution of the analytical example presented in the appendix C.2. We calculate the error $L^2$ named $E_{L^2} =||(\nabla_{\psi}H)_{exact} - (\nabla_{\psi}H)_{numerical}||_{L^2}$ for points which correspond to the spatial steps $dx = [10^{-5}, 10^{-4}, 10^{-3}, 10^{-2}, 0.1, 1, 10, 100]$ for a length $L = 1000$ m. We obtain the curves in figure 4.

<figure><img src="image_4.png" /><figcaption>**Figure 4:** A) Calculation of $\nabla_{\psi}H$ using Hadamard approximation with the following problem (see C.2): $\psi = ax + b$, $H = \cos(\psi)$, with $a = 0.02, b = -2$. B) L2 error and order of convergence for a comparison between the analytical solution of the simple problem described in C.2 with $H = \cos(\psi)$.</figcaption></figure>

We notice that the approximation becomes good very quickly. On the figure 4.A, we see that an increment $dx = 20$ m is enough to reach an almost perfect approximation. The figure 4.B shows that the error is very small and converges to the order $O(dx^{1/2})$. The sources of error for this calculation could be a) the approximation in the calculation of the vector $n$ (in this case, it is null because $\psi$ is linear), b) the computation of the gradient by finite differences
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Representation of two sea bottom profiles $\psi$ and $\psi + \varepsilon n$. To calculate the gradient, we need to calculate at all points the associated normal vector $n$.</figcaption></figure>

## 4.2 Numerical validation

The approximation (17) can be verified by calculating numerically the solution of the analytical example presented in the appendix C.2. We calculate the error $L^2$ named $E_{L^2} =||(\nabla_{\psi}H)_{exact} - (\nabla_{\psi}H)_{numerical}||_{L^2}$ for points which correspond to the spatial steps $dx = [10^{-5}, 10^{-4}, 10^{-3}, 10^{-2}, 0.1, 1, 10, 100]$ for a length $L = 1000$ m. We obtain the curves in figure 4.

<figure><img src="image_5.png" /><figcaption>**Figure 4:** A) Calculation of $\nabla_{\psi}H$ using Hadamard approximation with the following problem (see C.2): $\psi = ax + b$, $H = \cos(\psi)$, with $a = 0.02, b = -2$. B) L2 error and order of convergence for a comparison between the analytical solution of the simple problem described in C.2 with $H = \cos(\psi)$.</figcaption></figure>

We notice that the approximation becomes good very quickly. On the figure 4.A, we see that an increment $dx = 20$ m is enough to reach an almost perfect approximation. The figure 4.B shows that the error is very small and converges to the order $O(dx^{1/2})$. The sources of error for this calculation could be a) the approximation in the calculation of the vector $n$ (in this case, it is null because $\psi$ is linear), b) the computation of the gradient by finite differences
(order 1).

## 4.3 Validating the Hadamard solution

Historically, the OptiMorph model used the shoaling equation (1) and was based on the analytical differentiation of this equation (section 3.1). The Hadamard strategy allows us to obtain a calculation of $\nabla_{\psi} H$ in a numerical way, as with finite differences. To implement this approach practically,

we simply need to use the equation (17) with: $\nabla_X H = \begin{pmatrix} \frac{\partial H}{\partial x} \\ \frac{\partial H}{\partial \psi} \end{pmatrix}$ and $n = \frac{1}{\sqrt{d\psi^2+dx^2}} \begin{pmatrix} -d\psi \\ dx \end{pmatrix}$ and we obtain:

$$ \nabla_{\psi} H \approx \frac{\partial H}{\partial x} n_x + \frac{\partial H}{\partial \psi} n_z, \qquad (18) $$

with $n_x$ and $n_z$ the x and z component of $n$. In OptiMorph, we implement equation (18) and we compare the calculations of $\nabla_{\psi} H$ using the simple shoaling model presented in (1). The figure 5 shows a comparison of the Hadamard and exact solutions on a representative example: an offshore water level $H_0 = 2$ m, an offshore water depth $h_0 = 10$ m, a wave period $T_0 = 10$ s and a linear bottom profile $\psi$. The figure 5.A corresponds to a simple case and the figure 5.B to a case with small scales perturbations of the sea bottom.

We notice that the approximation is very good. There is still one point that has a defect in the non-linearity at $x = 670$ m. However, this does not alter the morphodynamic results. To be sure of the robustness, we add non-linearity with a random function that induces perturbations. These are composed of sinusoidal functions and random translations between [-0.2, 0.2], on the entire domain. We also set a hole at $x = 650$ m. We obtain the simulation figure 5.B. Even with all these perturbations, the Hadamard approximation remains very robust.
(order 1).

## 4.3 Validating the Hadamard solution

Historically, the OptiMorph model used the shoaling equation (1) and was based on the analytical differentiation of this equation (section 3.1). The Hadamard strategy allows us to obtain a calculation of $\nabla_{\psi} H$ in a numerical way, as with finite differences. To implement this approach practically,

we simply need to use the equation (17) with: $\nabla_X H = \begin{pmatrix} \frac{\partial H}{\partial x} \\ \frac{\partial H}{\partial \psi} \end{pmatrix}$ and $n = \frac{1}{\sqrt{d\psi^2+dx^2}} \begin{pmatrix} -d\psi \\ dx \end{pmatrix}$ and we obtain:

$$ \nabla_{\psi} H \approx \frac{\partial H}{\partial x} n_x + \frac{\partial H}{\partial \psi} n_z, \qquad (18) $$

with $n_x$ and $n_z$ the x and z component of $n$. In OptiMorph, we implement equation (18) and we compare the calculations of $\nabla_{\psi} H$ using the simple shoaling model presented in (1). The figure 5 shows a comparison of the Hadamard and exact solutions on a representative example: an offshore water level $H_0 = 2$ m, an offshore water depth $h_0 = 10$ m, a wave period $T_0 = 10$ s and a linear bottom profile $\psi$. The figure 5.A corresponds to a simple case and the figure 5.B to a case with small scales perturbations of the sea bottom.

We notice that the approximation is very good. There is still one point that has a defect in the non-linearity at $x = 670$ m. However, this does not alter the morphodynamic results. To be sure of the robustness, we add non-linearity with a random function that induces perturbations. These are composed of sinusoidal functions and random translations between [-0.2, 0.2], on the entire domain. We also set a hole at $x = 650$ m. We obtain the simulation figure 5.B. Even with all these perturbations, the Hadamard approximation remains very robust.
Comparison of numerical and analytical $\nabla_{\psi}H$ on the following problem: $\psi$ linear without / with perturbations - $H_0 = 2$ m - $h_0 = 10$ m - $T_0 = 10$ s

<figure><img src="image_5.png" /><figcaption>**Figure 5:** Comparison of numerical and analytical solution of $\nabla_{\psi}H$ using Opti-Morph model. Configuration without (A) / with (B) perturbations, $H_0 = 2$ m, $h_0 = 10$ m and the wave period $T_0 = 10$ s. In dodgerblue, the wave height $H$, in brown the bottom profile $\psi$, in red $\nabla_{\psi}H$ calculated analytically, in blue $\nabla_{\psi}H$ calculated by Hadamard strategy.</figcaption></figure>

# 5 Application of Hadamard Strategy

To go further, we can use the Hadamard strategy to couple any wave model to the morphodynamic model based on the gradient descent equation presented in (3). The figure 6 shows the detailed implementation of this coupling.

<figure><img src="image_6.png" /><figcaption>**Figure 6:** OptiMorph workflow coupled with wave model</figcaption></figure>

In this section, we perform Hadamard morphodynamics simulations forced by three distinct hydrostatic models: our extended shoaling model presented in appendix D, SWAN and XBeach. Simulations are performed on 5 different experimental data sets: (i) one configuration from the SANDS experience (Eichentopf et al. 2018); (ii) one configuration from the LIP 11D experience
Comparison of numerical and analytical $\nabla_{\psi} H$ on the following problem: $\psi$ linear without / with perturbations - $H_0 = 2$ m - $h_0 = 10$ m - $T_0 = 10$ s

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Comparison of numerical and analytical solution of $\nabla_{\psi} H$ using Opti-Morph model. Configuration without (A) / with (B) perturbations, $H_0 = 2$ m, $h_0 = 10$ m and the wave period $T_0 = 10$ s. In dodgerblue, the wave height $H$, in brown the bottom profile $\psi$, in red $\nabla_{\psi} H$ calculated analytically, in blue $\nabla_{\psi} H$ calculated by Hadamard strategy.</figcaption></figure>

# 5 Application of Hadamard Strategy

To go further, we can use the Hadamard strategy to couple any wave model to the morphodynamic model based on the gradient descent equation presented in (3). The figure 6 shows the detailed implementation of this coupling.

<figure><img src="image_7.png" /><figcaption>**Figure 6:** OptiMorph workflow coupled with wave model</figcaption></figure>

In this section, we perform Hadamard morphodynamics simulations forced by three distinct hydrostatic models: our extended shoaling model presented in appendix D, SWAN and XBeach. Simulations are performed on 5 different experimental data sets: (i) one configuration from the SANDS experience (Eichentopf et al. 2018); (ii) one configuration from the LIP 11D experience
flume experiment presented in table A1 (part of the XBeach benchmark (JA Roelvink et al. 1995)); (iii) three from open-sea configurations with linear, concave and convex bottom profiles.

# 5.1 Description of flumes experiments

In this section, we briefly present the LIP 11D (JA Roelvink et al. 1995) and SANDS (Eichentopf et al. 2018) experiments. These morphodynamic experiments are necessary to validate our model.

## 5.1.1 The SANDS experiments

The experimental setup for this study was conducted at the Canal d'Investigació i Experimentació Marítima (CIEM), a large-scale wave flume located within the Universitat Politècnica de Catalunya (UPC) in Barcelona, Spain. The CIEM is a large-scale wave flume of 100 m length, 3 m width and 4.5 m depth with a working water depth of 2.47 m and 2.5 m. Waves were generated using a hydraulic wave paddle positioned at the end of the deep-water section in the wave flume. The initial beach profile was carefully crafted by hand using well-sorted, commercial sand with a narrow grain size distribution ($d_{50}$ = 0.25 mm, $d_{10}$ = 0.154 mm, $d_{90}$ = 0.372 mm), resulting in a measured sediment fall velocity of $w_s$ = 0.034 m/s. The active portion of the beach profile featured a slope of 1/15. The experimental configuration of the SANDS project in Barcelona was meticulously documented in (Alsina et al. 2011).

This experiment (Eichentopf et al. 2018) is composed of two parts, an erosive part on a linear beach with slopes 1/15 with a forcing of $H_s$ = 0.53 m and $T_0$ = 4.14 s for an experiment duration of 23 hours and 30 minutes. An accretionary section on the final beach profile of the erosive section, with a forcing of $H_s$ = 0.32 m and $T_0$ = 5.44 s for an experiment duration of 20 hours and 25 minutes.

## 5.1.2 The LIP experiments

The Large Installations Plan (LIP) experiments were conducted in the Delta Flume of Delft Hydraulics (now Deltares) (JA Roelvink et al. 1995), which is a large-scale facility measuring 225 x 7 x 5 m. During these experiments, various parameters such as water levels, wave-averaged velocity and suspended concentration profiles, orbital velocities, and bed levels were measured.

Three types of experiments were carried out in LIP under different types of irregular waves, resulting in three distinct beach states: stable (LIP 1A), erosive (LIP 1B), and accretive (LIP 1C).

In LIP 1A, the initial profile was linear with a slope of 1/30 and a median grain size of 0.22 mm. This part of the experiment represented a pre-storm event with the creation of a sedimentary bar under moderate wave conditions ($H_s$ = 0.9 m, $T_0$ = 5 s).
flume experiment presented in table A1 (part of the XBeach benchmark (JA Roelvink et al. 1995)); (iii) three from open-sea configurations with linear, concave and convex bottom profiles.

# 5.1 Description of flumes experiments

In this section, we briefly present the LIP 11D (JA Roelvink et al. 1995) and SANDS (Eichentopf et al. 2018) experiments. These morphodynamic experiments are necessary to validate our model.

## 5.1.1 The SANDS experiments

The experimental setup for this study was conducted at the Canal d'Investigació i Experimentació Marítima (CIEM), a large-scale wave flume located within the Universitat Politècnica de Catalunya (UPC) in Barcelona, Spain. The CIEM is a large-scale wave flume of 100 m length, 3 m width and 4.5 m depth with a working water depth of 2.47 m and 2.5 m. Waves were generated using a hydraulic wave paddle positioned at the end of the deep-water section in the wave flume. The initial beach profile was carefully crafted by hand using well-sorted, commercial sand with a narrow grain size distribution ($d_{50}$ = 0.25 mm, $d_{10}$ = 0.154 mm, $d_{90}$ = 0.372 mm), resulting in a measured sediment fall velocity of $w_s$ = 0.034 m/s. The active portion of the beach profile featured a slope of 1/15. The experimental configuration of the SANDS project in Barcelona was meticulously documented in (Alsina et al. 2011).

This experiment (Eichentopf et al. 2018) is composed of two parts, an erosive part on a linear beach with slopes 1/15 with a forcing of $H_s$ = 0.53 m and $T_0$ = 4.14 s for an experiment duration of 23 hours and 30 minutes. An accretionary section on the final beach profile of the erosive section, with a forcing of $H_s$ = 0.32 m and $T_0$ = 5.44 s for an experiment duration of 20 hours and 25 minutes.

## 5.1.2 The LIP experiments

The Large Installations Plan (LIP) experiments were conducted in the Delta Flume of Delft Hydraulics (now Deltares) (JA Roelvink et al. 1995), which is a large-scale facility measuring 225 x 7 x 5 m. During these experiments, various parameters such as water levels, wave-averaged velocity and suspended concentration profiles, orbital velocities, and bed levels were measured.

Three types of experiments were carried out in LIP under different types of irregular waves, resulting in three distinct beach states: stable (LIP 1A), erosive (LIP 1B), and accretive (LIP 1C).

In LIP 1A, the initial profile was linear with a slope of 1/30 and a median grain size of 0.22 mm. This part of the experiment represented a pre-storm event with the creation of a sedimentary bar under moderate wave conditions ($H_s$ = 0.9 m, $T_0$ = 5 s).
The LIP 1B part of the experiment used the final profile from LIP 1A and represented a storm event with larger waves ($H_s$ = 1.4 m, $T_0$ = 5 s). The bar moves seaward under the action of large waves, highlighting the physical process of erosion.

Finally, the LIP 1C part of the experiment used the final profile from LIP 1B and represented a post-storm event with smaller waves ($H_s$ = 0.6 m, $T_0$ = 8 s). The bar moved back towards the coast asymmetrically, highlighting the physical process of accretion.

The bed profile was measured using a profile follower that used an automated sounding system. The LIP experiments provided valuable insights into the morphodynamic behaviour of sandy beaches under different wave conditions and have been widely used to validate numerical models of beach morphodynamics.

## 5.2 Hydro-Morphodynamic results on Flume experiment

To begin, we perform hydro-morphodynamic simulations with our morphodynamic approach using Hadamard's calculation of $\nabla_{\psi}H$. To highlight the phenomenological aspect of our model, we start by performing simulations on SANDS erosive experience (Eichentopf et al. 2018).

In this case, we set up the models as follows. We set a domain $\Omega$ of 53 m in length with a uniform subdivision of 530 cells. For XBeach and SWAN, the incoming wave boundary condition is provided using a JONSWAP wave spectrum (Hasselmann et al. 1973), with a significant wave height of $H_s$ = 0.53 m and a peak frequency at $f_p$ = 4.14 s⁻¹. For the extended shoaling model (appendix D), we use directly $H_s$ and a wave period $T_0$ = 4.14 s. The breaker model of XBeach uses the D. J. Roelvink (1993) formulation, with a breaker coefficient of $\gamma$ = 0.4, a power $n$ = 15, and a wave dissipation coefficient of 0.5. The breaker model of SWAN is based on the Battjes et al. (1978) breaking parametrization and the extended Shoaling model is simply based on a Munk breaking criterion $\gamma$ = 0.4. The mobility parameter $\Upsilon$ of our morphodynamic model has a value of $5 \times 10^{-3}$ m.s.kg⁻¹. The model is set to run 23.5 h using a coupling time of 42.3 s. We compare the numerical results to those experimental data. The hydro-morphodynamic results are presented in figure 7.A and the differences between the final bottom profile $\psi_f$ and initial bottom profile $\psi_0$ are presented in figure 7.B. The reference is the experimental curve in dark red.
The LIP 1B part of the experiment used the final profile from LIP 1A and represented a storm event with larger waves ($H_s$ = 1.4 m, $T_0$ = 5 s). The bar moves seaward under the action of large waves, highlighting the physical process of erosion.

Finally, the LIP 1C part of the experiment used the final profile from LIP 1B and represented a post-storm event with smaller waves ($H_s$ = 0.6 m, $T_0$ = 8 s). The bar moved back towards the coast asymmetrically, highlighting the physical process of accretion.

The bed profile was measured using a profile follower that used an automated sounding system. The LIP experiments provided valuable insights into the morphodynamic behaviour of sandy beaches under different wave conditions and have been widely used to validate numerical models of beach morphodynamics.

## 5.2 Hydro-Morphodynamic results on Flume experiment

To begin, we perform hydro-morphodynamic simulations with our morphodynamic approach using Hadamard's calculation of $\nabla_{\psi}H$. To highlight the phenomenological aspect of our model, we start by performing simulations on SANDS erosive experience (Eichentopf et al. 2018).

In this case, we set up the models as follows. We set a domain $\Omega$ of 53 m in length with a uniform subdivision of 530 cells. For XBeach and SWAN, the incoming wave boundary condition is provided using a JONSWAP wave spectrum (Hasselmann et al. 1973), with a significant wave height of $H_s$ = 0.53 m and a peak frequency at $f_p$ = 4.14 s⁻¹. For the extended shoaling model (appendix D), we use directly $H_s$ and a wave period $T_0$ = 4.14 s. The breaker model of XBeach uses the D. J. Roelvink (1993) formulation, with a breaker coefficient of $\gamma$ = 0.4, a power $n$ = 15, and a wave dissipation coefficient of 0.5. The breaker model of SWAN is based on the Battjes et al. (1978) breaking parametrization and the extended Shoaling model is simply based on a Munk breaking criterion $\gamma$ = 0.4. The mobility parameter $\Upsilon$ of our morphodynamic model has a value of $5 \times 10^{-3}$ m.s.kg⁻¹. The model is set to run 23.5 h using a coupling time of 42.3 s. We compare the numerical results to those experimental data. The hydro-morphodynamic results are presented in figure 7.A and the differences between the final bottom profile $\psi_f$ and initial bottom profile $\psi_0$ are presented in figure 7.B. The reference is the experimental curve in dark red.
<figure><img src="image_7.png" /><figcaption>**Figure 7:** A) Hydro-Morphodynamic results obtained with OptiMorph model using Hadamard strategy with wave models (Shoaling (green), SWAN (red) and XBeach (blue)). Bottom profile configuration from the SANDS erosive experience. Black: bottom profile, green: *H* and *ψ* from improved shoaling with Hadamard strategy, red: *H* and *ψ* from SWAN with Hadamard strategy, blue: *H* and *ψ* from XBeach with Hadamard strategy, dark red: *ψ* from experience. B) Morphodynamic ecarts of *ψ*<sub>*f*</sub> − *ψ*<sub>*i*</sub> obtained with the Shoaling, SWAN, XBeach models and experiment. Bottom profile configuration from the SANDS channel experiment. In green: morphodynamic differences from shoaling with Hadamard strategy (*ψ*<sub>RMSE</sub> = 11.7 cm), red: morphodynamic differences from SWAN with Hadamard strategy (*ψ*<sub>RMSE</sub> = 12.7 cm), blue: morphodynamic differences from XBeach with Hadamard strategy (*ψ*<sub>RMSE</sub> = 13.5 cm), dark red: morphodynamic differences from the experiment.</figcaption></figure>

In all three simulations and the experiment, a sedimentary bar is created over time and a trough is formed between the sandbar and the shore. These sedimentary bars are positioned below the breaking point of the wave. The sedimentary bars from the simulations have one main hump, whereas in the experiment there are two. In the simulations, the trough rises once the water has touched the shore (x = 1 m), while in the experiment, the trough continues afterwards (up to x = 7 m). The three simulations produce relatively similar results.

The next simulation from LIP - 1C flume experiment (JA Roelvink et al. 1995). In this other case, we set a domain Ω of 180 m in length with a uniform subdivision of 180 cells. For XBeach and SWAN, the incoming wave boundary condition is provided using a JONSWAP wave spectrum (Hasselmann et al. 1973), with a significant wave height of *H*<sub>s</sub> = 0.6 m and a peak frequency at *f*<sub>p</sub> = 8 s<sup>-1</sup>. For the extended shoaling model (appendix D), we use directly *H*<sub>s</sub> and a wave period *T*<sub>0</sub> = 8 s. The breaker model of XBeach, SWAN and
<figure><img src="image_8.png" /><figcaption>**Figure 7:** A) Hydro-Morphodynamic results obtained with OptiMorph model using Hadamard strategy with wave models (Shoaling (green), SWAN (red) and XBeach (blue)). Bottom profile configuration from the SANDS erosive experience. Black: bottom profile, green: *H* and *ψ* from improved shoaling with Hadamard strategy, red: *H* and *ψ* from SWAN with Hadamard strategy, blue: *H* and *ψ* from XBeach with Hadamard strategy, dark red: *ψ* from experience. B) Morphodynamic ecarts of *ψ*<sub>*f*</sub> − *ψ*<sub>*i*</sub> obtained with the Shoaling, SWAN, XBeach models and experiment. Bottom profile configuration from the SANDS channel experiment. In green: morphodynamic differences from shoaling with Hadamard strategy (*ψ*<sub>RMSE</sub> = 11.7 cm), red: morphodynamic differences from SWAN with Hadamard strategy (*ψ*<sub>RMSE</sub> = 12.7 cm), blue: morphodynamic differences from XBeach with Hadamard strategy (*ψ*<sub>RMSE</sub> = 13.5 cm), dark red: morphodynamic differences from the experiment.</figcaption></figure>

In all three simulations and the experiment, a sedimentary bar is created over time and a trough is formed between the sandbar and the shore. These sedimentary bars are positioned below the breaking point of the wave. The sedimentary bars from the simulations have one main hump, whereas in the experiment there are two. In the simulations, the trough rises once the water has touched the shore (x = 1 m), while in the experiment, the trough continues afterwards (up to x = 7 m). The three simulations produce relatively similar results.

The next simulation from LIP - 1C flume experiment (JA Roelvink et al. 1995). In this other case, we set a domain Ω of 180 m in length with a uniform subdivision of 180 cells. For XBeach and SWAN, the incoming wave boundary condition is provided using a JONSWAP wave spectrum (Hasselmann et al. 1973), with a significant wave height of *H*<sub>s</sub> = 0.6 m and a peak frequency at *f*<sub>p</sub> = 8 s<sup>-1</sup>. For the extended shoaling model (appendix D), we use directly *H*<sub>s</sub> and a wave period *T*<sub>0</sub> = 8 s. The breaker model of XBeach, SWAN and
extended Shoaling model are the same as in the previous simulation (still with $\gamma = 0.4$). The mobility parameter $\Upsilon$ of our morphodynamic model has a value of $5 \times 10^{-3}$ m.s.kg$^{-1}$. The model is set to run 13 h using a coupling time of 46.8 s.

We compare the numerical results to those experimental data. The hydro-morphodynamic results are presented in figure 8.A and the differences between the final bottom profile $\psi_f$ and initial bottom profile $\psi_0$ are presented in figure 8.B. The reference is the experimental curve in dark red.

<figure><img src="image_8.png" /><figcaption>**Figure 8:** A) Hydro-Morphodynamic results obtained with OptiMorph model using Hadamard strategy with wave models (Shoaling (green), SWAN (red) and XBeach (blue)). Bottom profile configuration from the LIP 1C channel experiment. Red points: measured HRMS, black: bottom profile, green: $H$ and $\psi$ from improved shoaling with Hadamard strategy, red: $H$ and $\psi$ from SWAN with Hadamard strategy, blue: $H$ and $\psi$ from XBeach with Hadamard strategy, dark red: $\psi$ from experience. B) Morphodynamic ecarts of $\psi_f - \psi_i$ obtained with the Shoaling, SWAN, XBeach models and experiment. Bottom profile configuration from the LIP 1C channel experiment. In green: morphodynamic differences from shoaling with Hadamard strategy, red: morphodynamic differences from SWAN with Hadamard strategy, blue: morphodynamic differences from XBeach with Hadamard strategy, dark red: morphodynamic differences from the experiment.</figcaption></figure>

In this experiment, the outer sedimentary bar is moving towards the coast. None of the simulations reproduces this behaviour: the outer bars remain in the same place ($x = 120$ m). In the experiment, the inner sediment bar grows ($x = 140$ m). Simulations show very similar behaviour. The three simulations produce relatively similar results.
extended Shoaling model are the same as in the previous simulation (still with $\gamma = 0.4$). The mobility parameter $\Upsilon$ of our morphodynamic model has a value of $5 \times 10^{-3}$ m.s.kg$^{-1}$. The model is set to run 13 h using a coupling time of 46.8 s.

We compare the numerical results to those experimental data. The hydro-morphodynamic results are presented in figure 8.A and the differences between the final bottom profile $\psi_f$ and initial bottom profile $\psi_0$ are presented in figure 8.B. The reference is the experimental curve in dark red.

<figure><img src="image_9.png" /><figcaption>**Figure 8:** A) Hydro-Morphodynamic results obtained with OptiMorph model using Hadamard strategy with wave models (Shoaling (green), SWAN (red) and XBeach (blue)). Bottom profile configuration from the LIP 1C channel experiment. Red points: measured HRMS, black: bottom profile, green: $H$ and $\psi$ from improved shoaling with Hadamard strategy, red: $H$ and $\psi$ from SWAN with Hadamard strategy, blue: $H$ and $\psi$ from XBeach with Hadamard strategy, dark red: $\psi$ from experience. B) Morphodynamic ecarts of $\psi_f - \psi_i$ obtained with the Shoaling, SWAN, XBeach models and experiment. Bottom profile configuration from the LIP 1C channel experiment. In green: morphodynamic differences from shoaling with Hadamard strategy, red: morphodynamic differences from SWAN with Hadamard strategy, blue: morphodynamic differences from XBeach with Hadamard strategy, dark red: morphodynamic differences from the experiment.</figcaption></figure>

In this experiment, the outer sedimentary bar is moving towards the coast. None of the simulations reproduces this behaviour: the outer bars remain in the same place ($x = 120$ m). In the experiment, the inner sediment bar grows ($x = 140$ m). Simulations show very similar behaviour. The three simulations produce relatively similar results.
## 5.3 Hydro-Morphodynamic results on Open-sea configurations

In this section, we perform simulations in open-sea configurations. The bottom profiles are linear, concave and convex shapes. These shapes are not directly observable in nature but representative of several typical settings (dissipative, reflexive). Forth, they allow to observe if the morphodynamic model is able to reproduce the phenomenology of sedimentary evolution of sand beaches. For these cases, we perform our morphodynamic model using waves from SWAN, XBeach and extended Shoaling models.

In this configuration, we use the same model settings as before with the exception of a $\Omega$ domain length resized at 1000 m with a uniform subdivision in 1000 cells. The forcing is no longer uniform but represents a storm event of 4 days with a peak wave height $H_s = 2$ m. The waves have a period of $T_0 = 12$ s and the water depth at $x = 0$ is $h_0 = 20$ m. The coupling time is set to 345 s. The results of these simulations are presented in figure 9.

<figure><img src="image_9.png" /><figcaption>**Figure 9:** Evolution of $\psi$ using Hadamard strategy with extended Shoaling (green), SWAN (red) and XBeach (blue) models. Simulation on open-sea configuration with linear, convex and concave configurations. Simulation parameters of $H_0 = 2$ m, $T_0 = 12$ s, $h_0 = 20$ m, $\Omega = 1000$ m.</figcaption></figure>

In all three cases, the simulations produce very similar results. Depending on the angle of the slope, a sedimentary bar is observed more or less far from the shore. For a steep angle (convex beach), the sediment bar is very close to the shore; whereas for a slight angle (concave beach), the bar is further from the shore. All these sedimentary bars are all followed by a trough and
## 5.3 Hydro-Morphodynamic results on Open-sea configurations

In this section, we perform simulations in open-sea configurations. The bottom profiles are linear, concave and convex shapes. These shapes are not directly observable in nature but representative of several typical settings (dissipative, reflexive). Forth, they allow to observe if the morphodynamic model is able to reproduce the phenomenology of sedimentary evolution of sand beaches. For these cases, we perform our morphodynamic model using waves from SWAN, XBeach and extended Shoaling models.

In this configuration, we use the same model settings as before with the exception of a $\Omega$ domain length resized at 1000 m with a uniform subdivision in 1000 cells. The forcing is no longer uniform but represents a storm event of 4 days with a peak wave height $H_s = 2$ m. The waves have a period of $T_0 = 12$ s and the water depth at $x = 0$ is $h_0 = 20$ m. The coupling time is set to 345 s. The results of these simulations are presented in figure 9.

<figure><img src="image_10.png" /><figcaption>**Figure 9:** Evolution of $\psi$ using Hadamard strategy with extended Shoaling (green), SWAN (red) and XBeach (blue) models. Simulation on open-sea configuration with linear, convex and concave configurations. Simulation parameters of $H_0 = 2$ m, $T_0 = 12$ s, $h_0 = 20$ m, $\Omega = 1000$ m.</figcaption></figure>

In all three cases, the simulations produce very similar results. Depending on the angle of the slope, a sedimentary bar is observed more or less far from the shore. For a steep angle (convex beach), the sediment bar is very close to the shore; whereas for a slight angle (concave beach), the bar is further from the shore. All these sedimentary bars are all followed by a trough and
are positioned below the breaking point of the wave. In the case of a convex profile, the sediment bar and the wave breaking, produced using SWAN, are closer to shore than the other two simulations.

# 6 Discussion

## 6.1 Computation time

This section is devoted to the analysis of simulation times of the wave models and our morphodynamic calculation. Table 1 corresponds to the computation times for the LIP 11D - 1C simulations (section 5.2) with Hadamard strategy using the SWAN, XBeach and extended Shoaling models. XBeach was used for providing wave calculation only but this model also calculates the circulation. Therefore, it was necessary to run it over a longer time than that required by morphodynamics to get the right significant wave height $H_s$.

<table><caption>**Table 1:** Computation time with 180 points calculated: LIP - 1C with different wave models. Simulations made with a 2.4 GHz computer using a single core on an Intel Xeon E5-2680 processor.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Hydrodynamic</th><th rowspan="2">Morphodynamic by gradient descent</th></tr><tr><th>Shoaling</th><th>SWAN</th><th>XBeach</th></tr></thead><tbody><tr><th>Simulation with 180 points</th><td></td><td></td><td></td><td></td></tr><tr><th>Computation time for 1 iteration (s)</th><td>0.004</td><td>0.278</td><td>7.372</td><td>0.012</td></tr><tr><th>Total computation time for 1000 iterations (mins)</th><td>0.26</td><td>4.83</td><td>123.06</td><td>0.2</td></tr></tbody></table>

We notice that the calculation time of the Shoaling model is very small (direct calculation in python); it is at least 50 times smaller than that of SWAN and XBeach. XBeach calculation times come from the circulation model, which has the advantage of giving the current $u$ (contrary to SWAN) and could be used for another definition of $J$ functional. The morphodynamic calculation time is very small and negligible compared to the hydrodynamic (except shoaling). By increasing the mesh size to 1000 points (5 times more), we obtain the table 2.

<table><caption>**Table 2:** Computation time with 1000 points calculated with different wave models. Simulations made with a 2.4 GHz computer using a single core on an Intel Xeon E5-2680 processor.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Hydrodynamic</th><th rowspan="2">Morphodynamic by gradient descent</th></tr><tr><th>Shoaling</th><th>SWAN</th><th>XBeach</th></tr></thead><tbody><tr><th>Simulation with 1000 points</th><td></td><td></td><td></td><td></td></tr><tr><th>Computation time for 1 iteration (s)</th><td>0.023</td><td>1.193</td><td>28.738</td><td>0.074</td></tr><tr><th>Total computation for with 1000 iterations (mins)</th><td>1.61</td><td>21.12</td><td>480.2</td><td>1.27</td></tr></tbody></table>

The calculation times are also multiplied by 5. To save computing time, we could use some interpolation strategy between grid of the wave tool and that of the morphodynamic model. This would allow performing wave computations on grids with fewer points with the same final results.
are positioned below the breaking point of the wave. In the case of a convex profile, the sediment bar and the wave breaking, produced using SWAN, are closer to shore than the other two simulations.

# 6 Discussion

## 6.1 Computation time

This section is devoted to the analysis of simulation times of the wave models and our morphodynamic calculation. Table 1 corresponds to the computation times for the LIP 11D - 1C simulations (section 5.2) with Hadamard strategy using the SWAN, XBeach and extended Shoaling models. XBeach was used for providing wave calculation only but this model also calculates the circulation. Therefore, it was necessary to run it over a longer time than that required by morphodynamics to get the right significant wave height $H_s$.

<table><caption>**Table 1:** Computation time with 180 points calculated: LIP - 1C with different wave models. Simulations made with a 2.4 GHz computer using a single core on an Intel Xeon E5-2680 processor.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Hydrodynamic</th><th rowspan="2">Morphodynamic by gradient descent</th></tr><tr><th>Shoaling</th><th>SWAN</th><th>XBeach</th></tr></thead><tbody><tr><th>Simulation with 180 points</th><td></td><td></td><td></td><td></td></tr><tr><th>Computation time for 1 iteration (s)</th><td>0.004</td><td>0.278</td><td>7.372</td><td>0.012</td></tr><tr><th>Total computation time for 1000 iterations (mins)</th><td>0.26</td><td>4.83</td><td>123.06</td><td>0.2</td></tr></tbody></table>

We notice that the calculation time of the Shoaling model is very small (direct calculation in python); it is at least 50 times smaller than that of SWAN and XBeach. XBeach calculation times come from the circulation model, which has the advantage of giving the current $u$ (contrary to SWAN) and could be used for another definition of $J$ functional. The morphodynamic calculation time is very small and negligible compared to the hydrodynamic (except shoaling). By increasing the mesh size to 1000 points (5 times more), we obtain the table 2.

<table><caption>**Table 2:** Computation time with 1000 points calculated with different wave models. Simulations made with a 2.4 GHz computer using a single core on an Intel Xeon E5-2680 processor.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Hydrodynamic</th><th rowspan="2">Morphodynamic by gradient descent</th></tr><tr><th>Shoaling</th><th>SWAN</th><th>XBeach</th></tr></thead><tbody><tr><th>Simulation with 1000 points</th><td></td><td></td><td></td><td></td></tr><tr><th>Computation time for 1 iteration (s)</th><td>0.023</td><td>1.193</td><td>28.738</td><td>0.074</td></tr><tr><th>Total computation for with 1000 iterations (mins)</th><td>1.61</td><td>21.12</td><td>480.2</td><td>1.27</td></tr></tbody></table>

The calculation times are also multiplied by 5. To save computing time, we could use some interpolation strategy between grid of the wave tool and that of the morphodynamic model. This would allow performing wave computations on grids with fewer points with the same final results.
## 6.2 Flume simulation

This section is devoted to the morphodynamic behaviour of our model using the Hadamard strategy on flume configuration (JA Roelvink et al. 1995; Eichentopf et al. 2018). The main question is to check whether the numerical model is capable of reproducing the morphodynamic behaviours measured experimentally.

In the SANDS results shown in figure 7, we can see that a sediment bar is created from a linear beach profile (1/15). Although the simulations do not reproduce the sedimentary bar exactly like the experiment, they show very similar results. The sedimentary bar in the simulations is much shorter (in the sense of $x$) than in the experiment. In the simulations, the pattern of troughs between the sediment bar and the shoreline is very similar to that in the experiment. However, in the experiment, the trough goes beyond the water level. This result cannot be observed in our model (except with a tide) as there is currently no mechanism to model this erosion beyond the water level. This induces errors in our model, which conserves the quantity of sand. This lack of sand could explain why our sedimentary bar is shorter than the experimental one.

In the LIP 1C results shown in figure 8, we notice that two main sandbars are observed. The inner one ($x=140$ m) seems to grow. The outer one ($x=120$ m) moves to the shore. The 3 simulations based on Hadamard strategy succeeded in reproducing the behaviour of the inner bar ($x=140$ m). XBeach model coupled to OptiMorph (blue) overestimates this sandbar and SWAN model coupled to OptiMorph (red) underestimates it. On this bar, there is a consequent loss of energy which induces a strong gradient and allows the bar to grow. However, none of the simulations has succeeded in reproducing the behaviour of the outer bar ($x=120$ m) moving towards the shore.

This is because, as mentioned in section 2.3, the parameterization of the sea bottom $\psi$ and equation (3) describing its dynamics only accounts for vertical variations using the gradient of the functional with respect to the sea bottom shape. Therefore, no lateral translation can be predicted by this model. To be able to account for lateral displacements, we need to introduce transport mechanisms, though, for instance, the following modification of the model:

$$ \left\{ \begin{array}{l} \psi_t = \Upsilon \Lambda d - V \nabla_s \psi \\ \psi(t=0) = \psi_0 \end{array} \right. , \qquad (19) $$

where we have introduced a transport operator in the right-hand side. $\nabla_s \psi$ is the spatial derivative of $\psi$ along the mean slope of the sea bottom and $V$ the velocity along this direction. We show the behaviour of the model using the following expression of $V$ involving, the amplitude of orbital velocity at bottom $U_b$ (Wiberg et al. 2008), and the significant wave height $H$:

$$ V = 0.01 U_b \left( \frac{H}{H_{max}} \right)^p \quad \text{with} \quad U_b = \frac{H \pi}{T_0 \sinh(kh)}. \qquad (20) $$
## 6.2 Flume simulation

This section is devoted to the morphodynamic behaviour of our model using the Hadamard strategy on flume configuration (JA Roelvink et al. 1995; Eichentopf et al. 2018). The main question is to check whether the numerical model is capable of reproducing the morphodynamic behaviours measured experimentally.

In the SANDS results shown in figure 7, we can see that a sediment bar is created from a linear beach profile (1/15). Although the simulations do not reproduce the sedimentary bar exactly like the experiment, they show very similar results. The sedimentary bar in the simulations is much shorter (in the sense of $x$) than in the experiment. In the simulations, the pattern of troughs between the sediment bar and the shoreline is very similar to that in the experiment. However, in the experiment, the trough goes beyond the water level. This result cannot be observed in our model (except with a tide) as there is currently no mechanism to model this erosion beyond the water level. This induces errors in our model, which conserves the quantity of sand. This lack of sand could explain why our sedimentary bar is shorter than the experimental one.

In the LIP 1C results shown in figure 8, we notice that two main sandbars are observed. The inner one ($x=140$ m) seems to grow. The outer one ($x=120$ m) moves to the shore. The 3 simulations based on Hadamard strategy succeeded in reproducing the behaviour of the inner bar ($x=140$ m). XBeach model coupled to OptiMorph (blue) overestimates this sandbar and SWAN model coupled to OptiMorph (red) underestimates it. On this bar, there is a consequent loss of energy which induces a strong gradient and allows the bar to grow. However, none of the simulations has succeeded in reproducing the behaviour of the outer bar ($x=120$ m) moving towards the shore.

This is because, as mentioned in section 2.3, the parameterization of the sea bottom $\psi$ and equation (3) describing its dynamics only accounts for vertical variations using the gradient of the functional with respect to the sea bottom shape. Therefore, no lateral translation can be predicted by this model. To be able to account for lateral displacements, we need to introduce transport mechanisms, though, for instance, the following modification of the model:

$$ \left\{ \begin{array}{l} \psi_t = \Upsilon \Lambda d - V \nabla_s \psi \\ \psi(t=0) = \psi_0 \end{array} \right. , \qquad (19) $$

where we have introduced a transport operator in the right-hand side. $\nabla_s \psi$ is the spatial derivative of $\psi$ along the mean slope of the sea bottom and $V$ the velocity along this direction. We show the behaviour of the model using the following expression of $V$ involving, the amplitude of orbital velocity at bottom $U_b$ (Wiberg et al. 2008), and the significant wave height $H$:

$$ V = 0.01 U_b \left( \frac{H}{H_{max}} \right)^p \quad \text{with} \quad U_b = \frac{H \pi}{T_0 \sinh(kh)}. \qquad (20) $$
The dimensionless morphodynamic factor 0.01 has been chosen in order to
make the ranges of lateral and vertical variations comparable.

Figure 10.A illustrates the behaviour of this model for different *p* values. The best choice appears to be *p* = 1 and as expected, the bar moves back towards the coast (asymmetrically) under smaller waves which shows that transport mechanisms were necessary. Velocity distribution figure 10.B shows that the transport mechanisms are mainly located around the two sandbars.

To go farther, and to make the approach generic, we should express this
velocity *V* using the gradient of the functional as done for the vertical mo-
tion using the gradient of the functional with respect to the shape. But this
development is by far beyond the scope of this paper.

<figure><img src="image_10.png" /><figcaption>**Figure 10:** A) Morphodynamic results by the OptiMorph model augmented by the transport mechanisms for *p* = 0, 1, 2 and the XBeach wave model, for the LIP 1C channel experiment. B) Velocity distribution for *p* = 1.</figcaption></figure>

6.3 Open-Sea simulation

This section is devoted to the morphodynamic behaviour of our model using
the Hadamard strategy on open-sea configuration. The two simulations figures
7 and 9, show that there is a creation of the sandbar at the wave breaking point.
Figure 9 shows that a slight pit is created before the sandbar and a trough
one after. These observations are providing because they represent the major
morphologic features along a typical sand bar profiles. Indeed, it is common to
observe sedimentary bars at the wave breaking point. Moreover, the steeper
the slope (convex), the later the breaking, the closer the sandbar is to the
shore (x = 950 m). Conversely, the gentler the slope (concave), the farther
the breaking point, the farther the sandbar is from shore (x = 650 m). These
types of beach profiles are usually observed in nature (Wright et al. 1984).
It highlights the fact that even with an unrealistic initial beach profiles, the
model can produce a realistic beach profile without any need in pre-nucleation
of the bottom perturbation. Forth, whatever the physics behind the waves,
the model is able to produce very similar morphodynamic results.
The dimensionless morphodynamic factor 0.01 has been chosen in order to make the ranges of lateral and vertical variations comparable.

Figure 10.A illustrates the behaviour of this model for different $p$ values. The best choice appears to be $p = 1$ and as expected, the bar moves back towards the coast (asymmetrically) under smaller waves which shows that transport mechanisms were necessary. Velocity distribution figure 10.B shows that the transport mechanisms are mainly located around the two sandbars.

To go farther, and to make the approach generic, we should express this velocity $V$ using the gradient of the functional as done for the vertical motion using the gradient of the functional with respect to the shape. But this development is by far beyond the scope of this paper.

<figure><img src="image_11.png" /><figcaption>**Figure 10:** A) Morphodynamic results by the OptiMorph model augmented by the transport mechanisms for $p = 0, 1, 2$ and the XBeach wave model, for the LIP 1C channel experiment. B) Velocity distribution for $p = 1$.</figcaption></figure>

## 6.3 Open-Sea simulation

This section is devoted to the morphodynamic behaviour of our model using the Hadamard strategy on open-sea configuration. The two simulations figures 7 and 9, show that there is a creation of the sandbar at the wave breaking point. Figure 9 shows that a slight pit is created before the sandbar and a trough one after. These observations are providing because they represent the major morphologic features along a typical sand bar profiles. Indeed, it is common to observe sedimentary bars at the wave breaking point. Moreover, the steeper the slope (convex), the later the breaking, the closer the sandbar is to the shore ($x = 950$ m). Conversely, the gentler the slope (concave), the farther the breaking point, the farther the sandbar is from shore ($x = 650$ m). These types of beach profiles are usually observed in nature (Wright et al. 1984). It highlights the fact that even with an unrealistic initial beach profiles, the model can produce a realistic beach profile without any need in pre-nucleation of the bottom perturbation. Forth, whatever the physics behind the waves, the model is able to produce very similar morphodynamic results.
## 6.4 Gamma sensibility

To highlight the creation of sandbars at the wave breaking point, we artificially change the breaking point by varying the Munk (1949) criterion $\gamma$ on hydro-morphodynamic simulations using Hadamard strategy. These simulations are performed with the SWAN wave model and the same wave parameters as the simulation 5.3 ($T_0$ = 12 s and $H_0$ = 2 m). By taking the Munk (1949) criteria at the values 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, we obtain the figure 11.

<figure><img src="image_11.png" /><figcaption>**Figure 11:** Hydro-morphodynamic results with different breaking criterion $\gamma$ - Simulation with OptiMorph (Hadamard strategy) using SWAN model - $H_0$ = 2 m, $T_0$ = 12 s, $h_0$ = 20 m.</figcaption></figure>

The figure 11 shows that the sandbars are formed systematically at the wave breaking point. The higher the coefficient $\gamma$, the closer to the coast the waves break and the closer the sandbar is to the coast. Moreover, troughs at the lee side of the sandbar like observed in the nature are systematically nucleated (Wright et al. 1984).

## 7 Conclusion

In this paper, we recall what is OptiMorph model and we have highlighted its limitations. Using Hadamard's approach, we calculate an approximation of the gradient $\nabla_{\psi} \mathcal{J}$ of the functional $\mathcal{J}$ with respect to the shape $\psi$ without any additional wave calculation. This study allows us to differentiate any functional $\mathcal{J}$ according to any input variable and relating to any wave model. The analytical and numerical comparisons performed prove that Hadamard strategy is accurate and robust. We applied this tool to realistic and idealized hydro-morphodynamic simulations. The morphodynamic results with
## 6.4 Gamma sensibility

To highlight the creation of sandbars at the wave breaking point, we artificially change the breaking point by varying the Munk (1949) criterion $\gamma$ on hydro-morphodynamic simulations using Hadamard strategy. These simulations are performed with the SWAN wave model and the same wave parameters as the simulation 5.3 ($T_0$ = 12 s and $H_0$ = 2 m). By taking the Munk (1949) criteria at the values 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, we obtain the figure 11.

<figure><img src="image_12.png" /><figcaption>**Figure 11:** Hydro-morphodynamic results with different breaking criterion $\gamma$ - Simulation with OptiMorph (Hadamard strategy) using SWAN model - $H_0$ = 2 m, $T_0$ = 12 s, $h_0$ = 20 m.</figcaption></figure>

The figure 11 shows that the sandbars are formed systematically at the wave breaking point. The higher the coefficient $\gamma$, the closer to the coast the waves break and the closer the sandbar is to the coast. Moreover, troughs at the lee side of the sandbar like observed in the nature are systematically nucleated (Wright et al. 1984).

## 7 Conclusion

In this paper, we recall what is OptiMorph model and we have highlighted its limitations. Using Hadamard's approach, we calculate an approximation of the gradient $\nabla_{\psi} \mathcal{J}$ of the functional $\mathcal{J}$ with respect to the shape $\psi$ without any additional wave calculation. This study allows us to differentiate any functional $\mathcal{J}$ according to any input variable and relating to any wave model. The analytical and numerical comparisons performed prove that Hadamard strategy is accurate and robust. We applied this tool to realistic and idealized hydro-morphodynamic simulations. The morphodynamic results with
SANDS (Eichentopf et al. 2018) are very encouraging because they succeed in reproducing the dynamics of the dominant sedimentary bar. However, the initial results on LIP 1C (JA Roelvink et al. 1995) failed to reproduce the displacement of the outer sedimentary bar, although the behaviour of the inner was well reproduced. The model still needs to be improved in order to fix the lateral displacement and erosion above the water level. Nevertheless, our model is of low-complexity and reproduces the phenomenology as shown by the open-sea and SANDS results where it creates a bar at the breaking point without the need of a priori nucleation or pre-location of the bar.

# 8 Declarations

## 8.1 Availability of data and material

All data, models, and code generated or used during the study appear in the submitted article.

## 8.2 Conflict of interest

The authors declare that they have no conflict of interest.

## 8.3 Acknowledgements

This work was conducted as part as M. Dupont's PhD studies which is funded by the CNRS with the MITI grant. We gratefully acknowledge funding from CNRS, OPTIBEACH projects and FEDER Europe. We also thank to GLADYS (www.gladys-littoral.org) for their continuous logistical and financial support of academic research and applications on coastal hydrodynamics in the South of France. We would also like to thank the CROCO team for providing us with the LIP data.

# Appendix

## A LIP data

<table><caption>Table A1: LIP flume experiment parameters</caption><thead><tr><td>Experiment</td><td>Initial Geometry</td><td>H<sub>S</sub> [m]</td><td>T<sub>p</sub> [s]</td><td>Duration [h]</td></tr></thead><tbody><tr><td>LIP-1A</td><td>Initial beach profile</td><td>0.9</td><td>5</td><td></td></tr><tr><td>LIP-1B</td><td>Result of 1A</td><td>1.4</td><td>5</td><td>18</td></tr><tr><td>LIP-1C</td><td>Result of 1B</td><td>0.6</td><td>8</td><td>13</td></tr></tbody></table>
SANDS (Eichentopf et al. 2018) are very encouraging because they succeed in reproducing the dynamics of the dominant sedimentary bar. However, the initial results on LIP 1C (JA Roelvink et al. 1995) failed to reproduce the displacement of the outer sedimentary bar, although the behaviour of the inner was well reproduced. The model still needs to be improved in order to fix the lateral displacement and erosion above the water level. Nevertheless, our model is of low-complexity and reproduces the phenomenology as shown by the open-sea and SANDS results where it creates a bar at the breaking point without the need of a priori nucleation or pre-location of the bar.

# 8 Declarations

## 8.1 Availability of data and material

All data, models, and code generated or used during the study appear in the submitted article.

## 8.2 Conflict of interest

The authors declare that they have no conflict of interest.

## 8.3 Acknowledgements

This work was conducted as part as M. Dupont's PhD studies which is funded by the CNRS with the MITI grant. We gratefully acknowledge funding from CNRS, OPTIBEACH projects and FEDER Europe. We also thank to GLADYS (www.gladys-littoral.org) for their continuous logistical and financial support of academic research and applications on coastal hydrodynamics in the South of France. We would also like to thank the CROCO team for providing us with the LIP data.

# Appendix

## A LIP data

<table><caption>Table A1: LIP flume experiment parameters</caption><thead><tr><td>Experiment</td><td>Initial Geometry</td><td>H<sub>S</sub> [m]</td><td>T<sub>p</sub> [s]</td><td>Duration [h]</td></tr></thead><tbody><tr><td>LIP-1A</td><td>Initial beach profile</td><td>0.9</td><td>5</td><td></td></tr><tr><td>LIP-1B</td><td>Result of 1A</td><td>1.4</td><td>5</td><td>18</td></tr><tr><td>LIP-1C</td><td>Result of 1B</td><td>0.6</td><td>8</td><td>13</td></tr></tbody></table>
## B Link with morphodynamic flux-based models and sediment characteristics Υ

In this section, we show how to link the bed receptivity coefficient in minimization-based to the bed porosity in classical flux-based morphodynamic models. The literature on morphodynamic models is vast (Nielsen 1992; Nielsen 2002; Rooijen et al. 2012; Chen et al. 2023). Modern numerical implementations rely on models which are in a divergence form. For instance, the Exner equation (Paola et al. 2005; Yang et al. 1996) describes the conservation of mass between sediment in the bed of a channel and sediment that is being transported. It states that bed elevation increases (the bed aggregates) proportionally to the amount of sediment that drops out of transport, and conversely decreases (the bed degrades) proportionally to the amount of sediment that becomes entrained by the flow. The model involves the local porosity of the bed $\lambda_p(x) \in [0, 1]$, a function in space $x$, but not in time. The model writes:

$$ \psi_t + \frac{1}{1 - \lambda_p(x)} \operatorname{div}(q(x, t)) = 0, $$

completed with initial and boundary conditions.

In the literature, the expression of $q$ is diverse. But our discussion remains the same whatever may be the formulation of $q$. For the sake of simplicity, we consider $q_x$ a flux in one dimension of space. Now, let us write the flux-based model and link it to our approach presented through the steepest descent formulation for simplicity:

$$ \psi_t = - \frac{1}{1 - \lambda_p(x)} q_x = -\Upsilon(x) \Lambda(x) \nabla_{\psi} \mathcal{J}(x, t). \quad (B1) $$

There is no explicit boundary condition in the minimization model. In this case, we consider $\Lambda(x) = 1$, the maximum disturbance. As we saw, global sand conservation, can be evaluated through a constraint. In the same way, the local maximum slope is expressed as a constraint. $\nabla_{\psi} \mathcal{J}(x, t)$ corresponds to the direction of the descent without constraint and $d$ with. The bed receptivity $\Upsilon(x)$ is a positive function which we link to the couple bed porosity $\lambda_p(x)$ and flux $q$ as follows.

Locally integrating in space equation (B1) over a small interval $[x - \varepsilon, x + \varepsilon]$ around $x$ we have:

$$ \int_{x-\varepsilon}^{x+\varepsilon} \Upsilon(s) \nabla_{\psi} \mathcal{J}(s,t) ds = \int_{x-\varepsilon}^{x+\varepsilon} \frac{1}{1-\lambda_p(s)} q_s(s,t) ds. $$

Assuming $\Upsilon$ and $\lambda_p$ constant over this small interval, which is physically realistic, we have:

$$ \Upsilon(x) \int_{x-\varepsilon}^{x+\varepsilon} \nabla_{\psi} \mathcal{J}(s,t) ds = \frac{1}{1-\lambda_p(x)} \int_{x-\varepsilon}^{x+\varepsilon} q_s(s,t) ds. $$
## B Link with morphodynamic flux-based models and sediment characteristics $\Upsilon$

In this section, we show how to link the bed receptivity coefficient in minimization-based to the bed porosity in classical flux-based morphodynamic models. The literature on morphodynamic models is vast (Nielsen 1992; Nielsen 2002; Rooijen et al. 2012; Chen et al. 2023). Modern numerical implementations rely on models which are in a divergence form. For instance, the Exner equation (Paola et al. 2005; Yang et al. 1996) describes the conservation of mass between sediment in the bed of a channel and sediment that is being transported. It states that bed elevation increases (the bed aggregates) proportionally to the amount of sediment that drops out of transport, and conversely decreases (the bed degrades) proportionally to the amount of sediment that becomes entrained by the flow. The model involves the local porosity of the bed $\lambda_p(x) \in [0, 1]$, a function in space $x$, but not in time. The model writes:

$$ \psi_t + \frac{1}{1 - \lambda_p(x)} \operatorname{div}(q(x, t)) = 0, $$

completed with initial and boundary conditions.

In the literature, the expression of $q$ is diverse. But our discussion remains the same whatever may be the formulation of $q$. For the sake of simplicity, we consider $q_x$ a flux in one dimension of space. Now, let us write the flux-based model and link it to our approach presented through the steepest descent formulation for simplicity:

$$ \psi_t = - \frac{1}{1 - \lambda_p(x)} q_x = -\Upsilon(x) \Lambda(x) \nabla_{\psi} \mathcal{J}(x, t). \quad (B1) $$

There is no explicit boundary condition in the minimization model. In this case, we consider $\Lambda(x) = 1$, the maximum disturbance. As we saw, global sand conservation, can be evaluated through a constraint. In the same way, the local maximum slope is expressed as a constraint. $\nabla_{\psi} \mathcal{J}(x, t)$ corresponds to the direction of the descent without constraint and $d$ with. The bed receptivity $\Upsilon(x)$ is a positive function which we link to the couple bed porosity $\lambda_p(x)$ and flux $q$ as follows.

Locally integrating in space equation (B1) over a small interval $[x - \varepsilon, x + \varepsilon]$ around $x$ we have:

$$ \int_{x-\varepsilon}^{x+\varepsilon} \Upsilon(s) \nabla_{\psi} \mathcal{J}(s,t) ds = \int_{x-\varepsilon}^{x+\varepsilon} \frac{1}{1-\lambda_p(s)} q_s(s,t) ds. $$

Assuming $\Upsilon$ and $\lambda_p$ constant over this small interval, which is physically realistic, we have:

$$ \Upsilon(x) \int_{x-\varepsilon}^{x+\varepsilon} \nabla_{\psi} \mathcal{J}(s,t) ds = \frac{1}{1-\lambda_p(x)} \int_{x-\varepsilon}^{x+\varepsilon} q_s(s,t) ds. $$
This leads to:

$$ \Upsilon(x) \int_{x-\varepsilon}^{x+\varepsilon} \nabla_{\psi} \mathcal{J}(s,t) ds = \frac{1}{1 - \lambda_p(x)} (q(x+\varepsilon, t) - q(x-\varepsilon, t)) $$

which we write as:

$$ \Upsilon(x) = F(x,t) \frac{1}{1 - \lambda_p(x)}, $$

where factor $F(x, t)$:

$$ F(x,t) = \frac{q(x + \varepsilon, t) - q(x - \varepsilon, t)}{2\varepsilon \overline{\nabla_{\psi} \mathcal{J}}|_{(x,t)}} $$

represents the ratio between the local flux difference and the local average shape gradient $\overline{\nabla_{\psi} \mathcal{J}}|_{(x,t)} = (1/(2\varepsilon)) \int_{x-\varepsilon}^{x+\varepsilon} \nabla_{\psi} \mathcal{J}(s,t) ds$ at point $x$.

If the bed porosity does not change in time, this evaluation is made only once at $t = 0$ and hence, given a flux and a bed porosity, the corresponding minimization-based procedure can receive an equivalent pointwise initialization (at the first iteration in an iterative time integration procedure).

In operational conditions, however, it is very unlikely to have a pointwise, even inaccurate, estimation of $\lambda_p(x)$. It is more reasonable to look for an 'equivalent' constant bed porosity for a given site knowing that what is important in coastal engineering is not the knowledge of the pointwise bed porosity, but the prediction of beach future behaviour based on this site macroscopic characteristics. Also, in homogeneous bed, as it is often the case in sandy beaches, $\lambda_p(x)$ is a constant. We therefore look for a constant bed receptivity $\Upsilon = \overline{F} \frac{1}{1-\lambda_p}$ over the domain of interest $]x_L, x_R[$ (L, R indicating Left and Right) given constant bed porosity $\lambda_p$ and flux $q$ with

$$ \overline{F} = \frac{q(x_R, 0) - q(x_L, 0)}{\int_{x_L}^{x_R} \nabla_{\psi} \mathcal{J}(s, 0) ds}, \quad (B2) $$

which is a scalar, and the ratio between flux variation over $]x_L, x_R[$ and the average of local shape gradients. Here we have defined $x_L$ and $x_R$ as the Left and Right extremities of the domain. So we have $q(x_R, 0)$ and $q(x_L, 0)$ the boundary conditions of the flux-based model. $\overline{F}$ is a measure of how the evaluation of local-based and optimization-based fluxes differs.

## B.1 Illustration using a simple model

Assuming that we are on a configuration of bed load transportation without suspended transport, we can calculate $q(x_R, 0)$ and $q(x_L, 0)$ by using a formula of the bed load transport rate $q$ with (Fredsøe et al. 1992):

$$ q = 10 \frac{\pi}{6} d_{50} p U'_{f} [1 - 0.7 \sqrt{\theta_c / \theta'}] \quad (B3) $$
This leads to:

$$ \Upsilon(x) \int_{x-\varepsilon}^{x+\varepsilon} \nabla_{\psi} \mathcal{J}(s,t) ds = \frac{1}{1 - \lambda_p(x)} (q(x+\varepsilon, t) - q(x-\varepsilon, t)) $$

which we write as:

$$ \Upsilon(x) = F(x,t) \frac{1}{1 - \lambda_p(x)}, $$

where factor $F(x, t)$:

$$ F(x,t) = \frac{q(x+\varepsilon,t) - q(x-\varepsilon,t)}{2\varepsilon \overline{\nabla_{\psi} \mathcal{J}}|_{(x,t)}} $$

represents the ratio between the local flux difference and the local average shape gradient $\overline{\nabla_{\psi} \mathcal{J}}|_{(x,t)} = (1/(2\varepsilon)) \int_{x-\varepsilon}^{x+\varepsilon} \nabla_{\psi} \mathcal{J}(s,t) ds$ at point $x$.

If the bed porosity does not change in time, this evaluation is made only once at $t = 0$ and hence, given a flux and a bed porosity, the corresponding minimization-based procedure can receive an equivalent pointwise initialization (at the first iteration in an iterative time integration procedure).

In operational conditions, however, it is very unlikely to have a pointwise, even inaccurate, estimation of $\lambda_p(x)$. It is more reasonable to look for an 'equivalent' constant bed porosity for a given site knowing that what is important in coastal engineering is not the knowledge of the pointwise bed porosity, but the prediction of beach future behaviour based on this site macroscopic characteristics. Also, in homogeneous bed, as it is often the case in sandy beaches, $\lambda_p(x)$ is a constant. We therefore look for a constant bed receptivity $\Upsilon = \overline{F} \frac{1}{1-\lambda_p}$ over the domain of interest $]x_L, x_R[$ (L, R indicating Left and Right) given constant bed porosity $\lambda_p$ and flux $q$ with

$$ \overline{F} = \frac{q(x_R, 0) - q(x_L, 0)}{\int_{x_L}^{x_R} \nabla_{\psi} \mathcal{J}(s, 0) ds}, \qquad (B2) $$

which is a scalar, and the ratio between flux variation over $]x_L, x_R[$ and the average of local shape gradients. Here we have defined $x_L$ and $x_R$ as the Left and Right extremities of the domain. So we have $q(x_R, 0)$ and $q(x_L, 0)$ the boundary conditions of the flux-based model. $\overline{F}$ is a measure of how the evaluation of local-based and optimization-based fluxes differs.

## B.1 Illustration using a simple model

Assuming that we are on a configuration of bed load transportation without suspended transport, we can calculate $q(x_R, 0)$ and $q(x_L, 0)$ by using a formula of the bed load transport rate $q$ with (Fredsøe et al. 1992):

$$ q = 10 \frac{\pi}{6} d_{50} p U'_{f} [1 - 0.7 \sqrt{\theta_c / \theta'}] \qquad (B3) $$
with $d_{50}$ the grain diameter, $p$ the fraction of bed surface particles in motion, $U'_f$ the skin friction velocity, $\theta_c$ the critical Shields parameter and $\theta'$ the Shields parameter. This formula has been chosen as one of the simplest. However, we can choose to take suspended sediment transport into account, simply by changing the expression of $q$ in our model. Combining the equations (B2) and (B3), we obtain the following expression of $\bar{F}$:

$$ \bar{F} = \pi d_{50} [10 - 7 \sqrt{\theta_c/\theta'}] \frac{p(x_L)U f'(x_L) - p(x_R)U f'(x_R)}{6 \int_{x_L}^{x_R} \nabla_{\psi}(s, 0) ds} $$

We have shown how a conjunct giving of a bed porosity and a flux permits the initialization of a minimization model according to the parameters of the chosen local flux-based model which is comforting for users familiar with such a more traditional approach. However, a same initialization does not mean that the two models will follow the same path, as the minimization-based approach introduces more physics. Indeed, in previous works, we have already shown how our minimization-based formulation can be seen as an Exner equation with a non-local flux (Mohammadi et al. 2011; Bouharguane et al. 2012) with terms similar to those encountered in Fowler-like models (Fowler 2001; Kouakou et al. 2006). Those terms bring the contribution of some non-local physics to the morphodynamics.

This formulation also permits the comparison of the bed $\psi$ evolution predicted minimizing different physical functionals $\mathcal{J}$. It is thus a very efficient exploratory model as defined by MURRAY (Murray 2007). However, it is not possible to find the functional $J$ associated to a given flux $q$ because this requires the mathematical concept of integration with respect to the shape to give sense to:

$$ \mathcal{J} = \frac{\Upsilon}{1 - \lambda_p} \int_{\psi} \nabla \cdot q \, d\psi. $$

Unfortunately, unlike differentiation with respect to the shape (Mohammadi et al. 2009), the concept of integration with respect to the shape does not exist as of today.

# C Analytical examples of Hadamard derivative

In this section, we illustrate analytical examples of derivation of $\psi$ on a quantity $A$; concretely, we calculate $\nabla_{\psi} A$.

## C.1 Flat form

We consider the relation $A = \psi^2$. We set in a general way $\psi = \{(x, y) \in \mathbb{R}^2 | y - f_{\psi}(x) = 0\}$ the space of $\psi$ with $f_{\psi}$ the function describing the bottom. The flat form $\psi = \{(x, y) \in \mathbb{R}^2 | y - c = 0\}$ deformed from $\varepsilon n$ is given by $\psi + \varepsilon n = \{(x, y) \in \mathbb{R}^2 | y - c - \varepsilon = 0\}$. It could be illustrated by the figure C1.
with $d_{50}$ the grain diameter, $p$ the fraction of bed surface particles in motion, $U'_f$ the skin friction velocity, $\theta_c$ the critical Shields parameter and $\theta'$ the Shields parameter. This formula has been chosen as one of the simplest. However, we can choose to take suspended sediment transport into account, simply by changing the expression of $q$ in our model. Combining the equations (B2) and (B3), we obtain the following expression of $\bar{F}$:

$$ \bar{F} = \pi d_{50} [10 - 7 \sqrt{\theta_c/\theta'}] \frac{p(x_L)U f'(x_L) - p(x_R)U f'(x_R)}{6 \int_{x_L}^{x_R} \nabla_{\psi}(s, 0) ds} $$

We have shown how a conjunct giving of a bed porosity and a flux permits the initialization of a minimization model according to the parameters of the chosen local flux-based model which is comforting for users familiar with such a more traditional approach. However, a same initialization does not mean that the two models will follow the same path, as the minimization-based approach introduces more physics. Indeed, in previous works, we have already shown how our minimization-based formulation can be seen as an Exner equation with a non-local flux (Mohammadi et al. 2011; Bouharguane et al. 2012) with terms similar to those encountered in Fowler-like models (Fowler 2001; Kouakou et al. 2006). Those terms bring the contribution of some non-local physics to the morphodynamics.

This formulation also permits the comparison of the bed $\psi$ evolution predicted minimizing different physical functionals $\mathcal{J}$. It is thus a very efficient exploratory model as defined by MURRAY (Murray 2007). However, it is not possible to find the functional $J$ associated to a given flux $q$ because this requires the mathematical concept of integration with respect to the shape to give sense to:

$$ \mathcal{J} = \frac{\Upsilon}{1 - \lambda_p} \int_{\psi} \nabla \cdot q \, d\psi. $$

Unfortunately, unlike differentiation with respect to the shape (Mohammadi et al. 2009), the concept of integration with respect to the shape does not exist as of today.

# C Analytical examples of Hadamard derivative

In this section, we illustrate analytical examples of derivation of $\psi$ on a quantity $A$; concretely, we calculate $\nabla_{\psi} A$.

## C.1 Flat form

We consider the relation $A = \psi^2$. We set in a general way $\psi = \{(x, y) \in \mathbb{R}^2 | y - f_{\psi}(x) = 0\}$ the space of $\psi$ with $f_{\psi}$ the function describing the bottom. The flat form $\psi = \{(x, y) \in \mathbb{R}^2 | y - c = 0\}$ deformed from $\varepsilon n$ is given by $\psi + \varepsilon n = \{(x, y) \in \mathbb{R}^2 | y - c - \varepsilon = 0\}$. It could be illustrated by the figure C1.
<figure><img src="image_12.png" /><figcaption>Figure C1: Illustration of $\psi$ and $\psi + \varepsilon n$ with the function $\psi: x \to c$</figcaption></figure>

Here we have:

$$ \psi = \{(x, y) \in \mathbb{R}^2 | y - c = 0\} \text{ and} \\ A = \psi^2 = \{(x, y) \in \mathbb{R}^2 | y - c^2 = 0\}, $$

then, we have, thank to the definition, on the one hand:

$$ \begin{aligned} \nabla_{\psi} A &= \lim_{\varepsilon \to 0} \left( \frac{1}{\varepsilon} [A(\psi + \varepsilon n) - A(\psi)] \right), \\ &= \lim_{\varepsilon \to 0} \left( \frac{1}{\varepsilon} [\psi^2 + 2\psi\varepsilon + \underbrace{\varepsilon^2}_{\to 0} - \psi^2] \right) && \text{because } A(\psi + \varepsilon n) = (\psi + \varepsilon)^2, \\ &= 2\psi, \end{aligned} $$

on the other hand:

$$ \nabla_X A.n = 2\psi\nabla_X(\psi).n = 2\psi \begin{pmatrix} 0 \\ 1 \end{pmatrix} \begin{pmatrix} 0 \\ 1 \end{pmatrix} = 2\psi, \quad (C1) $$

and therefore

$$ \nabla_X A.n = \nabla_{\psi} A. \quad (C2) $$

## C.2 Linear form

The linear form $\{\psi = (x, y) \in \mathbb{R}^2 | y - ax - b = 0\}$ deformed by $\varepsilon n$ is given by $\psi + \varepsilon n$. We consider the relation $A = \cos(\psi)$. We set in a general way $\psi = \{(x, y) \in \mathbb{R}^2 | y - f_{\psi}(x) = 0\}$ the space of $\psi$ with $f_{\psi}$ the function describing the bottom. It could be illustrated by the figure C2.

We know the point $P(0, b)$ is contained on the line. The point $P' = P + \varepsilon n$ is therefore contained on the new translated line. We deduce the equation rapidly:

$$ \psi + \varepsilon n = ax + \varepsilon \sqrt{a^2 + 1} + b. \quad (C3) $$
<figure><img src="image_13.png" /><figcaption>Figure C1: Illustration of $\psi$ and $\psi + \varepsilon n$ with the function $\psi: x \to c$</figcaption></figure>

Here we have:

$\psi = \{(x, y) \in \mathbb{R}^2 | y - c = 0\}$ and

$A = \psi^2 = \{(x, y) \in \mathbb{R}^2 | y - c^2 = 0\}$,

then, we have, thank to the definition, on the one hand:

$$ 
\begin{align*}
\nabla_{\psi} A &= \lim_{\varepsilon \to 0} \left( \frac{1}{\varepsilon} [A(\psi + \varepsilon n) - A(\psi)] \right), \\
&= \lim_{\varepsilon \to 0} \left( \underbrace{\frac{1}{\varepsilon} |\psi|^2 + 2\psi\varepsilon + \varepsilon^2}_{\to 0} - \underbrace{|\psi|^2}_{\to 0} \right) && \text{because } A(\psi + \varepsilon n) = (\psi + \varepsilon)^2, \\
&= 2\psi,
\end{align*}
$$

on the other hand:

$$ \nabla_X A.n = 2\psi\nabla_X(\psi).n = 2\psi \begin{pmatrix} 0 \\ 1 \end{pmatrix} \begin{pmatrix} 0 \\ 1 \end{pmatrix} = 2\psi, \quad (C1) $$

and therefore

$$ \nabla_X A.n = \nabla_{\psi} A. \quad (C2) $$

## C.2 Linear form

The linear form $\{\psi = (x, y) \in \mathbb{R}^2 | y - ax - b = 0\}$ deformed by $\varepsilon n$ is given by $\psi + \varepsilon n$. We consider the relation $A = \cos(\psi)$. We set in a general way $\psi = \{(x, y) \in \mathbb{R}^2 | y - f_{\psi}(x) = 0\}$ the space of $\psi$ with $f_{\psi}$ the function describing the bottom. It could be illustrated by the figure C2.

We know the point $P(0, b)$ is contained on the line. The point $P' = P + \varepsilon n$ is therefore contained on the new translated line. We deduce the equation rapidly:

$$ \psi + \varepsilon n = ax + \varepsilon \sqrt{a^2 + 1} + b. \quad (C3) $$
<figure><img src="image_13.png" /><figcaption>**Figure C2:** Illustration of $\psi$, $\psi + \varepsilon n$ and $A$ with the function $\psi: x \to ax + b$ and $A = \cos(\psi)$</figcaption></figure>

Let us check the equation (17) for $A = \cos(\psi)$. On the one hand, we have:

$$ 
\begin{aligned}
\nabla_{\psi} A &= \lim_{\varepsilon \to 0} \left( \frac{1}{\varepsilon} [A(\psi + \varepsilon n) - A(\psi)] \right), \\
&= \lim_{\varepsilon \to 0} \left( \frac{1}{\varepsilon} [\cos(ax + \varepsilon\sqrt{a^2+1} + b) - \cos(ax + b)] \right), \\
&= \lim_{\varepsilon \to 0} \left( \underbrace{\frac{1}{\varepsilon}[\cos(ax + b)[\cos(\varepsilon\sqrt{a^2+1}) - 1]] - \sin(ax + b)[\sin(\varepsilon\sqrt{a^2+1})]}_{\to \varepsilon\sqrt{a^2+1}} \right), \\
&= \lim_{\varepsilon \to 0} \left( \underbrace{-\cos(ax + b)\varepsilon(a^2+1)}_{\to 0} - \sin(ax + b)\sqrt{a^2+1} \right), \\
&= -\sin(\psi)\sqrt{a^2+1}.
\end{aligned}
$$

On the other hand, we have:

$$ \nabla_X A = -\sin(\psi)\nabla_X \psi = -\sin(\psi) \begin{pmatrix} -a \\ 1 \end{pmatrix}, \quad (C4) $$

and therefore:

$$ \nabla_X A.n = -\frac{\sin(\psi)}{\sqrt{a^2+1}} \begin{pmatrix} -a \\ 1 \end{pmatrix} \begin{pmatrix} -a \\ 1 \end{pmatrix} = -\sin(\psi) \frac{(a^2+1)}{\sqrt{a^2+1}} = -\sin(\psi)\sqrt{a^2+1}. \quad (C5) $$

The equality: $\nabla_{\psi} A = \nabla_{X} A.n$ is still verified.
<figure><img src="image_14.png" /><figcaption>**Figure C2:** Illustration of $\psi$, $\psi + \varepsilon n$ and $A$ with the function $\psi: x \to ax + b$ and $A = \cos(\psi)$</figcaption></figure>

Let us check the equation (17) for $A = \cos(\psi)$. On the one hand, we have:

$$ 
\begin{aligned}
\nabla_{\psi} A &= \lim_{\varepsilon \to 0} \left( \frac{1}{\varepsilon} [A(\psi + \varepsilon n) - A(\psi)] \right), \\
&= \lim_{\varepsilon \to 0} \left( \frac{1}{\varepsilon} [\cos(ax + \varepsilon\sqrt{a^2+1} + b) - \cos(ax + b)] \right), \\
&= \lim_{\varepsilon \to 0} \left( \underbrace{\frac{1}{\varepsilon}[\cos(ax + b)[\cos(\varepsilon\sqrt{a^2+1}) - 1] - \sin(ax + b)\sin(\varepsilon\sqrt{a^2+1})]}_{\to \varepsilon\sqrt{a^2+1}} \right), \\
&= \lim_{\varepsilon \to 0} \left( \underbrace{-\cos(ax + b)\varepsilon(a^2+1)}_{\to 0} - \sin(ax + b)\sqrt{a^2+1} \right), \\
&= -\sin(\psi)\sqrt{a^2+1}.
\end{aligned}
$$

On the other hand, we have:

$$ \nabla_X A = -\sin(\psi)\nabla_X \psi = -\sin(\psi) \begin{pmatrix} -a \\ 1 \end{pmatrix}, \quad (C4) $$

and therefore:

$$ \nabla_X A.n = - \frac{\sin(\psi)}{\sqrt{a^2+1}} \begin{pmatrix} -a \\ 1 \end{pmatrix} \begin{pmatrix} -a \\ 1 \end{pmatrix} = -\sin(\psi) \frac{(a^2+1)}{\sqrt{a^2+1}} = -\sin(\psi)\sqrt{a^2+1}. \quad (C5) $$

The equality: $\nabla_{\psi} A = \nabla_{X} A.n$ is still verified.
## D Extended Shoaling model

The last Shoaling model presented in Cook (2021) and Dupont et al. (2023) had certain limitations. This model was therefore improved to give birth to the extended model below:

$$H(x,t) = \begin{cases} H_0(x,t)K_S(x,t) & \text{for } x \in \Omega_S \\ \mathcal{F}(\gamma h(x,t)) & \text{for } x \in \Omega_B, \end{cases} \quad (D1)$$

where $\mathcal{F}$ is a numerical parameterization function of the breaking defines below (D2):

$$\mathcal{F}(\gamma h(x,t)) = H(x_{\text{start}}) + [H(x_{\text{stop}}) - H(x_{\text{start}})] \cdot f\left(\frac{x - x_{\text{start}}}{x_{\text{stop}} - x_{\text{start}}}\right) \cdot g\left(\frac{h_{\text{max}} - h}{h_{\text{max}} - h_{\text{min}}}\right), \quad (D2)$$

with $x \in \Omega_B = [x_{\text{start}}, x_{\text{stop}}]$, $h \in [h_{\text{min}}, h_{\text{max}}]$ and the notations on figure D1. $H_{\text{start}}$ and $H_{\text{stop}}$ are the wave heights at the beginning and the end of the surf
## D Extended Shoaling model

The last Shoaling model presented in Cook (2021) and Dupont et al. (2023) had certain limitations. This model was therefore improved to give birth to the extended model below:

$$ H(x,t) = \begin{cases} H_0(x,t)K_S(x,t) & \text{for } x \in \Omega_S \\ \mathcal{F}(\gamma h(x,t)) & \text{for } x \in \Omega_B, \end{cases} \qquad (D1) $$

where $\mathcal{F}$ is a numerical parameterization function of the breaking defines below (D2):

$$ \mathcal{F}(\gamma h(x,t)) = H(x_{\text{start}}) + [H(x_{\text{stop}}) - H(x_{\text{start}})] \cdot f\left(\frac{x - x_{\text{start}}}{x_{\text{stop}} - x_{\text{start}}}\right) \cdot g\left(\frac{h_{\text{max}} - h}{h_{\text{max}} - h_{\text{min}}}\right), \qquad (D2) $$

with $x \in \Omega_B = [x_{\text{start}}, x_{\text{stop}}]$, $h \in [h_{\text{min}}, h_{\text{max}}]$ and the notations on figure D1. $H_{\text{start}}$ and $H_{\text{stop}}$ are the wave heights at the beginning and the end of the surf
<figure><img src="image_14.png" /><figcaption>**Figure D1:** Illustration of notations.</figcaption></figure>

zone on the domain $\Omega_B = [x_{start}, x_{stop}]$. The first function $f$ gives an account of breaking without taking into account the bed shape. It simply gives the appearance of breaking. The second function $g$ takes into account the seabed and interacts with it. Note that if $f$ and $g$ are the identity functions, we find the linear breaking $\gamma h(x, t)$ illustrated on figure D1. We can present in figure D2 some of these functions that set the breaking.

<figure><img src="image_15.png" /><figcaption>**Figure D2:** Illustration of $f_1$, $f_2$ and $g$ defined in $[0, 1] \to [0, 1]$.</figcaption></figure>

These functions were chosen to try to capture a natural breaking. They have no physical meaning. It is necessary to stipulate that the model first locate all the $\Omega_B$ domains and then apply the equation (D2) on each of them. This type of model gives us the LIP - 1C simulation figure D3, which gives very similar results to those produced by SWAN and XBeach.
<figure><img src="image_15.png" /><figcaption>**Figure D1:** Illustration of notations.</figcaption></figure>

zone on the domain $\Omega_B = [x_{start}, x_{stop}]$. The first function $f$ gives an account of breaking without taking into account the bed shape. It simply gives the appearance of breaking. The second function $g$ takes into account the seabed and interacts with it. Note that if $f$ and $g$ are the identity functions, we find the linear breaking $\gamma h(x, t)$ illustrated on figure D1. We can present in figure D2 some of these functions that set the breaking.

<figure><img src="image_16.png" /><figcaption>**Figure D2:** Illustration of $f_1$, $f_2$ and $g$ defined in $[0, 1] \to [0, 1]$.</figcaption></figure>

These functions were chosen to try to capture a natural breaking. They have no physical meaning. It is necessary to stipulate that the model first locate all the $\Omega_B$ domains and then apply the equation (D2) on each of them. This type of model gives us the LIP - 1C simulation figure D3, which gives very similar results to those produced by SWAN and XBeach.
<figure><img src="image_16.png" /></figure>

**Figure D3:** Wave results obtained with the Shoaling, SWAN and XBeach models. Sea bottom configuration from the LIP 1C channel experiment. Black points, measured HRMS, blue H from SWAN ($H_{RMSE} = 0.033$ m), red H from XBeach ($H_{RMSE} = 0.028$ m).

## References

Alsina, José M. and Iván Cáceres (2011). “Sediment suspension events in the inner surf and swash zone. Measurements in large-scale and high-energy wave conditions”. In: *Coastal Engineering* 58.8, pp. 657–670. DOI: https://doi.org/10.1016/j.coastaleng.2011.03.002.

Bascom, Willard N (1951). “The relationship between sand size and beach-face slope”. In: *Eos, Transactions American Geophysical Union* 32.6, pp. 866–874.

Battjes, Jurjen A and JPFM Janssen (1978). “Energy loss and set-up due to breaking of random waves”. In: *Coastal engineering* 1978, pp. 569–587.

Blenkinsopp, Chris E, Paul M Bayle, Daniel C Conley, Gerd Masselink, Emily Gulson, Isabel Kelly, Rafael Almar, Ian L Turner, Tom E Baldock, Tomas Beuzen, et al. (2021). “High-resolution, large-scale laboratory measurements of a sandy beach and dynamic cobble berm revetment”. In: *Scientific Data* 8.1, p. 22.

Booj, N, LH Holthuijsen, and RC Ris (1996). “The” SWAN” wave model for shallow water”. In: *Coastal Engineering* 1996, pp. 668–676.

Bouharguane, Afaf, Pascal Azerad, Frédéric Bouchette, Fabien Marche, and Bijan Mohammadi (2010). “Low complexity shape optimization & a pos-
<figure><img src="image_17.png" /></figure>

**Figure D3:** Wave results obtained with the Shoaling, SWAN and XBeach models. Sea bottom configuration from the LIP 1C channel experiment. Black points, measured HRMS, blue H from SWAN ($H_{RMSE} = 0.033$ m), red H from XBeach ($H_{RMSE} = 0.028$ m).

## References

Alsina, José M. and Iván Cáceres (2011). "Sediment suspension events in the inner surf and swash zone. Measurements in large-scale and high-energy wave conditions". In: *Coastal Engineering* 58.8, pp. 657-670. DOI: https://doi.org/10.1016/j.coastaleng.2011.03.002.

Bascom, Willard N (1951). "The relationship between sand size and beach-face slope". In: *Eos, Transactions American Geophysical Union* 32.6, pp. 866-874.

Battjes, Jurjen A and JPFM Janssen (1978). "Energy loss and set-up due to breaking of random waves". In: *Coastal engineering* 1978, pp. 569-587.

Blenkinsopp, Chris E, Paul M Bayle, Daniel C Conley, Gerd Masselink, Emily Gulson, Isabel Kelly, Rafael Almar, Ian L Turner, Tom E Baldock, Tomas Beuzen, et al. (2021). "High-resolution, large-scale laboratory measurements of a sandy beach and dynamic cobble berm revetment". In: *Scientific Data* 8.1, p. 22.

Booj, N, LH Holthuijsen, and RC Ris (1996). "The" SWAN" wave model for shallow water". In: *Coastal Engineering* 1996, pp. 668-676.

Bouharguane, Afaf, Pascal Azerad, Frédéric Bouchette, Fabien Marche, and Bijan Mohammadi (2010). "Low complexity shape optimization & a pos-
teriori high fidelity validation”. In: *Discrete and Continuous Dynamical Systems* - B 13.4, pp. 759–772. DOI: 10.3934/dcdsb.2010.13.759.

Bouharguane, Afaf and Bijan Mohammadi (2012). “Minimisation principles for the evolution of a soft sea bed interacting with a shallow sea”. In: *International Journal of Computational Fluid Dynamics* 26.3, pp. 163–172. DOI: 10.1080/10618562.2012.669831.

Chen, W., J.J. van der Werf, and S.J.M.H. Hulscher (2023). “A review of practical models of sand transport in the swash zone”. In: *Earth-Science Reviews* 238, p. 104355. DOI: https://doi.org/10.1016/j.earscirev.2023.104355.

Cook, Megan (2021). “Calcul optimal pour la modélisation de la dynamique naturelle des plages sableuses et la conception d’ouvrages de défense du littoral à faible impact anthropique”. Theses. Université Montpellier.

Cook, Megan, Frédéric Bouchette, Bijan Mohammadi, Léa Sprunck, and Nico-las Fraysse (2021). “Optimal Port Design Minimizing Standing Waves with A Posteriori Long Term Shoreline Sustainability Analysis”. en. In: *China Ocean Engineering* 35.6, pp. 802–813. DOI: 10.1007/s13344-021-0071-7.

Dean, Robert and Robert Dalrymple (2004). “Coastal Processes With Engi-neering Applications”. In: *Coastal Processes with Engineering Applications, by Robert G. Dean and Robert A. Dalrymple*, pp. 487. ISBN 0521602750. Cambridge, UK: Cambridge University Press, March 2004.

Dupont, Ronan, Megan Cook, Frédéric Bouchette, Bijan Mohammadi, and Samuel Meulé (2023). “Sandy beach dynamics by constrained wave energy minimization”. In: *Ocean Modelling*, p. 102197. DOI: https://doi.org/10.1016/j.ocemod.2023.102197.

Eichentopf, Sonja, Iván Cáceres, and José M. Alsina (2018). “Breaker bar morphodynamics under erosive and accretive wave conditions in large-scale experiments”. In: *Coastal Engineering* 138, pp. 36–48. DOI: https://doi.org/10.1016/j.coastaleng.2018.04.010.

Fowler, A. C. (2001). “Dunes and Drumlins”. In: *Geomorphological Fluid Mechanics*. Ed. by N. J. Balmforth and A. Provenzale. Berlin, Heidelberg: Springer Berlin Heidelberg, pp. 430–454. DOI: 10.1007/3-540-45670-8_18.

Fredsøe, Jørgen and Rolf Deigaard (1992). *Mechanics of Coastal Sediment Transport*. WORLD SCIENTIFIC. DOI: 10.1142/1546.

Gâteaux, René (1913). “Sur les fonctionnelles continues et les fonctionnelles analytiques”. In: *CR Acad. Sci. Paris* 157.325-327, p. 65.

Griewank, Andreas and Andrea Walther (2008). *Evaluating derivatives: principles and techniques of algorithmic differentiation*. SIAM.

Hadamard, J (1914). “On the Maximum Module of a Function and Its Derivatives”. In: *CR Acad. Sci. Paris* 41, pp. 68–72.
teriori high fidelity validation”. In: *Discrete and Continuous Dynamical Systems* - B 13.4, pp. 759–772. DOI: 10.3934/dcdsb.2010.13.759.

Bouharguane, Afaf and Bijan Mohammadi (2012). “Minimisation principles for the evolution of a soft sea bed interacting with a shallow sea”. In: *International Journal of Computational Fluid Dynamics* 26.3, pp. 163–172. DOI: 10.1080/10618562.2012.669831.

Chen, W., J.J. van der Werf, and S.J.M.H. Hulscher (2023). “A review of practical models of sand transport in the swash zone”. In: *Earth-Science Reviews* 238, p. 104355. DOI: https://doi.org/10.1016/j.earscirev.2023.104355.

Cook, Megan (2021). “Calcul optimal pour la modélisation de la dynamique naturelle des plages sableuses et la conception d’ouvrages de défense du littoral à faible impact anthropique”. Theses. Université Montpellier.

Cook, Megan, Frédéric Bouchette, Bijan Mohammadi, Léa Sprunck, and Nico-las Fraysse (2021). “Optimal Port Design Minimizing Standing Waves with A Posteriori Long Term Shoreline Sustainability Analysis”. en. In: *China Ocean Engineering* 35.6, pp. 802–813. DOI: 10.1007/s13344-021-0071-7.

Dean, Robert and Robert Dalrymple (2004). “Coastal Processes With Engi-neering Applications”. In: *Coastal Processes with Engineering Applications, by Robert G. Dean and Robert A. Dalrymple*, pp. 487. ISBN 0521602750. Cambridge, UK: Cambridge University Press, March 2004.

Dupont, Ronan, Megan Cook, Frédéric Bouchette, Bijan Mohammadi, and Samuel Meulé (2023). “Sandy beach dynamics by constrained wave energy minimization”. In: *Ocean Modelling*, p. 102197. DOI: https://doi.org/10.1016/j.ocemod.2023.102197.

Eichentopf, Sonja, Iván Cáceres, and José M. Alsina (2018). “Breaker bar morphodynamics under erosive and accretive wave conditions in large-scale experiments”. In: *Coastal Engineering* 138, pp. 36–48. DOI: https://doi.org/10.1016/j.coastaleng.2018.04.010.

Fowler, A. C. (2001). “Dunes and Drumlins”. In: *Geomorphological Fluid Mechanics*. Ed. by N. J. Balmforth and A. Provenzale. Berlin, Heidelberg: Springer Berlin Heidelberg, pp. 430–454. DOI: 10.1007/3-540-45670-8_18.

Fredsøe, Jørgen and Rolf Deigaard (1992). *Mechanics of Coastal Sediment Transport*. WORLD SCIENTIFIC. DOI: 10.1142/1546.

Gâteaux, René (1913). “Sur les fonctionnelles continues et les fonctionnelles analytiques”. In: *CR Acad. Sci. Paris* 157.325-327, p. 65.

Griewank, Andreas and Andrea Walther (2008). *Evaluating derivatives: principles and techniques of algorithmic differentiation*. SIAM.

Hadamard, J (1914). “On the Maximum Module of a Function and Its Derivatives”. In: *CR Acad. Sci. Paris* 41, pp. 68–72.
Hascoet, L. and V. Pascual (2004). "Tapenade user's guide". In: *INRIA Technical report*. INRIA, pp. 1-31.

Hasselmann, Klaus, T.P. Barnett, E. Bouws, H. Carlson, D. Cartwright, K. Enke, J.A. Ewing, H. Gienapp, D. Hasselmann, P. Kruseman, A. Meerburg, Peter Muller, Dirk Olbers, K. Richren, W. Sell, and H. Walden (1973). *Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP)*. Tech. rep. Hamburg, Germany: Deutsches Hydrographisches Institut, pp. 1-95.

Hattori, Masataro and Ryoichi Kawamata (1980). "Onshore-Offshore Transport and Beach Profile Change". In: *Coastal Engineering* 1980, pp. 1175-1193. DOI: 10.1061/9780872622647.072.

Isèbe, Damien, Pascal Azerad, Frederic Bouchette, Benjamin Ivorra, and Bijan Mohammadi (2008a). "Shape optimization of geotextile tubes for sandy beach protection". In: *International Journal for Numerical Methods in Engineering* 74.8, pp. 1262-1277. DOI: 10.1002/nme.2209.

Isèbe, Damien, Pascal Azérad, Frédéric Bouchette, and Bijan Mohammadi (2014). "Design of Passive Defense Structures in Coastal Engineering". en. In: *International Review of Civil Engineering (IRECE)* 5.2, p. 75. DOI: 10.15866/irece.v5i2.2029.

Isèbe, Damien, Pascal Azérad, Bijan Mohammadi, and Frédéric Bouchette (2008b). "Optimal shape design of defense structures for minimizing short wave impact". In: *Coastal Engineering* 55.1, pp. 35-46. DOI: 10.1016/j.coastaleng.2007.06.006.

Kouakou, Kouamé Kan Jacques and Pierre-Yves Lagrée (2006). "Evolution of a model dune in a shear flow". In: *European Journal of Mechanics - B/Fluids* 25.3, pp. 348-359. DOI: https://doi.org/10.1016/j.euromechflu.2005.09.002.

Marchesiello, Patrick, Julien Chauchat, Hassan Shafiei, Rafael Almar, Rachid Benshila, Franck Dumas, and Laurent Debreu (2022). "3D wave-resolving simulation of sandbar migration". In: *Ocean Modelling* 180, p. 102127. DOI: https://doi.org/10.1016/j.ocemod.2022.102127.

Martins, Kévin, Philippe Bonneton, Paul M. Bayle, Chris E. Blenkinsopp, Arthur Mouragues, and Hervé Michallet (2020). "Surf Zone Wave Measurements from Lidar Scanners: Analysis of Non-hydrostatic Processes". In: *Journal of Coastal Research* 95.SI, pp. 1189-1194. DOI: 10.2112/SI95-231.1.

McCall, Robert T, JSM Van Thiel De Vries, NG Plant, AR Van Dongeren, JA Roelvink, DM Thompson, and AJHM Reniers (2010). "Two-dimensional time dependent hurricane overwash and erosion modeling at Santa Rosa Island". In: *Coastal Engineering* 57.7, pp. 668-683.

Mohammadi, Bijan (2007). "Global optimization, level set dynamics, incomplete sensitivity and regularity control". en. In: *International Journal of Computational Fluid Dynamics* 21.2, pp. 61-68. DOI: 10.1080/10618560701455855.
Hascoet, L. and V. Pascual (2004). "Tapenade user's guide". In: *INRIA Technical report*. INRIA, pp. 1-31.

Hasselmann, Klaus, T.P. Barnett, E. Bouws, H. Carlson, D. Cartwright, K. Enke, J.A. Ewing, H. Gienapp, D. Hasselmann, P. Kruseman, A. Meerburg, Peter Muller, Dirk Olbers, K. Richren, W. Sell, and H. Walden (1973). *Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP)*. Tech. rep. Hamburg, Germany: Deutsches Hydrographisches Institut, pp. 1-95.

Hattori, Masataro and Ryoichi Kawamata (1980). "Onshore-Offshore Transport and Beach Profile Change". In: *Coastal Engineering* 1980, pp. 1175-1193. DOI: 10.1061/9780872622647.072.

Isèbe, Damien, Pascal Azerad, Frederic Bouchette, Benjamin Ivorra, and Bijan Mohammadi (2008a). "Shape optimization of geotextile tubes for sandy beach protection". In: *International Journal for Numerical Methods in Engineering* 74.8, pp. 1262-1277. DOI: 10.1002/nme.2209.

Isèbe, Damien, Pascal Azérad, Frédéric Bouchette, and Bijan Mohammadi (2014). "Design of Passive Defense Structures in Coastal Engineering". en. In: *International Review of Civil Engineering (IRECE)* 5.2, p. 75. DOI: 10.15866/irece.v5i2.2029.

Isèbe, Damien, Pascal Azérad, Bijan Mohammadi, and Frédéric Bouchette (2008b). "Optimal shape design of defense structures for minimizing short wave impact". In: *Coastal Engineering* 55.1, pp. 35-46. DOI: 10.1016/j.coastaleng.2007.06.006.

Kouakou, Kouamé Kan Jacques and Pierre-Yves Lagrée (2006). "Evolution of a model dune in a shear flow". In: *European Journal of Mechanics - B/Fluids* 25.3, pp. 348-359. DOI: https://doi.org/10.1016/j.euromechflu.2005.09.002.

Marchesiello, Patrick, Julien Chauchat, Hassan Shafiei, Rafael Almar, Rachid Benshila, Franck Dumas, and Laurent Debreu (2022). "3D wave-resolving simulation of sandbar migration". In: *Ocean Modelling* 180, p. 102127. DOI: https://doi.org/10.1016/j.ocemod.2022.102127.

Martins, Kévin, Philippe Bonneton, Paul M. Bayle, Chris E. Blenkinsopp, Arthur Mouragues, and Hervé Michallet (2020). "Surf Zone Wave Measurements from Lidar Scanners: Analysis of Non-hydrostatic Processes". In: *Journal of Coastal Research* 95.SI, pp. 1189-1194. DOI: 10.2112/SI95-231.1.

McCall, Robert T, JSM Van Thiel De Vries, NG Plant, AR Van Dongeren, JA Roelvink, DM Thompson, and AJHM Reniers (2010). "Two-dimensional time dependent hurricane overwash and erosion modeling at Santa Rosa Island". In: *Coastal Engineering* 57.7, pp. 668-683.

Mohammadi, Bijan (2007). "Global optimization, level set dynamics, incomplete sensitivity and regularity control". en. In: *International Journal of Computational Fluid Dynamics* 21.2, pp. 61-68. DOI: 10.1080/10618560701455855.
Mohammadi, Bijan (2010). "Hadamard incomplete sensitivity and shape optimization". eng. In: *Control and Cybernetics* 39.3, pp. 615-626.

— (2017). "Uncertainty quantification in littoral erosion". In: *Computers & Fluids* 143, pp. 120-133. DOI: https://doi.org/10.1016/j.compfluid.2016.10.017.

Mohammadi, Bijan and Frédéric Bouchette (2014). "Extreme scenarios for the evolution of a soft bed interacting with a fluid using the Value at Risk of the bed characteristics". In: *Computers and Fluids* 89, pp. 78-87. DOI: 10.1016/j.compfluid.2013.10.021.

Mohammadi, Bijan and Afaf Bouharguane (2011). "Optimal dynamics of soft shapes in shallow waters". In: *Computers and Fluids* 40, pp. 291-298. DOI: 10.1016/j.compfluid.2010.09.031.

Mohammadi, Bijan and Olivier Pironneau (2009). *Applied Shape Optimization for Fluids*. Oxford University Press. DOI: 10.1093/acprof:oso/9780199546909.001.0001.

Munk, Walter (1949). "The solitary wave theory and its application to surf problems". In: *Annals of the New York Academy of Sciences* 51, pp. 376-424. DOI: 10.1111/j.1749-6632.1949.tb27281.x.

Murray, A. Brad (2007). "Reducing model complexity for explanation and prediction". In: *Geomorphology* 90.3, pp. 178-191. DOI: https://doi.org/10.1016/j.geomorph.2006.10.020.

Nielsen, Peter (1992). *Coastal bottom boundary layers and sediment transport*. Vol. 4. World scientific.

— (2002). "Shear stress and sediment transport calculations for swash zone modelling". In: *Coastal Engineering* 45.1, pp. 53-60. DOI: https://doi.org/10.1016/S0378-3839(01)00036-9.

Paola, C. and V. R. Voller (2005). "A generalized Exner equation for sediment mass balance". In: *Journal of Geophysical Research: Earth Surface* 110.F4. DOI: https://doi.org/10.1029/2004JF000274.

Pender, Doug and Harshinie Karunarathna (2013). "A statistical-process based approach for modelling beach profile variability". In: *Coastal Engineering* 81, pp. 19-29.

Quick, M. (1991). "Onshore-offshore sediment transport on beaches". In: *Coastal Engineering* 15, pp. 313-332.

Reineck, H.-E. and I. B. Singh (1973). *Depositional sedimentary environments; with reference to terrigenous clastics [by] H.-E. Reineck [and] I. B. Singh*. English. Springer-Verlag Berlin, New York, xvi, 439 p.

Roelvink, Dano J.A. (1993). "Dissipation in random wave groups incident on a beach". In: *Coastal Engineering - COAST ENG* 19, pp. 127-150. DOI: 10.1016/0378-3839(93)90021-Y.

Roelvink, Dano J.A., Ad Reniers, Ap van Dongeren, Jaap Thiel de Vries, Robert McCall, and Jamie Lescinski (2009). "Modelling storm impacts on
Mohammadi, Bijan (2010). "Hadamard incomplete sensitivity and shape optimization". eng. In: *Control and Cybernetics* 39.3, pp. 615-626.

— (2017). "Uncertainty quantification in littoral erosion". In: *Computers & Fluids* 143, pp. 120-133. DOI: https://doi.org/10.1016/j.compfluid.2016.10.017.

Mohammadi, Bijan and Frédéric Bouchette (2014). "Extreme scenarios for the evolution of a soft bed interacting with a fluid using the Value at Risk of the bed characteristics". In: *Computers and Fluids* 89, pp. 78-87. DOI: 10.1016/j.compfluid.2013.10.021.

Mohammadi, Bijan and Afaf Bouharguane (2011). "Optimal dynamics of soft shapes in shallow waters". In: *Computers and Fluids* 40, pp. 291-298. DOI: 10.1016/j.compfluid.2010.09.031.

Mohammadi, Bijan and Olivier Pironneau (2009). *Applied Shape Optimization for Fluids*. Oxford University Press. DOI: 10.1093/acprof:oso/9780199546909.001.0001.

Munk, Walter (1949). "The solitary wave theory and its application to surf problems". In: *Annals of the New York Academy of Sciences* 51, pp. 376-424. DOI: 10.1111/j.1749-6632.1949.tb27281.x.

Murray, A. Brad (2007). "Reducing model complexity for explanation and prediction". In: *Geomorphology* 90.3, pp. 178-191. DOI: https://doi.org/10.1016/j.geomorph.2006.10.020.

Nielsen, Peter (1992). *Coastal bottom boundary layers and sediment transport*. Vol. 4. World scientific.

— (2002). "Shear stress and sediment transport calculations for swash zone modelling". In: *Coastal Engineering* 45.1, pp. 53-60. DOI: https://doi.org/10.1016/S0378-3839(01)00036-9.

Paola, C. and V. R. Voller (2005). "A generalized Exner equation for sediment mass balance". In: *Journal of Geophysical Research: Earth Surface* 110.F4. DOI: https://doi.org/10.1029/2004JF000274.

Pender, Doug and Harshinie Karunarathna (2013). "A statistical-process based approach for modelling beach profile variability". In: *Coastal Engineering* 81, pp. 19-29.

Quick, M. (1991). "Onshore-offshore sediment transport on beaches". In: *Coastal Engineering* 15, pp. 313-332.

Reineck, H.-E. and I. B. Singh (1973). *Depositional sedimentary environments; with reference to terrigenous clastics [by] H.-E. Reineck [and] I. B. Singh*. English. Springer-Verlag Berlin, New York, xvi, 439 p.

Roelvink, Dano J.A. (1993). "Dissipation in random wave groups incident on a beach". In: *Coastal Engineering - COAST ENG* 19, pp. 127-150. DOI: 10.1016/0378-3839(93)90021-Y.

Roelvink, Dano J.A., Ad Reniers, Ap van Dongeren, Jaap Thiel de Vries, Robert McCall, and Jamie Lescinski (2009). "Modelling storm impacts on
beaches, dunes and barrier islands”. In: Coastal Engineering 56, pp. 1133–1152. DOI: 10.1016/j.coastaleng.2009.08.006.

Roelvink, J.A. (2006). “Coastal morphodynamic evolution techniques”. In: *Coastal Engineering* 53.2, pp. 277–287. DOI: https://doi.org/10.1016/j.coastaleng.2005.10.015.

Roelvink, JA and Ad Reniers (1995). *LIP 11D Delta Flume experiments: a dataset for profile model validation*. WL/Delft Hydraulics.

Rooijen, Arnold van, Ad Reniers, Jaap van Thiel de Vries, Chris Blenkinsopp, and Robert McCall (2012). “Modeling swash zone sediment transport at Truc Vert Beach”. In: *Coastal Engineering Proceedings* 1.33, sediment.105. DOI: 10.9753/icce.v33.sediment.105.

Schimmels, Stefan and Chris Blenkinsopp (2020). “DynaRev-Dynamic Coastal Protection: Resilience of Dynamic Revetments Under Sea Level Rise”. In.

Shafiei, Hassan, Julien Chauchat, Cyrille Bonamy, and Patrick Marchesiello (2023). “Adaptation of the SANTOSS transport formula for 3D nearshore models: Application to cross-shore sandbar migration”. In: *Ocean Modelling* 181, p. 102138. DOI: https://doi.org/10.1016/j.ocemod.2022.102138.

Short, Andrew D. (1996). “The role of wave height, period, slope, tide range and embaymentisation in beach classifications: A review”. In: *Revista Chilena de Historia Natural* 69, pp. 589–604.

Soulsby, R.L. (1987). “Calculating bottom orbital velocity beneath waves”. In: *Coastal Engineering - COAST ENG* 11, pp. 371–380. DOI: 10.1016/0378-3839(87)90034-2.

Sous, Damien, Guillaume Dodet, Frédéric Bouchette, and M Tissier (2020). “Momentum balance across a barrier reef”. In: *Journal of Geophysical Research: Oceans* 125.2, e2019JC015503.

Taylor, Jean E. (1976). “The Structure of Singularities in Soap-Bubble-Like and Soap-Film-Like Minimal Surfaces”. In: *Annals of Mathematics* 103.3, pp. 489–539.

Vos, Kilian, Mitchell D Harley, Kristen D Splinter, Andrew Walker, and Ian L Turner (2020). “Beach slopes from satellite-derived shorelines”. In: *Geophysical Research Letters* 47.14, e2020GL088365.

Wiberg, Patricia L. and Christopher R. Sherwood (2008). “Calculating wave-generated bottom orbital velocities from surface-wave parameters”. In: *Computers & Geosciences* 34.10, pp. 1243–1262. DOI: https://doi.org/10.1016/j.cageo.2008.02.010.

Wright, L.D and A.D Short (1984). “Morphodynamic variability of surf zones and beaches: A synthesis”. In: *Marine Geology* 56.1, pp. 93–118. DOI: https://doi.org/10.1016/0025-3227(84)90008-2.

Yang, Chih Ted et al. (1996). “Sediment transport: theory and practice”. In.
beaches, dunes and barrier islands”. In: Coastal Engineering 56, pp. 1133–1152. DOI: 10.1016/j.coastaleng.2009.08.006.

Roelvink, J.A. (2006). “Coastal morphodynamic evolution techniques”. In: *Coastal Engineering* 53.2, pp. 277–287. DOI: https://doi.org/10.1016/j.coastaleng.2005.10.015.

Roelvink, JA and Ad Reniers (1995). *LIP 11D Delta Flume experiments: a dataset for profile model validation*. WL/Delft Hydraulics.

Rooijen, Arnold van, Ad Reniers, Jaap van Thiel de Vries, Chris Blenkinsopp, and Robert McCall (2012). “Modeling swash zone sediment transport at Truc Vert Beach”. In: *Coastal Engineering Proceedings* 1.33, sediment.105. DOI: 10.9753/icce.v33.sediment.105.

Schimmels, Stefan and Chris Blenkinsopp (2020). “DynaRev-Dynamic Coastal Protection: Resilience of Dynamic Revetments Under Sea Level Rise”. In.

Shafiei, Hassan, Julien Chauchat, Cyrille Bonamy, and Patrick Marchesiello (2023). “Adaptation of the SANTOSS transport formula for 3D nearshore models: Application to cross-shore sandbar migration”. In: *Ocean Modelling* 181, p. 102138. DOI: https://doi.org/10.1016/j.ocemod.2022.102138.

Short, Andrew D. (1996). “The role of wave height, period, slope, tide range and embaymentisation in beach classifications: A review”. In: *Revista Chilena de Historia Natural* 69, pp. 589–604.

Soulsby, R.L. (1987). “Calculating bottom orbital velocity beneath waves”. In: *Coastal Engineering - COAST ENG* 11, pp. 371–380. DOI: 10.1016/0378-3839(87)90034-2.

Sous, Damien, Guillaume Dodet, Frédéric Bouchette, and M Tissier (2020). “Momentum balance across a barrier reef”. In: *Journal of Geophysical Research: Oceans* 125.2, e2019JC015503.

Taylor, Jean E. (1976). “The Structure of Singularities in Soap-Bubble-Like and Soap-Film-Like Minimal Surfaces”. In: *Annals of Mathematics* 103.3, pp. 489–539.

Vos, Kilian, Mitchell D Harley, Kristen D Splinter, Andrew Walker, and Ian L Turner (2020). “Beach slopes from satellite-derived shorelines”. In: *Geophysical Research Letters* 47.14, e2020GL088365.

Wiberg, Patricia L. and Christopher R. Sherwood (2008). “Calculating wave-generated bottom orbital velocities from surface-wave parameters”. In: *Computers & Geosciences* 34.10, pp. 1243–1262. DOI: https://doi.org/10.1016/j.cageo.2008.02.010.

Wright, L.D and A.D Short (1984). “Morphodynamic variability of surf zones and beaches: A synthesis”. In: *Marine Geology* 56.1, pp. 93–118. DOI: https://doi.org/10.1016/0025-3227(84)90008-2.

Yang, Chih Ted et al. (1996). “Sediment transport: theory and practice”. In.