1 Experimental Studies

In last few decades, a number of important experiments have been performed and fully analyzed to reveal the intrinsic physics of breaking waves. Different measurement methods including the high-speed cameras, bubble detection approach^[5], fiber optical probes^[6-7], Laser Doppler Velocimer^[8], Particle Image Velocimetry (PIV) technologies^[9], and acoustic measurement^[10] have been adopted.

By applying a visualization technique, Bonmarin^[11] analyzed the steep wave breaking evolution. The wave profile of a plunging crest, the related splash-up, and the air-entrainment phenomenon were observed. Qiao and Duncan^[12] performed detailed analysis for the unsteady evolution of moderate spilling wave breakings. Grue and Jensen^[13] measured the velocity field and corresponding accelerations of the steep plunging waves from experiment. They found that for the horizontal and vertical accelerations over 1.1g(g indicates the acceleration of gravity) and 1.5g, the wave breaking phenomenon will be induced.

The experiment of breaking waves over an underwater hydrofoil conducted by Duncan^[14-15] was the first attempt to study the energy dissipation caused by wave breaking. Rapp and Melville^[8] studied the dispersion of the deep-water wave breakings by using Laser Doppler measurement. They measured the loss of excess momentum and energy flux of breaking waves, and found over 10% flux loss for single spilling wave breakers to as much as 25% flux loss for plunging wave breakers. Chang and Liu^[16] investigated the turbulent flow originated from wave breaking with limited water depth. By applying the ensemble average approach, the mean velocity, Reynolds stress, and strain rate were obtained. It was concluded that the generation and evolution of turbulence flows including its dissipation are of equal importance, whereas the effect of turbulence diffusion occupies a small portion. Kimmoun and Branger^[17] performed the waves breaking test in water tank by using a wedge body. They calculated the spatial derivatives to acquire adequate information for describing the flow structures.

Many existing experimental studies suggest that the air fractions in wave-turbulence interaction process can induce low-frequency-noise in the ocean. For a better understanding of the air entrainment caused by wave breaking, Lamarre and Melville^[18] exploited the bubble plumes and found that the magnitude of containing air fractions is many orders larger than previous expectation. It was concluded that the contribution of bubbles to sound propagation might be underestimated to a very large extent if these plumes of large bubbles are out of consideration. Later on, Lamarre and Melville^[19]gave detailed measurements for the evolution of the void-fraction field in bubble plumes originated from large-scale deep water breakers. Various parameters of air fractions were calculated and comparison was made with previous 2-D laboratory breaking waves. In the experiment performed by Cartmill and Su^[20], bubble radii in the bounds of 34 to 1 200 μm were measured by an acoustic sensor at different locations and deepness in a water tank, and the scatters of bubble size are well conformed with -3 power-law from 10 to 1 200 μm. Loewen et al.^[10]applied an image processing technology to measure the size scatters of large bubbles generated by breaking waves in fresh and salt water. They found that the distributions of bubble size can be described by an exponential (or power law) equation, and concluded that the gas transportation resulted from large bubbles should be taken into consideration even for the low-speed wind. By measuring the bubble dimensions for breaking waves in open water and laboratory experiment, Deane and Stokes^[5] gave a quantitative illustration of the bubble generation principles. They found that for bubbles dimensions greater than about 1 mm, the bubble density conforms a -10/3 power law correspondent to the bubble radius, in which the bubble size distribution is dominated by turbulence intensity. In terms of smaller bubbles stemming from droplet and jet impact, a -3/2 power-law can be deduced. Air entrainment of typical gravity wave breakers on a slope-beach was studied experimentally by Mori and Kakuno^[21], in which the wave surface altitude, air fractions, and bubble transportation were analyzed in detail. Experimental results revealed a linear dependence between the air portions and turbulent intensity. Besides, the bubble distributions were observed to follow a -3.4 power law against the bubble diameters. It was noted that the bubble transportation has no significant relations with the breaking positions and free surface depth. Based on the experiment of Blenkinsopp and Chaplin^[22], further measurements were performed in the water tank to study the transportation of void fractions originated from the wave breaking of fresh water, salted water, and sea water. Experimental results suggested that the influence of water type might be less significant compared with the scale effect for investigating the evolution of bubble plumes entrained due to breaking waves.

