Irregular wave runup statistics on plane beaches: application of a Boussinesq-type model incorporating a generating-absorbing sponge layer and second order wave generation

Abstract Efficient absorption of reflected waves at the offshore boundary is a prerequisite for the accurate physical or theoretical modelling of long-duration irregular wave runup statistics at uniform, gently sloped beaches. This paper presents an implementation of the method suggested by Zhang et al. (2014) to achieve reflected wave absorption and simultaneous generation and propagation of incident waves in an existing numerical wave flume incorporating a moving boundary wavemaker. A generating–absorbing layer is incorporated within this 1DH hybrid Boussinesq-nonlinear shallow water equation model such that inshore-travelling incident waves, encompassing bound-wave structure approximately correct to second order, propagate unhindered while offshore-travelling reflected waves are absorbed. Once validated, the method is used to compile random wave runup statistics on uniform beach slopes broadly representative of dissipative, intermediate, and reflective beaches. Analyses of the individual runup time series, ensemble statistics and comparison to an empirical formula based on experimental runup data suggest that the main aspects of runup observed in the field are properly represented by the model. Existence of an upper limit on maximum runup is investigated using a simple extreme-value statistical analysis. Spectral saturation is examined by considering ensemble-averaged swash spectra for three representative beach slopes subject to incident waves with two different offshore significant wave heights. All spectra show f −4 roll-off at high frequencies in agreement with many previous field studies. The effect is also investigated of the swash motions preceding one particular extreme runup event on the eventual maximum runup elevation.

Abstract: Efficient absorption of reflected waves at the offshore boundary is a prerequisite for the accurate physical or theoretical modelling of long-duration irregular wave runup statistics concerning uniform, gently sloped beaches. This paper presents an implementation of the method, suggested by Zhang et al. (2014), to achieve reflected wave absorption and simultaneous generation and propagation of incident waves for an existing numerical wave flume incorporating a moving boundary wavemaker. A generating-absorbing layer is incorporated within this 1DH hybrid Boussinesq-nonlinear shallow water equation model such that inshore-travelling incident waves, encompassing bound wave structure correct to second order, propagate unhindered while offshore-travelling reflected waves are absorbed. Once validated, the method is used to compile random wave runup statistics on uniform beach slopes broadly representative of dissipative, intermediate, and reflective beaches. Analysis of the individual runup time series, ensemble statistics and comparison to an empirical formula based on experimental runup data suggests that the main aspects of runup observed in the field are properly represented by the model. Existence of an upper limit on maximum runup is investigated using a simple extremevalue statistical analysis. Spectral saturation is examined by considering ensemble-averaged swash spectra for three representative beach slopes subject to incident waves with two different offshore significant wave heights. All spectra show $f^{-4}$ roll off at high frequencies in agreement with many previous field studies. The effect of the swash motions preceding one particular extreme runup event on the eventual maximum runup elevation is also investigated.
The problem is that the model setup falls between two stools. On one hand, with the idea to mimic the circumstances of a physical wave flume, it includes weakly nonlinear wave generation based on physical wavemaker theory.
On the other hand, it also includes a wave absorption/generation sponge layer that has no equivalent in physical experiments and thus annihilates the rationale behind the use of physical wavemaker theory.
Nonlinear wavemaker theory basically addresses the Lagrangian nature of a moving paddle and the mismatch between paddle shape and wave velocity profile. While these issues are inevitable for physical wavemakers, they need not exist in Boussinesq models for which simpler and better wave generation alternatives are available. So why stick to generation using a wavemaker theory of questionable relevance, if comparison with physical experiments has been given up (and made impossible by introducing a sponge that cannot be physically realized)? A simple fixed-position boundary condition should suffice, especially if there is a generating-absorbing sponge layer available for fine tuning the incident field to the desired one. This would also allow the sponge to extend to the boundary for better performance.
It may be so that the precise nature of the wave generation at the boundary is not important due to corrections taking place in the generating-absorbing sponge. In that case there should be no need to change the model setup. However, the paper would then need to clearly identify the selfcontradictory nature of the setup, explain the historical reason for arriving at this setup and argue why it does not compromise the results. Response: see paragraph 3 (line 29 onward) on page 11 and the first line on page 12 of revised manuscript. We argue that the second-order generation achieved with the moving boundary wavemaker is of sufficient accuracy to justify inclusion. The performance of the compact sponge layer suggests our results will not be compromised. We do, however, point out the self-contradictory nature of the current implementation.
Despite the lengthy presentation, the section on sponge layer absorption provides only a vague impression of how well the sponge layer performs. A fundamental property of a wall-adjacent sponge layer is the associated reflection coefficient vs frequency as can be obtained by a single dedicated test run for a broad-band spectrum using reflection analysis. Response: see last paragraph (line 47) on page 13 extending to page 14 (lines 21 -41) and Figure 4. Note, we use a regular waves to study the wall adjacent sponge layer -reflections of outgoing long waves are likely to compromise an full repeat period irregular wave incidence run.
In general, the manuscript is quite wordy and would benefit from being shortened. Section 2 is reduced significantly in size and section 1 has also been shortened. More figures are present in the latest document with new material and so the manuscript's length remains the same.
Specific points: 1. Consider rephrasing the second bullet of the highlights -have tried to improve it. 2. Page 6, line 43: "approximant of the celerity" -check if it is not the squared celerity -Document shortened and discussion on page 6 removed Detailed Response to Reviewers 3. Bottom of page 6: Much has happened since that newest listed reference from 1998, see e.g. publications by Madsen et al. and Lynett et al. -Document shortened and discussion on page 6 removed 4. Page 8, equation (6) appears to be inconsistent with second order wavemaker theorydiscussion of second order wavemaker approach is now very brief in Section 2 (reference to the original work of Orszaghova is provided, at the end of paragraph 2 of Section 2). 5. Page 10, line 25-26: frequency dispersion is neglected, not just may be -see line 23 on page 6 of manuscript. 6. Page 10, line 31-36: How does this relate to equation (6) amplitudes are also taken to be random -all references to reference phase method have been corrected, now refer to random-phase/random-amplitude method. 13. Page 21, transition to page 22: fmax=2.08fp is a severe truncation of a PM spectrum. Specify if moments were computed from the full or the truncated spectrum -line 44 in first paragraph in 4.1.2, page 19 14. Page 23, line 19: Perhaps replace "submerged" by "wetted" -removed this sentence 15. Page 23, line 20: Therefore? Removed this sentence (repetition) 16. Page 24, line 6, we we -fixed 17. Page 28. Absorption and second order effects should be taken out separately to isolate their respective influence -I'd argue that this is a simple illustration of how different the experimentally generated wave fields may be from those generated in the numerical simulations. We've already isolated the effect of reflections and literature exists regarding the importance of second order corrections on runup -I don't think it's necessary to labour the point with further analysis. 18. Page 29, line 9 and on: It is always nice to know the consequences of an arbitrary choice and even without making new model runs it should be a fairly quick exercise to conduct an approximate analysis making powers-of-two changes to the record length by simple time series splitting or merging. See Figure 10, 11 and pages 28 -29 19. Page 30, line 37: "It is likely…", you could try it out (see previous point) -see previous response. 20. Page 34, line 2: I was actually surprised to see how small this influence is! Highlights of "Irregular wave runup statistics on plane beaches: application of a Boussinesq-type model incorporating a generating-absorbing sponge layer and second order wave generation" by Colm J. Fitzgerald, Paul H. Taylor, Jana Orszaghova, Alistair G. L. Borthwick, Colin Whittaker and Alison C. Raby  Simultaneous absorption of offshore-travelling reflected waves and propagation of inshoretravelling irregular incident waves on plane beaches is achieved by incorporating a generating-absorbing sponge layer in a 1DH hybrid Boussinesq-nonlinear shallow water equation  Irregular wave runup with outgoing wave absorption is simulated in a two-stage process at each time step: incident wave-fields encompassing bound harmonic structure to second order are computed in the incident wave propagation stage; this incident wave field is imposed on the sponge layer located offshore of the beach in the full runup stage damping any reflected waves  Long-duration irregular wave runup statistics are collated from the shoreline position timehistories for many realisations of one incident wave energy spectrum for three different plane beach types ranging from reflective to dissipative  Spectral saturation of swash is investigated on a reflective, intermediate and dissipative beach slope using ensemble averaged runup statistics  The effect of the swash motions preceding a large runup event on the eventual maximum runup elevation is also investigated 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 Irregular wave runup statistics on plane beaches: application of a Boussinesq-type model incorporating a generating-absorbing sponge layer and second order wave generation encompassing bound wave structure correct to second order, propagate unhindered while offshoretravelling reflected waves are absorbed. Once validated, the method is used to compile random wave runup statistics on uniform beach slopes broadly representative of dissipative, intermediate, and reflective beaches. Analysis of the individual runup time series, ensemble statistics and comparison to an empirical formula based on experimental runup data suggests that the main aspects of runup observed in the field are properly represented by the model. Existence of an upper limit on maximum runup is investigated using a simple extreme-value statistical analysis. Spectral saturation is examined by considering ensemble-averaged swash spectra for three representative beach slopes subject to incident waves with two different offshore significant wave heights. All spectra show f −4 roll off at high frequencies in agreement with many previous field studies.
The effect of the swash motions preceding one particular extreme runup event on the eventual maximum runup elevation is also investigated.

