Estimating nearshore infragravity wave conditions at large spatial scales

Leijnse, Tim W. B.; van Ormondt, Maarten; van Dongeren, Ap; Aerts, Jeroen C. J. H.; Muis, Sanne

doi:10.3389/fmars.2024.1355095

ORIGINAL RESEARCH article

Front. Mar. Sci., 13 March 2024

Sec. Coastal Ocean Processes

Volume 11 - 2024 | https://doi.org/10.3389/fmars.2024.1355095

Estimating nearshore infragravity wave conditions at large spatial scales

Tim W. B. Leijnse^1,2*

Maarten van Ormondt³

Ap van Dongeren^2,4

Jeroen C. J. H. Aerts^1,2

Sanne Muis^1,2

¹Department of Water and Climate Risk, Institute for Environmental Studies (IVM), Vrije Universiteit Amsterdam, Amsterdam, Netherlands
²Marine and Coastal Systems, Deltares, Delft, Netherlands
³Deltares USA, Silver Spring, MD, United States
⁴Coastal and Urban Risk and Resilience, Water Science and Engineering Department, IHE Delft Institute for Water Education, Delft, Netherlands

Infragravity waves may contribute significantly to coastal flooding, especially during storm conditions. However, in many national and continental to global assessments of coastal flood risk, their contribution is not accounted for, mostly because of the high computational expense of traditional wave-resolving numerical models. In this study, we present an efficient stationary wave energy solver to estimate the evolution of incident and infragravity waves from offshore to the nearshore for large spatial scales. This solver can be subsequently used to provide nearshore wave boundary conditions for overland flood models. The new wave solver builds upon the stationary wave energy balance for incident wave energy and extends it to include the infragravity wave energy balance. To describe the energy transfer from incident to infragravity waves, an infragravity wave source term is introduced. This term acts as a sink term for incident waves and as a complementary source term for infragravity waves. The source term is simplified using a parameterized infragravity wave shoaling parameter. An empirical relation is derived using observed values of the shoaling parameter from a synthetic dataset of XBeach simulations, covering a wide range of wave conditions and beach profiles. The wave shoaling parameter is related to the local bed slope and relative wave height. As validation, we show for a range of cases from synthetic beach profiles to laboratory tests that infragravity wave transformation can be estimated using this wave solver with reasonable to good accuracy. Additionally, the validity in real-world conditions is verified successfully for DELILAH field case observations at Duck, NC, USA. We demonstrate the wave solver for a large-scale application of the full Outer Banks coastline in the US, covering 450 km of coastline, from deep water up to the coast. For this model, consisting of 4.5 million grid cells, the wave solver can estimate the stationary incident and infragravity wave field in a matter of seconds for the entire domain on a regular laptop PC. This computational efficiency cannot be provided by existing process-based wave-resolving models. Using the presented method, infragravity wave-driven flooding can be incorporated into large-scale coastal compound flood models and risk assessments.

1 Introduction

Hundreds of millions of people live in low-elevation coastal zones (McMichael et al., 2020), of which there are more than 189 million in the 1-in-100-year floodplain (Neumann et al., 2015). The majority lives near deltas and estuaries, which are prone to compound flooding, where multiple coastal flood hazards such as storm surges from tropical and extra-tropical storms, spring tides, waves, high river discharges, and heavy rainfall events might occur at once with exacerbated flooding as a result (Mousavi et al., 2011; Wahl et al., 2015; Ward et al., 2018). While the number of people in the coastal zone will increase in the future due to socioeconomic trends (e.g., McGranahan et al., 2007; Jones and O’Neill, 2016; Merkens et al., 2016), the frequency and magnitude of the coastal hazards will also increase due to sea-level rise and possibly more intense storms due to climate change (e.g., Emanuel, 2013; Stocker, 2014; Tebaldi et al., 2021). It is, therefore, important to accurately assess flood risk to support the design of adequate flood protection measures (e.g., Diermanse et al., 2023) and flood early warning systems (e.g., Winter et al., 2020; De Kleermaeker et al., 2022).

In some coastal areas, waves can be the dominant driver of extreme water levels (e.g., Parker et al., 2023). Especially along steeper coasts without shelves, waves can be the dominant driver of coastal flooding. For instance, during Typhoon Haiyan (2013) in the Philippines, the village of Hernani experienced a storm surge of less than a meter. However, strong infragravity (IG) waves destroyed the entire village despite the relatively low surge and the fact that the village is sheltered behind a coral reef flat (Roeber and Bricker, 2015). IG waves are long waves forced by sea-swell wave groups, with periods between 25 and 250 s (Longuet-Higgins and Stewart, 1964; Bertin et al., 2018). In deep water, IG waves have amplitudes of mere centimeters (Bertin et al., 2018). Shoreward, as the incident waves shoal on the sloping bed, the bound IG wave is no longer in equilibrium with the incident wave groups, and a phase shift allows energy transfer from the incident wave band to the IG wave band (e.g., Janssen et al., 2003). This causes the IG wave to grow in amplitude depending on the normalized bed slope (Battjes et al., 2004), up to wave heights exceeding 1 m (Fiedler et al., 2015; Bertin et al., 2018). In the surf zone, the wave groups are destroyed and IG waves are released and propagate toward the shore as free waves (e.g., Masselink, 1995; Janssen et al., 2003). Nearshore, the IG waves can be dominant over the remaining incident wave signal (e.g., Guza and Thornton, 1982) and cause runup and flooding. Despite their importance, coastal compound flood risk modeling assessments and early warning systems do generally not include these IG wave-driven contributions, except at local (i.e., city) scales (e.g., inclusion of wave overtopping; Qiang et al., 2021).