Andreas^[23] was the first to analyze the transportation of seawater spray. In the study, previous droplets creation functions were summarized, and new formula for wind speeds over 32 ms^-1 were proposed. The droplet size distribution under strong wind conditions were also studied experimentally by Veron et al.^[24] by using high speed video camera. It showed that the droplet dimension distribution follows -3 to -5 power laws when the droplets diameters vary between 196 and 5 510 mm. The creation of droplets by two kinds of plunging waves was studied systematically by Towle^[25], and results revealed that the wave breakers induced by larger wave steepness generally lead to larger numbers of droplets with greater radius.

By using high-image-density PIV, Lin and Rockwell^[9] studied the evolution of quasi-steady breaker in different stages, including the capillary wave pattern and the fully developed breaking waves. For high Froude number, surface tension plays an important role to determine the breaker characteristics. In terms of intermediate Froude numbers, the evolution of free surface exhibits capillary wave pattern. Quasi-steady breaking waves accompany with flow separation may be generated by high Froude number at the breaker zone.

Some of the representative experimental investigations for breaking waves are summarized in Table 1.

2 Numerical Studies

Experimental research can only be done for usual wave conditions due to technical difficulties. Thus, it is difficult to acquire an intensive understanding of the multi-scale wave breaking phenomenon. Due to the existence of complex hydrodynamic characteristics(e.g., complex breaking patterns, coalescence and separation of bubbles/droplets, and multiscale liquid-gas turbulence), simulating the surf zone breaking waves is still a challenge for traditional CFD. With recent progresses of the numerical approaches, it becomes feasible to resolve the two-phase flow problems numerically even with the presence of the above difficulties. Because of the recent advances in computer architecture, detailed numerical studies of the intermediate scale wave breaking problems that takes account of water droplets and air bubbles become possible. Scardovelli and Zaleski^[26] provided a well-documented review on different numerical approaches.

In early numerical studies, Lemoes^[27]; Takikawa et al.^[28]; Lin and Liu^[29-30]; Bradford^[31] adopted the Reynolds Averaged Navier-Stokes (RANS) models for the analysis of breaking waves. Although the understanding of the fundamental wave breaking phenomenon arising in the surf zone can be improved accordingly, the turbulent intensity was over predicted at the breaking regions. Hieu et al.^[32] and Zhao et al.^[33] developed multi-scale method that can improve original RANS predictions. The 2-D flow structures can be fully resolved by RANS model, while the 3-D turbulence interactions are represented by the specified eddy viscosity model. Although air entrainment was not involved, the numerical results are generally consistent with experimental data.

Large eddy simulations (LES) and direct numerical simulations (DNS) are two major branches for the high-fidelity simulation of the turbulent wave breaking. In LES, the large scales eddies are modeled directly, while sub-grid viscosity model and bubbly-flow models are utilized for representing small scale eddy. Thus, Navier-Stokes equations can be solved using the computational mesh with moderate refinement. Shi et al.^[34], Liang et al.^[35-36], and Derakhti and Kirby^[37] developed turbulent bubbly flow models for the closure of governing equations. The numerical work of Zhao and Tanimoto^[38] was the first attempt to simulate breaking waves of a two-dimensional configuration by using LES method, in which good agreement was obtained compared with experimental measurements. Watanabe and Saeki^[39] presented the first numerical results of 3-D LES of a plunging breaking wave. Similar 3-D numerical experiment of plunging breakers were conducted by Mutsuda and Yasuda^[40], in which air entrainment and bubble interactions were discussed in detail. Christensen and Deigaard^[41] studied spilling breaking wave and plunging breaking wave with different intensity. Phenomenology study with a visualization of the internal velocity field were given. 3-D vortex structures created by the wave breaker in surf zone were investigated by Watanabe et al.^[42] with the adoption of single-phase model.

