Rapid Changes in Permeability: Numerical Investigation into Storm-Driven Pebble Beach Morphodynamics with XBeach-G

Abstract

This study delves into the morphodynamic changes of pebble beaches in response to storm events, employing a combination of observational and numerical approaches. This research focuses on three extreme events, meticulously examining morhological changes in intertidal topography on the beach of Etretat (Normandy, France). A robust dataset of daily beach topography, derived from video monitoring systems, validates a set of numerical simulations of cross-shore dynamics performed by the process-based model XBeach-G. Our study evaluates the model’s efficacy in estimating beach profile evolution under high-energy conditions and explores its sensitivity to the physical properties of pebbles, including permeability. The results underscore the significance of considering spatial and temporal variations in permeability during storms to enhance the numerical model’s accuracy in predicting pebble beach dynamics. Furthermore, this study advocates for the incorporation of grain size mapping techniques to refine numerical model implementations.

1. Introduction

Coarse-grained coastal systems, particularly gravel and pebble beaches, constitute a significant proportion of the world’s mid- to high-latitude shorelines [1]. These beaches play a crucial role in coastal protection, serving as natural buffers against sediment loss and energy dissipation during high-energy wave conditions and extreme storm events [2,3,4]. The morphodynamic evolution of these systems is intricately tied to the steep beach face’s highly reflective nature and the substantial variation in coarse-grain distribution, influenced by energetic swash motions generated by breaking waves. This complexity introduces heightened risks of flooding and coastal erosion.

Numerous studies spanning several decades have sought to understand and predict the behavior of these coastal systems in response to extreme hydrodynamic drivers [4,5,6,7,8]. Powell’s [8] seminal work in 1990 introduced a parametric model for predicting the evolution of gravel beaches’ cross-shore topography based on extensive physical model tests, proving effective under moderate-energy conditions. Similarly, Bradbury et al. [9,10] developed an empirical approach to examine the short-term morphodynamic response of UK barrier beaches to incident wave conditions. However, these approaches have limitations, particularly in accounting for sediment grain size and system permeability, crucial aspects in understanding swash dynamics on coarse-grained beaches.

In general, research studies on gravel beach morphodynamics rely on field observations, numerical modeling, and physical experimentations. For example, Almeida et al. [11] employed laser scanning for observational purposes, revealing asymmetrical morphological responses during tidal phases, impacting wave breaking and sediment transport. In the context of extreme storm events causing significant morphological changes on coarse-grained beaches, remote sensing methods, such as coastal video imaging, have emerged as effective solutions for monitoring temporal beach [12,13,14,15,16,17].

The development of numerical approaches for capturing gravel beach morphodynamics has also been pivotal. The XBeach model [18], originally dedicated to the simulation of sandy shorelines under storm conditions, was first adapted to gravel systems by Jamal et al. [6]. Masselink et al. [19] and McCall et al. [20] then introduced XBeach-G to address the process-based simulation of gravel beaches and barriers more specifically. The model was validated with extensive field data and used worldwide in numerous studies [5,21,22,23,24].

The sensitivity of XBeach-G is closely linked to extrinsic drivers (e.g., wave characteristics and tidal regime) and to the intrinsic properties of the system, such as its slope, sediment sizes, and related permeability. Indeed, the steeper slope of gravel beaches induces significant sediment transport in the swash zone, caused by a high infiltration rate through the permeable sediment porosity, contributing to the beach face’s morphology. The continuous re-organization of gravels interacting with incident waves suggests that permeability is likely to vary significantly in time and space on gravel beaches; however, such variability is not considered by current models.

Beyond a sediment grain size of 4 mm, coarse-grained coastal systems transition from gravel to pebble beaches, characterized by substantial sediments and a steeper slope. The evolution of such beaches should be very sensitive to the variability in sediment sizes and to the subsequent varying permeability. Understanding the evolution of pebble beaches, their response to extreme climate drivers, and making informed decisions on coastal management become crucial.

This study aims to answer pertinent questions related to XBeach-G’s reliability in capturing the morphological response of intertidal topography in pebble beaches, the dependence of pebble mobility on extreme energy drivers, and the sensitivity of morphological response to pebble distribution, sorting, and permeability. The investigation focuses on the beach of Etretat (Normandy, France), examining the intertidal topography changes before and during individual storms. First, high-resolution video images were used to estimate the pre- and post-storm cross-shore morphological profiles of the beach. Then, we analyzed the beach profile’s sensitivity to pebble permeability through a comprehensive ensemble of numerical simulations. Finally, limitations of the model are addressed, providing insights into the physical behavior of pebbles.

2. Study Site: Pebble Beach of Etretat (Normandy, France)

The bay of Etretat, nestled along the Normandy coasts in the southern part of the English Channel (Figure 1), is an approximately 1000 m long pocket beach framed by bedrock structures of indurated limestone [25]. The beach sediment supply predominantly consists of pebbles ranging from 2 to 10 cm [26]. A subtidal sandy substrate (D50 = 0.8 mm) occasionally appears on the beach face following intense successive storm conditions [4], which result in a temporary transition from a pure gravel beach to a sandy-dominated system, according to the classification of Jennings and Shulmeiser [27].

