Royal Society Publishing

Reynolds averaged Navier–Stokes modelling of long waves induced by a transient wave group on a beach

Javier L. Lara , Andrea Ruju , Inigo J. Losada


This paper presents the numerical modelling of the cross shore propagation of infragravity waves induced by a transient focused short wave group over a sloping bottom. A dataset obtained through new laboratory experiments in the wave flume of the University of Cantabria is used to validate the Reynolds averaged Navier–Stokes type model IH-2VOF. A new boundary condition based on the wave maker movement used in the experiments is implemented in the model. Shoaling and breaking of short waves as well as the enhancement of long waves and the energy transfer to low-frequency motion are well addressed by the model, proving the high accuracy in the reproduction of surf zone hydrodynamics. Under the steep slope regime, a long wave trough is radiated offshore from the breakpoint. Numerical simulations conducted for different bottom slopes and short wave steepness suggest that this low-frequency breakpoint generated wave is controlled by both the bed slope parameter and the Iribarren number. Moreover, the numerical model is used to investigate the influence that a large flat bottom induces on the propagation pattern of long waves.

1. Introduction

Sea surface records show the presence of groups of high and low waves as a result of the superposition of wave components with different frequencies. Since in sufficiently deep water, gravity waves are frequency dispersive, the groups are transient. Although the presence of transient groups in the wave field is the result of merely linear processes, nonlinear effects are not negligible in the focused wave groups. Local nonlinear interactions in conjunction with the focusing of wave energy lead to large waves with high and narrow crests that can break even in the deep ocean (Rapp & Melville 1990). Moreover, global nonlinear interactions produce low-frequency changes such as slow varying currents and oscillations of the mean water level known as infragravity waves.

In nature, typical periods of wind waves are in the order of 10 s, while infragravity motion evolves slowly with a typical period of several minutes. Two main generation mechanisms have been identified for these low-frequency waves. The first mechanism contemplates infragravity bound waves forced by short wave groups through the radiation stress gradient (Longuet-Higgins & Stewart 1962). The other main generation mechanism takes into account the radiation of free long waves both seaward and shoreward from a time varying breakpoint (Symonds et al. 1982; Schäffer 1993). Due to the low wave steepness, infragravity waves oscillate over the beach face, experiencing an almost complete reflection. It is recognized that the long period surface oscillation in the nearshore region correlated to short wave groups, called ‘surf beat’ (Munk 1949; Tucker 1950), has a strong influence on coastal processes such as beach erosion and bay resonance.

The cross shore propagation of infragravity waves over an uneven bottom has been extensively investigated by means of laboratory experiments. Baldock et al. (2000) and Baldock & Huntley (2002) corroborated the breakpoint forcing model first derived by Symonds et al. (1982). They also provided the physical conditions for which the breakpoint forcing should become ineffective. Janssen et al. (2003) investigated the phase lag between the short wave envelope and bound waves using high spatial resolution laboratory data. A theoretical linear model to predict this lag was presented, and model results showed good qualitative agreement with observations outside the surf zone. Battjes et al. (2004) studied the rate of energy transfer between primary waves and bound waves. According to Baldock et al. (2000) and Baldock & Huntley (2002), they differentiated between a mild slope regime, for which incident bound waves dominate over breakpoint generated long waves, and a steep slope regime in which the opposite is expected. van Dongeren et al. (2007) showed that in a mild slope regime, the long wave energy dissipation produced by breaking of long waves can overcome the dissipation owing to the bottom friction. Baldock (2006) presented a laboratory dataset where a single wave group was generated and focused close to the shoreline. The study of the transient wave group propagation over a steep slope enabled direct observation of the long wave breakpoint generated. He observed a negligible time domain lag between the short wave envelope and the bound wave.

In addition to laboratory experiments, the study of nearshore dynamics and long wave transformations has been approached with mathematical models. Traditionally, nonlinear shallow water (NSW) and Boussinesq-type equations models have been employed for the study of infragravity waves in the nearshore zone, providing satisfactory results. However, both approaches require setting the triggering wave breaking mechanism and the subsequent wave energy dissipation owing to wave breaking based on empirical formulations. A concise review in terms of their strengths and limitations can be found in Torres-Freyermuth et al. (2010).

In recent years, models based on Reynolds averaged Navier–Stokes (RANS) equations have been developed to describe nearshore processes with the aim to overcome the limitations related with NSW and Boussineq models. Due to the fewer assumptions involved in the governing equations, RANS models are able to simulate the frequency dispersive nature of gravity waves in deep waters as well as the nonlinear wave transformations in shallow waters over a sloping bottom. The refined information on the velocity, pressure and turbulence field makes them suitable to study surf zone hydrodynamics. Wave breaking and its evolution along the surf zone are directly solved without any imposed forcing. Although RANS models have focused traditionally on the study of wave and structure interaction (Lara et al. 2008; Losada et al. 2008), in recent years some attempts to model nearshore hydrodynamics on natural beaches (Torres-Freyermuth et al. 2007) have been carried out. More recently, Torres-Freyermuth et al. (2010) have analysed the performance of RANS equations to simulate long wave transformation on a barred beach. They found satisfactory agreement between numerical results and observations, although some tasks are still open related with the modelling of long wave transformation in the nearshore area.

It is the aim of this work to improve our understanding on the propagation pattern for long waves induced by a transient short wave group over a sloping bottom using a numerical model. In order to deal with infragravity motion in the nearshore zone, a high spatial and time resolution is required and a RANS model has been extended in order to test the capability in simulating hydrodynamic processes generated by a focused wave group. The use of a new laboratory dataset has been proposed to validate the model. This paper is organized as follows. Section 2 illustrates the laboratory experiments. Next, the numerical model is presented and described in §3. In §4, model validation is conducted by comparing the numerical results and the laboratory measurements. Long wave transformation analysis is carried out in §5. Finally, some conclusions are provided in §6.

2. Laboratory experiments

(a) Experimental set-up

The laboratory experiments were carried out in a wave flume at the University of Cantabria. The flume is 22 m long, 0.5 m wide and 0.8 m high. The still water depth during the experiments was set at 0.4 m at the wave maker position. Waves were generated by a hydraulically driven piston-type paddle located at one side of the flume. The paddle movement was controlled with a second-order signal wave generator in order to suppress spurious long waves and reproduce long bound waves (Ottensen-Hansen et al. 1980; Barthel et al. 1983). Wave absorption at the wave maker was turned off.