In contrast, as a promising technology, DNS is becoming prevalent in solving the multi-phase flow since no parametrizations are required. It is noted that the mesh resolution should be sufficiently refined to be taken into account of the small-scale flow structures. The early-stage numerical models by Lin and Liu^[29-30] as well as Bradford^[31] failed in considering air entrainment details until the studies of Lubin et al.^[43] and Biausser et al.^[44] DNS are often adopted to study the evolution of wave breakings to provide numerical data on wave dissipation and the splashing processes. Periodic wave conditions with relatively small wavelengths are usually adopted, and representative research works can be found in Refs.[45-49]. The physical processes of wave breaking, velocity profiles, and breaking intensities for plunging breakers were investigated by Chen et al.^[45] Iafrati^[50-51] provided quantitative statistics of the spray and air fraction in breaking waves. Numerical investigation of plunging breaking wave past an underwater bump were performed by Kang et al.^[52] and Koo et al.^[53], in which the whole plunging wave breaking procedure including the jet diving and secondary jet were analyzed.

Considering the computational cost for the 3-D cases, the simulations mentioned above were mainly conducted in 2-D domain. In addition, most of the DNS studies aimed at low-Reynolds-number problems. However, since 2-D simulations cannot capture the flow structure of Langmuir revolution in wave breaking proceedings^[54], 3-D simulations are necessary to fundamentally investigate the 3-D phenomenon including bubble and spray formation. The differences between 2-D and 3-D results were compared by Lubin et al.^[43], and numerical simulation revealed that the energy dissipation after the wave breaking can be enhanced by 3-D turbulent flow, while the enhancement is less than 5% for half-wave period. To alleviate such computational burdens of 3-D simulations, Sullivan et al.^[55-56] developed an approach with intermediate complexity. Different from the complete model of the breaking waves, a stochastic term was added in momentum equation with a linearizing free-surface condition to represent the existence of wave breaking. This model has been used to study the wave breaking phenomenon considering Langmuir circulation. However, forcing terms should be defined to keep consistence with the experimental measurements. Thus the stochastic model cannot capture the air entrainment, and the density changes and buoyancy effect are not involved in current stage (Sullivan et al.^[56]).

Recently, with the rapid development of numerical schemes and computer hardware, more investigations have been conducted to 3-D DNS focusing on the vortex structures, energy dissipation, and air entrainment, represented research works were given by Christensen^[57]; Fuster et al.^[58]; Brucker et al.^[59]; Lakehal and Liovic^[60]; Zhou et al.^[61]; Lubin and Glockner^[62], etc.

In the study of Brucker et al.^[59], up to 134 million grids were used and numerical tests showed that the parallel scalability remain acceptable. More recent numerical study presented by Lubin and Glockner^[62] used over 1 billion grids to investigate the vortex filaments created at an early stage after the occurrence of breaking. In the studies of Stern et al.^[63] and Wang et al.^[64], bow wave breakings adjacent to the wedge shape obstacle and breaking wave past a underwater bump were simulated using very large grids (up to 2.2 billion), which was the first endeavor to simulate wave-breaking problems with the mesh resolution scale around micrometers. Table 2 outlines the existing simulation approaches for the fundamental analysis of breaking waves.

3 Mechanism Analysis

The breaking process of complex multiphase turbulent flows is a challenging problem for analysis. General conclusions of the unsteady and statistical characteristics of the energy dissipation, air entrainment, bubble cloud, aerated vortex filaments, and the capillary effect are summarized and presented in this section. Corresponding statistical methodologies and important conclusions are also introduced and discussed.

The typical breaking procedure can be described in following aspects. As shown in Figs. 1(a)-(c) (reproduced from Wang et al.^[64]), a jet was first formed in the wave forward direction (Fig. 1(a)), and then it impacted the water surface due to gravity (Fig. 1(b)). In the initial breaking stage, the wave dynamics exhibited 2-D flow pattern. Afterwards, the jet impacting and entrapment of air cavity brought numerous bubbles and droplets, and the flow became fully three dimensional (Fig. 1(c)). As summarized by Deane and Stokes^[5] as well as Kiger and Duncan^[66], there are two major reasons for the generation of bubble clouds, i.e., jet impacting the free surface creates small scale bubbles, while the air-bubble collapsing accounts for the formulation of larger bubbles. The influence of initial slope for waves (S) was analyzed numerically by Deike et al.^[67], and they found that the wave evolution was close to the plunging breaking with initial slope S>0.42, while a spilling process was observed for 0.35 < S < 0.42, which was also conformed with the experiment done by Rojas and Loewen^[68].