The Channel, spanning about 560 km with varying widths ranging from around 240 km in the Atlantic Ocean to less than 35 km in the Strait of Dover, represents a shallow area, incorporating the continental shelf, with an average depth of 60 m in the central part and less than 20 m on the sides of the basin (Figure 1). The depth increases westward, exceeding 100 m but remaining distant from the continental margin.

The Normandy coasts are characterized by a large macrotidal range, averaging 9.5 m at spring tide [4]. Environments with tidal amplitudes greater than 8 m are classified as megatidal systems by Levoy et al. [28], which reflects the significant hydrodynamics associated with such extreme daily variations in the sea surface elevation. Average sea state conditions are characterized by waves from the west and northeast, with heights of approximately 1 m and periods of 5–7 s. Seasonal variations manifest high energy during winter (October through March) and low energy during the rest of the year [4].

A video monitoring system (VMS) was deployed in June 2018 on the beach of Etretat, consisting of three cameras with different fields of view. Each camera records 6 images per hour during daylight, providing high-resolution images with an exposure time of 10 min at a resolution of 3840 × 2160 pixels. The combination of cameras creates a panoramic view of the beach, covering the full bay. As the footprint of a pixel decreases with the distance from the camera, the region of interest was limited to a rectangular zone of approximately 200 × 500 m (Figure 1b). This ensures that the footprint remains under 1 m²/px and minimizes the vertical error of generated digital elevation models (DEMs).

According to Lopez et al. [29], dominant waves originate in the Atlantic Ocean, characterized by Gaussian behavior (Figure 2a) with a bimodal distribution, indicating the presence of sea and swell waves (Figure 2b). This behavior is consistent from deep water to shallow areas, with a reduction in average wave height ($H_s$) of around 50% and 25% of energy dissipation, respectively, from the Atlantic to the English Channel due to changes in water depth [29]. During storms, changes in wave direction occur, with an average clockwise variation of 18° at the peak of the storm, shifting from the sectors W-WNW to sectors WNW-NW. These changes are predominantly influenced by bathymetry variation and bottom friction, leading to wave refraction and energy dissipation processes [29].

3. Methodological Approach

3.1. Morphological Data

Three distinct data sources—SHOM [30], ROL (‘Réseau d’Observation du Littoral de Normandie et des Hauts-de-France’, [31]), and the video monitoring system (VMS)—contribute to morphological evolution observations across regional at local scales. The integration of these sources was crucial to establishing a comprehensive bathymetric cross-shore profile.

Offshore data from SHOM [30], provided on a regular grid with a resolution of 0.001° (~111 m), were refined to 3 m using linear interpolation for our study (Figure 3). Nearshore bathymetry (sub- to supratidal) was derived from aerial lidar surveys conducted by ROL [31], offering a resolution of 1 m with a vertical uncertainty of about 20 cm. The transition between SHOM and ROL (approximately 500 off the beach) data was smoothed using an interpolation method.

The intertidal topography was estimated daily using data from the video monitoring system (VMS) installed in Etretat, with a root mean square error (RMSE) of 22 cm on elevations [17] (Figure 3b). This section of the profile was measured by both ROL and VMS, with VMS-derived elevations used in the profile where available in the intertidal area. ROL data fill both sub- and supra-tidal parts of the profile. Notably, the extent of the VMS-derived DEM changes daily following tidal range variations and changes in the number of delineated shorelines after an image quality filtering process. Hence, there are two transitions between VMS and ROL data, at the upper and lower limits of the daily tidal range, which move from day to day. No smoothing was applied to these transitions, occasionally resulting in apparent jumps in the profile.

Coarse-grain size measurements were obtained using a deep learning approach developed by Soloy et al. [26]. Sediments, measured from the 9 June 2020 ortho-image, averaged D15, D50, and D90 of 4.59 cm, 5.99 cm, and 8.60 cm, respectively. While the embayment of Etretat tends to trap gravel within the bay [32], the overall grain size, unmeasured during storm dates, is assumed to remain close to the mentioned values. Past events of Etretat’s beach transitioning to mainly sandy conditions suggest that gravel may remain trapped in the subtidal part of the bay, although direct observations are lacking.

3.2. Hydrodynamic Data and Storm Identification

The hydrodynamic conditions of waves, winds, tides, and surges, close to Etretat, were assessed through a combination of in situ observations and numerical modeling.

The wave data were initially extracted from a wave model [29] as there is no active buoy close to the beaches of Normandy. Waves were then propagated to shore using the SWAN model [33]. Propagation was forced over a total of 26 months, between June 2018 and September 2020, on a regular rectangular grid of 20 m resolution. The model was validated against buoy data near Le Havre (~25 km south to southwest from Etretat) with an RMSE of 0.235 m and a correlation coefficient of 0.927. The time series of wave characteristics used for our XBeach-G models was extracted from the cell where the offshore end of the cross-shore bathymetric transect is located. A more precise description of SWAN model setting is provided in Solano et al. [29].