Introduction
Irregular wave runup on artificial and natural beaches is of significant interest to coastal engineers, scientists, and managers concerned with flood inundation, coastal erosion and maritime spatial planning. The maximum wave runup level defines the boundary of the region of wave action and is of obvious importance in planning coastal setback distances or assessing coastal flooding probabilities. Therefore, the understanding and prediction of extreme runup events during periods of energetic irregular wave incidence on beaches is crucial for the planning and protection of new and existing coastal communities. The extent of the wave action zone also directly affects onshore sand transport, deposition and erosion. However, difficulties persist in identifying the cause and effect of near-shore phenomena due to the inherently nonlinear physics of wave hydrodynamics and continuously changing environmental conditions in the surf and swash zones.
Investigations of irregular wave runup and its extreme events have been dominated to date by field studies at natural beaches. Runup, defined here as the time-varying shoreline vertical elevation, is typically separated into the wave setup component (time-averaged shoreline elevation relative to mean water level) and the swash oscillations about this setup level. Shoreline setup typically varies over time-scales of the order of 100 s at full scale and is rather difficult to measure accurately in situ where significant variations in beach morphology may occur on similar time scales. As an instantaneous quantity, swash is generally more straightforward to measure and a significant proportion of the literature on runup focuses on swash motions. Analysis of swash measurements obtained from field studies have been reported on intermediate beaches with local foreshore slopes between approximately 1/40 and 1/20 (cf. e.g. Guza and Thornton (1982) and Raubenheimer et al. (1995)), on dissipative beaches with foreshore slopes less than 1/40 (cf. e.g. Ruessink et al. (1998) and Ruggiero et al. (2004)) and on steeper, reflective beaches with slopes of order 1/10 (cf. e.g. Raubenheimer and Guza (1996)). The saturation of swash energy spectra at wind-wave frequencies and the absence of saturation at infragravity frequencies has been observed in many of these studies, e.g. Guza and Thornton (1982)and Raubenheimer and Guza (1996), although in highly dissipative conditions it has been noted that the saturation can occur for shorter infragravity frequencies (Ruessink et al., 1998;Ruggiero et al., 2004). This  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 indicates that swash oscillations arising from incident irregular waves involve competing processes whose magnitudes may vary depending on beach slope and incident conditions.
Numerical and laboratory studies of irregular wave runup on beaches are relatively rare compared with field studies. This is primarily due to difficulties in obtaining reliable data due to the re-reflection of otherwise outgoing waves at the offshore boundary (i.e. the wavemaker in the case of laboratory studies). A significant experimental study of irregular wave runup on beaches of different slopes was reported by Mase and Iwagaki (1984), with further analysis reported by Mase (1989), involving a linear displacement-controlled wavemaker but no attempt to address the issue of reflections was made. Reports of more recent experimental campaigns on irregular wave runup are difficult to find. For beaches with gentle slopes, much of the energy from the linear components of the incident wave train is dissipated through breaking and bed friction. However, both the second-order bound long waves, which are non-negligible for energetic wave groups, and the offshore propagating long free waves arising during breaking undergo significant reflection at the offshore boundary. Despite advances in active absorption wavemaker theory (Spinneken and Swan, 2009), it is still extremely difficult to absorb reflected long waves in laboratory flumes.
Furthermore, a numerical study by Orszaghova et al. (2014) recently emphasised the importance of generating the correct long wave bound harmonic structure for runup and overtopping estimates.
Long waves reflected at the beach and re-reflected at the wavemaker could also compromise runup measurements and computations. Thus, to obtain reliable estimates for runup it is necessary to conduct tests involving relatively short irregular wave trains to minimise interference from re-reflected waves.
In this paper, we seek to simulate long duration time series of irregular wave runup events using a hybrid Boussinesq-type model which describes the pre-breaking offshore waves using the enhanced dispersion Boussinesq equations (Madsen et al., 1991;Madsen and Sørensen, 1992) and the post-breaking inshore waves by the NLSW equations. This numerical wave flume, developed by Orszaghova et al. (2012), models wave evolution in one horizontal dimension and includes an in-built piston paddle wavemaker. In order to eliminate reflection of otherwise outgoing waves at the wavemaker boundary, we have modified the numerical model to incorporate sponge layer damping terms in the governing equations. Numerous applications of sponge layer absorption zones to Boussinesq-type numerical models appear in the literature, e.g. Larsen and Dancy (1983), Wei et al. (1999) and Madsen et al. (1997a), often in combination with an internal wave generator.
Moreover, Madsen et al. (1997b) describe numerical simulations of experimental tests in an irregular wave runup study using such a configuration. Here, a generating-absorbing sponge layer 3 similar to that described by Zhang et al. (2014) (also described by Siddorn (2012) in an ocean engineering context) is introduced. Irregular wave propagation in the Boussinesq region is solved separately and imposed as a solution on a sponge layer in the full simulation of the wave runup problem. Thus, reflected waves are damped out completely while incident waves are accurately propagated inshore. The distinct advantage of this method is its ability to propagate incident waves with correct second-order bound harmonics, using the wavemaker theory of Schäffer (1996), while fully absorbing the reflected waves. Combined generation and absorption using relaxation zones rather than sponge layers has previously been proposed for Boussinesq methods by Bingham and Agnon (2005).
An advantage of investigating irregular wave runup with a numerical model, rather than through field observations, is simplification through removal of temporal and spatial uncertainties arising from variations in offshore wave conditions and beach morphology, and non-uniformity in the transverse direction. This simplification enables the basic contributing factors to wave runup to be isolated. In order to illustrate the advantages of this approach, swash spectra for large ensembles of irregular wave runup simulations are obtained for two different significant wave heights on uniform bed slopes representing reflective, intermediate and dissipative beaches.
Analysis of these swash spectra confirms observed variations in swash saturation depending on beach geometry and the influence of beach slope on the swash spectra can clearly be seen.
Saturation of runup and the nature of the extreme runup is also investigated by examining the distribution of runup maxima with particular focus on the extreme event. Furthermore, an investigation into the effects of precursor bore motions at the shore on extreme runup events is also presented utilising the sponge layer in a slightly different manner than the standard irregular wave runup approach.
The structure of the paper is as follows. Section 2 discusses the governing equations and numerical solution method of the hybrid Boussinesq-type NLSW equation model. Section 3 presents the incorporation of the absorbing-generating sponge layer in the numerical model.
Predictions of long-duration irregular wave trains incident on a beach are used to demonstrate the 'coupled' or nested simulation approach. Section 4 presents a qualitative comparison of bulk runup statistics from the model with existing experimental data, distribution of runup maxima and averaged swash spectra from ensembles of simulations of irregular wave incidence on three beaches of different slopes. Conclusions are given in Section 5 with particular emphasis on future work on extreme wave runup modelling.  Madsen and Sørensen (1992), are employed to model the propagation and shoaling of weakly nonlinear, weakly dispersive waves from moderately deep water to the breaker line. When waves propagate into shallower water, the effects of nonlinearity become significantly greater than those of frequency dispersion and it is necessary to switch to the more appropriate nonlinear shallow water equations. The propagation of broken waves (modelled as bores) in the surf zone and the wetting and drying at the shoreline are naturally described within a Godunov-type finite volume framework. In particular, evolution of weakly discontinuous solutions of the NLSW equations (e.g. bores propagating up the beach) can be modelled using the shock-capturing capability of such finite volume schemes.
Wave generation is achieved by implementation of a moving boundary piston paddle wavemaker through the introduction of a local mapping of the time-varying physical domain in the vicinity of the paddle referred to as the paddle domain. A second-order wave generation methodology based on the theory derived by Schäffer (1996) provides the motion of the piston boundary. Finally, the switch between the Boussinesq and NLSW equations is determined from the current breaking location, identified from the most offshore wave where the local free-surface slope exceeds a prescribed threshold.
Numerical implementation of the hybrid solver is described in some detail by Orszaghova et al. An obvious limitation of this numerical model compared to other Boussinesq-type models is the weakly nonlinear approximation inherent in the enhanced dispersion equations. More sophisticated numerical models exist (e.g. Shi et al. (2012) and Tissier et al. (2012)) which utilise the 'fully nonlinear' Boussinesq equations (of Chen (2006) and Green and Naghdi (1976), respectively). Therefore, the accuracy of the model predictions of the propagation and shoaling of steep waves in the pre-breaking region will suffer due to the lack of sophistication of the underlying equations. A further limitation of the present numerical model is that all waves inshore of the most offshore breaking point are treated with the NLSW equations. Consequently, the influence of frequency dispersion on smaller waves inshore of a large breaker is neglected by the present breaking algorithm. Tissier et al. (2012) noted that many hybrid Boussinesq-NLSW models do not allow the termination of breaking, thus making them unsuitable for irregular wave transformations over uneven bathymetries. Both Tissier et al. (2012) and Tonelli and Petti (2012) introduce criteria to allow the governing equations to switch from Boussinesq to NLSW and also to revert from NLSW to Boussinesq so that the tracking of individual breaking wave fronts is possible. However, in surf zones of moderate length such as those considered herein the effect of neglecting frequency dispersion for small waves inshore of large broken waves on the shoreline dynamics is likely to be small. Therefore, the model described here provides a good balance between accuracy and computational speed -an important factor when collecting large sets of runup statistics.