In Fig. 1, Figs. 1(d)-(e) show the bottom view of the free surface for breaking waves, in which the air entrainment and cavity break-up process can be demonstrated clearly. The wave impact that induces large air cavity is accompanied by many small-size satellite bubbles. The thin filaments of air, which is an evident phenomenon in wave breaking, are under numerical investigation by Lubin and Glockner^[62].

3.1 Energy Dissipation

The total energy E=E_k+E_g+E_s for wave breaking can be decomposed into three components (the kinetic energy E_k, the gravitational potential energy E_g, and the capillary potential energy E_s), according to the formulations given by Chen et al.^[45], Iafrati^[51], and Deike et al.^[49] as follows:

${E_k} = \frac{1}{2}\int \rho {u^2}{\rm{d}}x{\rm{d}}y $

(1)

${E_g} = \int \rho gy{\rm{d}}x{\rm{d}}y $

(2)

${E_s} = \frac{{\gamma (\mathfrak{L} - 1)}}{{{\rho _w}g{\lambda ^2}}} $

(3)

where $\mathfrak{L}$ represents the arc length. Due to the viscosity of the fluid, the total wave energy exhibits descending trend in time. In the following part, the methodologies and preliminary conclusions for breaking energy dissipation will be introduced briefly. It is noted that in most of the theoretical and numerical studies, the surface tension effect E_s is ignored.

Fig. 2(a) shows the time evolution of wave energy, in which ▽ denotes the gravitational potential energy, ◇ represents the capillary potential energy, ○ indicates the kinetic energy, and □ denotes the total energy. Fig. 2(b) shows the wave breaking factors as a function of wave steepness, in which the semi-empirical results are represented by solid line, ▼ indicates the DNS by Deike et al.^[67], and other symbols are from the experimental measurements of Banner and Peirson^[69], Drazen et al.^[70] and Grare et al.^[71] Fig. 2(c) shows the total wave energy, kinetic energy, and gravitational potential energy as a function of time.

Followed the theoretical study of Phillips^[72], the energy dissipation rate caused by breaking waves can be described by ${E_{{\rm{diss}}}} = \int {{\mathit{\boldsymbol{\epsilon}}_l}} \Lambda (c){\rm{d}}c$, where c is the phase velocity, ${{\mathit{\boldsymbol{\epsilon}}_l}}$ indicates the dissipation rate of unit width of wave crest, and Λ(c) stands for the breaking distribution, which can be measured in experiments (Gemmrich et al.^[73], Kleiss and Melville^[74], and Sutherland and Melville^[75]). The dissipation rate of Duncan^[15] and Phillips^[72] can be adopted by defining ${{\mathit{\boldsymbol{\epsilon}}_l}}$=bρ c⁵/g, where b represents the non-dimensional breaking parameter. The semi-empirical relationships of the wave breaking factors derived by Romero et al.^[76] is applied as

${E_k} = b = 0.4{(S - 0.08)^{2.5}} $

(4)

where S=0.08 is an estimation of the threshold for breaking. The accuracy of the above inertial model was validated by numerous experimental studies (Deike et al.^[49], Romero et al.^[76], Garrett et al.^[77], Melville and Pizzo^[78], and Melville and Fedorov^[79]).

