Next Article in Journal
Multi-Step-Ahead Forecasting of Wave Conditions Based on a Physics-Based Machine Learning (PBML) Model for Marine Operations
Next Article in Special Issue
Observed Changes of a Mega Feeder Nourishment in a Coastal Cell: Five Years of Sand Engine Morphodynamics
Previous Article in Journal
Research on Real-Time Optimal Path Planning Model and Algorithm for Ship Block Transportation in Shipyard
Previous Article in Special Issue
Modeling Impact of Intertidal Foreshore Evolution on Gravel Barrier Erosion and Wave Runup with XBeach-X
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Université de Pau et des Pays de l’Adour, E2S UPPA, chair HPC-Waves, SIAME, 64600 Anglet, France
*
Author to whom correspondence should be addressed.
J. Mar. Sci. Eng. 2020, 8(12), 993; https://doi.org/10.3390/jmse8120993
Submission received: 15 October 2020 / Revised: 25 November 2020 / Accepted: 1 December 2020 / Published: 5 December 2020
(This article belongs to the Special Issue Observation, Analysis, and Modeling of Nearshore Dynamics)

Abstract

:
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 R 2 % 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 S I G is found to be sensitive to the initial phase distribution, likely because it is related to the short wave envelope.

1. 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.

2. Methodology

2.1. 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.
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.

2.2. 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:
ξ = β f H 0 / L 0
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.

2.3. 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 k h < π or a little beyond that. For k h > π , 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 x- and y-directions. In 1D, only the terms in the x-direction are retained:
H t + ( H u ) x + ψ C + ψ w m = 0
( H u ) t + H z α 2 2 u x x t + H z α ( h u t ) x x + ( H u 2 ) x + g H η x + u ψ C + τ 1 = 0
ψ C = z α 2 2 h 2 6 h u x x + z α + h 2 h h u x x x
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 ( H u ) instead of with only u:
H ( u t + u u x + g η x ) = ( H u ) t + ( H u 2 ) x + H g η x u [ H t + ( H u ) x ] = ( H u ) t + ( H u 2 ) x + H g η x + u [ ψ C ] ,
as H t + ( H u ) 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, ψ w m , 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:
η h > C b
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.

2.4. 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 d x = 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:
U = η . g h
where U is the horizontal particle velocity. Though the time series clearly shows dispersive waves, the selected conditions are of relatively low frequency dispersion ( k h < 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.

2.5. 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:
H s = 4 m 0
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:
S [ f l : f h ] = 4 f l f h P S D ( f ) d f
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:
S = S S W 2 + S I G 2
where S S W and S I G 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]:
R 2 % = 1.1 ( η s ¯ + S / 2 )

3. 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.

3.1. Wave Propagation

3.1.1. 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:
E ( % ) = 100 H s p H s m H s m
where H s p and H s m 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 S W and the infragravity significant wave height H I G are summarized in Table 2. The highest error is reached for case SBE2 at WG3 due to an overestimation of H I G , while it is low for the two other cases. The IG component H I G 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 I G 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.

3.1.2. 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:
d = 1 i = 1 n ( C i O i ) 2 i = 1 n ( | C i O ¯ | + | O i O ¯ | ) 2
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.
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.
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]:
S k = ( η η ¯ ) 3 ( η η ¯ ) 2 3 / 2
S k = H ( η η ¯ ) 3 ( η η ¯ ) 2 3 / 2
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.

3.2. Wave Runup

3.2.1. 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.

3.2.2. Runup Statistics

To quantify the ability of the model to compute the different contributions to R 2 % , SW and IG swash components ( S S W and S I G , 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 Table 3 and Table 4, suggesting that an overdissipation of the short wave energy in the surfzone results in underestimation of R 2 % values. The infragravity component S I G 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 S W . However, the relative errors are smaller than 10 % . Similarly to S S W , 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.

4. Discussion

4.1. Model Sensitivity to the Grid Size

A sensitivity analysis of the computed R 2 % and its components S S W , S I G , and η s ¯ to the grid size is conducted for all cases by varying the grid resolution from d x = 0.30 m to d x = 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 p 2 , this normalized grid size χ can be expressed as:
χ = H 0 L 0 ( d x . β f ) = g 2 π H 0 T p 2 ( d x . β f )
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 I G is slightly overestimated for coarse grids in the case SBA1—a reflective beach. In general, S I G 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 d x = 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.
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.
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%).