Generating-absorbing sponge layer -incorporation in existing model
In the present model, the in-built piston paddle wavemaker allows accurate reproduction of wave flume experiments involving piston-type wavemakers. In fact, if the experimental wavemaker moves subject to displacement control then the same paddle signal could be used as input to the numerical wave flume. Evanescent fields arising in the vicinity of a piston wavemaker will be differ significantly in the physical wave flume compared with the Boussinesq numerical wave flume where the assumption of fairly shallow water leads to small or negligible evanescent mode contributions.
Therefore, the second-order correction to the linear input piston-paddle signals for the numerical contributions, similar to the approach described by Barthel et al. (1983) rather than the theory derived by Schäffer (1996) which takes full account of evanescent modes effects. Nevertheless, the numerical wave flume suffers from the same disadvantages as its physical counterpart, in particular the difficulty in absorbing long waves reflected at the beach and those long free waves released at breaking. In a recent numerical investigation involving the model described in Section 2, Orszaghova et al. (2014) emphasised the importance of generating the correct long wave bound harmonic structure for runup and overtopping estimates. In certain laboratory flumes, effective absorption of reflected, otherwise outgoing waves in the linear incident frequency range has been achieved using active wave absorption (based on force feedback or free-surface elevation feedback control). Higuera et al. (2015) has recently implemented such active absorption in a numerical context. However, we seek to model irregular wave runup in cases where the incident irregular wave field is accurate to second-order so that outgoing reflected and free waves outside the linear frequency range must also be absorbed. Noting that we require wave absorption over the entire incident frequency spectrum, including sum and difference frequencies, we adopt a sponge layer damping zone.
The sponge layer approach was originally proposed by Israeli and Orszag (1981) for a general wave system, and has proved widely effective at absorbing outgoing waves of different celerities. Larsen and Dancy (1983) pioneered its use in shallow water Boussinesq-type modelling for the absorption of backward generated (from a free-surface source function) and reflected waves. A combined generating-absorbing sponge layer approach with second-order wave generation has not been implemented in a numerical wave flume before for the purposes of irregular wave runup modelling, as far as we are aware. By implementing these existing methods together in a Boussinesq-NLSW model, it is possible to conduct simulations of irregular wave runup at enhanced accuracy.
The traditional sponge layer damping condition is recovered by setting imposed free-surface elevation η 0 to the still water level, i.e. ζ 0 = 0 and q 0 = 0. However, if the imposed solution is an incident wave solution of the enhanced dispersion Boussinesq equations then the total wave will be damped towards this incident wave on the sponge layer. Therefore, in the absence of other wave disturbances the sponge layer will generate an incident wave whereas if, for example, a wave disturbance comprising incident and reflected components enters the sponge layer the reflected wave will be absorbed. The imposed solution can be linear or nonlinear, and can encompass regular or irregular waves. However, to avoid the propagation of error waves outside the sponge layer it is most useful to use solutions of the homogeneous (undamped) governing equations. In a numerical context, such solutions can be obtained directly from an existing numerical solution if no analytical expression exists.
The α 1 (x) terms in both mass and momentum conservation equations, dominate the damping control so that α 2 (x) is not strictly necessary, as will be demonstrated by means of numerical examples. The width of the sponge layer L s is typically assumed to be one or two wavelengths.
Although longer sponge layers are preferable to minimise reflections, it has proved necessary to specify a relatively compact sponge layer in the irregular wave simulations (discussed in Section 4).
The spatial profile f (x) of the sponge layer is location-dependent. If the sponge layer is situated Ls 0 f (x) dx = 1 so thatα 1 andα 2 denote the integrated sponge layer strengths.

