Jean-Baptiste Clément
email: jean-baptiste.clement@umontpellier.fr

Damien Sous
email: damien.sous@mio.osupytheas.fr

Frédéric Bouchette
email: frederic.bouchette@umontpellier.fr

Frédéric Golay
email: frederic.golay@univ-tln.fr

Mehmet Ersoy

A Richards' equation-based model for wave-resolving simulation of variably-saturated beach groundwater flow dynamics

Keywords: Swash zone, Barrier beach, Numerical modelling, Unsaturated porous media, Seepage, Adaptive discontinuous Galerkin method

This study introduces a model based on Richards' equation to describe variably-saturated beach groundwater flow. The surface wave propagation is computed by the phase-resolving non-hydrostatic SWASH code. The SWASH data are used to make a suitable dynamic boundary condition at the beach face to force Richards' equation. The latter is solved by a weighted discontinuous Galerkin method together with adaptive mesh refinement. The model is validated by comparison with a laboratory experiment of a transient water table recharge problem. Then, the BARDEX II prototype-scale experiment is considered to assess the model abilities for beach groundwater dynamics. The barrier beach is studied for three cases with different lagoon levels. Steady-state results with no-wave conditions show excellent agreement. Transient waves simulations are evaluated in terms of pressure heads, saturations, water table position and groundwater velocities for timeaveraged, swash-resolving and spectral analysis. Results bring interesting insights about beach groundwater modelling by comparison with the experimental data as well as a Darcy's equation-based model. A first investigation is carried out to assess the groundwater effect on the bed sediment dynamics through the modification of sediment relative weight.

Introduction

Beach groundwater dynamics has received a growing attention during the last two decades due to its importance in a number of physical, biological and chemical processes. The determination of groundwater transport through beaches or coastal barriers remains a key challenge in the prediction of salinity, contaminants and nutrients exchanges between the open sea, the lagoon and the aquifer which in turn control crucial biogeochemical cycles, see e.g. [START_REF] Robinson | Driving mechanisms for groundwater flow and salt transport in a subterranean estuary[END_REF][START_REF] Burnett | Quantifying submarine groundwater discharge in the coastal zone via multiple methods[END_REF][START_REF] Anschutz | Tidal sands as biogeochemical reactors[END_REF][START_REF] Sawyer | Stratigraphic controls on fluid and solute fluxes across the sediment-water interface of an estuary[END_REF][START_REF] Charbonnier | Aerobic respiration in a high-energy sandy beach[END_REF][START_REF] Mcallister | Dynamic hydrologic and biogeochemical processes drive microbially enhanced iron and sulfur cycling within the intertidal mixing zone of a beach aquifer[END_REF][START_REF] Heiss | Physical Controls on Biogeochemical Processes in Intertidal Zones of Beach Aquifers[END_REF]. The beach groundwater circulation is mainly controlled by the combination of the cross-barrier gradients due to tides or large-scale fluctuations of the mean water level and the wave forcing at the shoreline, including swell and long infragravity waves [START_REF] Nielsen | Tidal dynamics of the water table in beaches[END_REF][START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF][START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF][START_REF] Sous | Groundwater fluxes within sandy beaches: swash-driven flow vs cross-barrier gradients[END_REF][START_REF] Bertin | Infragravity waves: From driving mechanisms to impacts[END_REF]. The swash zone is the ultimate part of the beach exposed to ocean forcing, alternatively covered and uncovered by wave run-up (uprush) and run-down (backwash). It moves across the beachface following the fluctuations of the mean water level at the shoreline due to tides, atmospheric pressure variations, wind or wave setup. The swash zone hydrodynamics involves various processes characterized by wave dynamics, turbulence, sediment transport and interactions with the underlying porous bed. The percolation flows have early been suspected to affect swash zone sediment transport and beach morphodynamics. The first effect is the loss of volume of the swash tongue, only active on coarse-grained beach. The two other processes involved in sediment transport modification by in/exfiltration flows through the beach face are, on one hand, the modification of the boundary layer shear stress (the "ventilated bed" effect [START_REF] Chen | Velocity Distribution of Turbulent Open-Channel Flow with Bed Suction[END_REF][START_REF] Conley | Ventilated oscillatory boundary layers[END_REF][START_REF] Corvaro | Flow dynamics on a porous medium[END_REF][START_REF] Corvaro | Fluid-particle interaction and generation of coherent structures over permeable beds: an experimental analysis[END_REF]) and, on the other hand, the change of the relative weight of sediment. Significant research efforts have been engaged to estimate and model the net effect of these processes on sediment dynamics [START_REF] Turner | Swash infiltration-exfiltration and sediment transport[END_REF][START_REF] Butt | The influence of swash infiltration-exfiltration on beach face sediment transport: onshore or offshore?[END_REF][START_REF] Li | Wave-induced beach groundwater flow[END_REF][START_REF] Masselink | The role of swash infiltration in determining the beachface gradient: a numerical study[END_REF][START_REF] Li | BeachWin: modelling groundwater effects on swash sediment transport and beach profile changes[END_REF][START_REF] Karambas | Modelling of Infiltration-Exfiltration Effects of Cross-Shore Sediment Transport in the Swash Zone[END_REF][START_REF] Austin | Swash-groundwater interaction on a steep gravel beach[END_REF][START_REF] Bakhtyar | Wave-induced water table fluctuations, sediment transport and beach profile change: Modeling and comparison with large-scale laboratory experiments[END_REF][START_REF] Masselink | Large-scale laboratory investigation into the effect of varying back-barrier lagoon water levels on gravel beach morphology and swash zone sediment transport[END_REF][START_REF] Chardón-Maldonado | Advances in swash-zone research: Smallscale hydrodynamic and sediment transport processes[END_REF][START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF].

Overall, while the impact is major for coarse-grained beaches, it is now generally acknowledged that the ex/infiltration effects on sediment transport for fine sand beach should be overpowered by the other physical processes acting in the swash zone morphodynamics [START_REF] Sous | Swash-groundwater dynamics in a sandy beach laboratory experiment[END_REF][START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF]. This is in line with the mixed success met by beach drainage systems in trying to alter beach morphodynamics by controlling groundwater head [START_REF] Horn | Measurements and modelling of beach groundwater flow in the swash-zone: a review[END_REF][START_REF] Turner | Beach Dewatering as a 'Soft' Engineering Solution to Coastal Erosion-A History and Critical Review[END_REF].

Owing to the difficulty to measure experimentally the groundwater parameters, and in particular to monitor their spatio-temporal dynamics over a relevant area, a series of numerical approaches have been proposed during the last few decades. First models based on Darcy's theory and purely horizontal groundwater dynamics (Dupuit-Forchheimer assumption) [START_REF] Nielsen | Tidal dynamics of the water table in beaches[END_REF][START_REF] Kang | Watertable Dynamics in Coastal Areas[END_REF] allowed to describe the tidal dynamics. Similar assumptions of purely horizontal flows were made in the coupled approach using Boussinesq model for wave propagation and Darcy-Forchheimer within the porous medium [START_REF] Karambas | Modelling of Infiltration-Exfiltration Effects of Cross-Shore Sediment Transport in the Swash Zone[END_REF][START_REF] Karambas | Prediction of sediment transport in the swash-zone by using a nonlinear wave model[END_REF]. Finer vertically-resolving model have been proposed based on the Darcy-related 2D Laplace law to model groundwater dynamics forced by tides and waves [START_REF] Li | Numerical modelling of tide-induced beach water table fluctuations[END_REF][START_REF] Li | Wave-induced beach groundwater flow[END_REF] and further refined to describe in/exfiltration bed exchanges [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF]. The groundwater circulation cell highlighted by coupled VOF-RANS model [START_REF] Bakhtyar | Wave-induced water table fluctuations, sediment transport and beach profile change: Modeling and comparison with large-scale laboratory experiments[END_REF] has been later confirmed by laboratory [START_REF] Sous | Swash-groundwater dynamics in a sandy beach laboratory experiment[END_REF] and field [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF] experiments. In coarse-grained sand or gravel beaches, the Darcy theory is still used but extended beyond its laminar theoretical framework with the Darcy-Forchheimer approximation [START_REF] Van Gent | The modelling of wave action on and in coastal structures[END_REF][START_REF] Steenhauer | Numerical model of swash motion and air entrapment within coarse-grained beaches[END_REF][START_REF] Desombre | SIMULTANEOUS SURFACE AND SUBSURFACE AIR AND WATER FLOWS MODELLING IN THE SWASH ZONE[END_REF]. Quasi-3D models have been proposed, relying on an a priori parabolic profile of the non-hydrostatic groundwater pressure field [START_REF] Mccall | MODELLING OVER-WASH AND INFILTRATION ON GRAVEL BARRIERS[END_REF][START_REF] Mccall | Modelling storm hydrodynamics on gravel beaches with XBeach-G[END_REF].