On the local scale, wave-driven flooding can be modeled using advanced wave-resolving models such as NLSWE models (non-linear shallow water equation; e.g., SWASH; Zijlema et al., 2011, XBeach; Roelvink et al., 2009), Boussinesq models (e.g., BOSZ; Roeber et al., 2010), or RANS models (Reynolds-averaged Navier–Stokes; e.g., OpenFOAM^®). However, when modeling large coastal scales (>100 km), such models are too computationally expensive. Flood risk assessments and early warning systems increasingly cover national to continental scales (e.g., Leijnse et al., 2023a; Nederhoff et al., 2024b), but the current common practice is to either neglect wave-driven components (e.g., Muis et al., 2016), use assumptions such as a static wave setup of 20% of the offshore wave height (US Army Corps Of Engineers, 2002; Vousdoukas et al., 2016; Camus et al., 2021; Van Oosterhout et al., 2023), or use empirical (e.g., Stockdon et al., 2006) wave runup formulas (e.g., as in Kirezci et al., 2020).

Sufficiently accurate results in modeling compound flooding can be achieved using faster reduced complexity models such as LISFLOOD-FP (Bates et al., 2010) or SFINCS (Super-Fast INundation of CoastS) (Leijnse et al., 2021). SFINCS can be used at the scale of 100–1,000 km of coastlines (e.g., Eilander et al., 2022; Grimley et al., 2023; Leijnse et al., 2023b). For Typhoon Haiyan in the Philippines, Leijnse et al. (2021) have shown that when SFINCS is forced with water-level boundary conditions that include waves, flooding can be modeled with accuracy similar to a wave-resolving model (Roeber and Bricker, 2015). However, this approach still required accurate boundary conditions obtained from an expensive time-dependent two-dimensional (2D) XBeach two-layer non-hydrostatic model (de Ridder et al., 2021), making the method too expensive to be applied at large spatial scales (e.g., Gaido Lasserre et al., 2020; Nederhoff et al., 2024a).

Wave spectral models (e.g., Booij et al., 1999; Tolman, 2002) are also not yet suited to provide these nearshore conditions as they lack the physics to solve IG wave energy. Reniers and Zijlema (2022) recently extended a 2D spectral wave model to include bound IG waves. They showed that their method can provide improved IG boundary conditions in relatively shallow water, which can subsequently be used by other models to compute wave runup and overtopping. The applicability of their method is however limited to sandy coasts with mild slopes and gently varying alongshore bathymetry. Also, the transfer of wave energy to the IG waves is one way, and there is no sink term for the incident wave energy based on the energy transfer. Additionally, these 2D spectral wave models are too computationally expensive to be applied at high resolutions and large spatial scales.

On the other hand, one-dimensional (1D) derived methods can provide nearshore estimates very quickly, either based on process-based models (e.g., Pearson et al., 2017; van Ormondt et al., 2021; Athanasiou et al., 2022) and/or as data-driven methods (e.g., Stockdon et al., 2006; Dalinghaus et al., 2022). Although these are very fast, they are not applicable to complex coastlines. Also, often no data along the entire profile can be generated nor can they be forced by multidirectional offshore sea states and, therefore, do not suffice.

This study aims to overcome these limitations by presenting a new efficient wave solver that estimates the transformation of IG energy into the nearshore and can efficiently provide wave boundary conditions to drive flood models.

The new solver builds upon the stationary wave energy balance solver as used in Reyns et al. (2023), which simulates the transformation of incident wave energy. The formulations in this model are extended to also include the IG wave energy balance. For this purpose, we introduce an IG wave source term that describes the energy transfer from incident to IG waves. This term acts as a sink term for the incident waves and as a complementary source term for the IG waves. The IG wave source term is determined based on the steady-state wave energy balance of Schäffer (1993) and includes an IG wave shoaling parameter. The observed values of this parameter are obtained from the dataset of van Ormondt et al. (2021), which consists of XBeach simulation results over a wide range of wave conditions and beach profiles. These are used to derive an empirical relation that describes the shoaling parameter as a function of the local wave and bathymetrical conditions. With this parameterization, the IG wave source term can be included in the wave energy balance solver to compute the incident and IG wave transformation in the nearshore with offshore incident wave heights and local bathymetry as input. We demonstrate that, with the introduced wave solver, one can estimate incident and IG wave conditions in a computationally efficient way at large spatial scales.

2 Materials and methods

2.1 Methodology

The methodology starts with the formulations for the 2D wave solver (§2.1.1). In the next step, the 1D steady-state cross-shore IG wave energy balance of Schäffer (1993) is used to derive a description of the IG wave source term with a simple shoaling parameter $α_{i g}$ (§2.1.2). Then, a relation between this parameter and the local bed slope parameter $β$ and incident wave height over depth ratio $γγ$ is derived using 70% of the van Ormondt et al. (2021) (VO21) dataset as training profiles (§2.1.3). The introduced IG wave source term with parameterized shoaling parameter is validated for the unseen 30% validation profiles (§3.1). The IG wave transformation is verified in 1D for the Boers and GLOBEX laboratory flume tests (§3.2). Then, the performance in 2D for a real-world application is verified for the DELILAH field case (§3.3). Finally, the model is applied at the Outer Banks (North Carolina, USA) using ERA5 input data at a large spatial scale. This case shows how the model can be used to efficiently estimate large-scale nearshore boundary conditions of IG wave heights (§3.4).

2.1.1 Formulation of the 2D wave solver

The proposed wave solver is based on the 2D implicit stationary wave energy balance as described in Reyns et al. (2023). The model solves an incident band wave energy balance (Equation 1) using an implicit first-order upwind scheme with quadrant sweeping to ensure convergence of the solution. The model accounts for incident band wave-breaking dissipation and solves multiple directional bins and one frequency bin. Rather than solving the simplified 2D stationary incident wave energy balance over distances x and y as in, e.g., XBeach (Roelvink et al., 2009), it considers the distance x along the wave direction $θ$ . In the present formulation, an IG wave sink term $S_{i g}$ is added which describes the energy transfer from the incident to the IG waves:

\begin{array}{l} \frac{\partial c_{g} E_{i n c}}{\partial x} + \frac{\partial c_{θ} E_{i n c}}{\partial θ} + D_{w, i n c} + D_{b o t, i n c} - S_{i g} = 0 & (1) \end{array}