Astronomical tide, surge, and wind data were extracted from a Hycom2D model [34], built on a curvilinear grid with a resolution of 500 m near the coast. As for the SWAN data, the extraction point is located in the cell where the cross-shore transect ends. Regarding validation, Pasquet et al. [35] reported measuring a tidal phase shift of up to 12 min and an underestimation on the total water level of −10 cm at the peak of storm surges. This shift could induce a possible consequence on the estimation of the morphological elevation during storms. Nevertheless, this remains small compared to the errors generated by in situ measurements and corresponds to less than twice the median grain size.

The winters of 2018–2019 were characterized as the most energetic period in the last decade [36,37,38], causing substantial wave overtopping and seawall damage in Etretat [4] (Figure A1b). Storm events were identified in the hydrodynamic data record using an excessive wave significant height method. Identified storms were sorted according to their total energy. This method is often used to procedurally establish a sorted list of marine storm events [39,40,41]. In the present study, a marine storm was defined as an event during which the significant wave height ($H_s$) exceeds a threshold ($H_{s threshold}$) and lasts for at least 12 h. We introduced $H_{s threshold}$ as the mean value of $H_s$ plus two standard deviations (Equation (1)). In addition, consecutive events separated by an interval of less than 48 h were merged and considered as one longer event [42]. The list of storm events identified includes 30 storms distributed between June 2018 and September 2020. Sorting by decreasing accumulated energy ($E_{storm}$, Equation (3)) was then performed. $E_{storm}$ corresponds to the integral of the energy flux ($E_{flux}$, Equation (2)) over the duration of each event.

$$H_{s threshold}={\mu }_{H_s}+2{\sigma }_{H_s}^2$$

(1)

$$E_{flux}\left(t\right)=\frac{\varrho g^2}{64\pi }{H}_s^{2}T$$

(2)

$$E_{storm}={\int }_{t_{storm start}}^{t_{storm end}}{E}_{flux}\left(t\right) dt$$

(3)

where $\varrho$ is the sea water density (1025 kg/m³), $g$ is the acceleration of gravity (9.81 m/s²), $H_s$ is the wave significant height in meters, and $T$ is the wave period in seconds.

Once sorted, the three storms presenting both high intensity and sufficient bathymetry data were selected (Figure 4). Storm weather conditions tend to decrease the quality of VMS images (water drops, overtopping, depth of field reduced by rain and sea spray). This challenges the identification of an accurate coastline and creates gaps in our time series of intertidal bathymetry during storms. Nevertheless, DEMs could be reconstructed over a total of 11 days, including periods before and during the three selected storms. The selected events correspond to storms A, B, and C; their characteristics are summarized in

Table 1. Characteristics of modeled storms.

Storm Name	A	B	C
Storm Start Date	07 Dec 2018	13 Dec 2019	05 Dec 2019
Storm Duration	35 h	42 h	58 h
Simulation Start Date	06 Dec 2018	12 Dec 2019	04 Dec 2019
Simulation Duration	2 × 24 h	3 × 24 h	5 × 24 h
Wave Direction	104° N	104° N	109° N
Max Energy Flux	7.1 × 104 W/m/s	7.5 × 104 W/m/s	8.3 × 104 W/m/s
Storm Energy	8.9 × 109 W/m	7.0 × 109 W/m	5.4 × 109 W/m

. Figure A1a illustrates the wave heights of storm B, with a strong overtopping leading to the flooding of the backshore.

Figure 5 shows the evolution of $E_{flux}$ during the 3 selected storm events, as well as the dates on which the intertidal bathymetry data could be measured by the VMS (i.e., start, transition and end dates of the simulations). Only the dates where VMS-derived DEMs were wider than 35 m (an empirically chosen threshold) were considered usable to minimize the potential temporal decorrelation between our profile and the real bathymetry.

Storm A (Figure 5a) consisted of a series of four peaks increasing between 5 × 10⁴ and 7 × 10⁴ W/m/s, of which only the first (i.e., $E_{flux}$ = 5 × 10⁴ W/m/s) could be monitored and modeled. Since the calmer period preceding the storm was also covered by one VMS observation, three DEMs were used to simulate two days in using XBeach-G, each model lasting 24 h. Storm B (Figure 5b) consisted of three peaks varying between 5 × 10⁴ and 8 × 10⁴ W/m/s. It is covered from start to finish by four daily intertidal DEMs that allowed three consecutive days of simulation to be computed. Storm C (Figure 5c) consisted of only a single peak where $E_{flux}$ = 8 × 10⁴ W/m/s. Unfortunately, only observations prior to the peak are available, although the last available DEM covers the beginning of the event, when $E_{flux}$ reached 5 × 10⁴ W/m/s and was intensifying. However, the period preceding this maximum had another peak of lower intensity ($E_{flux}$ = 2 × 10⁴ W/m/s) with a good data temporal coverage. It was, therefore, decided to model this pre-storm period using the 6 DEMs available to simulate 5 days, as bathymetry was updated every 24 h.

3.3. XBeach-G Numerical Setup

XBeach is an open-source numerical model operating in the horizontal plane (2DH), addressing wave propagation, long waves, mean flow, sediment transport, and morphological changes [18]. Two main modes are utilized: nonhydrostatic (XB-NH) for resolving all wave motions (with higher computational expense) and Surfbeat (XB-SB) for resolving wave group scale motions more efficiently. The XB-SB mode, though not resolving sea-swell frequency motions, computes steady setup, (un)steady currents, and infra-gravity wave motions, especially relevant in dissipative events and fringing reef environments. XB-NH calculates depth-averaged flow due to waves and currents using non-linear shallow water equations, incorporating a nonhydrostatic pressure correction akin to a one-layer SWASH model.

