Estimation of Irregular Wave Runup on Intermediate and Reflective Beaches Using a Phase-Resolving Numerical Model

Accurate wave runup estimations are of great interest for coastal risk assessment and engineering design. Phase-resolving depth-integrated numerical models offer a promising alternative to commonly used empirical formulae at relatively low computational cost. Several operational models are currently freely available and have been extensively used in recent years for the computation of nearshore wave transformations and runup. However, recommendations for best practices on how to correctly utilize these models in computations of runup processes are still sparse. In this work, the Boussinesq-type model BOSZ is applied to calculate runup from irregular waves on intermediate and reflective beaches. The results are compared to an extensive laboratory data set of LiDAR measurements from wave transformation and shoreline elevation oscillations. The physical processes within the surf and swash zones such as the transfer from gravity to infragravity energy and dissipation are accurately accounted for. In addition, time series of the shoreline oscillations are well captured by the model. Comparisons of statistical values such as R2% show relative errors of less than 6%. The sensitivity of the results to various model parameters is investigated to allow for recommendations of best practices for modeling runup with phase-resolving depth-integrated models. While the breaking index is not found to be a key parameter for the examined cases, the grid size and the threshold depth, at which the runup is computed, are found to have significant influence on the results. The use of a time series, which includes both amplitude and phase information, is required for an accurate modeling of swash processes, as shown by computations with different sets of random waves, displaying a high variability and decreasing the agreement between the experiment and the model results substantially. The infragravity swash SIG is found to be sensitive to the initial phase distribution, likely because it is related to the short wave envelope.