Although significant effort has been deployed during the last two decades in the understanding and the modelling of beach groundwater dynamics, most of the studies have been (i) based on Darcy-Forchheimer approximation, i.e. totally ignoring the partial saturation and capillarity effect which are essential in finegrained sand beaches [START_REF] Heiss | Coupled surface-subsurface hydrologic measurements reveal infiltration, recharge, and discharge dynamics across the swash zone of a sandy beach[END_REF] and (ii) mainly focused on horizontal dynamics, although detailed field measurements demonstrated the presence of vertical head gradients and related circulation cells in sandy beaches [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF].

Aiming to tackle these two major challenges in the fine sand beach context, the present study aims to model wave-resolved variably-saturated groundwater flow based on Richards' equation. Generalized Richards' equation are generally used to model surface/subsurface flow in hydrological large-scale context (see e.g. [START_REF] Weill | A generalized Richards equation for surface/subsurface flow modelling[END_REF][START_REF] Pachepsky | Generalized Richards' equation to simulate water transport in unsaturated soils[END_REF]), but very few attempts have been made to implement Richards' equation in wave-forced environments [START_REF] Caviedes-Voullième | Numerical simulation of groundwater-surface interactions by external coupling of the 3D Richards equation and the full 2D shallow-water equations[END_REF].

Richards' equation is a nonlinear degenerate parabolic equation, since terms can vanish in the space and time derivatives, whose solving remains a challenge because of the nonlinearities and the degeneracies of hydraulic properties [START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF][START_REF] Zha | Review of numerical solution of Richardson-Richards equation for variably saturated flow in soils[END_REF]. In this regard, beach groundwater dynamics is expected to be a demanding problem because it holds a wide range of space/time scales and is forced by the fast dynamics of waves.

To solve efficiently Richards' equation, the work presented here uses a strategy based on a Discontinuous Galerkin method in combination with adaptive mesh refinement. This strategy has been introduced in Clément et al. [START_REF] Clément | An adaptive strategy for discontinuous Galerkin simulations of Richards' equation: Application to multi-materials dam wetting[END_REF] and successfully applied to a multi-materials dam wetting problem. Extensive test-cases are also available in [START_REF] Clément | Numerical simulation of flows in unsaturated porous media by an adaptive discontinuous Galerkin method: application to sandy beaches[END_REF]. The resulting code is called Rivage.

The first part of the manuscript introduces the groundwater flow model and the numerical strategy. The second part presents an experimental test-case to validate the numerical model and provide an overview of groundwater flows. Then, the numerical model is used to make a detailed analysis of a beach groundwater laboratory experiment. A discussion is handled to assess beach groundwater dynamics with Richards' equation. The last section is dedicated to conclusions and prospects.

The groundwater flow model

Richards' equation

Richards' equation [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF][START_REF] Richards | Capillary conduction of liquids through porous mediums[END_REF] is a degenerate nonlinear parabolic equation used to describe flow in variablysaturated porous media. It focuses on the water phase dynamics (the air phase is supposed to be continuouslyconnected with the atmosphere) and neglects inertial effects. In the present study, the mixed formulation of Richards' equation is chosen because it is mass conservative and can be applied to saturated/unsaturated and heterogeneous porous media [START_REF] Celia | A general mass-conservative numerical solution for the unsaturated flow equation[END_REF][START_REF] Farthing | Numerical Solution of Richards' Equation: A Review of Advances and Challenges[END_REF]. It is written here in terms of hydraulic head h = ψ + z [L]:

∂ t θ(h -z) -∇ • (K(h -z)∇h) = 0, ( 1
)
where ψ is the pressure head [L], θ the water content [-], K the hydraulic conductivity [L•T -1 ] and z the elevation [L]. The water table (ψ = 0) splits the porous medium between the saturated (ψ > 0) and the unsaturated zones (ψ < 0). Above the water table, the capillary fringe is here set wherever the water content is at least 95% of the saturated water content.

Solving Richards' equation requires two constitutive laws: one for water content and one for hydraulic conductivity. The water content is described in terms of effective saturation S e [-] and the tensor of hydraulic conductivity is supposed to react to pressure head identically for each space direction. This leads to write:

θ(ψ) = θ r + (θ s -θ r )S e (ψ), (2)
K(ψ) = K s K r (ψ), (3)
where θ s denotes the saturated water content [-], θ r the residual water content [-], K s the saturated hydraulic conductivity tensor [L • T -1 ] and K r the relative hydraulic conductivity [-]. The hydraulic properties present two different behaviours depending on whether the porous media is saturated (ψ ≥ 0) or not (ψ < 0):

S e (ψ) =      1 if ψ ≥ 0, S ⋆ e (ψ) otherwise, and K r (ψ) =      1 if ψ ≥ 0, K ⋆ r (ψ) otherwise. (4)
S ⋆ e and K ⋆ r are monotonic increasing functions of pressure head in the unsaturated zone. Relations used in the paper are compiled in Tab. 1.

Name Expression Parameters

Vachaud's relations

S ⋆ e (ψ) = (1 + (α|ψ|) n ) -m K ⋆ r (ψ) = S ⋆ e 0.5 (ψ) 1 -1 -S ⋆ e 1 m (ψ) m 2 α: parameter linked to air entry pressure inverse [L -1 ] n > 1: pore-size distribution [-] m = 1 - 1 n : pore-size distribution [-]
Table 1: Hydraulic relations used in this thesis.

Boundary conditions

Richards' equation is used with classical boundary conditions like Dirichlet boundary conditions for the hydraulic head and Neumann boundary conditions for the flux. But the interaction of the porous medium (beach) with the open surrounding medium (atmosphere and free surface water) has to be modelled by particular boundary conditions described in the sections below.

Seepage

The seepage boundary condition is used to model the interface between a porous medium and the atmosphere imitating an outflow condition [START_REF] Schweizer | Regularization of outflow problems in unsaturated porous media with dry regions[END_REF][START_REF] Scudeler | Examination of the seepage face boundary condition in subsurface and coupled surface/subsurface hydrological models[END_REF]. If the porous medium is saturated and outflow occurs, then water pours out at atmospheric pressure. Otherwise, the interface is an impervious boundary, i.e. there is no flux. This condition can be written under different forms but also simply viewed as a switch between Dirichlet and Neumann boundary conditions. Then, the seepage boundary condition may be interpreted as a nonlinear Robin boundary condition. In this study, it is treated inside the nonlinear iterative process according to the previous solution guess at a local level of the numerical scheme, see [START_REF] Clément | An adaptive strategy for discontinuous Galerkin simulations of Richards' equation: Application to multi-materials dam wetting[END_REF] for further insights.

This approach is very general since it does not need any assumption about the seepage face.

Dynamic waves forcing

Surface and subsurface flows are a unified system by nature. For technical reasons (space/time scales, computational solvability), they are often divided into two sub-models which work together through interface conditions. Numerically, different coupling schemes are possible, see the review of Furman [START_REF] Furman | Modeling Coupled Surface-Subsurface Flow Processes: A Review[END_REF] or the work of Caviedes-Voullième et al. [START_REF] Caviedes-Voullième | Numerical simulation of groundwater-surface interactions by external coupling of the 3D Richards equation and the full 2D shallow-water equations[END_REF]. The present study focuses on the groundwater dynamics in fine sand conditions and assumes that the infiltrated volume loss, pore velocities and boundary layer modification are weak enough to not affect the swash dynamics and therefore to neglect the groundwater feedback on the surface flow. The model coupling is therefore a pure one-way forcing: it only requires input data based on the solution from the surface water sub-model to prescribe boundary conditions for the groundwater sub-model.

To this end, a dynamic forcing boundary condition is built as a combination of a Dirichlet boundary and a seepage boundary condition. The former is applied on the flooded beach face while the latter is applied on the remaining part, i.e. beyond the swash tip. In this way, the model allows the swash tip and the water table to be disconnected. The forcing is one-way: the seepage water is supposed to not affect the surface flow.

The open-source SWASH code simulates the propagation of free surface waves according to selected input parameters at the offshore boundary. The governing equations are the nonlinear shallow water equations including non-hydrostatic pressure. The wave-resolving model is used here in a two-layer configuration, all parameters being kept by default [START_REF] Zijlema | SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[END_REF]. The bottom pressure computed by the wave model is supposed to transfer directly an instantaneous hydraulic head to the porous medium by the mean of the Dirichlet boundary condition. The latter is both space-and time-dependent, i.e. evolving with the waves elevation.

Technically, the data are interpolated linearly from the space-time grid of SWASH to the space-time grid of the groundwater sub-model. The position of the swash tip is treated specifically by extrapolation to correspond to the closest groundwater sub-model discretization.

Numerical strategy

The solution of Richards' equation is known to exhibit spurious oscillations, typically in the presence of sharp wetting fronts moving dynamically into initially dry soils. In addition, nonlinear convergence is difficult to achieve and can even fail. These drawbacks make Richards' equation costly to compute. In the case of wave-driven problems, rapid successive infiltrations are challenging to simulate. In order to alleviate these issues, a numerical strategy has been developed in [START_REF] Clément | Numerical simulation of flows in unsaturated porous media by an adaptive discontinuous Galerkin method: application to sandy beaches[END_REF] and presented in [START_REF] Clément | An adaptive strategy for discontinuous Galerkin simulations of Richards' equation: Application to multi-materials dam wetting[END_REF]. It is based on discontinuous Galerkin (DG) methods and adaptive mesh refinement. The main stages are summarized hereinafter.

Numerical methods

Space discretization of Richards' equation lies on DG methods for which extensive introduction can be found in Rivière [START_REF] Rivière | Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations[END_REF] or Dolejší and Feistauer [START_REF] Dolejší | Discontinuous Galerkin Method[END_REF]. In this work, the Incomplete Interior Penalty Galerkin (IIPG) under primal formulation is used because it keeps the formulation simple. For time discretization, the implicit Euler scheme is used. The implicit Euler scheme is simple and shows a wide region of stability suitable for the Richards' equation known to be stiff. The linearization of Richards' equation is crucial to obtain robust scheme for Richards' equation [START_REF] Paniconi | A comparison of Picard and Newton iteration in the numerical solution of multidimensional variably saturated flow problems[END_REF][START_REF] Lehmann | Comparison of Iterative Methods for Improved Solutions of the Fluid Flow Equation in Partially Saturated Porous Media[END_REF][START_REF] List | A study on iterative methods for solving Richards' equation[END_REF]. Here, the fixed-point method (also called Picard's iterations) is used. This is an iterative procedure based on an approximation of the Jacobian matrix from the Newton-Raphson method. Fixed-point iteration is only linearly convergent rather than quadratically as in the case for the Newton-Raphson method. On the other hand, the former is known to be more robust since it is not sensitive to the initial guess. To strengthen the solver robustness, the procedure is coupled with an adaptive time stepping [START_REF] Bergamaschi | Mixed finite elements and Newton-type linearizations for the solution of Richards' equation[END_REF][START_REF] Thoms | User's guide to the Variably Saturated Flow (VSF) process to MODFLOW[END_REF]. If it takes many iterations or fails to converge, then the time step is decreased. Inversely, if it converges in few iterations, then the time step can be increased.