The fixed bathymetry is made up of a horizontal bottom part and a slope-varying part that resembles a relatively steep beach profile. The part with a constant depth of 0.4 m begins at the wave maker and is 8.6 m long. A beach with three different slope regions extends up to the end of the flume. The first region has a constant slope of 1 : 20 and is 4 m long; the second region with a horizontal bottom and a depth of 0.2 m is 2 m long; finally the last region with the same slope, 1 : 20, as the first one reaches the end of the flume. The origin of the co-ordinate system is fixed at the intersection of the still water line with the bottom, with the xz axes positive shoreward and upward, respectively.

Thirteen surface-piercing resistance wire wave gauges mounted on a carriage above the flume were used to measure the fluid surface displacement at fixed spatial locations. The first three gauges, positioned between the wave paddle and the toe of the beach, were used as control gauges. The other 10 gauges were placed over the varying slope region. Data from the wave gauges were collected with a sampling frequency of 60 Hz to provide a high temporal resolution. Each experiment was 300 s long in order to achieve a high spectral resolution for frequency domain analysis. Video recordings of breaking waves and run-up were also carried out during the experiments. The bottom profile with the position of the gauges is shown in figure 1. The detailed gauge locations are summarized in table 1. Wave gauges are identified by the acronym ‘G’ followed by the gauge number ordered from the paddle to the shoreline.

Figure 1.

Wave flume and instrumentation.

View this table:
Table 1.

Location of surface gauges.

(b) Wave characteristics

A large transient group constituted by N=50 individual wave components of equal amplitude was generated. The primary short wave components were uniformly spaced over a specified band of the frequency domain in order to obtain a ‘top-hat’ spectral shape. Total wave group amplitude A, frequency bandwidth Δf and central frequency fc were defined as Embedded Image 2.1 where an is the amplitude of the nth frequency component. Since wave components had the same amplitude an=a, these parameters were varied in order to obtain wave groups with different characteristics. We proposed:

  • — Three group heights: H=(0.06,0.11,0.16) m.

  • — One frequency bandwidth: Δf=0.167 Hz.

  • — Two central frequencies: fc=(0.583,0.75) Hz.

Thus six cases for second-order wave generation were obtained; these are summarized with the corresponding Iribarren number Ir in table 2. The Iribarren number was determined as Embedded Image 2.2 where α=1:20 is the bottom slope, H=2A the total wave group height and Embedded Image the deep water wave length calculated with the central frequency fc. It is worth mentioning that for the cases with H=0.11 m (B cases) the Iribarren number calculation is theoretically inappropriate since in these cases the breakpoint is located above the horizontal bottom (see further description in the following). As a general rule, a lower limit for the plunging breaker type is defined as Ir≈0.4. However, owing to the strong nonlinearity occurring during focusing events, it seems reasonable to fix the lower limit to Ir≈0.3. In fact, during the experiments, short wave plunging-type breaking occurred in case A and a transition between spilling and plunging-type was observed in cases B and C.

View this table:
Table 2.

Wave groups conditions.

The period of the wave packets that comprise the surface elevation signal decreases downstream until reaching the minimum value of Tg=2/Δf at the focal point (Rapp & Melville 1990). The number of waves in the group was calculated as Nw=2fc/Δf. Table 3 summarizes wave group period and the number of waves in the group for the two different values of Δf/fc.

View this table:
Table 3.

Wave groups characteristics at the focal point.

For the condition considered, the wave-breaking location is controlled primarily by wave height. So breaking of short waves occurs at x≈−6.5,−5 and −2 m for cases C, B and A, respectively. This means that wave breaking occurs in intermediate depths (π/10<kh<π).

Following the Rapp & Melville (1990) procedure, the phase of each wave component at the wave paddle was determined so that, according to linear theory, constructive interference occurs at a point down the flume. Therefore, the first-order surface elevation at the wave paddle was calculated as Embedded Image 2.3 in which xp is the paddle location, kn and ωn are the wave number and radian frequency of the nth component wave, tf is time at which wave focusing occurs and xf is the focus position that was determined to coincide with the breakpoint location (xfxb). Identifying the exact location of the point where the individual wave components are brought into focus is not a trivial task (Baldock et al. 1996), since it is affected by both the bathymetry and the nonlinear interactions between wave components. A heuristic approach was used to focus the wave energy at x≈−5 m, close to the breakpoint of the short waves.

3. Numerical model

The use of a RANS model is motivated to study nearshore processes among other approaches because nonlinearity and wave breaking are directly solved with any imposed assumptions. In order to reproduce gravity and infragravity wave hydrodynamics, a detailed treatment of the boundary conditions has to be carried out. The simultaneous wave generation and absorption problem is of great importance for a detailed study of infragravity waves. The difficulty is in specifying incident waves through an inflow boundary with the presence of gravity wave reflection and offshore long wave radiation. Although several approaches have been focused on gravity waves, combining an internal wave maker and a sponge layer, first introduced by Lin & Liu (1999) and later enhanced by Lara et al. (2006), infragravity waves cannot be accurately generated or absorbed. More recently, Torres-Freyermuth et al. (2010) combined a method based on imposing the free surface and the velocity field with active absorption. Results showed a good reproduction of both gravity and infragravity wave transformations. However, the velocity field above the trough level must be extrapolated introducing additional uncertainties.

Since several physical wave studies (Baldock 2006 and others) have obtained successful results dealing with infragravity wave generation, a moving boundary is implemented in order to replicate wave maker movement during the experiments. No additional assumptions have to be made to define the wave velocity or the free surface, ensuring that the fluid is forced identically in the numerical model and in the experiments which is of great importance in proving the model performance. Active wave absorption is also incorporated to the numerical wave paddle. Numerical model details, wave generation and absorption methods are given next.

(a) General description

IH-2VOF model, an updated version of COBRAS-UC (Losada et al. 2008), solves the 2DV Reynolds averaged Navier–Stokes (RANS) equations based on the decomposition of the instantaneous velocity and pressure fields into mean and turbulent components, and the kϵ equations for the turbulent kinetic energy, k, and its dissipation rate, ϵ. The influence of turbulence fluctuations on the mean flow field is represented by the Reynolds stresses. The governing equations for kϵ are derived from the Navier–Stokes equations, and higher order correlations of turbulence fluctuations in k and ϵ equations are replaced by closure conditions. A nonlinear algebraic Reynolds stress model is used to relate the Reynolds stress tensor and the strain rate of mean flow. The free surface movement is tracked by the volume of fluid (VOF) method.