Introduction
Estimation of the total water level (TWL) at the shoreline is an important asset for coastal engineers and those involved in coastal zone management and engineering design. For instance, the TWL describes one of the key components in forecasting tools for the assessment of coastal flood risk or storm impact intensity [1][2][3]. The empirical formulae commonly used to design coastal structures, such as sea walls or rubble mound breakwaters, also rely on the determination of the maximum water level [4]. The calculation of the TWL has thus been the subject of many studies [5,6] aiming in particular to improve the estimation of the wave-induced runup [7][8][9][10][11][12], which is one of the primary contributions to the TWL with tide and atmospheric surge.
Wave runup is composed of a mean time component, the wave setup, and a time-varying component, the swash [13]. The setup depends on an increase in mean sea level at the wave period scale that balances the onshore component of the momentum flux of the waves in the breaking and surf zones [14]. The swash is composed of a short-wave (SW) or incident wave component, corresponding to high frequency oscillations of the water level in the frequency band between 0.04 and 0.25 Hz, and an infragravity (IG) component corresponding to the contribution of long waves with frequency ranging between 0.002 and 0.04 Hz. Therefore, the accurate determination of the runup contribution to the TWL requires tackling a series of challenges associated with the processes of transformation of short waves from intermediate to shallow waters together with the interaction of bound and free long waves. In addition, the respective contribution of SW and IG waves will depend on the type of beach. In the case of a dissipative beach, the dynamics of the swash zone will be dominated by IG waves, whereas for intermediate to reflective beaches both types of waves will contribute to the TWL at the shoreline [15].
One type of approach to estimate the runup consists of applying empirical formulations derived either from laboratory data [16][17][18][19][20] or field observations [13,[21][22][23][24]. These formulae have the advantage of providing an estimation of the runup essentially based on the knowledge of offshore bulk wave parameters, such as the significant wave height H s the peak wave period T p , and the beach geometry as the foreshore beach slope β f . This type of approach can easily be implemented into coastal risk forecasting tools based on fast and low cost computations. However, their application to beaches with complex 3D features is usually limited [7,10,23]. Indeed, comparing several empirical formulations Atkinson et al. [23] showed that the most accurate models give a relative error of R 2% of up to 25%. Thus, it is often necessary to develop site-specific runup formulations [25,26], which require a significant measurement effort to cover a wide range of oceanographic conditions at a given site [25]. Furthermore, it is often hazardous to collect data in natural environments, especially during extremely energetic events [10,27], which is probably the reason for the sparse existence of runup data from extreme events. Another limitation of empirical formulae is that they do not provide any information on the physical processes that control the wave-induced water level at the shoreline.
The limits of empirical formulae can be overcome through application of process-based deterministic numerical wave models. For instance, a phase-and-depth resolving model based on the Reynolds-averaged Navier-Stokes equations was recently used to study the sources of runup variability on planar beaches [28]. However, the application of this type of wave model is mainly intended for in-depth studies of physical processes that control wave transformations and their interactions with coastal structures [29]. Indeed, the high resolution and long computation time limit their application to real beach configurations. Phase-resolving and depth-integrated models offer a promising alternative. This type of approach allows to account for the main processes of wave transformation in intermediate and shallow waters, including dispersion and nonlinear effects, while requiring an acceptable computation time. For instance, the SWASH model [30], a widely used nonhydrostatic nonlinear shallow water model, was applied in 1D mode on an urbanized field site [11] and in 2D mode on a natural open sandy beach [10] to compute storm-induced runup. The COULWAVE model [31], a weakly dispersive and fully nonlinear Boussinesq-type model, was used to investigate wave processes in a fringing coral reef environment at two atoll sites in the western tropical Pacific [32]. The BOSZ model [33], a weakly dispersive and weakly nonlinear Boussinesq-type model was used to compute wave setup induced by energetic breaking waves at a fringing reef site in Hawai'i [34]. The model was incorporated into a full model suite for coastal inundation [35] and later used for probabilistic mapping of storm-induced coastal inundation under climate change scenarios [36]. Both studies involved large computational domains with millions of cells. Most of the previously cited studies demonstrate the ability of phase-resolving depth-averaged models to compute the cross-shore sea and swell waves, the IG waves, and wave-induced setup. This type of approach also succeeds in correctly estimating the 2% exceedance runup value (R 2% ) , which is usually used as an indicator of storm impact intensity. However, few studies have focused on the detailed computation of time-varying swash dynamics.
Accurate measurements of water level oscillation at the shoreline under real conditions are usually difficult to perform. It requires one to instrument a thin layer of water, usually during energetic wave conditions in a changing environment. Laboratory data can offer the advantage of providing synchronized high temporal and spatial resolution of wave transformations and wave-driven water level oscillations under controlled conditions. Free surface elevations can be measured using resistance or capacitance wave gauges distributed along a cross-shore transect. The runup oscillation on the beach face is usually measured using a long capacitance wire gauge mounted normal to the beach slope at a fixed height above the bottom. For example, the data collected during the GLOBEX project [37] was used to validate the application of SWASH to compute the runup variability under dissipative conditions corresponding to irregular waves breaking over a gentle slope [38] for three different incident wave conditions. The relative runup errors ranged between 1% and 11%. Laboratory data of free surface elevation for eight gauges and runup oscillations were used to provide a comprehensive and detailed methodology for sensitivity analysis, calibration, and validation of the SWASH model for its application to the computation of runup oscillations over fringing reefs [39]. The fully nonlinear and dispersive Boussinesq-type model FUNWAVE-TVD was tested with laboratory data to assess its ability to predict the cross-shore evolution of significant wave heights in the SW and IG frequency bands, and the runup spectrum for irregular waves propagating over a laboratory scale fringing reef [40]. Detailed data from a laboratory experiment for waves breaking over submerged reef [41] were also used to validate the computation of runup over a steep-sided coastal structure with a Boussinesq-type model [42]. As an alternative to commonly used measuring devices such as pressure sensors or runup wires, the use of LiDAR scanners in coastal research is becoming increasingly common in both field [43][44][45] and laboratory [46] experiments. The use of a LiDAR scanner provides a continuous description of the area of interest, as opposed to point-by-point measurements with previously cited measuring devices. A single instrument is required to cover a relatively large area. Moreover, LiDAR scanners allow for remote measurements, thus providing data in a nonintrusive way. This can be critical for studies of small-scale processes where surface piercing instruments can lead to obstruction or disturbance of the flow.
In the present work, the Boussinesq-type wave model BOSZ [33,34] is compared to an extensive runup laboratory data set based on LiDAR data obtained during the DynaRev large-scale experiment [46]. The study focuses on the ability of a depth-integrated phase-resolving model to accurately compute both the setup and the contribution of SW and IG waves of the time-varying water elevation at the shoreline in intermediate and reflective beach configurations. Furthermore, this work presents a detailed sensitivity analysis of the computed results to the model settings, including the influence of phase distribution of the incident short waves and the definition of the threshold value for runup determination.
The outline of the paper is as follows. In Section 2, the model used in this study is described together with the laboratory data set. The comparisons between the model computations and measurements of spectral wave characteristics, free surface elevation time series, and runup components are presented in Section 3. A discussion of best practices for proper setup of a phase-resolving wave model with the objective to compute runup oscillations over intermediate and reflective beaches is presented in Section 4. Finally, concluding remarks are drawn in Section 5.

Laboratory Experiment
The experimental runup data used in this study were collected during the DynaRev physical experiment that was carried out in the Großer Wellenkanal, GWK, (Large Wave Flume) in Hanover, Germany. The experiment originally focused on the investigation of the performance of specific revetments against erosion and runup excursions under varying wave conditions and water levels.
A series of tests was performed in the 309 m long, 7 m deep, and 5 m wide canal equipped with a combined piston-flap type wavemaker. A more comprehensive description of the laboratory experiments is presented in Blenkinsopp et al. [46].
In total, more than 130 runs were completed for a total of 141.6 h of testing. The wave conditions followed the distribution of a JONSWAP spectrum. The tested wave conditions varied from H s = 0.6 m to H s = 1.2 m and from T p = 6 s to T p = 12 s with the water level between 4.5 m and 4.9 m. In the present study, only the 20 min long runs are considered. The sandy bottom was initially composed of a 1/15 planar slope (≈6.6%) (Figure 1). After the first runs, the bed reached an equilibrium with the development of a stable inner bar. Measurements of the bottom profile before and after each 20 min run show no significant bed changes occurred during this time.  Table 1) with WG1 and WG2 wave gauge positions (red circles) as well as LiDAR scanner zone (grey area). WG3 and WG4 (red square) refer to virtual wave gauges derived from the LiDAR scanner measurements. The dashed line shows the profile for case SBE2 (see Table 1) indicating the inner bar.
Outside the surfzone, the free surface displacement was measured with two wave gauges WG1 and WG2 ( Figure 1) located at x = 50 m and x = 170 m, respectively. The free surface displacement in the entire surf and swash zones was measured at high resolution using an array of three SICK LMS511 2D LiDAR scanners mounted on the experimental roof. Those scanners allowed to measure the free surface elevation continuously over an extent of 65 m (grey area on Figure 1), starting from the cross-shore position of x = 215 m up to the upper part of the beach at x = 280 m with a resolution of 10 cm. For detailed comparisons between model results and measurements in the surf and swash zone, two virtual gauges WG3 and WG4 ( Figure 1) were set in order to extract water elevation time series from the LiDAR data. WG3 was located at a depth of 1 m and WG4 was positioned deep into the surfzone at a location where the profile is always submerged by at least 10 cm. For case DR0 (see Table 1), cross-shore locations of WG3 and WG4 correspond to x = 235 m and x = 245 m, respectively. Due to the nature of the LiDAR scanner vertical splashes occurring during the wave overturning process can potentially lead to abnormally high nonphysical values. After careful assessment of the data, these recording were removed from the samples. Wave gauges and LiDAR scanners were sampled at 25 Hz. Each run lasted 20 min and a spin-up time of 1 min was considered to allow the wave field to be fully saturated.