4.2. 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).
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].

4.3. 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 S W 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]).

4.4. 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 I G 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 S W 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.

5. 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.

Author Contributions

Conceptualization, J.P. and D.M.; methodology, J.P.; software and model development, V.R.; validation, J.P.; formal analysis, J.P.; writing—original draft preparation, J.P.; writing—review and editing, D.M. and V.R.; visualization, J.P.; supervision, D.M. and V.R.; project administration, V.R.; funding acquisition, D.M. and V.R. All authors have read and agreed to the published version of the manuscript.

Funding

The authors acknowledge financial support from the Isite program Energy Environment Solutions (E2S), the Communauté d’Agglomération Pays Basque (CAPB), and the Communauté Région Nouvelle Aquitaine (CRNA) for the chair HPC-Waves.

Acknowledgments

The authors express appreciation to Chris Blenkinsopp and Paul Bayle for providing data from the DynaRev project and guidance. The DynaRev project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 654110, HYDRALAB+. The authors gratefully acknowledge European POCTEFA Program funding under the research project MARLIT EFA344/19.

Conflicts of Interest

The authors declare no conflict of interest. The funding agencies had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
TWLTotal Water Level
SWShort Waves
IGInfra-Gravity
LiDARLight Detection And Ranging
NLSWENon Linear Shallow Water Equations
SWASHSimulating WAves till SHore
BOSZBoussinesq Ocean & Surf Zone
SWANSimulating WAves Nearshore

