Modelling storm hydrodynamics on gravel beaches with XBeach-G

In this paper we present a process-based numerical model for the prediction of storm hydrodynamics and hydrology on gravel beaches. The model comprises an extension of an existing open-source storm-impact model for sandy coasts (XBeach), through the application of (1) a non-hydrostatic pressure correction term that allows wave-by-wave modelling of the surface elevation and depthaveraged flow, and (2) a groundwater model that allows infiltration and exfiltration through the permeable gravel bed to be simulated, and is referred to as XBeach-G. Although the model contains validated sediment transport relations for sandy environments, transport relations for gravel in the model are currently under development and unvalidated. Consequently, all simulations in this paper are carried out without morphodynamic feedback. Modelled hydrodynamics are validated using data collected during a large-scale physical model experiment and detailed in-situ field data collected at Loe Bar, Cornwall, UK, as well as remote-sensed data collected at four gravel beach locations along the UK coast during the 2012-2013 storm season. Validation results show that the model has good skill in predicting wave transformation (overall SCI 0.14–0.21), run-up levels (SCI <0.12; median error <10%) and initial wave overtopping (85–90% prediction rate at barrier crest), indicating that the model can be applied to estimate potential storm impact on gravel beaches. The inclusion of the non-hydrostatic pressure correction term and groundwater model is shown to significantly improve the prediction and evolution of overtopping events. NOTICE: this is the author’s version of a work that was accepted for publication in Coastal Engineering. Changes resulting from the publishing process, such as peer review, editing, corrections, structural formatting, and other quality control mechanisms may not be reflected in this document. Changes may have been made to this work since it was submitted for publication. A definitive version was subsequently published in Coastal Engineering, Vol 91, September 2014 DOI: 10.1016/j.coastaleng.2014.06.007


Introduction
Gravel beaches and barriers occur on many high-latitude, wavedominated coasts across the world. Due to their natural ability to dissipate large amounts of wave energy, gravel coasts are widely regarded as an effective and sustainable form of coastal defence. However, during extreme events waves may overtop, overwash, and even lower, the crest of the gravel beach, flooding the hinterland. Although rare, such events can lead to loss of lives and significant damages to land and infrastructure.
In the UK, gravel is routinely used to nourish the coast (Moses and Williams, 2008). Despite this practice and previous research done to describe gravel barrier response to extreme storm events in a qualitative sense (e.g., Orford, 1977), coastal managers are currently forced to rely on simple empirical models to make quantitative predictions of gravel beach storm response and associated flooding risk (e.g., Bradbury, 2000;Powell, 1990). Although these empirical models have been applied with some success in the UK (e.g., Cope, 2005), they are inherently limited in their application by the range of conditions and data from which they are derived (cf. Bradbury et al., 2005;Obhrai et al., 2008). Since these models have been developed for relatively uncomplicated natural coastlines, managed coastlines (approximately 44% of the England and Wales coastline; DEFRA, 2010) containing man-made flood defence and beach regulation structures (seawalls, dikes, groynes) cannot be easily simulated using such models. More importantly, the application of these models outside their range of validity has been shown to underestimate the possibility of barrier overwash and breaching during storm conditions (cf. McCall et al., 2013;Van Rijn, 2010), leading to unsafe estimates of flooding. It is clear that these limitations will inhibit the use of such models to make accurate predictions of future storm impacts under changing environmental conditions. Process-based models offer an improvement over empirical models in that if the important underlying physics are understood and included in these models, they can be universally applied. In recent years advancements have been made in the development of process-based models for storm impact on sandy coasts (e.g., Roelvink et al., 2009;Tuan et al., 2006;van Thiel de Vries, 2009). In contrast, relatively few processbased models have been developed for gravel beaches. Due to the lack of measurement data collected under energetic to storm conditions on gravel beaches, existing process-based morphodynamic models for gravel beaches (e.g., Jamal et al., 2014;Pedrozo-Acuña et al., 2006;Van Rijn and Sutherland, 2011;Williams et al., 2012b) have been developed using data collected on natural or laboratory gravel beaches during low to moderate wave energy conditions, and may therefore not be representative of the physics and morphodynamics occurring during energetic storm events.
Additionally, while the initial results of existing process-based numerical models are promising, the validation of these models has thus far been limited to comparisons of morphological changes, rather than the hydrodynamic processes at the heart of the morphodynamic cycle. In particular, the implicit parametrisation (transfer of energy from wave-action balance to long waves; Jamal et al., 2014) and explicit parametrisation (effective onshore-directed swash velocity; Van Rijn and Sutherland, 2011) of the incident-band swash dynamics, which dominate wave run-up and overtopping motions on gravel beaches, require further validation before these models can be safely applied to simulate storm impacts on gravel beaches. Note that this point is also briefly referred to by Williams et al. (2012b), who applied both an incident-band phase-resolving shallow water wave approach, as well as an implicit wave-action balance type parametrisation, to simulate swash dynamics on a gravel barrier.
An accurate process-based model for storm impacts on gravel beaches would greatly increase the capacity of coastal managers to manage and plan for large storm events, such as those experienced in the UK in the winter of 2013-2014. Such a model could not only be used to provide early-warning of flooding events and assist emergency response coordination, but can also greatly improve the design of coastal defence structures and mitigation plans. The latter is particularly important when considering the large investments required in order to combat the potential effects of climate change and sea level rise on flooding (e.g., Environment Agency, 2009).
In this paper we attempt to address the need for a storm-impact model for gravel beaches by presenting a process-based model capable of simulating the hydrodynamics and hydrology on gravel beaches during storms. The model is validated using data collected during a large-scale physical model experiment (BARDEX; Williams et al., 2012a) and detailed in-situ field data collected at Loe Bar, Cornwall, UK (Poate et al., , 2014, as well as remote-sensed data collected at three other gravel beach locations along the UK coast during the 2012-2013 storm season, collected as part of the EPSRC-funded NUPSIG 1 project. The model is presented as a first step towards the development of a process-based morphodynamic model for storm impacts on gravel coasts.

Model description
In this paper we apply an existing open-source, process-based morphodynamic model for the nearshore and coast called XBeach 2 (Roelvink et al., 2009) to simulate the hydrodynamics on gravel beaches and barriers. The XBeach model has been shown to have good quantitative skill in hindcasting storm impact, overwash and breaching processes on sandy beaches Roelvink et al., 2009). Two modified versions of this model have been previously applied with relative success to model low wave-energy berm-building on a gravel beach (Milford-On-Sea; Jamal et al., 2014) and overwash on a gravel barrier (Slapton Sands; Williams et al., 2012b), although the physics included in the model were different in the two cases and both models effectively parametrised the incident-band wave run-up. In this paper, we use a onelayer, depth-averaged, non-hydrostatic extension to the XBeach model (Smit et al., 2010), similar to the SWASH model (Smit et al., 2013;Zijlema et al., 2011), that allows XBeach-G to solve wave-by-wave flow and surface elevation variations due to short waves in intermediate and shallow water depths. This is particularly important for application on gravel beaches, where due to steep slopes swash motion is mainly at incident wave frequencies, and infragravity wave motion, which dominates the inner surf and swash zone on sandy beaches during storms, is of secondary importance (e.g., Buscombe and Masselink, 2006). To correctly account for upper swash infiltration losses and exfiltration effects on lower swash hydrodynamics on gravel beaches, we compute groundwater dynamics and the exchange between groundwater and surface water using a newly developed groundwater model coupled to XBeach (McCall et al., 2012). Again, interaction between swash flows and the beach groundwater table are considered particularly important on gravel beaches due to the relatively large hydraulic conductivity of the sediment, while on sandy beaches this process is of significantly less importance (e.g., Masselink and Li, 2001).
In the following section we describe the central equations of the coupled surface water-groundwater model, which is termed XBeach-G in this paper (as in XBeach-Gravel). Although both surface water model and groundwater model are fully 2DH, in this paper we will restrict the description of the equations and application of the models to their 1D equivalent. We refer to Roelvink et al. (2009) andSmit et al. (2010) for a full description of the XBeach surface water model and its non-hydrostatic extension, and McCall et al. (2012) for a full description of the XBeach groundwater model.

