Estimation of incident and reflected components in nonlinear regular waves over sloping foreshores

The present paper deals with separation of long-crested regular waves into incident and reflected components. Such methods have been available for several decades for linear waves, but has recently been extended to cover nonlinear waves over horizontal foreshores. The overall goal of the present paper is to extend the separation method for nonlinear regular waves to cover also sloping foreshores. This requires the combination of the existing method with a nonlinear shoaling model. A nonlinear shoaling model was very recently found valid for the shoaling of the primary and bound components in regular waves when the slope angle is positive and mild. In the present paper this shoaling model is utilized and assumed valid also for the de-shoaling of the reflected waves, i.e. on a negative mild slope angle. However, if the reflected waves are nonlinear the de-shoaling process is much more complicated and will for example cause release of free waves. Interactions among those free reflected wave components may cause nonlinear interactions not included in the mathematical model. For that reason, the applicability range is limited to mildly nonlinear reflected waves. Using numerical model data with various foreshore slopes, wave nonlinearities and reflection coefficients the reliability of the developed model is examined in detail.


Introduction
To assess the performance of coastal structures in physical model tests, it is fundamental to know the incident and reflected waves present in the facility.This is because the target waves to be reproduced in the physical model is usually specified in terms of the incident wave parameters.Moreover, the response of coastal structures is dependent on the incident waves hitting the structure and not the total waves.However, only total waves can be measured, consisting of both the incident and reflected waves.Thus it is common practise to use either of the following methods to estimate the incident waves the structure is exposed to: • Method 1: Measure the waves without the model in place and assume that the total waves during these tests equal the incident waves with the model in place.This can be done by dividing the flume into two with a longitudinal guidewall or by repeating the tests without model in place.• Method 2: Measure in a single coordinate (x,y) several properties and calculate incident waves based on a mathematical model of these properties and a mathematical solver.Usually, the used co-located measurements quantities are velocity (u,w) and elevation (η) or pressure (p) which are used in various combinations.• Method 3: Similar to Method 2, but instead of co-located measurements, arrays of surface elevations (η) gauges are used.Thus the surface elevation is sampled in several (x,y) coordinates (η 1 , η 2 , …, η N ) and a mathematical model of the surface is used to separate into incident and reflected wave trains.
Method 1 is used by many hydraulic laboratories.The main advantage is that the incident waves can be measured exactly at the wanted position (usually at the toe of the structure) without any errors introduced by wrong mathematical models.The major disadvantage is that repeating the tests without the model doubles the testing time.This is avoided by using the guidewall, but this setup has less accuracy due to diffraction effects from the wall tip.Another disadvantage is that errors arise due to the reflection from the passive absorber during the tests without the model in place, which will typically be of 5-15 % of the incident wave height.Moreover, during the tests with the model in place, some degree of re-reflection at the wavemaker will always occur even for the best active absorption systems.This re-reflection is typically 0-10 %, according to Lykke Anderen et al. (2016).If the surface elevation gauge is placed at a sufficient distance to the reflective boundaries, then such small reflection will usually lead to very small errors on the incident spectral significant wave height (H m0 ) in irregular waves.This is because the nodes and antinodes do occur at different locations for the various frequencies in the spectrum.However, for regular waves as studied in the present paper it might be much more important as only one primary frequency is involved.The accuracy of the estimated incident waves thus depends on the reflection and re-reflection coefficients.
The main disadvantage of Methods 2 and 3 is that if the mathematical model does not reflect reality, errors are introduced.In most cases, the mathematical model used is based on linear wave theory.The consequence of using a mathematical model consisting of linear waves only is that errors occur if waves are nonlinear.
Method 2 with co-located measurements has the advantage over Method 3 that the spatial variations due to, for example, a sloping seabed or wave breaking is not influencing the performance of the methods significantly.However, similar restrictions on wave type are imposed for Method 2 where linear long wave theory was used by Guza and Bowen (1976), Guza et al. (1984), Tatavarti et al. (1988), Kubota et al. (1990) and Huntley et al. (1999).General linear theory was used by Kubota et al. (1990), Walton (1992); Hughes (1993) and McKee et al. (2018).Only the methods of Kubota et al. (1990) uses quasi-nonlinear wave theory, but only terms up to second order are considered and with long wave restriction.In all cases, the bathymetry slope is assumed so small that theory for a horizontal floor is valid.Method 2 is not considered further in the present study as the focus is on nonlinear and highly nonlinear waves on mild and steep foreshores.
From the 1960 ′ ties to the 1990 ′ ties the used mathematical models for Method 3 have assumed linear waves on horizontal bathymetry, cf.Kajima (1969), Thornton and Calhoun (1972), Goda and Suzuki (1976), Gaillard et al. (1980), Mansard and Funke (1980), Zelt and Skjelbreia (1992) and Frigaard and Brorsen (1995).These methods are to large extent still used today, but the last two decades also various nonlinear mathematical models have also been proposed for Method 3. The nonlinear methods are based on two paths of developments where the first is the methods based on local approximations using simulated annealing (LASA) proposed by Medina (2001), Figueres et al. (2003), and Figueres and Medina (2004).The LASA methods are local time domain solutions applying various wave theories for regular waves in local windows.Due to the time domain solution in windows the method is computationally heavy and lead not even to very accurate results, cf.Lykke Andersen et al. (2017) and Eldrup and Lykke Andersen (2019).The other path is based on frequency domain solutions with mathematical models that includes bound and free harmonics for the incident and reflected waves.This was proposed by Lin and Huang (2004) for regular waves, but they ignored amplitude dispersion and is thus only accurate to 2nd order waves.Lykke Andersen et al. (2017) included amplitude dispersion in their model which resulted in accurate separation of regular waves up to depth limitation.This method was extended to bichromatic waves by Lykke Andersen et al. (2019) using 3rd order wave theory.Finally, it was extended to irregular waves by Eldrup and Lykke Andersen (2019a) using a narrow band assumption for the primary spectrum.Based on second order theory the error for wider spectra was studied theoretically.This was supplemented by numerical and physical model data with highly nonlinear waves.Padilla and Alsina (2020) followed a similar path and they furthermore included sloping foreshores by linear shoaling of the free waves and a highly simplified shoaling model for the bound components.
In order to reflect the conditions in the prototype many physical model tests in coastal engineering involve a sloping bathymetry.If the mathematical model does not include a sloping seabed, then errors are introduced on the estimated incident and reflected waves.Baldock and Simmonds (1999) extended the linear method by Frigaard and Brorsen (1995) with a linear shoaling model (Burnside (1915)).They showed that for linear waves, the error on the estimated incident and reflected amplitudes might be up to 90 % when the sloping bathymetry is ignored.A similar study was performed by Chang and Hsu (2003) but using a frequency domain solution.In the present paper focus is on nonlinear waves and thus the linear shoaling used by Baldock and Simmonds (1999) and Chang and Hsu (2003) or the simple shoaling model used by Padilla and Alsina (2020) will not accurately describe the wave shoaling.This is because as waves travel into very shallow water, they become highly nonlinear and bound harmonics shoal significantly more than the components they are bound to.As highly nonlinear shallow water waves shoal, the bound harmonics thus increases much more than predicted with the linear shoaling coefficient.Energy from the primary component is transferred to the bound components.Thus the primary component might, in contradiction to the linear shoaling theory, actually decrease in amplitude.This was studied in detail by Eldrup and Lykke Andersen (2020) for regular waves.Eldrup and Lykke Andersen (2020) showed that the shoaling of a nonlinear wave and its components on gentle slopes can be described by stream function wave theory when the energy flux is conserved, cf.Rienecker and Fenton (1981).This was validated for both the total wave height and the amplitudes of the primary component and the bound superharmonics.Taking, for example, a regular wave of height 0.15 m and period 3 s in 0.7 m of water depth that shoals into 0.5 m depth.Then the linear shoaling coefficient is 1.06 while Rienecker and Fenton (1981) nonlinear shoaling coefficient is 1.10 for the wave height.The shoaling coefficients of the individual components vary tremendously as for example the primary component has a shoaling coefficient of 0.96, for the 2nd order component it is 1.67, for the 5th order component it is 7.9 and for the 10th order component it is 82.1.Such behaviour is not included in the simplified shoaling model used by Padilla and Alsina (2020).
The validity of the shoaling model studied by Eldrup and Lykke Andersen (2020) was given in terms of the water depth to wavelength ratio (h/L).For the slope angle of 1:100 the found validity range was h/L > 0.036, while for 1:50 it was h/L > 0.073 and finally for 1:30 it was h/L > 0.085.These numbers were found based on the validity of the component shoaling, but without studying the validity of the celerity which is very important for the reflection separation.
To avoid influence of evanescent waves caused by reflective structures Klopman and van der Meer (1999) recommends placing the wave gauge array at a distance of 40 % of the peak wavelength (0.4L p ) from the intersection between the waterline and the structure.In many cases the incident waves at the toe is needed and thus it is important to have a reliable separation method that is also able to extrapolate the results a distance that is up to 40 % of the wavelength.Thus, the lack of a reliable shoaling model has previously forced researchers to use Method 1 instead of Method 3 for steep sloping foreshores.This was, for example, done by Eldrup et al. (2018) and Eldrup and Lykke Andersen et al. (2019b) for their tests with steep foreshores.A reliable shoaling model to be used in the reflection analysis might make it possible to use Method 3, even for such cases.In the present paper, the conditions under which Method 3 is valid is extended and examined in detail.
Usually, energy loses are not included in the mathematical model and thus an example of a limitation is if significant energy loses are present.Energy losses due to bottom friction are usually small in wave flume experiments over the relevant distance for reflection separation.If experiments involve for example soft mud or vegetation this might not be the case and thus the energy loses should be included as for example proposed by Nouri et al. (2014).More relevant will usually be energy losses due to wave breaking.This is though difficult to include, and thus the present method is limited to non-breaking waves on the foreshore.Thus an upper wave steepness limit of the applicability might be estimated by, for example, the Miche (1944) breaker criteria.Unless the seabed slope is mild, the shoaling model though introduces significant errors before this theoretical limit is reached, cf.Eldrup and Lykke Andersen (2020).
Currents will change the celerity and shoaling of waves and thus if currents are present they should also be included in the mathematical model.Suh et al. (2001) presented a method for separation of linear waves on a known current by modifying the dispersion equation.This effect can also be included in the nonlinear method of Lykke Andersen et al. (2017) as the celerity of the bound components is calculated by stream function theory where a uniform and steady external current may be specified.For the free components, the method of Suh et al. (2001) applies.However, in the present paper only the case without currents will be considered.
The mathematical model is usually not including cross-modes.Thus, if cross-modes are present the incident and reflected wave predictions will have errors as studied in detail by Grønbech et al. (1996).Cross-modes should be avoided/minimized in 2-D physical model tests and will thus not be included in the mathematical model and is not considered further in the present paper.
The goal of the present paper is to include the nonlinear shoaling model by Eldrup and Lykke Andersen (2020) into the nonlinear separation method for regular waves by Lykke Andersen et al. (2017).Based on numerical model data the errors on the estimated incident waves is studied for the new method and compared with the predictions by the existing method based on horizontal foreshore.