Adaptive framework

The adaptive meshing strategy aims to optimize a compromise between feasibility, accuracy and computation time. It uses a block-based adaptive mesh refinement (AMR). Blocks can be viewed as macro-elements surrounding large regions of the domain. Provided that physical phenomena like wetting fronts stay within the same block, the block refinement can be kept, i.e. avoiding frequent remeshing and saving computational time. Local mesh modification can serve as a capturing technique to resolve physical phenomena. To do this, an error criterion has to be decided. In this study, the one chosen in [START_REF] Clément | Numerical simulation of flows in unsaturated porous media by an adaptive discontinuous Galerkin method: application to sandy beaches[END_REF] is used. Roughly speaking, it measures the numerical smoothness of the solution and it demonstrated its abilities to drive mesh adaptation to capture moving wetting fronts.

The solution may hold spurious oscillations at wetting fronts for vanishing or varying diffusion. Even though this phenomenon may be resolved by mesh refinement, the corresponding computational cost is substantial to maintain a physically acceptable solution. The weighted discontinuous Galerkin (WDG) framework [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF][START_REF] Proft | Discontinuous Galerkin methods for convection-diffusion equations with varying and vanishing diffusivity[END_REF] can help in this regard thanks to two key ingredients. The first one is the use of weighted averages instead of the standard arithmetic average inside the DG formulation. The second ingredient is to relax the face penalization used for continuity constraint by a coefficient. Such a penalty strategy turns out to tune automatically the amount of local penalty to regulate the degree of smoothness of the approximate solution. One should note that the WDG framekork was initially developed for convection-diffusion equations where the diffusivity variations fit the mesh. This is the case for heterogeneous media but it is not true for nonlinear diffusivity like in Richards' equation with moving wetting fronts. Thus, the WDG framework is not expected to perform just as well. However, a suitable adaptive mesh refinement like the one proposed here can be used as a capturing technique to work in synergy with the WDG framework. The details of this strategy are given in Clément et al. [START_REF] Clément | Numerical simulation of flows in unsaturated porous media by an adaptive discontinuous Galerkin method: application to sandy beaches[END_REF].

Model validation

The Vauclin et al.'s test-case [START_REF] Vauclin | Experimental and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem[END_REF] is used here to validate the model, as it is a classical but well-documented benchmark for Richards' equation-based models and because it can be viewed as a simplified version of the groundwater dynamics of a beach face subject to constant infiltration. The complete simulation of water table recharge by Vauclin et al. [START_REF] Vauclin | Experimental and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem[END_REF] was used by many other studies to assess their method, see [START_REF] Clement | A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow[END_REF][START_REF] Dogan | Saturated-Unsaturated 3D Groundwater Model. II: Verification and Application[END_REF][START_REF] Shen | A process-based, distributed hydrologic model based on a large-scale method for surface-subsurface coupling[END_REF][START_REF] Twarakavi | Evaluating Interactions between Groundwater and Vadose Zone Using the HYDRUS-Based Flow Package for MODFLOW[END_REF][START_REF] Xu | Integration of SWAP and MODFLOW-2000 for modeling groundwater dynamics in shallow water table areas[END_REF][START_REF] Zha | A generalized Ross method for two-and three-dimensional variably saturated flow[END_REF]. The validation of the code MODFLOW is partly done with this experiment dataset [START_REF] Thoms | User's guide to the Variably Saturated Flow (VSF) process to MODFLOW[END_REF]. One can notice that the majority of these simulations does not take into account the seepage boundary condition and prescribe an impervious wall instead which is in contradiction with the experimental set-up and numerical studies of Vauclin et al. In [START_REF] Clement | A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow[END_REF], the authors argue that a no-flow boundary can be used because the problem shows experimentally and numerically a steady-state solution without the development of a significant seepage face.

The simulation is carried for Richards' equation during T = 8 h. The mesh is a regular grid ∆x = ∆z = 5 cm. Monomial basis is used with p = 2. The penalty parameters are set to 100. The time step starts at 0.003 h and does not exceed 0.1 h.

Water table recharge

The transient position of the water table is reported every 2 hours, as shown in Fig. 2, and compared with experimental results taken from [START_REF] Vauclin | Experimental and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem[END_REF]. The results obtained from the numerical model closely agree with the experimental profile. The growth of water table mound is well simulated. Small discrepancies are observed in the middle of water table which is underestimated by numerical results. This observation is noticed by other study using different models [START_REF] Twarakavi | Evaluating Interactions between Groundwater and Vadose Zone Using the HYDRUS-Based Flow Package for MODFLOW[END_REF][START_REF] Vauclin | Experimental and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem[END_REF] and may be explained by the sandy soil which is non-perfectly isotropic and homogeneous as stated by hydraulic properties. Indeed, the recharge flux is redistributed horizontally by diffusion in the unsaturated zone before reaching the water table. Then, the unsaturated zone stores a portion of water inflow which moistens dry regions and does not contribute to the elevation of water table in this region. Similar observations are made in [START_REF] Xu | Integration of SWAP and MODFLOW-2000 for modeling groundwater dynamics in shallow water table areas[END_REF][START_REF] Shen | A process-based, distributed hydrologic model based on a large-scale method for surface-subsurface coupling[END_REF]. In addition, it may be difficult to identify experimentally the exact position of the water table because of the capillary fringe. One can see that seepage face length is very small and so, it may be disregarded for this experiment as stated in [START_REF] Clement | A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow[END_REF].