References

  1. Didier, D.; Baudry, J.; Bernatchez, P.; Dumont, D.; Sadegh, M.; Bismuth, E.; Bandet, M.; Dugas, S.; Sévigny, C. Multihazard simulation for coastal flood mapping: Bathtub versus numerical modelling in an open estuary, Eastern Canada. J. Flood Risk Manag. 2019, 12, e12505. [Google Scholar] [CrossRef]
  2. Xie, D.; Zou, Q.P.; Mignone, A.; MacRae, J.D. Coastal flooding from wave overtopping and sea level rise adaptation in the northeastern USA. Coast. Eng. 2019, 150, 39–58. [Google Scholar] [CrossRef]
  3. Sallenger, A. Storm impact scale for barrier islands. J. Coast. Res. 2000, 16, 890–895. [Google Scholar]
  4. Van der Meer, J.W.; Allsop, N.; Bruce, T.; De Rouck, J.; Kortenhaus, A.; Pullen, T.; Schüttrumpf, H.; Troch, P.; Zanuttigh, B. EurOtop 2018. Manual on Wave Overtopping of Sea Defences and Related Structures. An Overtopping Manual Largely Based on European Research, but for Worldwide Application, 2nd ed.; EurOtop: London, UK, 2018; p. 320. [Google Scholar]
  5. Ruggiero, P.; Komar, P.D.; McDougal, W.G.; Marra, J.J.; Beach, R.A. Wave Runup, Extreme Water Levels and the Erosion of Properties Backing Beaches. J. Coast. Res. 2001, 17, 407–419. [Google Scholar]
  6. Serafin, K.A.; Ruggiero, P. Simulating extreme total water levels using a time-dependent, extreme value approach. J. Geophys. Res. Ocean. 2014, 119, 6305–6329. [Google Scholar] [CrossRef] [Green Version]
  7. Stockdon, H.; Thompson, D.; Plant, N.; Long, J. Evaluation of wave runup predictions from numerical and parametric models. Coast. Eng. 2014, 92, 1–11. [Google Scholar] [CrossRef]
  8. Silveira, T.M.; Taborda, R.; Carapuço, M.M.; Andrade, C.; Freitas, M.C.; Duarte, J.F.; Psuty, N.P. Assessing the extreme overwash regime along an embayed urban beach. Geomorphology 2016, 274, 64–77. [Google Scholar] [CrossRef]
  9. Medellín, G.; Brinkkemper, J.; Torres-Freyermuth, A.; Appendini, C.; Mendoza, E.; Salles, P. Run-up parameterization and beach vulnerability assessment on a barrier island: A downscaling approach. Nat. Hazards Earth Syst. Sci. 2016, 16, 167–180. [Google Scholar] [CrossRef] [Green Version]
  10. Lerma, A.N.; Pedrero, R.; Robinet, A.; Sénéchal, N. Simulating wave setup and runup during storm conditions on a complex barred beach. Coast. Eng. 2017, 123, 29–41. [Google Scholar] [CrossRef]
  11. Fiedler, J.; Smit, P.; Brodie, K.; McNinch, J.; Guza, R. Numerical modeling of wave runup on steep and mildly sloping natural beaches. Coast. Eng. 2018, 131, 106–113. [Google Scholar] [CrossRef]
  12. Fiedler, J.; Young, A.; Ludka, B.; O’Reilly, W.; Henderson, C.; Merrifield, M.; Guza, R. Predicting site-specific storm wave run-up. Nat. Hazards 2020, 104, 493–517. [Google Scholar] [CrossRef]
  13. Holman, R. Extreme value statistics for wave run-up on a natural beach. Coast. Eng. 1986, 9, 527–544. [Google Scholar] [CrossRef]
  14. Longuet-Higgins, M.; Stewart, R. Radiation stresses in water waves; a physical discussion, with applications. Deep Sea Res. Oceanogr. Abstr. 1964, 11, 529–562. [Google Scholar] [CrossRef]
  15. De Bakker, A.; Tissier, M.; Ruessink, B. Beach steepness effects on nonlinear infragravity-wave interactions: A numerical study. J. Geophys. Res. Ocean. 2016, 121, 554–570. [Google Scholar] [CrossRef] [Green Version]
  16. Hunt, I. Design of seawalls and breakwaters, US Corps Eng. Lake Surv. Detroit 1959, 49, 123–152. [Google Scholar]
  17. Battjes, J. Surf similarity. In Proceedings of the 14th International Conference on Coastal Engineering, Copenhagen, Denmark, 24–28 June 1974; Volume 9, pp. 446–480. [Google Scholar]
  18. Mase, H. Random Wave Runup Height on Gentle Slope. J. Waterw. Port Coast. Ocean Eng. ASCE 1989, 115, 649–661. [Google Scholar] [CrossRef]
  19. Hedges, T.; Mase, H. Modified Hunt’s Equation Incorporating Wave Setup. J. Waterw. Port Coast. Ocean Eng. ASCE 2004, 130, 109–113. [Google Scholar] [CrossRef]
  20. Khoury, A.; Jarno, A.; Marin, F. Experimental study of runup for sandy beaches under waves and tide. Coast. Eng. 2019, 144, 33–46. [Google Scholar] [CrossRef]
  21. Ruggiero, P.; Holman, R.; Beach, R. Wave run-up on a high-energy dissipative beach. J. Geophys. Res. C Ocean. 2004, 109. [Google Scholar] [CrossRef] [Green Version]
  22. Stockdon, H.; Holman, R.; Howd, P.; Sallenger, A. Empirical parameterization of setup, swash, and runup. Coast. Eng. 2006, 53, 573–588. [Google Scholar] [CrossRef]
  23. Atkinson, A.; Power, H.; Moura, T.; Hammond, T.; Callaghan, D.; Baldock, T. Assessment of runup predictions by empirical models on non-truncated beaches on the south-east Australian coast. Coast. Eng. 2017, 119, 15–31. [Google Scholar] [CrossRef] [Green Version]
  24. Didier, D.; Caulet, C.; Bandet, M.; Bernatchez, P.; Dumont, D.; Augereau, E.; Floc’h, F.; Delacourt, C. Wave runup parameterization for sandy, gravel and platform beaches in a fetch-limited, large estuarine system. Cont. Shelf Res. 2020, 192, 104024. [Google Scholar] [CrossRef]
  25. Vousdoukas, M.; Wziatek, D.; Almeida, L. Coastal vulnerability assessment based on video wave run-up observations at a mesotidal, steep-sloped beach. Ocean Dyn. 2011, 62, 123–137. [Google Scholar] [CrossRef]
  26. Da Silva, G.; Gomes da Silva, P.; Araujo, R.; Klein, A.; Toldo, E. Wave run-up on embayed beaches. Study case: Itapocorói Bay, Southern Brazil. Braz. J. Oceanogr. 2017, 65, 187–200. [Google Scholar] [CrossRef] [Green Version]
  27. Senechal, N.; Coco, G.; Bryan, K.; Holman, R. Wave runup during extreme storm conditions. J. Geophys. Res. Ocean. 2011, 116. [Google Scholar] [CrossRef] [Green Version]
  28. Torres-Freyermuth, A.; Patiño, J.P.; Pedrozo-Acuña, A.; Puleo, J.; Baldock, T. Runup uncertainty on planar beaches. Ocean Dyn. 2019, 69, 1359–1371. [Google Scholar] [CrossRef]
  29. Martin-Medina, M.; Abadie, S.; Mokrani, C.; Morichon, D. Numerical simulation of flip-through impacts of variable steepness on a vertical breakwater. Appl. Ocean Res. 2018, 75, 117–131. [Google Scholar] [CrossRef]
  30. Zijlema, M.; Stelling, G.; Smit, P. SWASH: An operational public domain code for simulating wave fields and rapidly varied flows in coastal waters. Coast. Eng. 2011, 58, 992–1012. [Google Scholar] [CrossRef]
  31. Lynett, P.; Tso-Ren, W.; Liu, P.F. Modeling wave runup with depth-integrated equations. Coast. Eng. 2002, 46, 89–107. [Google Scholar] [CrossRef]
  32. Yao, Y.; Zhang, Q.; Becker, J.M.; Merrifield, M.A. Boussinesq modeling of wave processes in field fringing reef environments. Appl. Ocean Res. 2020, 95, 102025. [Google Scholar] [CrossRef]
  33. Roeber, V.; Cheung, K.; Kobayashi, M. Shock-capturing Boussinesq-type model for nearshore wave processes. Coast. Eng. 2010, 57, 407–423. [Google Scholar] [CrossRef]
  34. Roeber, V.; Cheung, K. Boussinesq-type model for energetic breaking waves in fringing reef environments. Coast. Eng. 2012, 70, 1–20. [Google Scholar] [CrossRef]
  35. Li, N.; Roeber, V.; Yamazaki, Y.; Heitmann, T.W.; Bai, Y.; Cheung, K.F. Integration of coastal inundation modeling from storm tides to individual waves. Ocean Model. 2014, 83, 26–42. [Google Scholar] [CrossRef]
  36. Li, N.; Yamazaki, Y.; Roeber, V.; Cheung, K.F.; Chock, G. Probabilistic mapping of storm-induced coastal inundation for climate change adaptation. Coast. Eng. 2018, 133, 126–141. [Google Scholar] [CrossRef]
  37. Ruessink, B.G.; Michallet, H.; Bonneton, P.; Mouazé, D.; Lara, J.; Silva, P.A.; Wellens, P. Globex: Wave dynamics on a gently sloping laboratory beach. In Proceedings of the HYDRALAB IV Joint User Meeting, Lisbon, Portugal, 1–9 July 2014; pp. 1351–1362. [Google Scholar]
  38. Ruju, A.; Lara, J.; Losada, I. Numerical analysis of run-up oscillations under dissipative conditions. Coast. Eng. 2014, 86, 45–56. [Google Scholar] [CrossRef]
  39. Pelaez-Zapata, D.; Montoya, R.; Osorio, A. Numerical study of run-up oscillations over fringing reefs. J. Coast. Res. 2018, 34, 1065–1079. [Google Scholar] [CrossRef]
  40. Liu, W.; Shao, K.q.; Ning, Y. Numerical Study of the Impact of Climate Change on Irregular Wave Run-up Over Reef-Fringed Coasts. China Ocean Eng. 2020, 34, 162–171. [Google Scholar] [CrossRef]
  41. Mase, H.; Miyahira, A.; Hedges, T.S. Random Wave Runup on Seawalls Near Shorelines with and without Artificial Reefs. Coast. Eng. J. 2004, 46, 247–268. [Google Scholar] [CrossRef]
  42. McCabe, M.; Stansby, P.; Apsley, D. Random wave runup and overtopping a steep sea wall: Shallow-water and Boussinesq modelling with generalised breaking and wall impact algorithms validated against laboratory and field measurements. Coast. Eng. 2013, 74, 33–49. [Google Scholar] [CrossRef]
  43. Blenkinsopp, C.; Mole, M.; Turner, I.; Peirson, W. Measurements of the time-varying free-surface profile across the swash zone obtained using an industrial LIDAR. Coast. Eng. 2010, 57, 1059–1065. [Google Scholar] [CrossRef]
  44. Brodie, K.; Raubenheimer, B.; Elgar, S.; Slocum, R.K.; McNinch, J. Lidar and pressure measurements of inner-surfzone waves and setup. J. Atmos. Ocean. Technol. 2015, 32, 1945–1959. [Google Scholar] [CrossRef]
  45. Fiedler, J.; Brodie, K.; McNinch, J.; Guza, R. Observations of runup and energy flux on a low slope beach with high-energy, long-period ocean swell. Geophys. Res. Lett. 2015, 42. [Google Scholar] [CrossRef] [Green Version]
  46. Blenkinsopp, C.; Bayle, P.; Conley, D.; Masselink, G.; Gulson, E.; Kelly, I.; Almar, R.; Turner, I.; Baldock, T.; Beuzen, T.; et al. Dynamic coastal protection: Resilience of dynamic revetments (DYNAREV). In Proceedings of the HYDRALAB+ Joint User Meeting, Bucharest, Romania, 21–25 May 2019. [Google Scholar]
  47. Wright, L.; Short, A. Morphodynamic variability of surf zones and beaches: A synthesis. Mar. Geol. 1984, 56, 93–118. [Google Scholar] [CrossRef]
  48. Nwogu, O. An alternative form of the Boussinesq equations for nearshore wave propagation. J. Waterw. Port Coast. Ocean Eng. 1993, 119, 618–638. [Google Scholar] [CrossRef] [Green Version]
  49. Wei, G.; Kirby, J.T.; Sinha, A. Generation of waves in Boussinesq models using a source function method. Coast. Eng. 1999, 36, 271–299. [Google Scholar] [CrossRef]
  50. McCowan, J. XXXIX. On the highest wave of permanent type. Lond. Edinb. Dublin Philos. Mag. J. Sci. 1984, 38, 351–358. [Google Scholar] [CrossRef] [Green Version]
  51. Emery, W.J.; Thomson, R.E. Time-series Analysis Methods. In Data Analysis Methods in Physical Oceanography; Emery, W.J., Thomson, R.E., Eds.; Elsevier Science: Amsterdam, The Netherlands, 2001; Chapter 5; pp. 371–567. [Google Scholar] [CrossRef]
  52. De Beer, A.; McCall, R.; Long, J.; Tissier, M.; Reniers, A. Simulating wave runup on an intermediate–reflective beach using a wave-resolving and a wave-averaged version of XBeach. Coast. Eng. 2020, 163, 103788. [Google Scholar] [CrossRef]
  53. Willmott, C.J. On the Validation of Models. Phys. Geogr. 1981, 2, 184–194. [Google Scholar] [CrossRef]
  54. Elgar, S. Observations of bispectra of shoaling surface gravity waves. J. Fluid Mech. 1985, 161, 425–448. [Google Scholar] [CrossRef]
  55. Hughes, M.; Aagaard, T.; Baldock, T.; Power, H. Wave runup (Swash) spectra on natural beaches: Morphodynamic controls. In Proceedings of the Coasts Ports 2013, Sydney, Australia, 11–13 September 2013; pp. 412–417. [Google Scholar]
  56. Tonelli, M.; Petti, M. Hybrid finite volume – finite difference scheme for 2DH improved Boussinesq equations. Coast. Eng. 2009, 56, 609–620. [Google Scholar] [CrossRef]
  57. Holland, K.; Raubenheimer, B.; Guza, R.; Holman, R. Runup kinematics on a natural beach. J. Geophys. Res. 1995, 100, 4985–4993. [Google Scholar] [CrossRef]
  58. Raubenheimer, B.; Guza, R.; Elgar, S. Swash on a gently sloping beach. J. Geophys. Res. 1995, 100, 8751–8760. [Google Scholar] [CrossRef]
  59. Romano, A.; Bellotti, G.; Briganti, R.; Franco, L. Uncertainties in the physical modelling of the wave overtopping over a rubble mound breakwater: The role of the seeding number and of the test duration. Coast. Eng. 2015, 103, 15–21. [Google Scholar] [CrossRef]
  60. Williams, H.E.; Briganti, R.; Pullen, T. The role of offshore boundary conditions in the uncertainty of numerical prediction of wave overtopping using non-linear shallow water equations. Coast. Eng. 2014, 89, 30–44. [Google Scholar] [CrossRef]
  61. McCabe, M.; Stansby, P.K.; Apsley, D.D. Coupled wave action and shallow-water modelling for random wave runup on a slope. J. Hydraul. Res. 2011, 49, 515–522. [Google Scholar] [CrossRef]
  62. Longuet-Higgins, M.S.; Stewart, R.W. Radiation stress and mass transport in gravity waves, with application to ‘surf beats’. J. Fluid Mech. 1962, 13, 481–504. [Google Scholar] [CrossRef]