Deike et al.^[49] compared the evolution of the whole wave energy and the kinetic and gravitational energy per unit length of breaking wave, as shown in Figs. 2(a)-(b). The total wave energy shows an abrupt decrease during the first wave period, which means the majority of the wave energy is dissipated. They also provided a fitted exponential function to predict energy decay which is written as E=E₀e^-ζt, where ζ represents the decay rate and E₀ stands for the primary wave energy per unit width of the wave crest. The breaking influencing factor b against the wave slope is shown in Fig. 2(b), in which the experimental data of Romero et al.^[76], Garrett et al.^[77], and Drazen et al.^[70] together with the 3-D DNS solutions (Deike et al.^[49]) were given. Deike et al.^[67] inferred that the initial setups for breaking cause the differences in dissipation rate between DNS and experimental result for S varying from 0.35 to 0.4. The kinetic, gravity, and potential energy components derived from 3-D NDS were also given by Wang et al.^[64], together with the 2-D DNS solutions by Chen et al.^[45] as shown in Fig. 2(c). Although the time evolution of energy exhibits similar trends for 2-D and 3-D simulations in general, a higher jump after jet plunges was observed for 3-D case, which is consistent with the numerical predictions by Lubin et al.^[43] and Lubin and Glockner^[62].

Deike et al.^[49] analyzed components of the breaking waves energy for various factors ${{\mathit{\boldsymbol{\epsilon}}}}$ and Bo (Bond number, Bo=Δρg/(γk²)), as shown in Figs. 3(a)-(d). Considering low surface tension (Bo=1 000) circumstance, the wave energy descends exponentially with constant decaying rate for the steepness below the breaking threshold. In the first few wave periods, the wave energy is shown to be dispersed very fast. Specifically, the gravitational energy is dissipated completely, while the kinetic fraction remains even in the immediate wave periods. They found that during the wave propagation, the surface tension potential energy could be neglected. Similar to the investigation by Iafrati^[50], the surface tension energy occupies less than 8% of the total energy dissipated.

When the surface tension effect becomes dominated (Bo=50, 100), it was concluded in Ref. [49] that if the Bond number is relatively low, more energy will be dissipated for non-breaking waves, while the total energy dissipation shows less relevance to the Bond number for the case of breaking wave.

The total energy dissipation affected by initial steepness at different Bond numbers was also provided in Ref.[49]. For low surface tension, the energy dissipation will be promoted for wave steepness larger than the threshold. In terms of strong surface tension, the existence of capillary effect will enhance the energy dissipation rate. It was observed that a larger initial steepness will further increase the total energy dissipation.

3.2 Air Entrainment and Void Fraction

Air entrainment occurs when the jet re-contacts the free surface with a large air pocket and bubbles entrained. It also occurs when the droplets drop into the water with high-speed, which was reviewed by Kiger and Duncan^[66]. However, the investigation for entrained air was not adequate in their study, so the available corresponding studies were performed by Deane and Stokes^[5], Lamarre and Melville^[18], and Blenkinsopp and Chaplin^[22].

The void fraction has been experimentally studied by repeating measurements at various locations by measuring electric conductivity^[19] and using optical fiber probes^{[6, 22]}. Deike et al.^[49] obtained void fraction variation through 3-D DNS, in which the predicted bubble flume evolution and the dynamics of void fraction are consistent with the available experimental measurements by Lamarre and Melville^[19] and Blenkinsopp and Chaplin^[22]. Both plunging case (S=0.55) and spilling breakers (S=0.43) were compared, and in the spilling case, bubble clouds generated by entrained air are accumulated adjacent to the free surface. In terms of plunging breaking waves, larger bubbles rise above the water surface into the air, while small bubbles remain in water phase for longer.

The entrained air evolving with time is presented in Fig. 4(a), which was reproduced from Deike et al.^[49] At t/T≈1, the first peak is observed since the jet impacts surface and a big pocket of air is entrapped. Then the measured air volume decreases quickly due to the collapse of the bubble and its re-contact into the free surface. After that, bubbles are created and a secondary air volume peak will be observed, which corresponds to the situation where large amount of bubbles are trapped beneath the water surface. Lamarre and Melville^[18] described the fast decay as an exponential function, V=V₀e^-Kt/T, Deike et al.^[67] deduced that the K factor varies from 2.5 to 4, and approximately, it relies on the initial wave steepness. K=3.9 and K=5 were deduced respectively by Lamarre and Melville^[19] and Blenkinsopp and Chaplin^[22] in their experiments.