Theory
In the present section first the method proposed by Lykke Andersen et al. ( 2017) for separation of nonlinear regular waves on horizontal foreshore is presented.Then the extension of the method to sloping foreshore is presented.

Original formulation (horizontal bathymetry)
The mathematical model used by Lykke Andersen et al. (2017) includes primary incident (η I (1) ) and reflected (η R (1) ) components as well as incident free (η I,F ) and bound (η I,B ) harmonics and reflected free (η R,F ) and bound (η R,B ) harmonics and finally also noise (e m ).The surface elevation at the x-coordinate of gauge m (x m ) is given by: η(x m , t) = a (1) where N is the highest order analysed and should be chosen high enough to include all the relevant wave energy.The wavenumbers of the free waves (k (n) ) are estimated from the linear dispersion equation (at the cyclic frequency nω).If the wavenumbers of the primary waves (k I and k R ) are also known, then all the amplitudes and phases in Eq. ( 1) can be obtained by a least-square approach at the primary frequency as well as at the higher harmonics (2 ≤ n ≤ N).Hereafter, the incident (H I ) and reflected (H R ) wave heights can be obtained from the amplitude and phases of the primary and bound components.However, as the wavenumbers of the primary components (k I and k R ) for nonlinear waves depends on the wave heights (H I and H R ) an iterative procedure or nonlinear solution method is needed.Lykke Andersen et al. (2017) proposed to use Rienecker and Fenton (1981) to calculate the wavenumbers (k I and k R ).An iterative procedure is used where the starting guess for the wavenumbers is based on H I = H R = 0 and based hereon the first estimate of the incident and reflected surface elevations including bound harmonics is obtained.This leads to new estimates of H I and H R and thus new k I and k R wavenumber may be calculated by Rienecker and Fenton (1981).This iterative procedure continues until convergence for k I and k R .Liu and Li (2016) worked simultaneously with Lykke Andersen et al. ( 2017) on the same problem.They suggested that k I and k R to be unknown parameters to be found by the mathematical solver.Thus the stream function wave celerity calculation is avoided, but on the other hand, it leads to a nonlinear mathematical solver being needed.The main purpose of the present paper is to study sloping foreshores and thus Liu and Li (2016) procedure for wavenumbers cannot be used and thus wavenumbers will be found by stream function wave theory as in Lykke Andersen et al. ( 2017).Furthermore, one should be aware that sometimes the bound and free waves at a given frequency cannot be reliable separated.That would be the case when the celerity of the two are almost identical.Lykke Andersen et al. ( 2017) thus included corrections for such cases and proposed for the problematic frequencies only to include the total amplitude instead of both the bound and free in the mathematical model.