Irregular wave modelling -coupled simulations
The generating/absorbing sponge layer is utilised to simulate long periods of irregular wave incidence on a beach using the following method, similar to nested domain approach demonstrated by Zhang et al. (2014). Assume the simulation domain consists of a piston wavemaker at the left end of the tank, a flat-bed section, wherein a generating-absorbing sponge layer is placed, followed by a plane beach of uniform bed slope extending from the beach toe through the still 9 3 4 5 6 7 8 9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  water level and beyond as illustrated in the lower diagram in Figure 2. An offshore incident wave field (η 0 (x, t), q 0 (x, t)) is imposed on the sponge layer offshore of the beach at every time step so that simultaneous absorption of reflected outgoing waves and propagation of incident waves is achieved. Any outgoing reflected wave components which interact with incident waves between the shore and the sponge layer are damped out within the sponge layer. The incident wave field offshore of the beach is determined by simulation of incident wave propagation in the absence of beach reflections. The incident wave propagation computational domain must have an identical flat-bed offshore geometry to the irregular wave runup domain which is continued inshore instead of the beach as represented in the upper and lower panels of Figure 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 be challenging to generate accurately the total incident wave disturbance from still water using proportional sponge-layer control. Instead, we seek solely to absorb the outgoing reflected wave component of the total wave signal. At each time step the incident wave solution is determined first, providing input for the sponge layer in the beach runup simulation, and then the beach runup solution is obtained.
The presence of the moving boundary wavemaker and the associated paddle domain complicate the imposition of the sponge layer boundary condition. Mappings to and from the physical paddle domain (modelling the expanding and contracting of the physical locations of the grid points) introduce extra terms in the enhanced dispersion Boussinesq equations (Orszaghova et al., 2012) and so it is more straightforward to restrict the sponge layer boundary conditions to the region beyond the paddle domain. Consequently, in the runup simulation the sponge layer is placed between the end of the paddle domain and the beach toe. The beach toe provides a natural shoreward limit for the sponge layer because the beach is omitted from the incident wave simulation. The absence of the beach in the incident wave domain is necessary to avoid any reflections from the beach occurring, as these would be re-reflected by the wavemaker causing contamination of the incident wave field. Therefore, the wave field obtained from the incident wave simulation beyond the beach toe cannot be utilised in the full simulation.
It is notable that a distinctly non-physical sponge layer is introduced to a model that incorporates a moving boundary piston wavemaker which mimics a laboratory wavemaker. Given that the sponge layer action, which cannot be realised physically, will prevent comparison of numerical simulation results with experimental measurement it seems contradictory to combine the two approaches. A fixed-position wavemaker which fluxes waves through the offshore boundary would suffice and, indeed, have the advantage of allowing an extended generating-absorbing sponge layer. However, given the accuracy of the second-order generation achieved using the moving boundary wavemaker developed by Orszaghova et al. (2012), it was considered expedient to retain the in-built wavemaker in the implementation described herein featuring the sponge layer terms. The generating-absorbing sponge layer is efficient even when quite compact and so any small outgoing long waves will be sufficiently damped in the region between the paddle domain and beach toe. An alternative formulation would see the analytical second order solution derived by Madsen and Sørensen (1993) for the enhanced dispersion Boussinesq equations directly imposed on the sponge layer thus utilising the generation capabilities of the sponge layer. For long irregular simulations, however, calculation of solution at every point on the sponge layer at every time step would prove costly in terms of computational times. Longer sponge layers, desirable for 11 the accurate generation of large waves from still water, would exacerbate this problem.
To verify the foregoing coupled incident wave propagation/irregular wave runup method for modelling irregular wave-incidence on a beach it is necessary to determine damping coefficients that minimise sponge layer reflections and transmissions. Absorption characteristics of the sponge layer are investigated for both focused wave group and regular waves incidence. The spectrum from which the components of the focused wave group are determined and from which the frequencies of the regular waves are selected is identical to the incident wave spectrum used in the irregular wave runup simulations, the results of which are presented in Section 4.

Sponge layer performance
In the simulations to be presented in Section 4, we seek to damp out reflections of broad-banded incident wave-trains while propagating the incident wave-trains over a compact sponge layer. In this context, it is useful to first consider reflected and transmitted waves due to focused wave group propagation over a sponge layer of width equal to two peak wavelengths. A typical value of the integrated damping strengthα 1 which yields good absorption is given by 4ω as observed in regular wave tests for incident waves of frequency ω for sponge layer widths of one or two wavelengths. Similar values forα 2 are also specified -although given the 1/h 2 term it is more difficult to specify a damping strength providing good absorption properties in general.
The computational domain geometry is chosen to be typical of shallow water flumes found in many experimental facilities. Thus, the still water depth at the paddle was chosen to be h SW = 0.5 m. The energy density spectrum of the incident focused wave group corresponds to that implemented by Hunt (2003) et al. (2012). Given that beach reflected waves comprise a fraction of the total incident wave energy to be damped, this moderately-sized incident wave is considered a reasonable test for  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 the damping effects of the sponge layer. A frictionless flat bed is considered, in order to model damping due solely to the sponge layer. Solutions were found to have converged for a uniform grid spacing ∆x of 0.02 m and a time step ∆t of 0.004 s.
Three different pairs of integrated damping strengths (α 1 ,α 2 ) are considered for a sponge layer width of approximately two peak wavelengths. The first two cases, (α 1 ,α 2 ) = (4ω p , 4ω p ) and (8ω p , 8ω p ), involve equal 'Newtonian cooling' and 'viscous' damping strengths. In the third case, only the proportional or 'Newtonian cooling' damping term is non-zero (α 1 ,α 2 ) = (4ω p , 0). The in a different simulation, from the total wave field. In the incident wave simulation, the central sponge layer is removed and replaced by a sponge layer extending four wavelengths inwards from the right hand boundary of the numerical wave flume. Integrated damping strengths of (α 1 , α 2 ) = (4ω p , ω p ) yield negligible reflection and effectively total absorption. Figure 3 shows the reflected and transmitted wave free-surface displacement time histories normalised by the focused wave envelope amplitude at the locations equidistant from the foucs (so that defocusing effects are separated from sponge layer absorption). The effect of the viscous damping term can be seen by comparing the relative amplitudes of the reflected and transmitted waves for the three pairs of damping strengths considered. Although inclusion of the viscous damping term does reduce wave transmission, it seems it is at the expense of increased reflection.
The optimal sponge layer specification from the three sets proposed is for (α 1 ,α 2 ) = (4ω p , 0) where both reflected and transmitted wave amplitudes are less than 0.5% of the focused wave amplitude. Further tuning is possible (α 2 = ω p yields slightly smaller reflections) but the benefits are insubstantial. In the coupled simulations, it is inevitable that any beach-reflected or offshoretravelling free waves transmitted through the sponge layer towards the wavemaker will be damped upon re-entering the sponge layer after reflection at the wavemaker boundary. Minimising sponge layer reflection is thus significantly more important than minimising sponge layer transmission. Therefore, the viscous damping term is set to zero in most cases considered hereon.
Absorption performance of wall-adjacent sponge layers is typically assessed by measuring the reflection coefficient, i.e. reflected wave amplitude normalised by incident wave amplitude, over a range of frequencies. A reflection analysis such as this is most efficiently achieved using a broadbanded irregular incident wave-field. However, re-reflections at the wavemaker inevitably become simultaneously. Therefore, regular waves at several different frequencies in the incident wave spectrum were utilised in the reflection analysis. Estimation of the reflected wave amplitudes can be obtained prior to the time of arrival of re-reflected waves (estimated from the group speed of the regular wave train) using the same approach to that outlined for the focused waves. Figure 4 shows the reflection coefficient versus frequency for two wall-adjacent sponge layers of sponge layer strength (α 1 ,α 2 ) = (4ω p , 0) and widths λ p and 1.5λ p in domains of similar geometry to the focused wave group analysis. (Note the distance through which the waves travel within in the wall-adjacent sponge layer is 2λ p and 3λ p .) Absorption performance improves considerably with sponge layer width. The minimum frequency (0.171 Hz) at which the regular wave reflection analysis was conducted is significantly shorter than the lower cut-off frequency f min = 0.330 Hz in order to assess long wave absorption. The sponge layer is still quite effective at absorbing these longer waves (reflected waves 2% of the incident waves) implying that, despite the compact nature of the sponge layer in the irregular wave simulations, the outgoing long waves will be damped.