The RANS equations are solved in a rectangular structured mesh. Pressure, p, and velocity, u, variables are defined in a staggered scheme. The model uses a cutting cell method first introduced by Clarke et al. (1986) in order to consider obstacles contained within the numerical domain. The basic idea behind this technique is that the obstacle can be modelled as a special case of the flow with an infinite density. This method improves the conventional approach of treating solid bodies defined by a sawtooth-shape. In order to describe the geometry of the mesh immersed bodies, a new set of coefficients needs to be defined, called openness functions, θ, which represent the fraction of volume of open space in this cell. Variables in the cell or on the cell faces are redefined by the product of the openness coefficients times the original variables. The governing RANS equations are redefined as follows: Embedded Image 3.1 and Embedded Image 3.2 where θ is the openness function, ρ is the fluid density, gi is the ith component of the gravitational acceleration and τij represents the sum of the viscous and Reynolds stresses. According to this definition, θ=0 is a ‘solid cell’ (entirely body occupied), θ=1 is a ‘fluid cell’ and 0<θ<1 is a ‘partial cell’. In the numerical implementation, partial cell coefficients need to be redefined at the cell centre and also at the cell boundaries, θr, θl, θb and θt, which correspond to the openness at the right, left, bottom and top, respectively.

Among the different numerical techniques used to simulate moving bodies within the computational domain, called ‘virtual force method’ (see Mittal & Iaccarino 2005 for a review), ‘direct forcing method’ (Mohd-Yusof 1997) is considered in the IH-2VOF model. This technique presents an advantage above others applied in computational fluid mechanics. The pressure field induced by the fluid and the moving boundary can be calculated by means of the Poisson pressure equation (PPE), avoiding iteration to close the velocity–pressure coupling, yielding a divergence-free velocity field. Moreover, no remeshing is required; the fluid is solved in a fixed-grid system and the method is also free of numerical diffusion (Lin 2008). The moving body velocity, ub,i, in this method needs to be known in advance, making this approach suitable to be applied to wave generation as discussed later.

The momentum equation (3.2) presents an additional body force, fb, which represents the virtual boundary force which acts only on the solid boundary of the partial cells (0<θ<1): Embedded Image 3.3 The two-step projection method (Chorin 1968) is applied. The first step consists in introducing an intermediate velocity in the momentum equation not considering the virtual boundary force: Embedded Image 3.4 where Embedded Image is the intermediate velocity, the superscript indicates the time level and Δt is the time step size corresponding to the (n+1)th time level. In equation (3.4), the openness function θ is written according to its time evolution. The second step can be written as follows: Embedded Image 3.5 Replacing Embedded Image by Embedded Image at the moving body boundaries, the virtual boundary force, fb, can be defined as Embedded Image 3.6 The pressure gradient can be evaluated at the (n+1)th time level, taking the divergence of equation (3.5) and applying the conservation of mass equation (3.1), yielding the PPE: Embedded Image 3.7

In the two-step projection method, the spatial derivations of the velocity components and the pressure field need to be expressed in finite-difference forms. The convection terms are discretized by the combination of the central difference method and upwind method. The combination of both is aimed at preventing their respective drawbacks of significant numerical damping and numerical instability. A weighting factor is introduced in the spatial derivative discretization expressions in order to adjust the influence of each one of the two schemes in the computation and to obtain stable and accurate solutions. Only the central difference method is employed to discretize the pressure gradient terms and the stress gradient terms.

(b) Wave maker implementation

Although piston-type wave maker implementation is explained in this section for the sake of space, its extension to flap type or mixed piston-flat type is straightforward. The wave maker movement is transferred into the model through the openness coefficients. Since wave maker position, X(t), and velocity, Ub(t), are known values, openness coefficients, θ,θr, θl, θb and θt, are then redefined according to the new wave maker position at any given time step. In the particular case of a piston-type wave maker position, which generates wave from left to right, θ=θb=θt, θl=0 and θr=1. A sketch of the openness functions in a piston-type wave maker is presented in figure 2 for a one-dimensional flow. According to the aforementioned classification, the left cell is a solid cell, the right cell is a fluid cell and the middle cell is a partial cell. Figure 2 represents the three adjacent cells in the wave maker front, where virtual force, fb is acting. fb is applied only in the middle cell.

Figure 2.

Sketch of openness function configuration in a piston-type wave paddle.

The right-hand side term of equation (3.7) is then discretized for the middle cell as Embedded Image 3.8 The moving wave maker appears as a source term in the momentum equation determined by the body surface velocity. This source term takes into account the velocity- and pressure-induced fields in the fluid.

(c) Active wave maker absorption

Active wave maker absorption is also included in the model using the same procedure followed in physical flumes. The numerical wavemaker is equipped with a control system for simultaneous wave generation and active wave absorption. The methodology proposed by Schäffer & Klopman (2000) is followed.

The target wave maker excursion, X(t), is corrected in time in order to absorb outgoing waves and to avoid reflection at the wave maker. The reflected wave is estimated comparing the target-free surface with the free surface recorded in front of the wave maker. A location three cells away from the wave maker location is chosen. The reflected free surface is then evaluated as Embedded Image 3.9 where ηr is the reflected wave, ηtarget is the target free surface and ηm is the measured wave in front of the wave maker. Wave maker velocity (U(t)=dX(t)/dt) has to be modified in order to match the velocity induced by the wave to be absorbed. In this case, linear wave theory is proposed, and the velocity correction owing to absorption can be written as follows: Embedded Image 3.10 To obtain the desired wave maker position, velocity has to be integrated considering both the target velocity and the correction by absorption Embedded Image 3.11

4. Model validation

The model validation procedure is conducted by comparing the numerical model predictions and the previously described laboratory measurements. The purpose is to identify the capabilities and the limitations of the IH-2VOF model in the simulation of the long period processes induced by a transient focused wave group over a sloping bottom. Since an accurate modelling of the short waves is required to predict infragravity wave motion, considerable attention is paid to short wave propagation.

(a) Numerical set-up

