A Model-Derived Empirical Formulation for Wave Run-Up on Naturally Sloping Beaches

Abstract

A new set of empirical formulations has been derived to predict wave run-up at naturally sloping sandy beaches. They are obtained by fitting the results of hundreds of XBeach-NH+ model simulations. The simulations are carried out over a wide range of offshore wave conditions (wave heights ranging from 1 to 12 m and periods from 6 to 16 s), and surf zone (Dean parameters a_D ranging from 0.05 to 0.30) and beach geometries (slopes ranging from 1:100 to 1:5). The empirical formulations provide estimates of wave set-up, incident and infragravity wave run-up, and total run-up R_2%. Reduction coefficients are included to account for the effects of incident wave angle and directional spreading. The formulations have been validated against the Stockdon dataset and show better skill at predicting R_2% run-up than the widely used Stockdon relationships. Unlike most existing run-up predictors, the relations presented here include the effect of the surf zone slope, which is shown to be an important parameter for predicting wave run-up. Additionally, this study shows a clear relationship between infragravity run-up and beach slope, unlike most existing predictors.

1. Introduction

Coastal zones are subject to wave and wind-induced hazards such as wave attack, flooding, and erosion. Worldwide, these coastal hazards affect millions of people living in low-lying coastal zones. With sea level rise and potentially an increase in storm and cyclone activity, this will increase substantially over the coming century.

It is therefore important to be able to reliably assess (in planning studies) and to forecast (in operational systems) the Total Water Levels (TWL) caused by wind- and wave-processes. Wind fields associated with storms exert a shear stress on the water surface which causes a water level gradient towards the coast called storm surge. This is exacerbated by the inverse barometer effect which causes an uplift of the water surface by the low-pressure fields near the eye of the storm. These wind-induced processes have been incorporated into hydrodynamic forecast and planning models such as Delft3D [1], ADCIRC [2], and SLOSH [3], which are routinely used.

The wave-driven component of the TWL consists of a mean set-up and swash run-up. The mean set-up is caused by wave breaking which exerts a wave-stress (the radiation stress [4]), which is balanced by a gradient in the water level towards the shore. The swash is the oscillatory motion of waves in the infragravity (0.004–0.04 Hz) and incident or sea-swell band (0.04–0.5 Hz) that run up and down the beach. The sum of tides, surge, mean set-up, and the crests of the oscillatory run-up motion are called the TWL. The wave contribution to the TWL may be in the order of 20% [5] on mildly sloping coasts and may be dominant on steeper coasts [6]. While the processes are largely understood and incorporated in numerical models which show considerable skill (e.g., [6,7]), process-based models to compute run-up are not routinely incorporated in forecast systems because of their computational expense. Rather, run-up formulas which are based on field observations (e.g., the Stockdon formula [8]) are typically applied [9]. The Stockdon formula [8] is the most widely used, but it is a class of empirical formulas which includes those in [10,11,12,13], among many others. A very thorough overview of existing run-up predictors and their performance is provided by Gomes da Silva [13].

While these formulas have been derived based on most of the available field observations, the overall set of observations is still small, which limits the applicability beyond the observed conditions (which have been mild mostly) and geomorphological setting (which have mostly been sandy shorelines). However, in absence of better widely applicable formulations, the empirical formulas are routinely used, which may result in an under- or overestimation of run-up levels in more extreme conditions. Both these deviations are undesired: they may result in false positives or negatives in forecasts, and in over- or under-dimensioning of coastal protection structures.

One way to extend the validity of the formulations is to actively observe run-up on a wider variety of coasts (e.g., on coral-reef lined coasts [14]) and to rapidly respond with quickly deployable field equipment to observe more extreme events, but this is costly and will increase the validity of new models only slowly in time. Physical laboratory tests of run-up can also be used to train empirical relations, although these may suffer from scaling effects. Another approach is to use numerical models in prototype scale to provide data on the basis of which more generalized formulas, or meta-models, can be derived. Stockdon [15] attempted to extend the validity of the empirical formula under more extreme conditions using the surfbeat version of XBeach [16,17], which is not capable of resolving incident swash motions. Furthermore, they found a relatively strong underestimation of modelled infragravity motions in 2D simulations.