New formulation (sloping bathymetry)
The aim of the present work is to modify the mathematical model given in Eq. ( 1) so that it includes shoaling.Shoaling means that the celerity and amplitude of the components changes with the x-coordinate.Because the model in Eq. ( 1) includes nonlinear regular waves, also the shoaling model should be valid for such waves.Thus, amplitude corrections due to shoaling are different for every component.Eldrup and Lykke Andersen (2020) showed that for gentle slopes, the nonlinear shoaling model of Rienecker and Fenton (1981) is accurately describing the shoaling of the incident primary and bound components.For steeper slopes, errors occur when the wave becomes highly nonlinear (small h/L and large H/L).This is because that Rienecker and Fenton (1981) method is based on horizontal seabed and thus the surface profiles of the waves are vertical symmetric around the crest.Observed waves on sloping bathymetries are vertical asymmetric and the asymmetry is increasing with the seabed slope.On steep slopes also reflection from the slope might occur and thus the assumption that the reflected wave is only changing by shoaling might not be fully correct.Thus, additional errors might occur for steep slopes and long wave gauge arrays due to reflection from the seabed within the array.Eldrup and Lykke Andersen (2020) provided limits of applicability of the shoaling model given in terms of minimum h/L for the seabed slopes 1:100, 1:50 and 1:30.However, these limits were based on reliable deep to shallow water shoaling coefficients.For reflection analysis it is only needed that the shoaling model is valid within the wave gauge array, but on the other hand accurate estimates of the celerity is at least as important as the shoaling coefficients.Thus, the limits of applicability of their shoaling model in reflection analysis will be studied in the present paper.
The Rienecker and Fenton (1981) shoaling methodology is to calculate the wave height that gives identical energy flux of stream function waves in the various depths (x-coordinates).Based on this wave height and related wavelength variation with x-coordinate the incident and reflected bound wavenumbers variation with x-coordinate are also known.From the wave height variation also the amplitude of all the bound waves components are found by Rienecker and Fenton (1981) and thus the variation of the shoaling coefficient of each component with x-coordinate is calculated by the ratio of the amplitude of each component at the given depth and the reference depth.
A similar model for the de-shoaling of the reflected primary and bound components is used.The celerity and amplitude change for the reflected waves is though different due to a different wave height.This was chosen even though Eldrup and Lykke Andersen (2020) showed that on slopes of 1:30 and 1:100 release of free harmonics from the reflected bound waves occur.They found that only for a very gentle slope of 1:1200 this did not occur.The released free harmonics interact and causes new nonlinear components at the same frequencies as the original wave (primary and bound components) but with different celerities.By assuming the de-shoaling model is identical to the shoaling model errors will thus occur unless the reflected waves are only mildly nonlinear or the seabed slope is very gentle (1:1200).The shoaling model for the free higher harmonics is based on linear wave theory and thus identical to Baldock and Simmonds (1999).
A reference coordinate (x r ) is defined and the wave heights of the primary wave including the bound higher harmonics components are here . The phase shifts due to x-coordinate shift from the reference coordinate x r to gauge m located at xcoordinate x m are defined as follows: where k I (x) is the wavenumber that depends on x-coordinate due to the changes in the depth as discussed above.This includes also variation in the celerity with x-coordinate caused by depth changes and wave height changes due to shoaling (amplitude dispersion).Thus, the variation of k I (x) is calculated by requiring conservation of the mechanical energy flux following Rienecker and Fenton (1981) procedure.Above integrals need thus to be evaluated by numerical integration.In Rienecker and Fenton (1981) either the current speed (U) or the discharge (q) may be provided.In Lykke Andersen et al. (2017) the current speed was assumed to be zero and thus a Stokes drift is present (U = 0, q > 0).However, after equilibrium has been established in a wave flume, the discharge would be zero and thus a return current is established to balance the stokes drift (q = 0, U < 0).Thus, both options were examined in the present paper, but the consequences on the obtained results were negligible for all tested cases and thus only results with U = 0 is presented.
The free higher harmonics are assumed of such a small amplitude that linear theory is valid.Otherwise, interaction terms between them and the other components would also be present.This assumption was already present in Eq. ( 1) for the method for horizontal floor.
The shoaling coefficients are evaluated according to: (5) where A (n) ) is the amplitude of the nth order incident and reflected components at the reference coordinate (x r ) and calculated from stream function waves with the depth at x r and the wave heights ) is calculated in the same manner, but using the depth and wave heights at x m .
K F is the shoaling coefficient of the free waves between the depth at x r and x m and is calculated by linear shoaling at the cyclic frequency nω (Eqs.( 7)-( 9)).K I,B and K R,B are evaluated for the primary frequency (n = 1) and bound superharmonics n > 1.If at the reference depth the energy is low at a given superharmonic frequency, then the shoaling coefficients at especially high order might be very large as the wave travel into shallower water and become nonlinear.The presence of noise on that frequency might thus be amplified significantly and especially if the waves are extrapolated to a shallower output location outside of the measurement array.Thus, an upper limit for the component shoaling coefficient needs to be included.K < 25 is used in the present paper, which seems like a reasonable compromise between the risk of amplifying noise and the distance that the waves can be extrapolated from the wave gauge array.For normal wave gauge arrays with a length up to 0.4 times the wavelength this limit has even on steep foreshores no or little effect.However, it might restrict the distance that the waves can be extrapolated from the array, especially for steep foreshores.For steep foreshores, the Eldrup and Lykke Andersen (2020) shoaling model anyway starts to become unreliable when the waves become highly nonlinear and thus above limit seems not to put significant additional constrains on the application area of the present method.