Here $E_{i n c}$ is the spectral energy density of the incident waves, $c_{g}$ the group velocity, $c_{θ}$ the refraction velocity, $θ$ the wave direction, $D_{b o t, i n c}$ the incident wave bottom friction as in XBeach (Roelvink et al., 2009), and $D_{w, i n c}$ the incident wave dissipation density following Baldock et al. (1998). For details, we refer to Reyns et al. (2023). The boundary conditions for Equation 1 are sea states prescribed at the offshore boundary using spectral parameters as the incident significant wave height $H_{s}$ , peak wave period $T_{p}$ , mean wave direction, and the wave spreading factor s.

The IG wave energy balance (Equation 2) is solved concurrently with Equation 1 and reads:

\begin{array}{l} \frac{\partial c_{g} E_{i g}}{\partial x} + \frac{\partial c_{θ} E_{i g}}{\partial θ} + D_{w, i g} + D_{b o t, i g} + S_{i g} = 0 & (2) \end{array}

where $E_{i g}$ is the spectral energy density of the IG waves. Dissipation of the IG wave energy $D_{w, i g}$ is included using the same wave dissipation formulation of Baldock et al. (1998), to be consistent with the incident wave energy balance. Also, the IG wave bottom friction $D_{b o t, i g}$ is included following Henderson and Bowen (2002) (as used previously in, e.g., van van Dongeren et al., 2007).

Together, Equations 1 and 2 describe our 2D wave solver to efficiently estimate incident and IG wave energy from offshore until nearshore. Hereby, a description to estimate the IG wave source term $S_{i g}$ efficiently is still needed, which is derived in Sections 2.1.2 and 2.1.3. Offshore IG wave boundary conditions are estimated using Herbers et al. (1994). A more comprehensive description of this method is given in Appendix B.

2.1.2 Infragravity wave source term

The starting point for deriving an IG wave source term $S_{i g}$ is the 1D steady-state wave energy balance for IG waves (i.e., averaged over the IG wave periods) in cross-shore direction, following Phillips (1977) and Schäffer (1993):

\begin{array}{l} \frac{\partial c_{g} E_{i g}}{\partial x} + U \frac{\partial S_{x x}}{\partial x} + D_{b o t, i g} = 0 & (3) \end{array}

which states that the IG energy flux (first term) changes due to the work (second term) that the incident waves (represented by the cross-shore radiation stress $S_{x x}$ ; Dean and Dalrymple, 1991) do on the IG waves (represented by its velocity amplitude $U$ ), with a sink term due to bottom friction ( $D_{b o t, i g}$ , third term). We only consider shoreward propagating IG wave energy, as we are interested in obtaining wave boundary conditions for nearshore wave transformation and flooding models.

The IG wave averaged work term can be written as the product of the velocity amplitude $U$ , the magnitude of the radiation stress gradient and an unknown factor (defined here as the shoaling parameter $α_{i g}$ ), which depends on the unknown phase between the velocity and the radiation stress. The velocity amplitude $U$ can be written in terms of wave height as in Svendsen (1984):

\begin{array}{l} U = \frac{H_{i g}}{2} \frac{c_{g}}{h} ~ \sqrt{E_{i g}} \frac{c_{g}}{h} & (4) \end{array}

where $H_{i g}$ is the incoming IG wave height, which is proportional to the square root of the IG wave energy $E_{i g}$ ; $c_{g}$ is the wave group speed; and $h$ is the water depth. The cross-shore radiation stress $S_{x x}$ can be written as (Dean and Dalrymple, 1991):

\begin{array}{l} S_{x x} = (2 n - 1 / 2) E_{i n c} & (5) \end{array}

With $n = c_{g} / c$ and $E_{i n c}$ is the incident wave energy. Substituting Equations 4 and 5 in Equation 3, we then obtain a description of the IG wave energy balance that depends on the known parameters $c_{g}$ , $n$ , $h$ , and $E_{i n c}$ and a yet unknown shoaling parameter $α_{i g}$ , to describe the IG wave energy $E_{i g}$ :

\begin{array}{l} \frac{\partial c_{g} E_{i g}}{\partial x} + α_{i g} \sqrt{E_{i g}} \frac{c_{g}}{h} \frac{\partial (2 n - 1 / 2) E_{i n c}}{\partial x} + D_{b o t, i g} = 0 & (6) \end{array}

The second term here is the simplified IG wave energy source term $S_{i g}$ , which can be written as Equation 7:

\begin{array}{l} S_{i g} = α_{i g} \sqrt{E_{i g}} \frac{c_{g}}{h} \frac{\partial S_{x x}}{\partial x} & (7) \end{array}

The third term is the dissipation due to bottom friction $D_{b o t, i g}$ , for which we use the formulation of Henderson and Bowen (2002):

\begin{array}{l} D_{b o t, i g} = f_{c w} ρ {(\frac{g}{h})}^{3 / 2} \frac{H_{r m s, i n c}}{\sqrt{8}} \frac{H_{r m s, i g}^{2}}{8} & (8) \end{array}

With $g$ the gravitational constant, $ρ$ the density of water, $H_{r m s, i n c}$ the incident root-mean-squared (rms) wave height, $H_{r m s, i g}$ the IG rms wave height, and $f_{c w}$ the friction coefficient, taken as $f_{c w} = 0.015$ following Swart (1974).

The parameterization to describe $α_{i g}$ is derived in the following subsection. With this, we obtain a full description of IG wave height growth over the cross-shore, with IG rms wave height as Equation 9:

\begin{array}{l} H_{r m s, i g} = \sqrt{8 E_{i g} / ρ g} & (9) \end{array}

IG wave breaking, which is not included in Equation 6, may become significant shoreward of the incident wave breaking point (e.g., van Dongeren et al., 2007). We therefore assume that Equation 6 is only valid from deep water until this point.

2.1.3 Parameterization of the shoaling parameter

