1. Introduction
Prediction of beach erosion and coastline changes requires understanding of the triggering mechanisms for sediment motion in surf and swash zones, which depend strongly on the flow characteristics (Chardón-Maldonado, Pintado-Patiño & Puleo Reference Chardón-Maldonado, Pintado-Patiño and Puleo2016). Many previous sediment transport studies have focused on the shear stress and pressure gradient in oscillatory flows and breaking waves (Cox, Kobayashi & Okayasu Reference Cox, Kobayashi and Okayasu1996; Ribberink Reference Ribberink1998; Nielsen Reference Nielsen2002; Sumer et al. Reference Sumer, Guner, Hansen, Fuhrman and Fredsøe2013; Frank et al. Reference Frank, Foster, Sou and Calantoni2015; Ting & Kern Reference Ting and Kern2022). Most numerical models can predict morphological changes with reasonable accuracy in the nearshore region where the water depth exceeds about 2 m (Ruessink et al. Reference Ruessink, Kuriyama, Reniers, Roelvink and Walstra2007; Brinkkemper et al. Reference Brinkkemper, Lanckriet, Grasso, Puleo and Ruessink2016), but the gap between model predictions and laboratory/field observations becomes significant in inner surf and swash zones, where the water depth is shallower. This is because the effects of broken bore generated turbulence on sediment suspension and transport are often ignored in the numerical models (Masselink & Puleo Reference Masselink and Puleo2006; Ruessink et al. Reference Ruessink, Kuriyama, Reniers, Roelvink and Walstra2007; Brinkkemper et al. Reference Brinkkemper, Lanckriet, Grasso, Puleo and Ruessink2016), which actually plays an important role for keeping the sediments suspended in the sediment transport processes (Puleo et al. Reference Puleo, Beach, Holman and Allen2000; Butt et al. Reference Butt, Russell, Puleo, Miles and Masselink2004). Some recent studies of a plunging solitary wave and plunging periodic waves indicated that the size of the eddies generated in the flow are proportional to the sand sediment suspension (Sumer et al. Reference Sumer, Sen, Karagali, Ceren, Fredsøe, Sottile, Zilioli and Fuhrman2011, Reference Sumer, Guner, Hansen, Fuhrman and Fredsøe2013). Moreover, the effects of wave–wave interaction on the sediment transport are also usually overlooked because of the lack of turbulence information. A better understanding of the turbulence structure during the dynamic interaction between uprush and downwash swash flows is essential in improving sediment transport modelling in the surf and swash zone.
Many laboratory experiments have been conducted to study the surf and swash flow induced by periodic waves. The typical phase-averaging method is applied to the time series of velocity measurements in a single experimental run so as to obtain the turbulent velocities based on the Reynolds decomposition concept (Ting & Kirby Reference Ting and Kirby1995, Reference Ting and Kirby1996; Cox & Kobayashi Reference Cox and Kobayashi2000; Petti & Longo Reference Petti and Longo2001; Cowen et al. Reference Cowen, Sou, Liu and Raubenheimer2003; De Serio & Mossa Reference De Serio and Mossa2006; Kimmoun & Branger Reference Kimmoun and Branger2007; Sou, Cowen & Liu Reference Sou, Cowen and Liu2010; Sou & Yeh Reference Sou and Yeh2011; De Serio & Mossa Reference De Serio and Mossa2019), the resulting turbulent intensities are likely overestimated if the wave record in a single run is not perfectly periodic and the errors are difficult to estimate (Sou et al. Reference Sou, Cowen and Liu2010). Moreover, if the wavemaker is not able to absorb all the reflected waves from the beach, the waves in the wave flume are no longer pure progressive waves, making the interpretation of the resulting wave–wave interactions in the swash flows difficult. Examining the turbulence and energy dissipation of spilling waves, Huang et al. (Reference Huang, Hsiao, Hwung and Chang2009) attempted to avoid the variations between waves by ensemble averaging one particular wave in a train of periodic waves from multiple test runs. They found that the variation of breaking point for the same wave caused by the long wave effect were problematic.
On the other hand, the methodology for generating a solitary wave has been greatly improved so that a sequence of solitary waves can be accurately generated with high repeatability (Lo, Park & Liu Reference Lo, Park and Liu2013; Pujara, Liu & Yeh Reference Pujara, Liu and Yeh2015b). Most recently, multiple solitary waves, up to nine, were generated in a wave flume for studying runup processes both experimentally and numerically (Wu, Higuera & Liu Reference Wu, Higuera and Liu2021). In this paper, the same experimental set-up is used with the new focus on measuring the swash flows. Since the consecutive solitary waves are highly repeatable, the same experiment can be repeated as many times as one wishes, allowing turbulence measurements with the ensemble-average procedure (Ting Reference Ting2006, Reference Ting2008), yielding more accurate turbulence measurements. This paper can be viewed as the sequel to Wu et al. (Reference Wu, Higuera and Liu2021) with the specific goal to examine the turbulence structures in terms of the length scale evolution of both weakly and strongly interacting cases as defined in Wu et al. (Reference Wu, Higuera and Liu2021).
In terms of velocity measurements within surf and swash zones, the laser Doppler velocimetry (LDV), featuring single-point measurements, was employed in many earlier works (Stive Reference Stive1984, Reference Stive1985; Svendsen Reference Svendsen1987; Nadaoka, Hino & Koyano Reference Nadaoka, Hino and Koyano1989; Ting & Kirby Reference Ting and Kirby1995, Reference Ting and Kirby1996; Pedersen, Deigaard & Sutherland Reference Pedersen, Deigaard and Sutherland1998; Cox & Kobayashi Reference Cox and Kobayashi2000; Petti & Longo Reference Petti and Longo2001; Stansby & Feng Reference Stansby and Feng2005; Shin & Cox Reference Shin and Cox2006; De Serio & Mossa Reference De Serio and Mossa2006, Reference De Serio and Mossa2019). However, the particle image velocimetry (PIV) measurement technique has become a favourable choice for measuring the whole-field flow velocity and turbulence structures because it can provide highly resolved temporal and spatial information (Adrian Reference Adrian1984, Reference Adrian1991). Over the past few decades, PIV has been used to study breaking waves and turbulence dissipation despite of the difficulty in imaging aerated regions associated with wave breaking (Lin & Rockwell Reference Lin and Rockwell1994, Reference Lin and Rockwell1995; Dabiri & Gharib Reference Dabiri and Gharib1997; Chang & Liu Reference Chang and Liu1998, Reference Chang and Liu1999; Melville, Veron & White Reference Melville, Veron and White2002; Cowen et al. Reference Cowen, Sou, Liu and Raubenheimer2003; Govender, Mocke & Alport Reference Govender, Mocke and Alport2004; Kimmoun & Branger Reference Kimmoun and Branger2007; Drazen & Melville Reference Drazen and Melville2009; Huang et al. Reference Huang, Hsiao, Hwung and Chang2009; Sou et al. Reference Sou, Cowen and Liu2010). To complement the PIV measurement in the aerated flow, Lim et al. (Reference Lim, Chang, Huang and Na2015), Na et al. (Reference Na, Chang, Huang and Lim2016) and Na, Chang & Lim (Reference Na, Chang and Lim2020) used a combination of PIV and bubble image velocimetry (BIV) technique to investigate the flow and turbulence properties of plunging and spilling breaking waves in deep water. They found that the turbulence properties related to the bubble-aerated region were overestimated when void fraction was not considered. Since BIV uses bubbles as the tracer particles for evaluating the flow velocity, BIV would work well if the flow velocity is larger than the bubble buoyancy velocity. During the uprush–downwash flows, the periodic forcing might not always satisfy this condition.
As mentioned earlier, the size of the eddies generated in the flow are correlated to the sand sediment suspension (Sumer et al. Reference Sumer, Sen, Karagali, Ceren, Fredsøe, Sottile, Zilioli and Fuhrman2011, Reference Sumer, Guner, Hansen, Fuhrman and Fredsøe2013). This motivates the examination of the integral length scale as it represents the largest eddy size in the turbulent flows (Pope Reference Pope2000). Many previous studies evaluated the integral length scale in various ways, e.g. Ting & Kirby (Reference Ting and Kirby1996), Petti & Longo (Reference Petti and Longo2001), Govender et al. (Reference Govender, Mocke and Alport2004), De Serio & Mossa (Reference De Serio and Mossa2019) and Sou et al. (Reference Sou, Cowen and Liu2010). One of the methods to estimate the integral length scale is to use the scaling relationship between the dissipation rate 
$\epsilon$ and the turbulent kinetic energy 
$k$ as 
$\mathcal {L}=c_dk^{3/2}/\epsilon$, in which 
$c_d$ is an empirical coefficient (Govender et al. Reference Govender, Mocke and Alport2004). The dissipation rate can be evaluated by fitting the 
$-$5/3 inertial subrange of the energy spectrum by assuming homogeneous isotropic turbulence (Pope Reference Pope2000; Govender et al. Reference Govender, Mocke and Alport2004). For swash flows, the assumption of homogeneous isotropic turbulence does not hold when the flow depth becomes shallow and when the location of the spectrum is close to the bed (Sou et al. Reference Sou, Cowen and Liu2010). Alternatively, the integral length scale can be estimated from the time-averaged integral length scale, which can be calculated from the measured time histories of velocity field (Ting & Kirby Reference Ting and Kirby1996; Petti & Longo Reference Petti and Longo2001; De Serio & Mossa Reference De Serio and Mossa2019). The integral time scale is then converted into the integral length scale by invoking the Taylor frozen turbulence hypothesis (Taylor Reference Taylor1938). The resulting integral length scale from this method is independent of time, which is not realistic because the eddy size actually varies in time within a wave cycle. On the other hand, using PIV spatial data for periodic waves, Sou et al. (Reference Sou, Cowen and Liu2010) evaluated the time-dependent integral length scale within a wave cycle directly from the turbulent kinetic energy spatial spectra. They found that the integral length scale was within the range of 10–20 % of the local water depth within a wave cycle in the surf zone, but the integral length scale increased up to 40 % of the local water depth in the swash zone. These results suggested that the vertical confinement due to the shallower water depth in the swash zone would constrain the large-scale turbulence structure (eddy structure) to horizontal motions.
In this study, a 5-m long-stroke wavemaker is used to generate six consecutive solitary waves. By repeating the same experiment 10 times, the ensemble averaging is then performed to determine the ensemble-averaged velocities and the turbulent velocities. A high-speed PIV (HSPIV) system is employed to obtain the two-dimensional (2-D) velocity measurements during the interactions of the uprush and downwash flows between consecutive solitary waves. The six consecutive solitary waves produce five uprush–downwash interactions in a single experimental test. Two solitary wave conditions were examined: (1) the leading solitary wave corresponds to a non-breaking wave with five subsequent surging breaking waves (NB case) and (2) the leading solitary wave is a surging breaking wave with five subsequent plunging breaking waves (SB case).
The HSPIV experimental set-up and the generation of solitary waves are presented in § 2. The distinctive flow behaviour is explained with the PIV raw images and the resolved instantaneous velocity fields in § 3. Then, the ensemble-averaged velocity and the variation of the Froude number are presented to gain insight of the processes during the uprush–downwash interaction in the surf zone. The differences and similarities of the flow and turbulence characteristics among the five interactions as well as between the two wave conditions are discussed in §§ 3.2 and 3.3. More discussions on the turbulent energy intensity are presented in § 3.4. The turbulence structures are investigated using the spatial spectral analysis in § 3.5. Lastly, the integral length scale and eddy viscosity are discussed in §§ 4.1 and 4.2. In § 5 concluding remarks are provided. The identification of free surface is described in Appendix A. The convergence test for the turbulence statistics and the uncertainty analysis are presented in Appendices B and C, respectively. Additional videos and information for the five interactions are presented in the supplementary material available at https://doi.org/10.1017/jfm.2023.484.
2. Laboratory experiments
2.1. Experimental set-up and wave conditions
Laboratory experiments were conducted in a wave flume (36 m long, 0.9 m high and 0.9 m wide) in the hydraulics laboratory at the National University of Singapore. The experimental set-up is the same as that reported in the companion paper, Wu et al. (Reference Wu, Higuera and Liu2021) (hereafter referred to as Reference Wu, Higuera and LiuWu21). For clarity the set-up is sketched in figure 1. The wave flume is equipped with a 5-m long-stroke wavemaker at one end and a 1/10 glass beach on the other end. The 
$x$ direction is positive in the on-shore direction along the beach face, the 
$z$ direction is pointing upwards normal to the beach face and the 
$y$ direction is set by the right-hand rule. Finally, the origin of the coordinate system is located at the still water shoreline (figure 1).
Figure 1. Experiment set-up and coordinate system (
$x$–
$z$ plane).
Following Goring (Reference Goring1978) and Lo et al. (Reference Lo, Park and Liu2013), successive solitary waves are generated by combining the trajectories of the piston-type wavemaker for six individual solitary waves. The third-order analytical solution for the free surface profile of a solitary wave (Grimshaw Reference Grimshaw1971) is applied. Using the offshore constant water depth, 
$h = 24$ cm, two solitary wave conditions are examined: the ratio of the incident wave height and water depth 
$H_o/h = 0.11$ (case 1) and 
$H_o/h = 0.22$ (case 2), representing weak and strong interaction, respectively, as discussed in Reference Wu, Higuera and LiuWu21. The separation time between two successive solitary wave crests are 
$T = 3.52$ and 2.65 s for the weakly and strongly interacting cases, respectively, which are approximately equal to the effective periods of individual solitary wave (Goring Reference Goring1978). Following the definition outlined in Lo et al. (Reference Lo, Park and Liu2013) and Pujara, Liu & Yeh (Reference Pujara, Liu and Yeh2015a), the corresponding effective wavelengths, 
$\lambda$, are 5.25 and 3.71 m, respectively. With the 5-m long stroke, the wavemaker can generate six successive solitary waves without the interference of the re-reflection from the wavemaker.
Figure 2 shows the time history of the six solitary waves at the wave gauge located at 10 m from the wavemaker. Multiple runs for each wave condition were performed and the free surface elevation measurements collapsed perfectly in a single line. It is apparent that the wavemaker is capable of generating repeatable successive solitary waves. The surf similarity parameter 
$\xi$ for periodic waves can be defined as 
$\xi =s/\sqrt {H_o/({gT^2}/{2{\rm \pi} })}$, where 
$s$ is the beach slope and 
$g$ is the gravity acceleration. The surf similarity parameters 
$\xi$ are found to be 2.74 and 1.47 for the weakly and strongly interacting cases, respectively. Using the solitary wave breaking criterion discussed in Grilli, Svendsen & Subramanya (Reference Grilli, Svendsen and Subramanya1997), the leading solitary waves of the two wave conditions correspond to non-breaking and surging breaking, respectively. As discussed in Reference Wu, Higuera and LiuWu21, for the NB case the subsequent waves become surging breakers. On the other hand, for the SB case, the subsequent waves become plunging breakers. These characterisations are consistent with the HSPIV measurements to be discussed herein.
Figure 2. Successive six solitary waves measured at the location 10 m from the wavemaker: (a) 
$H_o/h = 0.11$ (NB); and (b) 
$H_o/h = 0.22$ (SB). Green line, repeated measurements; black line, mean value. (They are not visually separable.)
While the overall wave shoaling, breaking and runup processes have been discussed in terms of the free surface profiles in Reference Wu, Higuera and LiuWu21, the focus of this paper is to discuss the corresponding turbulence evolution in the swash zone.
2.2. The HSPIV image acquisition
The 2-D cross-shore velocity measurements are acquired by a HSPIV system, consisting of an 8-W 532-nm continuous laser, a high-speed Phantom LAB340 camera and a Tokina AT-X Pro Macro 100 mm lens, laser sheet optics, a mirror to redirect the laser sheet and a data acquisition computer. The laser beam is shaped into a sheet by the laser optics and the light sheet is delivered to the measurement location from the bottom of the glass beach using a mirror beneath the beach. The high-speed camera is positioned by the wave flume taking images through the sidewall with the field of view (FOV) parallel to the sidewall, recording images at 1000 Hz. Because of the camera and camera lens combination, the optimal location for the measurement area is 15 cm from the wall. The entire water depth is captured throughout the measurements. The flow is seeded with hollow glass spheres (specific gravity 
$sg=1.10$, and mean particle diameter 
$D_p=10$ 
$\mathrm {\mu }$m).
The wavemaker and the HSPIV system are synchronised in such a way that the camera is triggered by the wavemaker with a predetermined time in order to skip unnecessary image recording before the arrival of wave. Due to the limit of the camera storage capability, the HSPIV images of the whole runup and rundown processes of all six consecutive solitary waves are captured in six separate image sets for each realisation, yielding five uprush–downwash interactions. Furthermore, each wave condition is repeated between 10 and 13 times and an ensemble average of the velocity fields is obtained from these realisations.
2.3. Data processing
The HSPIV velocity field is obtained from every two consecutive images of the whole image set (i.e. correlation between the first and second images, between the second and third images and so forth), yielding the 2-D velocity field, 
$U$ and 
$W$. The instantaneous velocity fields is evaluated using the free license software PIVlab (Thielicke & Stamhuis Reference Thielicke and Stamhuis2014). The PIV raw images are preprocessed using the high-pass filter before PIV processing. The sequence of each image set is initially interrogated with a 
$128\times 128$ pixel subwindow using a 50 % overlap. The results are further analysed with a 
$64\times 64$ subwindow with a 50 % overlap. The final interrogated subwindow is 
$32\times 32$ with a 50 % overlap, yielding a 
$159\times 55$ grid of velocity vectors for the NB case and a 
$159\times 74$ grid of velocity vectors for the SB case. A local median filter is applied to the final results to remove the invalid vectors and the missing data are obtained by interpolating the valid data points. In order to capture the entire water depth and maximise the camera spatial resolution, the PIV image size is 
$12.00\ {\rm cm}\times 35.11\ {\rm cm}$ for the NB case and 
$15.17\ {\rm cm}\times 32.83\ {\rm cm}$ for the SB case, yielding the spatial resolution of 0.222 and 0.208 cm for NB and SB cases, respectively.
Following Barranco & Liu (Reference Barranco and Liu2022), the instantaneous free surface is identified directly from the raw images. As shown in figures 3 to 4, the free surface profiles can be identified relatively easily during the later phase of the uprush flows and the beginning phase of the downwash flows as the free surface is relatively smooth (free of breaking and bubbles) and can be clearly identified in the raw images (panels a–f of the figures). When wave breaking occurs, i.e. the later phase of the downwash flows and the earlier phase of uprush flows, the free surface location is no longer easily identifiable in the raw images. The free surface shown in the instantaneous velocity vector fields (panels g–l of figures 3 to 4) is a smooth line due to the results of interpolating the identified points from the raw image in time and space. The blurry and violent aerated regions have been excluded in the HSPIV processing. However, the bubbles generated by wave breaking can still be seen underneath the free surface (see the HSPIV raw image videos in the supplementary material). Thus, the resulting velocity vector fields are biased toward the bubble vertical motions (Lim et al. Reference Lim, Chang, Huang and Na2015). As a result, the turbulent kinetic energy in the aerated region could be overestimated. We shall report the analysed results in the aerated region with cautions. Detailed discussions on the procedure of free surface identification are presented in Appendix A.
Figure 3. Raw HSPIV images (a–f) and instantaneous velocity fields (g–l) for NB45. The red vertical line indicates the location for the Froude number 
$Fr=1$. The velocity vector scale is the phase speed 
$c=\sqrt {g(h+H_o)}$.
Figure 4. Raw HSPIV images (a–f) and instantaneous velocity fields (g–l) for SB45. The red vertical line indicates the location for the Froude number 
$Fr=1$. The velocity vector scale is the phase speed 
$c=\sqrt {g(h+H_o)}$.
To synchronise the ensemble of measurements, the effective wave period 
$T$ is used for normalisation and the dimensionless time 
$t/T=0$ is determined for each interaction as the instant that the swash flow on the farthest offshore (left) side of the FOV just changes from the downwash phase to the uprush phase. Each interaction is examined starting with the downwash flow phase and ending with the uprush flow phase of the next incoming wave. Thus, the time sequence of each interaction is from the largest negative value of 
$t/T$ to the largest positive number, implying that the negative 
$t/T$ values represent the downwash flow phase of the previous wave and the positive 
$t/T$ values denote the uprush flow phase of the next incoming wave.
We adopt the Reynolds decomposition as
\begin{equation} U=\langle{U}\rangle+{u}, \quad W=\langle{W}\rangle+{w}, \end{equation}
where the ensemble-averaged velocity is obtained as the average of the instantaneous quantity over an ensemble of 
$M$ samples:
\begin{equation} \langle{U} \rangle = \frac{1}{M} \sum_{i=1}^M {U}_i, \quad \langle{W} \rangle = \frac{1}{M} \sum_{i=1}^M {W}_i, \end{equation}
and 
${u}$ and 
$w$ denote the turbulent fluctuating velocity components of the 
$i$th sample. And the magnitudes of the turbulent velocity intensity are defined as the root mean square of the turbulent fluctuating velocities 
$\sqrt {\langle u^2\rangle }$ and 
$\sqrt {\langle w^2\rangle }$ (Pope Reference Pope2000).
3. Results
3.1. Characteristics of the wave shoaling and breaking
Using HSPIV raw images and instantaneous velocity fields, general descriptions of the shoaling, breaking and swash flow processes will be provided in this section. To differentiate the five uprush–downwash interaction events, generated by two successive solitary waves, the following labels are used: NB
$IJ$ (SB
$IJ$) for the interaction process between the downwash of the 
$I$th solitary wave and the uprush of the 
$J$th solitary waves, where 
${I}=1,2,3,4,5$ and 
${J}=2,3,4,5,6$.
3.1.1. The NB case (weakly interacting case)
For the NB case, the uprush of the first wave does not break. Breaking occurs as a non-stationary hydraulic jump during the downwash flow for all six waves. In figure 3, snapshots of raw HSPIV images (a–f) and the corresponding instantaneous velocity vectors (g–l) during NB45 are plotted. The videos of the HSPIV images and the instantaneous velocity field for both NB12 and NB45 are provided in the supplementary material (movies 1–4). During the earlier downwash phase (at 
$t/T < -0.284$ and 
$-$0.171 in figure 3), the free surface is smooth and the velocity is accelerating on the landward (right) side of the FOV. Note that the downwash of the fourth wave moves in the negative 
$x$ direction along the beach face (from right to left). At 
$t/T=-0.114$ in figure 3, the breaking becomes visible on the free surface at around 
$x=-42$ to 
$-$30 cm (
$x/\lambda =-0.08$ to 
$-$0.057) on the raw image. The downwash flow (on the landward side of the FOV) becomes supercritical, 
$F_r >1$, as the water depth decreases and the velocity increases. The time-dependent Froude number is defined as
\begin{equation} Fr=\langle U \rangle_d/\sqrt{g h_L}, \end{equation}
where 
$\langle U \rangle _d$ is the depth- and ensemble-averaged velocity in the 
$x$ direction and 
$h_L$ is the time-dependent local water depth. The red vertical line shown in the instantaneous velocity vector field indicates the location of the critical flow at 
$Fr=1$ (figure 3), which is called the critical flow point (CFP). As the flow makes the transition from supercritical to subcritical, a hydraulic jump could occur. A detailed discussion on the temporal and spatial variations of the Froude number together with the ensemble-averaged velocity will be presented in § 3.3 in the context of possible occurrence of hydraulic jump. While the free surface continues to break (see 
$t/T=-0.057$ in figure 3), the flow velocity near the left (offshore) side of the FOV, which is increasingly under the influence of the incoming second solitary wave, becomes weaker and the water depth deeper, and the flow remains subcritical. During this period of time, the flow is moving in the offshore direction, but is transitioning spatially from supercritical to subcritical. As the supercritical flow depth subsides and the subcritical flow velocity decreases and reverses its direction, a broken bore is formed. At around 
$t/T=0.000$ when the subcritical flow velocity on the far left (offshore) side of the FOV is further reduced and is close to zero, the broken bore barely moves for a short period of time (the temporary ‘stationary’ broken bore is further discussed in § 3.3). The uprush flow from the second solitary wave then pushes the broken bore up the slope. The uprush flow eventually becomes unidirectional as shown in the instantaneous velocity field at 
$t/T=0.114$ (figure 3).
The eddies generated from the bed during the downwash flow can be traced by the white seeding particles in the HSPIV images between 
$x=-40$ and 
$-$35 cm or 
$x/\lambda =-0.08$ and 
$-$0.067 (see movie 1 in the supplementary material). The bed-generated eddies are eventually merged within the shoreward propagating surface roller.
The intensities of the first pair of uprush–downwash interactions, i.e. NB12 is strongest among all the interactions. The violent splash seen in NB12 (
$t/T=-0.028$ in movie 1) does not occur in the next four interactions even though the non-stationary hydraulic jump and broken bore also develop during the downwash phase (see § 3.3). More importantly, the HSPIV measurements reveal that the velocity fields of NB34, NB45 and NB56 are very similar. The videos of the velocity fields are included in the supplementary material. This is consistent with the observation made in Reference Wu, Higuera and LiuWu21 that the corresponding runup heights became the same after the third wave (see figure 17a in Reference Wu, Higuera and LiuWu21). For brevity, only the HSPIV raw images and the corresponding instantaneous velocity fields of NB45 are shown in figure 3.
3.1.2. The SB case (strongly interacting case)
For the SB case, the overall interaction dynamics among two successive downwash and runup flows are similar because the breakers are all plunging breakers, except that the wave breaking locations are different. Relatively speaking, SB12 breaks at the nearest onshore location, while SB23 breaks at the farthest location from the shoreline. However, SB34, SB45 and SB56 break roughly at the same location. This feature correlates well with the runup patterns reported in Reference Wu, Higuera and LiuWu21: the first wave produces the highest runup height, while the second wave generates the lowest runup height (refer to Reference Wu, Higuera and LiuWu21, figure 17e).
The HSPIV raw images and corresponding velocity fields during SB45 are shown in figure 4. Unlike the NB case, the non-stationary hydraulic jump is not observed before the subcritical flow near the left (offshore) side of the FOV reverses its direction (figure 4, 
$t/T=-0.057$). The second solitary wave arrives before the downwash flow of the first wave reaches the condition to develop an identifiable non-stationary hydraulic jump. This feature is further discussed in § 3.3 in terms of Froude number. An overturning jet is formed as the next uprush flow arrives (figure 4, 
$t/T = 0.000$). The clockwise overturning jet eventually merges with the counterclockwise bed-generated eddies. Finally, the clockwise overturning jet strikes on the beach slope, pushing the large broken bore up the slope (figure 4, 
$t/T=0.114$). The clockwise motion from the breaking wave dominates once the broken bore is fully formed. The results suggest that the time duration of downwash is not long enough for the development of the non-stationary hydraulic jump and the broken bore as observed in the NB case (see further discussions in § 3.3).
3.2. Ensemble-averaged velocity
In this section, the ensemble-averaged velocity fields are presented. For both NB and SB cases, although the overall processes are very similar among all five interactions, the detailed flow structures are different for the first two uprush and downwash interactions (i.e. NB12/NB23 and SB12/SB23) because the interactions have not reached the quasi-steady state until NB34 and SB34, respectively. After the third wave, the detailed flow structures become similar in both NB and SB cases. In addition to the velocity profiles to be shown in this section, the videos of the velocity vector fields for NB12, NB45, SB12 and SB45 are provided in the supplementary material. In the velocity profile plots (figures 5 and 6), the 
$z$ coordinate is normalised by the local water depth 
$h_L$ and the velocities 
$\langle U \rangle$ and 
$\langle W \rangle$ are normalised by the offshore constant water depth 
$h$ and the initial incident wave height 
$H_o$.
Figure 5. Ensemble-averaged velocity profiles 
$\langle U \rangle$ from (a) to (e) and 
$\langle W \rangle$ from (f) to (j) for the interactions NB12 (
$\triangleleft \triangleleft \triangleleft$), NB23 (
$\triangleright \triangleright \triangleright$), NB34 (
$\circ \circ \circ$), NB45 (
$+++$) and NB56 (
$\square \square \square$) at 
$x=-39.21$ cm: (a,f) 
$t/T=-0.114$, (b,g) 
$t/T=0.000$, (c,h) 
$t/T=0.057$, (d,i) 
$t/T=0.227$ and (e,j) 
$t/T=0.341$.
Figure 6. Ensemble-averaged velocity profiles 
$\langle U \rangle$ from (a) to (e) and 
$\langle W \rangle$ from (f) to (j) for the interactions SB12 (
$\triangleleft \triangleleft \triangleleft$), SB23 (
$\triangleright \triangleright \triangleright$), SB34 (
$\circ \circ \circ$), SB45 (
$+++$) and SB56 (
$\square \square \square$) at 
$x=-40.08$ cm. The time instants were the same as those shown in figure 5.
For the NB case, 10 snapshots of the ensemble-averaged velocity profiles, 
$\langle U \rangle$ and 
$\langle W \rangle$, in the water column at 
$x=-39.21$ cm (
$x/\lambda =-0.075$) are plotted in figure 5. Overall, the velocity profiles for NB34, NB45 and NB56 are similar in general. Relatively large deviations in both 
$\langle U \rangle$ and 
$\langle W \rangle$ appear during the downwash flow phase from 
$t/T=-0.114$ to 0.000 (figure 5a,b,f,g) when the influence of the non-stationary hydraulic jump and broken bore is significant. The downwash velocity profiles indicate that flows rotate in the counter-clockwise direction, which continue to evolve as the swash flow reverses its direction from downwash to uprush at 
$t/T=0.057$ (figure 5c,h). The flow becomes unidirectional, as the runup height is reaching its maximum (figure 5d,j). The flow decelerates and the velocity profiles becomes almost uniform in the water column. However, the bottom boundary layer flow structure is also clearly shown (figure 5d,e). Overall 
$\langle W \rangle$ is smaller than, but of the same order of magnitude of 
$\langle U \rangle$ during the downwash phases when the hydraulic jump and broken bore occur (figure 5a,b,f,g).
For the SB case, the 
$\langle U \rangle$ and 
$\langle W \rangle$ profiles in the water column at 
$x=-40.08$ cm (
$x/\lambda =-0.108$) for all five interactions are plotted in figure 6. The velocity profiles for SB12 and SB23 different from the other three, i.e. SB34, SB45 and SB56, because the first two solitary waves break at different locations. SB12 breaks at the farthest onshore (refer to figure S2 in the supplementary material), whereas SB23 breaks at the farthest offshore location. Comparing with the NB case, the flow depths in the NB case are thinner (roughly between 4 and 6 cm). The velocity profiles are mostly uniform during the downwash among SB34, SB45 and SB56 (figure 6a,b). Relatively larger deviation occurs when the incident wave breaks (figure 6c). During the early uprush phases for the SB case 
$\langle W \rangle$ is relatively large (figure 6h). The velocity profiles of SB34, SB45 and SB56 collapse into one when the flow becomes unidirectional during the uprush phase (figure 6d,e,i,j). The velocity profile of SB23 is close to those of SB34, SB45 and SB56 during the later uprush phases because the velocity profiles are almost uniform in the water column when the uprush flow decelerates.
At the quasi-steady state for the NB case, the maximum downwash 
$\langle U \rangle$ is 30 % of the phase speed 
$c$ at 
$t/T=-0.114$ and the maximum uprush 
$\langle U \rangle$ is 20 % of 
$c$ at 
$t/T=0.057$. For the SB case, the maximum downwash 
$\langle U \rangle$ is 40 % of 
$c$ at 
$t/T=-0.114$ and the maximum uprush 
$\langle U \rangle$ is 35 % of 
$c$ at 
$t/T=0.057$.
Both the wave-height-to-water-depth ratio 
$H_o/h$ and the separation time between two successive solitary waves play an important role in the development of the broken bore. In the present experiments the normalised separation time for both NB and SB cases is the same, 
$\tau /T \approx 1$, in which 
$\tau$ is the separation time and 
$T$ is the effective wave period. However, 
$H_o/h$ is 0.11 for the NB case and 0.22 for the SB case. Because of the relatively smaller 
$H_o/h$ for the NB case, the flow depth of the downwash flow diminishes before the next wave arrives, causing a short temporary ‘stationary’ broken bore in the NB case. In the SB case, the second solitary wave arrives not only before the downwash flow depth completely disappears, but also when the downwash flow is still strong. All the SB interactions (SB12 to SB56) show a similar pattern of flow field of a plunging breaker with the broken bore generated by wave breaking, i.e. the clockwise rotation inside the front roller fluid domain and counter-clockwise eddies generated along the bed.
3.3. Temporal and spatial distribution of Froude number
The temporal and spatial variation of the Froude number 
$Fr$, defined in (3.1), and the local depth, 
$h_L$, of NB45 and SB45 are plotted in figures 7 and 8, respectively. The sign of 
$Fr$ indicates the flow direction, positive being uprush and negative downwash. Therefore, the flow is supercritical for 
$|Fr|>1$, critical for 
$|Fr|=1$ and subcritical for 
$|Fr|<1$. The magnitude of 
$|Fr|$ for both NB and SB cases is within the range from 0 to 2.5. The flow changes spatially from supercritical to subcritical conditions at many instants (e.g. 
$t/T$ from 
$-$0.2 to 
$-$0.1 for NB45 shown in figure 7 and 
$t/T$ from 
$-$0.1 to 
$-$0.05 for SB45 shown in figure 8).
Figure 7. Temporal and spatial evolution of the non-dimensionalised local water depth (
$h_L/H_o$) and the Froude number (
$Fr$) times the sign of the depth-averaged velocity for NB45: (a) 
$h_L/H_o$, the vertical dashed lines A, B, C, D, E, and F indicate 
$t/T=-0.200$, 
$-$0.171, 
$-$0.140, 
$-$0.080, 
$-$0.038 and 0.038, respectively; (b) 
$Fr$ with the direction of the flow; the contour lines of 
$Fr=-1$, 
$-$1.5, 0 and 1 are traced as the black lines; (c) temporal evolution of 
$\mathcal {R}$ in (3.2). Interval from A to B: non-breaking jump; interval from B to C: weak jump; interval from C to D: weakening jump; interval from D to E: transition to broken bore; interval from E to F: momentarily ‘stationary’ broken bore.
Figure 8. Temporal and spatial evolution of the non-dimensionalised local water depth (
$h_L/H_o$) and the Froude number (
$Fr$) times the sign of the depth-averaged velocity for SB45: (a) 
$h_L/H_o$; (b) 
$Fr$ with the direction of the flow. The contour lines of 
$Fr=-1$, 
$-$1.5, 
$0$ and 
$1$ are traced as the black lines.
It is well known that for one-dimensional (1-D) rectangular open channel flows, a hydraulic jump occurs when the flow is changing spatially from supercritical flow to subcritical flow, while the upstream and downstream velocities are moving in the same direction. A stationary hydraulic jump is developed only if upstream Froude number 
$Fr_1$ and the downstream Froude number 
$Fr_2$ satisfy the following relation (Streeter, Wylie & Bedford Reference Streeter, Wylie and Bedford1998):
\begin{equation} \mathcal{R}=\frac{Fr_2}{Fr_1}\left (\frac{\sqrt{1+8Fr_1^2}-1}{2}\right)^{3/2} =1. \end{equation}
Note that when 
$\mathcal {R}$ is less than one, the hydraulic jump is non-stationary. Moreover, hydraulic jumps can have different free surface features depending on the upstream Froude number 
$Fr_1$ (Das Reference Das2008; Chanson Reference Chanson2009).
The magnitudes of 
$\mathcal {R}$ in (3.2) for NB45 are estimated so as to gain insights of swash flows. Since a hydraulic jump must occur from a depth less than the critical depth (the depth at 
$|Fr|=1$) to a depth larger than the critical depth, only the results during the downwash flow are shown in figure 7(c). The vertical dashed lines indicate the intervals of the downwash processes (also shown in figure 7a). The lines are labeled as A to F, being located at 
$t/T=-0.200$, 
$-$0.171, 
$-$0.140, 
$-$0.080, 
$-$0.038 and 0.038, respectively. The upstream and downstream Froude numbers 
$Fr_1$ and 
$Fr_2$ are estimated at 
$x=-25.67$ and 
$-$51.66 cm (
$x/\lambda = -0.049$ and 
$x/\lambda = -0.098$), respectively. As the upstream Froude number 
$Fr_1$ increases, 
$\mathcal {R}$ increases. Figure 7(c) shows that 
$\mathcal {R}$ is ranging from 0 to 0.5.
The early downwash flow is subcritical (
$|Fr| < 1$) across the 
$x$ direction within the measurement area at 
$t/T<-0.3$ (figure 7b). The variation of water depth 
$h_L$ is small across the 
$x$ direction at 
$t/T<-0.3$ (figure 7a). The free surface profile at the same time appears to be a straight line (refers to the supplementary raw image video for NB45 at 
$t/T<-0.3$). The CFP, defined as the point at 
$|Fr|=1$, begins at the far right of the measurement area (
$x=-17.2$ cm, 
$x/\lambda = -0.033$) at 
$t/T=-0.3$ (figure 7b). The 
$x$ location for the CFP has been identified with a vertical red line in the corresponding instantaneous velocity fields as shown in figure 3. At this instant, 
$\mathcal {R}$ is small (0.3 to 0.4) as shown in figure 7(c) and the velocity field is uniform in the water column (e.g. figure 3, 
$t/T=-0.284$).
As the downwash flow accelerates in time, the raw image at 
$t/T=-0.200$ shows that the water level has a smooth rise of 3 cm within 20 cm from 
$x\approx -45$ to 
$-$25 cm (the dashed line A in figure 7a). A non-stationary hydraulic jump begins at 
$t/T=-0.200$ (
$\mathcal {R}=0.42$ at A in figure 7c). Almost half of the downwash flow in the FOV is supercritical (
$|Fr| \approx 1.5$) on the onshore side and half is subcritical on the offshore side (figure 7b, 
$t/T=-0.200$). The smooth rise of water level continues until 
$t/T=-0.171$ (B) at which 
$\mathcal {R}\approx 0.5$. This initial hydraulic jump from 
$t/T=-0.200$ to 
$-$0.171 (the interval from A to B in figure 7c) is classified as a non-breaking jump with the upstream 
$|Fr|\approx 1.5$. The magnitude of 
$\mathcal {R}$ from A to B increases from 0.42 to 0.5 steadily as the non-breaking jump intensifies due to the increasing strength of the upstream flow as phase advances from 
$t/T=-0.20$ to 
$-$0.171.
The free surface shows initial breaking at 
$t/T=-0.171$ (figure 3). The non-breaking jump evolves into a breaking jump when the upstream 
$|Fr|$ increases from 
$\approx$1.7 to 2.5 between 
$t/T=-0.171$ and 
$t/T=-0.140$ (the interval from B to C in figure 7c). This breaking hydraulic jump is classified as a weak hydraulic jump as the 
$Fr$ is between 1.7 and 2.5 (Das Reference Das2008). The magnitude of 
$\mathcal {R}$ is at the maximum in the interval B to C (
$\mathcal {R}\approx 0.5$). The change in the local water depth 
$h_L$ between the upstream and downstream becomes more abrupt around 
$x=-38$ cm (
$|Fr|\approx 1$) from B to C (figure 7a,b from B to C).
The CFP moves to farther offshore direction until 
$t/T=-0.114$ at 
$x=-37$ cm (identified with the vertical line in the instantaneous velocity field 
$t/T=-0.114$ in figure 3). Then, the CFP remains roughly at the same 
$x$ location until 
$t/T=-0.08$. This is the farthest offshore 
$x$ location at which the CFP reaches.
In the interval between C and D (
$t/T=-0.140$ to 
$-$0.08), the 
$Fr$ estimation becomes noisy because the upstream flow depth is very shallow, roughly 1 cm (figure 7a, C to D). 
$\mathcal {R}$ decreases gradually as the upstream flow is weakening. The weakening jump evolves into a broken bore from 
$t/T=-0.08$ to 
$-$0.038 when the downwash flow depth diminishes (figure 7c, D to E). 
$\mathcal {R}$ values are noisy and close to zero because the downstream velocity is approaching to zero before flow reversal (
$Fr_2\approx 0$).
As shown in figure 7(b), the water depth is zero in the white area from 
$x=-17$ to 
$-$35 cm (
$x/\lambda =-0.032$ to 
$-$0.067) between 
$t/T=-0.038$ to 0.110. In the area where the water depth does not vanish, the 
$|Fr|$ values are almost zero from 
$t/T=-0.038$ to 0.038 (E to F) because the broken bore barely moves forward or backward when the flow reverses from downwash to uprush. This is the short duration, about 0.27 s, when the momentarily ‘stationary’ broken bore is identified (e.g. figure 3, 
$t/T=0.000$). In the same period, 
$t/T=-0.038$ to 0.038, the shoreline remains roughly at 
$x\approx -35$ cm (
$x/\lambda = -0.067$) as shown in figure 7(a, E to F).
The primary difference between the NB and SB cases is that the non-stationary hydraulic jump can not be identified during downwash in the SB case. As shown in figure 8(a), the variation of water depth across 
$x$ is small during the downrush phases and no sudden change of flow depth is observed for SB45. Therefore, the 
$\mathcal {R}$ values are not presented. Figure 8(a,b) shows that the flow is subcritical (
$|Fr|<1$) across the FOV when 
$t/T\approx {<}-0.151$ (e.g. figure 4, 
$t/T=-0.284$ and the supplementary raw image video SB45 from 
$t/T=-0.293$ to 
$-$0.151). The CFP moves from the onshore location (
$x=-17.2$ cm, 
$x/\lambda = -0.046$) rapidly to the offshore direction from 
$t/T=-0.151$ to 
$-$0.100, and the CFP is slowed down due to the influence from the incoming wave (figure 4, 
$t/T=-0.057$). The farthest location at which the CFP reaches is 
$x=-46$ cm (
$x/\lambda = -0.124$) roughly at 
$t/T=-0.05$ (figure 8b). The CFP stays at this location for a short time before the point of 
$|Fr|=0$ pushes the CFP back to the onshore direction rapidly from 
$t/T=0$ to 0.1 as shown in the line of 
$|Fr|=0$ in figure 8(b) (refers to figure 4, 
$t/T=0.000$ to 0.114). Even though the downwash flow on the right-hand side reaches supercritical condition (
$|Fr|=1.5$), the distance of the subcritical flow is not long enough for the hydraulic jump to develop (from 
$x=-47$ to 
$-$52 cm, 
$x/\lambda = -0.127$ to 
$-$0.140 at 
$t/T=-0.02$ shown in figure 8b). The leading flow during uprush reaches supercritical (
$|Fr|>1$) (
$t/T\approx 0.08$ to 0.12 and 
$x=-23$ to 
$-$32 cm, 
$x/\lambda = -0.062$ to 
$-$0.086). Two points of 
$|Fr|=0$ occur within the wave cycle when the flow reverses direction at 
${\approx }t/T=0$ as the downwash reverses to uprush and at 
$t/T=0.4$ as the uprush reverses to the downwash. These two points move up the beach slope when the phase advances as shown by the two 
$|Fr|=0$ contour lines in figure 8(b), indicating that the flow in the offshore region always changes direction first.
3.4. Turbulent velocity intensity
In this section, the turbulence characteristics of the fourth wave–wave interactions (NB45/SB45) are presented in figure 9 as they are representative of the quasi-steady state. In the supplementary material, the profiles of the turbulent velocity intensities 
$\sqrt {\langle {u^2}\rangle }$ and 
$\sqrt {\langle {w^2}\rangle }$ for all five interactions are provided at a fixed location 
$x=-39.21$ cm (
$x/\lambda =0.075$) for the NB case and at 
$x=-40.08$ cm (
$x/\lambda =0.108$) for the SB case (roughly in the mid-section of the FOV) in order to show the comparison among the five interactions. The locations (
$x=-39.21$ cm, 
$x/\lambda =0.075$) and (
$x=-40.08$ cm, 
$x/\lambda =0.108$) are chosen because the variation of the turbulent intensity is relatively large due to breaking.
Figure 9. The turbulent velocity intensities 
$\sqrt {\langle u^2 \rangle }$ and 
$\sqrt {\langle w^2 \rangle }$ profiles in the 
$z$ direction, (a) and (b) are for NB45, (c) and (d) are for SB45; (
$\circ \circ \circ$) 
$t/T=-0.114$; (
$+ + +$) 
$t/T=-0.057$; (
$\square \square \square$) 
$t/T=-0.028$; (
$\times \times \times$) 
$t/T=0.000$; (
$\triangleleft \triangleleft \triangleleft$) 
$t/T=0.028$; (
$\triangleright \triangleright \triangleright$) 
$t/T=0.057$; (
$\triangle \triangle \triangle$) 
$t/T=0.114$; (
$\triangledown \triangledown \triangledown$) 
$t/T=0.227$; (
$\ast \ast \ast$) 
$t/T=0.284$; (
$\diamond \diamond \diamond$) 
$t/T=0.341$.
The turbulent velocity intensities 
$\sqrt {\langle {u^2}\rangle }$ and 
$\sqrt {\langle {w^2}\rangle }$ for NB45 (figure 9a,b) are high during the downwash phases when the non-stationary hydraulic jump and the broken bore occur. The maximum 
$\sqrt {\langle {u^2}\rangle }$ is 15 % of the phase speed 
$c$. (figure 9a,b from 
$t/T=-0.114$ to 0.000). During the uprush phase, 
$\sqrt {\langle {u^2}\rangle }$ and 
$\sqrt {\langle {w^2}\rangle }$ decrease rapidly once the broken bore has been pushed to the onshore direction from 
$t/T=0.114$ to 0.341. The turbulent velocity intensity 
$\sqrt {\langle u^2 \rangle }$ is approximately 
$1\,\%$ of 
$c$ at 
$t/T=0.114$ to 0.341 because the velocity profiles of the uprush flow is quite uniform in the water column at 
$t/T=0.114$ to 0.341 (figure 5d,e). The turbulence production is insignificant as the velocity gradient is small and the turbulent kinetic energy cannot maintain itself.
The turbulent velocity intensities 
$\sqrt {\langle {u^2}\rangle }$ and 
$\sqrt {\langle {w^2}\rangle }$ for SB45 reveals the typical evolution of turbulence for plunging breaking wave as discussed in Sou et al. (Reference Sou, Cowen and Liu2010), i.e. the main turbulent kinetic energy generation happens during wave breaking with 
$\sqrt {\langle {u^2}\rangle } \approx 15\,\%$ of 
$c$ (figure 9c,d, 
$t/T=0.028$ to 0.057). The turbulent kinetic energy decays rapidly during the decelerating uprush phases with 
$\sqrt {\langle {u^2}\rangle } \approx 3\,\%$ of 
$c$ in the water column (figure 9c,d, 
$t/T=0.114$ to 0.341). During the downwash phase, the turbulent kinetic energy is generated from the bed but the 
$\sqrt {\langle {u^2}\rangle }$ in the water column decays to roughly 
$1\,\%$ of 
$c$. As mentioned, the downwash duration is too short for the non-stationary hydraulic jump to develop for SB45, such that the main source of turbulence during the downwash is the bed-generated turbulence for the SB case.
3.5. Spatial spectral analysis
To accurately calculate the full spatial spectra, the velocity vector field must cover a large spatial extent with a fine enough spatial resolution of the order of the Kolmogorov length scale. In the present experiments, the spatial resolution 
$\Delta x$ of the velocity vector field is 0.222 cm with 159 data points in 35.35 cm across the 
$x$ direction for NB45 and 0.207 cm with 159 data points in 33.05 cm across the 
$x$ direction for SB45. The largest length scale resolved in the spatial spectra depends on the data length in the 
$x$ direction and the smallest scale depends on the spatial resolution. Following directly from the Nyquist theorem, the largest wavenumber 
$\kappa _1$ that can be resolved equals to 
$2{\rm \pi} /(2\Delta x)={\rm \pi} /\Delta x$. Thus, the smallest scale resolved in the spatial spectra is twice of the spatial resolution 
$\Delta x$, which is 0.444 cm (
$\kappa _1=1.42\times 10^3$ rad m
$^{-1}$) and 0.420 cm (
$\kappa _1=1.50\times 10^3$ rad m
$^{-1}$) for the NB45 and SB45 cases, respectively. The largest length scale resolved in each spatial spectrum may change with the elevation 
$z$ as the number of data available across the 
$x$ direction in calculating the spectra varies with the local water depth. The full measurement length in the 
$x$ direction gives the largest length scale 
$=$ 35.35 cm (
$\kappa _1=17.77$ rad m
$^{-1}$) and 33.05 cm (
$\kappa _1=19.01$ rad m
$^{-1}$) for the NB45 and SB45 cases, respectively.
The present experimental data are capable of providing accurate spectra in the large-scale region because they accurately resolve the large-scale energy injection from the large eddies of the order of the integral length scales to be discussed in the following § 4.1. However, the spatial resolution of the PIV data is not fine enough to resolve the region of the Kolmogorov length scale in the spectra. The Kolmogorov length scale may be estimated using the dissipation rate of the turbulent kinetic energy as 
$\eta =(\nu ^3/\epsilon )^{1/4}$ , where 
$\nu$ is the kinematic viscosity and 
$\epsilon$ is the dissipation rate. The dissipation rate can be roughly estimated using the scaling method, 
$\epsilon =U_c^3/l_c$ (Tennekes & Lumley Reference Tennekes and Lumley1972). The characteristic velocity 
$U_c$ is of the order of the maximum uprush turbulent kinetic energy intensity in the 
$x$ direction 
$\sqrt {\langle u^2 \rangle }=20$ cm s
$^{-1}$ because the energy dissipated at the rate 
$\epsilon$ should be proportional to the rate of energy supply from the mean flow, and the maximum 
$\sqrt {\langle u^2 \rangle }$ gives the most conservative (smallest) estimate of the Kolmogorov length scale in the flow. The characteristic length scale 
$l_c$ is of the order of the integral length scale taken as 2 cm. Thus, the resulting dissipation rate 
$\epsilon$ is 0.4 m
$^2$ s
$^{-3}$ at the maximum uprush phase and the Kolmogorov length scale 
$\eta$ is 0.0398 mm, which is 56 times smaller than the spatial resolution of the present experimental data. Therefore, the spatial spectra are not capable of resolving the dissipation range and the present discussion will be focused on the large-scale region of the energy cascades.
Special treatment is made for determining the time-dependent ensemble-averaged spatial spectra in order to obtain a representative mean quantity as the integral length scale discussed in § 4.1 is estimated from the spatial spectrum. First, the ensemble-averaged turbulent kinetic energy spatial spectra 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$ are determined from the 10 1-D turbulent velocity spectra of the turbulent velocity field at a given elevation 
$z$ and phase 
$t/T$. Then, 10 individual ensemble-averaged spectra from the neighbouring 10 time steps are time-averaged, given the condition that the number of the data in the 
$\kappa _1$ direction (the wavenumber in the 
$x$ direction) for calculating the spectrum of the neighbouring time step remains the same. Since the HSPIV sampling rate is 1000 Hz, the total averaging in time is no more than 0.01 s. As a result, the final averaged spectrum is the average of 100 spatial spectra or less (10 runs 
$\times$ 10 time steps). 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$ are normalised so that their integration equals to the variances 
$\langle u^2 \rangle$ and 
$\langle w^2 \rangle$ of the turbulent velocity data for the spectra. Despite of averaging multiple time steps, the variances obtained from the single-phase-averaged spectra and the multiple-phase-averaged spectra are very close. The ensemble-averaged spatial spectra 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$ presented are the multiple-phase-averaged spectra. The spectra are non-dimensionalised with the maximum variance among the spectra 
$\langle u^2 \rangle _m$ (
$146$ cm
$^2$ s
$^{-2}$ for NB45 and 
$169$ cm
$^2$ s
$^{-2}$ for SB45) and the integral length scale 
$\langle L_{xx} \rangle$ (to be discussed in § 4.1). Each turbulent velocity spatial spectrum plot presented is labelled with the phase 
$t/T$ of the turbulent velocity field, the elevation 
$z$ of the turbulent kinetic velocity components in the 
$\kappa _1$ direction for the calculating the spectra of 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$, and the variances of the turbulent velocity components 
$\langle u^2 \rangle$ and 
$\langle w^2 \rangle$ accordingly. The vertical line 
$\langle L_{xx} \rangle /z$ represents the scale corresponding to the elevation 
$z$ for the individual spectrum.
3.5.1. Spatial spectra for NB45
The most interesting phases are at the non-stationary hydraulic jump (e.g. 
$t/T=-0.114$ in figure 3), the initiation of broken bore (e.g. 
$t/T=0.000$) and the uprush broken bore (e.g. 
$t/T=0.057$) for NB45. The spatial spectra and the corresponding raw image and ensemble-averaged velocity fields of the downwash phase at 
$t/T=-0.114$, the reversal phase at 
$t/T=0.000$ and the uprush phase at 
$t/T=0.057$ are shown in figures 10, 11 and 12, respectively.
Figure 10. The raw image, ensemble-averaged velocity field and spatial spectra of the downwash at 
$t/T=-0.114$ for NB45: (
$\circ$) 
${\langle S_{uu} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (
$+$) 
${\langle S_{ww} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (——-) 
$-$5/3 slope; (
$-----$) 
$-$1 slope; (
$-\cdot -\cdot -$) 
$-$3 slope; (vertical 
$-\cdot -\cdot -$) 
$\langle L_{xx} \rangle$/
$z$. Here 
$\langle u^2 \rangle _m = 146$ cm
$^2$ s
$^{-2}$ is the maximum 
$\langle u^2 \rangle$ among the spectra for NB45.
Figure 11. The raw image, ensemble-averaged velocity field and spatial spectra at 
$t/T=0.000$ for NB45: (
$\circ$) 
${\langle S_{uu} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (
$+$) 
${\langle S_{ww} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (——-) 
$-$5/3 slope; (
$-----$) 
$-$1 slope; (
$-\cdot -\cdot -$) 
$-$3 slope; (vertical 
$-\cdot -\cdot -$) 
$\langle L_{xx} \rangle$/
$z$. Here 
$\langle u^2 \rangle _m = 146$ cm
$^2$ s
$^{-2}$ is the maximum 
$\langle u^2 \rangle$ among the spectra for NB45.
Figure 12. The raw image, ensemble-averaged velocity field and spatial spectra of the uprush at 
$t/T=0.057$ for NB45: (
$\circ$) 
${\langle S_{uu} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (
$+$) 
${\langle S_{ww} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (——-) 
$-$5/3 slope; (
$-----$) 
$-$1 slope; (
$-\cdot -\cdot -$) 
$-$3 slope; (vertical 
$-\cdot -\cdot -$) 
$\langle L_{xx} \rangle$/
$z$. Here 
$\langle u^2 \rangle _m = 146$ cm
$^2$ s
$^{-2}$ is the maximum 
$\langle u^2 \rangle$ among the spectra for NB45.
As shown in figure 10 at 
$t/T=-0.114$, the energy level in the spatial spectra is at the highest near the bed at 
$z=0.444$ cm and decreases as the elevation 
$z$ increases. The energy injection around 
$\kappa _1 \approx 100$ to 500 rad m
$^{-1}$ (
$\kappa _1 \langle L_{xx} \rangle \approx 0.3$ to 2) near the bed is due to the eddies generated by the bed. At 
$z=4.22$ cm, both 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$ show the typical 
$-$5/3 slope without the energy injection as those near the bed at 
$z=0.444$ cm, as the turbulence is mostly homogeneous upstream of the flow field. These results indicate that the bed-generated eddies have not been evolved into the upper water body at the earlier stage of the downwash flow when the free surface breaks. The spectra energy level of 
$\langle S_{uu} \rangle$ is consistently higher than that of 
$\langle S_{ww} \rangle$ at 
$z=0.444$ and 4.22 cm but the turbulence is more isotropic in the mid depth at 
$z=1.33$ and 2.00 cm, indicating that the boundary effect from the bed and free surface are important during the downwash. The vertical turbulent motion is constrained near the free surface and bottom boundary during the downwash.
As shown in figure 11, at 
$t/T=0.000$, the large-scale energy level of the spectra at both the elevation 
$z=0.444$ and 1.33 cm increases compared with that of 
$t/T=-0.114$ (figure 10), yielding a 
$-$5/3 slope in the larger-scale region. The lower energy level in the vertical spectrum 
$\langle S_{ww} \rangle$ compared with the spectrum 
$\langle S_{uu} \rangle$ suggests that the boundary effect is still important near the bed at 
$z=0.444$ cm and near the free surface at 
$z=4.22$ cm even though the broken bore is momentarily ‘stationary’ at 
$t/T=0.000$. The energy is closer to isotropic with a 
$-$5/3 slope between the wavenumber 
$\kappa _1=2\times 10^2$ to 
$4\times 10^2$ rad m
$^{-1}$ (
$\kappa _1 \langle L_{xx} \rangle \approx 0.6$–1.3) when it is away from the bed at 
$z=1.33$ and 2 cm compared with the earlier downwash flow at 
$t/T=-0.114$, indicating that the broken bore generates more isotropic turbulence with the largest length scale in the flow for NB45.
As indicated in figure 12 during the uprush flow at 
$t/T=0.057$, the spectrum 
$\langle S_{uu} \rangle$ shows a 
$-$1 slope in the region of the smaller 
$\kappa _1$ as the larger-scale energy decreases when the broken bore is pushed forward by the next incoming wave. This is the moment when the ensemble-averaged velocity (
$\langle U \rangle \approx 60$ cm s
$^{-1}$) is at the maximum near the bed at 
$z=1$ cm. The results suggest that the turbulent kinetic energy level cannot maintain itself as the next incoming wave pushes the broken bore forward. The next incoming wave does not increase the turbulent kinetic energy production as the flow is relatively uniform behind the broken bore front and the region behind of the broken bore is free of air bubbles (i.e. the mean velocity gradient 
$\partial \langle U \rangle / \partial z\approx 0$). The majority of the turbulent kinetic energy comes from the non-stationary hydraulic jump and broken bore during the downwash phase for the NB45 case.
3.5.2. Spatial spectra for SB45
As shown in figure 13 during the downwash flow at 
$t/T=-0.114$ for the SB45 case, the energy levels in both 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$ are relatively low compared with the downwash flow at 
$t/T=-0.114$ for NB45. The 
$-$1 slope at the elevation 
$z=0.831$ and 1.87 cm is observed. Being consistent with the results discussed in Sou et al. (Reference Sou, Cowen and Liu2010), the turbulent kinetic energy cannot maintain itself during the later uprush flow and earlier downwash flow. Despite of the lower energy level in both 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$, the flow is still turbulent with the flow Reynolds number (
$Re=\langle U \rangle h_L/\nu$) of the order of 
$10^4$.
Figure 13. The raw image, ensemble-averaged velocity field and spatial spectra of the downwash at 
$t/T=-0.114$ for SB45: (
$\circ$) 
${\langle S_{uu} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (
$+$) 
${\langle S_{ww} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (——-) 
$-$5/3 slope; (
$-----$) 
$-$1 slope; (
$-\cdot -\cdot -$) 
$-$3 slope; (vertical 
$-\cdot -\cdot -$) 
$\langle L_{xx} \rangle /z$. Here 
$\langle u^2 \rangle _m = 169$ cm
$^2$ s
$^{-2}$ is the maximum 
$\langle u^2 \rangle$ among the spectra for SB45.
As displayed in figure 14 when the next incoming wave arrives at 
$t/T=0.000$, the spectrum 
$\langle S_{uu} \rangle$ shows a typical 
$-$5/3 slope as the larger-scale energy increases due to the large motion from the incoming breaking wave. The energy level of 
$\langle S_{ww} \rangle$ is almost at the same level as 
$\langle S_{uu} \rangle$ in a small range at 
$\kappa _1\approx 200$ rad m
$^{-1}$ (
$\kappa _1 \langle L_{xx} \rangle \approx 0.6$), indicating that there is energy injection from the vertical component of the length scale of 
$\kappa _1\approx 200$ rad m
$^{-1}$ (
$\kappa _1 \langle L_{xx} \rangle \approx 0.6$), possibly due to the large motion near the bed as the flow reverses direction from the downwash to uprush. This large eddy structure is visible in the ensemble-averaged velocity field at 
$x=-52$ to 
$-$45 cm and 
$z=2$ cm at 
$t/T=0.000$ in figure 14. The motion of the seeding particles at 
$x=-30$ to 
$-$25 cm and 
$z=2$ cm in the raw image at 
$t/T=0.057$ for SB12 shows a similar eddy structure (refer to movie 5 in the supplementary material). This particle motion can be thought of the sediment suspension in the flow to the extent proportional to the length scale of the order of 3 cm (
$\kappa _1\approx 200$ rad m
$^{-1}$, 
$\kappa _1 \langle L_{xx} \rangle \approx 0.6$).
Figure 14. The raw image, ensemble-averaged velocity field and spatial spectra at 
$t/T=0.000$ for SB45: (
$\circ$) 
${\langle S_{uu} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (
$+$) 
${\langle S_{ww} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (——-) 
$-$5/3 slope; (
$-----$) 
$-$1 slope; (
$-\cdot -\cdot -$) 
$-$3 slope; (vertical 
$-\cdot -\cdot -$) 
$\langle L_{xx} \rangle /z$. Here 
$\langle u^2 \rangle _m = 169$ cm
$^2$ s
$^{-2}$ is the maximum 
$\langle u^2 \rangle$ among the spectra for SB45.
During the uprush flow at 
$t/T=0.114$, the energy level of both 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$ increases as shown in figure 15. The turbulence becomes close to isotropic at 
$z=5.19$ cm due to wave surging into the water in front, generating a broken bore with high energy level covering the widest range (
$\kappa _1=70$ to 700 rad m
$^{-1}$, 
$\kappa _1 \langle L_{xx} \rangle \approx 0.2$–2.2 at 
$z=5.19$ cm) among all the spectra during the uprush–downwash interaction. Although the high energy level range at 
$z=0.831$ cm is not as wide as that at 
$z=5.19$ cm, the energy level of 
$\langle S_{uu} \rangle$ around 
$\kappa _1=200$ to 300 rad m
$^{-1}$ (
$\kappa _1 \langle L_{xx} \rangle \approx 0.6$–1) at 
$z=0.831$ cm is actually higher than that at 
$z=5.19$ cm. This is mainly due to the shear layer shown in the ensemble-averaged velocity field (
$x=-30$ to 
$-$20 cm at 
$z \approx 1$ cm) with the length scale of the order of 2 cm. The shear layer is generated when the uprush flow runs above the downwash flow from 
$t/T=0.05$ to 
$t/T=0.100$ (refer to movie 7 in the supplementary material). Due to the bubbles from the plunging breaker, the shear layer cannot be examined closely.
Figure 15. The raw image, ensemble-averaged velocity field and spatial spectra of the uprush at 
$t/T=0.114$ for SB45: (
$\circ$) 
${\langle S_{uu} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (
$+$) 
${\langle S_{ww} \rangle }/{({\langle u^2 \rangle }_{m} \langle L_{xx} \rangle )}$; (——-) 
$-$5/3 slope; (
$-----$) 
$-$1 slope; (
$-\cdot -\cdot -$) 
$-$3 slope; (vertical 
$-\cdot -\cdot -$) 
$\langle L_{xx} \rangle /z$. Here 
$\langle u^2 \rangle _m = 169$ cm
$^2$ s
$^{-2}$ is the maximum 
$\langle u^2 \rangle$ among the spectra for SB45.
4. Further discussions on turbulence characteristics
4.1. Integral length scales
In this section, the time-dependent turbulence structures in terms of the integral length scale are examined using the spatial turbulent velocity fields. The integral length scale represents the correlation distance of turbulent velocity in space and the largest eddy size in the flow. In this study, the integral length scale is calculated directly using the spatial spectrum of the spatial turbulent velocity vector field (Tennekes & Lumley Reference Tennekes and Lumley1972; Pope Reference Pope2000). Following Sou et al. (Reference Sou, Cowen and Liu2010), 
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle$ are evaluated from the 1-D spectra as
\begin{equation} \langle L_{xx} \rangle=\frac{\rm \pi}{2}\frac{\langle S_{uu}(\kappa_1=0)\rangle}{\langle u^2 \rangle},\quad \langle L_{zz} \rangle=\frac{\rm \pi}{2}\frac{\langle S_{ww}(\kappa_1=0)\rangle}{\langle w^2 \rangle}, \end{equation}
where 
$\kappa _1$ is the wavenumber in the 
$x$ direction. Since the turbulent energy level at the wavenumber 
$\kappa _1 = 0$ cannot be resolved with the data, the energy level at the smallest 
$\kappa _1$ resolved in the spectrum is used. The size of the HSPIV measurement area and the spatial resolution of the velocity vector field determine the range of the wavenumber 
$\kappa _1$ (or the length scale) that can be resolved in the spatial spectra.
An alternative approach for estimating the integral length scale using the spatial turbulent velocity field is to integrate the turbulent velocity autocorrelation function (Variano & Cowen Reference Variano and Cowen2008). Both the spectral method and the autocorrelation function method require a sufficiently large measurement length to ensure the inclusion of largest motions. The autocorrelation method also requires a large region of homogeneous turbulence. Otherwise, the resulting integral length scale will be biased towards the local turbulent length scale at the center point (Variano & Cowen Reference Variano and Cowen2008). Since the turbulence is inhomogeneous across the 
$x$ direction in our case, the integral length scale estimated directly from the spatial spectrum as shown in (4.1).
At each elevation 
$z$, 
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle$ are evaluated from 
$\langle S_{uu} \rangle$ and 
$\langle S_{ww} \rangle$, respectively, using the data points in the 
$x$ direction. A valid integral length scale estimate must satisfy three conditions. First, the minimum number of data points required for calculating the spectrum is chosen to be 
$2^4=16$. This gives a minimum data record length of 15 times the spatial resolution and 7 times the smallest length scale resolved in the spatial spectrum (refer to § 3.5). The minimum number of data points is chosen to be a power of two because of the fast Fourier transform (FFT) in the spectrum calculation. Second, the data length must be three times larger than the resulting integral length scale so that the largest scale motion is included (Variano & Cowen Reference Variano and Cowen2008). Lastly, the resulting integral length scale must be larger than the smallest possible length scale resolved from the spectra, i.e. twice the spatial resolution (
$2\times \Delta x=0.444$ cm for NB45 and 0.420 cm for SB45 as discussed in § 3.5). Most of the integral length scale estimate satisfies these conditions except a few estimates near the free surface. Only the valid integral length scales are presented.
The maximum integral length scales 
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle$, and their dimensionless values 
$\langle L_{xx} \rangle / h_l$ and 
$\langle L_{zz} \rangle / h_l$ at different phases for both NB45 and SB45 cases are listed in table 1. Here 
$h_l$ is the local water depth on the offshore (left) side of the measurement area. The maximum 
$\langle L_{xx} \rangle$ is mostly larger than the maximum 
$\langle L_{zz} \rangle$ except at 
$t/T=-0.114$ for NB45 (
$\langle L_{xx} \rangle \approx \langle L_{zz} \rangle$) and at 
$t/T=-0.057$ for SB45 (
$\langle L_{xx} \rangle < \langle L_{zz} \rangle$). At 
$t/T=-0.114$ for NB45, the maximum 
$\langle L_{zz} \rangle$ (
$\approx \langle L_{xx} \rangle$) occurs at 
$z=5$ cm at which the free surface breaks due to the hydraulic jump (figure 16a–e, 
$t/T=-0.114$). At 
$t/T=-0.057$ for SB45, 
$\langle L_{zz} \rangle$ is larger than 
$\langle L_{xx} \rangle$ because of the influence from the incoming wave as shown in figure 4 at 
$t/T=-0.057$. The vertical motion when the flow is about to reverse yields larger 
$\langle L_{zz} \rangle$. Interestingly, the maximum 
$\langle L_{xx} \rangle$ for NB45 and SB45 are very close at 
$t/T=0.000$ (
$\langle L_{xx} \rangle =2.21$ and 2.33 cm).
Figure 16. Spatial variation of the non-dimensionalised integral length scales at various times for NB45 (a–e) and SB45 (f–j): (
$\circ \circ \circ$) 
$\langle L_{xx} \rangle / h_l$; (
$+ + +$) 
$\langle L_{zz} \rangle / h_l$.
Table 1. Temporal variation of the maximum integral length scale and the non-dimensionlised maximum integral length scale 
$\langle L_{xx} \rangle$, 
$\langle L_{xx} \rangle /h_l$, 
$\langle L_{zz} \rangle$ and 
$\langle L_{zz} \rangle /h_l$ for NB45 and SB45.
An eddy of a certain wavenumber 
$\kappa _1$ generated in the flow can be thought of some additional energy around the length scales corresponding to 
$\kappa _1$. This eddy energy injection contributes a spike energy level to the spectrum at 
$\kappa _1$. Tennekes & Lumley (Reference Tennekes and Lumley1972) suggested that the eddy energy injection with the length scale of the wavenumber 
$\kappa _1$ contributes to a broad spike in the spectrum between 0.62
$\kappa _1$ to 1.62
$\kappa _1$ due to the interaction among the eddies of the same order of magnitude. This explains the consistent energy injection spike at the wavenumber around 
$\kappa _1=100$ to 400 rad m
$^{-1}$ (
$\kappa _1 \langle L_{xx} \rangle \approx 0.3$–1.3) in most of the spatial spectra because the characteristic large eddy size of the order of 
$\langle L_{xx} \rangle$ is around 1–2 cm (e.g. figure 10, 
$t/T=-0.114$). When the large-scale energy from the eddies becomes homogeneous due to the momentarily ‘stationary’ broken bore, a 
$-$5/3 slope is observed in the spectra at the smaller 
$\kappa _1$ region (e.g. figure 11, 
$t/T=0.000$).
In addition, the variation of the integral length scales in the 
$z$ direction at the same 
$t/T$ for both NB45 and SB45 is shown in figure 16. Since the local water depth within the FOV varies, the depth at the left side of the FOV, 
$h_l$, is chosen to non-dimensionalise the integral length scales in figure 16. Both 
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle$ are of the order of 1–2 cm. Overall, 
$\langle L_{xx} \rangle$ is roughly twice 
$\langle L_{zz} \rangle$ in the lower elevation 
$z$. Here 
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle$ are about the same in the higher elevation 
$z$ for most of the time instants for NB45 and SB45 (figure 16).
For NB45, 
$\langle L_{xx} \rangle$ is larger near the bed during the downwash phase at 
$t/T=-0.114$ due to the bed-generated eddies (figure 16a–e, 
$t/T=-0.114$). As the phase advances to 
$t/T=-0.057$, 
$\langle L_{xx} \rangle$ increases because the bed-generated eddies start growing and diffusing upwards into the water column (figure 16a–e, 
$t/T=-0.057$). At 
$t/T=0.000$, 
$\langle L_{xx} \rangle$ roughly maintains the same magnitude as that at 
$t/T=-0.057$ throughout the water column because of the energy supply from the momentarily ‘stationary’ broken bore. Then, 
$\langle L_{xx} \rangle$ decreases roughly from 2 to 1 cm near the bed when the next incoming wave arrives at 
$t/T=0.057$. As the broken bore is pushed onshore at 
$t/T=0.114$, 
$\langle L_{xx} \rangle$ is roughly 2 cm in the water column from 
$z=1$ to 5 cm.
For SB45, 
$\langle L_{xx} \rangle$ is large (2.77 cm) at 
$t/T=-0.114$ due to the relatively strong downwash flow (figure 6f–j, 
$t/T=-0.114$). It should be noted that 
$\langle L_{xx} \rangle =2.77$ cm is of the order of the local water depth at 
$t/T=-0.114$ as shown in figure 4. This implies that the eddy size can be as large as the local water depth when the flow becomes shallow during the downwash phase (e.g. the counter-clockwise eddies revealed by tracer particles during the downwash phase in the abstract figure). As the next breaking wave arrives at 
$t/T=0.000$, 
$\langle L_{xx} \rangle$ is larger near the bed as the flow is still downwash near the bed (figure 4, 
$t/T=0.000$). Between 
$z=2$ and 4 cm 
$\langle L_{zz} \rangle$ is larger than 
$\langle L_{xx} \rangle$ due to the strong vertical motion when the flow reverses (figures 16 and 4g–l, 
$t/T=0.000$). During the uprush phases 
$t/T=0.114$ and 0.057, 
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle$ in the upper water column increase due to the plunging broken bore.
In comparison, between NB45 and SB45, figure 16 indicates that the integral length scales 
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle$ are ranging from 
$0.05h_l$ to 
$0.32h_l$ with the maximum of 
$0.32h_l$ at the broken bore phase in the wave cycle for NB45. For SB45, the integral length scale is at the largest during the downwash phases, which is comparable to the local water depth (
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle \approx h_l$ at 
$t/T=-0.114$ and 
$-$0.057). The results of integral length scale from the previous studies of spilling and plunging waves show that the integral length scale is ranging from 0.05 to 0.4 of the local water depth from outer surf zone to the swash zone (Ting & Kirby Reference Ting and Kirby1996; Petti & Longo Reference Petti and Longo2001; Govender et al. Reference Govender, Mocke and Alport2004; Sou et al. Reference Sou, Cowen and Liu2010). The integral length scales reported in the previous studies were not larger than 0.4 of the local water depth. The present results illustrate that the integral length scale is comparable to the local water depth when the water depth is shallow during the downwash.
The turbulence Reynolds number can be defined as 
$Re_T=\langle u^2 \rangle ^{1/2} \langle L_{xx} \rangle / \nu$ (Pope Reference Pope2000), which yields 
$2.42 \times 10^3$ and 
$3.35\times 10^3$ for NB45 and SB45, respectively. The Taylor-microscale Reynolds number can also be estimated as 
$Re_{\lambda }\approx (20 Re_T)^{1/2}$ (Pope Reference Pope2000), resulting in 220 and 259 for NB45 and SB45, respectively. The variances of the turbulent velocities of the uprush flow from both cases are employed in these Reynolds numbers calculations. The turbulence Reynolds number is typically an order of magnitude smaller than the flow Reynolds number (
$Re=U h_L/\nu$), e.g. assuming 
$\sqrt {\langle u^2 \rangle }/U \approx 0.2$ and 
$L_{xx}/h_L\approx 0.5$. In general, 
$Re_{\lambda }=150$ is considered reasonably high in many laboratory flows (Pope Reference Pope2000). In addition, turbulence can be anisotropic even in flows with a reasonably high Reynolds number (e.g. George & Hussein Reference George and Hussein1991). The turbulence Reynolds numbers for both cases are of the same order of magnitude, although NB45 is weakly interacting and SB45 is strongly interacting. This implies that the energy levels of turbulence are of the same order of magnitude even though the sources of the turbulence generated are different for the NB and SB cases.
4.2. Eddy viscosity
A better understanding of the eddy viscosity 
$\nu _T$ for the breaking wave processes is essential for improving the numerical modelling of the swash flows. Following Schlichting & Gersten (Reference Schlichting and Gersten2000), the ensemble-averaged eddy viscosity 
$\langle \nu _T \rangle$ can be expressed as
\begin{equation} \langle \nu_{T}\rangle= c_p\langle L_{xx} \rangle\sqrt{\langle k \rangle}, \end{equation}
where 
$c_p$ is an empirical coefficient based on the Prandtl's mixing length model and 
$c_p \approx 0.55$ for turbulent boundary layer (Schlichting & Gersten Reference Schlichting and Gersten2000). Here 
$\langle k \rangle$ is defined as 
$\langle k \rangle = 1.33 (\langle u^2 \rangle + \langle w^2 \rangle )/2$. The factor 1.33 accounts for the turbulent kinetic energy in the 
$y$ direction similar to the plane wake (Svendsen Reference Svendsen1987). The turbulent length scale 
$\langle L_{xx} \rangle$ is used instead of 
$\langle L_{zz} \rangle$ as 
$\langle w^2 \rangle$ may be more affected by the bubbles in the aerated regions. Figure 17 shows the temporal and spatial variation of the ensemble-averaged eddy viscosity 
$\langle \nu _T \rangle$ calculated according to (4.2), corresponding to the same time instants shown in figure 16 for the integral length scales 
$\langle L_{xx} \rangle$ and 
$\langle L_{zz} \rangle$. The 
$x$ and 
$z$ coordinates are non-dimensionalised by the maximum 
$\langle L_{xx} \rangle$ at the corresponding 
$t/T$ listed in table 1. The eddy viscosity 
$\langle \nu _T \rangle$ is non-dimensionlised by the offshore constant water depth 
$h$ and the phase speed 
$c$. The eddy viscosity describes the spatial variation of the turbulent kinetic energy transfer at different time instants, ranging from 
$0.0005hc$ to 
$0.006hc$ (of the order of 
$10^{-4}$ to 
$10^{-3}$ m
$^2$ s
$^{-1}$) spatially and temporally during the uprush–downwash interaction. The evolution of the eddy viscosity shows that turbulent kinetic energy transfer occurs when the eddies are generated continuously from the bed, indicating the forming of turbulent boundary layer. The order of magnitude of the maximum eddy viscosity is consistent with that estimated by Svendsen (Reference Svendsen1987), i.e. 
$\nu _T=0.01 h \sqrt {gh} = 37$ cm
$^2$ s
$^{-1}$ (order of magnitude of 
$10^{-3}$ m
$^2$ s
$^{-1}$) and the same order of magnitude of 
$10^{-3}$ m
$^2$ s
$^{-1}$ found in Petti & Longo (Reference Petti and Longo2001).
Figure 17. Temporal and spatial variations of the non-dimensionalised eddy viscosity 
${\langle \nu _T \rangle }/{hc}$ for NB45 (a–e) and SB45 (f–j). The 
$x$ and 
$z$ coordinates are non-dimensionlised by the maximum 
$\langle L_{xx} \rangle$ at the corresponding phase 
$t/T$ listed in table 1.
For NB45, 
$\langle \nu _T \rangle$ is at the maximum (
$\approx$25 cm
$^2$ s
$^{-1}$, 
${\approx } 0.006hc$) at 
$z=1$ cm during the downwash and reversal phases (figure 17a–e, 
$t/T=-0.057$ and 0.000). The spatial distribution of 
$\langle \nu _T \rangle$ reveals the similar pattern of the eddies near the bed as shown in the HSPIV raw image in the abstract figure and the HSPIV raw image video (movie 1) in the supplementary material. The high-intensity pattern dimension of 
$\langle \nu _T \rangle$ is of the same order of magnitude of the integral length scale as shown at 
$z/\langle L_{xx} \rangle \approx 0$ to 1 when 
$t/T=-0.057$ and 0.000 (figure 17a–e). The magnitude is of the order of 5–10 cm
$^2$ s
$^{-1}$ during the later uprush phases and the earlier downwash phases (figure 17a–e, 
$t/T=-0.114$, 0.057 and 0.114). The range along the depth is narrower compared to the later downwash phase and the reversal phase (figure 17a–e, 
$t/T=-0.057$ and 0.000).
For SB45, the maximum 
$\langle \nu _T \rangle$ (
$\approx$25 cm
$^2$ s
$^{-1}$ around 
$z/\langle L_{xx} \rangle =0$ to 1) occurs during the uprush phases (figure 17f–j, 
$t/T=0.000$ and 0.114). The spatial distribution also shows the shear layer at 
$t/T=0.057$ when the uprush flow runs over the downwash flow (figure 17f–j, 
$t/T=0.057$). As the uprush phase advances to 
$t/T=0.114$, 
$\langle L_{xx} \rangle \sqrt {\langle k \rangle }$ is at the maximum near the bed when the plunging broken bore moves up the slope (figure 17f–j, 
$t/T=0.114$). The high-intensity pattern dimension corresponds to the same order of magnitude of the integral length scale 
$\langle L_{xx} \rangle$. During the downwash at 
$t/T=-0.114$, the integral length scale 
$\langle L_{xx} \rangle$ is approximately the same as the local water depth 
$h_L$ as 
$z/\langle L_{xx} \rangle \approx 1$ at the free surface (figure 17f–j at 
$t/T=-0.114$). Here 
$\langle \nu _T \rangle \approx 2$ cm
$^2$ s
$^{-1}$ even though the integral length scale (
$\langle L_{xx} \rangle =2.77$ cm) is large (figure 17f–j, 
$t/T=-0.114$). It is because the turbulent kinetic energy is small during the downwash phase as the turbulent energy production is small. The results suggest that the eddy viscosity between the uprush phase and the downwash phase can be different by an order of magnitude.
Large eddy viscosity value (
${\approx }0.006hc$) is associated with the bed-generated eddies when the hydraulic jump and broken bore occur with relatively large integral length scale and the turbulent kinetic energy for NB45. For SB45, high eddy viscosity occurs in the shear layer and the broken uprush bore front when turbulent kinetic energy production is high.
4.3. The 
$-$1 spectral law revisited
Nikora (Reference Nikora1999) suggested that the bed-generated eddies with the characteristic length scale of the of the elevation 
$z$ modify the slope of the inertial subrange in the 1-D spatial spectrum for the homogeneous turbulence from 
$-$5/3 to 
$-$1 slope as
\begin{equation} \langle S_{uu}(\kappa_1) \rangle \sim \epsilon^{2/3} \kappa_1^{{-}5/3} \sim u_*^2 \kappa_1^{{-}1}. \end{equation}
The increase in the dissipation rate 
$\epsilon$ because of the energy injection from the bed-generated eddies with the characteristic scale 
$z$ is of the order 
$u_*^3/z$, in which 
$u_*$ is the friction velocity. Studying swash flows generated by periodic waves, Sou et al. (Reference Sou, Cowen and Liu2010) demonstrated that the 
$-$1 slope region exists between the length scales of the order of the elevation 
$z$ and the boundary layer thickness in the energy spectra near the bed during the downwash phases at four locations covering from the inner surf zone to the swash zone.
In this study, the energy injection from the large-scale eddies of the order of the integral length scale 
$\langle L_{xx} \rangle$ is observed in the energy spectra during the uprush–downwash interaction (e.g. the spectra in the 
$x$ direction 
$\langle S_{uu}(\kappa _1) \rangle$ during the downwash phases for NB45 and SB45 shown in figures 10 and 13). The large-scale energy injection of the order of the integral length scale (
$\kappa _1 \langle L_{xx} \rangle \approx 1$) due to the bed-generated eddies at 
$t/T=-0.114$ is observed in the spectra for both NB45 and SB45 (figures 10 and 13). These results suggest that the 
$-$1 spectral law can be extended to the large-scale region. The dissipation rate 
$\epsilon$ is of the order of 
$\sqrt {\langle u^2 \rangle }^3/ \langle L_{xx} \rangle$ for the large-scale eddies (Tennekes & Lumley Reference Tennekes and Lumley1972). Thus, the superposition of the energy injection due to the length scales of the order of 
$\langle L_{xx} \rangle$ modifies the slope from 
$-$5/3 to 
$-$1 in the energy cascade as
\begin{equation} \langle S_{uu}(\kappa_1) \rangle \sim \epsilon^{2/3} \kappa_1^{{-}5/3} \sim \langle u^2 \rangle \langle L_{xx} \rangle \sim \langle u^2 \rangle \kappa_1^{{-}1}. \end{equation}
The turbulent energy injection is not limited within the boundary layer thickness during the uprush–downwash interaction. When the full spatial spectrum is measured accurately, it is possible to have multiple 
$-$1 slope regions in the energy cascade if the flow has multiple energy injections from the eddies with different dominant length scales of the spatial spectrum.
5. Concluding remarks
In this paper, the characteristics of the interaction between the uprush and downwash of consecutive solitary waves on a 1/10 beach slope have been investigated using the HSPIV measurements. The velocity fields for six consecutive waves were resolved for two wave conditions: weakly interacting case and strongly interacting case. The time separation between two consecutive waves is roughly the characteristic wave period of solitary wave for both cases. The first wave of the weakly interacting case does not break, while the first wave of the strongly interacting case is a surging breaker. Non-stationary hydraulic jump occurs during the interaction between the downwash and uprush flows for the weakly interacting case. The separation time between waves of the strongly interacting case is too short for the non-stationary hydraulic jump to happen before the next plunging breaking wave arrives.
Swash flow reaches the quasi-steady state after the third solitary wave for the two wave conditions. Flow similarity is observed during the third, fourth and fifth uprush–downwash interactions. The turbulence structures of the fourth interactions are examined using the spatial spectral analysis. The integral length scale is quantified directly from the spatial turbulent kinetic energy spectra without assuming homogeneous isotropic turbulence and invoking the Taylor frozen hypothesis because the PIV measurement provides large enough spatial velocity vector field to include the large-scale motion. The temporal and spatial evolution of the integral length scales show the bed-generated eddies evolve and diffuse from near-bed into upper water body during the downwash and reversal phases for the weakly interacting case. The largest integral length scale occurs during the flow reversal from uprush to downwash due to the momentarily broken bore for the weakly interacting case, while the largest integral length scale occurs during the downwash flow due to the strong velocity in a shallow flow depth for the strongly interacting case. The results of the eddy viscosity, estimated from the integral length scale and turbulent kinetic energy, suggest that the temporal and spatial variation of eddy viscosity within the wave cycle depends on the breaker types. The maximum eddy viscosity occurs during the later downwash and reversal phases for the weakly interacting case due to the non-stationary hydraulic jump and broken bore. For the strongly interacting case, the maximum eddy viscosity occurs during the uprush phase due to the plunging breaker. The maximum value of the eddy viscosity for both cases are of the same order of magnitude. The eddy viscosity between the uprush and the downwash can be different by an order of magnitude for both cases.
The spatial turbulent kinetic energy spectra for both cases show that the turbulence is close to isotropic in the water column. The strongly interacting case behaves similarly to the unsteady boundary layer flow during the downwash phases. The 
$-$1 slope is apparent due to the large-scale energy injection from the larger-scale eddies of the order of the integral length scale 
$\langle L_{xx} \rangle$. An eddy of a certain wavenumber 
$\kappa _1$ (or length scale) generated in the flow adds energy at the length scales corresponding to 
$\kappa _1$. Thus, the eddy energy injection at a certain wavenumber 
$\kappa _1$ (or length scale) contributes a spike energy level to the spectrum at 
$\kappa _1$, modifying the slope of the energy cascade as shown in the results of our spatial spectra. The typical 
$-$5/3 slope is modified to a 
$-$1 slope not only for the smaller scale region, but also for the larger-scale region corresponding to the integral length scales using the scaling method. Without the energy injection at a certain length scale (the integral length scale in our case), the spectrum energy cascades at 
$-$5/3 slope in the inertial subrange, i.e. the energy cascade for homogeneous turbulence in which the energy cascade in the inertial subrange of the spectrum depends on the dissipation rate. The results suggest that the 
$-$1 spectral law is not limited to the bed-generated eddies with the characteristic length scale of the boundary layer during the downwash flow. It can be applied to the large-scale energy injection of the eddies with the length scale of the order of the integral length scale due to wave breaking within the uprush–downwash interaction.
The separation time between the two consecutive waves and the initial wave-height-to-water-depth ratio are the key factors to determine the characteristics of the uprush–downwash interaction. Although the results in this work show that either non-stationary hydraulic jump and broken bore occur for the weakly interacting case or onshore wave broken bore occurs for the strongly interacting case, it is possible to generate a condition that the non-stationary hydraulic jump is followed by the next wave breaking within a single uprush–downwash interaction. By increasing the separation time between the two consecutive waves and the ratio of the initial wave height to water depth, such condition is observed in our wave flume. The longer separation time ensures that there is enough downwash duration for the forming of the non-stationary hydraulic jump before the next incident wave approaches. Larger ratio of the initial wave height to water depth is needed to ensure that the next incident wave height is large enough to break over the relatively deep subcritical flow depth of the hydraulic jump.
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2023.484.
Acknowledgements
Insightful discussions with Ignacio Barranco and Pablo Higuera during the early stage of the work are also appreciated. Helpful assistance from Yufei Wang with wave generation is greatly appreciated.
Funding
P.L.-F.L. would like to acknowledge the support from the National University of Singapore, Cornell University, and the Ministry of Education in Singapore through a research grant (MOE2018-T2-2-040). This research was also supported in part by the Yushan Program, Ministry of Education in Taiwan. Y.-T.W. would like to offer thanks for support from the Ministry of Science and Technology, Taiwan (MOST 108-2218-E-006-053-MY3) and, in part, by the Higher Education Sprout Project, Ministry of Education to the Headquarters of University Advancement at National Cheng Kung University.
Declaration of interests
The authors report no conflict of interest.
Author contributions
I.M.S.: conceptualisation; conducting experiments; data analysis and interpretation; writing. Y.-T.W.: conceptualisation; conducting experiments; data analysis and interpretation; writing. P.L.-F.L.: conceptualisation; writing – review and editing; supervision; project administration.
Appendix A. Free surface identification
The free surface profiles were identified by manually locating the 
$x$ and 
$z$ locations in pixels from the raw PIV images at least every 200 image frames and as many as every 5 frames when the free surface profile changed rapidly during wave breaking. Following Barranco & Liu (Reference Barranco and Liu2022), the full set of a single realisation of the free surface profile were interpolated in space and time with the selected points. As shown in figure 18, the green line is the estimated free surface profile and the red line is the beach bed. The bed was identified by manually locating the 
$x$ and 
$z$ points based on the reflection of the seeding particles near the glass bed. Only the region of tracer particles was included as the flow field. The blurry region in the raw image was out of the measurement plane and it was excluded. The free surface profiles during the downwash phases were relatively easy to be identified because they were typically smooth (figure 18a). The free surface profiles were not perfect when the wave broke especially during plunging breaking for the SB case (figure 18c,d). As one 
$x$ value in pixel could not have two values of 
$z$ for the interpolation, the free surface profiles had to include a portion of the blurry region within the plunger when the wave broke (figure 18c,d).
Figure 18. Examples of free surface and bed profiles during downwash and uprush: (a) downwash for the NB case; (b) uprush for the NB case; (c) early uprush for the SB case; (d) uprush for the SB case where the green line is the identified free surface and the red line is the bed.
Appendix B. Turbulence statistics convergence test
In order to verify whether the number of repeating realisations was sufficient to determine turbulence statistics, a sensitivity analysis was carried out. The effects of number of realisations on the root mean square of velocity fluctuations 
$\sqrt {\langle {u^2}\rangle }$ and 
$\sqrt {\langle {w^2}\rangle }$ were investigated. The normalised deviation 
$\Delta D_N$ of 
$N$ number of repeating realisations is defined as
\begin{equation} \Delta D_N = \frac{\sum_{k_xl_z}|D_N-D_{N+1}|}{\sum_{k_xl_z}D_{N+1}}, \end{equation}
where 
$D_N$ is the estimate of the turbulent velocity intensities 
$\sqrt {\langle {u^2}\rangle }$ and 
$\sqrt {\langle {w^2}\rangle }$. The sum is taken over all the 
$k_x \times l_z$ vectors of the measurement area. The total repeating realisations are 11 and 12 for NB45 and SB45, respectively. The convergence test was calculated on data during the most energetic time period from 
$t/T=-0.114$ to 
$t/T=0.341$. Figure 19 indicates that the values of 
$\Delta D_N$ for both 
$\sqrt {\langle {u^2}\rangle }$ and 
$\sqrt {\langle {w^2}\rangle }$ approach to 0.05 when 
$N$ increases to 10 and 11 for NB45 and SB45, respectively.
Figure 19. Convergence test of the turbulence statistics: (a) NB45; (b) SB45; (
$\circ \circ \circ$) 
$D=\sqrt {\langle {u^2}\rangle }$; (
$+ + +$) 
$D=\sqrt {\langle {w^2}\rangle }$.
Melville et al. (Reference Melville, Veron and White2002) evaluated the normalised error 
$\Delta E_N$ of the turbulence kinetic energy at each 
$N$ by defining 
$\Delta E_N$ equal to the difference between the 
$N$ number of repeating realisations and their total number of realisations (
$N=24$) as
\begin{equation} \Delta E_N = \frac{\sum_{k_xl_z}|E_N-E_{N=24}|}{\sum_{k_xl_z}E_{N=24}}. \end{equation}
Their normalised error 
$\Delta E_N$ converged to 0.16 with the total realisations of 
$N=24$. Although our normalised deviation 
$\Delta D_N$ is defined differently, the last value of 
$\Delta D_N$ can be compared with that of the normalised error 
$\Delta E_N$ in Melville et al. (Reference Melville, Veron and White2002) as the last normalised errors from both approaches are the difference between the last two 
$N$ numbers of realisations, i.e. 
$N$ and 
$N+1$. In addition, the sensitivity analysis indicates that our solitary waves are highly repeatable as 
$\Delta D_N$ reaches close to 0.16 when 
$N=3$ compared with 
$\Delta E_N \approx 0.16$ when 
$N=24$ in Melville et al. (Reference Melville, Veron and White2002).
Appendix C. Uncertainty analysis
The 
$95\,\%$ confidence intervals of the ensemble-averaged velocities and turbulent velocity intensities (
$\langle U \rangle$, 
$\langle U \rangle$, 
$\sqrt {\langle {u^2} \rangle }$ and 
$\sqrt {\langle {w^2} \rangle }$) are found by using the bootstrap uncertainty analysis (Efron & Tibshirani Reference Efron and Tibshirani1993). A single phase 
$t/T$ has only 10 realisations and the number of samples is insufficient for bootstrap analysis. Following the same approach used for calculating the spatial spectrum in § 3.5, the data samples from the 10 neighbouring phases of a single phase 
$t/T$ are used to obtain the measurement uncertainty (i.e. the data within 0.01 seconds). The 
$95\,\%$ confidence intervals for the downwash, reversal and uprush phases for NB45 and SB45 are presented in tables 2 and 3, respectively. The confidence intervals are calculated from the data points in the region with the higher turbulent kinetic energy level. For NB45, the 
$95\,\%$ confidence levels for all 
$t/T$ are based on the data in the region from 
$x=-45.88$ to 
$-$39.21 cm and 
$y=1.33$ to 
$3.11$ cm. For SB45, the region of data points is from 
$x=-40.07$ to 
$-31.75$ cm and 
$y=0.83$ to 
$1.26$ cm for 
$t/T=-0.114$ and 
$0.000$ at which the flow depth is shallow. For the strong uprush plunging broken bore at 
$t/T=0.114$, the region is from 
$x=-40.07$ to 
$-$31.75 cm and 
$y=0.83$ to 
$7.06$ cm. The typical uncertainty in the ensemble-averaged velocity and turbulent velocity intensity is approximately 
$1\,\%$. As expected, the uncertainty for 
$\langle W \rangle$ is typically higher than that of 
$\langle U \rangle$ since the dynamic range of pixels in the PIV processing for 
$\langle W \rangle$ is smaller compared with that for 
$\langle U \rangle$. For both NB45 and SB45, the maximum uncertainty for 
$\langle W \rangle$ is approximately 
$10\,\%$ when breaking occurs at 
$t/T=-0.114$ and 
$t/T=0.000$ for NB45 (table 2) and at 
$t/T=0.000$ and 
$t/T=0.114$ for SB45 (table 3). This is expected as the larger confidence intervals of the measured quantities are correlated to higher uncertainty due to the breaking phases of the non-stationary hydraulic for NB45 and the plunging breaking wave for SB45.
Table 2. The 95 % confidence intervals of 
$\langle U \rangle$, 
$\langle W \rangle$, 
$\sqrt {\langle u^2 \rangle }$ and 
$\sqrt {\langle w^2 \rangle }$ during the downwash (
$t/T=-0.114$), reversal (
$t/T=0.000$) and uprush (
$t/T=0.114$) phases for NB45.
Table 3. The 95 % confidence intervals of 
$\langle U \rangle$, 
$\langle W \rangle$, 
$\sqrt {\langle u^2 \rangle }$ and 
$\sqrt {\langle w^2 \rangle }$ during the downwash (
$t/T=-0.114$), reversal (
$t/T=0.000$) and uprush (
$t/T=0.114$) phases for SB45.
Open access