Results for wave-driven beach groundwater flow

The numerical model is applied to a series of large-scale experiments called BARDEX II [START_REF] Masselink | Large-scale Barrier Dynamics Experiment II (BARDEX II): Experimental design, instrumentation, test program, and data set[END_REF][START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF]. In these experiments, the beach groundwater dynamics induced by waves and cross-barrier gradients is investigated, with a focus on the swash zone.

Description of BARDEX II experiments

The BARDEX II experiment was performed in July 2012 at the Delta Flume in the Netherlands [START_REF] Masselink | Large-scale Barrier Dynamics Experiment II (BARDEX II): Experimental design, instrumentation, test program, and data set[END_REF][START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF].

One of the main goals of this large scale experimental project was to characterize, under controlled laboratory conditions, the effect of groundwater flow on the hydro-and morphodynamics in a sand barrier. In particular, we focus here on a series of experiments dedicated to the effect of waves and cross-barrier mean water level gradients on the groundwater dynamics. d 90 = 0.90 mm [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF]. Hydraulic conductivity at saturation is estimated 5

× 10 -4 m • s -1 ≤ K s ≤ 1 × 10 -3 m • s -1
[73] with an average value of 8 × 10 -4 m • s -1 [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF]. The porosity Φ [-] is evaluated 0.37 ≤ Φ ≤ 0.42 [START_REF] Masselink | Large-scale Barrier Dynamics Experiment II (BARDEX II): Experimental design, instrumentation, test program, and data set[END_REF] with Φ = 0.41 retained in [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF]. However, the relations and parameters governing the hydraulic conductivity and water content are not known in BARDEX II. For the present study, the van Genuchten-Mualem relations are chosen, see Tab. In the present confrontation tests with the model, we focus on six selected experimental cases allowing to differentiate or combine the effect of waves and cross-barrier MWL gradients on the groundwater field.

A6 and A6 w cases correspond to a zero cross-barrier MWL gradient (same elevation on both sides of the barrier), without and with wave forcing, respectively. A2 and A2 w cases correspond to a negative crossbarrier MWL gradient (lagoon MWL < sea MWL), without and with wave forcing, respectively. A4 and A4 w cases correspond to a positive cross-barrier MWL gradient (lagoon MWL > sea MWL), without and with wave forcing, respectively. For each wave case, the wave forcing is based on irregular random waves (JONSWAP spectrum) produced at the wave maker located at X = 0. The free surface forcing is provided by the SWASH model [START_REF] Zijlema | SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters[END_REF]. The BARDEX II simulation involves the main tools presented in this article: the WDG method in combination with the AMR technique. Fixed-point iteration is used together with adaptive time stepping.

Times steps are kept under 0.1 s. Mesh adaptation is done every five time iterations. Order of approximation is quadratic for space (p = 2). Penalty parameters are set to 100.

The actual initial groundwater state at the time periods selected for the experimental analysis of the BARDEX II cases is unknown. Turner et al. [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF] presented the experimental procedure which consisted in series of successive periods of wave action and calm conditions. In such circumstances, there is no guarantee for the attainment of a steady state for the barrier groundwater flow. In order to analyse conditions as close as possible to the steady state, they selected two characteristic time periods: one is taken before the first period of wave generation for the "no wave" conditions and the second one is taken at the end of the longest and last period of wave generation for the "wave" conditions. For the initialization of the numerical model: in the case of "no wave" conditions, the water table is assumed to be straight line which links the mean sea level and the lagoon level. The numerical initialization does not have much influence given that the steady-state equation is solved. Even though this is maybe not the case experimentally, it is important to look for a steady-state because the "no wave" conditions act like a controlled experiment;

in the case of "wave" conditions, the same initial condition can be used [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF] but it has been found that the Richards-based model does not provide good results after one run of 300 s because it is very sensitive to the initial unsaturated state. Thereby, the initial state of the BARDEX II time periods was reproduced numerically by launching successively 11 runs of 300 s. The final numerical solution of one run was kept to initialize the next one. The outcome of these 11 runs (roughly 1h) serves to initialize the analysed run of 300 s.

Comparison and assessment

The key numerical results for BARDEX II are presented in the next sections. They are gathered according to global and local points of view: time-averaged/resolved results and overall/swash dynamics. is made with the experimental BARDEX II results in Fig. 6a, 7a and 8a. Overall, very good agreement is found for the different cases. A discrepancy is observed on the lagoon side where the numerical water table is higher than the experimental one. But it seems that the experimental results in [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF] underestimate the lagoon water level because it is lower than the required experimental set-up (4.3 m). Small differences are observed for the case A4 but they stay within the range of the capillary fringe which can be misread by instruments. For the case A6, the experimental water table shows a bump near the swash zone whereas the numerical water table is an horizontal straight line. As no clear reason is seen for the presence of a bump in a fully balanced situation, one may assume that the experiment did not reach a true steady state. Contour plots of hydraulic head and flux vectors computed from numerical results are displayed in Fig. 6a, 7a and 8a.

Barrier groundwater dynamics

They show that the flow is completely driven by the lagoon water level in the absence of wave, i.e. by the cross-barrier difference of potentials induced by the sea and lagoon water levels. In the case A2, it leads to a small seepage face above the sea level. Time-averaged hydraulic head contours and flux vectors are also displayed. The model predicts very well the flow divide of BARDEX II experiments in Fig. 7b and 8b for the cases A4 and A6. The flow divergence takes place below the bump of the water table. This structure is built over time due to infiltration by the successive swash events. One should note that occasional strong swash event can have a lasting impact for the persistence of the bump-like structure because groundwater relaxation to equilibrium is slow. Groundwater pressure stays nearly "hydrostatic" in the dry area, in the upper beach.

Swash groundwater circulation cell

A focus is made on the area beneath the swash zone where the pore velocities give rise to a circulation cell.

Comparison is made in Fig. 9 to 11 between the numerical data (red vectors and contours) and experimental data (black vectors and grey contours) from Turner et al. [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF].

320

The cases with no wave conditions are presented in Fig. 9a, 10a and 11a. For the case A6, there should be no flow since the sea and lagoon levels are the same. In Fig. 11a, the numerical model shows no flow but a seaward flux is observed for the BARDEX II experiment. This confirms that the barrier did not attain a steady state despite the intention of experimenters. For the cases A2 and A4 both experimentally and numerically, the flow is mainly directed by the difference of potentials induced by the sea and lagoon levels which are different. The flux is seaward in Fig. 9a (raised lagoon level -A2) and landward in Fig. 10a (lowered lagoon level -A4). In the barrier zone between the swash tip and the lagoon, the flow is predominantly horizontal. However, as it comes closer to the swash tip, the velocity is decreasing and a vertical component is observed. The latter gets stronger below the surface water where it is upward in Fig. 9a and downward in Fig. 10a. This is observed both in numerical and experimental results but the numerical models tends to overestimate the vertical component compared to the experiments. The Darcy-based numerical model from Perera et al. [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF] led to the same discrepancy. They suggest that the experimental velocity component is not as strong because of localized unsaturated conditions which decrease the permeability near the beach surface. Even though this could be possible, we think it can be explained by the same reason for the case A6: the steady state was not reached yet by the experiments. Apart from this discrepancy, a good agreement is found in terms magnitude of flux between the present numerical model and the experiments.