We now seek a relation that describes the shoaling parameter $α_{i g}$ as a function of local conditions, based on observations of IG wave growth. For this, we use the VO21 dataset (see Section 2.2.1) which consists of 280 XBeach simulations covering a wide range of wave conditions and beach profiles (Figure 1). It contains values of incoming incident and IG wave height, bed slope, and water depth at 1-m depth contour intervals for each of the simulations. The training dataset, consisting of 70% of the 280 simulations, comprises a total of 4,441 training data points.

Figure 1

Figure 1 XBeach cross-shore profiles in the dataset of van Ormondt et al. (2021) for $H_{m 0} = 4 m$ , $T_{p} = 16 s$ , $\tan β_{b} = 0.10$ , and four variations of nearshore Dean slope $a_{D}$ (0.05, 0.10, 0.20, and 0.30) with indicated incident wave breaking point per profile in blue, adapted from van Ormondt et al. (2021).

For each point, the IG source term can be estimated using Equation 6, where the first term is obtained from the simulation results and the third term is estimated using Equation 8. Dividing the source term by the radiation stress gradient $\frac{\partial S_{x x}}{\partial x}$ yields the shoaling parameter at each point. Note that we use the conservative shoaling component of the radiation stress gradient rather than the observed value. This is done to avoid double counting as the observed gradient in the simulations already includes the effect of friction dissipation and energy transfer to the IG component.

Data points with a relative wave height $γ = H_{r m s, i n c} / h$ larger than 0.7 are excluded from the analysis. This was done because for these points IG wave breaking may be significant. Also, data points with a relative wave height $γ$ smaller than 0.1 are excluded, since for these points in deep water the increase in IG wave height is marginal and does not yield a clear trend. The result is a set of 2,778 observations of $α_{i g}$ for different combinations of $γ$ and corresponding local bed slope $β$ .

The observed shoaling parameter $α_{i g}$ for the training data are plotted against the local bed slope $β$ and color-coded with the relative wave height $γ$ (Figure 2A). The figure reveals a trend of decreasing values of $α_{i g}$ with increasing bed slopes. At mild bed slopes, the shoaling parameter is also dependent on the relative wave height $γ$ . Here, lower values of $γ$ , representing deeper water, result in higher shoaling parameter values. At relatively steep bed slopes, this dependency is less apparent. Therefore, the shoaling parameter $α_{i g}$ is hypothesized to be a function of the local values of $β$ and $γ$ :

Figure 2

Figure 2 Observed shoaling parameter $α_{i g}$ versus local bed slope $β$ and relative wave height $γ = H_{r m s, i n c} / h$ (colored dots) and determined fits for $α_{i g}$ depending on $γ$ (colored solid lines) (A) and scatterplot of predicted versus observed $α_{i g}$ for the training data (B) and the validation data (C).

\begin{array}{l} α_{i g} ~ f (β, γ) & (10) \end{array}

To parameterize Equation 10, we apply a non-linear least-square fitting procedure depending on $β$ and $γ$ and assume a negative exponential shape to represent the decay with increasing bed slope. The formulation distinguishes between shallower water depths ( $0.34 < γ < 0.7$ ) (green to red colors in Figure 2A) and deeper water depths ( $γ < 0.34$ ) (blue colors in Figure 2A), which have an extra term to represent the higher observed shoaling parameter for low values of $β$ . For very shallow water depths of $γ \geq 0.7$ (dark red colors in Figure 2A), $α_{i g}$ is set to 0 since those points were excluded, because IG wave dissipation is not accounted for in our simplified formulation. Furthermore, we assume that from there shoreward, there is no additional forcing to the IG waves, as the bound waves are released. For horizontal or negative bed slopes, we prescribe $α_{i g} = 0$ (Equation 12). In the case of negative bed slopes, bound IG waves may be released (Mei and Benmoussa, 1984), which we do not consider here. A maximum value of $α_{i g}$ = 1 is applied to prevent spurious values for very small $β$ in combination with small values of $γ$ . This gives:

\begin{array}{l} α_{i g} (β > 0, γ < 0.34) = 0.11 e^{- 17.7 β} (0.7 - γ) + (0.34 - γ) 0.017 / \sqrt{β} & (11) \end{array}

\begin{array}{l} α_{i g} (β > 0, γ \geq 0.34) = 0.11 e^{- 17.7 β} max (0.7 - γ, 0) & (12) \end{array}

The parameterization (colored lines in Figure 2A) produces an overall reasonable agreement with the observed shoaling parameter. The presence of the considerable scatter in the data introduces an inherent variability that the model may not fully capture. Furthermore, imperfections in the parametrization itself contribute to deviations from an idealized fit. Regression analysis shows an R² coefficient of determination of 0.566, a root-mean-squared error (RMSE) of 0.03, and a bias of 0.001 for the training data (Figure 2B). The skill is similar to the validation data, which includes all validation data points including the deepest water depths, with an R² value of 0.58 (Figure 2C).

2.2 Data

2.2.1 Synthetic data