At every time step the model provides the updated dependent variables of the governing equations. The time step is automatically adjusted during the computation to fulfil the stability conditions. Free surface elevation at gauge locations is sampled with a 125 Hz frequency. The frequency output for dependent variables at every cell location is set to 25 Hz. Since a high spectral resolution is required for frequency domain analysis, each numerical simulation is 300 s long, the same as the laboratory experiment. The computational mesh is 21.6 m long and 0.6 m high. This extends from x=−19.6 m up to x=2 m with respect to the reference system introduced in §2. Therefore, the zero position of the ‘numerical’ wave maker lies at 1 m from the start of the mesh.

The dependence of the numerical accuracy on the spatial grid resolution is studied through a mesh sensitivity analysis. A uniform grid system in the vertical direction is used with Δz=0.5 cm. In the horizontal direction, the grid cell size Δx is allowed to vary between 0.8 cm in the swash zone and 1.2 cm close to the wave maker. The total number of cells in the numerical domain is 2203×121. The order of magnitude of the time step lies between 10−4 and 10−2 s in the simulated cases. The computational time is in the order of 30 h in a standard PC.

(b) Analysis techniques

Once the time series of surface elevations from laboratory experiments and model simulations at every gauge locations are obtained, a resampling procedure is needed. Since the two data sources considered have different sampling frequencies, an interpolation procedure is used to homogenize them with the same frequency equal to 25 Hz. Distribution of wave energy within the frequency domain as well as the spatial evolution can be visualized by an amplitude spectra plotted at different locations. Amplitude spectra are obtained via a fast Fourier transform of the surface elevation data series. The surface elevation record is separated into long wave and short wave components, denoted by ηlf and ηhf, respectively. Long wave components are obtained from low-pass filtered surface elevation time series with a cut-off frequency of fc/2. Then, a high-pass filter with the same cut-off frequency is applied to the surface elevation series in order to obtain short wave components. For the purpose of describing spatial and time evolution of wave groups, the short wave envelope is determined by means of the Hilbert transform operator of the surface elevation data: Embedded Image 4.1 in which Embedded Image denotes the Hilbert transform operator and ∥lf represents the low-pass filter operator. Cross correlation analysis between two time series is a commonly applied technique to find out their relationship. The normalized cross-correlation function between two time series a(t) and b(t) is determined as Embedded Image 4.2 where 〈〉 denotes time averaging over the time series length, τ is the time shift, σa and σb are the standard deviations of the two time series, and −1≤Rab≤1 owing to the normalization procedure.

(c) Free surface displacement

Figure 3 shows the measured (solid line) and predicted (dashed line) free surface elevation time series at different cross shore locations for case B2. Growth of wave amplitude due to the shoaling process is observed up to gauge G9 where breaking of the two largest waves in the group occurs. The free surface elevation shape becomes more peaked and asymmetric as waves approach breaking. Groupiness appears to be destroyed in the surf zone where broken waves show a conoidal shape which tends to a saw-tooth shape in the inner surf zone (gauge G13). The model seems to be able to reproduce at every location the propagation features of all waves present in the time series. The changes in amplitude and in shape owing to the shoaling process are reproduced satisfactorily by the model. Breakpoint location as well as wave amplitude damping within the surf zone is correctly predicted. Moreover, wave celerity along the flume appears to be well captured.

Figure 3.

Free surface elevation time series, case B2. Solid line: laboratory observations; dashed line: model predictions.

Figure 4 shows the measured (thick solid line) and predicted (thick dashed line) low-frequency component of the free surface time series at different cross shore locations for case B2. In order to simplify the interpretation of the figure, the measured (thin solid line) and predicted (thin dashed line) total free surface elevation are also illustrated, and the scale of the low-frequency motion is magnified ×10. At gauge G1, the bound long wave trough appears under the peak of the short wave group, consistent with the constant depth solution of Longuet-Higgins & Stewart (1962). A smooth positive forward long wave is observable ahead of the group. Further shoreward, the dynamic set-down is delayed with respect to the short wave group. The time delay starts in the shoaling zone over the sloping bottom (gauge G4) and increases as the waves approach the shore. Both the long wave trough and the positive surge in advance of the group are strongly enhanced, reaching their maximum values inside the surf zone (gauge G13). The model appears to be capable of capturing the long wave evolution within the nearshore region. The amplitude and phase of incident long waves, which are the forward long wave and long wave trough behind the group, are correctly reproduced. Only a slight sub estimation of the long wave trough at gauge G10 is noticeable.

Figure 4.

Long wave elevation, case B2. Thin solid line, total free surface elevation measured; thin dashed line, total free surface elevation calculated; thick solid line, long wave elevation measured (right scale); thick dashed line, long wave elevation calculated (right scale).

Besides incident long waves, in figure 4 long waves propagating offshore can be recognized. Focusing the attention on gauge G4, the positive pulse at t≈112 s is identified as a long wave propagating offshore owing to the reflection of the positive forward long wave in front of the group at the shoreline (x=0 m). This long wave crest is preceded and followed by weaker long waves with negative elevation. The long wave trough centred at t≈114.5 s is the dynamic set-down reflected from the shoreline. The negative pulse at t≈108 s is identified as the long wave radiated from the breakpoint. In the model results, the amplitude of the long wave crest reflected from the shoreline is slightly smaller than in the laboratory data. A possible explanation is that in the model, the beach is more reflective with a partial reflection of the positive surge noted before the shoreline (x<0 m). This hypothesis is endorsed by the fact that the amplitude of the long positive pulse radiated from the breaking zone is slightly higher than that measured in the laboratory (gauge G4, the crest at t≈105 s). That is, part of the energy of the incident long wave crest is reflected before the shoreline. As a consequence, the run-up predicted by the model is slightly lower than that observed during the experiments. In general terms, the model seems to be able to correctly reproduce the breakpoint generation and propagation of outgoing long waves.