The cases with waves action are presented in Fig. 9b, 10b and 11b for time-averaged results. Contrary to the cases without wave, the flow is now directed seaward, no matter the lagoon conditions are. For the cases A4 and A6, the hump-like structure of the water table built by the successive infiltrations forms a water potential higher than the lagoon and sea levels, see Fig. 7b and8b. Thereby, a circulation is established with infiltration in the upper part of the swash zone and exfiltration in the lower part. This confirms the observations of Turner et al. [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF], supported by Perera et al. [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF], which show that wave action is powerful enough to rule the groundwater dynamics in this area (see also Sous et al. [START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF] for similar observations on the field). As noted experimentally by Turner et al. [START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF] and simulated by the present model, the cases A4 and A6 show greater seaward hydraulic gradients than the case A2, see the tightened equipotential lines in Fig. 10b and11b. It means that local gradients driven by waves dominate the overall gradients induced by the difference of sea and lagoon levels. Concerning all the cases, while experimental flow is mainly horizontal, the numerical flow holds an upward vertical component near the beach surface in the left part of the domain.

This discrepancy is also stated by the Darcy-model from Perera et al. [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF]. The fact this is observed with waves action for every case may lead to think that the modelling in this zone needs further improvements.

It can be because of either the surface water modelling which overestimates the arch of the swash tongue, or the interface condition between the surface or the groundwater models which misestimates vertical flows.

Apart from this discrepancy, a good agreement is found in terms magnitude and direction of flux between the present numerical model and the experiments. In particular, no major deviance has been found for the case A4 with waves contrary to the numerical simulation of Perera et al. [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF].

Spectral analysis

The analysis of spectral properties of groundwater dynamics is performed for the numerical results of BARDEX II cases in terms of hydraulic head at different locations under the swash zone. The locations are depicted in Fig. 5. Results are displayed in Fig. 12 to 14. Firstly, one can note the transformation from the wave-maker location (blue line) to the lower part of the swash zone (orange line) in Fig. 14a. The energy 360 of generated waves is divided into a lower gravity band (0.1 ≤ f ≤ 0.2 Hz for A2 and A4, 0.06 ≤ f ≤ 0.18 Hz for A6) and a smaller component in the infragravity band (f ≤ 0.06 Hz for A2 and A4, non-existent for A6). At the swash zone, the energy of the gravity band has been reduced and transferred partly to the infragravity band. Besides, a small energy peak is appearing around f = 0.21 Hz for the cases A2 and A4 and f = 0.16 Hz for the case A6 which is probably related to swash-swash interactions.

365

In Fig. 12b, 13b and 14b, one can observe that the groundwater hydraulic head field fluctuates in the same frequency range as the free surface in the overlying swash zone. This is expected because the porous medium is almost always saturated underneath the swash zone. Then, Richards' equation is purely diffusive which means that information propagates at infinite speed. Yet, it does not mean that flow velocity is the same since pore velocity is driven by gradient of hydraulic head throughout the porous medium. Note also that spectral analysis of groundwater head fluctuations helps to understand how free surface wave energy propagates into the soil but does not give a direct insight on the groundwater flows since these latter result also from phase shift of pressure waves.

The spectra of hydraulic head highlights the filtering effect of the porous medium in the depth and cross-shore direction: (i) the deeper and (ii) the more landward the point in the swash zone, the smaller the remaining energy is. This damping is frequency-dependent, with stronger for the frequencies above the infragravity band, revealing the low-pass filtering role played by the sand soil. Energy attenuation is much more effective horizontally than vertically. Indeed, if the beach is saturated, the pressure of the water column is transferred directly to the pore fluid vertically with little attenuation. Horizontally, different mechanisms are involved but they are mainly controlled by diffusion which depends on the groundwater conditions. Thus, in this direction, energy is more attenuated. In addition, it should be noted that the fluctuations of heads stop quickly beyond the swash zone since the frequencies drop below f < 10 -2 at X = 100 m (data are not showed) for all cases. As a consequence, the barrier groundwater flow can be divided into two zones with different dynamics: the first one corresponding to the surf and swash zones is very active while the second in relation with the back-barrier remains relatively calm. For the cases A4 and A6, the delimitation of these two dynamics zones nearly matches the flow divide line.

Time-resolved swash event

Instrumental observations generally lack resolutions to measure correctly the groundwater response to a single swash event, especially the local patterns or the rapid water table fluctuations. Few experimental descriptions have been carried out in [START_REF] Sous | Swash-groundwater dynamics in a sandy beach laboratory experiment[END_REF][START_REF] Sous | Field evidence of swash groundwater circulation in the microtidal rousty beach, France[END_REF] but it stands out that numerical models are needed to bring A swash event is selected for the cases A2, A4 and A6 whose snapshots are respectively gathered in Fig. 15 395 to 17.

The overall trend for the different cases of BARDEX II is as follows:

A first observation concerns wave propagation. When a wave is moving, the groundwater flow un-derneath in the saturated zone reacts immediately. Exfiltration occurs in front of the bore, due to a low-pressure zone induced by the trough, while infiltration occurs under the bore due to a high-pressure zone induced by the crest. It is particularly visible in Fig. 15b and16b;

Velocities magnitude is rather small everywhere, below 5 × 10 -5 m/s, except below the propagating wave and swash tip at the end of uprush, where there is strong infiltration (between 5 × 10 -5 m/s and 8 × 10 -4 m/s). Moreover, flow divides into seaward and landward components near the uprush limit as seen in Fig. 15d, 16d and17d;

A global pattern can be depicted for swash cycles. Two recirculation zones are attached respectively to the back (low-pressure) and the front (high-pressure) of the swash tongue. These two groundwater cells are interacting together and are evolving during uprush and backwash, growing or diminishing.

Sometimes the low-pressure cell vanishes completely (Fig. 15d) and the high-pressure cell may be fixed by the recharge induced by the water table elevation (Fig. 16c);

During uprush, while the swash tongue is inundating the beach, downward infiltration occurs. It is mainly concentrated in the zone over the fully saturated part of the beach (Fig. 16b). The uprush makes the water table rise by (i) raising the exit point of the water table and (ii) the induced-infiltrations percolating in the unsaturated zone. Rapid uprush can create trapped unsaturated pockets as in Fig. 17b which obstruct the downward infiltrations;

During backwash, a low hydraulic head zone emerges at the base of the swash zone where exfiltration occurs. Velocities are mainly seaward under the swash tongue moving back. Water brought by uprush in the unsaturated zone is infiltrating, until reaching the hump-like structure of the water table. While the swash toe is retreating to the sea, the water table exit point is following with a small delay (Fig. 15d). For strong and fast backwash, a wide seepage face may form (Fig. 17d).

The water table's back and forth movement car create locally trailing fluctuations in its topology (two hump-like structures), see Fig. 16a and17a.

It should be noted that this swash event has been selected to be the clearest possible but the swash-swash interactions create more complex patterns which interfere with or delay these general observations. This is why for example the swash event of the case A6 in Fig. 17 is longer than the two others.

Infiltration and exfiltration involve greater pore velocities than in other parts of the beach. Like it is depicted by Perera et al. [START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF], gradients brought by the wave are much larger than the average values in the beach due to larger differences in water depth across the wave. The generalised observation that net infiltration predominates in the upper swash zone, in contrast to net exfiltration across the lower beach face, confirms the results of numerical simulation presented in Bakhtyar et al. [START_REF] Bakhtyar | Wave-induced water table fluctuations, sediment transport and beach profile change: Modeling and comparison with large-scale laboratory experiments[END_REF] and match the experimental results of Sous et al. [START_REF] Sous | Swash-groundwater dynamics in a sandy beach laboratory experiment[END_REF]. Such distinction must be considered in the conceptual framework summarized by Turner and Nielsen [START_REF] Turner | Beach Dewatering as a 'Soft' Engineering Solution to Coastal Erosion-A History and Critical Review[END_REF]. They showed that pressure fluctuations under the swash zone are related to oscillations of the water table. At the swash motion, the groundwater dynamics is mainly attributed to the displacements of the water table inside the capillary fringe thickness.

Discussion

435

This section is organised in three main topics: the potential groundwater effect on sediment dynamics, the overall drivers of beach groundwater dynamics and the numerical strategy.

Sediment mobility by swash groundwater flows

