Extreme Inundation Statistics on a Composite Beach

The runup of initial Gaussian narrow-banded and wide-banded wave fields and its statistical characteristics are investigated using direct numerical simulations, based on the nonlinear shallow water equations. The bathymetry consists of the section of a constant depth, which is matched with the beach of constant slope. To address different levels of nonlinearity, the time series with five different significant wave heights are considered. The selected wave parameters allow also seeing the effects of wave breaking on wave statistics. The total physical time of each simulated time-series is 1000 hours (~360000 wave periods). The statistics of calculated wave runup heights are discussed with respect to the wave nonlinearity, wave breaking and the bandwidth of the incoming wave field. The conditional Weibull distribution is suggested as a model for the description of extreme runup heights and assessment of extreme inundations.


Introduction
Estimating extreme runup events in coastal zones is an important task. Flood prediction received a lot of attention in recent decades, in order to reduce hazard risks in coastal zones [1][2][3][4]. The statistical distribution of wave runup characteristics is influenced by many factors, such as topography and coastline, nonlinearity and wave breaking [5][6][7].
Also, some individual waves at the coast may be unexpected, extreme and hazardous. This regards sneaker waves or freak wave runups [8][9][10]. Such extreme events at the coast often lead to human injuries and fatalities, when people are washed off to the sea from a gentle beach or from coastal rocks or sea walls, and damage of coastal structures. During the period of 2011-2018, there were cases when freak wave runups (unrelated to tsunami) washed cars and motorcycles into the sea and damaged houses and buildings in the coastal zone [10]. These events correspond to the very tails of the statistical runup height distribution and their analysis requires extremely large datasets.
Previous studies have employed different methods to study the statistics of long wave runup, including numerical models, experiments, and field measurements.
Theoretically, [11] studied the statistical characteristics of long waves on a beach of constant slope using an analytical solution of the nonlinear shallow water theory. The study revealed that the runup height was distributed according to the Rayleigh distribution, if the incident wave elevation was described as having a normal distribution and a narrow-band spectrum. In terms of the statistical moments of the moving shoreline on a beach of constant slope, this study asserts that the kurtosis is positive for weak amplitude waves and negative for strongly nonlinear waves, whereas the skewness is always positive. Later [12] showed that for the description of even non-breaking waves the Gaussian distribution is inappropriate. Both theoretical studies had a number of assumptions, which were putting in question the applicability of these results.
Experimentally, [13] tried to reproduce theoretical results of [11] in the wave flume at Warwick University. However, they could not generate "pure" Gaussian wave field. Moreover, the generated waves were affected by capillary effects. Thus, the only result [13] could reproduce regarded an increase in the mean sea-level elevation with an increase in wave nonlinearity attributed to the known phenomenon of wave set-up. They also found that the values of the statistical moments of wave runup (skewness and kurtosis) were similar to those of the incident wave field. [14] studied statistics of narrowband and wide-band wave runups in the Large Wave Flume of the University of Hannover, Germany. They found that wave fields with a narrow-band spectrum were associated with a higher loss of the wave energy compare to the waves with a wide-band spectrum. However, their experimental records were not long enough to discuss freak runups.
In the field measurements, [15] studied runup heights, measured on a wide spectrum of sandy beaches in New South Wales; they found that runup was distributed according to the Rayleigh distribution. [16,17] studied wave runup at Canadian and Australian coasts and demonstrated that wave runup deviates from the Gaussian distribution. Although some of these conclusions were similar to those of [11][12][13], it was not possible to put direct correspondence between these works due to a number of reasons. First, the field measurement studies lacked information about an incident wave field. Second, they had a different bathymetry and coastal topography, deviating from the ideal plane beach. Third, the data included an error associated with measurement techniques.
However, the main issue, which complicates comparison of theoretical [11,12] and experimental [13,14] results, is the insufficient length of the experimental time-series, which do not support analysis of extreme runup statistics. Potentially, this issue can be overcome nowadays with the use of IP highresolution cameras permanently installed on a beach and associated techniques [18][19][20][21][22][23][24]; however, we have not seen such works yet.
In this paper we cover the existing gap in long-term experimental records by using digital data obtained with intensive numerical computations. This approach has clear advantages. It gives control on the initial wave field offshore and allows checking the applicability of the approximated analytical results by [11,12] to a more realistic bathymetry: plane beach merged with the flat bottom.
The paper is organized as follows. In section 2, the numerical model, based on nonlinear shallow water equations is described. The statistical moments and the distribution functions of random wave and runup fields, as well as distribution functions of wave and runup heights, are described in detail in Section 3. Then, the main results are summarized in Section 4.