Model coordinate system and grid
XBeach-G uses a coordinate system where the computational x-axis is orientated in the cross-shore direction, positive towards the coast, and a staggered grid system in which bed levels, surface water levels, groundwater levels, dynamic pressure, groundwater head and vertical fluxes are defined in cell centres, and horizontal fluxes are defined at cell interfaces. In 2DH applications, the computational grid may be curvilinear ; however, in this paper we apply only 1D rectilinear, non-equidistant grids. Since incident-band wave motions are resolved explicitly in the XBeach-G model, the grid resolution for an XBeach-G model is higher than for a regular XBeach model in which wave motions are computed on the wave group scale. In a 1D application of XBeach-G, this increase in model resolution leads to approximately 2-3 times greater computation times than a coarser resolution 1D XBeach model. The simulations presented in this paper have a simulation to computation time ratio of approximately 1:1-2:1 on a standard desktop PC, although higher ratios of approximately 5:1 can be achieved at the expense of detailed model resolution on the foreshore and backshore.
Surface water and groundwater dynamics are both computed using one layer in the vertical each, with the computational surface water layer located above the groundwater layer. Although the model equations are depth-averaged, two quasi-3D models are used to compute vertical velocities and pressures at the surface and bottom of the surface water and groundwater layers, in order to approximate the non-hydrostatic pressure distribution. The principal equations of the XBeach-G model are described in the following sections in their 1D-form.
Surface water flow is solved using a limited MacCormack (1969) predictor-corrector scheme that is second-order accurate in areas where the solution is smooth, and first-order accurate near discontinuities (Smit et al., 2010). The scheme is mass and momentum conserving following Stelling and Duinmeijer (2003), allowing for the correct representation of drying and flooding, as well as the capture of sub-and supercritical flows and shock-like features. Groundwater flow is solved using first-order central differences, which is considered sufficient to describe the inherently dissipative groundwater dynamics.

Surface water flow
Depth-averaged flow due to waves and currents are computed using the non-linear shallow water equations, including a non-hydrostatic pressure term and a source term for exchange with the groundwater: where x and t are the horizontal spatial and temporal coordinates respectively, ζ is the free surface elevation above an arbitrary horizontal plane, u is the depth-average cross-shore velocity, h is the total water depth, S is the surface water-groundwater exchange flux, v h is the horizontal viscosity, ρ is the density of water, q is the depth-averaged dynamic pressure normalised by the density, g is the gravitational constant and c f is the bed friction factor. Note that the exchange of horizontal momentum between the surface water and groundwater layer is assumed negligible. The horizontal viscosity (v h ) is computed using the Smagorinsky (1963) model to account for the exchange of horizontal momentum at spatial scales smaller than the computational grid size, which under assumption of longshore uniformity in flow and absence of longshore current is given as: where c s is the Smagorinsky constant, set at 0.1 in all model simulations and Δx is the computational grid size.
The bed friction factor (c f ) is computed using the Chézy equation for turbulent flow: where C ¼ ffiffiffiffiffiffiffiffi ffi 32g p log 12h k À Á is the Chézy bed friction coefficient and k is the characteristic roughness height. In this paper we assume k = 3D 90 in order to estimate bed friction in the swash zone on gravel beaches, which is the focus of this paper. Since this estimate of the roughness height is only valid for flat beds, the bed friction may be underestimated in the shoaling and surf zone. The depth-averaged normalised dynamic pressure (q) is derived in a method similar to a one-layer version of the SWASH model (Zijlema et al., 2011), in which the depth-averaged dynamic pressure is computed from the mean of the dynamic pressure at the surface and at the bed, assuming the dynamic pressure at the surface to be zero and a linear change in the dynamic pressure over depth. In order to compute the normalised dynamic pressure at the bed, the contributions of advective and diffusive terms to the vertical momentum balance are assumed to be negligible: where w is the vertical velocity and z is the vertical coordinate. The vertical velocity at the bed is set by the kinematic boundary condition: where ξ = ζ − h is the elevation of the bed and the subscript b refers to the location at the bed.
Combining the Keller-box method (Lam and Simpson, 1976) as applied by Stelling and Zijlema (2003) for the description of the pressure gradient in the vertical and Eq. (5), the vertical momentum balance at the surface can be described by: where the subscript s refers to the location at the surface. The dynamic pressure at the bed is subsequently solved by combining Eq. (7) and the local continuity equation: Smit et al. (2010) have shown that the inclusion of the dynamic pressure described above reduces the relative dispersion and celerity errors in the non-linear shallow water equations of XBeach to less than 5% for values of kh ≤ 2.5 and allows for accurate modelling over wave transformation on dissipative beaches. In order to improve the computed location and magnitude of wave breaking, we apply the hydrostatic front approximation (HFA) of Smit et al. (2013), in which the pressure distribution under breaking bores is assumed to be hydrostatic. Following the recommendations of Smit et al. (2013), we consider waves to be hydrostatic bores where δζ δt N0:6 and to reform if δζ δt b0:3 . Although this method greatly oversimplifies the complex hydrodynamics of plunging waves on gravel beaches, we show in this paper that the application of this model provides sufficient skill to describe dominant characteristics of the flow, without requiring computationally-expensive high-resolution discretisation of the vertical and surface tracking of overturning waves.