K
(1) ) are the shoaling coefficients between x r and x m of the incident primary and bound components calculated with Rienecker and Fenton (1981) as described above.
) is similar but for the reflected wave components.
The applied mathematical model using above definitions is: The amplitudes (a) and phases (ϕ) in above model are given at the reference coordinate x r .
An iterative least-square procedure similar to the one used for Eq. ( 1) may be used for Eq. ( 10).Note also that just as for Eq. ( 1) there might in shallow water cases be problems in separating the bound and the free harmonics.That would be the case if the celerity of the two components is almost identical.In the case of sloping foreshores this has to take into account that the celerity of both the bound and the free varies through the array, and thus it is evaluated according to: where θ is evaluated with x r at the first gauge and x m at the last gauge.
Mathematically α = 1.0 corresponds to a singularity in the mathematical T. Lykke Andersen and M.R. Eldrup problem at frequency nω.However, it is not sufficient that α > 1.0 as otherwise small errors on the bound or free celerity or gauge positions will influence results significantly.A safe value for α seems to be 1.15, while for α = 1.05 to 1.15, reliable results are sometimes obtained, but the sensibility of results must be checked.Additionally, the following criteria is used to assess if the separation might be unreliable: where x r and x m again should be taken at the first and last gauge.This criterion is used as it is not sufficient that the celerity has a certain difference if the array is too short to lead to sufficiently large phase differences.A safe value seems to be β = 0.08π while for β = 0.03π to 0.08π sensibility of results has to be checked.Usually, H I > H R and thus the problem occurs only for the reflected waves (Eqs.( 12) and ( 14)) and in some cases for both reflected and incident waves (also Eqs. ( 11) and ( 13)).Also, the problem occurs only for shallow water waves and usually only for the first bound harmonics (usually only for n = 2).The corrections performed are the same as used by Lykke Andersen et al. ( 2017).
The new method has been implemented in the wave analysis software WaveLab3 (2020).

Numerical model data and analysis methodology
As in Eldrup and Lykke Andersen (2020), the test cases are generated by the MIKE 3 Wave Model FM by DHI.This model is a phase-resolving, fully nonlinear numerical model that is based on the Reynolds Averaged Navier-Stokes equations with a k-ε turbulence model.The free surface is described by a height function.
Regular waves were generated by stream function wave theory inside a relaxation zone also used to absorb the reflected waves reaching the deep end of the model.The waves were generated on a constant depth (h) of 3 m using a wave height (H) of 0.1 m.The large depth leads to almost linear waves being generated by the wavemaker.Furthermore, the long relaxation zone ensures re-reflections are insignificant.The overall conclusion is thus that no free incident higher harmonics is being generated in the model.Three different wave periods were used (T = 2.4, 3.7 and 5.0 s).Foreshore slopes of 1:100, 1:50 and 1:30 were used to shoal the waves until being close to breaking.These nine test combinations were performed with a sponge layer located after the slope and on a horizontal plateau.Thus only minor reflections from the foreshore slope should be present in the model.The purpose of these nine tests is to study the effect of including the nonlinear shoaling model or not on the estimated wave parameters.The mathematical model, as well as the numerical model, is only valid for non-breaking waves.Thus the test cases were planned so that the water depth at the end of the foreshore slope (h s ) for every test was set to 115 % of the depth where the shoaled incident wave is expected to break (h b ) according to the method by Goda (2010).For this breaking calculation, the shoaled wave height was calculated by Kweon and Goda (1996), but the waves are shoaling more in the numerical model than calculated by Kweon and Goda (1996).Thus, the waves reach for the three test cases on the 1:100 slope between 100 % and 106 % of the breaking wave height, while for the 1:50 slope they reach between 99 % and 102 % and for the 1:30 slope between 89 % and 92 %.For both the slope 1:100 and 1:50 slopes, there is one test where the waves exceed the Goda (2010) breaking criteria in the last wave gauge array and for the above reasons, these two arrays are not used in the present paper.
Additional setups where the sponge is replaced by an impermeable structure with slopes 1:3, 1:6 and 1:9 are used to demonstrate the performance of the separation method when reflected waves are present.To avoid breaking caused by the interaction of incident and reflected waves, the depth at the toe of the impermeable structure was selected to 0.6 m and thus larger than the minimum depth in the models with the sponge.The setups used are shown in Fig. 1.For the setup with the impermeable structure six selected test cases out of 27 combinations were selected.The tests selected were all three slope angles for the lowest period (T = 2.4 s) on the steep 1:30 foreshore slope.For the other two periods only the 1:9 impermeable structure was tested and again mainly for the 1:30 foreshore, but the longest periods was repeated also for the 1:100 foreshore.
The horizontal grid cell size used was L 1 (h s )/200 where L 1 is the wavelength calculated by linear theory using the primary frequency and h s is the smallest water depth in the model domain for the cases with the sponge absorber.40 vertical layers were used in all models.The time step was selected based on a CFL number of 0.4.The cases with the impermeable structures applies the same discretization as the identical test with the sponge absorber.
Outputs were given by surface elevation gauges placed with a distance of 0.02L 1 (h s ) and sampled at 50 Hz.These gauges were used to define 100 wave gauge arrays for reflection separation.The method requires four gauges, but five is used in the present study to have an overdetermined system.Gauges for each array were selected to have gauge distances of approximately of x 12 = 0.05L 1m , x 13 = 0.12 L 1m , x 14 = 0.27 L 1m and x 15 = 0.45 L 1m , where L 1m is the linear wavelength at the middle of the array, i.e. at (x 1 +x 5 )/2.In most cases, the incident and reflected wave trains were calculated at x = x 1 , but few tests are used to demonstrate the extrapolation to for example the toe of a structure.This is to demonstrate the calculation of incident waves at the toe if the array is placed with the distance recommended by Klopman and van der Meer (1999) to a reflective structure.The analysed time frame is always after a stationary wave pattern is developed from the incident and reflected waves, and thus the stationarity requirement assumed in Eqs. ( 1) and ( 10) is fulfilled.
The data are analysed with the present method and with Lykke Andersen et al. ( 2017) method.The Lykke Andersen et al. ( 2017) method is for horizontal foreshores only, and thus the water depth used in that method was selected as the depth in the middle of each wave gauge array.

