1. Introduction
Mixed-phase turbulence is a common phenomenon in nature and engineering applications. In contrast to two-phase turbulent flow without surface breaking (Brocchini & Peregrine Reference Brocchini and Peregrine2001), mixed-phase turbulence is accompanied by violent surface deformations and surface breakups, leading to complex mass and momentum transfers across the interface.
Water jet plunging into a pool is a typical process that induces mixed-phase turbulence (Kiger & Duncan Reference Kiger and Duncan2012; Delacroix et al. Reference Delacroix, Germain, Gaurier and Billard2016). Plunging jets were investigated early in chemical engineering applications, such as wastewater treatment (van der Lans, Donk & Smith Reference van der Lans, Donk and Smith1979) and the mixing and reacting of liquids and gases (McKeogh & Ervine Reference McKeogh and Ervine1981). The research in this area has focused on the mechanism of air entrainment and the characteristics of the resulting bubble flows. It was found that the impact velocity of the liquid jet plays a dominant role in the air entrainment inception conditions. It was reported that the impact velocity is a key criterion of air entrainment (Lorenceau, Quéré & Eggers Reference Lorenceau, Quéré and Eggers2004). When the critical condition was reached, a minimum amount of energy was available in the flow to do work against the surface tension and/or the potential energy of gravity to entrap air. Laboratory experiments were conducted to analyse the critical entrainment velocity for both high-viscosity liquids (Joseph et al. Reference Joseph, Nelson, Renardy and Renardy1991; Jeong & Moffatt Reference Jeong and Moffatt1992; Biń Reference Biń1993; Eggers Reference Eggers2001; Lorenceau, Restagno & Quéré Reference Lorenceau, Restagno and Quéré2003; Lorenceau et al. Reference Lorenceau, Quéré and Eggers2004) and low-viscosity liquids (Sene Reference Sene1988; Biń Reference Biń1993; Lin & Reitz Reference Lin and Reitz1998; Chirichella et al. Reference Chirichella, Gomez Ledesma, Kiger and Duncan2002; El Hammoumi, Achard & Davoust Reference El Hammoumi, Achard and Davoust2002). Clanet & Lasheras (Reference Clanet and Lasheras1997) proposed a model to predict the penetration depth of the air bubbles entrained by a water jet impacting onto a flat water surface. This model shows that the penetration depth is determined by the initial jet momentum and the bubble terminal velocity as a function of the bubble size. Comprehensive reviews of the related research can be found in Biń (Reference Biń1993) and Kiger & Duncan (Reference Kiger and Duncan2012).
The aforementioned studies focused on the vertical plunging jet, and recently, horizontal and shallow-angle plunging jets have also been investigated. Deshpande et al. (Reference Deshpande, Trujillo, Wu and Chahine2012) studied a shallow-angle plunging jet and revealed a periodic pattern of air entrainment, which does not occur when the impingement angle is steep. Later, Deshpande & Trujillo (Reference Deshpande and Trujillo2013) corroborated that the periodicity scaled linearly with the Froude number through numerical simulations. They also studied, both computationally and analytically, the underlying causes responsible for large cavity formation at shallow angles. They found a strong stagnation pressure region that deflects the entire incoming jet flow radially outwards, producing a large cavity and subsequently creating splashing events. Hsiao et al. (Reference Hsiao, Wu, Ma and Chahine2013) studied a stationary and moving horizontal jet plunging into a quiescent water pool. Their numerical and experimental results showed that the frequency of air entrainment depended on the jet diameter and relative velocity with respect to the free surface.
In addition to the plunging jet, mixed-phase turbulence with air entrainment also occurs in many other flows. Chachereau & Chanson (Reference Chachereau and Chanson2011) conducted experimental studies on the air entrainment in hydraulic jumps and observed that air entrainment takes place as the Froude number exceeds a critical value. They also found that the volume of air entrainment increases with an increasing Froude number. Garrett, Li & Farmer (Reference Garrett, Li and Farmer2000), Deane & Stokes (Reference Deane and Stokes2002), Wang, Yang & Stern (Reference Wang, Yang and Stern2016) and Deike, Melville & Popinet (Reference Deike, Melville and Popinet2016) studied the bubbles induced by breaking waves, and observed that the bubble size spectrum is proportional to 
$r^{-10/3}$, where 
$r$ is the effective radius of the bubbles. This power law of bubble size spectrum was explained by Garrett et al. (Reference Garrett, Li and Farmer2000) using an energy cascade theory. Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) performed direct numerical simulations (DNS) of a canonical three-dimensional two-phase viscous turbulent flow, with the turbulent kinetic energy (TKE) supplied by an underlying near-surface shear flow. They investigated the dependence of the air entrainment and bubble size on the Froude number and Weber number. Moreover, they proposed a heuristic model that qualitatively matched and explained the salient evolution behaviour of the bubble size spectrum. For some advanced ships, the plunging jet adds another source of air entrainment in the transom region and beyond (Hsiao et al. Reference Hsiao, Wu, Ma and Chahine2013). Hendrickson et al. (Reference Hendrickson, Weymouth, Yu and Yue2019) performed implicit large-eddy simulations (LES) with a conservative volume-of-fluid interface-capturing method to investigate the mixed-phase turbulent wake behind a dry transom stern of a surface ship. They conducted a detailed analysis on the air entrainment to investigate air entrainment parameters such as the entrainment rate and bubble size spectrum.
The above reviews indicate that the mechanism and physical process of air entrainment and bubble generation have been studied extensively in mixed-phase turbulence. However, studies on the statistical characteristics and transport mechanics of turbulence are limited in the literature. As the most common phenomenon that induces two-phase turbulence in nature, breaking waves modulate the transfer of mass, momentum and energy between the ocean and atmosphere. Deike (Reference Deike2022) summarized the recent research on breaking waves. Based on canonical experiments (Rapp & Melville Reference Rapp and Melville1990; Melville Reference Melville1994; Melville, Veron & White Reference Melville, Veron and White2002; Banner & Peirson Reference Banner and Peirson2007; Drazen, Melville & Lenain Reference Drazen, Melville and Lenain2008; Tian, Perlin & Choi Reference Tian, Perlin and Choi2010), DNS (Chen et al. Reference Chen, Kharif, Zaleski and Li1999; Iafrati Reference Iafrati2009; Deike, Popinet & Melville Reference Deike, Popinet and Melville2015; Deike et al. Reference Deike, Melville and Popinet2016; Wang et al. Reference Wang, Yang and Stern2016; Yang, Deng & Shen Reference Yang, Deng and Shen2018; Chan et al. Reference Chan, Johnson, Moin and Urzay2021; Mostert, Popinet & Deike Reference Mostert, Popinet and Deike2022) and LES (Lubin & Glockner Reference Lubin and Glockner2015; Derakhti & Kirby Reference Derakhti and Kirby2016), the dynamics of wave breaking were investigated. However, since the turbulence induced by breaking waves is a statistically unsteady process, it is infeasible to apply time averaging to define turbulent statistics. To date, research on turbulent statistics corresponding to breaking waves is limited to its impact on wind turbulence (Yang et al. Reference Yang, Deng and Shen2018), and the study on the mixed-phase region is limited. In terms of statistically stationary mixed-phase turbulence, Deshpande et al. (Reference Deshpande, Trujillo, Wu and Chahine2012) investigated the mean velocity, mean volume of fluid and TKE of a plunging jet. Yu et al. (Reference Yu, Hendrickson, Campbell and Yue2019) investigated the characteristics of free-surface turbulence and elucidated the qualitatively distinct characteristics of strong versus weak free-surface turbulence. Strong anisotropy was found in weak free-surface flow, while mixed-phase turbulence became almost isotropic in strong free-surface turbulence. Regarding the statistical averaging of variable density fluid motion equations, an additional unclosed term appears (Taulbee & Vanosdol Reference Taulbee and Vanosdol1991), that is, the turbulent mass flux (TMF). Modelling of TMF was investigated previously in compressible flow (Jones Reference Jones1979; Grasso & Speziale Reference Grasso and Speziale1989; Nichols Reference Nichols1990). Comprehensive reviews can be found in Chassaing (Reference Chassaing2001). Recently, Hendrickson & Yue (Reference Hendrickson and Yue2019) studied the incompressible highly variable density turbulence in the wake behind a three-dimensional dry transom stern. They developed an explicit algebraic closure model for the TMF term.
The objective of the present study is to investigate the statistical characteristics of the mixed-phase turbulence with air entrainment caused by a plunging jet. Specifically, we aim to study in detail the basic statistics reflecting the turbulent characteristics, such as TKE, TMF and their transport equations. The Froude number effect is examined through different cases. The closure problem of TMF is also discussed. The remainder of this paper is organized as follows. In § 2, the numerical method and physical set-up of the present simulation of the plunging jet is described. Then the results are presented and discussed in § 3. The conclusions are provided in § 4.
2. Details of the numerical simulations
2.1. Numerical method
To study a water jet plunging into a quiescent pool, we performed high-resolution LES by solving the three-dimensional two-phase incompressible Navier–Stokes equations. The coupled level-set (LS) and volume-of-fluid (VOF) method was used to capture the air–water interface on a Cartesian grid. The three-dimensional incompressible Navier–Stokes equations with varying densities and viscosities are expressed as
$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0 , \end{gather}$$
$$\begin{gather}\frac{\partial(\rho \boldsymbol{u})}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho \boldsymbol{u} \boldsymbol{u}) =-\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot}(2 \mu \boldsymbol{S}) + \rho \boldsymbol{g} + \boldsymbol{f}_s , \end{gather}$$
where 
$\rho$ and 
$\mu$ are the density and dynamic viscosity, respectively, 
$\boldsymbol {u} = [u, v, w]$ is the velocity, 
$p$ is the pressure, 
$\boldsymbol {S} = (\boldsymbol {\nabla } \boldsymbol {u} + \boldsymbol {\nabla } \boldsymbol {u}^{\rm T}) / 2$ is the strain-rate tensor, and 
$\boldsymbol {g} = [0, -g, 0]$ is the gravitational acceleration. The last term in (2.2) is the surface tension force, which is defined as
\begin{equation} \boldsymbol{f}_s = \sigma \kappa\,\boldsymbol{\nabla} H(\phi), \end{equation}
where 
$\sigma$ is the surface tension, 
$\kappa = \boldsymbol {\nabla }\boldsymbol {\cdot } (\boldsymbol {\nabla } \phi )$ at 
$\phi = 0$ is the interface curvature, with 
$\phi$ being the LS function, and 
$H(\phi )$ is the Heaviside function, defined as
\begin{equation} H(\phi ) = \begin{cases} 0, & \phi \leq 0, \\ 1, & \phi > 0 . \end{cases} \end{equation}In LES, the dynamic viscosity is expressed as
\begin{equation} \mu = \rho (\nu + \nu_t), \end{equation}
where 
$\nu$ is the kinematic viscosity, and 
$\nu _t$ is the subgrid-scale (SGS) eddy viscosity. In the present work, the SGS model proposed by Vreman (Reference Vreman2004) is used to determine 
$\nu _t$ for cases with high Reynolds numbers. At low Reynolds numbers, 
$\nu _t$ is omitted, following Hendrickson et al. (Reference Hendrickson, Weymouth, Yu and Yue2019).
The interface between the two fluid phases is captured using the coupled LS and VOF method. The following convection equations of the LS function 
$\phi$ and the VOF function 
$\psi$ are solved:
$$\begin{gather} \frac{\partial \phi}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\phi \boldsymbol{u})=0, \end{gather}$$
$$\begin{gather}\frac{\partial \psi}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\psi \boldsymbol{u})=0. \end{gather}$$
The LS function 
$\phi$ is defined as the signed distance from the fluid phase to the interface. Its sign is negative and positive for air and water, respectively. The VOF function 
$\psi$ is defined as the volume fraction of water in a grid cell ranging from 0 to 1. Subscripts 
$a$ and 
$w$ denote air and water, respectively. The density and viscosity are determined using the LS function as follows:
$$\begin{gather} \rho=\rho_a+\left(\rho_w-\rho_a\right) H(\phi), \end{gather}$$
$$\begin{gather}\mu=\mu_a+\left(\mu_w-\mu_a\right) H(\phi). \end{gather}$$As noted in the literature (Rudman Reference Rudman1998; Arrufat et al. Reference Arrufat, Crialesi-Esposito, Fuster, Ling, Malan, Pal, Scardovelli, Tryggvason and Zaleski2018; Nangia et al. Reference Nangia, Griffith, Patankar and Bhalla2019; Yang, Lu & Wang Reference Yang, Lu and Wang2021), if the density is determined solely using (2.8), then the simulation is unstable for two-fluid flows with a high-density contrast. A useful approach for improving the numerical stability is to calculate the mass and momentum fluxes using a consistent scheme and evolve the density by solving the following convection equation:
\begin{equation} \frac{\partial \rho}{\partial t}+\boldsymbol{\nabla} \boldsymbol{\cdot}(\rho \boldsymbol{u})=0. \end{equation}In the present study, (2.8) is used to determine the density at the beginning of each time step, while (2.10) is evolved together with the momentum equation to provide the density within a time step. Therefore, the interface is captured accurately, and the simulation is stable. More details regarding the numerical method can be found in Yang et al. (Reference Yang, Lu and Wang2021).
2.2. Physical set-up
Figure 1 shows schematically the computational domain and definition of the key parameters. As shown, a water jet spurts out horizontally and plunges into a quiescent water pool. Except for the water column and pool, other spaces of the computational domain are initially filled with air, some of which tends to be entrained into the water with the jet and then evolves into air cavities and bubbles beneath the free surface.
Figure 1. Definition sketch of the computational model.
We use the jet diameter 
$D = 0.04\ \text {m}$ as the characteristic length scale, and the horizontal outlet velocity of the jet 
$U = 4\ \text {m}\ \text {s}^{-1}$ as the characteristic velocity scale. The density of the water 
$\rho _w = 1.0 \times 10^3\ \text {kg}\ \text {m}^{-3}$ is used as the density unit. Hereafter, all variables are non-dimensionalized using 
$U$, 
$D$ and 
$\rho _w$, unless stated otherwise. The computational domain size is 
$L_x \times L_y \times L_z = 75.0 \times 27.0 \times 13.0$, where 
$x$, 
$y$ and 
$z$ represent the streamwise, vertical and spanwise directions of the domain, respectively. The water depth is 
$h = 18.75$. The distance between the centre of the jet orifice and the free surface is 
$\delta = 1.0$. The boundary condition indicates no penetration at the top, bottom, front and back of the domain, while a constant inlet velocity is applied in the jet orifice area of the left boundary, and a zero gradient condition is specified at the right outlet boundary.
The entire computational domain is discretized using a uniform Cartesian grid. The number of grid points is 
$1125 \times 405 \times 195$, giving grid resolution 
$\Delta x = \Delta y = \Delta z = 0.067$. The jet diameter is resolved by 15 grid points. Because the strong shear around the plunging points tears the interface, droplets and bubbles at very small sizes can be generated. To capture droplets and bubbles at all scales using grids is infeasible with the current computer power. There is no guarantee that the instantaneous dynamics of small droplets and bubbles is realistic in the present simulations. However, although small-scale dynamics of bubbles is not captured fully, turbulence statistics dominated by large-scale motions can be obtained accurately. This point is verified in numerical studies of single-phase turbulence. To ensure the reliability of the simulation results, this grid resolution is chosen through a careful mesh-convergence test. It is found that a grid resolution of 
$1/10$ jet diameter is sufficient to make accurate prediction of the mean velocity and TKE, while the TMF requires a finer resolution of 
$1/15$ jet diameter to converge. The details of the mesh-convergence test are given in the Appendix. Furthermore, we have also run all cases using a lower resolution of 
$1/10$ jet diameter. We find that the simulation with lower resolution reaches the same conclusions of this paper. We present the results obtained from the finer resolution of 
$1/15$ jet diameter to provide more accurate reference data for future research.
To consider the effect of the Reynolds number 
${Re} = {\rho _w UD}/{\mu _w}$ and the Froude number 
${Fr} = \sqrt {{U^2}/{gD}}$, we conducted nine cases in the present study. The air–water density ratio and viscosity ratio are 
$\rho _a/\rho _w = 1.2\times 10^{-3}$ and 
$\mu _a/\mu _w = 1.54\times 10^{-2}$, respectively. The Weber number defined based on the water density and the surface tension between the air and the water, 
${We} = \rho _w U^2 D / \sigma$, is set to a constant. The key parameters are listed in table 1. Cases 1, 2 and 6 were performed to examine the effect of the Reynolds number, and cases 3–9 were conducted to investigate the Froude number effect. The Reynolds number and Froude number for case 1 remained the same as the case considered in Deshpande et al. (Reference Deshpande, Trujillo, Wu and Chahine2012) to facilitate validation. At the lowest Reynolds number, 
$Re = 1.6\times 10^3$, the SGS eddy viscosity is smaller than the fluid viscosity except for a thin layer around the interface, where a cubic upwind interpolation (CUI) scheme is used to calculate the momentum flux to ensure the stability of the simulation (Yang et al. Reference Yang, Lu and Wang2021). The CUI scheme provides additional numerical dissipation near the interface. If the SGS eddy viscosity is added in this region, then the simulation is found to be overly dissipative. Therefore, the SGS eddy viscosity is omitted in cases 3–9, following Hendrickson et al. (Reference Hendrickson, Weymouth, Yu and Yue2019) to perform implicit LES.
Table 1. Parameters in the simulations of the plunging jet.
3. Results and discussion
3.1. Formation of the mixed-phase region
An instantaneous flow field under the statistically steady state is shown in figure 2, in which the interface between air and water is visualized using the iso-surface 
$\phi = 0$. This figure illustrates entrained air pockets and bubbles under the surface, and water splash and droplets above the surface. A large air pocket is entrained into the water near the plunging point, and it breaks into small air pockets and bubbles downstream. Above the free surface, water splashes and droplets hit the downstream surface, causing secondary plunging. The plunging event results in violent breaking of the free surface and highly mixed air–water turbulent flow. Movies for cases 3, 6 and 9 are provided as supplement materials, which are available at https://doi.org/10.1017/jfm.2023.1081.
Figure 2. Simulation result of instantaneous air–water interface visualized using the iso-surface 
$\phi = 0$. The result for case 6 is shown.
To better understand the mixed-phase turbulence in the near-surface region, we followed Hendrickson & Yue (Reference Hendrickson and Yue2019) to define a mixed-phase region as the variable density region, in which the mean volume of the fluid 
$\bar {\psi }$ satisfies 
$0.05 \leq \bar {\psi } \leq 0.95$. Here, the overline defines the time averaging, which is performed over a time duration 
$T = 1200.0$, after the turbulence is fully developed. The sampling rate is 
$\Delta T = 1.0$, which provides 1200 samples for time averaging. Figure 3 shows the mixed-phase region and mean free surface in case 6. As the streamwise coordinate 
$x$ increases, the size of the mixed-phase region increases and reaches a peak shortly after the jet plunging point, and then decreases downstream. The mean free surface with 
$\bar {\psi } = 0.5$ is also shown in figure 3(a) using the dash-dotted line. There exist a hollow mean free surface near the jet plunging point and a hump shortly downstream. They correspond to air entrainment and water splash-up, respectively.
Figure 3. Mixed-phase region in case 6: (a) contour of the mean volume of the fluid 
$0.05 < \bar {\psi } < 0.95$ at the mid-span; (b) transverse cuts of the mixed-phase region in different streamwise positions. The dash-dotted line in (a) and iso-surface in (b) represent the mean free surface with 
$\bar {\psi } = 0.5$.
3.2. The Reynolds number effect
The Reynolds number is an important parameter in turbulent flow. However, in many previous studies of mixed-phase turbulence (Brocchini & Peregrine Reference Brocchini and Peregrine2001; Deike et al. Reference Deike, Popinet and Melville2015; Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019), it was found that the Reynolds number effect was less significant than the Froude number effect. To minimize the effect of LES modelling, it is common to reduce the Reynolds number. In this study, we tested three Reynolds numbers to examine their effects on turbulent statistics. Meanwhile, we compared these results with previous experimental and numerical studies to validate our simulations.
Figure 4 shows the vertical profiles of the mean volume of the fluid 
$\bar {\psi }$ and mean streamwise velocity 
$\bar {u}$ at the mid-span and different streamwise locations for the three Reynolds numbers. As shown in figure 4(a), the mean volume of the fluid 
$\bar {\psi }$ varies mainly inside the mixed-phase region. The results for 
$\bar {\psi }$ for different Reynolds numbers are close to each other. They also agree with the numerical results of Deshpande et al. (Reference Deshpande, Trujillo, Wu and Chahine2012). Figure 4(b) shows that the results for the mean velocity for different Reynolds numbers are also close to each other. This observation indicates that variation in the Reynolds number (from 
$1.6 \times 10^3$ to 
$1.6 \times 10^5$) does not impose significant effects on the mean flow. The results for mean velocity also agree with the experimental and numerical results of Deshpande et al. (Reference Deshpande, Trujillo, Wu and Chahine2012).
Figure 4. Vertical profiles of (a) the mean volume of the fluid and (b) the mean streamwise velocity at the mid-span and different streamwise locations for different Reynolds numbers. Solid line indicates 
${Re} = 1.6 \times 10^5$, dashed line indicates 
${Re} = 1.6 \times 10^4$, dotted line indicates 
${Re} = 1.6 \times 10^3$; circles and triangles represent the numerical and experimental data, respectively, of the same case with 
${Re} = 1.6 \times 10^5$ in Deshpande et al. (Reference Deshpande, Trujillo, Wu and Chahine2012). Dash-dotted lines represent the edge of the mixed-phase region. The results for case 6 are shown.
We also calculated the time-averaged bubble-size density spectra 
$\bar {N}(r_{eff})$, defined as
\begin{equation} \bar{N}\left( r_{eff} \right) = \frac{1}{T} \int_t^{t+T} \frac{n\left(r_{eff}, t ; b\right)}{V \times b}\,\mathrm{d} t, \end{equation}
where 
$n(r_{eff}, t ; b)$ is the number of bubbles, whose effective radii fall between 
$r_{eff}$ and 
$r_{eff} + b$ in a given fluid volume 
$V$ at time 
$t$. In the present work, 
$b=0.001$ is chosen. The fluid volume 
$V$ for bubble statistics is a cuboid in the computational domain 
$x \in [15, 45]$, 
$y \in [0, h]$, 
$z \in [0, L_z]$, which contains most of the air cavities and bubbles beneath the interface. The effective spherical radius is defined as
\begin{equation} r_{eff} = \left[(3 / 4 {\rm \pi}) v_e\right]^{1 / 3}, \end{equation}
where 
$v_e$ is the volume of an individual bubble. To determine the number and volume of each bubble, a connected component algorithm (Samet & Tamminen Reference Samet and Tamminen1988) is used to identify and label the entrained air cavities. Figure 5 shows the results for the bubble-size density spectra 
$\bar {N}(r_{eff})$ for different Reynolds numbers. It can be seen that the results for different Reynolds numbers are close to each other, indicating that the Reynolds number effect on the bubble-size density spectra is negligible. The solid line in figure 5 represents the 
$r^{-10/3}$ power law, which is satisfied in cases at different Reynolds numbers.
Figure 5. Bubble-size density spectra for cases at different Reynolds number. The Froude number is 
${Fr} = 6.4$.
The results for mean velocity, volume of the fluid, and bubble-size density spectra for different Reynolds numbers (ranging from 
$1.6 \times 10^3$ to 
$1.6 \times 10^5$) indicate that the Reynolds number imposes a limited impact on the turbulence statistics. We also examined the Reynolds number effects on other turbulent statistics, including TKE and TMF. The impact of the Reynolds number is found to be less significant than that of the Froude number. Therefore, in the following context of this paper, we focus on the effect of the Froude number.
3.3. Cross-sectional area of the mixed-phase region
We start the analyses of the Froude number effect from the size of the mixed-phase region. Figure 6 shows the streamwise variation of the cross-sectional area 
$S^R(x)$ of the mixed-phase region. The streamwise location of jet plunging point 
$x^p$ varies with the Froude number. We use 
$x - x^p$ as an independent variable to facilitate comparisons among different cases. As shown in figure 6, the maximum cross-sectional area increases as the Froude number increases. Figure 7 shows the variation in the maximum cross-sectional area of the mixed-phase region, 
$\max (S^R)$, with respect to the Froude number. It is seen that 
$\max (S^R)$ increases linearly with the Froude number. The observations from figures 6 and 7 indicate that as the Froude number increases, the mixing of air and water is enhanced. A similar conclusion was drawn in previous studies of plunging jets (Chirichella et al. Reference Chirichella, Gomez Ledesma, Kiger and Duncan2002; Kiger & Duncan Reference Kiger and Duncan2012) and other mixed-phase turbulent flows, such as hydraulic jumps (Chachereau & Chanson Reference Chachereau and Chanson2011; Ma et al. Reference Ma, Oberai, Drew, Lahey and Hyman2011) and mixed-phase turbulence induced by shear near the interface (Yu et al. Reference Yu, Hendrickson, Campbell and Yue2019).
Figure 6. Streamwise variation of the mixed-phase region area for cases at different Froude numbers.
Figure 7. Maximum mixed-phase region area at different Froude numbers, and the linear fitting line.
3.4. Mean velocity
Figure 8 shows the contours of the mean streamwise and vertical velocities at the mid-span for case 6 at an intermediate Froude number 
${Fr} = 6.4$. The dotted lines show the upper and lower edges of the mixed-phase region. Figure 8(a) shows that the jet plunging induces a mean streamwise velocity 
$\bar {u}$ near the free surface, and the large magnitude of 
$\bar {u}$ is collocated with the mixed-phase region. As shown in figure 8(b), the direction of the mean vertical velocity 
$\bar {v}$ varies along the streamwise direction in the mixed-phase region. Near the jet plunging point, the fluid around the jet moves downwards with it. Shortly downstream, the pool water moves upwards, and droplets are generated. Meanwhile, the air cavities under the surface move upwards under buoyancy. As a result, the vertical velocity is positive in this region. After the droplets reach the peak, they fall back to the pool and form a secondary plunging, which causes another region with a negative vertical velocity.
Figure 8. Contours of time-averaged velocities (a) 
$\bar {u}$ and (b) 
$\bar {v}$ at the mid-span. Dotted lines represent the edges of the mixed-phase region. Here, 
${Re} = 1600$ and 
${Fr} = 6.4$.
Because the size of the mixed-phase region varies with the Froude number, to facilitate comparison of the results for different cases, we followed Hendrickson & Yue (Reference Hendrickson and Yue2019) to define the conditioned average in the mixed-phase region as
\begin{equation} \langle \,f^R \rangle _{y z}(x) = \frac{1}{S^R(x)} \int_{y z} f^R(y, z ; x) \,\mathrm{d} S^R. \end{equation}
Here, 
$f^R$ represents the variable 
$f$ inside the mixed-phase region. The integration denoted by 
$\langle \cdot \rangle _{y z}$ is performed over a cross-stream section 
$(y$–
$z$ plane), and 
$S^R(x)$ represents the area of the mixed-phase region in the corresponding plane.
Figure 9 compares the mean velocities 
$\langle \bar {u}_i^R \rangle _{yz}$ averaged in the mixed-phase region for different Froude numbers. We note here that the results in the flow region for 
$x - x^p = 0.0\unicode{x2013}4.0$ show some uncertainty because the area of the mixed-phase region 
$S^R(x)$ is small. Therefore, the following analyses focus mainly on the remaining flow region for 
$x - x^p > 4.0$, where the sampling number is sufficiently large to provide more reliable statistics. This does not influence our understanding of the statistical properties of the mixed-phase turbulence induced by the plunging jet, because active turbulence occurs mainly downstream, where air and water are sufficiently mixed.
Figure 9. Streamwise variation of the mean velocities averaged in the mixed-phase region for different Froude numbers: (a) 
$\langle \bar {u}^R \rangle _{yz}$, (b) 
$\langle \bar {v}^R \rangle _{yz}$.
Figure 9(a) shows that 
$\langle \bar {u}^R \rangle _{yz}$ near the jet plunging region (
$x - x^p = 4.0\unicode{x2013}8.0$) remains almost unchanged at different Froude numbers. Noting that the outlet horizontal velocity of the jet remain the same in all cases, the above observation indicates that the mean flow is determined by the inflow in this region, while turbulent motion does not impose significant influences. At approximately 
$x-x^p=10.0$, the magnitude of 
$\langle \bar {u}^R \rangle _{yz}$ reaches a valley due to the drag effect of the pool water. Downstream, 
$\langle \bar {u}^R \rangle _{yz}$ shows a complex non-monotonic response to the increase in the Froude number. At lower Froude numbers 
${Fr} \leq 5.3$, 
$\langle \bar {u}^R \rangle _{yz}$ increases with the Froude number. Regarding 
${Fr} \geq 6.4$, there exists intense vertical motion that expands the size of the mixed-phase region, and as a result, the streamwise momentum is diffused and 
$\langle \bar {u}^R \rangle _{yz}$ decreases as the Froude number increases.
Figure 9(b) shows that a primary negative peak of 
$\langle \bar {v}^R \rangle _{yz}$ occurs at approximately 
$x-x^p=4.0$ shortly downstream of the plunging point. The magnitude of this negative-valued peak decreases as the Froude number increases because of the reduction of the gravitational potential energy of the injected water, which is proportional to 
${Fr}^{-1}$. Downstream of the plunging region, the vertical motion is weak at low Froude numbers 
${Fr} \leq 5.3$, resulting in a small magnitude of 
$\langle \bar {v}^R \rangle _{yz}$. Regarding 
${Fr} \geq 6.4$, the water plunging generates a large amount of droplets, which induce the first positive peak of 
$\langle \bar {v}^R \rangle _{yz}$. When the droplets reach the highest altitude, the vertical velocity becomes zero. Moreover, at the corresponding streamwise location, the size of the mixed-phase region reaches a maximum (figure 6). The secondary negative peak of 
$\langle \bar {v}^R \rangle _{yz}$ represents the downward motion of droplets, which leads to secondary plunging. Downstream of the secondary plunging, there is still a small magnitude of 
$\langle \bar {v}^R \rangle _{yz}$ at high Froude numbers. This indicates that the increasing Froude number results in more intense vertical motions of the surface.
To perform a statistical study of turbulent properties, an instantaneous variable 
$f$ is decomposed as 
$f(x, y, z; t) = \bar {f}(x, y, z) + f'(x, y, z; t)$, where 
$f'$ is the fluctuation. The mean momentum equation of incompressible variable-density flow can be expressed as
\begin{equation} \frac{\partial (\bar{\rho} \bar{u}_i)}{\partial t} = 0 = CM_i + GM^p_i + GM_i + DM^v_i + DM^t_i + AM_i. \end{equation}
Because all statistics are calculated at a statistical stationary stage of the flow, it is assumed that the time derivative vanishes owing to a sufficiently long time duration used for performing the time averaging. The budget terms on the right-hand side of (3.4) include the convection term 
$CM_i$, pressure gradient term 
$GM^p_i$, gravity term 
$GM_i$, viscous diffusion term 
$DM^v_i$, Reynolds stress term 
$DM^t_i$, and TMF term 
$AM_i$. These terms are defined as
$$\begin{gather} CM_i =- \frac{\partial\left(\bar{\rho} \bar{u}_i \bar{u}_j\right)}{\partial x_j}, \end{gather}$$
$$\begin{gather}GM^p_i =- \frac{\partial \bar{p}}{\partial x_i}, \end{gather}$$
$$\begin{gather}GM_i =-\frac{\bar{\rho}}{{Fr}^2}\,\delta_{i2}, \end{gather}$$
$$\begin{gather}DM^v_i = \frac{\partial \bar{\tau}_{i j}}{\partial x_j}, \end{gather}$$
$$\begin{gather}DM^t_i =- \frac{\partial \overline{\rho u_i^{\prime} u_j^{\prime}}}{\partial x_j}, \end{gather}$$
$$\begin{gather}AM_i =- \frac{\partial \overline{\rho u_i^{\prime}}}{\partial t}- \frac{\partial}{\partial x_j}\left(\overline{\rho u_i^{\prime}}\,\bar{u}_j + \overline{\rho u_j^{\prime}}\, \bar{u}_i\right), \end{gather}$$
where 
$\tau _{ij} = \mu ( u_{i, j} + u_{j, i} ) / {Re}$ is the viscous stress tensor. There are unclosed terms on the right-hand side of (3.4), namely the Reynolds stress 
$\overline {\rho u_i^{\prime } u_j^{\prime }}$ and TMF 
$\overline {\rho u_i^{\prime }}$.
Figure 10 shows the streamwise variation of each budget term of (3.4) averaged in the mixed-phase region. As shown in figure 10(a), the convection term 
$CM_1$, Reynolds stress term 
$DM_1^t$ and TMF term 
$AM_1$ make dominant contributions to the transport of the mean streamwise momentum 
$\bar {\rho } \bar {u}$. The convection term 
$CM_1$ and Reynolds stress term 
$DM_1^t$ balance each other near the jet plunging point. They decay downstream, and the TMF term 
$AM_1$ becomes a dominant term. It is seen from figure 10(b) that among the budget terms of 
$\bar {\rho } \bar {v}$, the summation of the gravity term 
$GM_2$ and mean pressure gradient term 
$GM^p$ makes significant contributions. The TMF term 
$AM_2$ is important near the jet plunging point and decays to a small magnitude for 
$x - x^p \geq 10.0$. The Reynolds stress term 
$DM_2^t$ plays an important role for 
$x - x^p \geq 10.0$. The results shown in figure 10 indicate that the closure of both the Reynolds stress and TMF is important in the mixed-phase turbulence induced by jet plunging.
Figure 10. Streamwise variation of the budget term of the mean momentum equation averaged in the mixed-phase region. The results for 
${Re} = 1600$ and 
${Fr} = 6.4$ (case 6) are shown.
3.5. Turbulent kinetic energy
There are different strategies for closing the Reynolds stress 
$\overline {\rho u_i^{\prime } u_j^{\prime }}$. In single-phase flows, an important strategy is to use the dynamic equation of TKE 
$k = \frac {1}{2} \overline {\rho u_i^{\prime } u_i^{\prime }}$ for closure, such as the 
$k$–
$\varepsilon$ model (Chien Reference Chien1982; Kaul Reference Kaul2010, Reference Kaul2011) and the 
$k$–
$\omega$ model (Wilcox Reference Wilcox1988, Reference Wilcox2008; Menter Reference Menter1994; Spalart & Rumsey Reference Spalart and Rumsey2007). In the following context, we first analyse the effect of the Froude number on TKE, followed by some discussions on its closure model in the mixed-phase turbulence.
Figure 11 displays the contours of TKE at the mid-span for 
${Fr} = 6.4$. The figure demonstrates a collocation between the mixed-phase region and large magnitude of TKE. The highest TKE is observed near the jet plunging point (
$x - x^p = 4.0\unicode{x2013}12.0$), where the shear between the jet and the pool water is strong. At approximately 
$x - x^p = 24.0$, there is a secondary peak of TKE caused by the secondary plunging.
Figure 11. Contours of TKE at the mid-span. Dotted lines represent the edges of the mixed-phase region. The results for 
${Re} = 1600$ and 
${Fr} = 6.4$ (case 6) are shown.
Figure 12 compares the streamwise variation of TKE averaged in the mixed-phase region 
$\langle k^R \rangle _{yz}$ for different Froude numbers. At all Froude numbers, a primary peak occurs at approximately 
$x-x^p = 6.0$. Regarding large Froude numbers (
${Fr} \geq 6.4$), a secondary peak of TKE occurs. Figure 13 compares the magnitudes of the two peaks at different Froude numbers. At low Froude numbers, the magnitude of the primary peak increases slightly as the Froude number increases. As the Froude number exceeds 
${Fr} = 6.4$, the magnitudes of both primary and secondary peaks decrease linearly with an increasing Froude number. The observations from figure 13 indicate that the Froude number imposes dual effects on TKE. At low Froude numbers, the entrained air volume increases, and the shear is enhanced near the interface as the Froude number increases. As a result, the magnitude of TKE increases with the Froude number for 
${Fr} \leq 5.3$. At larger Froude numbers, 
${Fr} \geq 6.4$, the increase in the Froude number leads to a decrease in the gravitational potential energy. Consequently, the vertical velocity and plunging angle decreases when the jet hits the free surface. Furthermore, at higher Froude numbers, more droplets are generated, resulting in secondary plunging. These effects lead to the distribution of TKE in a wider streamwise range at higher Froude numbers, resulting in a smaller peak value.
Figure 12. Streamwise variation of TKE averaged in the mixed-phase region for various Froude numbers.
Figure 13. Variation of the magnitudes of the primary peak (red circle symbols) and secondary peak (blue square symbols) of TKE with respect to the Froude number.
The transport equation of TKE in variable-density flows is expressed as (Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002)
\begin{equation} \frac{\partial k}{\partial t} = 0 = CK + DK^t + PK + DK^p + DK^v + \varepsilon K + AK. \end{equation}
The budget terms on the right-hand side of (3.11) include the convection term 
$CK$, turbulence diffusion term 
$DK^t$, production term 
$PK$, pressure diffusion term 
$DK^p$, viscous diffusion term 
$DK^v$, dissipation term 
$\varepsilon K$, and the TMF correlation term 
$AK$. These terms are defined as follows:
$$\begin{gather} CK =-\frac{\partial\left(k \bar{u}_j\right)}{\partial x_j}, \end{gather}$$
$$\begin{gather}DK^t =- \frac{\partial\left(\tfrac{1}{2}\,\overline{\rho u_i^{\prime} u_i^{\prime} u_j^{\prime}}\right)}{\partial x_j}, \end{gather}$$
$$\begin{gather}PK =-\overline{\rho u_i^{\prime} u_j^{\prime}}\,\frac{\partial \bar{u}_i}{\partial x_j}, \end{gather}$$
$$\begin{gather}DK^p =- \frac{\partial \overline{p^{\prime} u_i^{\prime}}}{\partial x_i}, \end{gather}$$
$$\begin{gather}DK^v = \frac{\partial \overline{\tau_{i j}^{\prime} u_i^{\prime}}}{\partial x_j}, \end{gather}$$
$$\begin{gather}\varepsilon K =- \overline{\tau_{i j}^{\prime}\,\frac{\partial u_i^{\prime}}{\partial x_j}}, \end{gather}$$
$$\begin{gather}AK = \frac{\overline{\rho u_i^{\prime}}}{{Fr}^2}\,\delta_{i2} - \overline{\rho u_i^{\prime}}\left(\frac{\partial \bar{u}_i}{\partial t} + \bar{u}_j\,\frac{\partial \bar{u}_i}{\partial x_j}\right). \end{gather}$$ Figure 14 shows the streamwise variation of the budget terms of 
$\langle k^R \rangle _{yz}$ averaged in the mixed-phase region. The results for 
${Re} = 1600$ and 
${Fr} = 6.4$ are shown. Similar to many single-phase flows, the balance between the production term 
$PK$ and dissipation term 
$\varepsilon K$ dominates the transport of TKE. The convection term 
$CK$ is positive upstream and becomes negative downstream, indicating the energy convection from upstream to downstream. The magnitude of the TMF correlation term 
$AK$ is smaller than the production term 
$PK$. This indicates that the closure model of single-phase turbulence can be referenced by the mixed-phase turbulence induced by a plunging jet.
Figure 14. Streamwise variation of the budget terms in the transport equation of TKE averaged in the mixed-phase region. The results for 
${Re} = 1600$ and 
${Fr} = 6.4$ (case 6) are shown.
3.6. Turbulent mass flux and its transports
It is understood from (3.4)–(3.10) that there are two unclosed terms in the mean momentum equation of the mixed-phase turbulence. The first term is the Reynolds stress tensor, and it is usual to consider its isotropic part (i.e. TKE) for closure problems. The analysis of the TKE transport equation in § 3.5 shows that its budget remains similar to many single-phase turbulent flows. The energy generated by the TMF correlation term is relatively small. Thus the closure model single-phase flow can be adopted. However, § 3.4 shows that the TMF term plays an important role in the transport of the mean momentum. Therefore, its closure model also needs to be considered for industrial applications. Hendrickson & Yue (Reference Hendrickson and Yue2019) developed an algebraic model for TMF based on their LES data of the wake of a three-dimensional dry transom stern. Their a priori tests showed that it is challenging to obtain a high correlation coefficient between the model and LES data. Another strategy to close the TMF term is to develop a dynamic model, which requires an analysis of its transport equation. In this subsection, we investigate the effects of the Froude number on TMF, followed by its transport equation.
Figures 15(a) and 15(b) show the contours of 
$\overline {\rho u^\prime }$ and 
$\overline {\rho v^\prime }$, respectively, at the mid-span for 
${Re} = 1600$ and 
${Fr} = 6.4$. The upper and lower dotted lines represent the edges of the mixed-phase region. The dash-dotted line in figure 15(a) represents the mean location of the free surface corresponding to 
$\bar {\psi } = 0.5$. As shown in figure 15(a), 
$\overline {\rho u^\prime }$ is positive above the mean free surface, indicating the downstream transport of water droplets. Below the mean free surface, 
$\overline {\rho u^\prime }$ is negative valued, corresponding to the downstream transport of bubbles. Figure 15(b) shows that near both the primary and secondary plunging points, 
$\overline {\rho v^\prime }$ is positive, indicating air entrainment. Shortly downstream of the plunging points, 
$\overline {\rho v^\prime }$ changes its sign to negative, corresponding to bubble detrainment.
Figure 15. Contours of (a) streamwise TMF component 
$\overline {\rho u^\prime }$ and (b) vertical TMF component 
$\overline {\rho v^\prime }$ at the mid-span. The upper and lower dotted lines represent the edges of the mixed-phase region. The dash-dotted line in (a) represents the mean location of the free surface. The results for 
${Re} = 1600$ and 
${Fr} = 6.4$ (case 6) are shown.
Figure 16 compares the streamwise variation of TMF averaged in the mixed-phase region for different Froude numbers. Figure 16(a) shows that the negative value of 
$\overline {\rho u^\prime }$ dominates in the mixed-phase region. Recalling that negative and positive 
$\overline {\rho u^\prime }$ correspond to downstream convection of bubbles and droplets, respectively, the negative value of 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$ indicates that the downstream transport of air beneath the free surface is dominant. The primary peak of the negative-valued 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$ occurs downstream of the plunging point. Its magnitude increases as the Froude number increases. This indicates that more bubbles are convected downstream at higher Froude numbers. When 
${Fr} \geq 6.4$, there exists a secondary peak in 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$. However, it shows a different trend of the secondary peak as the Froude number increases. This is because the convection of droplets above the surface balances part of the bubble motion beneath the surface. At large Froude numbers, water splash-up induces droplets, resulting in the decrease in the magnitude of the secondary negative-valued peak of 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$.
Figure 16. Streamwise variation of TMF 
$\langle \overline {\rho u_i^\prime }^R\rangle _{yz}$ averaged in the mixed-phase region for various Froude numbers: (a) streamwise component and (b) vertical component.
Figure 16(b) shows that 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$ is positive near the plunging point, indicating air entrainment in this region. The magnitude of 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$ is small at 
${Fr} = 3.2$ and 4.2. It increases as the Froude number increases from 4.2 to 6.4. As the Froude number continues to increase, its magnitude decreases slightly, but the streamwise range with positive 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$ expands. This indicates that the air entrainment takes place in a larger streamwise region at a higher Froude number. A negative peak of 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$ occurs for 
${Fr} \geq 6.4$, caused by the bubble detrainment after the secondary plunging. The magnitude of this negative-valued peak of 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$ increases as the Froude number increases from 6.4 to 9.6, indicating that more bubbles are detrained at higher Froude numbers.
To investigate the closure of TMF, we examined the following transport equation of TMF:
\begin{equation} \frac{\partial \overline{\rho u_i^{\prime}}}{\partial t} = 0 = CT_i + PT^{(1)}_i + PT^{(2)}_i + DT_i + ET_i. \end{equation}
The budget terms on the right-hand side include the convection term 
$CT_i$, the production terms 
$PT^{(1)}_i$ and 
$PT^{(2)}_i$ corresponding to the velocity gradient and density gradient, respectively, the turbulent diffusion term 
$DT_i$, and a combining term 
$ET_i$. The definitions of these terms are given as follows:
$$\begin{gather} CT_i =-\frac{\partial \left(\overline{\rho u_i^{\prime}}\,\bar{u}_j \right)}{\partial x_j}, \end{gather}$$
$$\begin{gather}PT^{(1)}_i =-\overline{\rho u_j^{\prime}}\,\frac{\partial \bar{u}_i}{\partial x_j}, \end{gather}$$
$$\begin{gather}PT^{(2)}_i =-\overline{u_i^{\prime} u_j^{\prime}}\,\frac{\partial \bar{\rho}}{\partial x_j}, \end{gather}$$
$$\begin{gather}DT_i = \frac{\partial \bar{\rho}\,\overline{u_i^{\prime} u_j^{\prime}}}{\partial x_j} - \frac{\partial \overline{\rho u_i^{\prime} u_j^{\prime}}}{\partial x_j}, \end{gather}$$
$$\begin{gather}ET_i = \left(\bar{\rho}\,\overline{\frac{1}{\rho}\,\frac{\partial p}{\partial x_i}} - \frac{\partial \bar{p}}{\partial x_i}\right) + \left(\frac{\partial \bar{\tau}_{i j}}{ \partial x_j} - \bar{\rho}\,\overline{\frac{1}{\rho}\,\frac{\partial \tau_{i j}}{\partial x_j}}\right). \end{gather}$$
The combining term 
$ET_i$ consists of a pressure gradient part and a viscous stress part. From its expression, it is understood that this term is essentially the difference between their time-averaged values (i.e. 
${\partial \bar {p}}/{\partial x_i}$ and 
${\partial \overline {\tau }_{i j}}/{\partial x_j}$) and their density-weighted time-averaged values (i.e. 
$\bar {\rho } \overline{{\rho }^{-1}\,{\partial p}/{\partial x_i}}$ and 
$\bar {\rho }\, \overline {{\rho }^{-1}\,{\partial \tau _{i j}}/{\partial x_j}}$). To perform density-weighted time averaging, the instantaneous density needs to be interpolated. Due to the use of sharp-interface treatment in the present simulation, the instantaneous density 
$\rho$ varies sharply across the interface. As a result, the interpolation of 
$\rho$ causes oscillations near the interface, resulting in non-physical values. Therefore, in the present study, the pressure gradient and viscous stress terms are combined into one term, 
$ET_i$, and its value is determined indirectly as the opposite number of the summation of other terms. This treatment is supported by the assumption that 
${\partial \overline {\rho u_i^{\prime }}}/{\partial t} = 0$ is attained, because a sufficiently long time duration is used for performing time averaging.
Figure 17 shows the streamwise variation of the budget terms of TMF averaged in the mixed-phase region for 
$Re=1600$ and 
$Fr=6.4$. Figure 17(a) shows that all budget terms of 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$ are active near the jet plunging point. Downstream, the convection term 
$CT_1$, production term corresponding to the velocity gradient 
$PT^{(1)}_1$ and turbulent diffusion term 
$DT_1$ decay to a relatively small value, and the budget is balanced mainly by the production term corresponding to the density gradient 
$PT^{(2)}_1$ and combining term 
$ET_1$. From figure 17, it is seen that the balance between the production term corresponding to the density gradient 
$PT^{(2)}_2$ and combining term 
$ET_2$ also dominates the budget of the vertical component of TMF 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$.
Figure 17. Streamwise variation of the budget terms in transport equation of (a) the streamwise component and (b) the vertical component of TMF averaged in the mixed-phase region. The results for 
${Re} = 1600$ and 
${Fr} = 6.4$ (case 6) are shown.
Regarding the closure of (3.19), the convection term 
$CT_i$ and production term corresponding to the velocity gradient 
$PT^{(1)}_i$ do not require modelling. The turbulent diffusion term 
$DT_i$, which can be seen as a diffusion effect of velocity fluctuations on TMF, is important near the plunging point, where the turbulent fluctuations are strong. This term can be modelled by estimating a characteristic diffusion velocity through TKE. The combining term 
$ET_i$ is an important transportation term, but due to the lack of reliable data and deeper understanding, it is currently difficult to develop a rational closure model for this term. Since both the pressure and viscous stress fluctuations are induced by the velocity fluctuations, 
$ET_i$ can be seen as a passive response of the flow field to the change in the other budget terms of TMF, and possibly it can be modelled as a diffusion effect corresponding to an artificial viscosity. The production term 
$PT^{(2)}_i$ is related to the density gradient caused by the two-phase mixture, which is an important source of TMF. Therefore, it is crucial to close 
$PT^{(2)}_i$ in a dynamic model of TMF.
Figure 18 compares the streamwise variation of 
$PT^{(2)}_i$ averaged in the mixed-phase region at different Froude numbers. At lower Froude numbers (
${Fr} \leq 5.3$), 
$\langle PT_1^{(2) R} \rangle _{yz}$ shows a sharp peak near the jet impact point, and its absolute value increases as the Froude number increases. At higher Froude numbers (
${Fr} \geq 6.4$), 
$\langle PT_1^{(2) R} \rangle _{yz}$ shows two peaks corresponding to the two plunging events, and their magnitudes both decrease with an increasing Froude number. From the comparison between figures 18(a) and 18(b), it is understood that despite the opposite sign, 
$\langle PT_2^{(2) R} \rangle _{yz}$ is similar to 
$\langle PT_1^{(2) R} \rangle _{yz}$ in terms of both the variation in the streamwise direction and the Froude number effect.
Figure 18. Streamwise variation of the production term corresponding to the density gradient 
$PT_i^{(2)}$ averaged in the mixed-phase region for various Froude numbers: (a) streamwise component and (b) vertical component.
The consistency between the production term 
$PT^{(2)}_i$ and TKE is better than that between TMF and TKE. Specifically, TKE shows two peaks at large Froude numbers (
${Fr} \geq 6.4$), and the magnitudes of both peaks decrease as the Froude number increases (figure 12). Although TMF also shows two peaks at large Froude numbers (
${Fr} \geq 6.4$), their magnitudes show different trends as the Froude number increases (figure 16). Furthermore, the secondary peak value of 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$ is negative, while both peaks of 
$\langle PT_2^{(2) R} \rangle _{yz}$ are positive. This indicates that 
$PT_i^{(2)}$ can be modelled through TKE.
Based on the above analyses, we propose a closure model for 
$PT^{(2)}_i$ as
\begin{equation} PT_i^{(2)} =-C_i\,\overline{u_k^{\prime} u_k^{\prime}}\,\frac{\partial \bar{\rho}}{\partial x_i}, \end{equation}
where 
$C_i$ is the model coefficient. A linear least squares fit between 
$PT_i^{(2)}$ and 
$-\overline {u_k^{\prime } u_k^{\prime }}\,{\partial \bar {\rho }}/{\partial x_i}$ in the mixed-phase region is used to determine the model coefficient. The model coefficients and conditioned correlation coefficients 
$R_i$ between the two sides of (3.25) at different Froude numbers are listed in table 2. The coefficients of the streamwise and vertical components of (3.25) are displayed, while the spanwise component is not shown because the spanwise component of TMF is much smaller than the other two components. The correlation coefficient of the vertical component exceeds 0.89, and the model coefficient 
$C_2$ varies little with the Froude number. This indicates that 
$PT_2^{(2)}$ can be well estimated by the proposed model. Although the accuracy of the streamwise component estimated by the model is not as good as that of the vertical component, the correlation coefficients of the streamwise components are all above 0.55, which is overall higher than the correlation coefficient of the model that directly fits TMF using TKE in Hendrickson & Yue (Reference Hendrickson and Yue2019).
Table 2. Model coefficients in (3.25) and conditioned correlation coefficients between the two sides of (3.25) in the mixed-phase region. The two columns of 
$C_i$ and 
$R_i$ correspond to the streamwise and vertical components, respectively.
In the above discussion, the Reynolds stress (without density) 
$\overline {u_i^{\prime } u_j^{\prime }}$ is modelled using TKE (without density) 
$\overline {u_k^{\prime } u_k^{\prime }}$. In single-phase turbulence, the Reynolds stress is modelled conventionally by the strain-rate tensor 
$\bar {S}_{i j} = ({\partial \bar {u}_i}/{\partial x_j} + {\partial \bar {u}_j}/{\partial x_i})/2$ via an eddy viscosity. In view of this, we also examine the feasibility of incorporating an eddy viscosity model in the closure of 
$PT^{(2)}_i$. In this situation, the closure model for 
$PT^{(2)}_i$ is expressed as
\begin{equation} PT_i^{(2)} = 2\nu_{e} \bar{S}_{i j}\,\frac{\partial \bar{\rho}}{\partial x_j}, \end{equation}
where 
$\nu _e$ is the eddy viscosity. We calculate the correlation coefficients between the two sides of (3.26), which are listed in table 3. By contrasting tables 2 and 3, it is understood that the correlation coefficients between the two sides of (3.26) are smaller than those of (3.25), indicating that (3.25) is more appropriate for modelling 
$PT^{(2)}_i$ in the mixed-phase turbulence generated by a plunging jet.
Table 3. Correlation coefficients between the two sides of (3.26) in the mixed-phase region. The two columns of 
$R_i$ indicate the streamwise and vertical components, respectively.
4. Conclusion
In the present study, we performed high-resolution interface-resolved LES to study the mixed-phase turbulence induced by a water jet plunging into a quiescent pool. Among all simulation cases, the Reynolds number ranged from 
$1.6 \times 10^3$ to 
$1.6 \times 10^5$, and the Froude number ranged from 
$3.2$ to 
$9.6$. Because the effect imposed by the Reynolds number on the turbulent statistics was found to be less significant than that imposed by the Froude number, we focused mainly on the effects of the Froude number on the turbulent statistics in this paper. The simulation results showed that increasing the Froude number led to an increase in the area of the mixed-phase region. To facilitate direct comparison of turbulent statistics among different cases, a conditioned average inside the mixed-phase region was adopted. The main findings of this study are summarized below.
The mean velocity averaged in the mixed-phase region varies non-monotonically with the Froude number. Among all Froude numbers under investigation, the magnitude of the mean streamwise velocity reaches its maximum at an intermediate value 
$Fr=6.4$. The mean vertical velocity shows a single negative-valued peak for low Froude numbers with 
$Fr \geq 5.3$. As the Froude number increases to 
$Fr\ge 6.4$, water splash-up and secondary plunging occur, causing a positive-valued peak and a secondary negative-valued peak. The complex behaviour of the mean velocity is correlated to the nonlinear effects associated with the mixed-phase turbulence, which requires an accurate closure model in industrial applications.
To discuss the closure problem of the Reynolds-averaged mean momentum equation, we analysed TKE, TMF and their transport equations. Our simulations show that TKE also varies non-monotonically with an increasing Froude number. At low Froude numbers, TKE shows a single peak near the plunging point. The magnitude of this peak increases as the Froude number increases from 3.2 to 5.3. As the Froude number increases to 
$Fr\ge 6.4$, the secondary plunging induces a secondary peak in TKE. The magnitudes of both the primary and secondary peaks of TKE decrease as the Froude increases from 6.4 to 9.6. The analysis of the TKE transport equation shows that it is dominated by the balance between the production and dissipation terms. The convection and turbulent diffusion terms are mainly responsible for the spatial transport of the TKE. Although there exists a TMF correlation term in the TKE transport equation, its magnitude is smaller than the above dominant terms, indicating that the closure model of TKE for single-phase turbulence can be used in the mixed-phase turbulence induced by the plunging jet.
Compared to the single-phase turbulence, the additional TMF term occurs in not only the TKE transport equation, but also the mean momentum equation. Although the importance of the TMF correlation term is found to be insignificant in the TKE transport equation, its influence on the mean momentum transport is found to be important. The streamwise component of TMF is positive above the mean water elevation, and negative below the mean water elevation, corresponding to the streamwise convection of droplets in the air and bubbles in the water, respectively. When the streamwise component of TMF is averaged in the mixed-phase region, its value is negative, indicating that the convection of bubbles dominates in the streamwise direction. As the Froude number increases, the magnitude of the streamwise component of TMF increases, corresponding to enhanced downstream convections of bubbles. The positive and negative vertical components of TMF occur alternately along the streamwise direction. Due to air entrainment, the vertical component of TMF is positive near the plunging points. In the downstream, bubbles are detained, causing a negative vertical component of TMF. As the Froude number increases, the air entrainment is enhanced. This is characterized by the expansion of the streamwise region with positive vertical TMF. Meanwhile, the downstream air detrainment is also enhanced, as reflected by the magnitude of the increment in the negative peak of the vertical TMF. In a further analysis of the TMF transport equation, it is discovered that the production term corresponding to the density gradient shows high consistency with the TKE in terms of both the streamwise variation and the response to the change in the Froude number. Based on this finding, a model of this production term is proposed. The a priori test shows a satisfactory correlation between the modelled value and LES data.
To finalize this paper, we compare the main findings of the present study on the mixed-phase turbulence generated by the plunging jet with the results of Hendrickson & Yue (Reference Hendrickson and Yue2019) on the mixed-phase turbulence in a wake flow. The purpose of providing this comparison is to evaluate whether the main conclusions of the present study are potentially common for different types of mixed-phase turbulence, or they are special in plunging jets. In the present study, it is discovered that the TKE transport is dominated by the balance between the production and dissipation terms. This is similar to the wake flow. However, the TKE generated by the TMF correlation term shows less significance in the plunging jet than in the wake flow. As noted by Hendrickson & Yue (Reference Hendrickson and Yue2019), the model for TMF is needed to close the TKE transport equation in the mixed-phase wake turbulence. In terms of the closure model of TMF, the results of both the present study and Hendrickson & Yue (Reference Hendrickson and Yue2019) show that the correlation between TMF and TKE is weak. We further investigated the transport equation of TMF, and reached a new finding that the production term of TMF corresponding to the density gradient is in good agreement with TKE, leading to a closure model for the production term of TMF, which is potentially useful for the future development of a dynamic model of TMF.
Supplementary movies
Supplementary movies for cases 3, 6 and 9 are available at https://doi.org/10.1017/jfm.2023.1081.
Funding
This research is supported by the NSFC Basic Science Centre Program for ‘Multiscale problems in nonlinear mechanics’ (no. 11988102), and NSFC projects (nos 11972038, 12272357).
Declaration of interests
The authors report no conflict of interest.
Appendix. Mesh-convergence test
In order to ensure that the results reported in this paper are independent of grid resolution, case 6 (see table 1) is simulated on four different grids. The information about grids used for conducting the mesh-convergence test is summarized in table 4. As shown, the grid resolution is gradually refined from grid 1 to grid 4.
Table 4. Information about grids used for resolution test.
Figure 19 compares the results obtained from different grid resolutions. Figures 19(a–f) show the results for the area of mixed-phase region 
$S^R$, mean streamwise velocity 
$\langle \bar {u}^R \rangle _{yz}$, mean vertical velocity 
$\langle \bar {v}^R \rangle _{yz}$, TKE 
$\langle k^R \rangle _{yz}$, streamwise TMF 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$, and vertical TMF 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$, respectively. It is observed from figures 19(a–d) that the area of mixed-phase region, mean velocity and TKE obtained from grids 2, 3 and 4 collapses. In comparison, the TMF requires a finer grid resolution to converge. As shown in figures 19(e,f), the magnitude variations of 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$ and 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$ are respectively 
$25\,\%$ and 
$20\,\%$ between grids 2 and 3, while they reduce to 
$5\,\%$ and 
$4\,\%$ as the resolution is further refined from grid 3 to grid 4. This comparison indicates that grid 3 is needed to obtain accurate results for 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$ and 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$. Therefore, grid 3 is used for the simulations of other cases in the present study.
Figure 19. Streamwise variation of (a) area of mixed-phase region 
$S^R$, (b) mean streamwise velocity 
$\langle \bar {u}^R \rangle _{yz}$, (c) mean vertical velocity 
$\langle \bar {v}^R \rangle _{yz}$, (d) TKE 
$\langle \langle k^R \rangle _{yz}$, (e) streamwise TMF 
$\langle \overline {\rho u^\prime }^R \rangle _{yz}$, and (f) vertical TMF 
$\langle \overline {\rho v^\prime }^R \rangle _{yz}$, obtained from different grids. Information about the different grids is given in table 4.