Deike et al.^[67] further deduced the formulation for predicting the time-averaged air entrainment volume against the dissipation rate per unit width of wave crest or the slope, given as, V=B_l${{\mathit{\boldsymbol{\epsilon}}_l}}$L_c/ρWg · Apply the semi-empirical function for ${{\mathit{\boldsymbol{\epsilon}}_l}}$ as (Romero et al.^[76])

$\bar V = 0.4B{(S - 0.08)^{2.5}}{L_c}{c^5}/W{g^2} $

(5)

Fig. 4(b) shows the total air volume V against ${{\mathit{\boldsymbol{\epsilon}}_l}}$L_c/ρWg in which a good agreement was obtained between the theoretical and numerical scaling (Eq. (3)), and the CFD solution is also consistent with available experimental measurement by Deane and Stokes^[5], Duncan^[14], and Lamarre and Melville^[19]. Total volume of air against wave steepness S is presented in Fig. 4(c), in which a reasonable match between the theoretical model (Eq. (5)) and the DNS data is found.

3.3 Bubble and Droplet Size Distribution

The identification of the droplets and bubbles was implemented with a neighbor searching algorithm utilizing a non-recursive strategy. Herrmann^[80] proposed an efficient identification approach for classifying droplets or bubbles shared by different blocks^{[49, 67]}.

Previous experimental investigation by Deane and Stokes^[5] indicated that the size distribution of bubbles is dependent with different power law scaling. A theoretical scaling equation accounting for the turbulent fragmentation was proposed by Garrett et al.^[77], which is written as

$N(r) \propto Q{\mathit{\boldsymbol{\epsilon}}^{ - \frac{1}{3}}}{r^{ - \frac{{10}}{3}}} $

(6)

where N(r) stands for the bubble numbers in a unit mm, Q and ${{\mathit{\boldsymbol{\epsilon}}}}$ represent the air volume entrained per volume and energy dissipation respectively, and r indicates the bubble dimensions. It was found that when the bubble dimension is smaller than Hinze scale, the bubble shape can be held by surface tension effect. In this situation, the influence of turbulence is expected to be unimportant. Hinze scale was defined by Deane and Stokes^[5] as α_H=2^-8/5${{\mathit{\boldsymbol{\epsilon}}}}$^-2/5(σ We_c/ρ)^3/5, where We_c denotes the critical Weber number and σ is the surface tension coefficient. The scaling equation of small bubbles was given by Deane and Stokes^[5] and is expressed as

$N(r) \propto Q{(\sigma /\rho )^{ - \frac{3}{2}}}{v^2}{r^{ - \frac{3}{2}}} $

(7)

where v indicates the jet falling speed. Soloviev and Lukas^[81] derived a scaling equation for bubble size distribution based on the analysis of buoyancy forces dominated bubble fraction. Figs. 5(a)-(c) illustrate that bubbles created at the early stages of breaking waves are relatively small while larger bubbles are formed in the later stages of after breaking. Figs. 5(a)-(b) present the bubble size distributions for fine-grid simulations^[64], in which the bubble numbers (time-averaged) is close to the experimental measurements (Deane and Stokes^[5], Mori and Kakuno^[21], and Tavakolinejad^[82]). In the simulation work of Wang et al.^[64], the two-phase flow is not as violent as those in the experiments, so the slopes are slightly lower than the theoretical and experimental results for large bubbles.

In terms of droplet size distribution, a similar scaling equation was given by Garrett et al.^[77] and Deane and Stokes^[5], which is written as

$N(r) \propto {v^{ - 1}}{(\sigma /\rho )^{\frac{1}{2}}}{r^{ - \frac{9}{2}}} $

(8)