The relative errors are calculated as the difference between predicted and experimental values, divided by the experimental value: Embedded Image 4.3 The errors are expressed in percentage; the modulus operator returns only positive values. The variables chosen for the model validation are the root mean square height, Hrms, of the focused group waves and the incident bound long wave trough. For case A1, characterized by a small wave group amplitude (A=0.03 m) and a high-wave frequency (fc=0.75 Hz), the group-induced long waves are significantly smaller than in the other cases. Due to the scarce relevance of long waves on the total energy balance, case A1 is not accounted for in the validation analysis. Hence, the comparison is conducted for 65 values of each variable (13 gauges × 5 cases). Figure 5 shows the comparison between predicted and measured values including all the gauges both inside and outside the surf zone. Figure 5a illustrates the comparison of root mean square wave height calculated from the Nw−2 largest waves of the focused group, Nw being the number of waves per group (Nw=7 and 9 waves for the cases with fc=0.583 and 0.75 Hz, respectively). The predicted Hrms presents a mean relative error of 5.11 per cent. Figure 5b shows the comparison between predicted and measured values of the incident bound long wave. The predicted value of the dynamic set-down presents a mean relative error of 9.20 per cent.

Figure 5.

(a) Observed versus predicted Hrms and (b) long bound wave. Solid line, perfect agreement; dashed line, +10% error bound.

In general terms, the error of the long wave component is higher than that of the short wave. The main reason is because the ratio between long wave amplitude and cell size is really small. The dynamic set-down amplitude is in the order of 1 cm (≈ 2 cells), while the short wave amplitude is in the order of 10 cm (≈ 20 cells). Despite that, the mean relative error of the long wave component is lower than 10 per cent, which represents a satisfactory result, attesting the model ability to study infragravity wave transformation over a sloping bottom.

(d) Amplitude spectra

Figure 6 shows the surface elevation amplitude spectra at different cross-shore locations for case C1. The locations selected here correspond to gauges G1, G6, G10 and G13. The amplitude spectra obtained from the free surface elevation data of the numerical model (blue dashed line) and laboratory experiments (red solid line) are displayed. Moreover, amplitude spectra of free surface elevation predicted at the wave paddle location by second-order (black line) and first-order theory (green line) are also represented in figure 6. These theoretical spectra are represented at gauges G1 and G6, G10 and G13 in order to emphasize the spatial redistribution of wave energy during the shoaling and breaking processes. Amplitudes are presented on a logarithmic scale in order to visualize in the same figure infragravity and primary bands that cover a large range of values.

Figure 6.

Amplitude spectra, case C1. Red line, observed; blue line, predicted by the model; green line, predicted by linear theory at wave paddle location; black line, predicted by second-order theory at wave paddle location. (a) Gauge 1; (b) gauge 6; (c) gauge 10; (d) gauge 13.

At the most seaward location, the ‘top hat’ spectral shape of measured data for the frequency range corresponding to primary short wave components is properly reproduced by the model. The predicted amplitude of the first superharmonic (in the frequency range 1.4<f<1.6 Hz for case C1) agrees well with laboratory data. Furthermore, long wave components are satisfactorily reproduced by the model and the additional energy not predicted by theory is attributable to long wave reflection at the shoreline and wave maker. At shoreward locations, the original ‘top hat’ shape is modified and energy is redistributed between wave components. The energy is transferred into components at frequencies slightly higher than the upper limit of the primary wave band, consistent with the transient wave group behaviour described by Rapp & Melville (1990). Moreover, amplitude increases in the superharmonic band owing to wave shoaling is clearly visible (gauge G6). Inside the surf zone (gauges G10 and G13), wave breaking and consequently enhanced turbulence cause a decrease of primary components as well as superharmonics. Energy dissipation of the primary components is stronger for the frequencies close to the upper limit of the primary wave band. Low frequency amplitudes increase across the whole nearshore region as a result of shoaling of long waves and energy transfer from higher frequency components. Inside the surf zone, the energy concentrated in low-frequency components represents a significant part of the total wave energy. Summarizing, we can say that at every location and within the entire frequency domain, amplitudes from numerical results are in very good agreement with those obtained from laboratory data. Only for the punctual frequency of 0.04 Hz (=25 s) there is a sensible mismatch at gauges G1, G10 and G13. At this frequency, which corresponds to the oscillation period of the flume, numerical results are higher than laboratory data. So, while the seiche is dissipated in the real flume by friction forces during the experiment, in the numerical flume the dissipation is less effective. This can be a result of neglecting friction forces on the side-walls of the flume by a 2DV model, while the model employs a no-slip boundary condition on the bottom. Since gauge G6 is located close to the centre of the flume close to the node of the seiche, no amplitude peak is observed at f=0.04 Hz.

5. Wave transformation analysis

In the validation procedure, the active absorption at the wave maker has not been implemented in order to correctly reproduce the laboratory conditions. In this section, the active absorption at the wave maker is turned on. This allows the outgoing long waves to leave the computational domain avoiding the interaction with the long wave reflected at the wave maker.

The analysis techniques used in this section are the same as those described in §4. Since the spatial laboratory sampling does not enable an exhaustive plot construction in the 〈τ,x〉 or 〈t,x〉 planes, the very high spatial resolution numerical data are used to overcome this limitation. Therefore, contour plots of cross correlation function in the 〈τ,x〉 plane as well as contour plots of surface elevation or short wave envelope in the 〈t,x〉 plane are constructed from numerical data.

(a) Free surface and short wave group evolution

The contour plot of surface elevation for case C1 is illustrated in figure 7. Note that owing to the colour scale chosen, the surface elevations lower than 0 m are coloured dark. This allows a direct observation of the long wave trough reflected from the shore at t≈110 s. In the upper portion of the figure the run-up is specified. The propagation path defined by the group velocity calculated from linear theory using the central frequency, fc, is also illustrated by a white dashed line, its mathematical expression being: Embedded Image 5.1 The values of x0 and t0 are chosen in order to superpose the dashed line over the propagation path of the short wave envelope peak. Breaking of the largest waves in the group occurs at x≈−6.5 m, just before the horizontal bottom region with 0.2 m depth. It is possible to note that the celerity of the short waves ahead of the group is higher than predicted by linear theory in the inner surf zone for x≥−1.5 m (linear theory predicts that in shallow water wave celerity tends to group celerity indicated by the white dashed line). On the contrary, waves that reach the inner surf zone, corresponding to the group peak predicted path, interact with the long wave trough and propagate slower than expected from theory. Swash zone motion appears to be driven by two short wave bores leading the short wave group and enhanced by the long wave crest, a large run-up event is generated.

Figure 7.

Space–time evolution of surface elevation, case C1. Dashed line illustrates the group path propagation predicted by linear theory. Horizontal lines indicate the bottom slope changes. Solid line in the upper portion indicates the run-up on the beach face. Elevation in m.