XBeach dataset: The synthetic dataset of wave transformation and wave runup (van Ormondt et al., 2021) consists of 280 1D simulations using the incident wave-resolving XBeach non-hydrostatic two-layer model (Roelvink et al., 2009; de Ridder et al., 2021) for a large range of beach profiles and hydrodynamic conditions. The profiles are representative of sandy coasts and consist of a constant linear beach slope with angle $\tan β_{b}$ above the mean water level. Between mean water and a depth of $4 H_{m 0}$ , the profile consists of a Dean profile (Dean, 1991) with $z = - a_{D} x^{2 / 3}$ (Figure 1), with $x$ the cross-shore distance and $z$ the elevation of the profile. In deeper water, the profile has a constant 1:50 slope up until the depth where $n = c / c_{g} = 0.8$ , as recommended in de Ridder et al. (2021), with $c$ and $c_{g}$ representing the individual wave and the wave group velocity, respectively. Beyond that limit, the depth is constant. The 280 simulations consist of combinations of an offshore wave height ( $H_{m 0}$ ; 1, 2, 4, 8, 12 [m]), peak wave period ( $T_{p}$ ; 6, 8, 12, 16 [s]), Dean slope ( $a_{D}$ ; 0.05, 0.10, 0.20, 0.30 [–]), and beach slope ( $\tan β_{b}$ ; 0.01, 0.02, 0.05, 0.10, 0.20 [–]). Combinations with unrealistic (too steep) wave steepness values were discarded in the dataset. At the offshore boundary, a default JONSWAP spectrum with a peak enhancement factor of 3.3 is used in XBeach. The output consists of values of incoming incident and IG rms wave heights and wave periods at selected points in the cross-shore at every 1-m depth contour interval (e.g., 1, 2, 3 m, etc.). The incoming wave signal was obtained in van Ormondt et al. (2021), using the method of Guza et al. (1984). The definition of the IG wave cutoff frequency is half the peak frequency ( $f_{p} = 1 / 2 T_{p} f_{p} = 1 / 2 T_{p}$ ), the definition of which is also used throughout this study. For more details on the dataset, see van Ormondt et al. (2021). The IG wave heights are used for derivation of the parameterization describing the IG wave growth. The performance of XBeach in describing IG wave growth compared with observations is already shown in, e.g., de Ridder et al. (2021) and Roelvink et al. (2009), and for field cases, we refer to existing validations in the literature (e.g., Lashley et al., 2018; Roelvink et al., 2018).

ERA5 wave dataset: EMCWF’s ERA5 reanalysis dataset (Hersbach et al., 2020) is used for retrieving wave conditions on 13 October 1990 around the Outer Banks, USA, as used in the application case in Section 3.4. We use the significant wave height, peak wave period, and mean wave direction.

2.2.2 Laboratory data

GLOBEX dataset: The data of the GLOBEX experiments (Ruessink et al., 2013) consist of experimental laboratory experiments of wave transformation over a gently sloping linear bed with a 1:80 slope. We focus on cases A1, A2, and A3, the conditions of which can be found in Table 1 and are representative for real-world conditions of H_s = 2 m and T_p = 7 s for test A1, H_s = 4 m and T_p =10 s for A2, and H_s = 2 m and T_p =10 s for A3. From these data, we use the measured incident and IG wave heights over the profile. The incoming wave signal was obtained from de Bakker et al. (2015) using the method of Guza et al. (1984).

Table 1

Table 1 The GLOBEX and Boers datasets detailing the (average) slope, wave conditions, and skill metrics R², RMSE, and relative bias.

Boers dataset: The data of the experiments of Boers (1996) consist of laboratory experiments of wave transformation over a barred beach, with an average slope of about 1:38. We focus on cases 1A, 1B, and 1C, the conditions of which can be found in Table 1. From these data, we use the measured incident and IG wave heights over the profile. For both GLOBEX and Boers, to overcome data gaps and sudden jumps in the incident wave height measurements, the data are smoothed with a 10- and 5-point moving mean, respectively, to obtain useful data input as incident wave energy. The incoming wave signal was obtained in Boers (1996) using the method of Hughes (1993).

2.2.3 Field data

DELILAH dataset: The data of the DELILAH field campaign (Birkemeier et al., 1997) consist of alongshore varying 2D bathymetry data at the research facility in Duck, NC, representing the state on 13 October 1990, as, e.g., modeled by Van Dongeren (2003). The bathymetry consists of a breaker bar and small alongshore variations. From 13 October 1990, from 16:00 to 17:00, wind-still swell conditions with a significant wave height H_s of 1.81 m and a peak wave period T_p of 10.6 s were observed in 8 m water depth. These conditions are also used in other studies (e.g., Van Dongeren, 2003; Roelvink et al., 2009, 2018; Reyns et al., 2023) and are a result of Hurricane Lili. The incident wave angle is 88° from the north, corresponding to 16° from shore normal. The directional spreading factor s of the spectrum is approximately 6 (Roelvink et al., 2018), corresponding to 30.6 degrees. The local wind conditions consist of a wind speed of only 2 m/s (Birkemeier et al., 1997) and are therefore excluded in the model, as in Reyns et al. (2023). The offshore water level is 0.69 m+NGVD (National Geodetic Vertical Datum).

Topography and bathymetry datasets of North Carolina: For the application case of the Outer Banks, USA, we use the 1-m grid resolution Coastal National Elevation Database (CoNED) topographic model (Danielson et al., 2016; Tylor et al., 2022). For regions not covered by these data, the Continuously Updated Digital Elevation Model (CUDEM; CIRES, 2014), Coastal Relief Model (NOAA National Geophysical Data Center, 2001), and GEBCO (Becker et al., 2009) datasets are used.

3 Results

3.1 Validation of the infragravity wave source term against synthetic data

Validation of the prediction of IG wave transformation using the IG wave source term is performed for the synthetic profiles of VO21. For this, we include 1) a validation of the IG wave height at the incident wave breaking point and 2) a validation of the IG growth over the profile, through determining the shoaling rate.

With Equations 6, 8, 11, and 12, the IG wave growth from deep water toward nearshore can be computed. This is illustrated for one profile (Figure 3), where the shoaling parameter along the profile (Figure 3B) is determined based on the observed incident wave height (Figure 3A; black dots) and bathymetrical conditions (Figure 3C). The resulting computed IG wave height, obtained by numerically integrating Equation 6 (Figure 3A; blue line), compares well with the observed IG wave heights (Figure 3A; black triangles). Further visual inspection of the IG wave growth for more profiles and wave conditions is done in Appendix A.

Figure 3

Figure 3 Example of the observed incident and infragravity wave heights (in black) and predicted infragravity wave height (in blue) (A), predicted shoaling parameter $α_{i g}$ (B), and bed level with indicated incident wave breaking point (C) for XBeach simulation $H_{m 0, i n c} = 4 m$ , $T_{p} = 16 s$ , $a_{D} = 0.2, \tan β_{b} = 0.10$ in the dataset of van Ormondt et al. (2021).