In the present approach study, we apply the XBeach-NH+ model [18] which is a quasi-two-layer dispersive model which resolves both the incident and infragravity waves up to a relative depth kh of 4. This model has shown good skill in predicting the wave hydrodynamics in field conditions [6] and wave run-up in laboratory conditions [7]. This model, using calibrated settings, is applied to a wide range of geomorphological parameters/coastal types and for a wide range of wave forcing conditions, including obliquely incident waves and wave directional spreading. The computed run-up levels are treated as observed data and are used to fit a new empirical “VO21” formulation for the individual set-up z_mean, infragravity run-up S_IG, incident run-up S_inc, and total run-up R_2% components. Furthermore, reduction coefficients are derived to account for the effects of incident wave angle and directional spreading, using the results of a limited set of 2D XBeach-NH+ simulations.

The new formulations are validated against the Stockdon dataset [8]. The main aim of this study is not to explain the physics behind the derived relationships, but rather to define a set of empirical relations that can be used in engineering and forecasting applications. The implications for future use, the limitations of applicability, and next steps are discussed.

2. Methods

The empirical relations presented in this paper were derived from many 1D XBeach-NH+ model runs [18], in which a wide range of offshore boundary conditions and cross-shore geometries were simulated. Time series of vertical run-up were extracted from the results, and from these the mean set-up (z_mean), and the incident (S_inc) and infragravity (S_IG) components were computed. Next, empirical relations were defined that describe z_mean, S_inc, and S_IG as a function of the offshore wave height, and dimensionless surf zone and beach Irribarren numbers [19]. A relationship was then found that describes the R_2% run-up as a function of z_mean, S_inc, and S_IG. Finally, a set of reduction factors was derived to account for the effects of incident wave direction and directional spreading from a limited set of 2D simulations. The synthetic model data have not been divided here into training and testing data. The entire synthetic dataset was used for training. The Stockdon field dataset [8] was used for testing.

2.1. Model Set-Up

For the new empirical relations to be valid over a wide range of conditions, we chose a large set of offshore wave conditions (ranging from long-distance swell to locally generated storm waves) and morphologic parameters (from very mild to very steep sloping beaches). Each of the runs was carried out with varying offshore wave height (H_m0) and period (T_p), surf zone steepness (Dean parameter a_D [20]) and beach slope (tan β_b) (see

Table 1. Offshore wave and cross-shore geometry parameters.

Variable	Unit	Range
Wave height Hm0	m	1, 2, 4, 8, 12
Peak wave period Tp	s	6, 8, 12, 16
Dean slope aD	-	0.05, 0.10, 0.20, 0.30
Beach slope tan βb	-	0.01, 0.02, 0.05, 0.10, 0.20

), resulting in an initial set of 5 × 4 × 4 × 5 = 400 model runs. Unlikely (but not impossible) combinations with either very mild foreshores and very steep beaches or very steep foreshores and very mild-sloping beaches were also included in the set of simulations. Conditions, however, with physically unrealistic values for deep-water wave steepness (>0.04) were removed (e.g., H_m0 = 16 m and T_p = 6 s). This results in a final set of 280 simulations.

An example of a cross-shore profile is given in Figure 1. At the dry beach, the profile has a slope of tan β_b. The profile in the surf zone is schematized with a Dean profile [20] up to a depth of 4 times the offshore wave height. Beyond the surf zone, the bathymetry gradually drops to a depth where n (C/C_g) equals 0.8 [18]. The cross-shore grid resolution varies along the profile, such that individual waves are covered by approximately 80 grid cells. The minimum and maximum cross-shore grid spacing were set to 0.2 m and 5.0 m, respectively.

All model simulations were carried out for a 3 h period during which the offshore boundary conditions were kept constant. Deep-water waves (from Table 1) were first shoaled to the depth of the model boundary. At the boundary, the incoming waves were generated by XBeach-NH+ assuming a JONSWAP spectrum with a peak enhancement factor of 3.3.

A constant Chezy roughness (C = 70) was applied. Vertical run-up time series z_t were obtained from a run-up gauge in the model which tracks the most landward wet point. The threshold water depth for determining whether a grid cell is wet was set to 0.01 m. Tests with this parameter did not show strong sensitivity to the results in the range between 0.01 and 0.05 m. The maximum breaker steepness was set to 0.6, as recommended by [18].