The high spatial resolution of the numerical results enables the construction of the short wave envelope in the 〈t,x〉 plane. Figure 8 shows the contour plot of the space and time evolution of the short wave envelope A for case A2 (figure 8a) and C1 (figure 8b). The velocity propagation appears to be equal to the group velocity as predicted by linear theory. A strong decrease in peak amplitude owing to the short wave energy dissipation is observed at the breakpoint. For case C1, the collapse of the envelope peak leads to the formation of two separated little peaks that propagate ahead of and behind the original envelope ridge, pointing out a weak groupiness inversion in the surf zone for x>−5 m. Total dissipation of short wave energy in the inner surf zone cancels any leftover of wave envelope; no reflection is observed at the shoreline.

Figure 8.

Space–time evolution of short wave envelope. (a) Case A2 and (b) case C1. Dashed line illustrates the group path propagation predicted by linear theory. Horizontal lines indicate the bottom slope changes. Elevation in m.

(b) Evolution of the low-frequency component

The evolution in the spatial and time domain of the low-frequency component is shown through a contour plot in figure 9. The position where a slope change occurs is indicated in the figure by a black horizontal line. The gauge location is similarly indicated by horizontal numbered white lines. The contour plot extends seaward from x=0 m. On the upper portion the low-frequency component of run-up along beach face is displayed. The shoreward propagation path of the long bound wave predicted by linear theory using group velocity of the central primary wave (fc) is indicated by the white dashed line. In the same way, seaward propagation of the reflected long wave calculated with shallow water wave celerity (Embedded Image) is shown. The mathematical expression of the reflected predicted path is: Embedded Image 5.2 For case A2 (figure 9), the long wave trough appears to propagate with group celerity until strong shoaling and breaking of short waves take place over the sloping bottom. Then a time shift is observed and in the surf zone the bound long wave celerity appears slower than the theoretical group velocity. The long wave trough is amplified by the increasing radiation stress gradient in the shoaling zone and the intensified forcing owing to decreasing depth over the beach profile, reaching a minimum value inside the surf zone. A sudden growth of the long wave crest preceding the long wave trough is observed just offshore of the breakpoint (x=−2 m). This positive surge drives a low frequency run-up on the beach face, then a positive pulse radiates offshore. Interaction between the positive pulse owing to the reflection of the long wave crest at the shoreline and the incident long wave trough generates a transient standing wave node at x≈−1.5 m that lasts approximately 5 s. Seaward propagation of the long wave reflected at shoreline seems to be well described by linear shallow water wave theory, pointing out that outgoing waves behave as free waves. The positive pulse owing to the reflection of the long wave crest is followed by a negative pulse owing to the reflection of the bound long wave trough at the shoreline. The amplitude of the reflected negative pulse is smaller than that of the incident bound wave trough in the inner surf zone. The process of reduction of the negative pulse appears to be linked to the superposition of long waves generated within the surf and the swash zone (Baldock 2006) rather than to the dissipation owing to bottom friction (Henderson & Bowen 2002). Therefore, the presence of the outgoing positive pulse which propagates from the shoreline at t≈115 s can be explained as a low-frequency wave generated within the surf and swash zone.

Figure 9.

Evolution of low-frequency surface elevation, case A2. Dashed line illustrates the group path propagation predicted by linear theory. Black horizontal lines indicate the bottom slope changes. White horizontal lines show the location of wave gauges. Solid line in the upper portion indicates the low-frequency component of the run-up on the beach face. Elevation in m.

Additionally, a long wave trough is observed to lead the positive pulse owing to the reflection of the long wave crest at the shoreline. This negative pulse appears to be generated from the breakpoint and it propagates offshore with shallow water wave celerity. Direct radiation of a long wave trough from the breakpoint was observed by Baldock (2006) and is consistent with theoretical models (Symonds et al. 1982; Schäffer 1993). Moreover, a weak positive pulse propagating offshore leads the negative pulse from the breakpoint. This outgoing long wave crest can be considered as a free wave owing to the reflection of the forward positive surge over the sloping bottom for x≥−4 m, enhanced by a long wave radiated during the amplification of radiation stress forcing as the short wave group shoals in intermediate depths.

(c) Cross-correlation analysis

The relation between the short wave envelope and long wave motion is investigated through a cross correlation analysis of the numerical results. Figure 10 illustrates the cross correlation function RηA between low-frequency motion and squared short wave envelope for case B2, measured at the same location. The white dashed line indicates the travel time required for the groups of waves to reach the shore with group velocity and for long waves to travel back to the specified location with shallow water wave celerity: Embedded Image 5.3 A strong negative correlation at zero time lag is observed near the wave generation region (x<−14 m). This is consistent with bound long wave trough propagating under the peak of the short wave envelope. However, as the wave group approaches the shore, the time lag tends to increase assuming small positive values indicating a time delay of the low-frequency motion with respect to the short wave group, as observed by Janssen et al. (2003). The positive correlation at τ≈−6 s indicates the presence of a positive surge shoreward propagating in front of the short wave group. Outside the surf zone, the positive surge appears to travel at the group celerity suggesting that it is driven by local forcing of the short waves. Due to the increasing long wave crest during shoaling the positive correlation also increases towards the shore until the breakpoint is reached. Then a small discontinuity between gauges G8 and G9 is observed. This can be ascribed to the sudden shoaling and the subsequent breaking of the greater waves in the group that produces a displacement of the envelope peak ahead of the group. Moving further towards decreasing depth, the positive correlation ridge experiences a strong discontinuity and the time lag changes from τ≈−6 to 0 s very close to the shoreline (x≥−1.5 m). The positive correlation between long wave motion and short wave group at near zero time-lag is owing to the fact that in very shallow water low-frequency motion strongly influences the propagation pattern of short wave groups. Short waves above the crest of the long wave experience a shoreward current and propagate in deeper water than short waves just in front of them: a convergence of wave bores can take place.

Figure 10.

Cross correlation between low-frequency motion and |A|2, case B2. Black horizontal lines indicate the bottom slope changes. White horizontal lines indicate the location of wave gauges. White vertical line indicates the zero time lag and dashed white line indicates summed travel times determined from the wave celerity.