Groundwater flow
Horizontal groundwater flow in the aquifer is computed assuming incompressible flow and the Law of Darcy (1856): where u gw is the depth-averaged horizontal groundwater velocity, h gw is the height of the groundwater surface above the bottom of the aquifer, w gw,s is the vertical groundwater velocity at the groundwater surface, K is the hydraulic conductivity of the aquifer and H is the depth-averaged hydraulic head. Note that the bottom of the aquifer is assumed impermeable and the vertical groundwater velocity at the bottom of the aquifer is zero. Since Darcy's Law is only strictly valid for laminar flow, we approximate turbulent groundwater flow conditions using a modification of the laminar hydraulic conductivity similar to Halford (2000): where K lam is the laminar hydraulic conductivity, Re crit is the critical Reynolds number for the start of turbulent flow, Re ¼ u gw j jD50 nν is the Reynolds number of the groundwater flow in the pores of the aquifer, D 50 is the median grain size, v is the kinematic viscosity of water and n is the porosity.
In order to compute the non-hydrostatic groundwater pressure, the groundwater head is approximated by a parabolic curve in the vertical, which is bound by a zero vertical velocity condition at the impermeable bottom of the aquifer, the imposed head at the groundwater surface, and an assumption of a constant gradient in the vertical groundwater velocity over the vertical (McCall et al., 2012): in which H is the groundwater head, varying over the vertical, σ is the vertical coordinate above the bottom of the aquifer, β is the parabolic curvature coefficient and H bc is the head imposed at the groundwater surface. In the case of hydrostatic pressure, β reduces to zero. The depth-average value of the groundwater head is found by integrating Eq. (12) over the groundwater column: The vertical velocity at the groundwater surface is computed from the gradient of Eq. (12) at the surface and the hydraulic conductivity: Eqs. (9) and (10) form a coupled system that is solved by substitution of Eqs. (13) and (14). The groundwater level is subsequently computed as: 2.4. Groundwater-surface water exchange The groundwater and surface water are said to be in a connected state where and when the groundwater level reaches to the top of the bed and surface water exists above the bed. In this case the rate of exchange between the surface water and groundwater, defined positive from surface water to groundwater, is determined by the vertical groundwater velocity at the interface between the groundwater and surface water: Infiltration and exfiltration occur in locations where the groundwater and surface water are not connected. Infiltration takes place where surface water covers an area in which the groundwater level is lower than the bed level. The flux of surface water into the bed is related to the pressure gradient across the wetting front in a manner similar to Packwood (1983): where p b = ρ(q b + gh) is the total surface water pressure at the bed and d i is the thickness of the wetting front, which increases over time during the infiltration event according to the infiltration velocity: Since the groundwater model has one vertical layer and cannot track multiple layers of groundwater infiltrating into the bed, the wetting front thickness is reset to zero when the surface water cell becomes dry, or the groundwater and the surface water become connected. All infiltrating surface water is instantaneously added to the groundwater volume, independent of the distance from the bed to the groundwater table. Since the groundwater model neglects the time lag between infiltration at the beach surface and connection with the groundwater table a phase error may occur in the groundwater response to swash dynamics. However, this phase error is expected to be small on permeable gravel beaches where the distance between the waterline and the groundwater table is generally small, as also shown by McCall et al. (2012), and does not affect the modelled infiltration velocities at the beach face.
Exfiltration occurs where the groundwater and surface water are not connected and the groundwater level exceeds the bed level. The rate of exfiltration is related to the rate of the groundwater level exceeding the bed level: 3. Measurement data and model setup The data used in this paper to set-up and validate the XBeach-G model have been collected during the BARDEX large-scale physical model experiment in the Deltaflume, The Netherlands, and at four gravel beach locations along the coast of the UK as part of the NUPSIG project ( Fig. 1 and Fig. 2).
During the BARDEX physical-model experiment, the hydrodynamics and morphodynamics of a 4-metre high and 50-metre wide gravel barrier were measured under varying hydraulic boundary conditions, ranging from wave run-up to wave overtopping and overwash (see Williams et al., 2012a for details). The morphodynamic response of the gravel barrier was measured by a mechanical roller and actuator following the bed profile from an overhead carriage before and after each 5-20-minute wave sequence. Measurements of the groundwater head in the gravel barrier were made using 15 pressure transducers buried in the bed beneath the gravel barrier. Wave transformation across the foreshore was measured using three wave gauges located c. 40 m offshore of the beach and one nearshore pressure transducer near the toe of the gravel beach. Wave run-up and overtopping levels were measured using a cross-shore array of 45 acoustic bed level sensors (BLS; cf. Turner et al., 2008) that spanned the entire subaerial portion of the gravel barrier. Poate et al. (2013) collected in-situ and remote-sensed hydrodynamic and morphodynamic data on a fine gravel barrier (Loe Bar, Cornwall, UK) over a period of four weeks. Two energetic events occurred during this period on 8 March 2012 (LB1) and 24 March 2012 (LB2) with offshore significant wave heights of 1.6-2.3 m. Offshore wave conditions were measured by a directional wave buoy in 15-20 m water depth maintained by the Channel Coastal Observatory (CCO). Tide and surge levels were measured by a pressure transducer located in Porthleven harbour, approximately 2 km from the field site ( Fig. 1, Site B). Wave transformation across the beach face was measured by a cross-shore array of five pressure transducers (PTs), as shown in Fig. 3. Wave run-up time series were extracted from water level and bed level data collected by a crossshore array of 45 BLS that spanned the beach face from MHWS-level to the barrier crest. Bed levels along the main instrument cross-shore transect were measured every low tide using Real Time Kinematic GPS (RTK-GPS). During wave event LB2, high-frequency (2 Hz) and horizontal resolution (0.05-0.20 m) bed level and water level data were collected continuously by a tower-mounted cross-shore laser scanner (cf., Almeida et al., 2013).
Wave run-up data were collected during storm conditions in the winter of 2012-2013 at three gravel beaches along the coast of the UK (Loe Bar, Slapton Sands and Chesil Beach) and one composite beach with a gravel upper beach fronted by a sandy low-tide terrace (Seascale) ( Fig. 1 and Fig. 2). At Loe Bar and Slapton Sands, offshore wave data were provided by directional wave buoys maintained by CCO, located approximately 500 m from the study site in 15-20 m and 10-15 m water depth respectively. At Chesil Beach offshore wave data were provided by a directional wave buoy maintained by CCO, located approximately 7 km from the study site in 12-15 m water depth. Wave data at Seascale were obtained from an offshore wave buoy maintained by the Centre for Environment, Fisheries & Aquaculture Science (CEFAS), located 50 km from the study site in 20 m water depth, supplemented by a nearshore PT in 0-4 m water depth (depending on tide), which is used to scale the offshore wave height to account for offshore wave refraction and sheltering. Tide and surge data at Loe Bar, Slapton Sands and Seascale were provided by PTs located approximately 2 km, 1 km and 500 m from the main instrument transect lines, respectively, whereas tide data at Chesil Beach were derived from tidal predictions. At all four gravel beach field sites, shoreline position time series were extracted along 4-6 cross-shore transects from digital video camera pixel time stacks collected at 3.75 Hz, following the method described by Poate et al. (2013). Pre-and post-event bed levels were measured using RTK-GPS at all four locations.
A summary of the measured, or estimated, median grain diameter (D 50 ), hydraulic conductivity (K) and beach slope (tan(β)) at all four gravel beach field sites and the BARDEX experiment, as well as a summary of the instruments used to collect hydrodynamic and morphodynamic data presented in this paper, is given in Table 1. The table furthermore lists the maximum hydrodynamic forcing conditions (offshore significant wave height (H m0 ), peak wave period (T p ), relative freeboard (R c /H m0 ) and wave angle relative to shore normal (θ rel )) during each of the storm events discussed in this paper, and an overview of the use of the measurement data in the validation of the XBeach-G model.

