Numerical modelling of intra-wave sediment transport on sandy beaches using a non-hydrostatic, wave-resolving model

The mutual feedback between the swash zone and the surf zone is known to affect large-scale morphodynamic processes such as breaker bar migration on sandy beaches. To fully resolve this feedback in a process-based manner, the morphodynamics in the swash zone and due to swash-swash interactions must be explicitly solved, e.g., by means of a wave-resolving numerical model. Currently, few existing models are able to fully resolve the complex morphodynamics in the swash zone, and none is practically applicable for engineering purposes. This work aims at improving the numerical modelling of the intra-wave sediment transport on sandy beaches in an open-source wave-resolving hydro-morphodynamic framework (e.g., non-hydrostatic XBeach). A transport equation for the intra-wave suspended sediment concentration, including an erosion and a deposition rate, is newly implemented in the model. Two laboratory experiments involving isolated waves and wave trains are simulated to analyse the performance of the model. Numerical results show overall better performance in simulating single waves rather than wave trains. For the latter, the modelling of the morphodynamic response improves in the swash zone compared with the existing sediment transport modelling approach within non-hydrostatic XBeach, while the need of including additional physical processes to better capture sediment transport and bed evolution in the surf zone is highlighted in the paper.


Introduction
Sandy beach evolution plays a key role in coastal vulnerability, influencing the stability of ecosystems and coastal communities' economy and safety. Their morphodynamical evolution and response to drivers such as increased storminess remain difficult to predict (Wong et al. 2014).
The exchange of sediments between the swash zone and the surf zone determines the evolution of the beach et al. 2011, Kim et al. 2017). While hydrostatic NLSWE do not allow frequency dispersion, thus limiting their application, the other two aforementioned sets of flow equations are widely used for wave propagation from intermediate waters to the shoreline. The non-hydrostatic NLSWE are used in the open-source Non-Hydrostatic version of the XBeach (XBNH) model (Smit et al. 2010;Roelvink et al. 2018). To enable the computation of morphodynamics, these models would require sub-models that compute suspended and bed load sediment transport based on intra-wave hydrodynamics, from which in turn bed level change can be computed.
In the framework of non-hydrostatic NLSWE, formulations of sediment transport have been extensively tested for gravel beaches (e.g., McCall et al. 2015), but not for sandy beaches. Ruffini et al. (2020) showed that application of a wave-averaged sediment transport equation (Van Thiel de Vries 2008;Van Rijn 2007) in XBNH led to inaccurate simulated beach morphodynamics, which was related to inaccuracies in modelled sediment concentrations, particularly during flow reversal.
The aim of this work is to improve the modelling of intra-wave sediment dynamics in XBNH. To this end, the Pritchard and Hogg (2003) and Meyer-Peter and Müller (1948) formulations are chosen because of their good performance in the swash zone (see Zhu and Dodd 2015). The present study focuses on modelling the dynamics of sediment in the nearshore zone in all stages of the flow both for solitary waves, i.e., isolated waves, and wave trains, with significant swash-swash interactions. First, the proposed hydro-morphodynamic model is verified against a high-resolution numerical solution of an idealised bore generated by a solitary wave over an erodible sloped beach. Subsequently, two experimental case studies are simulated involving bichromatic wave groups and consecutive noninteracting solitary waves over sandy beaches. For the former, a sensitivity analysis of the results to the parameters used in the Pritchard and Hogg (2003) sediment transport equation is also carried out. This paper is organised as follows: XBNH, with a focus on the sediment transport and bed-updating modelling, is described in Section 2; verification of the proposed model is illustrated in Section 3; the performance of the model against two laboratory experiments is presented in Sections 4 and 5, respectively; discussion of results and concluding remarks with recommendations for future works are given in Sections 6 and 7, respectively.

Model description
XBNH consists of coupled wave-resolving hydrodynamics and morphodynamics equations. The existing approach for suspended and bed load transport in XBNH were originally developed for the Wave-Averaged Sediment Transport (WAST) modelling (see Appendix 1). Therefore, in the present work, a wave-resolving formulation for the suspended sediment transport is implemented in the model. To this end, the transport advection equation of Pritchard and Hogg (2003) is chosen. XBNH computes the bed load transport using the Meyer-Peter and Müller (1948) expression. The combined use of the newly implemented Pritchard and Hogg (2003) equation with the Meyer-Peter and Müller (1948) formula within XBNH for the Intra-Wave Sediment Transport (IWST) modelling is herein referred as XBNH-IWST. In this work, only the cross-shore direction is considered; Fig. 1 shows a schematisation of a typical cross-shore profile with the main variables used.
The model performance is quantified by computing the normalised Root-Mean-Square Error (nRMSE), defined as: where y m,i is the i-th sample of the modelled quantity y, and y ref,i is the i-th sample of the corresponding reference variable (e.g., semi-analytical, experimental); N is the number of samples; s y ref is the standard deviation of the reference quantity, y ref , and it is defined as: Following Bosboom et al. (2020), the Root-Mean-Square Transport Error (RMSTE) is also computed to quantify the model performance in terms of final bed changes. The RMSTE (m 2 ) measures the mismatch between the predicted final bed level, z b f,m , and the reference one, z b f,ref , in terms of the minimum (i.e., optimal) quadratic sediment transport cost required to transform the predictions into the reference field, and it is computed as: where Q i is i-th sample along the cross-shore coordinate, x, of the sediment volume, Q, required to transform z b f,m into z b f,ref . The conservation of mass is satisfied so that ∂Q/∂x = z b f,m − z b f,ref and Q i=1 = 0 is assumed (with i = 1 referring to the onshore boundary of the x− domain, located landward of the maximum run up).
To further assess the correlation between time series of modelled and reference quantities, Pearson's cross-correlation coefficient, −1 ≤ ρ mr ≤1, is used and it is defined as: where cov(y ref , y m ) is the covariance of the time series of modelled and reference quantities and it is computed as: whereȳ m = (1/N) N i y m,i is the mean value of y m , and s y m is the standard deviation of the time series of the predicted quantity.