To quantify this overall 280 profiles, the observed and predicted IG wave heights at the incident wave breaking point $h_{b r}$ are compared, as illustrated in Figure 4. The incident wave breaking point here is defined as the location where the incident wave height $H_{m 0}$ over water depth ratio equals 0.7 (dashed line in Figure 3C). For both the training data (Figure 4A) and the validation data (Figure 4B), the R² values are similar with a low RMSE and relative bias, indicating that the results are not overfitted and the parameterization has predictive skill. For the validation profiles, this yields an R² value of 0.964, an RMSE of 0.054 m, and a bias of −0.011 m. Overall, 94.4% of the training data and 92.9% of the validation data are within the +−25% error margins.

Figure 4

Figure 4 Comparison between predicted and observed infragravity wave heights at the wave breaking point $h_{b r}$ (in blue) for the training data (A) and validation data (B), as well as a comparison between the predicted and observed shoaling rate over the entire profile (red), as fitted following the method in van Dongeren et al. (2007) for the training data (C) and the validation data (D). The dashed line shows the perfect fit line, while the solid black lines indicate the +−25% error margins.

A second validation considers the IG wave shoaling rate $m_{i g}$ , defined as $H_{r m s, i g} ~ h^{- m_{i g}}$ . This is a single parameter used to evaluate the skill of predicting the IG wave evolution. The parameter $m_{i g}$ can be determined based on known values of $H_{r m s, i g}$ and $h$ following the method as in van Dongeren et al. (2007). These values are taken from deep water up to the incident wave breaking point. This gives one value for the shoaling rate $m_{i g}$ per profile, for both the observed and predicted data, not to be confused with the local shoaling parameter $α_{i g}$ as defined in Section 2. For both the training (Figure 4C) and the validation data (Figure 4D), the shoaling rate $m_{i g}$ is between 0.5 and 1.5. The performance of the training and validation data is comparable in terms of patterns and skill score, indicating that the results are not overfitted and the IG wave source term has predictive skill. The RMSE of the validation data is 0.123, and the relative bias indicates a minor underestimation of 8.4%. Overall, 96.4% of the training data and 92.9% of the validation data are within the +−25% error margins.

3.2 Validation of infragravity wave transformation in 1D against laboratory data

In this section, we compare the results of the introduced wave solver for predicting IG wave transformation (Equation 2) against laboratory data. The data are obtained from three test cases of the GLOBEX experiments (Ruessink et al., 2013) and three test cases of the Boers experiments (Boers, 1996) (see Figure 5). For each case, the IG wave transformation is computed using the observed incident wave height as input to the IG source term, rather than the values computed with Equation 1. By constraining the model in this manner, we ensure that the validation focuses solely on the IG source term and is not influenced by possible errors in the computed incident wave transformation. The wave conditions and profile characteristics are summarized in Table 1, and the model settings are summarized in Appendix B.

Figure 5

Figure 5 Observed incident and infragravity wave heights (in black) and predicted infragravity wave height (in blue) for GLOBEX tests A1 (A), A2 (B), and A3 (C) including the bed level profile (D) and for Boers tests 1A (E), 1B (F), and 1C (G) including the bed level profile (H).

The evolution of the IG wave heights for the GLOBEX test cases is quite well predicted (blue line compared with black dots), which is reflected in the skill scores (Table 1). The R² is 0.91–0.97, RMSE 0.7–4.6 mm, and relative bias 0.7%–22.4%. The IG wave transformation is especially well predicted up until the surf zone, within the region where the incident waves have started breaking. In test A3, the IG wave breaking nearshore is underestimated.

The development of IG wave heights for the Boers tests with a barred beach profile turns out to be more difficult to model correctly (Figure 5H). For test cases 1A and 1B (Figures 5E, F), the general trend and order of magnitude are still reasonably well predicted. In test case 1A, the nearshore IG wave dissipation is underestimated. In test case 1C (Figure 5G), there is an overestimation of the predicted IG wave growth. This case combines the relatively steep profile (compared with GLOBEX), with a significantly lower wave steepness than the other Boers tests (Table 1). The overall R² is 0.78–0.86, RMSE 4.0–9.3 mm, and relative bias −12.3% to +49.9% (Table 1).

3.3 Validation of the wave solver in 2D against field data

The prototype scale validation of the wave solver involves the DELILAH field experiment (Birkemeier et al., 1997) at Duck, NC. This site includes a non-uniform alongshore bathymetry, energetic wave conditions, non-oblique mean wave direction, and directional wave spreading. The 2D stationary wave solver (Equations 1, 2) is used to compute the incident and IG wave transformation at a spatially uniform 5-m grid resolution (as the finest resolution as used in Reyns et al., 2023). For more details on the settings of the wave solver, see Appendix B.

Spatial variations in bathymetry (Figures 6A, E) lead to alongshore and cross-shore variations in both the modeled incident (Figures 6B, D) and the modeled IG wave heights (Figures 6C, D). While the modeled incident and IG wave heights differ from the observations, the results predict the overall magnitude of wave heights reasonably well (Figures 6B–D). In between 600< x< 700 m, the incident wave height is underestimated, which can also be observed in the results of XBeach SurfBeat (Roelvink et al., 2009) and XBeach non-hydrostatic (Roelvink et al., 2018). Consequently, the IG wave height is also underestimated in that region. The offshore wave height evolution cannot be validated because of a lack of measurements. In shallow water depths landward of the breaker bar (x > 750 m, Figure 6E), the breaking of the IG waves and conservative shoaling is well captured. The RMSE over all eight measuring locations combined is 0.108 m for the incident and 0.042 m for the IG wave heights.

Figure 6

Figure 6 DELILAH case: bed level elevation including observation stations (A), observed (circles) and modeled incident wave heights (B), observed (circles) and modeled infragravity wave heights (C), observed and modeled incident and infragravity wave heights along the transect of measurements (D), and cross-shore elevation including observation stations (E).

Overall, the results have similar accuracy as Reyns et al. (2023) and Roelvink et al (Roelvink et al., 2009, 2018), who modeled the IG wave heights using a dynamic wave-resolving model.