Model validation
This section describes the comparison of the model simulation data and data collected during the BARDEX experiment and the field measurements at the four UK gravel beaches. The model results are split into five categories: (1) groundwater dynamics, (2) wave transformation, (3) wave set-up, (4) wave run-up, and (5) wave overtopping. For comparison with the measurements, cross-shore transect models are set up in XBeach-G for all five gravel beaches (BARDEX, Loe Bar, Chesil Beach, Slapton Sands and Seascale). In each model, the bed level is set to the bed level measured along the main instrument array (for Loe Bar simulations LB1 and LB2, and BARDEX), or along the main crossshore video image pixel time stack transect (for Loe Bar simulation LB3, Slapton Sands, Chesil Beach and Seascale) for the low tide prior to the storm event. The models for Loe Bar, Chesil Beach, Slapton Sands and Seascale are forced using wave spectra measured at the nearest wave buoy, described in the previous section, and measured (Loe Bar, Slapton Sands, Seascale) or predicted (Chesil Beach) tide and surge levels. The XBeach-G model uses the input wave spectrum to generate a random time series of incident waves and bound low-frequency second order waves at the model boundary. In the BARDEX simulations, measured time series of the water elevation at the flume wave generator are used to force the XBeach-G model. The hydraulic conductivity of the beach used by the groundwater component of the XBeach-G model and grain size properties are derived from in-situ measurements, literature or estimates (Table 1). The cross-shore resolution of the models is set to vary gradually in the cross-shore direction, from Lm 25 ≈2-3 m at the offshore boundary of the model, where L m is the wave length related to the mean wave period, to 0.1 m near the waterline in order to correctly capture wave breaking and wave run-up in the model. In the case of the BARDEX simulations, the resolution has been increased to 0.25 m at the wave generator and 0.05 m at the beach.
Since not all types of measurement data are available at all five beaches, the validation of the model results will be restricted to certain datasets, as outlined in Table 1. Multiple simulations are carried out at all five gravel beaches, representing periods of different wave or tidal forcing.

Groundwater dynamics
Groundwater pressure data collected at 13 of the buried PTs during one measurement series of the BARDEX experiment (B-E10; Table 1), which was characterised by significant overwash activity, are used to validate the groundwater component of the model. For this comparison, groundwater pressure at the PTs is converted to groundwater head as H ¼ p ρg . The groundwater component of the model is initialised using the measured groundwater head at the start of the series. The surface water component (see following sections) provides all the boundary conditions for the groundwater component of the XBeach-G model during the simulation. For modelling purposes, the gravel barrier is assumed to be homogeneous, with a constant hydraulic conductivity of 0.16 m s −1 and porosity of 0.32, based on analysis by Turner and Masselink (2012).
Time series of measured and modelled groundwater head at four PTs under the gravel barrier, as well as the groundwater head variance density at the four PTs, are shown in Fig. 4. The figure shows a gradual decrease in the measured amplitude of the groundwater variation and a transition of the peak of the groundwater head variance from the incident wave frequency (0.13 Hz) to lower frequencies, with increasing distance from the front to the back of the gravel barrier. Both phenomena are well described by the groundwater model. The groundwater head at all 13 buried PTs is simulated with a combined RMSE of 0.064 m and bias of −0.012 m (Table 2), which is appropriate for the application purpose of this model, especially considering the inherent uncertainties in the hydrological and geotechnical properties of gravel barriers. An earlier study using the groundwater model in XBeach provided similar accuracy in the prediction of the groundwater head for three other measurement series of the BARDEX experiment (McCall et al., 2012).