In addition to the 1D profile simulations, a limited number of 2D runs were carried out to investigate the effects of the incident wave angle and directional spreading, as, by definition, the 1D profile models cannot properly take these into account. The 2D models apply the same cross-shore bathymetry profile and grid spacing. They each cover a 1 km stretch of coast with an alongshore grid spacing of 1 m, using cyclical lateral boundary conditions.

2.2. Extraction of Run-Up Results from Simulation Output

An example of XBeach-NH+ results is given in Figure 2, which shows a snapshot of cross-shore wave evolution for one of the simulations and the vertical run-up time series z_t. Note that, whereas the offshore peak wave period T_p here is 12 s, the interval between the run-up peaks in Figure 2b is typically in the order of one minute, indicating that in this simulation, the run-up signal is dominated by lower-frequency infragravity wave motions.

For each simulation, the set-up z_mean (Equation (1)) was calculated as the 3 h mean of the run-up elevation z_t.

$$z_{mean}=\overline{z_t}$$

(1)

Like Stockdon [8], we calculated the significant swash height S (Equation (2)) from the run-up elevation z_t similar to the calculation of significant wave height [21].

$$S=4\sqrt{\sum PSD\left(f\right)df }$$

(2)

where PSD is the power spectral density and f is frequency.

The significant swash height in the incident band S_inc was calculated by summing only over frequencies greater than a cut-off frequency f_c. Similarly, the swash height in the infragravity band S_IG was determined by summing over frequencies below f_c. Here, we used a cut-off frequency of half the peak frequency (Equation (3)):

$$f_c=\frac{1}{2T_p}$$

(3)

Note that, in this regard, our method differs from [8] who use a constant cut-off frequency f_c of 0.05 Hz.

The R_2% run-up height (defined as the value exceeded by two percent of the run-up events) was determined by taking the 98th percentile of the peaks in the vertical run-up time series.

The derivation of the empirical relation for S_IG was carried out in two steps (as shown in Section 3.3) and first requires an estimate of the incoming infragravity wave height H_in,IG and the infragravity period T_m-10,IG close to the shoreline. These were determined with the method presented in [21].

The aim of this study was to define z_mean, S_inc, and S_IG as the product of the deep-water wave height and some function of the dimensionless surf zone and beach Irribarren numbers, i.e., $H_{m0}f\left({\xi }_s,{\xi }_b\right)$. The procedure for determining the Irribarren numbers is as follows.

For each simulation, we first estimated the breaking wave depth with Equation (4):

$$h_{br}=\frac{H_{m0}}{{\gamma }_{br}}$$

(4)

In this study, we assumed a constant wave breaking parameter γ_br = 1.0. From there, the width W_s and mean slope of the surf zone tan β_sz are determined with Equations (5) and (6):