Results for cases with sponge absorber
In Figs.3-5, the estimated wave heights at the reference coordinate x r = x 1 are presented for the nine test cases with the sponge absorber.Each of the defined 100 wave gauge arrays were analysed and thus the amplitudes and wave heights of the incident and reflected wave trains may be plotted as a function of the depth.The wave heights shown are the height of the primary wave including its bound components (referred to as H i,b and H r,b ) the wave height of the sum of the free higher harmonics (referred to as H i,f and H r,f ) and the height of the total waves (H i and H r ) for incident and reflected waves respectively.
In these figures also, the validity ranges of the used shoaling model are included.Eldrup and Lykke Andersen (2020) provided the validity ranges of the shoaling model, but this was based only on the shoaling of the wave height and the individual components from deep water.For reflection analysis, it is the shoaling within a shorter distance that is important and thus larger errors on the component shoaling can be accepted.However, the celerity is very important, and examination of Eldrup and Lykke Andersen (2020) results show that celerity error is acceptable (1-2%) if the limits given in Table 1 are fulfilled.As the depth given in the x-axis of Figs.3-5 is at the first gauge in each array, the validity limits in Table 1 are corrected for the depth difference between x 1 and x 5 in the validity limits plotted in these figures.In order to apply the provided validity limits in other cases, a dimensionless plot of the validity ranges is provided in Fig. 2. In that figure also the breaking limitation according to Miche (1944) is shown.This demonstrates that the validity is quite a bit from the breaking limitation and thus it is interesting to observe how large errors are introduced in the reflection separation from end of shoaling model validity and until breaking limitation.
As only minor reflections from the slope should be present in the performed tests, the total measured wave height at the reference coordinate might be used as the target value for the incident wave height.This target is also plotted as a function of depth in Figs.3-5.Oscillations with x-coordinate occur though in the total wave height and indicate thus that reflections are present.These oscillations match in wavelength with a partial node and antinode pattern due to incident waves and reflected waves from a partial reflective structure.Their amplitude also increases with the foreshore slope angle as expected when caused by reflections.The true incident wave height in the model will thus be the total measured but excluding these oscillations.
The results demonstrate that the incident total wave height is correctly estimated with the new shoaling model up to the limit of validity of the shoaling model (the minimum depth indicated by the black dashed line).It also indicates that outside the validity range the errors on the incident wave height are minor (less than 1-2%) unless the seabed is steep and the waves are highly nonlinear where the error increases to up to 4 %.In some cases, the last few arrays have so nonlinear waves that stream function theory did not converge through the entire array and thus the wave height had to be reduced for convergence and indicate waves are very close to breaking limitation (on a mild slope).This part is shown with dotted line in the figures and demonstrates the error on the incident wave height and the free wave amplitude in many cases increase for that reason.This is especially visible for the wave period of 3.7 s and foreshore slope 1:100, cf.Fig. 4.
Compared to the Lykke Andersen et al. ( 2017) method, improvements are seen in the incident total wave height, but even by using the method for horizontal seabed, small errors are observed unless the waves are highly nonlinear.For highly nonlinear waves the error on the incident wave height with existing method might be up to 12 % and even for mild slopes large errors might occur.It could also be observed that when waves become highly nonlinear Lykke Andersen et al. ( 2017) always overestimates the incident wave height.This is because that method does not consider shoaling and thus present an average over the gauges in the array.Thus the incident wave height in the first gauge (as shown in Figs.3-5) is overpredicted.If instead the last gauge had been shown it would have underpredicted.These errors would be even larger if the waves had to be extrapolated outside of the array.
The differences on the wave height of the incident free higher harmonics are much larger and especially for the steepest foreshore slope (1:30).No free incident higher harmonics are present in the model, and this is also well predicted by the new method within the range of validity of the shoaling model.This is not the case for the original formulation with horizontal seabed where incident free higher harmonics are predicted, and they increase in height with seabed slope and wave nonlinearity.Thus, the separation into bound and free components on sloping foreshores is only reliable when the nonlinear shoaling model is included.
The actual reflected wave heights are unknown, but they should be small and without oscillations with depth.The new method can be seen to provide predictions having much smaller oscillations than the original method for horizontal seabed.This is especially visible for the steepest foreshore and within the validity range of the shoaling model.Outside the validity range spurious reflections are predicted even with the new method, but they are much smaller than predicted by the method assuming horizontal foreshore.To further demonstrate these differences, arrays at two depths are selected for showing a comparison between the estimated incident time series and the measured total.Note as discussed earlier that the total includes a small degree of reflection and especially for the steepest foreshore.Thus, small differences should be acceptable, but large differences indicate errors in the reflection separation.The two depths selected are corresponding to the lowest depth where the shoaling model is valid, and the array at the end of the foreshore slope, respectively.The slopes 1:100 and 1:30 are selected and results are shown in Figs.6-8.For the array placed at the depth where the shoaling model is valid (sub figures a and c), the predicted incident total wave profile is accurately matching the measured total profile no matter which of the two methods are used.However, significant differences exist in the estimated free higher harmonics.Using the present method, the free higher harmonics are correctly estimated to be insignificant in all three cases and minor differences occur only for the most nonlinear wave and especially if the seabed is also steep, cf.Fig. 8c.For the method based on horizontal seabed assumption (Lykke Andersen et al. ( 2017)) significant free higher harmonics are predicted already at this depth and the errors increase with increasing seabed slope and wave nonlinearity.Correctly assessing the height of the free waves thus requires the nonlinear shoaling to be included in the mathematical model.For the arrays at the end of the foreshore slopes (sub figures b and d) significant errors are present for both methods.The new method gives acceptable results when the foreshore slope and wave nonlinearity is mild.For highly nonlinear waves on steep foreshores very significant errors are though present.In T. Lykke Andersen and M.R. Eldrup all cases, the new method leads though to significantly smaller errors than the method for horizontal seabed.Further improvements in predictions require a shoaling model valid in shallower water.
It is quite surprising that even for quite steep foreshores, the Lykke Andersen et al. ( 2017) method has quite small errors on the total incident wave height.This is though caused by free waves being predicted that compensates for the shoaling within the length of the array.Due to the recommendations by Klopman and van der Meer (1999) it is important that it is possible to extrapolate the measurements to the toe of the structure.In Figs.9-10 such extrapolation is examined, and these figures demonstrate that it is essential to use the present method and the nonlinear shoaling model when extrapolating waves outside of the array.Thus, the errors in the mathematical model of Lykke Andersen et al. (2017) given in Eq. ( 1) becomes very evident when extrapolating outside of the array.For the 1:100 foreshore, the present method is very accurately predicting the waves one peak wavelength ahead of the array.This applies both when extrapolating within the validity range of the shoaling model (Fig. 9a) and when the array is placed outside of the validity range and extrapolating until waves are almost breaking (Fig. 9b).In the latter case, there is though a small error in the predicted celerity that cause a slight time shift.In this situation, the method for horizontal seabed underestimates the wave height by up to 20 % and causes a very significant time shift.For the steep 1:30 slope, the present method is only accurately extrapolating the waves inside the validity ranges of the shoaling model.The predictions outside of the validity range of the shoaling model are though significantly better than with the original method, but not accurate enough to avoid repeating tests without the model in place.
For all test cases, significantly improved estimates are found by using the new method including nonlinear shoaling compared to the existing method without shoaling.Inside the validity range of the shoaling model the estimates by the new method are highly reliable for all cases and this applies also when extrapolating wave trains outside of the measurement array.Outside the range of validity of the shoaling model errors are  present, but for mild slopes, the errors are acceptable over the distance needed to fulfil Klopman and van der Meer (1999) requirements.For steep slopes, errors are more significant and for the most nonlinear wave, extrapolating wave trains outside of the measurement array lead to significant errors.This is because steep slopes lead to vertical asymmetric waves and thus the energy flux estimation and shoaling of the components calculated by stream function theory for horizontal seabed introduce errors.Even in such cases, significantly improved estimates are provided compared to the original reflection separation method without shoaling.However, outside of the shoaling model validity, it is not possible to extrapolate waves outside of the measurement array for such cases.Thus it might be needed to calibrate waves without the model in place if highly nonlinear shallow water waves on steep foreshores should be extrapolated outside of the measurement array.Accurate reflection analysis for such cases requires improvements in the nonlinear shoaling model and thus requiring a wave model for nonlinear waves on sloping foreshores.Some efforts on this was done by Swart and Crowley (1989) and Chen et al. (2012).Swart and Crowley (1989) developed a covoidal wave theory for regular waves with a gentle slope of the seabed, where a covocoidal mathematical expression describes the surface elevation.Swart and Crowley (1989) found that the solution did not accurately describe measured wave parameters.Chen et al. ( 2012) used a third-order asymptotic solution to derive a wave theory for regular waves on a sloping seabed based on a Lagrangian description.The method is valid for linear and mildly nonlinear waves only, and thus it is not applicable for the present study.Overall, the results validate the new method and further improvements are only possible by improved understanding of shoaling of nonlinear waves on steep foreshores.