Wave transformation
Wave transformation from offshore to the gravel barrier toe and the lower swash is compared in the model to data collected during the BARDEX experiment and to data collected at Loe Bar. In the BARDEX experiment, surface water pressure was measured by a shallow water PT near the toe of the barrier (bottom panels in Fig. 5). In this analysis, we convert the pressure measurements at the toe of the barrier to surface elevation time series using the local approximation method of Nielsen (1986). XBeach-G is used to simulate the wave transformation during two measurement series of the BARDEX experiment with Table 1 Overview of beach characteristics and data collection instrumentation at each of the gravel beach sites and maximum hydrodynamic forcing conditions and model-data comparisons for each of the storm events. In model-data comparison, measurement data are used to validate groundwater dynamics (G), wave transformation ( Heijne and West (1991). In the case of Loe Bar § was determined for the beach face by in-situ falling head tests, and at Seascale, † is estimated for the gravel section of the beach and ‡ is relative to top of gravel beach.  Table 1). In these simulations, the model is forced at the offshore boundary using time series of the water surface elevation measured at the wave-maker and an estimate of the intra-wave depth-average cross-shore velocity at the boundary based on linear wave theory.
To validate the transfer of the incident-band wave energy to higher and lower harmonics across the barrier foreshore, the wave spectrum at the model boundary is compared to the computed and measured wave spectrum at the location of the shallow water PT (Fig. 5). The figure shows a transfer of wave energy from the peak of the wave spectrum (0.23 Hz) to lower frequencies (0.05 Hz) in B-C2 (left panels) and from the peak of the spectrum (0.13 Hz) to lower frequencies (0.02 Hz) as well as higher frequencies (0.25 Hz and 0.36 Hz) in B-E10 (right panels), representing the transfer of energy to higher and lower harmonics of the peak frequency band. The results of the model simulations show that XBeach-G is capable of reproducing this transfer across the frequency band relatively well, although the energy in the upper and lower frequency bands appears to be under-predicted somewhat by the model. Since the measured surface water elevation at the toe of the gravel barrier contains both incident and reflected waves, this under-prediction may be both due to lower energy transfer rates in the incident wave components, as well as an incorrect representation of the amplitude or phase of the reflected wave components. Despite the under-prediction in the high and low frequency components, the overall spectral significant wave height at the shallow water PT is predicted well by the model with an RMSE of 0.034 m (4.6%) and 0.028 m (3.1%) in the B-C2 and B-E10 measurement series, respectively. Comparable model results and accuracy were found for the simulation of B-BB1 and B-C1 (not shown).
To determine whether the model is also capable of predicting wave transformation well on natural beaches, the XBeach-G model is used to simulate wave transformation at Loe Bar. During this field experiment, five PTs were mounted near bed level to a cross-shore scaffold instrument frame spanning the upper inter-tidal (see Fig. 3 for an overview of the location of the PTs and Poate et al. (2013) for further details). As in the case of the BARDEX pressure data, water surface elevation time series were derived from the measured pressure data using the local approximation method of Nielsen (1986). An XBeach-G model was set up for two high-energy wave events on 8 March 2012 (LB1) and 24 March 2012 (LB2) with offshore significant wave heights 1.6-2.3 m, as discussed in Section 3. The XBeach-G model is forced using directional wave spectrum time series measured by the CCO  A comparison of measured and modelled wave heights, split into high-, mid-and low-frequency components at the five cross-shore PTs at Loe Bar is shown in Fig. 6 for LB1, and in Fig. 7 for LB2. Fig. 6 shows that for LB1, little wave height transformation takes place between the nearshore wave buoy and the most seaward pressure transducer (PT9), except for an increase in the low-frequency band. The wave height in the high-frequency band gradually decreases in the cross-shore direction between PT9 and PT6, whereas the wave height in the mid-frequency band shows relatively little decay compared to LB2, which is likely due to the slightly reflective state of the beach for the long-period waves of LB1 (Table 1). Note that water depths at PT5 during LB1 are too small to compute wave statistics during any part of the tide. Fig. 7 shows a strong increase in the measured low-frequency wave height from the offshore boundary of the model to the most seaward pressure transducer (PT9) in LB2. During this event, wave heights in the mid-and high-frequency components of the wave spectrum are generally lower at PT9 than offshore. In the cross-shore direction, all measured wave heights are modulated by the tide level. Both figures show that wave heights in the low-, mid-and high-frequency bands are generally predicted well in the model. In contrast to the results of the BARDEX simulations, the highand low-frequency components of the wave spectrum are slightly overpredicted during the LB1 and LB2 (positive bias), instead of underpredicted. During LB2, the accuracy of the model predictions of the wave height decreases over time at the most landward pressure transducers (in particular PT5, PT6 and PT7), which may be due to the lack of morphological updating in the model. Notwithstanding these errors, the quantitative model skill in predicting wave height transformation across the foreshore and gravel beach is good, with RMSE in the high-, midand low-frequency band b0.24 m for LB1 and b 0.30 m for LB2, which is approximately 15% and 13% of the total offshore wave height of the two wave events, respectively. The SCI of the model wave height prediction is low (SCI b 0.26) for all frequency bands at the two most offshore pressure transducers (PT8 and PT9), and reasonable (SCI b 0.57) at the three landward pressure transducers (PT5, PT6 and PT7). The overall RMSE for the integrated wave height is 0.11 m during LB1 and 0.28 m during LB2, corresponding to a SCI of 0.14 and 0.21, respectively ( Table 3).
The evolution of the wave spectrum from offshore to the five crossshore PTs is shown in Fig. 8 at four stages of the tide during LB2. The figure shows a distinct drop in wave energy at the peak of the spectrum across the PT array, caused by depth-induced wave breaking, and transfer of wave energy to lower and higher harmonics of the peak frequency band. Both phenomena are represented well by the XBeach-G model, indicating that the model skill is not restricted to ensemble wave heights and the total wave energy, but can also be used to study wave spectrum transformation on gravel beaches.
Finally, the transformation of the wave shape is examined in terms of wave skewness (Sk) and wave asymmetry (As). In this analysis, both parameters are computed from a low-pass (f ≦ 5f p ) filtered time series, where f p represents the offshore spectral peak frequency, of the modelled water surface elevation and the water surface elevation derived from the measured pressure time series (ζ lpf ) as follows: As Modelled and measured wave skewness and wave asymmetry at the five cross-shore PTs are shown in Fig. 9. The figure shows that wave skewness and asymmetry are predicted relatively well by the model at the two most offshore pressure transducers (PT8 and PT9), but that in general wave asymmetry is slightly overpredicted by the model, particularly at the three most shoreward pressure transducers (PT5, PT6 and PT7). The overprediction of the wave asymmetry in the model may be the result of the simplified method in which the model attempts to simulate the complex hydrodynamics of breaking waves using the hydrostatic front approximation (HFA), as also found to lesser extent in the SWASH-model under narrow-banded wave conditions (Smit et al., 2014). However, since wave skewness and asymmetry are sensitive to water depth, changes in the wave asymmetry due to errors in the imposed bed level may also contribute to the differences found between the model and measurements. It should be noted that although the error in the predicted wave shape is sufficiently small for the purpose of the current hydrodynamic model, differences in the wave skewness and asymmetry may have more adverse consequences if the model is used to compute sediment transport and morphology.