A weak negative correlation is observed along the white dashed line. This indicates that the bound wave trough propagating with the wave group is released as the short waves break and, after reflection at the shoreline, propagates offshore with shallow water wave celerity. Inside the surf zone (for x>−5 m), a possible superposition of long waves generated inside the surf and the swash zone in conjunction with the deformation of the group structure owing to the short wave breaking cause a negligible correlation along the white dashed line. The positive correlation due to the reflection of the forward positive surge leads the negative correlation path of the long wave trough already commented. A long wave trough preceded by a long wave crest appears to be radiated directly seaward from the breaking zone.

(d) Radiated free waves

Here, the numerical model is applied for the study of the bottom slope and the short wave steepness influence on the radiation of free long waves from the breakpoint. New simulations are conducted for the cases with a relatively low central frequency (fc=0.583 Hz). Hence, cases A2, B2 and C2 are simulated keeping the same wave characteristics already described but changing the beach slope. Three different profiles with a uniform bottom slope, hx, equal to 1:20, 1:30 and 1:40 are analysed. The toe of the slope is kept at 8.6 m from the zero position of the wave maker. A mesh grid with the same spatial resolution of the grid described in §4 is introduced. The grid is uniform in the vertical direction with Δz=0.5 cm. In the horizontal direction, grid cell size Δx varies between 0.8 cm in the swash zone and 1.2 cm close to the wave paddle. New simulated cases with the corresponding normalized bed slope parameter are summarized in table 4.

View this table:
Table 4.

New simulated cases.

The normalized bed slope parameter, β, was introduced by Battjes et al. (2004) to differentiate between the mild slope regime and the steep slope regime: Embedded Image 5.4 where h and hx are representative values of depth and bottom slope in the shoaling zone; k and ω are the long wave number and angular frequency. It is worth mentioning that the parameter β is closely related to the parameter previously introduced by Baldock et al. (2000) and Baldock & Huntley (2002) who examined the influence that the relative magnitude of the breakpoint excursion induces on the breakpoint forcing mechanism. Taking hx(eq)=1:25 and h=0.3 m as representative bed slope and depth in the shoaling zone, we find β≈0.44 indicating a steep slope regime (we refer to an equivalent bottom slope, hx(eq), to indicate that the slope is not constant along the profile). The new simulations conducted changing the beach slope bring values of β between 0.55 (for hx=1:20) and 0.27 (for hx=1:40).

The presence of long waves radiated directly from the breakpoint, as noted above and as observed by Baldock (2006), is consistent with theoretical models (Symonds et al. 1982; Schäffer 1993). Using the linear, constant depth solution, Nielsen et al. (2008) suggested that a sudden growth of the radiation stress forcing in the shoaling zone can generate a positive pulse propagating ahead of the group. Nielsen & Baldock (2010) identified the strong long wave crest outside the surf zone observed in previous experiments (Baldock 2006) as a free wave propagating onshore. Since the radiation stress growth occurs in the shoaling zone over a sloping bottom, the constant depth solution remains only qualitative.

The negative pulse radiated offshore from the breakpoint, observed in the present experiments and reproduced by the numerical model, can be explained as a breakpoint generated free wave. In agreement with previous experiments of Baldock et al. (2000) and Baldock & Huntley (2002), Battjes et al. (2004) noted that the relative importance of the breakpoint generated waves over the incident bound waves increases with increasing bed slope. They stated that breakpoint generation is ineffective in the mild slope regime, for say βb<0.3 (βb is calculated with depth at breaking). The amplitude of radiated long waves, for weakly modulated bichromatic wave groups, varies owing to the changing position of the breakpoint. For a transient focused group with the breakpoint located in intermediate depths, this phenomenon seems to be strictly related to the breaking type and then to the rate of decrease of radiation stress forcing in the surf zone. As short waves break in intermediate depths, the long wave released from the breakpoint is not a shallow water wave and a transformation of the propagation pattern (in essence a reduction in the incident trough amplitude and the generation of free long waves) is needed as the radiation stress forcing decays. So, it appears that the short wave breaking type and the radiation stress dissipation rate in the surf zone strongly influence the long wave released at the breakpoint (Baldock & O’Hare 2004). It is worth mentioning that the bed slope parameter does not bring information on the local characteristics of short waves. Since the breaking type is mainly controlled by the bed slope and the short wave steepness, the importance of the breakpoint generated free wave is expected to increase with the Iribarren number (Embedded Image).

The relation between the breakpoint generated outgoing wave and the incoming bound wave at gauge G1 location is investigated using the simulations conducted with different beach slopes. For each simulation, gauge G1 lies at the same distance from the wave paddle and thus the distance from the origin of the horizontal co-ordinate increases for decreasing beach slope. Figure 11 shows the low-frequency free surface elevation time series at gauge G1. The negative pulse breakpoint generated is observable at t≈108,112 and 115 s in figure 11ac, respectively. No breakpoint generated wave is noted in figure 11d.

Figure 11.

Long wave elevation at gauge G1. solid line, total free surface elevation; thick line, long wave elevation (right scale). (a) case A2-1, hx=1:20, (b) case A2, hx(eq)=1:25, (c) case A2-2, hx=1:30, (d) case A2-3, hx=1:40.

The ratio (γ=Arad/Ainc) between the breakpoint generated outgoing wave and the incoming bound wave both at gauge G1 is related to both the Iribarren number, Ir, and the normalized bed slope parameter, β, in table 5. From the table, it can be noted that for a given β the relative importance of the outgoing radiated wave, γ, increases with Ir and as the wave steepness decreases. Furthermore, for decreasing beach slope, and thus decreasing β and Ir, the relative importance of the outgoing radiated wave also seems to decrease. Finally, for β=0.27 no significant breakpoint generated wave is detected. This is consistent with Battjes et al. (2004) who stated that breakpoint generation is ineffective under a mild slope regime.

View this table:
Table 5.

Relative importance of breakpoint generated waves.

(e) Influence of a flat bottom

Finally, we are mainly interested in the analysis of the effects that the presence of a large flat bottom induces on the propagation pattern of long waves. In particular, attention is focused on the propagation of the long wave crest in front of the group and the subsequent run-up produced on the beach face. The flat bottom region is introduced since it resembles a morphological feature such as a submerged marine terrace, often present on continental shelves. A new simulation, called A2-flat, is conducted for case A2 in which the bottom profile is obtained using the same geometry of the laboratory but extending the 0.2 m deep horizontal bottom section from 2 to 12 m. The mesh grid resolution is similar to that described previously: the grid is uniform in the vertical direction with Δz=0.5 cm, while in the horizontal direction Δx varies between 0.8 cm in the swash zone and 1.2 cm close to the wave maker. We choose the wave condition of case A2 since breaking occurs at x≈−2 m and therefore the flat bottom part lies in the shoaling zone.