Selected Test Cases
Three test cases were selected from the DynaRev data set as they are representative of different beach states. Usually, the beach state is defined by the Irribaren number [17] or surf similarity parameter given by: where β f is the foreshore slope and H 0 and L 0 are the deep-water significant wave height and wavelength, respectively. In the present study, the foreshore slope is calculated as the average slope over the swash zone as defined in Stockdon et al. [22]. Traditionally, a low Iribarren number (ξ < 0.3) indicates a dissipative beach. For higher values, ranging between 0.3 and 1.25, the beach is classified as intermediate, and for ξ > 1.25, it is classified as reflective [47]. According to this classification, the three selected cases allow to study runup for intermediate and reflective conditions. They are summarized in Table 1. The Ursell number, U r = Hλ 2 /h 3 , where H is the total wave height, λ the wavelength, and h the still water depth, expresses the degree of nonlinearity of the waves.

Numerical Model
For the present study, the Boussinesq-type wave model BOSZ [33,34] is used to compute the runup for irregular waves propagating over intermediate and reflective beaches. This phase-resolving depth-integrated numerical model is based on a conserved variable formulation of Nwogu's equations [48]. Contrary to the nonlinear shallow water equations (NLSWE), Boussinesq equations naturally include dispersion terms to account for the nonhydrostatic pressure effects of periodic waves. To account for frequency dispersion, Nwogu [48] expressed the vertical gradient of the horizontal velocity at an arbitrary depth z α through a truncated Taylor series expansion in combination with the irrotationality condition u z = w x . This allows for an approximation of the horizontal velocity's vertical variation in terms of only the horizontal velocity components at z α . The third momentum equation in the z-direction vanishes from depth integration and a pseudo-3D solution is obtained in the 2D horizontal plane. The resulting set of equations agrees well with the Airy theory in terms of its dispersion properties for kh < π or a little beyond that. For kh > π, the dispersion error increases, which causes an overestimation of the wavelength in deep water typical for most equations of this type. The value of z α can be adjusted for an optimal compromise between linear dispersion and shoaling properties. For most applications, z α = −0.531 · h works reasonably well and is used throughout this study. The set of equations, expressed in conservative form, consists of a continuity equation and two momentum equations in the xand y-directions. In 1D, only the terms in the x-direction are retained: The second and third term in the momentum equation arise from the Boussinesq-type approximation derived by Nwogu [48] and account for nonhydrostatic pressure correction. The term uψ C is not part of the original equation by Nwogu [48] and due to conserved variable formulation of Roeber et al. [33]. The term results from the following expression of the momentum equation of the NLSWE with the local acceleration term expressed with the conserved variable (Hu) instead of with only u: as H t + (Hu) x = −ψ C from the continuity equation. The momentum equation therefore includes information from the continuity equation that supports the correct representation of a shock front in terms of flow depth, speed, and dissipation. The mass source term, ψ wm , serves as the internal wavemaker for the generation of periodic waves (see Wei et al. [49]). Here, h denotes the water depth, η denotes the free surface elevation, and H = h + η denotes the total flow depth. The variables are defined on Figure 2. The subscripts x and t stand for partial derivatives with respect to space and time, g is the gravitational acceleration, u is the horizontal flow velocity at the reference depth z α , and τ 1 is the Manning roughness term. Since the equations are depth-integrated, the solution cannot describe wave overturning, which essentially requires more than one value of the solution in the vertical direction. Therefore, a breaking wave is approximated as bore or hydraulic jump with a discontinuous profile. Since the governing equations contain the elliptic dispersion terms, discontinuous solutions are not directly possible, but they require special treatment along the breaking wave front. One possibility is the deactivation of the dispersion terms locally and momentarily based on particular conditions queried in each time step. The BOSZ code offers several options to identify individual cells along the wavefront where the dispersion terms can be ignored and the solutions from the subset of the NLSWE can be used instead. These options include criteria based on geometry and kinematics. Here, the criterion based on a free surface elevation to water depth ratio is used, expressed as follows: If the ratio of free surface to water depth is larger than a given value C b , the dispersion terms are deactivated. It is important to note that no additional term is strictly required to account for the energy dissipation under breaking waves, since the shock-solution of the underlying NLSWE in combination with a shock-capturing numerical scheme properly accounts for the dissipation rate. The approach of deactivation of dispersion terms has been used in multiple previous studies and provides a robust solution as long as the grid spacing is not excessively fine.
Theoretical calculations from McCowan [50] showed that the highest ratio of η/h in shallow water prior to breaking was 0.78, which is the value adopted here. However, the study by McCowan [50] refers to the ratio of wave height to water depth that includes the leading trough, which is difficult to compute on the fly. The sensitivity of the results to the value of C b is addressed in the discussion.