The development of XBeach-G, a one-dimensional extension of XBeach dedicated to gravel barriers, was incentivized by the diverging results of XBeach when applied to gravel systems, suggesting that physical forcings needed to be considered differently in such contexts [19]. XBeach-G uses the NH mode to solve the non-linear shallow water equations (NLSWE) and calculate the depth-averaged flow due to waves and currents, taking groundwater exchange into account (Equation (4)).

$$\begin{gathered} \frac{\delta \zeta }{\delta t}+\frac{\delta hu}{\delta x}+S=0 \\ \frac{\delta u}{\delta t}+u\frac{\delta u}{\delta x}-\frac{\delta }{\delta x}\left({\nu }_h\frac{\delta u}{\delta x}\right)=-\frac{1}{\rho }\frac{\left(\rho \overline{q}+\rho g\zeta \right)}{\delta x}-\frac{{\tau }_b}{\rho h} \end{gathered}$$

(4)

where $t$ and $x$ represent the temporal and spatial (horizontal) coordinates, respectively, $\zeta$ is the elevation of the free surface, $h$ is the water depth, $u$ is the depth-averaged cross-shore velocity, $S$ is the exchange flux between surface water and groundwater (proportional to the permeability $k$), ${\nu }_h$ is the water viscosity (horizontal), $\rho$ is the water density, $\overline{q}$ is the depth-averaged dynamic pressure (normalized by the water density), $g$ is the acceleration due to gravity, and ${\tau }_b$ is the bed shear stress.

The aim of the present numerical investigation is to assess beach permeability evolution by modeling storms with varying permeability ($k$) within the range of 0 to 0.6 m/s. This range aligns with the expected permeability for the measured surficial grain size in Etretat [7,43,44]. Storms are segmented and modeled in 24 h intervals, approximately matching the daily frequency of intertidal bathymetry estimation from VMS observations.

The 3500 m long bathymetric profile (Figure 6a) is distributed on a 1D mesh of 1232 nodes with a variable spacing ranging between 3 m offshore and 1 m close to the beach. An empirically defined non-erodible layer immobilizes sections corresponding to the seawall and below elevation −6 m (Figure 6b) and the hardground bedrock. The seawall profile is defined as a vertical plane extending from −6 m to +9 m, although it is most likely an inclined plane (the precise morphology of the seawall below the sediment is unknown). Thus, only the section of the profile between −6 and +9 m in elevation is composed of material that can be mobilized by waves, although sediment could theoretically be deposited and then removed outside these limits, depending on hydrodynamic conditions.

A sensitivity analysis on permeability $k$ involves running 10 to 13 model scenarios for each modeled day. Three metrics—the root mean square error (RMSE, Equation (5)) between measured and modeled profiles, the coefficient of determination (R², Equation (6)) on the relative change in elevation, and the Brier Skill Score ($BSS$, Equation (7))—are computed to evaluate model performance, focusing on the intertidal zone where VMS-derived elevation data are available.

$$RMSE=\sqrt{{\left(z^{\prime }-z\right)}^2}$$

(5)

where $z$ is the observed elevation along the profile, and $z\prime$ is the modeled one.

$$R=\frac{Cov\left(X,Y\right)}{{\sigma }_X\times {\sigma }_Y}$$

(6)

where $Cov\left(X,Y\right)$ is the covariance between two vectors $X$ and $Y$, as calculated by $Cov\left(X,Y\right)=E\left[\left(X-\overline{X}\right)\times \left(Y-\overline{Y}\right)\right]$.

The equation of $\mathrm{B}\mathrm{S}\mathrm{S}$ considers the initial observed and both final observed and modeled profiles. The main advantage of this index is that it only takes into account the active part of the profile [45], thus limiting potential bias such as the presence of a non-erodible profile section. $BSS$ values are generally considered good for values superior to 0.6;

Table 2. BSS interpretation table according to Van Rijn et al. [46].

BSS Range	Interpretation
BSS ≤ 0	Bad
0 < BSS ≤ 0.3	Poor
0.3 < BSS ≤ 0.6	Fair
0.6 < BSS ≤ 0.8	Good
0.8 > BSS	Excellent

provides an interpretation chart based on the work of Van Rijn et al. [46].

$$BSS=1-\frac{\frac{1}{N}\sum _{i=1}^N{\left(z_{final, i}-{z\prime }_{final, i}\right)}^2}{\frac{1}{N}\sum _{i=1}^N{\left(z_{initial, i}-{z\prime }_{final, i}\right)}^2}$$

(7)

The cross-shore profile starts offshore at a depth of 24.6 m and extends 3490 m up to the seawall at +9 m. Resolution is gradually refined from 3 m offshore to 1 m near the seawall to capture wave processes adequately. The selected transect minimizes alongshore transport effects, which are not considered by XBeach-G, due to its location around the beach cell rotation pivot center, despite significant transport demonstrated by shoreline rotation mechanisms [4].