where v represents the jet velocity of the spray. Fig. 6 shows the theoretical and experimental results (Lhuissier and Villermaux^[83], Veron et al.^[24], and Towle^[24]), as well as the numerical predictions of Wang et al.^[64]for time-averaged droplet size distribution. It is noted that there are significant differences between the experimental studies and simulations for analyzing the spray formation mechanisms, (e.g., the plunging wave breaking^[26]), the high wind speed generated spray spume drops (Veron et al.^[24]), and bubble bursting (Lhuissier and Villermaux^[83]), however, the distribution of drop size also conforms a scaling power-law. Based on the DNS results of Deike et al.^[49], the time-averaged bubble size distribution N(r) against the wave breaking time is presented in Fig. 6, in which N(r) is found to obey a power law with $N(r) \propto {r^m}$, and m falls into the range of 3-3.5 according to the numerical results. The size distribution predicted by the above formula is consistent with many experimental measurements (e.g., Deane and Stokes^[5], Rojas and Loewen^[6], Blenkinsopp and Chaplin^[7], and Terrill et al.^[84]).

3.4 Capillary Effects

Few researchers have paid attention to the micro-mechanism inside the wave breaking processes, considering the influence of capillary effects. The surface tension effects were firstly studied by the experiments of Liu and Duncan^[85-87]. They found that the breaking kinematic energy dissipation is able to be modified through adding surfactants. Capillary effects on the wave breaking measured in experiment were discussed extensively by Duncan^[87]. The numerical investigation of capillary effects was accomplished in the studies of Liu and Duncan^[86] and Deike et al.^[49]

Deike et al.^[49] classified various small-scale wave patterns into four groups, including the gravity induced parasitic capillary waves, non-breaking waves, and spilling and plunging breaking waves. By changing the initial steepness of the wave and the Bond number, a wave-state diagram was summarized, as shown in Fig. 7(a), where PB indicates plunging breakers, SB represents spilling breakers, PCW represents parasitic capillary waves, and NB stands for non-breaking gravity waves. From Fig. 7(b), the critical steepness ${{\mathit{\boldsymbol{\epsilon}}_c}}$ for this case is around 0.32. Thus when $\mathit{\boldsymbol{\epsilon}} < {\mathit{\boldsymbol{\epsilon}}_c}$, no wave breaking occurs. However, as pointed by Deike et al.^[49], since the critical steepness predicted by numerical simulation only gives an ideal situation, ${{\mathit{\boldsymbol{\epsilon}}_c}}$ in experiment or open sea is strongly corresponding to the initial conditions, therefore the discrepancies of ${{\mathit{\boldsymbol{\epsilon}}_c}}$ between numerical and experimental predictions are unavoidable. Accurate evaluation of critical steepness is challenging even in experimental studies, and the corresponding experiments for critical steepness were performed by Romero et al.^[76], Drazen et al.^[70], and Perlin et al.^[88]

As shown in Fig. 7, when surface tension becomes dominated, two types of wave breaking will be induced, i.e., the gravity-capillary breakers and the spilling breakers, and both are scaled as $\mathit{\boldsymbol{\epsilon}} \propto \left({1 + B{o^{ - 1/3}}} \right)$, which is consistent with the theoretical study of Longuet-Higgins^[89].

3.5 Other Research Findings Related to Wave Breaking

Recent progress in computational capabilities offers the possibility to study multiscale phenomenon existing in wave breaking. In this subsection, recent development in the simulation of vortex filaments, modulational instability, and wind turbulence over breaking waves are introduced briefly.

The vortex filaments generated in the primary stage of the wave breaking are often observed in laboratories or open sea. Lubin and Glockner^[62] confirmed the existence of the unusual vorticial structures under breaking waves and investigated the mechanisms through 3-D DNS. From numerical results, the streamwise vortex filaments were detected with an upstream obliqueness (approximately 50°) under plunging breaking waves. However, Lubin and Glockner^[62]pointed out that there is no contribution on the dissipation process due to the presence of the vortex filaments. The formation of large-scale vortex structures in the air induced by the plunging wave and its influence to the interaction between the ocean and atmosphere were investigated numerically by Iafrati et al.^[90] The energy dissipated by the wave breaker both in fluid and gas was considered. It is particularly important to note that the energy dissipation in air is larger than that in water. Recently, the wind turbulence over breaking waves was analyzed by Yang et al.^[91] in detail. Effects of wave age and wave steepness were investigated through statistics of turbulent airflow over breaking waves.