Demonstration of an irregular wave runup simulation
This section provides a practical demonstration of the effectiveness of the coupled irregular wave propagation and beach interaction simulations described in Section 3.1. Rather than investigating the amplitudes of reflected and transmitted waves outside the sponge layer (both of which would be difficult to isolate in a domain featuring incident and beach-reflected waves) the behaviour of the shoreline is used to measure the success of the generation-absorption zone.
The energy spectrum for the irregular wave train is again based on one used in the UKCRF tests conducted by Hunt (2003). The spectrum is discretised into 157 components corresponding to a repeat period of approximately 245.76 s. The mean zero-crossing period of the sea state is approximately 1.75 s so that on average each repeat period comprises 140 waves. The significant wave height of the offshore spectrum is chosen to be H s = 0.10 s corresponding to a steepness coefficient of (H s /λ z ) 0.03 where λ z is the wavelength associated with the mean zero-crossing period. Amplitudes and phases of the linear components of the random realisation of the sea state are obtained using the random-amplitude/random-phase method (Tucker et al., 1984), i.e. the amplitudes are determined from the Rayleigh distribution with scale parameter S(f )∆f and the phases from the uniform distribution on (0, 2π), respectively. The second-order wavemaker signal comprising sub-and super-harmonics is computed based on the wavemaker theory of Schäffer (1996).
The computational domain geometry is based on the UKCRF wave basin: a beach with a 1:20 slope is present in the full runup simulation, extending from the beach toe located 8.33 m from the equilibrium position of the paddle to beyond the still water level. In the region between the paddle and the beach toe, the still water depth is 0.5 m. The paddle domain extent is approximately ten times the maximum paddle sweep (equal to approximately 0.25 m from Figure 5). Thus, the 15  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 Figure 5 shows both the first-order (erroneous) and second-order (corrected) wavemaker paddle signals over a single repeat period. It is evident that the paddle must undergo large excursions subject to the second-order wavemaker correction in order to eliminate the spurious long error waves arising from first order theory, and so obtain the correct bound subharmonic structure.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63 17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62 63 64 65 Figure 7 shows the difference in shoreline displacement obtained over the first and second repeat periods as shown in Figure 6 (a). The discrete nature of the predicted shoreline motion is particularly evident in the shoreline difference time-history because the magnitude of the runup difference is not much larger than the shoreline step size ∆x = 0.02 m and so is affected by the wetting and drying algorithm. Beyond the initial start-up transient for the first repeat period, the shoreline positions are generally no further apart than two or three grid cells, thus revealing a satisfactory agreement between the simulation predictions during the first and second repeat period. It should be noted that discrete changes in the shoreline position time-history over time increments of ∆t = 0.004 s cannot be represented adequately in Figure 7 and give the time-history plot a block-like appearance. Based on this comparison of the shoreline displacement differences for the sponge layer strengths α = 4ω p and α = 6ω p , it appears the stronger sponge layer gives better repeatability and is used as standard for the simulations in the following section.

Irregular wave runup simulations -results and analysis
Multiple numerical simulations of irregular wave runup on beaches are now used to assemble a comprehensive set of runup records for a one incident spectral shape (with various target significant wave heights). Simple bulk and extreme-value statistical analyses are conducted on the ensembles of runup records. Qualitative comparison of the results with measured laboratory data should at least partially validate the numerical model. Although the 1DH hybrid model has several theoretical and empirical deficiencies (discussed in Section 2) particularly with regard to nonlinear wave behaviour during shoaling and the mechanics of wave breaking, which might lead to inconsistencies between numerical predictions and physically-measured data, the model is fast and robust allowing large datasets for statistical analysis to be constructed efficiently.
It should be emphasised that the primary interest of the present investigation is not to model conditions and runup in the field as accurately as possible -instead it is to understand the processes that affect wave runup probability distributions (extreme or otherwise) under idealised conditions. Such conditions, where the beach slope is uniform in space and constant in time, are more representative of shallow water flume laboratory than field conditions. However, both existing field data and laboratory data are used as benchmarks for comparison with model results.  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 beyond which there is a beach of uniform slope. Irregular wave runup on the beach is characterised by the offshore wave spectrum (specified by its shape and the significant wave height H s and peak period T p parameters), offshore depth and beach slope. No other length-scales are needed.

Beach slopes
Uniform bed slopes of 1/10, 1/20 and 1/40 are considered with the intention that the slopes represent reflective, intermediate and dissipative beach types, respectively. The bed slope is sometimes defined as tan β where β is the angle of inclination of the beach. In the literature, there is no precise consensus on what beach slopes correspond to what types; steep reflective beaches are sometimes considered to have slopes greater than 1/10. However, generally speaking, dissipative gently sloping beaches have gradients less than 1/20 and reflective slope beaches are of order 1/10 and greater, as stated by Baldock et al. (1998). Care should be taken not to choose a very steep slope as the assumption of a slowly-varying bathymetry underlying the Boussinesq governing equations may be invalidated. The intermediate slope of 1/20 is often specified in laboratory wave flumes (Borthwick et al., 2006;Hunt-Raby et al., 2011). In addition to these three representative slopes, a slope of 1/30 is also considered when analysing the bulk statistics of runup in Section 4.2.1 although only a small ensemble of irregular wave runup records are gathered in this case.

Offshore spectrum and simulation times
The spectral shape, cut-off frequency, and peak frequency are identical to those specified in