$$W_s={\left(\frac{H_{m0}}{{\gamma }_{br}{a}_D}\right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$

(5)

$$\tan {\beta }_{sz}=\frac{h_{br}}{W_s}$$

(6)

The slope of the beach tan β_b is one of the input variables and does not require calculation. The dimensionless surf similarity numbers for the surf zone ξ_sz and beach ξ_b are given by Equations (7) and (8).

$${\xi }_{sz}=\frac{\tan {\beta }_{sz}}{\sqrt{\frac{H_{m0}}{L_0}}}$$

(7)

$${\xi }_b=\frac{\tan {\beta }_b}{\sqrt{\frac{H_{m0}}{L_0}}}$$

(8)

Here, L₀ is the deep-water wave length (Equation (9)):

$$L_0=\frac{g{T}_p^2}{2\pi }$$

(9)

3. Results

3.1. Set-Up z_mean

Qualitative analysis of the computed set-up z_mean reveals a clear pattern (Figure 3) where the relative dimensionless set-up z_mean/H_m0 always increases with increasing ξ_sz. The relation between z_mean/H_m0 and ξ_b is more complicated. It shows that under dissipative beach conditions (i.e., low values of ξ_b), the relative set-up first increases with increasing ξ_b, and then, beyond a certain tipping point where beach conditions become reflective, starts to decrease. The model results indicate that the location of this tipping point is not fixed but depends on the surf zone conditions. Under dissipative conditions in the surf zone (ξ_sz < 0.1) the tipping point is located around ξ_b ≈ 0.05. For more reflective surf zone conditions (ξ_sz > 1.0), it is located around ξ_b ≈ 0.6.

It was experimentally found that this behavior can be described in the following form:

$$\begin{array}{c}{z}_{mean}=H_{m0}\left(A+B{\xi }_b^C{e}^{-D{\xi }_{sz}^E\sqrt{{\xi }_b}}\right)\end{array}$$

(10)

where A, B, C, D, and E are constants yet to be determined. After applying a non-linear least-square fitting procedure to determine their optimal values, the formulation for z_mean is given by Equation (11):

$$z_{mean}=H_{m0}\left(0.099+3.05{\xi }_b^{0.66}{e}^{-1.77{\xi }_{sz}^{-0.36}\sqrt{{\xi }_b}}\right)$$

(11)

The comparison between simulated set-up and Equation (11) (Figure 4) reveals very good agreement between the empirical fit and the XBeach-NH+ results.

3.2. Incident Run-Up S_inc

The relation between ξ_sz, ξ_b, and the relative dimensionless incident run-up S_inc/H_m0 is somewhat similar to that of the wave set-up (see Figure 5) described in the previous section. S_inc/H_m0 also increases with increasing ξ_sz values, and under dissipative beach conditions (i.e., low values of ξ_b), S_inc/H_m0 increases with increasing ξ_b, albeit much more rapidly than the relative wave set-up. Beyond the tipping point, S_inc/H_m0 remains fairly constant. Under very reflective conditions (ξ_b > ~3 and ξ_sz > ~1), the incident run-up approaches the standing wave limit (S_inc/H_m0 = 2).

It was again experimentally found that a good prediction of S_inc can be obtained when written in the following form:

$$S_{inc}=A{H}_{m0}{\xi }_b^2\tanh \left(\frac{B{\xi }_{sz}^C}{{\xi }_b^2}\right)$$

(12)

The non-linear fitting procedure is now used to determine optimal values for the constants A, B, and C, resulting in the following empirical relation for S_inc:

$$S_{inc}=0.986H_{m0}{\xi }_b^2\tanh \left(\frac{1.79{\xi }_{sz}^{0.62}}{{\xi }_b^2}\right)$$

(13)

The comparison between simulated set-up and Equation (13) (Figure 6) reveals generally very good agreement between the empirical fit for S_inc and XBeach-NH+ results.

3.3. Infragravity Run-Up S_IG

Inspection of the model results reveals that the frequency of infragravity waves at the shoreline strongly depends on the surf zone slope and, to a much lesser extent, the offshore wave steepness. For relatively steep surf zones (tan β_s > 0.04), the T_m-10,in,IG wave period of the infragravity part of the power density spectrum is typically in the order of 5–7 times the deep-water peak wave period T_p. At milder sloping surf zones however, the model shows a shift of infragravity wave energy towards lower frequencies (Figure 7b). We believe that this shift is related to infragravity wave–wave interactions, but the exact physical process behind it remains to be investigated.

It is expected that the infragravity wave period at the shoreline is an important parameter for predicting the S_IG run-up at the beach. A different strategy is therefore chosen to derive the formula for S_IG in which this assumption is reflected. In the first step, we derive a relation for the incoming infragravity wave height H_in,IG at a location close to the shoreline, as well as the corresponding infragravity wave period T_m-10,in,IG. In the second step, we attempt to derive a relation for S_IG, based on H_in,IG, T_m-10,in,IG, and the beach slope.

The nearshore location in the first step was chosen such that the bulk of the incident wave energy has dissipated, yet it is deep enough that it does not become dry at any time (and therefore yields a continuous time series of water levels and velocities). The latter is a requirement for accurately determining H_in,IG and T_m-10,in,IG. Here, we chose a point at which the still water depth equals one third of the offshore wave height H_m0. For a Dean profile and γ_br = 1.0, the distance from the shore to this point equals 19.2 percent of the surf zone width.

After some experimentation, the following fits are found to estimate the incoming infragravity wave height (Equation (14)) and period (Equation (15)).

$$\frac{H_{in,IG}}{H_{m0}}=2.29\sqrt{{\beta }_{sz}}{e}^{-17.5{\beta }_{sz}^2}+0.183e^{-2798{\left(\frac{H_{m0}}{L_0}\right)}^2}$$

(14)

$$\frac{T_{m-1,0,in,IG}}{T_p}=18.1\frac{{\left(\frac{H_{m0}}{L_0}\right)}^{0.07}}{{\beta }_s^{0.43}}+20.9{\beta }_s$$

(15)

The comparison between simulation results and Equations (14) and (15) (Figure 8) reveals generally good agreement.

Rather than estimating the infragravity run-up with the original beach Irribarren number ξ_b which includes the deep-water steepness, we use an Irribarren number ξ_b,IG that is computed with the infragravity wave steepness near the shore (Equation (16)):

$${\xi }_{b,IG}=\frac{{\beta }_b}{\sqrt{\frac{H_{in,IG}}{L_{in,IG}}}}$$

(16)

where the infragravity wave length L_in,IG is computed with Equation (17).

$$L_{in,IG}=T_{m-1,0,in,IG}\sqrt{\frac{1}{3}{H}_{m0}g}$$

(17)

The response of the relative dimensionless infragravity run-up S_IG/H_in,IG shows a clear distinction between dissipative and reflective beach conditions (Figure 9). Under dissipative beach conditions (i.e., low values of ξ_b,IG), the relative run-up first increases with increasing ξ_b,IG, and then, beyond a tipping point ${\xi }_{b,IG}^m$, starts to linearly decrease with increasing ξ_b,IG. We also see that S_IG cannot be properly described as a function of H_in,IG and ξ_b,IG alone. For more reflective conditions on the beach, a dependency on the surf zone slope β_sz is found which affects the boundary ${\xi }_{b,IG}^m$ between the dissipative and reflective regimes. This dependency can likely be attributed to the fact that in steeper surf zones, the incident wave energy is relatively higher at the location where the incoming infragravity wave height is determined (i.e., at a depth of one third the offshore wave height). Shoreward of this location, a portion of the remaining short-wave energy is still converted to infragravity energy, thereby increasing the infragravity run-up at the beach. Under very reflective conditions (ξ_b,IG > ~4), the infragravity run-up approaches the standing wave limit (S_IG/H_in,IG = 2).

The behavior seen in Figure 9 can be reasonably well described when written in the following form:

$$\frac{S_{IG}}{H_{in,IG}}=A {\xi }_{b,IG}^B \mathrm{for} {\xi }_{b,IG}\le {\xi }_{b,IG}^m$$

(18)

$$\frac{S_{IG}}{H_{in,IG}}=max\left(A {\xi }_{b,IG}^m{}^B-C\sqrt{tan {\beta }_{sz}}\left({\xi }_{b,IG}-{\xi }_{b,IG}^m\right), 2.0\right) \mathrm{for} {\xi }_{b,IG}>{\xi }_{b,IG}^m$$

(19)

$${\xi }_{b,IG}^m=D\tan {\beta }_{sz}{}^E$$

(20)

After determining the optimal values for the constants A, B, C, D, and E, with a non-linear fitting procedure, the full set of equations is given by Equations (21)–(23).

$$\frac{S_{IG}}{H_{in,IG}}=3.48 {\xi }_{b,IG}^{0.59} \mathrm{for} {\xi }_{b,IG}\le {\xi }_{b,IG}^m$$

(21)

$$\frac{S_{IG}}{H_{in,IG}}=max\left(3.48 {\xi }_{b,IG}^m{}^{0.59}-0.512\sqrt{tan {\beta }_{sz}}\left({\xi }_{b,IG}-{\xi }_{b,IG}^m\right), 2.0\right) \mathrm{for} {\xi }_{b,IG}>{\xi }_{b,IG}^m$$

(22)

where:

$${\xi }_{b,IG}^m=2.32\tan {\beta }_{sz}{}^{0.24}$$

(23)

The comparison of S_IG obtained from XBeach-NH+ versus the outcome of Equations (14)–(23) reveals good agreement between the numerical simulation and the empirical formulations (Figure 10).

3.4. Run-Up R_2%

With satisfactory estimates for z_mean, S_IG, and S_inc, it is now possible to derive a formulation for the total run-up R_2%. Initially, we attempted to use a formulation similar to the one used by Stockdon [8]:

$$R_{2\%}=z_{mean}+A\sqrt{B{S}_{IG}^2+C{S}_{inc}^2}$$

(24)

Here, R_2% is the run-up computed by the XBeach-NH+ model and z_mean, S_IG, and S_inc are computed with Equations (11), (13) and (14)–(23), respectively. The optimal values for the constants A, B, and C can be determined using a non-linear fitting procedure.

It was however found that the inclusion of ξ_sz and ξ_b in the estimation of R_2% yielded a significant performance improvement. Eventually, the best fit that was found for R_2% is given by Equation (25).

$$R_{2\%}=z_{mean}+0.834\sqrt{\left(0.873S_{IG}^2+0.725S_{inc}^2\right)} {\xi }_{sz}^{-0.094}{\xi }_b^{0.140}$$

(25)

There is generally good agreement between run-up values from the XBeach-NH+ simulations and those computed with Equation (25) (Figure 11).

3.5. Angle of Incidence

A limited set of 2D simulations were carried out to determine the reduction factors for the offshore wave angle φ with respect to the coast. Due to the strongly increased computational effort required to run the 2D simulations, the 2D runs were only performed with on offshore wave height H_m0 of 2.0 m. Furthermore, the beach slope tan β_b was kept constant at 0.05. Model input parameters were kept identical to those in the 1D simulations. The 2D model was forced at the boundaries with waves coming in at angles φ of 0, 20, 40, and 60 degrees. The alongshore extent and resolution of the 2D model was 1 km and 1 m, respectively. The reduction factors were determined with:

$$f_{\phi }=\frac{S_{2D}}{S_{2D, 0}}$$

(26)

where S_2D and S_2D,0 are z_mean, S_inc, or S_IG from the 2D simulations with oblique (φ > 0), and normal incident wave (φ = 0), respectively. The reduction factors are to be applied to Equations (11), (13), (21) and (22) before computing the total run-up with Equation (25).

The incident wave run-up height S_inc was found to be insensitive to the incident wave angle. We therefore used f_φ,inc = 1.0. The set-up z_mean and infragravity run-up height S_IG, however, showed a clear decrease with increasing wave angle φ (Figure 12). The reduction effect of the wave angle on z_mean was stronger for steeper offshore waves (Figure 12a), whereas the reduction of the infragravity run-up S_IG was found to be greater for relatively steep surf zones (Figure 12b).

A good match for f_φ,zmean and f_φ,inc (see Figure 13) was found when the reduction factors for set-up and S_IG are estimated with Equations (27) and (28):

$$f_{\phi ,z_{mean}}=1-0.0035{\left(\frac{H_{m0}}{L_0}\right)}^{0.32}{\phi }^{1.43}$$

(27)

$$f_{\phi ,IG}={\left(\cos \phi \right)}^{0.40+2.7\tan {\beta }_{sz}}$$

(28)

3.6. Directional Spreading

A similar approach was used to determine reduction factors for directional spreading D (here defined as the directional standard deviation in degrees), but rather than varying the angle of incidence, we varied the directional spreading D in the offshore wave spectra that make up the offshore boundary conditions, using values of 0, 15, and 30 degrees. Similar to the angle of incidence, the incident wave run-up height S_inc was found to be insensitive to the directional spreading. We therefore used f_D,inc = 1.0. The set-up and infragravity run-up, however, both decrease with increasing directional spreading (Figure 14). The reduction effect for both z_mean and S_IG was enhanced under more dissipative conditions in the surf zone (lower values of ξ_sz). Qualitatively similar results were obtained by [22].

The largest effect of directional spreading was found on the infragravity run-up, where S_IG can be reduced by over 40 percent. The best fit is found when the reduction factors are estimated with Equations (29) and (30) (see also Figure 15).

$$f_{D,set-up}=e^{-\left(0.031D^{0.69}(1-\tanh \left(0.67{\xi }_s\right)\right)}$$

(29)

$$f_{D,IG}=e^{-\left(0.048D^{0.67}(1-\tanh \left(0.51{\xi }_s\right)\right)}$$

(30)

3.7. Validation

Until this point, the data to derive the empirical relations have all been obtained from hydrodynamic model results. Here, we now validate the results against field observations, using the video-derived run-up data collected by Stockdon [8], hereafter referred to as S06. As the original S06 data do not include information on offshore wave angle and directional spreading, these were obtained from WAVEWATCH III re-analysis data by NOAA [23]. Surf zone profiles were obtained by digitizing the bathymetry profiles provided by S06. Tide levels at each data point were computed from the TPXO 7.2 database [24].

To compute the run-up, the time-varying surf zone slope was computed separately for each data point in the S06 dataset. This was conducted with Equations (4)–(6), using the surf zone depth profile and including the tide level. The beach steepness β_b was taken from the S06 dataset.

The scatter plot in Figure 16a shows the comparison between observed and computed R_2% values (using Equation (25)). For reference, the results of the S06 formulations are shown in Figure 16b. The results with the empirical relations from the present study (VO21) generally compare well with the observations. The VO21 error statistics are significantly better than those from S06. In particular, the VO21 formulations do not suffer from the negative bias seen in S06. Furthermore, the R² correlation coefficient of VO21 (0.772) is higher than that of SO6 (0.715), and the root-mean-square error is lower (0.304 vs. 0.372 m). It should be noted here that there is a considerable amount of scatter in the S06 run-up data. The scatter can to some extent be explained by the relatively short 17 min video observation intervals which introduce uncertainty in the run-up estimates. The uncertainty in infragravity run-up estimates may be as high as 65% [25].

The performance of the VO21 results and comparison to S06 is further discussed in the next section.

4. Discussion and Conclusions

The empirical relations described in Section 3.1, Section 3.2, Section 3.3 and Section 3.4 are able to accurately reproduce the run-up characteristics computed by the XBeach-NH+ 1D profile models. Furthermore, when combined with the reduction factors for incident wave angle (Section 3.5) and directional spreading (Section 3.6), the use of these new relations results in improved error statistics compared to the S06 relations, using the same dataset. In particular, the new relations do not suffer from the negative bias that is seen in the S06 relations. On the one hand, this result is not entirely surprising, as the method presented here takes in more forcing data (tides, surf zone slope, incident angle, and directional spreading) than was used in the S06 formulations. On the other hand, it must be noted that the new VO21 relations were derived from an entirely different dataset obtained from numerical hydrodynamic model results, whereas the S06 formulations are the result of a direct fit of the field observation data. Furthermore, the model results were obtained from a wider range of wave and topographic conditions (see

Table 2. Range of hydrodynamic and topographic conditions in the VO21 dataset and the Stockdon field data [8].

	VO21 (-)	Stockdon [8]
Wave height Hm0	1.0–12.0	0.4–4.1
Wave period Tp	6.0–16.0	4.0–17.0
Surf zone slope βs	0.007–0.097	0.007–0.083
Beach slope βb	0.01–0.20	0.01–0.16

), with the exception of somewhat lower minimum wave heights and periods in the field data. It is therefore to be expected that the VO21 formulations are valid over a wider range of conditions than S06, and in particular, under extreme wave events which are not included in the field data.

It should be noted here that many alternative run-up formulations have been derived in recent years. For a thorough overview, we refer to [13]. A number of these formulas have been parameterized with larger and/or different field datasets than S06, and some have shown better performance than S06. The main reason for quantitively comparing VO21 only to the S06 relations lies in the fact that S06 is the only formulation that has been trained solely with our validation dataset (i.e., the video-derived field data collected by [8]). Furthermore, S06 remains at the time of writing the most known and applied formula in the engineering community.

Whereas most existing predictors in [13] include only one slope (typically the swash zone slope), VO21 considers both the surf zone slope (which varies with the tide) and the swash slope. Another important feature of VO21 is that it distinguishes between dissipative and reflective conditions in both the surf zone and on the beach. The fact that these effects are not included in many other predictors may explain why previous studies did not see an improvement of their infragravity run-up predictions by including beach slope (e.g., [8,12,22]). In the present study, there is a clear response of the infragravity run-up height to changes in the beach slope. We see a positive relation between Irribarren numbers and both set-up and wave run-up under dissipative conditions at the beach. Beyond a certain tipping point, however, set-up and infragravity run-up start to decrease with increasing Irribarren numbers, and incident run-up remains constant. Other parameters that have been shown to be important in this study, and which are typically not included in existing relations, are the incident wave angle and directional spreading of the offshore wave conditions.

The effect of excluding the above-mentioned parameters in our validation against the field data have been investigated. The results of this sensitivity test are presented in

Table 3. Effects of omitting input parameters on error statistics.

	R2 (-)	RMSE (m)	Bias (m)
Full	0.772	0.304	0.006
No tide/surge	0.756	0.315	−0.071
φ = 0	0.755	0.323	0.101
D = 0	0.750	0.508	0.334
tan βs = 0.046	0.706	0.387	0.047

. Exclusion of the vertical tide results in slightly poorer performance (with slightly decreased R² and increased RMSE). Not including offshore wave angle (φ = 0) and particularly the directional spreading (D = 0) leads to slightly lower correlation R² with the data and introduces a strong positive bias. The greatest effect is however seen when excluding the effect of the surf zone slope. This was investigated by maintaining a constant surf zone slope of 0.046 (the mean of all values computed for the S06 dataset). Not including the effect of the surf zone slope significantly reduces the R², and increases the RMSE by 27 percent, again highlighting the relative importance of this parameter.

In Figure 17, R_2% predictions from VO21 are shown for mild (18a), intermediate (18b), and steep (18c) beach and mild (red), intermediate (blue) and steep (green) underwater slopes. The S06 formula results, which do not include the surf zone slope, are shown in black. The figure highlights the importance of the effect of the underwater steepness on run-up (compare the red, blue, and green lines). Furthermore, it suggests that S06 yields significantly lower run-up during large wave events (compared to VO21) at intermediate and steep underwater slopes, except for the case of very steep beach slopes.

The fact that VO21 requires more input data can of course also be a disadvantage, in particular when such data are not available. Here, we list recommendations when some of the data are missing. Typically, tide data can always be obtained, either from nearby tide gauges, or from global tide models, such as [24]. The incident wave angle and directional spreading can be obtained from nearby directional wave buoys, or from global wave re-analyses or forecasts. Setting either parameter to 0.0 typically leads to an overestimation of the run-up and is therefore not recommended. Tests with the validation dataset reveal a zero bias between VO21 predictions and observations when incident wave direction and spreading are set to 29 degrees and 25 degrees, respectively. When the surf zone profile is not known, one can estimate a Dean slope parameter a_D by, e.g., estimating the distance L₁₀ between the shoreline and the 10 m depth contour from a global bathymetry dataset, such as [26]. The Dean parameter a_D can then be computed with $a_D=10.0/L_{10}^{\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}$. From there, the surf zone slope can be computed with Equations (4)–(7).

The VO21 formulations are not valid along reef-fringed coastlines. Nor are they likely valid on gravel beaches where permeability may have a significant effect on run-up elevations [27].

Some, potentially important, parameters have not been taken into account in VO21. These include the effects on vertical run-up of bed roughness, the presence of breaker bars, and alongshore non-uniformities (e.g., rip channels). Furthermore, the effect of the offshore wave frequency spectrum has not been investigated. Guza and Feddersen [22] showed that R_2% wave run-up tends to increase with increasing frequency spread. Here, however, we used a JONSWAP spectrum with a constant peak enhancement factor of 3.3. Further work is required to include the above-mentioned effects.

It is recommended that the VO21 relations are further validated against more extensive datasets, in particular datasets with more extreme conditions. Although its predictions now closely match the numerical model results under extreme conditions (H_m0 > 4 m), we have not performed a validation of VO21 against field data under such conditions, as these are not available in the validation dataset [8].

Finally, the present study shows a relation between the surf zone slope on the infragravity wave period at the shoreline. For relatively steep surf zones, the infragravity wave period T_m-10,IG is typically in the order of 5–7 times the deep-water peak wave period T_p. At milder sloping surf zones, however, this study shows a shift of infragravity wave energy towards lower frequencies. The cause of this phenomenon remains to be further investigated.