Hydrodynamics
The hydrodynamics in XBNH is similar to the one-layer version of SWASH (Zijlema et al. 2011). The depthaveraged flow is computed using the non-hydrostatic 1D (one-dimensional) NLSWE: where t is time, η is the water surface elevation from the Still Water Level (SWL), u is the depth-averaged crossshore velocity, h is the total water depth, ν h is the horizontal viscosity, ρ is the density of water, p nh is the depth-averaged dynamic pressure normalised by the density, g is the gravity acceleration constant and τ b is the total bed shear stress, which is computed as: where c f is the dimensionless friction coefficient. Relatively simple unsteady Bottom Boundary Layer (BBL) models, such as the momentum integral method (Sumer et al. 1987) used also in NLSWE solvers (e.g., Briganti et al. 2011), could be considered. However, the results in terms of τ b are comparable with simpler formulations, such as the one considered in this study (see e.g., Briganti et al. 2018). Also, phase differences could be significant and more complex BBL models should be used (e.g., Rijnsdorp et al. 2017). Nevertheless, the detailed modelling of the BBL is outside the scope of the present work. p nh allows to account for wave dispersion with similar accuracy to that of weakly non-linear Boussinesq-type models (Bai et al. 2018). In the cases analysed in this study, the dispersivity parameter, k w d is lower than 0.5, where k w is the wave number (defined as k w = 2π/L, with L being the local wave length) and d is the still water depth (see Fig. 1). Therefore, the celerity error in the description of frequency dispersion is of the order of 1% (Bai et al. 2018).
Wave breaking in XBNH is modelled by using the Hydrostatic Front Approximation (HFA) of Smit et al. (2013), in which the non-hydrostatic pressure term is set to 0 when ∂η ∂t > 0.4c (with c = √ gh the wave celerity in shallow water). After this condition is reached, waves propagate as hydrostatic bores. The reader is referred to Smit et al. (2010) and McCall (2015) for a full description of the XBNH model.

Intra-wave sediment transport modelling
In this study, the Pritchard and Hogg (2003) advection equation for the intra-wave suspended sediment transport is included in XBNH: where C is the depth-averaged suspended sediment concentration; its maximal value, C max , is herein considered the higher physically possible sediment concentration for a fluidised bed and defined as C max = 1−n p,d , with n p,d = 0.6, the porosity for a fluidised bed; S sl represents the bed slope effects computed following Deltares (2015): where α sl = 1.6 according to Deltares (2015) and z b is the bed level. Therefore, the suspended sediment transport rate, q s , is defined as q s = huCS sl . m e is the mobility parameter, which determines the erodibility of the sediment as suspended load, τ b,cr is the critical bed shear stress, τ ref is the reference bed shear stress, R > 0 is a dimensionless exponent (Pritchard and Hogg 2003), w s is the sediment settling velocity and C nb is the near-bed suspended sediment concentration at a small near-bed reference height, d nb , above z b . The two terms on the right side in Eq. (9) represent the erosion rate, E, and the deposition rate, D, respectively. C nb in D is computed as: where the shape factor, K C , represents the relative importance of sediment settling and mixing. When good mixing is assumed K C = 1. Since the suspension is assumed to be sufficiently diluted, there is no feedback between the suspended sediment and the vertical distribution of turbulence. Therefore, K C depends only on sediment properties and the depth-averaged hydrodynamics. Consequently: where B is the Rouse number defined as: where κ = 0.4 is the von Karman constant and u * is the friction velocity: u * = (τ b /ρ) 1/2 . d nb is the dimensionless near-bed reference height, given by a simplified form of Van Rijn formula as shown in Soulsby (1997) with the relationship: in which D 50 is the median grain diameter, and λ is a reference length scale. The model described above allows expressing the vertical distribution of suspended sediment by a power-law profile as in Soulsby (1997): where d z is the vertical elevation from z b (see Fig. 1). The concentration profile described by Eq. (15) corresponds to a linearly increasing eddy diffusivity of the sediment with the height above the bed. Note that C z is not a model output, but it is herein computed in the aftermath of the numerical simulations.
The bed load transport rate, q b , is calculated using the equation derived by Meyer-Peter and Müller (1948). The reader is referred to McCall (2015) for the existing implementation of this formula in XBNH, and it is summarised here, following with the main equation: where θ and θ cr are the Shields and critical Shields parameters, respectively; = (ρ s − ρ)/ρ, in which ρ s is the sediment density. θ is computed as θ = τ b /( ρgD 50 ), and θ cr is given by Soulsby (1997).
The Meyer-Peter and Müller (1948) formula is considered appropriate for the swash zone according to previous studies (see Chardón-Maldonado et al. 2016 among others) and variations of the formula have been tested, for example in Postacchini et al. (2012) for sand and Briganti et al. (2018) for coarse sand. When compared with the original Meyer-Peter and Müller (1948) formula, the Postacchini et al. (2012) formulation showed very similar results in terms of net bed changes (see Briganti et al. 2016). Therefore, in this study we did not test other formulas because of the limited differences in the swash zone shown in the literature.