Collation of runup statistics
Random realisations of irregular sea states are first obtained by assuming the sea surface to be a Gaussian random variable. A finite wave record is constructed at the paddle by a Fourier summation of wave components using the random-phase random-amplitude method described in Section 3.3. The corresponding linear paddle signal can be obtained using the Biésel transfer function, and corrected by applying the second-order transfer function derived by Schäffer (1996).
For each random realisation of the spectrum, the paddle domain width (defined as ten times the maximum paddle excursion) and hence the sponge layer width vary.
Free-surface elevation and shoreline position time-histories resulting from the propagation of the incident wave up the beach are computed from the second-order wave generation irregular wave runup simulations. For a given shoreline elevation time-history, the elevation maximareferred to as runup maxima or crest -are determined as follows. As shown in Figure 8, the minima of the shoreline elevation are first identified (marked with circles in the figure). In most  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  For every random realisation of the prescribed offshore conditions we can determine the maximum runup from the shoreline record and so construct a distribution of extreme runup elevations. Each extreme runup value corresponds to the maximum shoreline elevation that occurs during the incidence of (on average) 280 random waves over an interval of 491.52 s. Of course, the number of runup crests from which we choose the maximum value is smaller than the number of offshore waves due to surf zone interactions.

Statistical analysis of beach runup
The statistical properties of the ensembles of shoreline elevation records are investigated next, with the bulk statistics and certain extreme value statistics of the collated runup maxima quantified for each beach slope. Coastal engineers and managers are interested in the extreme runup value in a given period, e.g. the largest runup during a storm, and so of most relevance is the nature of the extreme value distribution of runup maxima. Offshore data is often provided in three hour segments because the wave statistics, characterised by a certain significant wave height and peak period, tend not to vary much over such an interval. A three hour wave record from a storm with a typical 10 s zero crossing period will comprise approximately 1000 waves and so the largest runup in 1000 waves is often of interest despite the arbitrary threshold. By considering the ensemble of incident sea-state realisations, and the associated shoreline motions, it is hoped that we can determine (qualitatively at least) how often these extreme runup events 21 occur. However, in order to assess the influence of beach slopes on shoreline motion and the driving physical processes, it is useful first to consider the bulk statistics of the local shoreline maxima or runup crests.

Ratio of number of runup crests to incident-wave crests
In a laboratory study of random wave runup Mase (1989) calculated from the measured shoreline and free-surface elevation time series, among other runup quantities, the ratio of the number of runup crests to incident wave crests for a wide range of beach slopes and incident wave conditions. Mase then proposed an empirical formula for this ratio α in terms of the surf similarity parameter or Iribarren number ξ = tan β/ H s /L 0 , where L 0 is the peak offshore wavelength, having previously shown (Mase and Iwagaki, 1984) that α is governed by this surf similarity parameter alone. A useful preliminary test for the irregular wave runup model described herein is to show that the ratio α of runup crest to offshore incident wave crest occurrences can be described solely in terms of the surf similarity parameter. Thereafter, we compare the α values obtained from numerical simulation to the empirical formula proposed by Mase (1989). Mase (1989) estimated the ratio α for beach slopes of 1/5, 1/10, 1/20, and 1/30 using a similar flume geometry and incident wave conditions (water depth 45 cm, tank length 27 m, Pierson-Moskowitz spectrum, peak frequencies ranging from 0.4 Hz to 1.2 Hz) to those described in Section 4.1.2. Time series of offshore random waves were generated for different peak spectral frequencies, with observed significant wave heights ranging from 5 cm to 10 cm. The length of the runup and free-surface elevation records in each random wave test corresponded to 650-900 individual waves. The empirical formula for the ratio α, expressed in terms of the surf similarity parameter is: height H s = 0.025 m on a 1/10 beach slope are also simulated in order to obtain a surf similarity parameter ξ > 0.91 for comparison with the second case in the empirical formula in Equation (3).
The final case ξ > 3.57 corresponds to highly reflective conditions where every incident wave causes a shoreline motion and is therefore not investigated. The full range of incident wave conditions, beach slopes and ensemble averaged α values are presented in Table 1. The incident free-surface elevation offshore of the beach toe is obtained as part of the coupled runup simulation and so the number of waves N W in the offshore region is straightforward to determine using a zero up-crossing method. In particular, the number of incident waves is counted at one location x = 8.0 m just offshore of the beach toe. Given the discrete nature of the domain at the shore, care must be taken in counting the runup crests or shoreline elevation maxima. However, adopting the definition of a runup crest from Section 4.1.3 it is straightforward to compute the number of runup crests and the associated runup crest to wave crest ratio for each sea-state realisation.
Empirically predicted values for the ratio of runup crest to wave crest occurrences α EMP are also included in Table 1 for completeness. Uncertainties are estimated using a single standard deviation from the ensemble of values.
Comparison of the ratio of the number of runup crests to incident wave crests, α, computed for case V and case V I in Table 1 confirms that the shoreline motions predicted by the model depend on the surf similarity parameter alone (within the statistical margins of uncertainty) .  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 This is despite the significant wave height for the waves incident on the 1/20 beach being twice that of the waves on the 1/30 beach. In both cases, the number of runup crests is slightly less than a half the number of incident waves counted offshore of the beach toe. This confirms that our irregular wave runup model captures the surf zone and swash processes observed by Mase and Iwagaki (1984) in a qualitatively correct manner. Physically, it implies that the width of the surf zone is a significant factor in determining the nature of the shoreline motions. A wider surf zone can be expected for larger waves (where breaking will occur further offshore) and/or gentler beach slopes both of which yield smaller surf similarity parameters. Bore capture, whereby a larger bore overtakes and subsumes a smaller bore, will occur on more occasions in a wider surf zone leading to fewer swash motions at the shore. The converse is true for narrower surf zones. To ensure consistency between the target incident wave spectrum and the randomly realised incident Despite inevitable differences in incident irregular wave field generation (and the absence of reflected wave absorption in the laboratory flume) we compare the ratio of the number of runup crests to that of the incident waves obtained from numerical simulation to the empirical formula derived by Mase (1989) in Table 1. Precise details of the wavemaker configuration (paddle geometry and paddle motions) are not provided by Mase (1989). It is assumed that Mase used linear superposition of sinusoidal motions with amplitudes at each frequency determined from the Pierson-Moskowitz spectrum. Rather than prescribing the desired incident (significant) wave height a priori, Mase (1989) specified different input levels for the electrical signals supplied to the wavemaker and estimated the significant wave height based on wave gauge measurements.
Thus, reflected wave contributions were assumed negligible, which may not have been the case for the steepest beaches considered. In contrast, in the numerical simulations which contribute data to Table 1 the significant wave height H s is prescribed beforehand and the corresponding paddle motions determined using a combination of a linear transfer function and a second-order correction.
Identification of runup crests from numerical and experimental shoreline elevation time- histories  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 might also differ. The numerical results in Table 1 overestimate the ratio of the number of runup crests to offshore incident wave crests compared to the empirical formula based on the experimental data by approximately 20% for all but the gentlest beach slope. Several factors may underlie the disparity in the number of runup maxima observed. Small runup crests identified from simulated runup time histories using the criterion illustrated in Figure 8 may not have been identified or be visible in the laboratory measurements of Mase (1989). Furthermore, absence of a second-order correction to the laboratory paddle motion and the presence of re-reflected waves from the paddle in the laboratory flume mean numerical and experimental incident wave fields differ.