Numerical Model
In this section, the 1D nonlinear shallow water model, which represents the mass and momentum conservation, is briefly described: Here D = h + η is the total water depth, η (x, t) is the water elevation, with respect to the still water level, x is the coordinate directed onshore, and t is time, h(x) is the unperturbed water depth, u(x, t) is the depthaveraged water flow velocity, and g is the gravitational acceleration. The dimensionless formulation can be obtained by choosing a typical water depth h0 as the length scale (in this problem, the depth of the constant section can be taken as h0√ ℎ 0 , 0 gh as the velocity scale and ℎ 0 /√ ℎ 0 00 / h gh as the time scale. The dimensionless equations take the form of equations (1), (2) with h0 = 1 and g = 1. All computations reported in this study were performed in the dimensionless formulation.
The modelling is performed in the framework of equations (1), (2), which are solved using a modern shock-capturing finite volume method. Although the shallow water model does not pursue the wave breaking and undular bore formation in a general sense (including the water surface overturning), it allows shock-wave formation and propagation with the speed given by Rankine-Hugoniot jump conditions, which, to some extent, approximates wave breaking. The numerical scheme is second order accurate, thanks to the spatial reconstruction (UNO2). For details, see [25].
In this simulation, the corresponding bathymetry ( Figure 1) set-up is used: the flat part of the flume matches the beach of constant slope: where h0 is the constant water depth, kept at 3. The spatial grid step is, therefore, 25 cm, which corresponds to 4 cm vertical resolution for runup height. This was done in order to limit simulation time, when running 10 000 hrs of physical time of wave propagation. However, this also implies that we have a low resolution and not so reliable statistics especially for small amplitude waves Hs = 0.1 m. In a similar manner to the significant wave height, Hs, the significant runup height, Rs, is introduced as an average of one third of the largest runup heights in the time-series. The significant runup height for this small amplitude case is Rs = 0.23 m, so even in this case the resolution is low, but considerable. Of course, the number of extreme runups in this resolution is also somehow underrepresented, however all qualitative and comparative conclusions of this study still hold on.

Boundary Condition
On the left extremity x = a of the computational domain, the Dirichlet boundary condition on the total water depth component D(a, t) = h0 + η0(t) of the solution (D, Du) is imposed. Namely, the free surface elevation function, η0, is drawn from a narrow-or wide-band Gaussian signal depending on the experiment. This data turns out to be enough to obtain a well-posed initial boundary-value problem provided that the flow is subcritical at the point ,, which is always the case for Riemann waves (see [26] for the rigorous mathematical justification of this fact in case of transparent boundary conditions). The boundary conditions are implemented in the finite volume scheme according to the method described in [27], see also [28] for more details on the application to the nonlinear shallow water equations). The considered boundary condition (wave field offshore) is distributed by the Gaussian distribution: where, σ is a standard deviation, and μ is a mean value of the distribution. To ensure this, all individual time-series have been verified by the Kolmogorov-Smirnov test [29].
The spectrum of the generated waves is ( ) where f is the wave frequency, Δf is the frequency band, f0 = 0.1 Hz is the central frequency, and S0 is the constant, which is calculated in order to achieve the desired Hs.
In this work, the case with Δf/f0 = 0.1 is referred to narrow-band spectrum, while the case with Δf/f0 = 0.4 is attributed to the wide-band spectrum. In order to study the influence of wave nonlinearity during wave propagation to the coast, waves of different significant wave heights, which is calculated as averaged of one third of the largest wave heights in the time-series (Hs = 0.1 m, 0.2 m, 0.3 m, 0.4 m, and 0.5 m) are considered. The calculated time-series for each Hs is 1000 hours (360 000 wave periods). Parallel computations have facilitated the calculation of the statistics of wave runup characteristics for 5000 hours, for each bandwidth, and 10000 hours in total. The numerical computations have been carried out in MATLAB and run on a cluster containing 28 cores.
Parameter of the nonlinearity for generated waves is estimated as Hs/h0 and is changing from 0.03 to 0.14. The characteristic parameter kh0 = 0.38 is at the border of validity of the shallow water theory taking into account the horizontal extent of the wave tank. The phase velocity relative error committed by nondispersive theory for kh0 = 0.38 is only 2.3%. Thus, at the end of the numerical wave tank the difference between wave crest positions (between dispersive and non-dispersive models) is less than 10%. Since the focus of this paper is on wave runup, the choice of the theory is justifiable. The choice of wave parameters allows us to see the effects of wave breaking on statistics of their runups. The type of wave breaking is defined by the Iribarren number [5]: where H is the wave height, and L is the characteristic wavelength offshore. It is surging or collapsing for