Bed-updating modelling
The bed-updating is modelled using the Exner-type equation: where n p is the bed porosity; E and D are formulated as shown in Eq. (9), and q b is computed as in Eq. (16).

Numerical scheme
XBNH uses a staggered grid where depth, water level and sediment concentration are defined in the cells centres, and velocity and sediment flux at the cells interfaces. The hydrodynamic equations are solved by applying a limited version of the McCormack (1969) predictor-corrector scheme, which is second-order accurate where the solution is smooth and reduces to first-order accuracy in proximity of discontinuities. The method is mass and momentum conservative (Smit et al. 2010). For the sediment transport and bed-updating modelling, a finite volume approach is applied, where upwind approximations are used (Deltares 2015). XBNH uses a dynamically adjusted time step. Thus, a value for the maximum Courant Number, CN, is defined and the program in turn adjusts the time step, t, in order to guarantee that u t/ x < CN, where x is the computational grid size. In this study, CN = 0.7 was used.

Model verification
In this section, the performance of the hydro-morphodynamic model proposed in this study is verified against the highresolution numerical solution obtained by Zhu and Dodd (2015), in which an idealised bore generated by a solitary wave over an erodible sloped beach was simulated. Figure 2 shows the model domain in Zhu and Dodd (2015). In the region x < −10 m the bed is flat, while for x ≥ −10 m an erodible 1:15 sloped beach is considered. The initial shoreline position is located at x = 5 m. The initial conditions of η and u along the cross-shore direction were given by Mei (1989) and the hydrodynamic Riemann condition, respectively. As shown in Fig. 2, at the initial condition, the wave crest is located at x = −22 m. The wave height, H , is equal to 0.60 m. The governing equations in Zhu and Dodd (2015) were solved using the Method Of Characteristics (MOC), and the hydrodynamics were solved using the hydrostatic NLSWE, which included bed shear stress. The suspended sediment transport was computed using the Pritchard and Hogg (2003) transport equation, assuming a well-mixed condition (i.e., K C = 1). The bed load was given by the Meyer-Peter and Müller (1948) formula.

Model set-up and parameters
The model set-up and physical parameters followed closely those used in Zhu and Dodd (2015). Time series of η and u were provided by Zhu and Dodd (2015) at x 0 = −20 m, i.e., where the wave does not propagate as a bore. Thus, as shown in Fig Initial (Zhu and Dodd, 2015) Fig. 2 Model domain and initial condition in Zhu and Dodd (2015) and upstream boundary location in XBH-IWST model domain (reddashed line) time was approximately 33 s. Table 1 shows a summary of the main parameters included and conditions assumed in the Pritchard and Hogg (2003) equation. Similarly to Zhu and Dodd (2015), the well-mixed condition was assumed (i.e., K C = 1 was set). Consistently with the cited study, bed slope effects in the computation of q s and q b were not taken into account. To compare the response of the model proposed in this study with Zhu and Dodd (2015) in terms of intra-wave sediment transport, the hydrostatic approach was considered by turning off the non-hydrostatic pressure term. This model configuration is herein referred to as XBH-IWST.

Comparison between model predictions and Zhu and Dodd (2015)
Comparison of model predictions with Zhu and Dodd (2015) for the hydrodynamics is shown in Appendix 2; here, only the modelling of beach morphodynamics is discussed. Figure 3 shows the final bed profiles, z b f (a), and the final bed changes, is the time at the end of the simulation. Despite the height of the bed step being underestimated by 20% with respect to Zhu and Dodd (2015), XBH-IWST captures the erosion and deposition well. nRMSE = 0.0085 for z b f and RMSTE = 0.003 m 2 , showing good performances of XBH-IWST in simulating the solution of Zhu and Dodd (2015). Similarly to Briganti et al. (2012), z b f was post-processed by using a moving average. Spurious oscillations are shown in the region x > 2 m during the backwash bore, which runs down the beach and generates a sharp deposition at x 2 m. Note that Zhu and Dodd (2015) used a shock-fitting scheme. Increasing x by an order of magnitude reduces the oscillations; however, as expected, it was found that the much lower resolution would lead to the underestimation of the height of the backwash step by 50%. While not being the aim of this study, the implementation of a shock-capturing numerical scheme could help overcome this issue. Table 1 Main parameters and conditions in the Pritchard and Hogg (2003) transport equation Pritchard and Hogg (2003) equation