Results for cases with reflection from impermeable slope
Above results considered tests with a highly efficient passive absorber and thus only the reflection from the seabed is present.For such cases, reflection separation is not needed, but may anyway be used to evaluate the performance of the separation method.However, more realistic tests cases are also performed with impermeable slopes of 1:3, 1:6 and 1:9.Mainly the foreshore slope 1:30 was tested, but a single test with the most nonlinear wave is also performed for the 1:100 foreshore.The slope angle is only varied for the shortest period of 2.4 s for which reflection coefficients from approximately 10 %-90 % are found depending on the slope angle.As the numerical model does not well model the energy loses due to wave breaking, these numbers are higher than what would be observed in a physical model with the same structure.Thus the three slopes actually cover the typical range of reflectivity observed for coastal structures including vertical breakwaters.For the other wave periods only the slope angle 1:9 is tested, and this leads to reflection coefficients (H r /H i ) of approximately 30 % for the middle wave period (T = 3.7 s) and 60 % for highest wave period (T = 5 s).This range includes the typically obtained reflection coefficients in hydraulic model tests with rubble mound breakwaters exposed to long waves.
As only part of the primary wave is reflected the reflected wave cannot bind as much superharmonic energy as the incident.If the reflection coefficients is low, the reflected harmonics will thus mainly be free while for the higher reflection coefficients both free and bound reflected energy is present.If the reflection coefficient is low, further release of free harmonics from the bound components during the deshoaling will not be significant and thus the shoaling model must be expected to be valid.However, for medium and high reflection coefficients this is not the case as free energy is released during deshoaling.Moreover, if the free harmonics are not of small amplitude they will interact with the other free harmonics and also with the primary component and lead to sub and superharmonics.Such interactions do not fit in the mathematical model as only two components are assumed at each frequency and such interaction may not fit with the celerity of either of those.Thus, errors in the separation must be expected when the reflected waves cause significant nonlinear interactions.
The tests have been carried out in such a way that the previous tests with the sponge absorber may be used as reference cases.This is because identical incident waves are being generated, and identical measurement positions are used.Thus, the correct incident waves are known if no interaction of incident and reflected waves occur.
For every test, all arrays along the foreshore slope have been analysed, but significant differences in the accuracy of the predicted incident wave trains do not exist between arrays at the toe and those placed for example half a wavelength further offshore.Thus, only the array with the last gauge placed at or closest to the structure toe is considered in the following.This is also the location with the most nonlinear waves and according to above findings it also leads to the highest inaccuracy of the shoaling model.
Fig. 11 shows the results for the tests with the wave period T = 2.4 s.These waves are only mildly nonlinear and there are no significant differences between using the method without or with the shoaling model.However, the present method, including the shoaling model, shows a better match between the target profile and the primary wave including its bound components.For that test, the predicted incident and reflected wave amplitudes are also correctly predicted to be very small with the present method, while they are a little larger for the method without the shoaling model.When reflections are present, a small timeshift in the incident wave is predicted.This might be caused by a small change in celerity over a long distance caused by interaction of incident and reflected waves or small changes in water levels and return currents compared to the tests with the sponge.This timeshift increases with the amount of reflection and also the predicted incident wave profile changes slightly.Thus, the target profile is unknown when reflections are present and thus the results may only be qualitatively assessed based on comparison to the target with the sponge.Instead, the amount of incident free higher harmonics may be used as a quantitative estimate of the accuracy of the separation.This is because that these should be zero in all cases.When the 1:9 slope is considered (approximately 10 % reflection), the free incident higher harmonics are still small and the present method shows a clear improvement over the method without the shoaling model.Moreover, the reflected waves are, as expected, predicted to mainly consist of free higher harmonics and are not significantly different if Lykke Andersen et al. (2017) or present method is used.When the reflection increases to approximately 40 % Fig. 11.Incident and reflected waves at first gauge (x 1 ) for the array closest to the structure toe.The test with a wave period T = 2.4 s and foreshore slope 1:30 is considered.Both a highly efficient absorber and 1:3, 1:6 and 1:9 impermeable slopes with the toe at a depth of 60 cm are shown.
(1:6 slope) the primary reflected wave component is significant and bound harmonics start to occur leading to a higher and narrower crest than the trough, which again seems realistic.For this test, both significant bound and free reflected harmonics are present.Incident free waves are still predicted very small and especially with the present method.Also the total incident wave profile is a reasonable match with the target profile without the structure in place except for a time shift.Thus, for this test, the results seem reliable, but there are only minor improvements by including the shoaling model.For the 1:3 slope the reflection coefficient is around 90 %, and thus significant bound reflected harmonics occur that will be released to free waves as the waves de-shoal.Thus, this test will include high free harmonics that will interact and cause bound sub-and superharmonics by many interactions (not only by self-interaction).This was also shown by Eldrup and Lykke Andersen (2020) that described the de-shoaling process of nonlinear waves as very complicated unless the seabed slope is very mild.Thus, the reflected waves are not behaving as assumed in the mathematical model, which is only valid if free incident and reflected superharmonic wave components are so small that they can be described as linear waves.This is expected to be the main reason for the free incident harmonics being predicted a bit larger for this test.Moreover, significant differences between the predicted incident primary wave including the bound components and the profile observed without the structure in place is predicted.Despite this the total incident wave profile seems not to include significant errors neither with nor without the shoaling model included.However, the separation of the bound and free reflected harmonics seems not correct as the reflected wave profile including the bound components seems unrealistic.The overall conclusion from the tests with the lowest period is thus that the separation is reliable if free harmonics are so small that linearity assumption for these is reasonable.Thus, this puts an upper limit to the magnitude of the reflection coefficient.For this wave period, only small improvements are found by Fig. 12. Incident and reflected waves at first gauge (x 1 ) for the array closest to the structure toe.The test with a wave period T = 3.7 s and foreshore slope 1:30 is considered.Both a highly efficient absorber and a 1:9 impermeable slopes with the toe at a depth of 60 cm is shown.
including the shoaling model.This may also be attributed to the rather small errors found in the tests without the structure in place for an array placed at a depth of 0.6-0.7 m, cf.Fig. 3. Fig. 12 present the result of the wave period T = 3.7 s with a foreshore slope of 1:30 and a structural slope of 1:9.Without the model in place, very little free incident higher harmonic energy is predicted with the present method, but significant spurious free higher harmonic energy is predicted with the method for horizontal seabed.For the method with horizontal seabed, large deviations for the primary incident wave including all the bound harmonics is observed, but the total is almost the same as with the present method.Free waves are thus erroneously predicted when not including the shoaling in the mathematical model.This is because the free waves will partly cancel some of the errors from not including the shoaling model and thus reduce the errors on the total waves.As demonstrated in Figs. 9 and 10 this is though only the case inside the measurement array and extrapolating outside the measurement array will introduce large errors when the shoaling model is not included.With the structure in place, approximately 30 % reflection is found and both methods lead to a small amount of free incident higher harmonics being predicted and in general almost identical results.However, the reflected primary wave including bound components show a slightly unrealistic form indicating errors are present in both methods.This might again be caused by free reflected higher harmonics being significant and causing sub-and superharmonic interactions.
For the highest wave period T = 5 s both the foreshore slope 1:30 and 1:100 is tested with the structure slope 1:9, cf.Figs. 13 and 14 respectively.For the 1:30 slope, the present method is a huge improvement over the method for horizontal bottom by Lykke Andersen et al. (2017).This can be seen by the height of the free incident higher harmonics and also the form of the predicted incident and reflected primary wave including the bound components.However, both methods lead to Fig. 13.Incident and reflected waves at first gauge (x 1 ) for the array closest to the structure toe.The test with a wave period T = 5.0 s and foreshore slope 1:30 is considered.Both a highly efficient absorber and a 1:9 impermeable slope with the toe at a depth of 60 cm are shown.accurate total incident waves being predicted, and also almost identical total reflected waves.As mentioned previously, this would though not be the case if waves are extrapolated outside of the array.It is interesting to note that for this highly nonlinear wave, the present method is as accurate as for the mildly nonlinear waves This is not the case for the method for horizontal seabed which observe a significant reduction in accuracy for the highly nonlinear wave.For the 1:100 foreshore, the improvement is not as significant and both methods lead to quite small incident free higher harmonics.This is because the wave is shoaling less over the length of the wave gauge array than for the steeper foreshore.The wave on the 1:100 foreshore also includes more higher order energy as they are largest for a gentle foreshore, cf.Eldrup and Lykke Andersen (2020).Thus, the reflected wave train also contains higher free harmonics that are not well described by linear theory.Thus, the separation is less reliable, which is also shown by the form of the reflected primary wave including the bound harmonics.
The overall conclusion from the tests that include reflected waves is that accurate separation of bound and free harmonics on a sloping foreshore requires a mathematical model including shoaling.Thus, the separation is always better with the present method including the Eldrup and Lykke Andersen (2020) nonlinear shoaling model than with the existing method for horizontal seabed.The improvement is largest for steeper foreshore and not very significant on a 1:100 foreshore unless waves are extrapolated outside of the array.The separation is very reliable for all cases where free higher harmonics are so small that they may be accurately described as linear waves.Larger deviations occur if free higher harmonics become so large that they cause nonlinear interactions.For coastal structures the requirement of small reflected free harmonics would usually be fulfilled unless the incident waves are highly nonlinear and the slope angle of the tested structure is large.The deviations in the total incident and reflected wave trains is though always small independent of which method is used.This would however not be the case if waves are extrapolated outside of the measurement array.For such cases, the only accurate method to estimate the incident waves would be by performing calibration tests without the structure in place.This is though also far from a trivial task as it requires highly effective passive and active absorption systems.