Model Settings
The upstream boundary of the numerical domain is lined up with the first gauge WG1 of the physical experiment. The total length of the numerical domain is 222 m long. The grid size of the computational domain is constant and set to dx = 0.5 m. For each test case, the free surface elevation time series measured at WG1 is prescribed at the offshore wave boundary, allowing the model to be forced with the exact incident wave phases. Since the incident wave velocity imposed at the offshore boundary is not available from the experiment measurements, it is computed using the shallow water approximation based on Airy wave theory that relates wave speed to the local water depth: where U is the horizontal particle velocity. Though the time series clearly shows dispersive waves, the selected conditions are of relatively low frequency dispersion (kh < 0.25π), which reasonably justifies the long wave approximation. In order to match the conditions of the experiment, a constant Manning friction coefficient of n = 0.02 sm −1/3 is used to account for sand of medium grain size. The time step is adaptive to ensure stability under the Courant-Friedriechs-Lewy (CFL) condition of a fixed value of 0.5. The free surface was saved to output files at a frequency of 10 Hz.

Data Analysis
The power spectral density (PSD) of the free surface elevation is computed by applying a fast Fourier transform (FFT) algorithm to 5 min segments of the 19 min time series using a 50% overlapping Hanning window, resulting in about 20 degrees of freedom and a frequency resolution of d f = 0.0033 Hz. Associated 95% confidence intervals are calculated according to Emery and Thomson [51]. While 19 min is relatively short for low frequency analysis, it is a common duration in runup studies [7,10,22,52], as the results over this time span are hardly influenced by tidal fluctuations or offshore conditions. Nonetheless, given the short duration of the experimental data, the results of low frequency components, i.e., periods longer than 1 min, are subject to uncertainty. The significant wave height is then calculated as: where m 0 denotes the zero-order spectrum moment. Instantaneous free surface elevation is averaged over the 19 min at each cross-shore location to provide spatial distribution of the setupη. The time-varying shoreline elevation η s (t) is tracked as the last cell in the shoreward direction where the water depth is greater than a threshold depth δ. The value of δ was fixed to 10 cm as recommended in previous studies [7]. This value was applied to process both numerical and physical data. An example of the results of shoreline tracking carried on scanner data is shown in Figure 3.
In the numerical experiment, all calculations are performed using fixed bathymetry. Comparisons of the bathymetric profiles before and after each 20-min trial showed that no significant changes were observed, validating the use of fixed bathymetry. Alternative methods that take small-scale bed variations into account [44] are not reproducible with a numerical model such as BOSZ . For the sake of consistency in the analyses, the data from the numerical model and from the experiment were processed in an identical way, i.e., over a fixed profile. The swash time series is obtained by removing the shoreline setupη s from the shoreline elevation time series η s (t). The swash power spectral density is derived in a manner similar to that of the free surface spectra. The significant swash height S [ f l : f h ] of a given frequency band obtained between f l and f h is calculated by computing the power spectral density of the swash time series according to: where f l and f h denote the low and high frequency cutoff, respectively, and d f is the frequency resolution (first non-null value of the frequency vector). The total significant swash height is then computed according to: where S SW and S IG denote the significant swash heights computed in the SW frequency band ( f l = f p /2 and f h = 3 f p ), and in the IG frequency band ( f l = d f and f h = f p /2), respectively. Finally, the 2% runup exceedence (R 2% ) corresponding to the maximum water elevation reached by 2% of the highest runup is computed according to Stockdon et al. [22] consistently with other studies [11,28,52]:

Results
The BOSZ model computes wave transformation processes including shoaling, breaking, and energy transfer between SW and IG waves. These quantities are compared with measurements for the three test cases from Section 2.1. Moreover, the detailed scanner data are used to examine the model for the calculation of wave-induced oscillations of the shoreline water elevation and runup statistics.

Spectral Wave Characteristics
The depth-induced wave breaking between gauges WG2 and WG3 results in a significant energy dissipation that is well captured by the model for the three cases ( Figure 4). In addition, the transfer of energy from the SW frequency band to the IG band caused by time variations of the breaking point position is also accurately computed. Measurements carried out deep in the surf zone, at gauge WG4 for the intermediate beach cases DR0 and SBE2, reveal that most of the peak energy has been dissipated and the energy is now distributed between the IG and SW frequency bands (Figure 4c,f). For the reflective beach case SBA1 (Figure 4i), the energy at WG4 is distributed between the two frequency bands with a clear peak in the SW band. The detailed wave transformation patterns from the surf zone up to the swash zone for the three beach state configurations are well reproduced by the model. The model's accuracy to compute wave energy transfers in the surf and swash zones is evaluated by computing the relative error of H s at the different gauge locations that is defined as: where H p s and H m s stand for the predicted and the measured significant wave heights H s , respectively. The errors of the total significant wave height H s , the short-wave significant wave height H SW and the infragravity significant wave height H IG are summarized in Table 2. The highest error is reached for case SBE2 at WG3 due to an overestimation of H IG , while it is low for the two other cases. The IG component H IG appears to be slightly overestimated in all cases. Deeper into the surfzone, at WG4, the discrepancy reduces. The absolute errors, the difference between the modeled and measured setup, are less than 3.1 cm, which is low given the overall scale of the test. The errors are of the same order of magnitude as the precision of the initial still water level measurements.
Continuous cross-shore LiDAR measurements of the evolution of the significant wave heights associated with the different frequency bands are compared to the numerical results on Figure 5 for all cases. As Table 2 showed, H IG is slightly overestimated in the outer surfzone but reduces in the inner surfzone. The free surface profiles modeled are in good agreement with the LiDAR measurements as shown in the bottom panels of Figure 5.

Free Surface Elevation
Comparisons of the free surface time series are shown for case DR0 at the four wave gauges ( Figure 6). The agreement between physical and numerical model results is fairly good. The wave amplitudes and phases computed with the model generally match the measurements. The asymmetry of the waveform that resembles a sawtooth wave as a result of wave shoaling and the amplitude attenuation by friction and breaking are well reproduced. The results for cases SBE2 and SBA1 are of comparable agreements and are therefore not shown for brevity. The ability of the model to reproduce the evolution of the free surface along the flume as well as the shoreline motion is evaluated through the root mean square error (RMSE), bias, coefficient of determination R 2 , and Willmott's index of agreement. Willmott's index [53] is computed as: where C and O denote the computed and observed values, respectively, and n the total number of points. The agreement index has values between 0 and 1, d = 1 meaning a perfect agreement between the observed and computed values and, conversely, d = 0 indicating no agreement at all between the values. The results for the three wave gauges are summarized in Table 3. Table 3. RMSE, bias, R 2 , and Willmott's index of agreement d at the three wave gauges. The Willmott's index of 0.86 and higher confirms the satisfying agreement between computed and measured time series for all cases and at all locations. The coefficient R 2 shows a decreasing correlation between the model and the experiment in the surfzone at WG3 and WG4. A continuous cross-shore analysis of the coefficient R 2 in the LiDAR scanner area (not shown here) shows that high values of R 2 are observed prior to breaking and that the values decrease continuously throughout the surfzone. The limitation of depth-integrated equations in the breaking zone is a possible explanation for this trend. Moreover, LiDAR scanners are known to capture large vertical splashes occurring during energetic breaking events, locally leading to overestimation in the free surface measurements. A negative bias is observed for all experiments at WG3 and WG4. A possible explanation, consistent with the analysis of R 2 , is a slightly excessive dissipation of wave energy around the onset of wave breaking due to the selected condition of deactivation of the dispersion terms in the governing equations of the model.