Extreme runup distributions
Each shoreline elevation record corresponds to the propagation of approximately 280 incident waves from the paddle to beach, based on the mean zero up-crossing period of the underlying spectrum. A reasonable estimate for the mean zero up-crossing period at full scale is T z = 8.0 s and thus the time record corresponds to about 37 minutes of data. This duration is a similar order of magnitude to those of the free-surface elevation wave records used to characterise a sea state and similar to the runup record length considered in the extreme value statistical analysis by Holman (1986). Note that this record length, which equals one repeat period of the discretised offshore wave energy density spectrum, has been chosen somewhat arbitrarily. A different discretisation of the energy spectrum would yield a different repeat period. Nevertheless, the records may be sub-divided or combined to form smaller or larger sample sizes from which to obtain the maximum runup elevation. We analyse the statistics for the maximum predicted runup elevation in a single wave record and sets of two combined wave records.
In the following we examine distributions of runup maxima for evidence of an upper limit for runup. For each bed slope, simulations of wave runup have been conducted for 100 random realisations of the incident sea state. The maximum runup elevations from each of the shoreline records of one repeat period in duration are then collated to yield an extreme-value type distribution comprising 100 samples. Note that the number of runup maxima (from which the extreme value is selected) in each shoreline elevation record varies according to the beach slope -for the 1/10 slope, 1/20 slope and the 1/40 slope the number of runup maxima will be approximately 80%, 50% and 30% of the number of incident waves, respectively. Thus, in terms of extreme value statistics, the distributions are not exactly equivalent. Nevertheless, in the context of coastal (and ocean) engineering the return period of a runup maximum for a given set of offshore conditions will be of substantial interest, and so it is more consistent to consider the same time interval of wave incidence for each beach slope. As noted by Holman (1986), the fitting of an extreme value frequency distribution to a relatively small set of data points creates large uncertainties in the details of the upper tail. The return period of a 'rare event' can be derived from the tail of an extreme value distribution; however, it is very sensitive to the tail shape and so we do not seek to obtain such an estimate. Instead, a comparison between the distribution shapes of the maximum runup crest occurring in the shoreline elevation time records is presented for each of the three beach slopes. Figure 10 shows the distribution of the maximum runup crests obtained from shoreline elevation records of one full repeat period in length for each beach slope.
In each case there appears to be a smaller secondary peak in the upper tail of the distributions.  Figure 11. It might have been hoped that evidence of an upper limit to runup would become clearer for ensembles of longer shoreline records; however, the smoothness of the distributions is then limited by the reduced number of samples.

Saturation in swash spectra
Evidence of runup saturation in certain frequency bands has been observed by many authors when analysing swash spectra on gently sloping beaches. Raubenheimer and Guza (1996) 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  investigated swash saturation on beaches with intermediate slopes while Ruessink et al. (1998) and Ruggiero et al. (2004) investigated swash saturation observed in data sets collected on dissipative natural beaches. Saturation was assessed by examining the properties of swash or run-up energy density spectra obtained from the field measurements and by investigating the correlation of infragravity or sea swell swash heights with beach slope and/or offshore significant wave height. Ruessink et al. (1998) averaged the measured runup data in the alongshore direction whereas Ruggiero et al. (2004) considered measurements along each transect of the beach independently giving a range of runup data for different beach slopes. In many of the field studies reported in the literature, it is difficult to isolate the dependence of swash on one particular environmental condition (e.g. beach slope or offshore significant wave height) due to continuous changes in these conditions.
Using the numerical wave flume modified to model irregular wave runup described in this paper, we are able to constructing large runup data sets from which reliable swash statistics can be obtained for any permutation of environmental conditions. Therefore, we can investigate the dependence of swash on any given factor by holding all other conditions constant. To demonstrate this capability, swash spectra on each slope are determined from ensembles of 16 runup records on each of three representative beach slopes (1/10, 1/20, 1/40) for two different offshore significant wave heights (0.05 m and 0.10 m) and compared in Figure  spectra are obtained in two ways: first, using a direct FFT of the shoreline elevation time histories and secondly, using Welch's power spectral density estimate with zero overlap and an averaging factor of 8 (see documentation of the pwelch function in Matlab R ). Thick lines in Figure 12 indicate the averaged power spectral density estimate and are smoother than the direct FFT approach (thin lines). The dashed black line in each log-log plot indicates the estimate for the division between the infragravity and sea swell frequency bands at 'experimental' scale. At full scale, this threshold is 0.05 Hz and by assuming offshore depths of 15 m in the field (Ruessink et al., 1998), the corresponding value for a numerical wave flume of depth 0.5 m can be obtained.
A variety of interesting swash behaviours are evident in Figure 12. Spectral roll-off proportional to f −4 is clearly observed for all beach slopes. However, as noted by Mase (1988) this may be largely due to the parabolic nature of the runup time histories as illustrated in Figure 8. Swash spectra for steeper beach slopes are shifted towards higher frequencies than spectra for gentler slopes, consistent with observations in the field by Ruggiero et al. (2004). Saturation of runup, identified in Figure 12 as where the swash spectra for the two different offshore significant wave heights are equal, is observed to occur for lower frequencies for the gentler beach slopes. In particular, on the dissipative beach the saturated portions of the spectra extend a significant distance into the infragravity band, as observed by Ruessink et al. (1998) and Ruggiero et al. (2004).
The occurrence of saturation in the sea swell band for the reflective beach is notable in Figure 12. Two peaks in the swash spectra are observed for this relatively steep beach slope ,  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 one in the sea swell band and one in the infragravity band. The spectral peak in the sea swell band for offshore wave heights H s = 0.05 m and H s = 0.10 m occurs at approximately at 0.40 Hz and 0.36 Hz, corresponding to swash periods of 2.5 s and 2.8 s respectively. Swash periods can be estimated from Table 1 using T p /α (where T p = 1/f p = 2.16 s) and are approximately 2.5 s and 2.8 s for the relevant cases II and III. Saturation of sea swell swash has occurred on the reflective beach. However, larger swash excursions due to increased infragravity energy remain possible for larger incident significant wave heights. Some evidence of a secondary peak in the sea swell band is present for the intermediate beach. Only for the dissipative beach does all incident wave energy appear to have been down shifted to the infragravity band.