Numerical modelling of bichromatic wave groups on an intermediate beach
In this section, the performance of XBNH-IWST is assessed against the experiments conducted within the Hydralab IV  Alsina et al. 2016). The experiments studied the hydromorphodynamics of bichromatic wave groups on a 1:15 sloped beach built at prototype scale with commercial sand characterised by D 50 = 0.25 mm, w s = 0.034 m/s and n p = 0.36, which showed clearly swash-swash interactions. Two bichromatic wave group conditions with the same energy content were generated in the flume: BE1 2 (broadbanded wave condition) and BE4 2 (narrow-banded wave condition), respectively, with varying wave group period, T g , and repeat period, T r . For BE1 2, T g = 15 s and T r = 195 s, whereas for BE4 2, T g = T r = 27.7 s (see also Alsina et al. 2018). T g = 1/f g , where f g is the group frequency, which is defined as the difference of the primary frequencies, f 1 and f 2 . A summary of the bichromatic wave groups is shown in Table 2, where H 1 and H 2 are the wave heights of the primary components.
For each wave condition, starting from the same initial z b (1:15 uniform sloped bed), eight successive bichromatic wave sequences, from SEG1 to SEG8, each of 1800 s duration, were generated. Figure 5 shows the initial z b and the location of the instruments in the wave flume and selected for comparison. Wave Gauges (WG) and Acoustic Wave Gauges (AWG) measured η; Acoustic Doppler Velocimeters (ADV) measured the local flow velocity. Optical Back-Scattering sensors (OBS) and Conductivity Concentration Measurements (CCM + ) tanks measured C z and time-dependent z b in the swash zone, respectively. d in the horizontal part of the domain was 2.48 m for BE1 2 and 2.46 m for BE4 2.
The beach is classified using the dimensionless settling velocity, Ω, defined as Gourlay and Meulen (1968): Ω = H rms /(w s T p ), where H rms is the root-mean-square wave height computed at WG3 (i.e., H rms = 0.39 m for BE1 2 and H rms = 0.40 m for BE4 2) and T p = 1/(f 1 + f 2 )/2) = 3.70 s is the mean primary wave period for both wave conditions. It is computed that 1< Ω <6, which indicates an intermediate beach.
The reader is referred to Alsina et al. (2016) for a detailed description of the experimental procedure.

Model set-up and parameters
The model domain is the same as Ruffini et al. (2020). As shown in Fig. 5, the upstream boundary in the model is located at x 0 = 30.55 m, where the WG3 was installed. Thus, the model domain extended from WG3 to the end of the beach, located at x = 85.05 m, and like the cited study x = 0.1 m. Time series of η and u, updated as offshore forcing, were the same as those used in Ruffini et al. (2020) as boundary conditions. For the computation of c f , a slightly lower value of the Manning coefficient, n, than in Ruffini et al. (2020) was used. n = 0.018 s/m 1/3 was calibrated considering the best compromise between the accuracy of maximum run-up and morphological evolution. The value chosen still reflects the characteristics of the considered sandy beach. Model parameters, which are not mentioned herein, were set to their default values defined in Deltares (2015). The calibration of the sediment transport model in XBNH-IWST was carried out by varying m e , R and λ. Table 3 summarises the main parameters included and conditions assumed in the Pritchard and Hogg (2003) transport equation. This set of parameters was chosen as it provided the best modelling in the sensitivity analysis shown in Section 4.2.
Only the first two segments, SEG1 and SEG2, were simulated for both BE1 2 and BE4 2 because those segments showed larger morphological changes than the subsequent ones. For BE1 2, the experimental bed evolution reached an equilibrium more rapidly compared with BE4 2. However, SEG1 and SEG2 were far from equilibrium for both cases.

Sensitivity analysis for the Pritchard and Hogg (2003) transport equation
The sensitivity analysis of the results to the parameters used in the Pritchard and Hogg (2003) transport equation was carried out for SEG1 of BE1 2 wave condition. The aim of the sensitivity analysis is to show the relative effects of these parameters in terms of the modelled C and z b f = Table 3 Main parameters and conditions in the Pritchard and Hogg (2003) transport equation Pritchard and Hogg (2003) equation The parameters considered are m e , R and λ. Note that λ also affects τ ref (see Table 3) and K C following Eqs. (12) and (14). According to Zhu and Dodd (2015), τ b,cr is not analysed because the effect of a threshold for suspended load is negligible for fine sand; hence, τ b,cr = 0 N/m 2 .
Each parameter is varied by keeping the others to their reference values as in Zhu and Dodd (2015) (i.e., m e = 0.002 m/s, λ = 1 m and R = 1). m e is the least welldetermined parameter, due to the lack of data to provide its estimates. Since m e = 0.002 m/s is found to underestimate both z b f and C, the sensitivity analysis for m e was carried out by increasing it by up to two orders of magnitude with respect to the reference value. R > 0 is a numerical parameter and it was increased and decreased with respect to R = 1 considering R = 0.25, 0.5, 1.5. Values of λ were chosen to be physically representative of the Alsina et al. (2016) configuration. Therefore, λ = H rms = 0.39 m and λ = d = 2.48 m at WG3 were selected. Figure 6 shows the sensitivity analysis for the parameters considered in terms of z b f (Fig. 6a, b and c) and C (Fig. 6d, e and f), respectively. The corresponding nRMSE, ρ mr and RMSTE are presented in Table 4. Note that the experimental C was computed as the average of the observed C z at OBS4 and OBS7, which were located at two different d z above the initial z b at AWG7 (OBS4 at d z = 0.03 m and OBS7 at d z = 0.08 m). The OBSs did not measure when the free surface was lower than the instrument sensor. Therefore, the corresponding nRMSE and ρ mr were computed when at least one of the two OBSs was submerged.
The sensitivity analysis reveals that m e is the most influencing parameter within the ranges considered for both the predicted C and z b f . Variations of λ and R affect C more than z b f in terms of nRMSE. By increasing λ or R, C decreases and the nRMSE increases due to the increasing underestimation of C. For both λ and R, the maximum relative difference in terms of nRMSE is 14%. Instead, for m e , the difference between the maximum and minimum nRMSE is 73%. The sensitivity analysis shows that the variation of the parameters included in the Pritchard and Hogg (2003) model leads to a variability of its peaks and magnitude, as a consequence of the variability of C. However, for all the parameters considered, the low values of ρ mr highlight a poor correlation between the modelled and experimental C. For z b f , variations of λ and R lead to negligible differences in terms of the corresponding nRMSE; differences in terms of RMSTE are lower than 0.05% and 0.8% for λ and R, respectively. z b f is quantitatively more sensitive to the variation of m e , with the difference between the maximum and minimum RMSTE being 25%. From a qualitative point of view, XBNH-IWST is able to capture the peak of the accretion in the upper swash zone if m e is increased by two orders of magnitude with respect to the values suggested in Zhu and Dodd (2015). In particular, the predicted erosion pattern in the upper swash region evolves into a deposition one for m e ≥ 0.05 m/s. Also, by increasing m e , the erosion in the lower swash region increases and the peak of deposition in the surf zone moves shoreward. Different combinations of values for m e , λ and R were selected within the ranges considered, for the model calibration. Results are shown in Fig. 7 and the corresponding nRMSE, ρ mr and RMSTE are presented in Table 5.