3.4 Application of the wave solver in 2D for a large-scale domain

We now demonstrate the application of the wave solver on large domains. To this end, a model was constructed to compute nearshore IG wave conditions for 450 km of coastline. It covers the entire Outer Banks in North Carolina and most of the Virginia coastlines. We force the model with the wave conditions on 13 October 1990, using the wave reanalysis data of ERA5 (Hersbach et al., 2020), where the mean wave direction is from the east. The area model has an offshore boundary at the 200-m water depth contour at the edge of the continental shelf. A spatially varying grid resolution is used, containing 1,000 m grid cells in deep water, increasing in resolution in steps of a factor of 2, up to the finest grid cells of 31.25 m in the surf zone. The total domain has 4,493,839 active grid cells.

Various bathymetric features that affect the wave field are discernible in the bathymetry (Figure 7A). These include barrier islands, breaker bars, and shoals. Behind the Cape Lookout Shoals (orange triangle in Figure 7B), a sheltered area is visible where the wave height is lower. The computed IG wave height field (Figure 7C) shows that the highest IG wave heights occur close to the shore, where the incident wave shoaling is the strongest. The shoaling is particularly strong just north of Duck (red triangle in Figure 7C), where a range of nearshore bathymetric features lead to large changes in radiation stress, leading to stronger IG wave growth. At the Cape Lookout Shoals, the waves from the east shoal up and generate IG waves. In general, the patterns of the IG wave fields are consistent with the bed level and incident wave height changes.

Figure 7

Figure 7 Outer Banks case: bed level elevation (A), modeled incident wave heights (B), and modeled infragravity wave heights (C) for the Outer Banks upscaling demonstration model. Duck, NC, is indicated in the red triangle and Cape Lookout Shoals, NC, in the orange triangle.

The computation of a single stationary wave field (combining both incident and IG waves) of 4.5 million grid cells takes mere seconds on a regular laptop PC. This is significantly faster than computations by advanced wave-group-resolving numerical models that would take hours to days to calculate this for a domain of this size.

4 Discussion

The presented wave solver provides reasonable to good estimates of the nearshore IG wave conditions for large spatial scales and with a low computational expense. Since a reduced complexity approach is taken, it is not a full description of the physics involving the generation of IG waves (as is the case for most numerical models). The main limitations, assumptions, and discussion points of both the IG wave source term and the full 2D wave solver are outlined below:

4.1 Limitations of the infragravity wave source term

The proposed IG wave source term using a parameterization for the shoaling parameter $α_{i g}$ , based on the local bed slope and wave height over water depth ratio, is a simplified physical description. Therefore, the goodness of the fit of the parameterization of the shoaling parameters is not very high with an R² score of 0.580 (Section 2.1.3). However, when predicting the IG wave development from offshore to nearshore, the uncertainty in the exact value of $α_{i g}$ matters less, which is indicated by an R² value of 0.964 in estimating the IG wave height at the incident wave breaking point and a small bias of 1.1 cm (Section 3.1).

The dataset of VO21 used for deriving and validating the parameterization consists of 1D XBeach non-hydrostatic two-layer model results. Shore normal 1D models overestimate the IG wave growth compared with the 2D XBeach models (van Ormondt et al., 2021; McCall et al., 2023), because all energy is forced into the shore normal wave direction and directional spreading is not considered. Furthermore, the parameterization assumes a JONSWAP spectral shape in the incident wave band, as is used in the VO21 dataset, which is applicable to nearshore conditions. The evolution of IG wave heights and the corresponding IG wave source term may differ for variations of this spectral shape and subsequent groupiness. The effect of directional spreading and wave groupiness on IG wave shoaling requires further investigation.

In the IG source term, the shoaling parameter (and therefore the source term itself) becomes zero for negative bed slopes (Equation 12). These occur where the water depth increases in the wave propagation direction, e.g., on the landward side of a submerged breaker bar or when IG waves enter a deeper tidal channel after passing a shallow shoal (Reniers and Zijlema, 2022). These conditions are not present in the VO21 dataset from which the parameterization was derived. The choice to set the source term to zero for these cases was made to keep the IG source term simple and not have to account for the possibility that bound IG waves may become released. This aspect of the parameterization is however likely inaccurate and should be further investigated.

No clear correlation was found in the VO21 dataset between incident wave steepness and the IG source term. The steepness is therefore not explicitly considered in our parameterization. However, some cases with low wave steepness, such as Boers test 1C, showed an overestimation of the predicted IG wave growth. It is recommended that the effect of wave period and wave steepness on IG wave transformation is further investigated.

4.2 Limitations of the wave solver

In the wave solver, IG wave-breaking dissipation was included using the formulation of Baldock et al. (1998). This formulation was derived for sea swell and has not been validated for IG waves. Further investigation into the validity of the formulation is needed and might also improve results, e.g., GLOBEX test A3 and Boers tests 1A and 1C, which would benefit from introducing more IG wave breaking.

Break point forcing may contribute significantly to IG wave generation at steep rocky coasts and coral-reef-lined coastlines (Battjes et al., 2004; Van Dongeren et al., 2013). A source term to represent this process has not yet been implemented in the wave solver. As a result, the solver’s applicability is currently limited to coasts with relatively mild slopes.

As with other numerical models, the present wave solver requires a relatively good estimate of the nearshore bathymetry, resolving, e.g., breaker bars. This in contrast with other approaches that use a certain mean (surf zone or beach) bed slope to represent the coastal profile. These bathymetric details are not covered in coarser global datasets (e.g., GEBCO) or derived profile datasets (e.g., Athanasiou et al., 2019, 2023). However, the sensitivity between estimates of IG wave growth derived from using a detailed or coarser bathymetry has not been assessed and could be considered for future work.

The presented wave solver has not yet been validated in 2D during extreme conditions as the significant wave height of 1.8 m for Duck was relatively limited. However, because the IG wave source term has been derived and validated based on conditions in XBeach with significant wave heights up to 12 m, the derived wave solver is assumed to be valid for extreme events along sandy beach coastlines. Field validation could be considered as a direction for future research.