Wave set-up
Steady wave set-up at the five cross-shore PTs at Loe Bar is extracted from the measured pressure records for LB1 and LB2 by subtraction of the tide and surge level measured by the harbour tide gauge, from 15minute averaged water levels measured at the PTs. Time series of the steady wave set-up for both wave events are shown in Fig. 10. The figure shows little measured wave set-up at the most offshore cross-shore pressure transducer (PT9), where set-down dominates during LB1, and set-up is less than 0.5 m during LB2. For both events, wave set-up increases in shoreward direction across the PT-array, and reaches a minimum at all PTs at high tide (16:45 and 18:00 for LB1 and LB2, respectively). Wave set-up at all cross-shore PTs is predicted reasonably well, with RMSE b0.10 m (approximately 6% of the tidal amplitude) for LB1 and b 0.25 m (approximately 13% of the tidal amplitude) for LB2. The larger error in the steady wave set-up during LB2 than LB1 is primarily due to an under-estimation (negative bias) of the measured wave set-up at the most landward pressure transducers (PT5, 0.25 m; PT6, 0.20 m). This may partly be explained by the lack of morphological updating in the model, also noted in the discussion of the wave height transformation in Fig. 7, and is addressed in Section 5.1. It should be noted that although the SCI is included in Fig. 10 for reference, the value at the most seaward pressure transducers (PT8 and PT9) are poor in the case of LB1 due to the very low value of the denominator in the SCI calculation, rather than to particularly large errors in the predictions.