Figure 1. Bathymetric profile for case DR0 (see 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.
Figure 1. Bathymetric profile for case DR0 (see 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.
Jmse 08 00993 g001
Figure 2. Definition sketch of the free surface flows.
Figure 2. Definition sketch of the free surface flows.
Jmse 08 00993 g002
Figure 3. LiDAR scanner flow depth measurement using a threshold depth of 10 cm.
Figure 3. LiDAR scanner flow depth measurement using a threshold depth of 10 cm.
Jmse 08 00993 g003
Figure 4. Comparisons between the observed (red dashed line; mostly obstructed) and computed (blue solid line) spectra at the three wave gauge locations for case DR0 (ac), SBE2 (df), and SBA1 (gi). The vertical dotted line represents the peak frequency and the black dashed line represents the frequency cutoff between the SW and IG band. Shaded areas define the 95% confidence intervals. The bathymetric profile was assumed fixed in the numerical computation for each of the three scenarios.
Figure 4. Comparisons between the observed (red dashed line; mostly obstructed) and computed (blue solid line) spectra at the three wave gauge locations for case DR0 (ac), SBE2 (df), and SBA1 (gi). The vertical dotted line represents the peak frequency and the black dashed line represents the frequency cutoff between the SW and IG band. Shaded areas define the 95% confidence intervals. The bathymetric profile was assumed fixed in the numerical computation for each of the three scenarios.
Jmse 08 00993 g004
Figure 5. Cross-shore evolution of computed (blue solid line) and observed (red dashed lines) significant wave heights. Values at WG2 are shown as circles at x = 170 m. (ac): H S W of cases DR0, SBE2, and SBA1. (df): H I G of cases DR0, SBE2, and SBA1. (gi): bed profiles with modeled (blue) and LiDAR measurements (red) of the instantaneous free surface at t = 961 s.
Figure 5. Cross-shore evolution of computed (blue solid line) and observed (red dashed lines) significant wave heights. Values at WG2 are shown as circles at x = 170 m. (ac): H S W of cases DR0, SBE2, and SBA1. (df): H I G of cases DR0, SBE2, and SBA1. (gi): bed profiles with modeled (blue) and LiDAR measurements (red) of the instantaneous free surface at t = 961 s.
Jmse 08 00993 g005
Figure 6. Comparisons between observed (red dashed line) and computed (blue solid line) free surface time series at different cross-shore locations for case DR0.
Figure 6. Comparisons between observed (red dashed line) and computed (blue solid line) free surface time series at different cross-shore locations for case DR0.
Jmse 08 00993 g006
Figure 7. Cross-shore evolution of the computed (blue solid line) and observed (red circles) skewness S k and asymmetry A s for all cases. Values at WG2 are shown as circles at x = 170 m. (ac): S k of cases DR0, SBE2, and SBA1. (df): A s of cases DR0, SBE2, and SBA1.
Figure 7. Cross-shore evolution of the computed (blue solid line) and observed (red circles) skewness S k and asymmetry A s for all cases. Values at WG2 are shown as circles at x = 170 m. (ac): S k of cases DR0, SBE2, and SBA1. (df): A s of cases DR0, SBE2, and SBA1.
Jmse 08 00993 g007
Figure 8. Comparisons between the observed (red) and simulated (blue) continuous shoreline elevation η s and swash spectra for cases DR0 (a,b), SBE2 (c,d), and SBA1 (e,f). Blacked dashed line: IG cutoff frequency ( f p / 2 ). Black dotted line: peak frequency f p . Shaded areas define the 95% confidence intervals.
Figure 8. Comparisons between the observed (red) and simulated (blue) continuous shoreline elevation η s and swash spectra for cases DR0 (a,b), SBE2 (c,d), and SBA1 (e,f). Blacked dashed line: IG cutoff frequency ( f p / 2 ). Black dotted line: peak frequency f p . Shaded areas define the 95% confidence intervals.
Jmse 08 00993 g008
Figure 9. Relative error of the runup components between the numerical results and the laboratory experiment in function χ for case DR0 (blue stars), SBE2 (green plus signs), and SBA1 (black circles). The corresponding χ to d x = 0.5 m is shown in red. Note that χ increases for decreasing values of d x .
Figure 9. Relative error of the runup components between the numerical results and the laboratory experiment in function χ for case DR0 (blue stars), SBE2 (green plus signs), and SBA1 (black circles). The corresponding χ to d x = 0.5 m is shown in red. Note that χ increases for decreasing values of d x .
Jmse 08 00993 g009
Figure 10. Self-contained deviation of the runup components from the results obtained with the smallest grid size as a function of χ . For markers and curves colors, refer to Table 6.
Figure 10. Self-contained deviation of the runup components from the results obtained with the smallest grid size as a function of χ . For markers and curves colors, refer to Table 6.
Jmse 08 00993 g010
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).
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).
Jmse 08 00993 g011
Figure 12. Experimental (dashed line) and numerical (solid line) runup components depending on the threshold depth δ for the DR0 (blue stars), SBE2 (green plus signs), and SBA1 (black circles) cases.
Figure 12. Experimental (dashed line) and numerical (solid line) runup components depending on the threshold depth δ for the DR0 (blue stars), SBE2 (green plus signs), and SBA1 (black circles) cases.
Jmse 08 00993 g012
Figure 13. Runup components normalized by the mean value for the 10 random tests for case DR0 (a), case SBE2 (b), and case SBA1 (c). Red: R 2 % , green: η s ¯ , blue: S S W , and black: S I G .
Figure 13. Runup components normalized by the mean value for the 10 random tests for case DR0 (a), case SBE2 (b), and case SBA1 (c). Red: R 2 % , green: η s ¯ , blue: S S W , and black: S I G .
Jmse 08 00993 g013
Table 1. Three beach state cases (I: intermediate, R: reflective). WL: water level.
Table 1. Three beach state cases (I: intermediate, R: reflective). WL: water level.
Case Name H s (m) T p (s) β f (-) ξ (-)Initial Still WL (m)Ursell Number U r (-)State
DR00.866%0.504.511.7I
SBE21.2810%0.914.928.2I
SBA10.61210%1.934.933.7R
Table 2. Relative error of H s , H S W , and H I G , and absolute errors of the setup η ¯ at the different wave gauges locations. The still water level is indicated for reference.
Table 2. Relative error of H s , H S W , and H I G , and absolute errors of the setup η ¯ at the different wave gauges locations. The still water level is indicated for reference.
Case Name WG2WG3WG4
DR0depth (m)410.30
H s (%)−1.95%−2.52%−1.36%
H S W (%)−1.70%−6.70%−3.17%
H I G (%)+5.24%+23.3%+3.25%
η ¯ (cm)−1.0−2.0−3.1
SBE2depth (m)4.510.45
H s (%)1.28%8.51%6.90%
H S W (%)+0.58%−1.23%−2.47%
H I G (%)+13.4%+19.5%+16.6%
η ¯ (cm)2.7−0.3−1.4
SBA1depth (m)4.510.45
H s (%)−1.14%−0.55%−0.91%
H S W (%)−2.01%−3.06%−6.06%
H I G (%)+10.4%+23.9%+9.58%
η ¯ (cm)0.08−1.5−2.6
Table 3. RMSE, bias, R 2 , and Willmott’s index of agreement d at the three wave gauges.
Table 3. RMSE, bias, R 2 , and Willmott’s index of agreement d at the three wave gauges.
RMSE (m)Bias (m) R 2 (-)Willmott’s d (-)
DR0
WG20.06−0.010.930.98
WG30.09−0.020.670.90
WG40.09−0.040.600.87
SBE2
WG20.110.030.900.97
WG30.17−0.0030.550.86
WG40.13−0.010.600.88
SBA1
WG20.040.010.950.99
WG30.1−0.010.740.93
WG40.1−0.030.680.89
Table 4. RMSE, bias, R 2 , and Willmott’s index of agreement d for η s .
Table 4. RMSE, bias, R 2 , and Willmott’s index of agreement d for η s .
Case NameRMSE (m)Bias (m) R 2 (-)Willmott’s d (-)
DR00.07−0.010.710.87
SBE20.13−0.020.790.91
SBA10.08−0.010.890.96
Table 5. Relative errors of the R 2 % and the different runup components. Values from the numerical model are shown after the vertical bar.
Table 5. Relative errors of the R 2 % and the different runup components. Values from the numerical model are shown after the vertical bar.
Case Name R 2 % S SW S IG η s ¯
DR0−3.4%0.44 m−1.2%0.28 m−3.3%0.42 m−4.7%0.15 m
SBE2−6.0%0.84 m−9.8%0.63 m−1.7%0.86 m−8.7%0.23 m
SBA1−3.9%0.65 m−5.9%0.77 m+2.5%0.47 m−4.4%0.14 m
Table 6. Synthetic tests and their graphic markers. WL: water level.
Table 6. Synthetic tests and their graphic markers. WL: water level.
H s (m) T p (s)Intial Still WL (m) β f (-) ξ (-)Marker
L1.1313205%0.47blue stars
L1.2515205%0.42blue circles
L2.13132010%0.94green stars
L2.25152010%0.84green circles
S1.10.152.915%0.47black stars
S1.20.253.3515%0.42black circles
S2.10.152.9110%0.94red stars
S2.20.253.35110%0.82red circles
Table 7. Maximum and mean relative errors of the computed scenarios in comparison with the laboratory experiment.
Table 7. Maximum and mean relative errors of the computed scenarios in comparison with the laboratory experiment.
R 2 % S SW S IG η s ¯
DR0Maxmean rel. err.+35.0%+22.5%+6.8%+3.5%+71.0%+47.3%+3.9%−0.3%
SBE2Maxmean rel. err.+20.0%+9.8%+12.9%+2.8%+40.0%+24.8%−10.4%−6.0%
SBA1Maxmean rel. err.+8.5%+4.32%−8.9%−3.2%+62.8%+30.4%−6.3%−0.01%
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Pinault, J.; Morichon, D.; Roeber, V. Estimation of Irregular Wave Runup on Intermediate and Reflective Beaches Using a Phase-Resolving Numerical Model. J. Mar. Sci. Eng. 2020, 8, 993. https://doi.org/10.3390/jmse8120993

AMA Style

Pinault J, Morichon D, Roeber V. Estimation of Irregular Wave Runup on Intermediate and Reflective Beaches Using a Phase-Resolving Numerical Model. Journal of Marine Science and Engineering. 2020; 8(12):993. https://doi.org/10.3390/jmse8120993

Chicago/Turabian Style

Pinault, Jonas, Denis Morichon, and Volker Roeber. 2020. "Estimation of Irregular Wave Runup on Intermediate and Reflective Beaches Using a Phase-Resolving Numerical Model" Journal of Marine Science and Engineering 8, no. 12: 993. https://doi.org/10.3390/jmse8120993

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop