Radiation stress introduced by Longuet-Higgins and Stewart (1964) appear in the phase-averaged momentum equation of motion describing oceanic fluid flow and can successfully explain the phenomena such as wave set-up and longshore currents observed in the surf zone (Svendsen,1984; Xia et al., 2004). The variation of the radiation stress due to shoaling and breaking is balanced in the on-offshore direction by the slope of the mean water level (set-up and set-down). In the case of wave arriving at an angle to the shore, the variation across the surf zone of the wave-induced longshore momentum flux in a direction perpendicular to the coast provides a driving force for the longshore current. With the development of the concept of radiation stress, it has also been used widely in rip current (Bowen, 1969), modeling wave-induced currents (Bao and Nishimura, 2000; Yao et al., 2020; Xia et al., 2020) and wave-current coupling (Cao and Wang, 1993; Song et al., 2020). For the calculation of radiation stress under the action of regular waves, many researchers (Longuet-Higgins and Stewart, 1964; Feddersen, 2004; Mellor, 2011) have given mature results. For the radiation stress under the action of irregular waves, Li et al. (2006) gave a comprehensive consideration of multiple deformation factors mathematical model of near-shore multi-directional irregular wave propagation. Tang et al. (2008) considered the randomness of incident waves and used JONSWAP wave spectrum to discretize the unidirectional irregular wave elements at the incident. Based on the wave potential function and other parameters of the parabolic gentle slope equation, the wave radiation stress is calculated by Zheng et al. (2000) and the numerical simulation results are verified. Compared to regular waves, irregular waves are more reflective of actual wave motion. By considering irregular waves as the superposition of a series of regular waves (small-amplitude waves), this paper proposes exact and approximate methods of calculating irregular wave radiation stress. They are results through vertical integration, which are different from the forms of vertical distribution of radiation stress (Xia et al., 2004; Mellor, 2011) and general forms based on regular radiation stress (Longuet-Higgins and Stewart, 1964; Zheng et al., 2000). The two calculation methods are analyzed and compared. On the basis of further experimental results (Shen et al., 2016), the key variable wave energy in the approximate calculation method of irregular wave radiation stress is verified. The results show that in the case of narrow-band spectrum, the approximate calculation method has good accuracy, which can save a lot of calculation time and improve the calculation efficiency.

The paper is made up as follows. Firstly, the exact calculation method of irregular wave radiation stress is derived in Section 2. Secondly, we give a approximate solution of irregular wave radiation stress in Section 3. Thirdly, the comparison between exact solution and approximate solution of irregular wave radiation stress is given in Section 4. Subsequently, we verify the key variable wave energy of radiation stress through measured wave set-up time series in Section 5. Finally, conclusions are given in Section 6.

The irregular wave radiation stress is calculated by treating the irregular wave as a series of regular waves (small-amplitude waves) in this section. First consider the one-way irregular wave, whose wave setup $\eta $ can be expressed as:

where ${\theta _i} = {k_i}x - {\omega _i}t + {\varepsilon _i}$, ${A_i}$ is the amplitude of the constituent waves, ${k_i}$ and ${\omega _i}$ are the wave numbers and wave frequencies, respectively, ${\varepsilon _i}$ is a random phase (evenly distributed between 0 and $2{\rm{\pi }}$), x is the lactation of wave surface, and t is time.

Because of $\overline {\tilde u\tilde w} \ne 0$ at the two vertical sides of control body, the average compressible stress of irregular waves should be corrected, there is

The radiation stress in the x and y directions can be obtained based on Longuet-Higgins and Stewart (1964) as follows:

where $\tilde u$, $\tilde v$ and $\tilde w$ are the components of wave water particle velocities along the three directions of x, y and z, respectively; h is water depth.

From the potential function of the wave water particle velocity

yields

From the formula of the sum of two complex numbers ${\beta _1}$ and ${\beta _2}$ (where “$*$” means complex conjugate)

yields

Substituting Eqs (9)–(11) into Eqs (3) and (4) gives

where

For the traveling waves incident at an angle $\alpha $, the corresponding radiation stress tensor can be obtained by a two-dimensional tensor coordinate transformation (Zou, 2009).

The above exact calculation method can consider each component wave of the irregular wave separately, but it will take a large amount of calculation because of cumulative summation of each component wave for each time step. Here, it is simplified and approximated in order to achieve a result with suitable calculation speed and accuracy. This method is for narrow-band wave energy spectrum, that is, the spectrum $S(\omega)$ has a large value only near the peak frequency ${\omega _{\rm{p}}}$, as shown in Fig. 1, which has obvious wave group characteristics for wave surface.

The amplitude ${A_i}$ and frequency ${\omega _i}$ of the constituent waves can be determined by the wave energy spectrum $S(\omega)$,

The frequency and amplitude of each component wave are determined by the center frequency and area of strips with a width of $\Delta \omega $ divided from the wave energy spectrum, as shown in Fig. 1. The division here is equidistant, $\Delta \omega $ is a constant. In order to avoid periodic repetition of the frequency of the simulated component waves, the frequency is randomly selected as the representative frequency of the interval within each frequency division interval.

After the water depth h of the wave area is known, the wave number k of the component wave can be calculated from the following dispersion equation:

The irregular wave is regarded as the superposition of a series of regular waves (small-amplitude waves). From Eq. (1), the time series of wave setup can be obtained, and the envelope expression of the time series of wave setup $\eta (x,t)$ can be obtained as:

in which

From Eq. (20),

the accurate values of wave energy at various locations on flat-bottomed terrain over time can be obtained. However, when the topography changes, the wave energy at each position in the above formula will no longer maintain the form of the analytical solution. At this time, the wave energy ${E_{\rm{w}}}(t)$ calculated at the origin of the above formula with time can be used as the incident boundary condition of the wave energy. The wave group propagates forward at the propagation speed ${c_{\rm{g}}}$ corresponding to the peak frequency (because of the narrow-band wave energy spectrum), and the wave energy conservation Eq. (21) can be used to obtain the approximate value of the wave energy ${E_{\rm{w}}}(x,t)$ at any time and position.

where ${c_{\rm{g}}}$ is the propagation speed of the wave group, $\alpha $ is wave incident angle, and ${S_{\rm{d}}}$ is wave energy dissipation.

The radiation stress at each position at each time can also be calculated approximately according to the results of regular waves. For traveling waves incident at an angle $\alpha $, there is

where $n = (1 + 2{k_{\rm{p}}}h/\sinh 2{k_{\rm{p}}}h)/2$, is wave energy transfer rate. k_p is the wave number corresponding to the peak frequency.

From the formula for calculating the irregular wave radiation stress, compared to the exact calculation method of the irregular wave radiation stress, there are three approximations to the results of radiation stress obtained by this approximation method: one approximation from Eq. (22) using the regular wave result itself, one approximation from a calculation of wave energy ${E_{\rm{w}}}(x,t)$ at each position and each time in Eq. (22), it is assumed that the wave energy at the initial time propagates forward at the propagation speed ${c_{\rm{g}}}$ corresponding to the peak frequency, another place comes from the wave energy transfer rate $n$ in the expression of radiation stress calculation in Eq. (22), which is approximately calculated from the wave number corresponding to the peak frequency. This method replaces the frequency of each component wave with the peak frequency, and has better accuracy only in the case of a narrow-band wave energy spectrum. The advantage of this method is that the irregular wave is calculated approximately as regular wave, the original regular wave result can be directly applied, and the calculation speed is much faster than the exact calculation method.

In order to analyze the calculation accuracy and applicable range of the approximation method of irregular wave radiation stress, an improved JONSWAP type spectrum is used to generate irregular waves, and its expression (Hasselmann et al., 1973; Goda, 1999) is as follows:

where $f$ is wave frequency, ${f_{\rm{p}}}$ is peak frequency, ${T_{\rm{p}}}$ is peak period, ${T_{H_{1/3}}}$ is effective period, ${H_{1/3}}$ is effective wave height, ${\text{γ}} $ is peak factor, and $\sigma $ is peak shape factor ($f \leqslant {f_{\rm{p}}}$, $\sigma $=0.07; $f > {f_{\rm{p}}}$, $\sigma $=0.09).

The advantage of this spectral formula is that once the value ${\text{γ}} $ is selected, the spectral shape can be determined by the designed wave element. Here, the commonly used spectral peak factor ${\text{γ}} = 3.3$ is selected for further analysis and comparison. Other parameters ${T_{H_{1/3}}}$ is 1.5 s, ${H_{1/3}}$ is 4.49 cm, water depth $h$ is 0.045 m and the angle between the wave incident direction and the normal of shoreline $\alpha $ is 30°. The corresponding spectrum chart is shown in Fig. 2. The influence of the peak factor on the shape of the spectrum is also shown in Fig. 2 at the results of ${\text{γ}} $=5.0 and 8.0.

Figure 3 shows the comparison between the accurate calculation results and the approximate calculation results of the wave energy of the irregular waves generated by the above JONSWAP spectrum at x=0 m, x=5 m, and x=10 m, respectively. Figure 4 shows the comparison between the exact calculation results and the approximate calculation results of the radiation stress (S_xx, S_yy, S_xy) of the irregular wave in the three directions at the corresponding position. In order to show the main source of the error, Fig. 4 also gives the results of irregular wave radiation stress calculated using accurate wave energy which can be obtained by Eq. (20).

To measure the closeness between the approximate calculation result and the accurate calculation result, the figures also show the statistical parameter d proposed by Willmott (1981) to judge, the value of parameter d is between 0 and 1, and the value equal to 1 means that the two are completely matching, the closer the value is to 1, the closer the two results are. The calculation expression is as follows:

where $\alpha (i)$ and $\beta (i)$ are the corresponding results of accurate calculation and approximate calculation method respectively, $\bar \alpha $ represents the average value of the results obtained by the accurate calculation.

It can be seen from Fig. 3 that wave energy of irregular wave calculated by approximate method is close to that by accurate method, especially not far away from the initial propagation position. The greater the distance, the greater the relative error. From Fig. 4, we can see that the wave radiation stresses of irregular wave in three directions calculated by approximate method is close to that by exact method. The main error comes from the wave energy obtained by approximate method, which can be easily seen from the statistical parameter d. When obtained the accurate wave energy, the corresponding results agree much better. Further comparison will find that the exact calculation method can more accurately reflect the fluctuation of the radiation stress at each position at each time, which is of great significance for the further study of the unstable movement of longshore currents.

Since we cannot directly compare the irregular wave radiation stress through experiments, here we can verify the reliability of the results obtained by the wave energy approximation method through experiments, so that the reliability of the irregular wave radiation stress calculated by the approximation method can be deduced by the above section’s conclusion.

The measurements were performed in the wave basin of the State Key Laboratory of Coastal and Offshore Engineering in Dalian University of Technology, with the dimensions of 55 m in length, 34 m in width, and 0.7 m in depth. Slope of 1:40 was selected in the present study to verify the reliability of the results obtained by the wave energy approximation method.

Figure 5 showed the experimental setup. The origin of the coordinate system (x, y) located at the right side of the upstream. The offshore direction was set as the positive direction of the x-axis, and the downstream direction was the positive direction of the y-axis. The centreline of the beach model was tilted at 30° in respect to the wave-maker board. In order to ensure that the submerged slope length of both slopes was equal to 18 m, and the corresponding water depths in front of the slopes was set as 0.45 m. The distances between the slope toe and wave maker board ranged from 8.0 m to 22.5 m. Two waveguide walls were installed between the wave maker and slope toe, in order to prevent the influence of currents outside the experimental area. The wave absorbers, which was constructed with net boxes filled with plastic scraps, were attached to the inner side of the waveguide walls to prevent wave reflection. An additional 3 m long slope section was reserved above the water surface, thereby allowing sufficient space for wave run-up. A distance of 4.4 m was reserved between the model and walls to eliminate the boundary effects. During the experiment, currents were generated by the water level difference caused by the wave motion, which finally formed a circular flow in the basin.

The wave height was measured by 60 capacitive gauges placed in three different cross-sections (I to III) normal to the coastline (as shown in Fig. 5). The distance from the origin O to cross-section I was 7.0 m, while the interval between every two adjacent sections was 5 m. Both cross-sections I and III had 14 wave gauges which were allocated with an interval of 2.0 m. However, cross-section II had 32 wave gauges which were arranged with smaller intervals (0.5 m and 1.0 m for the ranges of 0.5–10 m and 10–22 m, respectively) to acquire the location of breaking point more accurately. Figure 6 showed the configuration of the velocity meters and wave gauges.