Influence of 'precursor' waves on total runup
It is well known that the properties of preceding waves and the size of the resultant backswash can have a substantial influence on subsequent swash motions. For example, such swash interactions have been observed by Mase and Iwagaki (1984) and Mase (1989) in laboratory random wave runup tests and more recently Erikson et al. (2005) modified a semi-empirical model based on the NLSW equations to account for such interactions. Therefore, a brief investigation is now described into the number of preceding waves required in order to yield accurate shore conditions preceding the arrival of a prescribed incident wave group. Given the importance of extreme runup events, a particularly large wave runup crest from the ensemble of irregular wave runup simulations for the 1/20 beach slope is first considered. In Figure 13, the origin of the time axes is specified to coincide with the occurrence of the extreme runup event and so waves, bores, and shoreline motions can be identified relative to this extreme event; however, in the actual simulation the event occurs after more than 450 s. The analysis is also repeated for large runup events on the beaches with steeper (1/10) and gentler (1/40) slopes.
To investigate how preceding waves affect wave setup and, ultimately, the extreme runup crest, a series of simulations are conducted with progressively fewer waves preceding the arrival of the main group yielding the large runup event, and the corresponding shoreline elevations compared. A particular difficulty in assessing the effect of incident waves on runup arises from the bore-bore interactions in the surf zone. Interactions such as large bore 'capture' of smaller bores preclude establishing a one-to-one correspondence between offshore waves and runup crests (as observed in Section 4.2.1). Thus, the precise effect of each wave on shoreline elevation is not  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 as to how many preceding waves are necessary to represent fully an extreme runup event. In order to control the number of 'precursor' waves arriving at the shore immediately before an extreme runup event, an absorbing sponge layer spanning the entire beach is included in the runup simulation. With reference to Figure 2, an absorbing sponge layer similar to that in the incident wave simulation is imposed from beach toe to shore in the runup simulation. Incident waves steadily decay to zero over the sponge layer and the initial damping strength is specified so that no shoreline motions occur. To allow a certain number of waves to arrive at the shore prior to the extreme runup event we simply ramp the sponge layer strength to zero at the appropriate time. Figure 13 presents time-histories for free-surface displacement at the mid-point of the beach, at the still water shoreline, and shoreline elevations (with respect to SWL) for one simulation involving a full repeat period of incident waves and three further simulations where the sponge layer strength is set to zero several periods before the occurrence of the extreme runup event.
Convergence of the shoreline motions arising from short periods of irregular wave incidence on the shore to the fully established shoreline motions depends on the number of waves reaching the shore prior to the arrival of the set of waves causing the large runup event. The number of undamped waves arriving at the shore in turn depends on time of flight of the waves from the beach toe to the (still water) shoreline. The first bore to arrive at the shore after the sponge layer is 'removed' may incorporate some of the damped, reduced amplitude waves that are present in the absorbing sponge layer prior to its removal. Consequently, to select the time at which to 'remove' the absorbing sponge layer so that the first waves to reach the shore are those forming the main bore that leads to the extreme runup event requires some subjective judgement (based on the time of flight estimate). (Note that the absorbing sponge layer is 'removed' by ramping the strength to zero over half a mean zero-upcrossing period T z /2.) By ramping the sponge layer strength to zero 9.0 s before the maximum runup event, the first significant waves to arrive at the shore form the main bore, yielding the extreme runup event, as shown by the red-dashed line in Figure 13. This simulation is referred to as runup case A. Compact focused wave groups typically comprise three significant crests flanked by two smaller crests and so the duration of significant wave motions in a typical compact wave group in such a sea state may be roughly approximated as being equal to three mean zero-upcrossing periods (3T z = 5.25 s). Thus, by reducing the sponge layer strength to zero at times 14.25 s (runup case B) and 19.5 s (runup case C) before the extreme runup event, we can allow approximately three and six waves or one and two wave groups, respectively, to arrive at the shore just before the main bore. In this manner, the effect 33  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64 3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  also. However, the rundown or backswash occurring after the runup event at t = −5.0 s in case B appears to impede the propagation of the main bore up the beach because the predicted maximum runup elevation is smaller than for the fully-established shoreline motions. Nevertheless, the free-surface elevation at the still water shoreline is still over-predicted relative to the free-surface elevation of the fully-established irregular wave train at the shore. Runup case C, allowing approximately six fully-formed waves to propagate past the mid-point of the beach leading to three bore motions at the still water shoreline, provides the best agreement with the shoreline dynamics -both in terms of the free-surface elevation of the main bore at the still water shoreline and the shoreline elevation.
Figures 14 and 15 compare moving shoreline elevations and still water shoreline free-surface displacements around an extreme runup event for different numbers of precursor waves on the 1/10 and 1/40 beaches, respectively. Shoreline motions for case A comprising the extreme bore motion only with no precursor waves approximate the actual extreme shoreline excursion reasonably well with relative differences compared with fully established shoreline excursions of just less than 10%. Improved agreement is achieved in cases B and C where one and two wave groups are allowed to precede the wave group causing the extreme bore motion. This is particularly evident on 1/10 beach slope where both the extreme bore motion and the preceding bore motion are 35  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 almost exactly reproduced by Case C where two precursor wave groups are allowed to arrive at the shore. In fact, on the reflective 1/10 beach slope only one precursor wave group is necessary to precondition the shore so that that the following bore motion is accurately reproduced. On the 1/40 beach, it appears that one precursor wave group (case B) is also sufficient although the reproduction of the bore motions at the still water shoreline is not as accurate as for the 1/10 beach.
In all cases, the maximum shoreline elevation is reasonably well approximated by incident wave trains which, after passage through the surf zone, yield just the main bore. Better agreement is achieved by including a precursor wave group (approximately three waves) with the main incident wave train which preconditions the shoreline to resemble the fully established shoreline conditions. Even on the gentlest slopes where the surf zone is longest and swash zone interactions may be significant, a single precursor wave group appears to lead to an excellent approximation of the extreme shoreline motion.

Conclusions
A coupled incident-runup model based on a hybrid 1DH Boussinesq-NLSW solver incorporating a generating-absorbing sponge layer was used to create long duration simulations of undirectional irregular wave runup on uniform, gently sloping beaches. The generating-absorbing sponge layer facilitated simultaneous propagation of incident irregular waves and absorption of outgoing reflected waves. The offshore incident wave field was first simulated in a domain without the beach present. Incident and runup solutions were then coupled by imposing the incident wave solution on the sponge layer in the runup simulation at each time-step. An in-built moving piston paddle wavemaker allowed accurate generation of incident waves to second order; it is essential to include correctly the subharmonic component for accurate runup and overtopping modelling (as previously demonstrated by Orszaghova et al. (2014)).
Before application to long duration irregular wave runup tests, the damping strength of the sponge layer was tuned to absorb focused wave groups (having already obtained preliminary estimates for suitable strengths with regular waves). The damping strength was then further tuned to absorb irregular waves reflected from a beach, using a sponge layer of width equal to one peak wavelength width in a coupled irregular wave propagation and runup simulation. Convincing evidence of the need to absorb reflected waves was demonstrated by comparing shoreline motions for two coupled simulations involving the tuned sponge layer strength and negligible sponge layer strength (approximating zero absorption) .  1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60  61  62  63  64  65 Statistics of runup maxima for irregular wave incidence (with correct second-order bound wave structure) on beaches of three different slopes were compiled from a large ensemble of long irregular wave runup simulations. In each case the input offshore wave field was taken as a random unidirectional sea state determined from a prescribed wave energy density spectrum. The beach slopes corresponded to reflective, intermediate, and dissipative beach types. In qualitative terms, the runup statistics for the various beach slopes indicated that the 1DH hybrid model was capable of reproducing efficiently the important surf zone interactions for runup on plane beaches, thus allowing compilation of large sets of runup data. Comparison between the predicted ratio of number of wave runup crests to number of offshore wave crests and corresponding results from an empirical formula derived by Mase (1989) for random wave runup in a laboratory wave flume revealed certain discrepancies due in part to differences in the incident wave fields. However, this ratio was shown to be governed by the Iribarren number alone in accordance with Mase's findings.