The role played by swash-induced groundwater flows on sediment transport is a debated question [START_REF] Elfrink | Hydrodynamics and sediment transport in the swash zone: a review and perspectives[END_REF][START_REF] Bakhtyar | Modeling sediment transport in the swash zone: A review[END_REF][START_REF] Chardón-Maldonado | Advances in swash-zone research: Smallscale hydrodynamic and sediment transport processes[END_REF] to assess the beach erosion/accretion. The coupled sediment-groundwater-surface flow dynamics is 440 a very complex issue, which involves a wide range of processes outside the scope of the present research, such as granular mechanics and groundwater-surface flows interactions at very fine scale. The development of a comprehensive modeling tools will remain a challenge for the incoming years. However, considering their importance in terms of beach morphodynamics, we place our discussion in the simplified framework proposed these last decades to estimate the net effect of swash groundwater dynamics on the beach face sediment transport. For sandy beaches, two mechanisms are supposed to be in competition to explain why infiltration/exfiltration affect sediment stability: the relative weight modification (RWM) of sediment particles by seepage forces and the alteration of bed shear stress because of the boundary layer. The Shields parameter can account for the balance of these two processes [START_REF] Turner | Swash infiltration-exfiltration and sediment transport[END_REF]. Butt et al. [START_REF] Butt | The influence of swash infiltration-exfiltration on beach face sediment transport: onshore or offshore?[END_REF] show that the effect of through-bed flows on the swash boundary layer can be neglected compared to the RWM of sediments. Thus, for this study, the aim is to give a first insight about the question by investigating the effective weight of the sediment which leads to bed destabilization/stabilization. Following Butt et al. [START_REF] Butt | The influence of swash infiltration-exfiltration on beach face sediment transport: onshore or offshore?[END_REF], the RWM, denoted here by µ RWM , is computed according to:

µ RWM = - β s -1 w K , (5)
where

K [L•T -1
] is the (variably-saturated) hydraulic conductivity in the direction of the beach bed, w

[L•T -1
] is the flux normal to the beach bed, s [-] is the specific gravity of the sediment (relative density) and β is a dimensionless coefficient giving the strength of the downward drag. It is assumed that s = ρ sand ρ water ≈ 1.65 and β = 0.5 [START_REF] Turner | Swash infiltration-exfiltration and sediment transport[END_REF].

Computation of the RWM is performed on the numerical results of the different cases of BARDEX II for five points (S1 to S5) lying on the beach face and whose X-position is determined to be evenly spread in the swash zone. The latter is delimited by the lowest and the highest point reached by the swash tip over time. Then, it has been found that the swash zone stretches from X = 86 m to X = 98 m for the case A2, from X = 85 m to X = 99 m for the case A4 and from X = 84 m to X = 103 m for the case A6. As a consequence, the S1 to S5 points have different absolute position but the same relative position in the swash zone. The probability density of the RWM during the 300 s run is showed in Fig. 18 for the S1 to S5 points in each BARDEX II case. A first observation is that the RWM is small for the considered period (300 s), mostly under 0.1%. Besides, the distribution of RWM appears to be quite symmetric and only slightly deviated from zero for each point. This means lots of swash cycles are needed to potentially affect the sediment transport by groundwater effects. The overall tendency according to the location in the swash zone follows the observations previously performed on groundwater flows. On the lower part of the swash zone, distribution of RWM is deviated towards negative values, meaning that exfiltration lightens the effective weight of sediments because of wave-induced pressure gradients and groundwater discharge/seepage during backwash. In the upper swash zone, the dynamics is different. The RWM is now increased in relation with infiltration which occurs during uprush events.

The influence of the elevation of the water table due to the lagoon level is evaluated between the different BARDEX II cases. A lower water table (cases A4 and A6) tends to widen and to shift towards greater values the distribution of the RWM, as illustrated in Fig. 18b and18c. This means that (i) infiltration/exfiltration events in the case of lower water table are more spread for a given location in the swash zone, i.e. pore velocities and saturation show larger range of values, and (ii) the effective weight of sediments is increased.

These observations may be explained by a medium more often unsaturated where downward infiltrations take place during uprush but seepage exfiltration is prevented during backwash. Figure 18a shows narrow distributions of RWM due to dominant influence of the upper water table to drive the groundwater flows compared to swash. This statement is supported by the absence of flow divide in this case. Finally, one can also remark there is no noticeable difference between the lowered (case A4) and medium (case A6) elevation of the water table. The computation of RWM was also carried out for the same abscissas of S1 to S5 points but at a depth of 5 cm under the beach face. Results were very similar to the ones presented in Fig. 18 for the different cases, with slightly more centered and narrowed RWM distributions. We did not observe any particularly intense hydraulic gradient near the beach surface as experimentally measured by Baldock et al.

[80] and suspected to lead to bed fluidization.

For operational purposes, the present observations may suggest that lowering the beach water table should help to prevent erosion in the case of beaches for which the water table is known to be high most of the time. However, it should be kept in mind that the RWM values observed here are small, and therefore that the groundwater control on sediment dynamics is expected to be overpowered by many physical processes acting in the swash zone. BARDEX II. S1 to S5 are five points lying on the beach face and whose X-position is selected to be evenly spread within the swash zone defined for each BARDEX II case.

Drivers of beach groundwater dynamics

The numerical data of BARDEX II point out some of the drivers of beach groundwater flow identified by Santos et al. [START_REF] Santos | The driving forces of porewater and groundwater flow in permeable coastal sediments: A review[END_REF]: the wave setup, the water level difference across the barrier, the flow/topography-induced pressure gradients, the wave pumping and the shear. The following observations can be summarized from the present modeling study.

The wave setup results in a mean recirculation cell with infiltration and exfiltration respectively in the upper and lower parts of the swash zone, as showed in Section 4.2.2;

The water level difference across the barrier plays a significant role in the no-wave cases (see Fig. 6a, 7a and 8a);

The bottom topography interacts with flows causing local differences of pressure on the bed. When the flow meets an obstacle (like a bump), the pressure is increased, causing infiltration. Conversely, when the flow accelerates passing over the obstacle, the pressure is decreased, causing exfiltration. In the case of wave-induced oscillatory flow, a symmetric groundwater pattern takes place at the ripples.

The numerical simulations of BARDEX II show this phenomenon for the bars developed on the beach face. The groundwater flow pattern is particularly relevant in Fig. 8b for the bar of case A6. However, it should be noted that the involved time-averaged velocity is low and so it concerns mainly pressure gradients;

The difference of hydrostatic pressure underneath wave crests and troughs produces a localized groundwater pumping. Along its propagation, the wave carries a low-pressure zone at the front and a highpressure zone underneath, resulting respectively in exfiltration and infiltration whose magnitude decreases with increasing depth. This process is presented in Section 4.2.4 for time-resolved events with It is interesting that BARDEX II barrier beach can be divided into three parts where different drivers seem to be dominant for the groundwater flow: (i) the toe and surf zones with topography-induced flow and wave pumping, (ii) the swash zone for much complex processes with wave setup, shear, seepage and water table fluctuations, and (iii) the back-barrier with lagoon-induced potential and capillary forces. In addition, to the groundwater dynamics drivers identified by Santos et al. [START_REF] Santos | The driving forces of porewater and groundwater flow in permeable coastal sediments: A review[END_REF], the unsaturated mechanisms addressed by the present study should be considered, including the role played by capillary-raised water and the swash-induced mound in the water table built by infiltration of successive swash cycles.

Wave-resolved and variably-saturated numerical modelling

The present approach is innovative as it is the first wave-resolved and variably-saturated numerical model used to simulate beach groundwater dynamics induced by swash. Compared to previous studies mostly based on Darcy's law, the present model is based on the solving of Richards' equation which extends the study of beach groundwater to the unsaturated part. In particular, the capillary fringe -the porous zone almost entirely filled by water but with negative pressure head -can be predicted. This is important because many instrumental setups are likely to misread the capillary fringe and its effect on seepage, beach water storage and water table position. In addition, Richards' equation is a major asset compared to Darcy's law to get access to the flux in the unsaturated zone, which is of primary importance for swash zone dynamics exposed to cyclic water level fluctuations. Richards' equation is solved at the wave event scale, allowing a variety of time/space/spectral/distribution analysis.

The Richards' equation has some limitations. In terms of modelling, it cannot be used for rapid groundwater flows and when air phase phenomena with dynamic pressure are involved. This framework remains reasonable for fine sand beaches, but would certainly become irrelevant for larger sediment. Besides, hydraulic properties (hydraulic conductivity and water content) need to be defined for the unsaturated part.