4. Results

4.1. XBeach-G Storm Simulations

The morphological response of the intertidal topography was numerically investigated for three distinct storm events (A, B, and C), each characterized by different severity and duration. The validity of these investigations was established through video observations, utilized for calibrating the initial profile extracted before each event.

Numerical simulations of the selected beach profile were conducted for thirteen permeability values, ranging from 0 to 0.6, during each storm event (A, B, and C). For each storm, XBeach-G simulations were performed every 24 h (twice for storm A over 35 h, 3 times for storm B over 42 h, and 5 times for storm C over 58 h), resulting in a total of 26, 39, and 56 numerical simulations for storms A, B, and C, respectively. Selected simulations are illustrated in Figure 7, Figure A2 and Figure A3. Subsequently, error calculations were performed for each simulation, encompassing varying permeability values and time intervals corresponding to the progression of storms, with assessments conducted every 24 h.

Table 3. Validation results of the XBeach-G simulations of storm A. Metrics include the Brier Skill Score (BSS), the root mean square error (RMSE) and the coefficient of determination (R²). Red values in bold highlight scenarios with BSS > 0.6, RMSE < 0.25 m, or R² > 0.9. Yellow values in italic indicate BSS > 0.3, RMSE < 0.5 or R² > 0.5.

Parameter	Dates	$\mathbf{Permeability} \mathit{k}$ (m/s)
		0	0.01	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5	0.6
BSS	06 Dec 2018	−36.19	0.43	0.50	0.56	0.60	0.63	0.63	0.63	0.63	0.60	0.57	0.54	0.49
	07 Dec 2018	−27.65	0.68	0.69	0.68	0.63	0.58	0.52	0.46	0.38	0.31	0.25	0.16	0.02
RMSE (m)	06 Dec 2018	2.86	0.34	0.32	0.30	0.28	0.27	0.26	0.26	0.27	0.28	0.28	0.29	0.31
	07 Dec 2018	2.43	0.28	0.26	0.26	0.30	0.32	0.34	0.37	0.39	0.42	0.44	0.47	0.51
R2	06 Dec 2018	0.21	0.56	0.63	0.69	0.71	0.73	0.72	0.71	0.70	0.66	0.62	0.59	0.53
	07 Dec 2018	0.72	0.75	0.73	0.74	0.73	0.72	0.70	0.69	0.64	0.56	0.47	0.32	0.12

presents a summary of BSS, RMSE, and R² results for various permeability scenarios during storm A. Further results of storms B and C are presented in

Table A1. Validation results of the XBeach-G simulations of both storms B and C. Metrics include the Brier Skill Score (BSS), the root mean square error (RMSE) and the coefficient of determination (R²). Red values in bold highlight scenarios with BSS > 0.6, RMSE < 0.25 m, or R² > 0.9. Yellow values in italic indicate BSS > 0.3, RMSE < 0.5 or R² > 0.5.

Parameter	Storm	Dates	k (m/s)
			0	0.01	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5	0.6
BSS	B	12 Dec 2019	−12.54	−0.53	−0.48	−0.48	−0.50	−0.48	−0.50	−0.49	−0.50	−0.48	−0.48	−0.47	−0.44
		13 Dec 2019	−76.51	−0.68	−0.53	−0.41	−0.31	−0.20	−0.13	−0.04	0.05	0.12	0.21	0.22	0.33
		14 Dec 2019	−12.68	−3.03	−2.75	−2.55	−2.42	−2.29	−2.20	−2.13	−2.04	−1.93	−1.87	−1.78	−1.55
	C	04 Dec 2019	-	−0.05	−0.19	−0.33	−0.41	−0.49	−0.52	−0.52	−0.52	−0.51	-	-	−0.37
		05 Dec 2019	-	0.32	0.14	−0.07	−0.23	−0.37	−0.51	−0.59	−0.66	−0.64	-	-	−0.51
		06 Dec 2019	-	−1.69	−2.12	−2.52	−2.85	−3.17	−3.45	−3.74	−3.95	−4.10	-	-	−4.11
		07 Dec 2019	-	−0.30	−0.41	−0.49	−0.55	−0.58	−0.59	−0.60	−0.56	−0.55	-	-	−0.41
		08 Dec 2019	-	0.29	0.08	−0.14	−0.40	−0.67	−0.94	−1.21	−1.47	−1.72	-	-	−2.56
RMSE	B	12 Dec 2019	2.74	0.90	0.87	0.87	0.88	0.89	0.90	0.90	0.91	0.90	0.91	0.90	0.90
		13 Dec 2019	5.65	0.78	0.75	0.73	0.70	0.68	0.65	0.62	0.59	0.56	0.53	0.53	0.49
		14 Dec 2019	1.40	0.81	0.80	0.78	0.77	0.76	0.75	0.74	0.72	0.71	0.70	0.68	0.64
	C	04 Dec 2019	-	0.22	0.23	0.24	0.24	0.25	0.25	0.25	0.24	0.24	-	-	0.23
		05 Dec 2019	-	0.11	0.13	0.14	0.15	0.15	0.16	0.17	0.17	0.17	-	-	0.16
		06 Dec 2019	-	0.40	0.43	0.46	0.48	0.49	0.51	0.52	0.53	0.53	-	-	0.52
		07 Dec 2019	-	0.35	0.37	0.38	0.39	0.40	0.39	0.39	0.38	0.38	-	-	0.34
		08 Dec 2019	-	0.44	0.48	0.53	0.59	0.65	0.70	0.75	0.79	0.83	-	-	0.97