Wave run-up
Data on wave run-up levels were collected using a cross-shore array of bed-level sensors during the BARDEX experiment (Table 1;  Significant offshore wave height time series (*), significant wave height time series measured by five nearshore pressure transducers (○) and significant wave height time series modelled at the location of the nearshore pressure transducers (□) during LB1, separated into three frequency bands, where f p represents the offshore spectral peak frequency. The locations of the five nearshore pressure transducers in the cross-shore are shown in Fig. 3.

Sands (SS) and Seascale (SE).
For the purpose of this study, the shoreline derived from the pixel time stacks is assumed to correspond to a water depth of 0.01 m, and this value is used as a depth criterion to determine the shoreline time series in the bed-level sensor data and XBeach-G model results. The 2%, 5%, 10% and 20% (R 2% , R 5% , R 10% , R 20% ) run-up exceedence levels are computed from 15 to 20-minute sections of the shoreline time series above Still Water Level (cf., Stockdon et al., 2006).
To compare predicted and measured run-up levels, XBeach-G models are set up for the three measurement series of the BARDEX experiment and the six storm events discussed above (cf. Table 1). Each BARDEX series simulation is run for one measured wave sequence of approximately 20 min. In the case of the storm events, one simulation is run for every 1-3 sequential daytime high-tides of the storm event. Each high-tide simulation is run for the duration of maximum tide levels and contiguous camera or bed-level sensor data, which was generally in the order of 0.5-1 h. Run-up exceedence levels are computed from the modelled shoreline time series using identical methods and computation periods as used in the derivation of the measured run-up levels. To investigate the sensitivity of the modelled run-up levels to the selection of random wave components at the model boundary, each XBeach-G simulation is run ten times using a new random wave time series of the imposed offshore wave spectrum.
Mean measured and modelled run-up levels computed for every 15-20-minute section of shoreline time series data at all sites are Significant offshore wave height time series (*), significant wave height time series measured by five nearshore pressure transducer (○) and significant wave height time series modelled at the location of the nearshore pressure transducers (□) during LB2, separated into three frequency bands, where f p represents the offshore spectral peak frequency. The locations of the five nearshore pressure transducers in the cross-shore are shown in Fig. 3. shown in Fig. 11. Vertical error bars in the figure represent variations in the modelled run-up levels due to variations in the random wave times series applied at the model boundary. Horizontal error bars represent the variation in measured run-up data across the multiple cross-shore camera pixel stacks, cf. Poate et al. (2014). The figure shows very good correspondence and little scatter between measured and modelled runup levels for all exceedence probabilities and at all five gravel beaches. Importantly, the model shows practically no systematic relative bias (defined as the absolute bias, normalised by the measured run-up) in the computation of the extreme run-up levels, and only a very small negative bias (under-prediction) of the 10% and 20% run-up exceedence levels.
Variations in modelled and measured run-up levels due to variations in the imposed wave time series and cross-shore camera pixel stack locations are up to 1 m (20%) for run-up levels over 5 m. The model skill is further examined in Fig. 12, which presents histograms of the absolute relative error in the mean run-up level prediction for all 15-20-minute sections of shoreline time series at all five gravel beaches. The figure shows that the majority of absolute relative run-up level prediction errors are in the order of 0-15% of the measured runup. The empirical relative error exceedence function in the same figure shows that the median (50% exceedence) relative error for R 2% is less than 10%, and the maximum relative run-up error for R 2% is 29.4%. These values indicate that even without morphological updating, the model can potentially be applied to investigate extreme run-up levels and the possibility of wave overtopping under energetic wave conditions.

Wave overtopping
Time series of overtopping waves were measured by a cross-shore array of 45 bed-level sensors during the BARDEX experiment. Data provided by these instruments are the level of the bed directly below the ultrasonic sensor (when the bed is dry), or the water level below the ultrasonic sensor (when the bed is covered with water). To study the applicability of the XBeach-G model to predict overtopping waves on gravel barriers, XBeach-G simulations are set up of BARDEX measurement series B-E9 and B-E10, during which wave overtopping of the barrier crest took place. Due to lowering of the crest during the experiment, the relative freeboard of the barrier is higher in B-E9 than in B-E10, which in combination with a slight change in the beach slope results in more overtopping waves in B-E10 than in B-E9. Since considerable bed level change occurred during both measurement series, the XBeach-G simulations are limited to the first 10 min of overtopping waves during which the crest level was lowered less than 0.15 m from the level at the start of each series.
Comparisons of modelled and measured time series of the bed level and water level at three locations across the gravel barrier in B-E9, and four locations in B-E10 shown in Figs. 13 and 14 respectively. Data at the most landward sensor (BLS45) are not shown in the comparison of B-E9 due to the lack of reliable measurement data. The figures show a reduction in the number of waves, described by spikes in the time series, and their amplitude, from the most seaward sensor (BLS30) to the most landward sensors (BLS40 and BLS45). This reduction in the number and the size of overtopping waves is due to infiltration of the swashes on the back barrier. Periods in which the dry bed is measured by the sensor are indicated by the horizontal sections in the time series between waves. The measurements of the dry bed show that the bed at BLS30 erodes approximately 0.15 m in the first 10 min of B-E9 and B-E10, and that some accretion takes place at BLS40 in both series. Figs. 13 and 14 show that the XBeach-G model is able to reproduce the time series of overtopping waves at most locations on the gravel barrier well. At the locations of BLS30 and BLS35, the model correctly  Fig. 8. Offshore wave spectra (black dashed), wave spectra measured by five nearshore pressure transducers (black) and wave spectra modelled at the location of the nearshore pressure transducers (orange) at 15:00 (first row), 16:00 (second row), 17:00 (third row) and 18:00 (fourth row) at Loe Bar during LB2. Note that PT9 did not record any data at 18:00. The locations of the five nearshore pressure transducers in the cross-shore are shown in Fig. 3.
predicts more than 78% of the overtopping wave occurrences that exceed the initial bed level (Table 4). The wave height of the majority of these overtopping events is also predicted well by the model, although accuracy of the wave height predictions at BLS30 is strongly reduced by the erosion of the bed. Wave overtopping at BLS40 is poorly predicted by the model during B-E10, where only 28% of overtopping waves are correctly reproduced by the model, however at BLS45 in the same series, the model skill improves by correctly predicting the four largest of six overtopping wave events. The reason for the improvement in the model skill from BLS40 to BLS45 is not clear. However, the approximation of the infiltration velocity in the groundwater component of the XBeach-G model, the lack of morphological updating in the XBeach-G model, and possible longshore non-uniformities in the barrier response of the BARDEX physical model, may all be considered sources for discrepancies between the measurements and modelled results.
The results of simulations B-E9 and B-E10 show that the XBeach-G model is well capable of predicting initial wave overtopping at the crest of the gravel barrier. The model also correctly predicts the evolution of most initial overtopping waves across the back barrier. These results show that the model may be considered a useful tool with which to estimate the potential for overtopping on gravel barriers. However, since much bed level change is expected during overtopping and overwash events, the addition of morphodynamic feedback in XBeach-G is considered necessary in order to properly predict the development of overtopping and overwash discharge during these events.

Effect of morphological updating on computed wave setup, wave transformation and wave run-up
As discussed in Section 4.3 and shown in Fig. 10, application of the XBeach-G model to LB2 underestimates the measured wave set-up at the most landward pressure transducers (PT5 and PT6) by as much as 0.35 m. This underestimation of the set-up is mainly attributed to the absence of morphodynamic updating in the XBeach-G model, specifically ignoring the fact that the high tide beach morphology is significantly different from that during low tide. Here, the effect of including morphodynamic updating on the predicted set-up is investigated. The mean bed level position during LB2 is derived every 15 min from high-frequency (2 Hz) laser data along the model cross-shore transect from the wave run-down level to the barrier crest (cf. Almeida et al., 2013). The model is then re-run using the laser-derived time series of the bed level elevation as a time-varying bed boundary condition. Note that because no laser-derived bed elevation data exist below the wave run-down level, we assume for the purpose of this sensitivity analysis that the bed level below the wave run-down level remained constant during the event.
The results of the wave set-up predicted by the XBeach-G model for LB2 with, and without, measured bed level updating are shown in the top panel of Fig. 15 at the moment of maximum wave set-up (18:00). The figure shows that the build-up of the gravel step (derived from The bottom panel of Fig. 15 shows the measured wave height at PT5-PT8 at 18:00, as well as the wave height computed at the same PTs by the XBeach-G simulations with, and without, bed level updating. The figure shows that although bed level updating does modify the computed wave height, the model prediction of the nearshore wave height is less sensitive to bed level updating than the computed setup. Wave heights computed by the XBeach-G simulation with bed level updating are 5-10% lower than those in the XBeach-G simulation without bed level updating, leading to lower model bias and RMSE, and slightly lower SCI values in the simulation with bed level updating. The overprediction of the wave asymmetry at PT5-7 discussed in Section 4.2 is not reduced significantly by the application of bed level updating in the model, indicating that a modification of the HFA-model may be necessary, alongside more accurate bed level information below the waterline, in order to correctly predict the wave skewness and asymmetry in the lower swash and inner surf zone. The application of bed level updating in the XBeach-G model affects the computed wave run-up levels to a similar magnitude as the wave height (not shown in Fig. 15). In the case of run-up however, the computed value is 5-10% higher in the simulation with bed level updating compared to the simulation without bed level updating, leading to slightly better predictions of the maximum run-up extent during LB2. The increase in the run-up height is explained to a great extent by the large increase in the nearshore wave set-up, in combination with relatively little wave height reduction, in the model simulation with bed level updating relative to the model simulation without bed level updating. This model observation appears contrary to measurement data presented by Poate et al. (2013), who show a reduction in the run-up height due to the development of a step during LB2. This difference between the observed and modelled behaviour may indicate a limitation of the XBeach-G model, but may also be the result of the lack of updated bed level information below the wave run-down level. It should also be noted that the difference in run-up height between both model simulations is of the same order as the model prediction error and the natural spread in run-up heights described in Section 4.4. Due to the absence of overtopping during LB2, the effect of the morphological development of the beach on overtopping discharge is not examined. However, the effect of the morphodynamic changes on the beachface are expected to be small compared to the effect of morphodynamic changes at the beach or barrier crest during overtopping and overwash events.

Effect of non-hydrostatic wave and groundwater model components
The version of the XBeach-G model discussed in this paper has been modified from the standard version of XBeach for sandy coasts (e.g., Roelvink et al., 2009) through two extensions to the XBeach model: (1) the application of a non-hydrostatic pressure correction term (Smit et al., 2010) that allows wave-by-wave modelling of the surface elevation and depth-averaged flow due to the incident-band short waves, instead of the use of the standard wave-action balance (surf beat) approach to model short waves; and (2) the application of a groundwater model (McCall et al., 2012) that allows infiltration and exfiltration through the permeable gravel bed to be simulated.
To study the effect of these two extensions for the purpose of simulating storm impact on gravel barriers, and to access the improvement to the model performance, we re-simulate B-E9 and B-E10, previously discussed in Section 4.5, with two variants of the XBeach-G model. In the first variant groundwater interactions are included, but the short waves are modelled using the standard wave-action balance approach (Variant 1). In the second variant waves are modelled using the nonhydrostatic wave-by-wave approach, but groundwater interactions are excluded (Variant 2). The models are all forced using the same model grid and wave boundary condition information. However, since Variant 1 uses a wave-action balance approach to model the incident waves, the total incident wave signal for this Variant is split into a high-frequency wave energy part (f ≥ 0.5f p ) varying on the wavegroup time scale, which is used as a boundary condition for the waveaction balance, and a low-frequency flux component (f b 0.5f p ) that is imposed as a boundary condition in the hydrostatic non-linear shallow water equations (cf. Roelvink et al., 2009;van Thiel de Vries, 2009). Measured wave overtopping time series, and wave overtopping time series modelled by XBeach-G and the two Variants are shown in Figs. 16 and 17 for B-E9 and B-E10, respectively.
The results of the simulations using Variant 1 show that the simulation of the incident waves using the non-hydrostatic wave-by-wave method greatly increases the model skill in predicting overtopping waves compared to the wave-action balance method. This effect is particularly clear in the case of B-E9 (Fig. 16), in which Variant 1 does not predict any of the 63 wave overtopping events at the crest of the gravel barrier (BLS30), whereas the XBeach-G version correctly predicts 90% of the overtopping wave events (see Table 4). The improvement of the XBeach-G model over Variant 1 is less pronounced in the case of B-E10 (Fig. 17), which has a lower relative freeboard than B-E9, causing almost every wave to overtop. In this simulation, Variant 1 predicts wave overtopping events that are generally lower in amplitude than the measured overtopping events, and have a duration of several incident waves, corresponding to low-frequency wave-group motions. Interestingly, the model skill of Variant 1 is comparable to that of XBeach-G at the back of the gravel barrier (BLS45) Fig. 13. Measured (black) and modelled (orange) time series of overtopping waves (spikes) and bed levels (horizontal sections) at the locations of three acoustic bed level sensors during B-E9. The locations of the bed level sensors relative to the barrier profile (black) and water levels (grey) are shown in the bottom panel. Note that the gradual erosion at BLS30 and BLS35, and accretion at BLS40 found in the measurements is not accounted for in the XBeach-G model. Sparse data were collected at BLS45 due to the proximity of the water level to the instrument and are therefore not shown in this figure. that the large swash events that reach the back of the barrier are related to low-frequency motions on the wave group time scale. Simulations of B-E9 and B-E10 using Variant 2 show that the inclusion of infiltration and exfiltration through the groundwater component does not significantly alter the prediction of overtopping waves at the barrier crest (BLS30 and BLS 35 in B-E9 and B-E10, respectively). However, the XBeach-G model shows substantially better model skill in predicting overtopping time series at the back of the barrier compared to Variant 2. In these locations, Variant 2 greatly overpredicts the number, and the magnitude, of the overtopping swash events compared to the measured time series, whereas XBeach-G correctly predicts 74% and 67% of the overtopping swashes for B-E9 and B-E10 respectively (see Table 4).  From these results it can be concluded that the non-hydrostatic wave-by-wave modelling of the incident wave field is necessary to predict run-up levels and the start of overtopping on gravel beaches, and can only partially be replaced by a wave-action balance approach in case of very low relative freeboards and large infragravity motions. Groundwater interaction is required in order to correctly model the Cross shore distance (m) Fig. 15. Effect of bed level updates on computed mean water levels and wave heights at Loe Bar. Top panel: measured 20-minute mean water at 18:00 during LB2 (□), modelled mean water level without updated bed levels (black solid line) and modelled mean water level with updated bed levels (black dashed line). Bed levels corresponding to the period of the simulation without and with updated bed levels are indicated by the solid and dashed grey lines respectively. Bottom panel: measured wave height at 18:00 during LB2 (□), modelled wave height without updated bed levels (○) and modelled wave height with updated bed levels (grey △).  evolution of overtopping waves across the gravel barrier. It should be noted that the gravel barrier in the BARDEX experiment is exceptionally thin, and that even larger prediction errors will occur on wider gravel barriers if groundwater interaction is excluded. Additionally, the importance of accounting for groundwater interactions becomes increasingly important as the hydraulic conductivity (i.e., sediment size) of the barrier material increases.

Conclusions
This paper presents an extension of the XBeach numerical model to simulate hydrodynamics on gravel beaches under energetic wave conditions. The model is modified from the standard XBeach model for sandy beaches by the inclusion of (1) a non-hydrostatic pressure correction term (Smit et al., 2010) that allows wave-by-wave modelling of the surface elevation and depth-averaged flow, and (2) a groundwater model (McCall et al., 2012) that allows infiltration and exfiltration through the permeable gravel bed to be simulated, and is referred to as XBeach-G. The model does not include sediment transport formulations for gravel and is therefore run without morphodynamic updating. The XBeach-G model is applied to simulate groundwater dynamics, wave transformation, wave set-up, wave run-up and wave overtopping on one large-scale physical model dataset (BARDEX; Williams et al., 2012a), one large-scale field experiment dataset (Loe Bar; Poate et al., 2013) and storm wave run-up measurements at three gravel beaches and one composite beach with a gravel upper beach fronted by a sandy low-tide terrace. A comparison between modelled and measured hydrodynamics shows that the model is capable of reproducing groundwater dynamics, wave height, wave spectrum transformation and wave run-up well. Wave shape transformation is predicted reasonably well by the model, although it is shown that the model does overestimate the wave asymmetry in the lower swash and inner surf zone region. Model results of wave overtopping and local gradients in wave set-up are shown to be accurate if the correct bed level development is imposed on the model, or short sections of the dataset are analysed in which little bed level change takes place. Sensitivity studies showed that modelling of the incident-band wave motion, instead of the wave group motion, was essential in predicting wave overtopping on a gravel barrier, and that groundwater interaction was required to correctly model the evolution of overtopping waves across a gravel barrier.
The results of this paper show that XBeach-G can be applied to estimate the potential storm impact on gravel barriers through a prediction of wave height transformation, wave run-up levels and initial wave overtopping discharge on gravel and composite beaches. However, relevant aspects of the storm response of a gravel barrier, including the development of wave overtopping and wave overwash during a storm, cannot be successfully simulated without morphodynamic updating of the bed level. The inclusion of morphodynamic updating represents the next stage of the XBeach-G model development.
The XBeach-G model, the XBeach-G model source code (Fortran95) and a graphical user interface for the XBeach-G model are available for download on the XBeach project website (www.xbeach.org).