Comparison between model predictions and observations
The hydrodynamics response of XBNH-IWST is very similar to that presented in Ruffini et al. (2020). For this reason only the nRMSE and ρ mr for η at selected crossshore locations are shown in Table 6. Figures 8 and 9a show the variation of H rms along the beach profile during SEG2 for BE1 2 and BE4 2, respectively. For both wave conditions, the model is able to capture the evolution of H rms across the domain.
(with t f being the time at the end of SEG2) are illustrated in Figs. 8 and 9 (b and c, respectively); see Table 7 for the corresponding computed nRMSE and RMSTE. Numerical results show a better performance for BE1 2 than BE4 2, in both the surf and swash zones. This is indicated by the lower nRMSE for BE1 2 compared with BE4 2. For the former, XBNH-IWST can capture the deposition in the upper swash zone and the erosion in the lower swash region, whereas the development of the breaker bar is not accurately simulated for both wave conditions. Consequently, H rms is more underestimated in the shoaling zone for BE4 2 than in BE1 2. Indeed, the experimental results suggest that reflection occurred in the shoaling zone due to the bar; thus, H rms increased more than the predicted one.
The net sediment transport rate,q sed , over SEG1 and SEG2 is shown in Figs. 8 and 9d for BE1 2 and BE4 2, respectively. This was computed using a sediment balance, which was numerically integrated over the x-domain between the start of SEG1 and the end of SEG2: where q sed is the instantaneous sediment transport and the bar refers to the averaging over the duration of the two segments, t SEG1−2 ; the subscript i refers to the ith point along the x-domain for both the numerical mesh and the experimental domain, where z b is available. Thus, i = 1, ...N, with i = 1 at the onshore boundary of the domain (i.e., landward of the maximum run-up limit), where q sed = 0 is assumed, and i = N at the offshore start of the beach. z b SEG1−2 is the difference between z b at the end of SEG2 and at the start of SEG1. Figure 8 (d) highlights that XBNH-IWST is able to simulate the magnitude of the onshore-directed sediment transport in the upper swash zone and the offshore-directed one in the lower swash region and surf zone up to the crest of the bar, located at x = 65 m (see also Appendix 3). For BE4 2, the model can capture the sign ofq sed up to the bar at x = 63 m (Fig. 9d), but the magnitude is underestimated. This might be explained by the more prominent bar observed in BE4 2 than BE1 2, which XBNH-IWST cannot reproduce. Therefore, the exchange of sediment between the swash and surf zones is not well simulated, resulting in a deterioration  in the overall modelling ofq sed . For both wave conditions, some limitations are visible in the shoaling region and surf zone up to the bar crest, where the experimental onshoredirectedq sed is not predicted by the model. Indeed, when the experimentalq sed changes in sign, the modelled one continues being negative for both wave conditions. Note that for BE4 2 the observedq sed goes to zero at the offshore boundary, which is not shown in Fig. 9d. However, the positive and quasi-uniform value of the observedq sed in the shoaling zone is most likely affected by measurement effects due to the mechanical wheel profiler used to measure the bed level. This instrument has a wheel that is too large to detect individual ripples. Therefore, the change in the bed level is below the sensitivity of the instrument. Moreover, the modelled θ at the offshore side of the bar is larger than θ cr for the most part of the event.
A detailed analysis of the local sediment transport dynamics within SEG2 at AWG7 (x = 75.81 m) is shown in Figs. 10 and 11 for BE1 2 and BE4 2, respectively. Note that the observed u is a local measurement and herein assumed depth-uniform. The experimental timedependent z b recorded by the CCM+ tank at the same x was used to compute C z for XBNH-IWST with Eq. (15). For BE1 2, only two groups are selected from the sequence of groups within T r over SEG2. Figure 10b and c show the time series of C z for both the model and observations at OBS4 and OBS7, respectively. Results show that the initial peak of the intra-wave C z corresponding to the first bore of both groups (at t/T r 0.005 and t/T r 0.084, respectively) is captured by XBNH-IWST, as well as the order of magnitude of C z corresponding to the peaks of u. However, the correlation between the predicted and observed C z is affected by the underestimation of C z close to flow reversal, as confirmed by the corresponding nRMSE and ρ mr , shown in Table 7, which reflect a lower performance of XBNH-IWST compared with the hydrodynamics modelling. Figure 10d shows the modelled q s and q b ; q s being always higher than q b . In the broadbanded wave condition, large backwashes are allowed to develop. Consequently, the suspension of sediment particles is dominant with respect to the sediments settling. BE4 2 allows analysis of results within one group over SEG2, since T r = T g . Note that the observed u was filtered with a low-pass filter (cut-off frequency set to 3 Hz) to remove the noise in the measurements. Similarly to BE1 2, XBNH-IWST is able to capture the order of magnitude of the observed C z after the first event due to swash-swash interactions present in the group at OBS4, whereas, except for the first uprush within the group, C z is underestimated at OBS7 (Fig. 11b and c, respectively). The peak in the modelled C z corresponding to the first wave (at t/T r 0.095) might be the result of a larger bore-induced advection than observations. The lower model performance for BE4 4 than BE1 2 is reflected by the higher nRMSE for C z , while the values of ρ mr are of the same order of magnitude of those for the broad-banded wave condition (see Table 7). In BE4 2, a higher number of swash-swash interactions occurred within the group than in BE1 2. Therefore, backwashes corresponding to subsequent events were allowed to develop for a shorter duration compared with the broad-banded wave condition. Consequently, the observed and predicted C z are lower and the difference between the modelled q s and q b (Fig. 11d) is smaller than in BE1 2. Indeed, q b > q s during the backwash events within the group.