Here we compare the wave energy result converted from the wave setup measured by the wave gauges at the four positions (x=22 m, 20 m,10 m, 3 m) in the middle row of the experiment with the approximate wave energy calculation result. Among them, x=22 m is the location of the first wave gauge in the middle column, which is located at the flat bottom of the wave basin. The wave energy at this location will also be used as the incident boundary condition of the approximate calculation method of wave energy; x=20 m is located at the foot of the slope; x=10 m is located outside the wave breaking zone of the slope; x=3 m is located in surf zone of the slope. Figure 7 gives the results of wave setup time series at the above four locations for irregular wave case with T_p=1.5 s and H_rms=4.49 cm. It should be pointed out that the irregular wave in the experiment is generated by JONSWAP spectrum, which is a narrow-band wave energy spectrum mentioned above.

Before wave energy comparison between experimental results and approximate calculation results, a method of determining the wave energy corresponding to the wave setup time series $\eta (t)$ obtained from experimental measurements is given. $\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\eta } (t)$ can be obtained by Hilbert transform of the above $\eta (t)$. The wave envelope can be got by the relation of $A(t) = {({\eta ^2}(t) + {\overset{\lower0.2em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\eta } ^2}(t))^{1/2}}$. The wave energy changing with time ${E_{\rm{w}}}(t)$ at the location of the wave gauge can be deduced by ${E_{\rm{w}}}(t) = \dfrac{1}{2}\rho g{A^2}(t)$.

In the approximate numerical calculation, the incident boundary is directly taken from the time-varying wave energy ${E_{\rm{w}}}(t)$ obtained by converting the wave setup time series $\eta (t)$ measured at x=22 m. According to the approximate processing method in Section 2, it is assumed that the wave group propagates forward at peak-frequency ${\omega _{\rm{p}}}$, the wave energy ${E_{\rm{w}}}(t)$ can be calculated at any position at any time from the wave energy conservation Eq. (21).

Figure 8 shows the comparison of wave energy time series after filter (low-pass filter cutoff frequency is 0.1 Hz) of the above-mentioned irregular wave conditions (T_p=1.5 s, H_rms=4.49 cm) at different positions between approximate method and experimental results deduced by the wave setup measured by the experiment after Hilbert transformation.

The first line in Fig. 8 is the wave energy at 22 m from the shoreline, and is also the wave energy condition at the incident boundary, and the two coincide. The second line is the wave energy at 20 m from the shoreline, which is in a flat terrain area. Comparing the experimental results and numerical results, it is found that the two are in good agreement, the corresponding peak frequency positions are basically coincident, and the energy amplitude is also in good agreement; the third line is the wave energy at a distance of 10 m from the shoreline, which is in the plane slope terrain field, and the wave is not broken at this location. Comparing the experimental results and numerical results shows that the agreement between the two is also good, but the agreement is slightly worse than the flat-bottomed terrain field, which shows that the corresponding peak frequency position is slightly offset, and the energy magnitude deviation is also larger than the flat-bottomed terrain field; the fourth line is the wave energy at a distance of 3 m from the shoreline, which is also in the plane slope terrain field, but the wave has been broken at this location. Comparing the experimental and numerical results, it is found that the deviation between the two is large at this time. The high-frequency peak position and energy magnitude of the wave energy in the experimental results have been difficult to correspond with each other in the numerical calculation results, but the numerical calculation results can basically reflect the average wave energy in the experimental results. In general, the numerical calculation results are in good agreement with the experimental results, especially before the wave breaking, which indicates that the wave energy approximate calculation method assuming that the wave group propagates forward at the wave group velocity corresponding to the peak frequency is reasonable. Figure 8 also gives the calculation result of the corrected propagation frequency $\,{f_{\rm{p}}} = 0.5$ Hz to reduce the error by correcting the wave group velocity to propagate toward the shore at a wave group velocity close to the peak frequency. By comparison, it can be found that the value of d after correction before wave breaking is larger than the result before correction, which makes the result after correction slightly improved, but it fails to improve after wave breaking. This is because the wave breaking is more complicated, and the energy dissipation model fails to accurately reflect the wave breaking progress, so that the numerical calculation results after the wave breaking are not accurate enough.

In this paper, the exact calculation method and approximate calculation method of the radiation stress under the action of irregular waves are given. The results of the two calculation methods for wave energy and irregular wave radiation stress are analyzed and compared. The experimental results are used to further verify the wave energy calculation in the approximate calculation method. The results show that in the cases of narrow spectrum, the approximate calculation method has good enough accuracy, which can save a lot of calculation time, thereby greatly improving the calculation efficiency. However, the exact calculation method can more accurately reflect the fluctuation of the radiation stress at various positions at all times, which is of great significance for the further study of the instability of longshore currents.