Conclusion
In the present paper a new method to separate incident and reflected nonlinear regular waves over sloping bathymetries is presented.The method is a combination of the nonlinear method for separation of nonlinear regular waves over horizontal foreshores (Lykke Andersen et al. (2017)) and a nonlinear shoaling model (Eldrup and Lykke Andersen (2020)).In order to validate the method, numerical test data is generated by the MIKE 3 FM numerical model.Foreshore slopes 1:100, 1:50 and 1:30 were tested.The present separation method is first validated on tests with incident nonlinear waves only and thus a highly efficient passive absorber is present at the end of the foreshore slope.For these test cases, excellent agreement between estimated and the theoretical incident wave profile is found within the validity range of the nonlinear shoaling model.The validity ranges were quantified as functions of the wave steepness and the relative water depth.Outside of the validity range and up to breaking, the error on the incident wave height increases but is typically below 1-2%, but may be up to 4 % for highly nonlinear waves on steep foreshores.If shoaling is disregarded significant spurious free incident waves and/or reflected waves are predicted when waves become nonlinear.With the existing method, the error on the incident wave height increases to up to 12 % when waves are close to breaking and is high even for mildly sloping foreshores.These tests show that the inclusion of the nonlinear shoaling model is important, especially when considering highly nonlinear waves close to breaking.When extrapolating the measurements outside of the measurement array this becomes even more visible.For waves that are only mildly nonlinear the existing method for horizontal floor provides acceptable results.However, even for mildly nonlinear waves it is essential to include shoaling when extrapolating waves outside of the measurement array.
The new method is also validated based on test cases with structures giving from 10 % to 90 % reflection.In each of these tests, an improvement over the method for horizontal seabed is also found.Especially when the seabed is steep and waves are highly nonlinear a significant improvement is found.However, the accuracy of the present method is dependent on the amount of reflection and highly accurate results are only found when reflection is so low that reflected waves are linear or mildly nonlinear.This is because when reflection is so high that nonlinear free superharmonics are present, they will interact with each other and with the primary component and such interactions are not included in the mathematical model.To include those correctly in the mathematical model would also require many more wave gauges in the array and thus it becomes impractical.Thus further improvements are only expected to be possible by an improved shoaling and de-shoaling model.As most hydraulic model tests are performed with irregular waves, future work will of course be to implement a nonlinear shoaling model also in the method of Eldrup and Lykke Andersen (2019a).