Numerical modelling of consecutive non-interacting solitary waves over a sloped beach
In this section, the performance of XBNH-IWST is further tested against the experiments of an erodible sloped beach exposed to nine consecutive non-interacting solitary waves presented in Young et al. (2010). Similarly to Zhu and Dodd (2015), this case allows an individual swash event and the evolution of the beach to be examined without the presence of swash-swash interactions occurring in wave groups. Figure 12 shows the experimental set-up and the location of the instrumentation installed in Young et al. (2010) and considered for comparison in this study. For η, WG8 (x = 23 m) and the Distance Sonic, DS2 (x = 29 m) sensors were considered. For u, ADV8 (x = 23 m) and ADV5 (x =

Model set-up and parameters
The model domain was set up following the experimental settings, described in detail in Young et al. (2010). The computational domain was x 0 = 0 ≤ x ≤ 40 m and x = 0.05 m. The initial wave-modified z b was used as the initial z b in the simulations. As shown in Fig. 12, time series of η and u were updated at the upstream boundary located at x 0 , and were given by Titov and Synolakis (1995), considering a solitary wave over the initial d = 1 m: where H = 0.60 m and u is:  Fig. 11 Time series of flow and sediment transport variables at AWG7 over SEG2 for BE4 2: u a; C z at OBS4 b; C z at OBS7 c; q s and q b d; reference line: grey dashed line Moreover, to be consistent with the observations, a reflection was taken into account at the upstream boundary. The sediment transport model was set up similarly to the Alsina et al. (2016) test case (see Table 3), with λ = H . c f was modelled using n = 0.025 s/m 1/3 , which was chosen matching the simulated maximum run up with the observations. The maximum excursion point was observed at x = 38.5 m. Parameters that were not mentioned in this study were set to their default values defined in Deltares (2015). Numerical simulations were performed for the first three waves of the nine experimental runs. Following the   Figure 13 shows the numerical and experimental time series of η and u within one wave at selected x. In general, XBNH-IWST is able to capture the hydrodynamics. However, a shift of 4 s in the prediction of the reflected wave due to the finite size of the flume is also visible, which affects the corresponding nRMSE and ρ mr (see Table 8). At WG8 (Fig. 13a), the average overestimation of η after the wave run-down, including the reflected wave (i.e., for t > 22 s), is equal to 8%. Comparison at ADV5 (Fig. 13d) is affected by some noise in the collected signal at t = 10 s and when the water level dropped down the sensor (17 < t < 37 s); hence, no signal was recorded. Figure 14 shows z b f (a) and z b f (b) after 3 waves (z b f = z b (x, t f ), with t f = time at the end of the third run. Despite the overestimation of the deposition in the upper swash zone and the erosion in the lower swash region for x > 30.5 m, XBNH-IWST is able to reproduce the observed   Table 8. An analysis of the intra-wave sediment transport within one wave is shown in Fig. 15 for x = 23 m. Time series of C are shown in Fig. 15b. Note that the experimental C was computed as the average of the observed C z at OBS3 and OBS4. The corresponding computed nRMSE and ρ mr are shown in Table 8. Numerical results show that the largest peak of C corresponds to the backwash phase of the flow (Fig. 15a) (at t 20 s). According to Young et al. (2010), the large peak of C observed between t = 20 s and t = 25 s was possibly due to the dispersion of sediment generated by the hydraulic jump observed during the experiments at x = 24 m.

Comparison between model predictions and observations
Time series of the modelled q s and q b are shown in Fig. 15c. At the early stage of the uprush q b q s . This is consistent with the observations, where the experimental C is almost zero in the uprush, which means that sediment motion mainly occurred as near-bed sediment transport. The reason might be addressed by the location of the OBS sensors, which was seaward of the wave plunging point. As the stirred up sediments are entrained in the water column, the contribution of the modelled q s also increases. Close to flow reversal, sediment settling occurs and both predicted q s and q b decrease. During the backwash, q b and q s increase until q s > q b in the stage of the run-down.

Discussion
First, the XBNH-IWST modelling of sediment transport and morphodynamics is discussed by comparing the results of numerical simulations of the Alsina et al. (2016) case obtained in the present study with those of Ruffini et al. (2020) using the XBNH-WAST approach. Figure 16 shows the XBNH-IWST and XBNH-WAST predictions for SEG2, with a focus on the swash zone for BE1 2 and BE4 2, respectively. Comparison between the modelled and experimental C is shown in Fig. 16d and h and the corresponding nRMSE, ρ mr and RMSTE are presented in Table 9. The accuracy of XBNH-IWST in the prediction of C in terms of nRMSE and ρ mr is similar to that of XBNH-WAST (see Table 9). XBNH-ISWT, however, is able to better describe the sediment suspension observed after the first bore generated by swash-swash interactions within the group and C close to flow reversal. Also, the accuracy of XBNH-IWST in terms of nRMSE and ρ mr is higher when the comparison is carried out for C z (see also Table 7).
Differences in the predictions of C for the two approaches lead, in turn, to differences in the simulated z b f (Fig. 16a and e). XBNH-IWST shows a better performance in the prediction of z b f than XBNH-WAST; for BE1 2, the RMSTE of the former approach is lower by 21% than that of the latter one, while for BE4 2 the difference is 10%. XBNH-IWST better simulates the deposition in the upper swash zone and the erosion in the lower swash region. In fact, z b f predicted with XBNH-WAST diverges from observations, especially in the upper swash zone. This might be explained by the behaviour of the sediment transport model in XBNH-IWST near flow reversal (i.e., u = 0 m/s), when sediment particles are allowed to settle. The better performance of the present approach is confirmed by the modelling ofq sed (Fig. 16b and f). Unlike the model proposed in this study, XBNH-WAST does not capture the sign ofq sed in the upper swash zone. For the Young et al. (2010) experiments, the predicted z b f shows lower nRMSE and RMSTE than for the Alsina et al. (2016) case (see Table 8) by 60%, and by an order of magnitude respectively. Regarding the intrawave sediment transport, the experimental C was obtained by averaging the measurements performed by OBS3 and OBS4 in the surf zone, where, as indicated in Young et al. (2010), C z was observed to vary a lot along the depth. Therefore, discrepancies between observations and numerical predictions are due to some limitations of the model performance, but uncertainty in the comparison exists because of the low resolution of the measurements.
For both the simulated laboratory experiments, the performance of the Meyer-Peter and Müller (1948) formulation is consistent with results of previous studies (e.g., Xiao et al. 2010 andPostacchini et al. 2012). The modelled θ for Young et al. (2010) at x = 23 m is higher than θ cr for almost the whole duration of both the uprush and backwash events (with maximum values of approximately 15 and 4, respectively). showed that m e is the parameter that mostly influences the prediction of z b f and C (see Table 4 and Fig. 6). However, because no field data are available to provide its estimates, the only comparisons with the findings of the present study are previous numerical studies, such that of Zhu and Dodd (2015). The value of m e that showed best modelling in this study is an order of magnitude larger than that used by the cited study for a single solitary wave when a well-mixed condition is assumed (i.e., K C = 1). The optimal m e found was used in combination with R = 1.5 (R being an arbitrary numerical parameter), λ equal to the considered wave height at the upstream boundary, and by considering K C ≥ 1. To simplify the choice of the values corresponding to the parameters considered in the Pritchard and Hogg (2003) expression, τ b,cr was set equal to 0. This choice was justified following Zhu and Dodd (2015) who pointed out that the effect of a threshold for suspended load is not significant for sandy beach morphodynamics. λ is an arbitrary length scale, which in turn, affects the other scale parameter, τ ref . Therefore, representative values of λ were selected for the model calibration (see Section 4.2). τ ref was computed similarly to Zhu and Dodd (2015) for all the simulations carried out.
Although the proposed parameter set allows obtaining a lower RMSTE and a better prediction of erosion and deposition patterns than the other combinations tested, the prediction of C is still inaccurate in XBNH-ISWT as well as the process of bar formation, indicating the need of explicitly modelling further important physical processes. The missing explicit representations of processes such as wave breaking-induced turbulence in the HFA model, phasing effects of the velocity in the BBL, such as the earlier flow reversal near the bottom with respect to the rest of the water column in the swash zone (see e.g., Zhang and Liu 2008), are certainly factors that concur to this inaccuracy.

Conclusions
In this study, the implementation of the Pritchard and Hogg (2003) intra-wave sediment transport equation in XBNH was presented, together with the verification of the model against a high-resolution numerical solution. A sensitivity analysis allowed choosing the set of parameters included in the Pritchard and Hogg (2003) equation that provided the best modelling and the performance of the resulting model was assessed with two laboratory test cases. Verification of XBH-IWST against Zhu and Dodd (2015) highlighted that the Pritchard and Hogg (2003) transport equation performs qualitatively and quantitatively well when compared with a high-resolution numerical solution of NLSWE. Therefore, this modelling approach is suitable for solving the intra-swash sediment transport in the context of wave-resolving models.
Numerical simulations of the Alsina et al. (2016) and Young et al. (2010) laboratory experiments showed that XBNH-IWST is able to simulate the beach evolution driven by isolated waves, whereas the exchange of sediments between the swash and surf zones under wave trains is not accurately predicted yet. The sediment transport modelling approach herein used is able to improve the simulation of the intra-swash sediment transport within XBNH, especially in terms of C at flow reversal, and in turn, the prediction of the beach accretion in the upper swash zone for the Alsina et al. (2016) case. However, it is not sufficient to accurately capture the mutual feedback between the swash and surf zones when swash-swash interactions are present.
In evaluating the comparisons with the laboratory data, consideration was given to the associated uncertainty. In the Alsina et al. (2016) experiments, a state of the art measurement system was used with high resolution in space and time for the suspended sediment transport. The OBS sensors inevitably provide low-resolution data in the water column. On the other hand, the previous Young et al. (2010) experiments were conducted with fewer instruments. This provides limited resolution for the data. High-resolution experiments of solitary waves on mobile bed are particularly valuable because the intra-wave dynamics could be investigated in depth. Therefore, such experiments are highly recommended for future studies.
At present, there are still limitations in the qualitative and quantitative representation of the morphodynamics and sediment transport with XBNH-IWST. The inclusion of the wave-resolving Pritchard and Hogg (2003) sediment transport equation in XBNH represents a first step to improve the morphodynamics prediction of the model in the context of sandy beaches evolution. Future work is necessary for an accurate prediction of the mutual feedback between the surf and swash zones. Further work will be addressed to the improvement of the intra-wave sediment transport modelling within XBNH-IWST. Additional physical processes, such as BBL effects and turbulence, and the inclusion of the vertical structure of sediment concentration and flow, need to be included in a formulation suitable for the class of models at hand. A simplified wave breakinginduced turbulence model, similar to those used by Alsina et al. (2009) and Reniers et al. (2013), could be considered to take into account the turbulence diffusion in the computation of C(z) and, in turn, in the prediction of C with XBNH-IWST.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.

Appendix 1: Sediment transport modelling in XBNH-WAST
In XBNH-WAST, the advection equation for the sediment transport is formulated as: where T s is the adaptation time, which represents the entrainment of the sediment, depending on h and w s , and C eq is the total equilibrium sediment concentration. XBNH-WAST calculates the equilibrium concentration using the the Van Rijn et al. (2007), Van Thiel de Vries (2008) formulations. For both the bed load and the suspended load, referred to as b and s, respectively, the equilibrium concentrations are expressed in function ("φ") of u and the wave orbital velocity, u rms , as: C eq,b,s = A b,s φ (u, u rms ) , where A b,s represents, alternatively, the bed load and suspended load coefficients, depending on sediments' grain size and flow properties. The Van Rijn et al. (2007) and Van Thiel de Vries (2008) sediment transport formulations were originally developed for the wave-averaged version of XBeach, where u represents the wave-averaged crossshore flow velocity, computed by means of NLSWE, and u rms is computed as a parameterisation of the short-wave energy (Roelvink et al. 2009). XBNH solves the intra-wave flow through the extended NLSWE, and the total intra-wave cross-shore velocity is included in the term "u." Consequently, in XBNH-WAST, the term "u rms " in Eq. (22) is equal to 0. For a more detailed implementation of XBNH-WAST, the reader is referred to Deltares (2015). (2015) for the hydrodynamics Figure 17 shows the time series of h and the scaled velocity, u/ √ gλ, at two different cross-shore locations. Overall, results show a very good agreement in terms of

Appendix 3: Representation of the net return flow in XBNH-IWST
Since XBNH is based on the one-layer version of the SWASH model, no vertical discretisation of the velocity profile is available. Therefore, only a depth-and phaseaveraged net current can be obtained from the model. Results show that the predicted net current is directed offshore. Figure 18 shows the modelled and experimental u for SEG1 of BE1 2 for the Alsina et al. (2016) case and the wave-averaged velocity over the wave group period, T g = 15 s, u mean ; u mean = − 0.46 m/s for XBNH-IWST and u mean = − 0.48 m/s for the experiments. Figure 18 also shows the modelled q tot = q s + q b for the same segment. q tot is also depth-averaged, and the presence of a net depth-averaged velocity leads to net mean offshore transport, q tot,mean = -0.001 m 3 /m/s, even if entrainment is considered.

Fig. 18
Time series of u and q tot , and u mean and q tot,mean at x = 75.81 m over T g for SEG1 of BE1 2 (Alsina et al. 2016)