Data analysis and results
Figure 2 shows probability density functions (PDF) of narrow-band and wide-band wave fields for different nonlinearities, Hs/h0. The data of the narrow-band spectra, Δf/f0 = 0.1 are shown by triangles (different colors correspond to different nonlinearities), while the corresponding Gaussian distribution ( = 0, = 0.25) is shown by the black solid line. The data of the wide-band spectra, Δf/f0 = 0.4 are shown by pluses, and the corresponding Gaussian distribution ( = 0, = 0.27) is shown by the red solid line. It can be seen that the generated waves are well described by the Gaussian distribution, which has zero mean, skewness and kurtosis for all nonlinearities, Hs/h0.
To describe the wave statistics in Figure 2, the Rayleigh distribution, which is well used for this type of problem [5], is applied: where ξ is a data vector, is the scale parameter. For a better fit, a two-parameter Weibull distribution is also considered: where is the scale parameter, and k is the shape parameter.
The wave height distributions of both narrow-band and wide-band wave fields are shown in Figure  5. As expected, the narrow-band data are well described by a Rayleigh distribution ( = 0.5), although a Weibull distribution gives a slightly better fit ( = 0.74, k = 2.27). The data of wide-band spectra tend to be distributed according to a Weibull distribution ( = 0.71, k = 2.06).
The waves which are twice higher than the significant wave height (H/Hs ≥ 2) are the so-called freak waves. It can be seen from Figure 3 that the probability of the freak wave occurrence in the initial wave field is higher for narrow-band signals than for wide-band ones.
The calculated significant runup heights Rs for narrow-band and wide-band signals are shown in Figure 4. Interesting to see that Rs for wide-banded waves is always higher than for narrow-banded waves, which can be explained by higher variability in wave periods for wide-banded waves. Also, figure  4 indirectly shows us how many of our waves are breaking. The wave runup height, at which the first wave breaking occurs in the wave trough can be estimated as Rcr /h0 = g(αT/(2π)) 2 /h0 = 0.2, see for details [30]. This means that our case of "small" nonlinearity Hs/h0 = 0.03 is very little affected by wave breaking (< 1% according to Iribarren criterion). The case of Hs/h0 = 0.06 is affected by wave breaking only for extreme runups (32-35% according to Iribarren criterion). In the case of Hs/h0 = 0.09, more than a half of waves are breaking (61-65% % according to Iribarren criterion). However, in the cases of Hs/h0 = 0.11 and Hs/h0 = 0.14 the majority of waves are breaking. Figure 5 shows the probability distribution functions of runup oscillations, r/Rs for initial Gaussian narrow-banded and wide-banded wave signals. It can be seen from Figure 5a, that runup oscillations of narrow-banded waves are no longer distributed by a normal distribution, and are slightly shifted to the right towards larger positive values with an increase in nonlinearity. Partially this effect was observed both theoretically for an infinite plane beach [11,12] and experimentally [13,14]. What is interesting and peculiar is a strong deformation of the distribution itself. In addition, the tails of these distributions are much thinner than of Gauss, and reflect a relatively weak probability of extreme floods for narrowbanded waves.
The distributions of runup oscillations of initial wide-band signal are also shifted to the right towards higher runups with an increase in nonlinearity, but this shift is much larger compared to the one of the narrow-band signal. Moreover, the tails of these distributions are much thicker than those for narrowband data, and are rather close to the normal distribution, which corresponds to a relatively large probability of extreme floods for wide-banded waves.
It can also be seen that for both narrow-banded and wide-banded waves, the probability of large waves decreases with an increase in wave nonlinearity, which can be explained by wave breaking.    These effects can also be seen in Figure 6, which shows the statistical moments of narrow-banded and wide-banded waves offshore, normalized by Hs, and the corresponding runup oscillations on a beach, normalized by Rs. The statistical moments, mean, variance, skewness, and (normalized) kurtosis are calculated as: where ξ is a data vector, and n is its length. Noteworthy, the mean, skewness and kurtosis of both narrow-banded and wide-banded wave fields are zero, which provide the desired Gaussian statistics. Regarding runup oscillations, one can see that for both narrow-and wide-banded waves the mean of runup oscillations rises with the nonlinearity, which reflects the known effect of wave set-up on a beach. For small-amplitude waves, the set-up for narrowbanded waves is larger than for wide-banded ones, while for large amplitude waves, affected by wave breaking, it is the opposite. For wide-banded waves, the variance decreases with an increase in nonlinearity, while for narrow-banded waves it changes non-monotonically. The higher moments, skewness and kurtosis of runup oscillations for waves with narrow-band spectrum are negative, while for waves with wide-band spectrum they are sign-variable. Also for the narrow-banded waves, the skewness decreases with an increase in wave nonlinearity, while kurtosis changes non-monotonically with an increase in wave nonlinearity. Moreover, for wide-banded waves, both skewness and kurtosis change non-monotonically with an increase in nonlinearity. This somehow only partially corresponds to the theoretical findings in [11], where the kurtosis was positive for weak amplitude waves and negative for strongly nonlinear waves, while the skewness was always positive. However, in the experimental study of [13], the skewness was both positive and negative. It is also important to say that for all four moments one can see different dynamics for small-amplitude non-breaking or almost non-breaking waves and large-amplitude waves, strongly affected by wave breaking. Runup oscillations deviate from the Gaussian distribution even for weak-amplitude waves (see Figure 6). With an increase in nonlinearity, all statistical moments of runup oscillations change. It can be seen that statistical moments of narrow-banded irregular waves (except kurtosis) change with Hs monotonically, while for the wide-banded waves, they vary non-monotonically (except mean values).
The large (extreme) wave runup heights, Rextrm = R/Rs ≥ s, where s is some threshold value, somehow behave similar to a conditional Weibull law whose density is given by Eq. (11): A conditional Weibull law is characterized by three parameters: the shape k, the scale λ and the threshold s. Given the data (Ri extrm) = 1...n, s is fixed and k and λ are computed by maximum likelihood estimator. The scale parameter, λ can be obtained from Eq. (12): where n is a number of extreme wave runups. In order to obtain the shape parameter, k, one should solve Eq. (13): Similarly to freak waves, the waves on a beach, whose runup height is twice larger than the significant runup height (R/Rs ≥ 2), we call freak runups. On gentle beaches, such freak runups are manifested as sudden floods and may result in human injuries and fatalities [8][9][10]. Figure 7 shows probability distribution functions of large runup heights (R ≥ 0.7 Rs), for narrowband and wide-band spectra, for different nonlinearities. It can be seen in Figure 7, the tails of distributions for runup heights corresponding to freak events for narrow-banded waves decay much faster than those for wave heights offshore (except waves of weak amplitude with Hs/h0 = 0.03), which means that for narrow-banded waves the probability of freak runup occurrence on a beach is less than the probability of freak wave occurrence in the sea coastal zone and a gentle beach works as some kind of "filter" for narrow-banded freak events. This is also manifested in the numbers of actual freak events, given in Table 1. It can be seen that for non-breaking waves of the smallest amplitude Hs/h0 = 0.03, the number of freak events on a beach was reduced twice compared to the original number of freak waves offshore, while for waves of larger amplitude, which were affected by the wave breaking, there were no freak runups at all.
In contrast, for wide-banded waves the probability of freak events on a beach is more or less the same as in the sea coastal zone and may even be higher (Figure 7). The number of freak runups for small non-breaking wide-banded waves increased twice as compared to the original number of freak waves offshore (see Table 1). With an increase in wave amplitude (and consequently, wave breaking), the number of freak runups on a beach decreases, however for waves of moderate amplitude, the number of freak runups is still larger than the number of freak waves offshore, while for waves strongly affected by the wave breaking (Hs/h0 = 0.11 and 0.14), the number of freak runups on a beach suddenly drops down (see Table 1). Probability of extreme wave runups on a beach is noticeably higher for waves with wide-band spectrum, than for waves with narrow-band spectrum (see Figure 7), although the probability of extremes wave heights in the wave field offshore is significantly higher for narrow-banded waves than for widebanded ones (see Figure 3, Table 1).
The probability of extreme runup formation changes with the wave nonlinearity. It is decreasing with an increase in wave nonlinearity for wide-banded waves and changes non-monotonically with nonlinearity for narrow-banded waves. It is also interesting to see that the tails of distributions in Figure 9 are somehow gathered into clusters and can be separated in two groups for "relatively large Hs" and "relatively small Hs", where the "small Hs" group is always higher than the "large Hs" group. The latter holds for both narrow-banded and wide-banded waves and can be explained by the wave breaking.
The corresponding data of wave runup heights are also approximated by a conditional Weibull distribution [Eq. (11)], which gives reasonable results and can be used to evaluate the probability of freak runups. Here the threshold s is selected as 0.7 and the calculated parameters k and λ are given in Table 2.