In figure 12, the time history of low-frequency free surface is shown at gauge G13 for the simulation A2-1, A2 and A2-flat. We recall that for these simulations the wave characteristics at the wave maker are the same, and the bottom profile differs only in the flat bottom part of 0.2 m depth (see figure 13 where the upper part of the profile is drawn). For each simulation, gauge G13 is kept at the same distance from the shoreline (x=−2.91 m). As can be seen in figure 12 a slight increase of the incident bound long wave trough takes place. More evident is the strong growth of the positive surge in front of the group. At gauge G13 the positive surge amplitude is 0.8, 1.4 and 3.4 mm in figure 12ac, respectively. The strong increase of the long wave crest in front of the group can be explained considering this positive surge constituted by a superposition of the positive part of the bound long wave (Baldock 2006) and free waves radiated ahead of the group as short waves shoal in intermediate depths (Nielsen & Baldock 2010).

Figure 12.

Long wave elevation at gauge G13 (x=−2.91 m). Thin line, total free surface elevation; thick line, long wave elevation (right scale). (a) Case A2-1, hx=1:20, (b) case A2, hx(eq)=1:25, (c) case A2-flat, hx(eq)=1:50.

Figure 13.

Spatial evolution of long wave elevation. Solid line, total free surface elevation; grey dashed line, wave envelope; black dashed line, long wave elevation (right scale). (a) Case A2-1, hx=1:20, t=96.7 s; (b) case A2, hx(eq)=1:25, t=98.4 s; (c) case A2-flat, hx(eq)=1:50 t=107 s.

Figure 13 illustrates the spatial evolution of the low-frequency component and the short wave envelope taken at the time instant when maximum amplitude positive surge occurs at gauge G13 (see figure 12). For case A2-flat (figure 13c), the entire group structure is propagating above the constant depth (0.2 m deep) region. As a consequence of the smaller depth, all the short waves in the group have already shoaled and are steeper than in the case with constant slope (figure 13a). The marked radiation stress gradient in conjunction with the intensified forcing owing to decreasing depth enhances the long bound wave crest in front of the group. Furthermore, free long waves radiated ahead of the group from the shoaling zone can play an important role in determining the positive surge shape. It turns out that the distance from the toe of the beach, which corresponds to the point where the shoaling process starts, to gauge G13 is larger for the case with a flat bottom. Along this distance, it is expected that free long waves radiated ahead of the group during the shoaling of short waves are able to partially separate from the bound wave trough. The long wave development is the result of the time Tf spent by the group of length Lg over the flat bottom in comparison to the time Ts required for complete separation (Nielsen 2009). With the flat bottom geometry, Tf=9.8 s and Embedded Image s; this means that although the 12 m flat section appears sufficiently long to allow a significant positive surge enhancement by free waves, it is not long enough to obtain a complete separation from the long wave trough.

In figure 13a the envelope appears almost in anti-phase with the low-frequency motion indicating that local forcing governs the long wave propagation. While, a significant phase lag between the envelope peak and the long wave trough is observed in figure 13c. This suggests that, with a long flat bed, the long wave shape is strongly affected by free waves which play a prominent role in the enhancement of the total positive surge observed at gauge G13.

Since short waves are allowed to propagate in very shallow water on the crest of the low-frequency motion (Janssen et al. 2003; note the positive correlation for x≥−1.5 m at near zero time lag observed in figure 10 for case B2), the strong long wave crest observed with a flat bottom is expected to enhance the run-up along the beach face. In fact, the maximum run-up predicted by the model is 2.5 cm for the case A2-flat, and 1.8 and 1.5 cm for cases A2 and A2-1, respectively (not shown). Therefore, with the same wave characteristics in the open sea, the flat bottom appears to be responsible for the generation of a significant increase of the run-up on the beach face which for the simulated cases is of 67 per cent.

6. Conclusions

In this work, cross shore propagation of long waves induced by a transient short wave group has been studied by means of numerical modelling. Dataset obtained by a new laboratory experiments have been used to validate the RANS type model IH-2VOF. Short wave groups with different central frequency and short wave steepness have been generated and focused at the breakpoint of the short waves. Free surface elevation data have been measured simultaneously at different cross shore locations both outside and inside the surf zone.

Model validation has been conducted to test the ability of a RANS type model IH-2VOF to study low-frequency processes in the nearshore zone. A great effort has been paid to the generation and absorption procedure. Numerical wave generation has been carried out simulating the wave maker movement by means of a moving body algorithm. Comparisons between numerical and laboratory free surface time history have proven the high accuracy of the model in reproducing both gravity and infragravity wave transformation. The model is able to reproduce shoaling and breaking of short waves. Moreover, the enhancement of long waves and the energy transfer to low-frequency motion is well addressed by the model.

In the steep slope regime, a gentle positive pulse radiated directly from the shoaling region, owing to a superposition of a partial reflection of the long wave crest in front of the group and outgoing free long waves, is noted. Moreover, a long wave trough radiated directly from the breakpoint is detected in the laboratory and predicted by the model. By changing the bottom profile, the model is applied to study the influence of the Iribarren number and the bed slope parameter β (Battjes et al. 2004) on the breakpoint long wave generation. For a given β, an increase of relative importance of the breakpoint generated long wave over the incident bound long wave as a result of the increase of the Iribarren number is observed. For gentle slope, no significant radiated long wave is detected. This is consistent with Battjes et al. (2004). So, especially for transient wave groups that break in intermediate depths, it appears that the outgoing long wave released at the breakpoint is controlled by both the bed slope parameter and the Iribarren number.

Finally, a new simulation is conducted to study the influence of a large flat bottom region on the propagation pattern of incident long waves. It is observed that, as the group propagates above the flat bottom located outside the surf zone, a strong long wave crest in front of the group is generated. The flat bottom region seems to be responsible for the development of a strong radiation stress forcing and at the same time enables a partial separation ahead of the group of incident free long waves radiated from the shoaling region. Furthermore, the marked long wave crest observed with the flat bottom significantly enhances the subsequent run-up on the beach face.

  • Received June 23, 2010.
  • Accepted October 14, 2010.


View Abstract