R2	B	12 Dec 2019	0.68	0.24	0.18	0.14	0.10	0.07	0.05	0.03	0.02	0.00	0.00	0.00	0.00
		13 Dec 2019	0.24	0.02	0.02	0.01	0.01	0.01	0.00	0.00	0.00	0.01	0.03	0.03	0.07
		14 Dec 2019	0.01	0.09	0.07	0.06	0.05	0.04	0.04	0.03	0.03	0.02	0.01	0.01	0.01
	C	04 Dec 2019	-	0.02	0.13	0.19	0.19	0.20	0.20	0.19	0.18	0.17	-	-	0.11
		05 Dec 2019	-	0.57	0.59	0.59	0.59	0.59	0.57	0.55	0.53	0.52	-	-	0.42
		06 Dec 2019	-	0.11	0.13	0.13	0.14	0.14	0.15	0.17	0.19	0.19	-	-	0.17
		07 Dec 2019	-	0.47	0.55	0.60	0.64	0.62	0.58	0.53	0.48	0.45	-	-	0.19
		08 Dec 2019	-	0.56	0.45	0.35	0.25	0.18	0.11	0.07	0.04	0.02	-	-	0.00

. Among the modeled dates, only four achieved BSS values higher than 0.3 (threshold of the “fair” category Van Rijn et al. [46]), with two exceeding 0.6 (“good” category). The highest BSS values occurred during the most intense storm (storm A), particularly during the peak between 07 Dec 2018 and 08 Dec 2018 ($E_{flux max}$ = 7.1 × 10⁴ W/m/s). However, storm B, with a higher-intensity peak ($E_{flux max}$ = 7.5 × 10⁴ W/m/s), did not yield BSS values above 0.3, indicating no clear correlation between BSS and wave condition intensity. The lowest BSS values were associated with permeability value $k$ = 0, equivalent to the deactivation of the groundwater flow process.

BSS values covaried with permeability, suggesting that the maximum BSS value for one date indicates the optimal permeability value. Two consecutive dates with BSS exceeding 0.6 were identified: on 16 Dec 2018, optimal $k$ was between 0.2 and 0.35 m/s under relatively calm conditions ($E_{flux max}$ = 5 × 10³ W/m/s), and on 07 Dec 2018, optimal $k$ ranged between 0.01 and 0.15 m/s (maximum for $k$ = 0.01 m/s) during a storm peak ($E_{flux max}$ = 5 × 10⁴ W/m/s). Both dates exhibited a minimum RMSE of 0.26 m and maximum R² of 0.73 and 0.75, respectively.

The two dates with BSS classified as “fair” (13 Dec 2019, BSS = 0.33, and 05 Dec 2019, BSS = 0.32) remained in the lower limit (0.3) of the fair class. These dates corresponded to storm conditions ($E_{flux max}$ = 7 × 10⁴ W/m/s) and calm conditions ($E_{flux max}$ = 1 × 10² W/m/s), respectively, with associated RMSE values of 0.49 m and 0.11 m, and R² of 0.07 and 0.57. Results from the XBeach-G model were more significant for the storm A compared to storms B and C (Table 3 and Table A1).

Detailed analysis of the intertidal profile changes on 06 Dec 2019 revealed the accurate reproduction of accretion phenomena at the top of the beach for various permeability values. However, simulated movements in the lower part of the beach were underestimated, particularly with lower permeability values, excluding $k$ = 0 m/s. On 07 Dec 2019, where hydrodynamic conditions were stronger, the model accurately reproduced accretion at the bottom of the beach for different $k$ values, indicating lower sensitivity to permeability. However, erosion of the upper beach was only reproduced for low $k$ values. Notably, for permeability values exceeding 0.3 m/s, the model predicted accretion of the upper beach instead of erosion, highlighting greater sensitivity to permeability in this region of the profile. The scenario with $k$ = 0 m/s significantly overestimated both lower-beach accretion and upper-beach erosion. Modeled movements outside the area monitored by VMS remained relatively small.

4.2. Morphological Observations

Three indices of morphological variability on the beach in response to the different storm events were extracted from the VMS-derived DEMs (Figure 8). The beach width (BW) was estimated as the average horizontal distance separating the DEM’s contour lines of various elevations from −2 m to 3 m. The beach slope (BS) was computed as the ratio between the tidal range and BW. The beach orientation is calculated by approximating the shoreline to a parabola. The orientation angle (BOA) is then calculated as the angle between the seawall and the parabola’s tangent at the transect location (Figure 3b). The results are presented in Figure 8 for storm A, and Figure A2 and Figure A3 for storms B and C, respectively.

5. Discussion

5.1. Influence of Permeability on Morphodynamics