Conclusions
In this paper, irregular waves runup on a plane beach is studied by means of direct numerical simulations. The numerical model is based on the nonlinear shallow water equations and is of the second order of accuracy. The corresponding bathymetry consists of the section of constant depth, which is matched with the beach of a constant slope. The irregular waves are represented by the Gaussian wave field with spectra of two different bandwidths, which are referred to as narrow-banded and wide-banded waves. To address different levels of wave nonlinearity, time-series with five different significant wave heights are considered. The selected wave regimes represent (i) non-breaking waves, (ii) waves slightly affected by wave breaking, (iii) moderate wave breaking and (iv) significant wave breaking, when the majority of waves are breaking. Each of these time-series has a duration of 1000 hours (360 000 wave periods).
The heights of narrow-banded waves are well described by Rayleigh distribution, while heights of wide-banded waves are described by Weibull distribution irrespective of the wave nonlinearity. However, for wide-banded waves the very tails of these distributions show larger variability than for narrow-banded ones.
As expected, the runup oscillations are not Gaussian, which confirms results of many previous studies, both theoretical [11,12] and experimental [13,16,17]. For both narrow-band and wide-band cases, one can observe the effect of wave set-up (increase in the mean value of runup oscillations), which increases with an increase in wave nonlinearity. However, for wide-banded waves this increase is significantly stronger than for narrow-banded ones.
What regards extreme, so-called "freak events", their statistics in the initial narrow-banded wave signal offshore is more representative, than on the beach ("freak runups") even for non-breaking waves. Therefore, for narrow-banded waves, gentle beaches reduce the number of freak events as compared to the sea coastal zone, and work as a 'low-pass filter' for extreme wave heights. This may explain why freak events on a beach are so unexpected [8][9][10]. However, for wide-banded waves, such an effect has not been observed and the probability of freak events on a beach was similar to or even larger than the one in the sea coastal zone.
The number of freak events in wide-band and narrow-band cases varies, so that increase in the bandwidth leads to a substantial increase in the number of freak events. This can be explained by higher variability in wave periods for wide-banded waves, and wave runup height is rather sensitive to these variations. In addition, the number of freak waves decreases with an increase in wave amplitude and consequently, wave breaking. The largest number of freak waves was observed for non-breaking widebanded waves, which almost doubled the number of freak waves in the boundary condition wave record.
Finally, to describe statistics of extreme wave runup heights on a gentle beach, a conditional Weibull distribution is suggested. It gives reasonable results and may be used for assessment of extreme inundations on a beach (freak runups). In addition, in future applications the statistical analysis hereby provided might also be useful in the study of the wave run-up phenomenon in other applications, e.g. in structures placed in shallow water conditions [31,32].
The limitation imposed by the resolution of the numerical simulations should also be taken into account. Although the number of freak waves on the beach may be somehow reduced by a coarse model resolution, the qualitative and comparative conclusion of this study should not be affected. This point will be improved in our future studies.