Relations exist in the literature but outside the beach groundwater context, calling for dedicated experimental field work to reduce the uncertainty associated to the lack of knowledge of hydraulic parameters. In addition, dynamic effects and hysteresis of hydraulic properties are not taken into account into the present model due to the lack of experimental knowledge, even though it would be possible from theoretical and numerical aspects. The expected gain in this direction seems however to be limited compared to the complexity it would bring. Finally, the one-way forcing implies that no groundwater-to-surface flow feedback is accounted for. The successful comparison with the BARDEX II experiments indicates that this framework is valid for the fine sand situation, at least for the parameters monitored here. Some discrepancies are observed on the near-surface groundwater flow in the swash zone. At such fine scale, corrections of Richards' equation can be envisioned to take into account fast flow dynamics for example. This would also enlarge the scope of the targeted types of beaches (larger sediment) and waves which can be considered. It is also expected that a full coupling will provide increased performance, together with an implementation of inter-granular stresses model to reach a comprehensive view of the sediment matrix interaction with groundwater and surface flow.

In a longer term prospect, a natural extension would be to integrate a transport model to study the transport of sediments (beach erosion/accretion), nutrients or pollutants within the swash-forced beach to assess the biogeochemical cycles in this important exchange area between land and ocean.

The numerical model gives promising results and offers lots of opportunities to be further improved to make accurate, efficient and robust wave-resolving simulations of Richards' equation. In particular, dedicated linearization schemes would help the nonlinear convergence and specific stabilization methods would prevent systematically oscillations in the solution. From a computational view, the adaptive meshing strategy can be automated to help the user with the parameters and the DG methods can be quite easily parallelized to save time. For long term prospects, the high-order and local abilities of DG methods allow other developments like order adaptation and local solvers to reach high accuracy and efficiency. Finally, DG methods are also suitable to solve surface wave equation. Then, the latter can be done along with Richards' equation computation to make a full coupling.

Conclusion

This paper addresses the challenge of the numerical modelling of wave-driven beach groundwater dynamics accounting for partial saturation conditions. The resulting code, called Rivage, is based on Richards' equation. The model is forced at the seabed by the pressure field computed by the non-linear and nonhydrostatic shallow water SWASH model. This one-way forcing relies on the assumption that groundwater feedback on surface water dynamics can be neglected because pore velocity is much slower than wave velocities, restricting the application domain to fine-grained beach.

The Rivage code has been first validated on the classical Vauclin's experiment and then used to simulate BARDEX II large-scale experiments of a sand barrier beach exposed to irregular waves and various mean water elevation gradients. Based on the numerical results, the following hints and remarks can be outlined.

While a satisfactory overall agreement is found between the BARDEX II experiments and the Rivage model predictions, the experimental data show some limitations which make sometimes difficult the direct comparison with numerical data. Steady state seems to be hardly achieved experimentally. The lack of precision/resolution for the measurements is problematic to reconstruct small velocities from pressure head field (just beneath the beach surface) and to locate the water table in relation with the capillary fringe. In addition, Richards' equation needs data which are often omitted by the experiments and are important to initialize the numerical model and choose hydraulic properties for the unsaturated part. Sediments heterogeneity/anisotropy may also be important at small scales.

Confirming and refining the experimental observations, the numerical predictions highlight the role played by waves in the beach groundwater dynamics. In the absence of waves, the groundwater head field is governed by the difference of potentials induced by sea and lagoon water levels. By contrast, the BARDEX II cases with waves show the creation of a bump-like structure in the water table where infiltration occurs during uprush. This leads to a flow divide of the groundwater with a very active part seaward (driven by waves and swash) and a relatively calm part landward (controlled by the lagoon elevation).

Swash motion induces a mean groundwater circulation cell with infiltration in the upper part and exfiltration in the lower part, no matter the lagoon elevation is.

At the swash event scale, two localized groundwater circulation zones evolve and interact. One, associated with low-pressure (exfiltration), is attached to the swash base during backwash and to the forefront of the next wave during uprush. The other one, associated with high-pressure (infiltration), is attached to water table recharge during backwash and to the wave crest during uprush. This latter elevates the water table next to beach face. Strong occasional events of uprush can have a (i) delaying effect due to water partially stored in the unsaturated zone and a capillary barrier effect but (ii) a lasting impact because the groundwater relaxation time to equilibrium is slow. If the uprush is fast enough, it creates local unsaturated zones. In the upper unsaturated part of the beach, the uprush results in sharp wetting fronts. They propagate downward until reaching the water table. Thus, water table elevates locally to shape a bump-structure which induces a high-head zone from which the groundwater flow divides into landward and seaward components. During backwash, the overall groundwater trend is seaward under the swash tongue moving back to the sea. Water table exit point follows the swash tip with delay. Seepage face can develop only when strong and fast backwash happens.

The beach sediments act as a low-pass filtering, with only the low band of infragravity waves being able to deeply penetrate the barrier. Close to the swash zone sand surface, the groundwater fluctuations are strongly connected to the swash dynamics.

The relative weight modification (RWM) of sediments remains small in the swash zone which means lots of swash cycles are needed to show a potential cumulative effect on beach morphology, if there is any. The overall trend is that the effective weight of sediments decreases in the lower swash zone whereas it increases in the upper swash zone, which is consistent with the groundwater circulation cell.

A lower water table seems to shift slightly the RWM towards positive values by providing unsaturated conditions more often. This might help to prevent erosion, although the groundwater effect should likely be overpowered by many morphodynamical processes acting in the swash zone.

The BARDEX II simulations highlight the spatialization of the groundwater flow drivers. In the toe and surf zones, topography-induced flow and wave pumping are dominant, in the swash zone, wave setup/shear/seepage/water table fluctuations show a complex groundwater pattern shaping globally a circulation cell, and in the back-barrier, lagoon-induced potential and capillary forces drive the flow. This general beach groundwater scheme is consistent with observations made by the experimental and numerical studies in the literature. The numerical model allows a comprehensive and refined analysis of the beach groundwater, out of the reach for usual experimental approaches. Additional developments and investigations are needed to further validate and to improve the present model. On the one hand, the selection of the unsaturated hydraulic properties needs to be improved as well as the initialization of the saturation state of the beach, both issues being of primary importance for the model accuracy and difficult to perfectly control experimentally. On the other hand, if one wants to catch more complex processes and study more types of waves/sediments, the model should be completed in priority with a fully coupled boundary layer model for the interface between surface and subsurface models where high pore velocities can happen, an evolving beach topography describing erosion/accretion including intra-granular stresses and a heterogeneous density for water (fresh/salt). Such features would provide a complete model to assess the transport of sediments, contaminants or nutrients within the beach.

  : empirical shape parameters [L B;D ] B, D: empirical shape parameters [-]Van Genuchten-Mualem relations (1980)[START_REF] Mualem | A new model for predicting the hydraulic conductivity of unsaturated porous media[END_REF][START_REF] Van Genuchten | A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils[END_REF]

3. 1 .

 1 Description of Vauclin's benchmark Vauclin et al. [66] experimented and simulated water flow recharge through a slab of soil. It consists in a box of 6 m by 2 m whose only one half is simulated because of symmetry. Sides are impervious on the left, at the top for x > 50 cm and at the bottom: ∇h • n = 0. A constant flux is prescribed on the top for x ≤ 50 cm: -K(ψ)∇h • n = 14.8 cm/h. The water level is kept constant in the ditch on the right for z ≤ 65 cm: h = 65 cm. The remaining boundary on the right for z > 65 cm accounts for a seepage boundary condition. Initial state is at hydrostatic equilibrium with water table at z = 65 m. Test-case configuration is depicted in Fig. 1. Hydraulic properties are taken from Vachaud's relations (see Tab. 1) with A = 2.99 × 10 6 , B = 5.0, C = 40000, D = 2.90, K s = 35 cm/h, θ s = 0.30 and θ r = 0.

Figure 1 :

 1 Figure 1: Watertable recharge test-case configuration.

Figure 2 :

 2 Figure 2: Numerical simulation of two-dimensional recharge with water table position compared to experimental data from Vauclin et al. [66].

Figure 3 :

 3 Figure 3: Distribution of hydraulic head h (turned upside down) and flux, location of water table (white line) and capillary fringe (magenta line) at t = 3 h (top) and t = 8 h (bottom) for water table recharge. Hydraulic head contours are in red.