During storm A, the results highlight the sensitivity of pebble beach morphodynamics to permeability. The absence of groundwater flow ($k$ = 0 m/s) exaggerates sediment erosion and accumulation phenomena significantly, emphasizing the crucial role of considering groundwater flow in estimating these dynamics. This is consistent with observations from Buscombe and Masselink [43], Jennings and Shulmeister [27], and Mason and Coates [47], and was one of the main incentives for the adaptation of XBeach for gravel environments [7]. This sensitivity is explained by the importance of water seepage in the rapid loss of wave energy during uprush on gravel beaches.

Generally, the results show that for non-zero values, $k$ controls the amplitude of the profile changes, while the position along the profile of the maximum and minimum changes depends on the antecedent topography and hydrodynamic conditions (tidal amplitude, wave characteristics). Lower $k$ values in calm conditions tend to smoothen the profile, while higher $k$ values favor berm buildup in the upper intertidal zone. Under energetic conditions, low permeability fosters profile erosion and accumulation at the beach step, while high values tend to build a berm by eroding the beach step. This aligns with findings by She et al. [48], who observed increased offshore sediment transport and reduced onshore transport with decreasing permeability in a mixed sand and gravel beach physical model.

5.2. Limitations of the Modelling Approach

The low number of models meeting the BSS ≥ 0.6 condition underscores the limitations of the approach. Challenges include the reduced quality of VMS images in winter conditions, especially during storms, which affects the reconstruction of DEMs in terms of number, extension, and potentially quality. Additionally, the model updates only the intertidal section of the profile visible to cameras from one storm to the next, neglecting potential subtidal changes and features like the beach step, known for its importance in gravel beach dynamics [43,49,50]. Similarly, supratidal morphology, constituting a sedimentary stock accessible to storm overtopping waves, is not updated. Despite these limitations, the model results for storm 2 align well with observed morphological changes, suggesting relative closeness to reality for that specific date.

Other studies at different sites reveal similar model performances. For example, Pollard [51] modeled the storm response of a gravel barrier in the UK, where the D50 varied between 1 and 2 cm. Forty-five calibration runs, including permeability variations, yielded a maximum BSS of 0.71 for $k$ = 0.2 m/s. The model was then applied to seven different profiles, among which only two resulted in BSS values greater than 0.6. Bergillos et al. [22] applied XBeach-G to a mixed sand and gravel deltaic coastline with a D50 of 2 cm. Permeability was not reported by the authors, but their results show BSS values greater than 0.89. Brown et al. [23] tested the model on a mixed sand and gravel beach in New Zealand, with a D50 of 1.36 cm and a permeability $k$ = 0.01 m/s. Only 7 of the 30 scenarios obtained a BSS greater than 0.6, with a maximum of 0.79.

5.3. Implications of Grain Size for Modelling Pebble Beaches

The observed morphological changes during the calm period involve the construction of a horizontal berm, while the initial peak of storm A results in berm erosion and profile smoothing. These behaviors align with findings by McCall et al. [20] and correspond to the rollover mechanism [4]. McCall et al. [20] defined rollover as the cross-shore landward displacement of the entire profile of a barrier without a change in height, width, or slope. Although the beach of Etretat is not precisely a barrier due to its fixed seawall constraint, the rollover manifests as an oscillation of the slope around a pivot axis in the system at the x = 3450 m coordinate and between 0 and +1 m in elevation (Figure 7 and Figure 8a,c).

Furthermore, the position along the profile and the width of the horizontal berms vary from one date to another and appear to be correlated with the permeability value and the intensity of hydrodynamic conditions. Specifically, for any given date, a larger $k$ value tends to shape the berm, while a lower $k$ value is likely to erase it and smoothen the profile shape (Figure 7, Figure A2 and Figure A3). The study of the morphological characteristics of berms as a proxy for permeability during storms emerges as a relevant research axis for future investigations.

Considering a well-sorted sediment, a larger sediment size implies a grater porosity and increased seepage velocities under swash (i.e., an enhanced ability to dissipate wave energy) [43]. The observed significant granulometric sorting on gravel beaches with strong spatial and temporal variability [4,43], as well as the mobilization of a substantial amount of sediment during storm events [52], suggests that the permeability of a beach may vary on a daily scale. In the case of a gravel beach, it is reasonable to expect greater permeability during sustained calm conditions over a longer period, allowing finer sediments to migrate through the porosity. Conversely, under more energetic conditions, the sudden mobilization of sediment can mix grain sizes and potentially decrease permeability.

Indeed, studies on mixed gravel and sand beaches highlight that the presence of sand plays a crucial role in influencing the response of the beach profile to wave impact. Permeability, in such cases, is controlled by the finest 10% of the sediment [53,54]. Mason et al. [55] demonstrated that when the surface sedimentary layer (depth < 1 m) contains more than 25% of its weight in sand, the beach morphodynamics become more similar to that of sandy beaches than mixed ones.

Although Etretat is primarily considered a pebble, there exists a significant amount of underlying sand in the system. In certain instances, after storms, this sand can become the predominant substrate of the beach. While this occurrence is rare, it demonstrates that finer materials, such as sand, could fill up the pebble’s porosity and potentially reduce the permeability of the beach.