RMSE (m) Bias (m) R 2 (-) Willmott's d (-)
The ability of the numerical model to capture second order statistics is evaluated with the comparisons of wave skewness S k and asymmetry A s , calculated as [54]: where H denotes the Hilbert transform and . denotes the time-average operator. These two quantities provide information on the wave-by-wave shape in contrast to the computed and measured comparison of the free surface elevation. Increasing values in wave skewness mean more asymmetry in the vertical direction, i.e., the wave crests increase proportionally more than the troughs decrease. Negative values in the waves' asymmetry indicate a steepening of the wave face toward the beach with a sharper leading edge and a flattened back. Especially for breaking waves, where the local wave dissipation can influence the shape of the individual waves significantly, it is difficult to obtain close agreement between measured and computed skewness and the asymmetry value. The cross-shore evolution of S k and A s are shown for all cases on Figure 7. The sharp decrease of the wave asymmetry A s , due to the steepening of the face as the depth decreases, is well captured by the model. The vertical deformation characterized by increasing values in skewness is observed similarly in both numerical and laboratory data. Considering the challenging conditions in the surf zone, the agreement between model and laboratory data is satisfying.
Overall, these comparisons show that the wave transformation processes in the surf and swash zones on intermediate and reflective beaches can be computed with a high degree of accuracy and confidence by a phase-resolving depth-integrated model such as BOSZ .

Shoreline Elevation Oscillations
The LiDAR scanner allows for accurate measurements of the shape of the runup tongue from which the water line position can be inferred. The measured and computed shoreline elevation time series η s (t) are compared on Figure 8. The model succeeds in reproducing the amplitude and phases of the shoreline oscillations. The swash energy distribution is also well captured by the model. For instance, the model shows that the swash spectrum at high frequencies exhibits a spectral decay of f −4 , which is consistent with the measured spectra and other studies [21,55]. Furthermore, for the two intermediate beach state cases DR0 ( Figure 8b) and SBE2 (Figure 8d), the swash is mainly dominated by low-frequency or infragravity oscillations, whereas for the reflective case SBA1 (Figure 8f), the SW contributions to the shoreline oscillations become more important. For all cases, the results show good agreement between the measured and modeled swash spectra. Following the statistical analysis of the wave gauges from Table 3, the RMSE, bias, R 2 coefficient, and Willmott's index for the shoreline motion η s are presented in Table 4. Willmott's index indicates a close match between the numerical and experimental time series, with values higher than 0.87. A negative bias is observed for all cases, suggesting that the model underestimates the amplitude of the oscillations. It is consistent with Table 3 where all bias at WG4 are negative. A possible explanation is a slight overdissipation of energy in the surfzone due to the breaking parametrization. Interestingly, the values presented in this table indicate a better match of η s than for η at WG4. Overall, the statistics show that the model is capable of capturing the time series oscillations with a satisfying accuracy.

Runup Statistics
To quantify the ability of the model to compute the different contributions to R 2% , SW and IG swash components (S SW and S IG , respectively) and the shoreline setupη s are computed and compared to the measured values. The results are displayed in Table 5. The relative discrepancy between the observed and computed R 2% is low, ranging between −3.4% and −6.0% which is comparable to the SWASH model performance for dissipative beach [38]. The underestimation of R 2% is consistent with Tables 3 and 4, suggesting that an overdissipation of the short wave energy in the surfzone results in underestimation of R 2% values. The infragravity component S IG is well reproduced by BOSZ for the three cases with a maximum relative error of the order of 3%. Discrepancies are more pronounced for the computation of S SW . However, the relative errors are smaller than 10%. Similarly to S SW , the shoreline setupη s displays error smaller than 10%. Despite the minor discrepancies with the LiDAR data, the model satisfactorily computes the dynamics of the shoreline elevation oscillations, including interaction between incoming bore and the receding runup. The results attest the ability of the model to compute the different components of the runup with a high degree of accuracy for the intermediate and reflective beach states considered in this study.