Figure 3

 3 Figure 3 illustrates a typical field distribution of fluxes and hydraulic head as well as the location of water table and capillary fringe (set at θ = 0.29), at time t = 3 h and t = 8 h. Hydraulic head values have been turned upside down to correspond with a downward vertical axis and to facilitate comparison with the numerical and experiment results from Vauclin et al. [66]. At t = 8 h, Vauclin et al. [66] stated that a 195

Figure 4 :

 4 Figure 4: Steady-state simulation of water table recharge with distributions of hydraulic head h and flux. Hydraulic head contours are in red. Calculated water table and capillary fringe are respectively represented by a blue solid line and a blue dotted line. Those of Vauclin et al. [66] are represented by an orange solid line and an orange dotted line.

Figure 5

 5 Figure 5 depicts the experimental set-up. The sand barrier was 4.5 m high and 5 m wide. The barrier profile comprised five successive sections: (i) a concrete toe of slope 1:10 (24 < X < 29 m), (ii) a horizontal section (29 < X < 49 m), (iii) a beach face of seaward slope of 1:15 (49 < X < 109 m), (iv) a horizontal crest of height 4.5 m (109 < X < 114 m) and (v) a landward slope of 1:5 (114 < X < 124 m). The barrier is separated from a 10 m long lagoon by a permeable retaining wall. The water level in the lagoon is controlled by a pump system. A series of pressure transducers were buried in the barrier to monitor the groundwater pressure head. Beach sand of BARDEX II has a sorting of a medium-sized sand: d 10 = 0.26 mm, d 50 = 0.42 mm,

1 .

 1 The parameters are taken from the study of Benson et al. [74] in which the Vail sand (Winsconsin Geotechnics Laboratory measures) shows sorting and porosity properties similar to the sand used in BARDEX II: d 10 = 0.22 mm, d 60 = 0.55 mm, θ s = 0.412 and θ r = 0.030 (Φ = θ s -θ r = 0.382). Van Genuchten-Mualem parameters of Vail sand are α = 4.6 m -1 and n = 5.14. For the numerical model, θ s = 0.41, θ r = 0.03 and K s = 8 × 10 -4 m • s -1 .

Figure 5 :

 5 Figure 5: Description of experimental set-up for BARDEX II benchmark cases. Irregular waves are generated from a wavemaker at X = 0 m. The mean water level (MWL) at sea boundary is kept constant while the lagoon MWL can be maintained to a higher (A2), lower (A4) or same (A6) level. The diamonds indicate the points chosen for spectral analysis of Section 4.2.3. The boundary conditions prescribed for the simulations are shown with different colours.

  The model is implemented along the studied flume profile (2DV simulations) with two vertical layers and a horizontal resolution of 0.1 m. The model is forced by 5-min irregular waves synthetic time series built from a Jonswap spectrum based on the experimental wave parameter H s and T p . Second-order bound infragravity waves are added at the offshore boundary of the model based on weakly non-linear finite-depth wave theory [75, 76]. The resulting bottom pressure data are used to enforce a forcing boundary condition along the seaward beach face of the barrier (X = 29 -109 m). Boundary conditions for the simulations are depicted in Fig. 5. Since the concrete toe is impervious, the present groundwater numerical model does not take into account this structure where a zero-valued Neumann boundary condition is prescribed. The flume floor (X = 29 -124 m) is impervious and there is no flux at the barrier crest (X = 109 -114 m) which means a zero-valued Neumann boundary condition can be prescribed for both. To model the lagoon, a constant Dirichlet boundary condition (equal to lagoon water height) is prescribed at the lagoon side both along the vertical wall and the landward-sloping back barrier (X = 114 -124 m). A seepage boundary condition is prescribed for the remaining part above the lagoon water level because outflow can happen according to the groundwater conditions.

Figures 6 to 8

 8 Figures 6 to 8 represent the global results of groundwater dynamics obtained by the Richards-based numerical model for BARDEX II. Groundwater flow is simulated for the test series A2, A4 and A6 without wave by solving the steady-state Richards' equation, i.e. without considering the time derivative. Comparison of the water table position

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Comparison of the barrier groundwater dynamics between the numerical model and the run A2 of BARDEX II experiments. For the numerical model, hydraulic head distribution is showed together with contours in grey. Vectors show the flux from the numerical model. There is no flow divide.

( a )Figure 9 :Figure 10 :Figure 11 :

 a91011 Figure 9: Comparison of the circulation cell beneath the beach face between the numerical model and the case A2 of BARDEX II experiments. Red (numerical model) and black (experiments) arrows indicate the flux at their origin. They are scaled by their magnitude. Red (numerical model) and grey (experiments) equipotential lines show hydraulic head spaced by 1 cm. The blue line is the surface water computed by the numerical model and the yellow region represents the barrier.

( a )Figure 12 :Figure 13 :Figure 14 :

 a121314 Figure 12: Power spectral density (periodogram) for the simulation of the case A2 from BARDEX II.

390

  insights about such events. Perera et al.[START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF] used their Darcy-based model to look into groundwater flows for a swash event in BARDEX II. Clément et al.[START_REF] Clément | Wave-driven Groundwater Flows in Sandy Beaches: A Richards Equation-based Model[END_REF] used a Richards-based model to investigate the groundwater response to the action of swash in the case of an idealized beach. In the continuity of such studies, some assessments of time-resolved dynamics are here considered for the BARDEX II experiments.

Figure 15 :Figure 16 :Figure 17 :

 151617 Figure 15: Snapshots during a swash event for the simulation of the case A2 from BARDEX II. Blue and white solid line show respectively the surface water level and the water table. Magenta dotted line shows the capillary fringe. Grey contours show the equipotential lines of hydraulic head.

Figure 18 :

 18 Figure 18: Probability density function for the relative weight modification of sediments in the swash zone computed for

Fig. 15 to 17 .

 17 Fig. 15 to 17. Wave pumping should result in no net flux along the beach face which is supported by Fig. 6b, 7b and 8b where time-averaged groundwater shows in/outflow-symmetric cell-like patterns of very low velocity in the barrier toe and face (before the swash zone);

Table 2 :

 2 The presented results are obtained from 300 s model runs. Significant wave height H s and peak period T p are 0.8 m and 8 s, respectively. BARDEX II experimental cases.

	Experimental parameters are summarized in Tab. 2.

  The flow divide determined by Turner et al.[START_REF] Turner | Groundwater fluxes and flow paths within coastal barriers: Observations from a large-scale laboratory experiment (BARDEX II)[END_REF] from experiments (A4 and A6) is also plotted.For every case, waves raise the water table near the sea side, eventually shaping a bump-like structure in the cases A4 and A6. This is also observed by BARDEX experiments, confirming previous observations[START_REF] Sous | Swash-groundwater dynamics in a sandy beach laboratory experiment[END_REF]. However, some discrepancies are observed between numerical and experimental data about the exact location of the water table. One should note the Darcy-based model of Perera et al.[START_REF] Perera | Surface-groundwater flow numerical model for barrier beach with exfiltration incorporated bottom boundary layer model[END_REF] provided results closer to experimental data. For the present numerical model, one can notice that the experimental water table is almost always in the range of the numerical capillary fringe or in the unsaturated area under strong infiltrations. This may suggest that, either the parameters of unsaturated hydraulic properties are not good, and/or the measurement of BARDEX II partly misread the capillary fringe. In the model of Perera et al., the water table follows the swash tip and no flow occurs in the unsaturated zone. As a consequence, the position of water table should be overestimated because the unsaturated beach, particularly the capillary fringe, should store a substantial quantity of water, as observed in the present study. In support of this comment, the present numerical model shows a disconnection between the water table and surface water elevation. Besides, occasional event of strong uprush can have a delayed impact on the water table because infiltration in the unsaturated zone is slow (and even stopped by a capillary barrier effect if the medium is

The numerical model based on Richards' equation is computed for the test series A2, A4 and A6 with waves. Results are time-averaged over the 300 s run. Comparison of water table location is shown in Fig.

6b

, 7b and 8b. very dry) due to water being partially stored. A better understanding on these issues would likely require further measurements with fine control of local groundwater properties.

Acknowledgments

The authors are grateful to Deltares and the involved persons for the provision of the BARDEX II dataset.

Exchanges between the authors were facilitated through the GLADYS group (www.gladys-littoral.org).