CRediT authorship contribution statement
Thomas Lykke Andersen and Mads Røge Eldrup contributed equally to all parts of the paper.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Numerical model setups.L 1 (h) is the wavelength calculated by linear wave dispersion using the local water depth (h).

Fig. 2 .
Fig. 2.Validity limits of the shoaling model for reflection analysis given in terms of maximum acceptable wave steepness (H/L).Breaking limitation according toMiche (1944) is also shown.

Fig. 3 .
Fig. 3.Estimated incident and reflected wave heights for bound and free wave trains along the foreshore for the period T = 2.4 s using methods of Lykke Andersen et al. (2017) and present method.Vertical dashed grey lines show identified minimum depth of validity of the Eldrup and Lykke Andersen (2020) shoaling model according to Fig. 1.Dotted part for present method shows depths for which wave height had to be reduced for convergence in the stream function theory calculation of the shoaling.

Fig. 6 .
Fig. 6.Estimated incident and reflected time series for the case with T = 2.4 s for 1:100 and 1:30 foreshore slopes for two water depths.The two depths correspond respectively to an array where the shoaling model is just valid (sub figure a,c) and the array at the end of the slope (sub figure b,d).

Table 1
Lowest depth in metre where the Eldrup and Lykke Andersen (2020) shoaling model is valid.