4 Conclusions and Future Prospect

In this review, previous studies for wave breaking phenomenon by laboratory experiments, numerical simulations, and semi-empirical models are summarized. The methodologies for analyzing characteristic variables during wave breaking, including the void fraction, the bubble size distribution, the energy dissipation of the breaking wave, and the time evolution of the total volume of air are illustrated.

In the last few decades, numerous experimental, theoretical, and numerical methods have been proposed for analyzing the breaking waves and their effects in different scales. For experimental and open ocean investigation, extensive measurements for the bubble density, size distribution (Lamarre and Melville^[18]; Cartmill and Su^[20]; Lin and Rockwell^[9]; Deane and Stokes^[5]; Mori and Kakuno^[21]; Blenkinsopp and Chaplin^[7]) under various initial wave steepness were performed. Acoustic resonator, fiber optical probes, and high speed camera were often used for those analyses. Two experiments (Lin and Rockwell^[9]; Qiao and Duncan^[12]) were designed aiming at analyzing the surface tension effects, in which the PIV technology was adopted for visualizing the velocity field of spilling breakers. Kinetic energy dissipation and turbulence intensity during the breaking procedure were well studied by Rapp and Melville^[8], Chang and Liu^[16], and Mori and Kakuno^[21] with PIV and high speed camera. There were also quantitative statistics for droplets distribution (Veron et al.^[24]; Towle^[25]) and air entrainment (Bonmarin^[11]; Rapp and Melville^[8]) during the violent breaking. Compared with experiment, numerical simulations gave more details and deep insights into the complex two-phase flow phenomenon. Among the prevalent numerical methods mentioned in this review, it was noticed that the RANS simulations (e.g., Takikawa et al.^[28]; Lin and Liu^[29-30]; Bradford^[31]) were only conducted in early years. Recently, the LES (e.g., Watanabe et al.^[42]; Lubin et al.^[62]; Shi et al.^[34]), especially the DNS (e.g., Lakehal and Liovic^[60]; Wang et al.^[64]; Lubin and Glockner^[62]; Deike et al.^[49]) were more preferred in the fundamental study of wave breaking problems. The statistic of air entrainment (Lubin et al.^[43]) and droplet/bubble distribution (Wang et al.^[64]), as well as energy dissipation (Deike et al.^[49]) were more simple and straightforward from a volume of fluid (VOF) field than from experimental photograph. In terms of semi-empirical method, Romero et al.^[76]proposed breaking parameter model, which can be used for estimating dissipation rate. Deike et al.^[67] deduced the formulation for predicting the volume of entrained air. Garrett et al.^[77] developed bubble size scaling equation accounting for the turbulent fragmentation. Deane and Stokes^[5] further provided a scaling equation of small bubbles and droplet size distribution.

Although existing theoretical, experimental, and numerical tools are capable to provide majority characteristics, such as the mean void fraction, and droplet/bubble distribution, many fundamental physical processes of breaking wave are lack of quantitative analysis. For experimental studies, tracking and recording evolution of single droplet/bubble are challenging, especially for fully 3-D experiment. In addition to the reliable approaches for accurate statistics of air entrainment, droplets distributions need further development, and existing PIV technologies cannot provide adequate details for intensive investigations of the turbulent two-phase flows, while the promising X-ray micro-computed tomography technology^[92] is considered to be extended for this purpose. For numerical approaches, DNS provides an important alternative in resolving wave breaking mechanism in micro-scale. Due to the high requirement of the computing resources in DNS, adaptive mesh (Liu and Hu^[93-96]), which is capable to resolve the minimum bubble/droplet by localized mesh refinement, shows obvious advantage in large scale parallel computation. Besides, the fast algorithm for two-phase flow problems, which can greatly improve the computational efficiency, can be considered in future work. In terms of theoretical methods, up to now, all the semi-empirical formulae mentioned in this review were validated by model experiment and numerical simulations, while the scaling effect is rarely considered. In the future, more experimental and numerical investigations are required to bridge the gap between model scale analysis and real scale applications.