5 Conclusions

In this study, we introduce an efficient stationary wave energy solver designed to estimate the evolution of incident and IG waves from offshore to the nearshore, specifically for mildly sloping coasts. This wave solver can be used to generate nearshore IG wave boundary conditions to drive overland flood models at large scales and with little computational expense.

The new wave solver builds upon the stationary wave energy balance for incident wave energy. The formulations are extended to also include the IG wave energy balance. For this purpose, we introduce an IG wave source term that describes the energy transfer from incident to IG waves. This term acts as a sink term for the incident waves and as a complementary source term for the IG waves. The IG wave source term is simplified using an IG wave shoaling parameter $α_{i g}$ . Observed values of this parameter are obtained from a synthetic dataset of covering a wide range of wave conditions and beach profiles. These are used to derive an empirical relation, which describes the shoaling parameter as a function of the local bed slope $β$ and the incident wave height over water depth ratio $γγ$ . While the parameterization of $α_{i g}$ exhibits some scatter relative to observed values (R² value of 0.566), it captures the general trend without overfitting. Most importantly, employing this parameterization accurately predicts the IG wave source term and subsequent IG wave height evolution across all profiles in the validation dataset. The shoaling rate over the cross-shore profile is predicted with an RMSE of 0.123 with 92.9% of the validation data falling within the +−25% error margins. IG wave heights at the incident wave breaking point are predicted with high accuracy, indicated by an R² value of 0.96 and an RMSE of 0.05 m.

With this parameterization, the wave solver can estimate IG wave transformation based on offshore incident wave heights and local bathymetry as input. The wave solver is shown to be able to predict the evolution of IG wave heights for three laboratory tests of GLOBEX (RMSE between 0.7 and 4.6 mm, R² of 0.91–0.97) and for the more complex bed level of three laboratory tests of Boers (RMSE between 4.0 and 9.3 mm, R² of 0.78–0.86). Additionally, validation of the wave solver in 2D using the 1990 DELILAH field case observations at Duck, NC, USA, yields results that align well with those obtained with computationally expensive dynamic wave models. Predictions of incident wave height yield an RMSE of 0.11 m, while IG wave height estimates exhibit an RMSE of 0.04 m.

The large-scale applicability of the wave solver is further demonstrated by modeling the entire Outer Banks, North Carolina, coastline in the USA, covering 450 km of coastline from 200 m water depth up to the coast. For this model, consisting of 4.5 million grid cells, the wave solver can calculate IG wave conditions for the entire domain within seconds.

The presented wave solver makes it possible to include IG wave-driven flooding in large-scale compound flood modes. Its computational efficiency, with results in seconds, stands in contrast to process-based wave-resolving models, which require hours to days to run. In future work, this wave solver will be used to force a reduced complexity model (SFINCS), for IG wave-resolving flood hazard modeling at the national to continental scales.

Data availability statement

The developed wave solver has been integrated in the fast compound flood model SFINCS, which is openly available from https://www.deltares.nl/en/software-and-data/products/sfincs. The source code of the coupled model with wave solver and infragravity source term as developed in this paper is available from https://github.com/Deltares/SFINCS/releases/tag/v2.0.5-alpha-branch-sfincs_snapwave_ig_2024_02_10 and will be incorporated in future releases of SFINCS. The data analyzed in this study are subject to the following licenses/restrictions: The used dataset of XBeach simulations was not made available open source by van Ormondt et al. (2021); however, the data are generated using the open-source XBeach model. Access to the dataset can be requested on a personal basis. Requests to access these datasets should be directed to bWFhcnRlbi52YW5vcm1vbmR0QGRlbHRhcmVzLXVzYS51cw==.

Author contributions

TL: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft. MO: Conceptualization, Methodology, Investigation, Software, Writing – review & editing. AD: Conceptualization, Methodology, Supervision, Writing – original draft, Writing – review & editing. JA: Supervision, Writing – review & editing, Funding acquisition. SM: Supervision, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This project has been funded by EU ERC COASTMOVE nr 884442 and by Deltares SITO-IS research funding under Moonshot 2 – Flooding.

Acknowledgments

We would like to thank Johan Reyns for providing the input data as used in the Duck field case, as well as Carola Seyfert for making the initial version of the model setup. Permission to use the data provided by the Field Research Facility of the U.S. Army Engineer Waterways Experiment Station’s Coastal Engineering Research Center is greatly appreciated.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Athanasiou P., Van Dongeren A., Giardino A., Vousdoukas M., Antolinez J. A.A., Ranasinghe R. (2022). Estimating dune erosion at the regional scale using a meta-model based on neural networks. Natural Hazards Earth Sys. Sci. 22, 3897–39155. doi: 10.5194/nhess-22-3897-2022

ORIGINAL RESEARCH article

Estimating nearshore infragravity wave conditions at large spatial scales

1 Introduction

2 Materials and methods

2.1 Methodology

2.1.1 Formulation of the 2D wave solver

2.1.2 Infragravity wave source term

2.1.3 Parameterization of the shoaling parameter

2.2 Data

2.2.1 Synthetic data

2.2.2 Laboratory data

2.2.3 Field data

3 Results

3.1 Validation of the infragravity wave source term against synthetic data

3.2 Validation of infragravity wave transformation in 1D against laboratory data

3.3 Validation of the wave solver in 2D against field data

3.4 Application of the wave solver in 2D for a large-scale domain

4 Discussion

4.1 Limitations of the infragravity wave source term

4.2 Limitations of the wave solver

5 Conclusions

Data availability statement

Author contributions

Funding

Acknowledgments

Conflict of interest

Publisher’s note

References

Appendix A: further inspection of infragravity wave growth for profiles

Appendix B: boundary conditions and model settings of the wave solver

Boundary conditions of infragravity waves

Model settings of the GLOBEX and Boers lab tests

Model settings of the Duck and Outer Banks case studies