Model Sensitivity to the Grid Size
A sensitivity analysis of the computed R 2% and its components S SW , S IG , andη s to the grid size is conducted for all cases by varying the grid resolution from dx = 0.30 m to dx = 2 m by a 0.1 m step. Previous studies have shown a strong correlation between R 2% and the quantity √ H 0 L 0 , where H 0 is the deep-water significant wave height and and L 0 is the wavelength [22]. In order to relate the grid size to the incident offshore wave parameters, a normalized grid size based on this parameter is proposed. Considering that, according to the linear theory, L 0 = g 2π T 2 p , this normalized grid size χ can be expressed as: The model's performance is evaluated using the relative error given by Equation (12). The results are displayed on Figure 9. It is worth noting that the higher the normalized grid size χ, the finer the grid. Overall, the model's accuracy increases consistently with the grid resolution. For large grid sizes, R 2% and its components are underestimated by up to 40% for the majority of the beach states and even up to nearly 60% for the shoreline setup computed for case DR0. Only the IG swash component S IG is slightly overestimated for coarse grids in the case SBA1-a reflective beach. In general, S IG has a lower sensitivity to the grid size, which is consistent with the previous numerical study carried out with SWASH on a fringing reef [39]. The relative error curve shows an asymptotic shape from χ ∼ 200, which is close to the value corresponding to the grid size dx = 0.5 m used in the present study. To further verify the relevance of the parameter χ, synthetic cases are carried out at two different scales. A TMA spectrum is propagated over a straight slope with different H s , T p , and still water levels. The conditions are summarized in Table 6. Table 6. Synthetic tests and their graphic markers. WL: water level. Still WL (m) β f (-) ξ (-  The computed cases exhibit fairly different wave conditions to represent a wide range of scenarios and to generalize the previous finding. The evolution of R 2% and its components are displayed on Figure 10. Similarly to what is observed in Figure 9, the values tend to converge as χ increases. Again, χ ∼ 200 appears to be a reasonable value to ensure correct computation of the total runup and its components. Though this value is model-dependent to a certain extent, other numerical models of a similar kind should not behave entirely differently.

H s (m) T p (s) Intial
This result can serve as a recommendation for properly setting up a phase-resolving model for runup computation. For instance, a runup computation along a 1D transect over an intermediate to reflective beach with incident irregular waves of H 0 = 3 m and T p = 13 s offshore would require a grid size of around 2 m (foreshore slope β f ∼ 7%).

Runup Sensitivity to the Wave Breaking Detection Criterion
The simulation of wave breaking in depth-integrated wave models is generally a challenging task. Indeed, this type of model cannot solve for the free surface overturning of a breaking wave and does not include the 3D turbulent dissipation process. To overcome this limitation, the dispersive term of the governing equations can be deactivated once an onset breaking criterion is reached. Thus, the set of equations reduces to the NLSWE, which are a subset of the Boussinesq-type equations. The NLSWE have the advantage that they can describe a discontinuity in the free surface and implicitly treat the dissipation in the hydraulic jump. The solution benefits from a finite volume method such as it is the case in BOSZ . In the BOSZ version used in the present study, the onset breaking criterion is based on the free surface height to water depth ratio C b given by Equation (6). A sensitivity analysis of the runup computation to the value of C b is conducted by running the model for the three cases with C b values varying from 0.4 to 1.2 ( Figure 11). Figure 11. Evolution of the relative errors of R 2% and its components depending of the breaking coefficient C b for case DR0 (blue stars), SBE2 (green plus signs), and SBA1 (black circles).
For case DR0, the intermediate beach configuration without inner bar, R 2% and its components are more sensitive to C b than for the two other cases. For this case, relative errors increase with C b except for the IG swash component, for which the relative error decreases. In the SBA1 and SBE2 cases, the swash components are insensitive to C b . In contrast, the relative error of the mean shoreline setup increases linearly with C b , a trend that can be seen in the relative error evolution of R 2% . Overall, the best results for all of the components are obtained for C b ranging between 0.6 and 0.8, which is consistent with the value used in this study (C b = 0.78), based the theoretical work from McCowan [50], and with other studies [42,56].

Sensitivity of the Runup Determination to the Threshold Depth
The determination of the leading edge of runup requires to define an ad hoc criterion. In the present study, the threshold depth δ used to track the limit between dry and wet cells was set to 10 cm according to recommendations from previous studies based on field data [57,58] and numerical results [7,12]. In this section, the influence of the value of δ on R 2 % and its components is investigated for both the numerical results and the laboratory data by varying δ from 4 cm up to 20 cm ( Figure 12). In general, a threshold depth resembles a low-pass filter. Overall, low values of δ result in higher runup, swash, and shoreline setup values regardless of the beach type, even if this trend is less pronounced for the determination of S SW in case DR0. Moreover, for small δ lower than 6 cm, it is worth noting that runup components computed from laboratory data are particularly variable, especially for the SBE2 case. It is hypothesized that for this case, small changes in the measured bed profile of the order of 5 cm have an influence on the determination of the runup tongue. This tendency is not observed in the numerical results. In fact, all the computations were carried out using a fixed bottom. Additional tests (not shown here) reveal that when the experimental shoreline position is extracted using a variable profile, the runup components show a behavior similar to that of the numerical results, confirming the influence of small changes of the foreshore slope when using low values of δ to identify the runup limit. Finally, the 10 cm threshold value used in this study provides similar runup values between experimental data and numerical results, which is consistent with other studies (e.g., [22,57]).

Influence of Phase Distribution of Incident Waves on Runup
For the comparisons between measured and computed runup, the BOSZ model was forced using the free surface times series measured at the wave gauge WG1. The results confirm the ability of the model to accurately compute the oscillations of the shoreline elevation in case the phase distribution is known (Figure 8). However, for practical applications of the model, an empirical offshore wave spectrum or a spectrum from a phase-average model such as SWAN or WaveWATCH III is commonly used to generate the input free surface elevation time series. In this case, the phases are not known and a distribution of random phases is used instead. Obviously, these random phases remain fixed throughout the computation.
The influence of incident phase distribution on the assessment of overtopping volumes has been highlighted in previous studies [42,59,60], showing that the overtopping volume can be overestimated by more than 100%. Wave runup sensitivity to phasing of incident waves was also studied for sloping beach [28,61] and fringing reef environments [32], showing significant differences between R 2% computed with measured phases and random phases. In particular, the significant IG wave height computed near the shoreline of an idealized fringing reef profile was overestimated by 20% using randomly distributed phases.
The influence of the phase distribution on the shoreline elevation is investigated here for the three test cases. Ten runs are computed with the same input energy spectrum as the one measured, but with different uniform random phase distributions within the interval [−π : +π].
The runup components, normalized by the mean value of the ten runs, are compared in Figure 13 to assess the variability of the runup in dependence of the random phase distribution. The results show that S IG is the most sensitive component of the runup with significant variations in all cases. The total IG signal is composed of both bound and free long waves, where bound waves are phased-locked to the wave group traveling at the group speed [62]. The wave groups are the result of the effect of the dispersion relation on the initial superposition of the different spectral components; their shape and appearance in time depend essentially on the initial phase distribution. The bound components in the IG band are directly affected by the wave phases. Consequently, changes in the phases lead to some variability of the bound waves, which are released as free IG waves after wave breaking. Moreover, Yao et al. [32] suggested that the short waves envelope, which is highly dependent on the interaction of the swell waves with each other, partially controls the IG wave transformations. On the other side, S SW and the averaged quantity of setup at the shorelineη s change only slightly with the phase angles. The SW swash and setup are mainly controlled by the depth-induced breaking of individual short waves and might be less sensitive to the SW envelope and, therefore, to the phase distribution. As reported by Torres-Freyermuth et al. [28], the variability increases with decreasing values of ξ, suggesting that the runup variability is larger on dissipative beaches. One possible reason for that is the increasing contribution of the IG band to the runup on dissipative beaches.
The relative errors between R 2% and its components from the random phases and the laboratory data are summarized in Table 7. The errors in R 2% are much larger than when the measured time series is used directly as input in the model (see Table 5). On average, the use of random phases overestimates the measured runup, with a maximum overestimation of 35% for the DR0 case. This highlights the aleatory nature of the runup and the notorious difficulties of comparing numerical results to measured laboratory or field data when no phase information is available. Furthermore, the IG swash is largely overestimated in all cases.

Conclusions
As phase-resolving depth-integrated models are gaining importance for runup studies, precise validations under controlled conditions and recommendations for best practice are required. In this work, a high-quality laboratory data set is used to investigate the capability and sensitivity of the Boussinesq-type model BOSZ for the computation of nearshore wave transformations, including swash processes over intermediate and reflective beaches. The data set includes accurate LiDAR measurements of the free surface elevation in the surf and swash zone as well as shoreline elevation oscillations.
Wave transformations are accurately captured with low cross-shore errors for both the significant wave height H s and the wave setupη. Time series from the numerical model output of shoreline elevation oscillations as well as swash spectra show a satisfying agreement with the laboratory data. The statistical runup quantity R 2% is successfully computed with relative errors of less than 6%. The IG swash is well predicted with errors smaller than 3.5%. The SW swash and shoreline setup are reasonably well predicted with errors of less than 10%.
The discussion evaluates the sensitivity of the results to the model settings for general numerical computations of wave runup by depth-integrated phase-resolving models. Multiple computations with different breaking indices show that in the range of C b = 0.6 to 0.8, the numerical results show little variability overall. Some parameters, such as the grid size or the threshold depth defining the edge of the runup tongue, are found to have a significant impact on the results, and thus the performance of the model. A nondimensional parameter is proposed to find the optimal grid size to improve numerical accuracy. A depth threshold of 10 cm, consistent with other studies [11,22], is found to be the most appropriate value for systematic comparisons of numerical and laboratory data to prevent small changes in the beach profile from having a disproportional impact. For model/data comparison, a free surface time series is used as the boundary condition for the numerical model, thus providing information of both amplitude and phase angles. Computations with different sets of random phases demonstrate that accurate replication of the laboratory data can only be achieved when the exact phases are not known. Moreover, a significant variability is observed among the runs with different random phases due to the sensitivity of the IG swash to the initial phase distribution. For beaches with high influence of IG energy, i.e., intermediate to dissipative beaches, the variability of the runup can be significant. If the goal is to reproduce particular runup values, which were previously measured in the laboratory or field, the lack of information of the incident wave phases can contribute to substantial uncertainty. It should also be noted that special attention is necessary when data from laboratory experiments such as the one in this study are used. The generation of irregular waves in a laboratory environment essentially relies on the same technique as that which is used in phase-resolving models, i.e., a wave spectrum is decomposed into a series of individual waves. Laboratory data are inevitably subject to the same problem related to the waves' phases as phase-resolving wave models.
Overall, the phase-resolving depth-integrated BOSZ model shows satisfying capabilities in modeling irregular wave runup on intermediate and reflective beaches. It proves that this type of numerical model can be a powerful tool for coastal risks assessment and hazard mitigation projects. The sensitivity analysis performed provides guidelines on how to utilize the model and, more generally, any phase-resolving depth-integrated model to find the best accuracy at the lowest computational cost and ensure quality results for runup modeling studies.