Under this hypothesis, the impact of successive storms of moderate intensity in a cluster could be greater than stronger but isolated events. This is because clustered storms would occur on a beach with a permeability lowered by antecedent storms, hence with a reduced ability to absorb wave energy through infiltration/exfiltration processes.

Storm A was preceded by a month of calm conditions, with the previous storm occurring on November 8. Before the first peak of storm A, the model results indicate that the permeability optimum was between 0.20 and 0.35 m/s in calm conditions, later reduced to around 0.05 m/s under more energetic conditions. These values approximately correspond to well-sorted coarse and fine gravel sizes, respectively [43,56]. These results suggest strong and rapid temporal variability in permeability in pebble systems, with a decrease by a factor of 4 at least within a day.

The strict character, “purely gravel”, of the beach (the word gravel is here referring to particles coarser than 2 mm, i.e., sand) could shift occasionally to a “mixed sand and gravel”-type system. Considering that storms are seasonal events, the transition from pure gravel to mixed sand and gravel behavior (i.e., not necessarily resulting in surficial sand) could be seasonal as well. An analogous seasonal transition between beach types was documented by Casamayor et al. [56] on a composite beach. The authors found that their composite study site behaved as a pure gravel system in winter and as a composite system during summer due to wave climate variability and associated morphological changes. They proposed adding a seasonal component to the composite beach type of Jennings and Shulmeister [27], and the results here tend to support that it could be a relevant addition for the pure gravel type as well.

Moreover, the low BSS results obtained with models of storms B and C could be partly due to the beach behaving more like a sandy one. In such cases, the sediment transport equations used by XBeach-G could become obsolete, highlighting the need for a more nuanced understanding of the beach composition and dynamics in different seasonal and storm contexts.

This analysis provides valuable insights into the complex interplay between storm dynamics, permeability, and morphological changes in pebble beaches, highlighting the need for further exploration in this field.

5.4. Insights on the Complex Evolution of Pebble Beaches

In summary, it was possible to successfully model a portion of storm A in Etretat, despite the inherent limitations of the dataset. However, the modeling capabilities remain limited by the previously mentioned data accuracy limitations and by the conceptual simplifications of Xbeach-G. The reliance on 1D model computations leads to ignoring potentially important longshore exchanges (Soloy et al., [4], especially in the case of an oblique storm swell such as in Etretat. Notably, storm A stood out as the only modeled event where BOA remained consistently aligned at all elevations (from −2 to −3 degrees), with only 1 degree of clockwise rotation during the modeled period. Conversely, other storms exhibited broader variations in shoreline orientations (Figure A4 and Figure A5), pointing to localized alongshore processes introducing significant complexity not accounted for in Xbeach-G.

The findings also underscore notable daily variability in permeability, aligning with recent observations on the spatiotemporal variability in the surface grain size of pebble beaches. However, k is imposed as a spatially and temporally uniform constant, due to the inherent challenges in accurately measuring this parameter, which may limit model convergence and, by extension, the ability to model gravel coasts. Moreover, although gravel size strongly influences permeability and is known to vary over time and space (Horn and Walton, [54]), it is described in Xbeach as spatially and temporally uniform, relying on single values of D50 and D90. Advancements in coastal monitoring and grain size mapping techniques proposed in recent studies [17,26] hold the promise of providing high-resolution observational datasets. Utilizing such data could significantly enhance models, enabling a more accurate consideration of the spatiotemporal variability in permeability during storm events.

6. Conclusions

This study employed Xbeach-G to simulate the impact of three intense storms occurring between 2018 and 2020 on the Etretat beach in Normandy, France. The modeling aimed to calibrate the permeability value (k) to understand its day-to-day and inter-storm variations. Among the modeled events, satisfactory results were obtained for storm A, with maximum BSS values of 0.63 and 0.69 on the first and second modeled days, respectively. These results indicated permeability values of k = 0.20 m/s to 0.35 m/s and k = 0.05 m/s, consistent with expectations for well-sorted coarse and fine gravel sizes, respectively.

However, storms B and C yielded unsatisfactory results, likely due to limitations in the approach, including imprecise subtidal and supratidal topography, lack of updating in the subtidal area, and the neglect of longshore sediment transport. Indeed, the time series of beach orientation angle suggested significant longshore transport, which the model could not capture due to the one-dimensional nature of Xbeach-G.

This study emphasized the sensitivity of the beach profile to permeability, where high values led to berm formation, and low values resulted in erosion. The intensity of wave conditions influenced the amplitude of topographic changes. The storm A results indicated a potential fourfold decrease in permeability within 24 h, emphasizing the dynamic nature of pebble beaches during storm events, especially regarding their ability to dissipate wave energy through infiltration/exfiltration processes.

The findings challenged the strict characterization of Etretat’s beach as “purely gravel”, suggesting occasional transitions to a “mixed sand and gravel” system, possibly seasonal. The spatial and temporal uniformity of k, grain size, and porosity were identified as limitations in the Xbeach-G model. Future developments could involve incorporating granulometric observations and addressing these limitations for a more accurate representation of gravelly systems. This study underscores the complexity of pebble beaches, advocating for innovative techniques, model improvements, and integrative approaches, combining numerical models, physical experiments, and field observations to enhance the understanding of morphological responses and uncertainties.