1. Introduction
Flows involving complex, dynamic free-surface motion are found widely in industry and nature, with fuel sloshing in aircraft wings and wave impacts on coastal and offshore structures being prime examples. For waves in particular, the violent motion of the free surface often results in the entrainment of bubbles at the free surface, which can have significant effects on the overall dynamics and peak loads, and plays a major role in the exchange of gas between the ocean and atmosphere.
With a greater understanding of the effect of bubbles in breaking waves as our motivation, we seek improved approaches to their numerical simulation. The topic of bubble entrainment in breaking waves has been the subject of considerable experimental (e.g. Rapp, Melville & Longuet-Higgins Reference Rapp, Melville and Longuet-Higgins1990; Deane & Stokes Reference Deane and Stokes2002) and numerical (Chen et al. Reference Chen, Kharif, Zaleski and Li1999; Ma, Shi & Kirby Reference Ma, Shi and Kirby2011; Derakhti & Kirby Reference Derakhti and Kirby2014; Deike, Popinet & Melville Reference Deike, Popinet and Melville2015; Deike, Melville & Popinet Reference Deike, Melville and Popinet2016) research, and with advances in computational resources and mesh adaptivity in recent years, researchers have begun to conduct multi-phase simulations with the aim of resolving even the smallest bubble and droplet scales (Deike et al. Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016; Mostert, Popinet & Deike Reference Mostert, Popinet and Deike2022). This work has elucidated fundamental aspects of the process of wave breaking, including the degree of three-dimensionality in the flow (Mostert et al. Reference Mostert, Popinet and Deike2022), non-locality in the bubble breakup cascade (Chan, Lozano-Durán & Moin Reference Chan, Lozano-Durán and Moin2020; Chan et al. Reference Chan, Johnson, Moin and Urzay2021b), and the underlying physical mechanisms controlling bubble breakup (Rivière et al. Reference Rivière, Mostert, Perrard and Deike2021; Ruth et al. Reference Ruth, Aiyer, Rivière, Perrard and Deike2022). Whilst high-fidelity simulations are desirable for obtaining fundamental insight, their applicability in more industrially relevant settings is limited due to computational costs. Even with adaptive mesh refinement, the computational cost of multi-phase simulations as in Mostert et al. (Reference Mostert, Popinet and Deike2022) is significant: of the order of a month, and half a million CPU hours across several hundred cores.
An alternative numerical approach is to model the presence of bubbles (rather than resolving individual bubbles), with some form of population balance equation and, typically, an assumption that the bubbles are spherical. With this approach, there are two options for the treatment of the dispersed phase. First, it may be modelled as a continuum, through a bubble volume fraction, or number density field, which is subject to an evolution equation. Here, the evolution equations for the dispersed phase are partial differential equations, which are integrated in the same numerical framework as the equations of motion of the continuum liquid phase. In the mesh-based literature, such schemes are referred to as Eulerian–Eulerian, and have been developed for the simulation of bubble plume dynamics (Becker, Sokolichin & Eigenberger Reference Becker, Sokolichin and Eigenberger1994; Sokolichin & Eigenberger Reference Sokolichin and Eigenberger1994; Pfleger & Becker Reference Pfleger and Becker2001), air entrainment in breaking waves (Ma et al. Reference Ma, Shi and Kirby2011; Derakhti & Kirby Reference Derakhti and Kirby2014), and liquid jet breakup (Edelbauer Reference Edelbauer2017). Since the dynamics of bubbles dispersed in a liquid has a strong dependence on the bubble size, the treatment of polydisperse bubble distributions here is problematic. Typically, a population of bubbles is segregated into groups of similar sizes, each group requiring two additional evolution equations. Bubble breakup and coalescence may then be incorporated by source and sink terms exchanging mass (or number density) between bubble groups. An additional limitation is that this approach does not allow the tracking of individual bubbles over their lifetimes, but only statistical averages, and this constrains potential additional physical models, such as those for bubble break up. Although numerically this approach allows bubbles larger than the discretisation length scale of the continuous liquid phase (but with small number density), the assumptions used in the derivation of such models limit bubble sizes to smaller than the resolution of the continuous phase (Lakehal, Smith & Milelli Reference Lakehal, Smith and Milelli2002).
The second option for treating the dispersed phase is to model the bubbles individually, representing each bubble as a Lagrangian particle. Again, such approaches have been developed widely for mesh-based schemes, for example the work by Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016) focusing on bubble plumes, or the studies on turbulence–bubble interactions (Mazzitelli, Lohse & Toschi Reference Mazzitelli, Lohse and Toschi2003a,Reference Mazzitelli, Lohse and Toschib; Pozorski & Apte Reference Pozorski and Apte2009; Breuer & Hoppe Reference Breuer and Hoppe2017; Olsen, Skjetne & Johansen Reference Olsen, Skjetne and Johansen2017), and cavitation bubble clouds (Fuster & Colonius Reference Fuster and Colonius2011; Maeda & Colonius Reference Maeda and Colonius2018). In the mesh-based community, these methods are described as Eulerian–Lagrangian schemes. Each individual bubble interacts with the continuous phase through exchanges of momentum (and sometimes volume, as in Finn, Shams & Apte (Reference Finn, Shams and Apte2011), although in many cases where the concentration of the dispersed phase is small, volume exchanges are neglected). Schemes of this type allow for the tracking of individual bubbles, and a continuous polydisperse bubble distribution poses no additional challenge. However, the resolution of the liquid phase imposes an upper limit on the maximum bubble size that can be represented (Fraga et al. Reference Fraga, Stoesser, Lai and Socolofsky2016), and for very small bubbles, the computational costs increase with (a) the increasing number of bubbles, and (b) the increasing stiffness of the equation of motion of small bubbles due to the closure model for the drag force.
Temporarily setting aside the presence of bubbles, mesh-free methods have shown significant promise for simulations of breaking waves in recent decades. Smoothed particle hydrodynamics (SPH) is one mesh-free method, developed originally for astrophysical simulations (Gingold & Monaghan Reference Gingold and Monaghan1977; Lucy Reference Lucy1977), and since applied with considerable success to a range of terrestrial flows, including those with dynamically evolving free surfaces (Monaghan Reference Monaghan2012). The fluid is discretised by a set of Lagrangian particles, and spatial derivatives are approximated by weighted sums of fluid properties at neighbouring particles. Whilst tracking a deforming surface undergoing topological changes is a complex task in mesh-based methods, for SPH, little additional effort is required. There are now a wide variety of SPH schemes and related methods capable of simulating breaking waves, including weakly compressible SPH (Domínguez et al. Reference Domínguez2022), incompressible schemes (Lind et al. Reference Lind, Xu, Stansby and Rogers2012; Chow et al. Reference Chow, Rogers, Lind and Stansby2018; Guo et al. Reference Guo, Rogers, Lind and Stansby2018), and the moving particle semi-implicit method (Khayyer & Gotoh Reference Khayyer and Gotoh2009). Although multi-phase SPH schemes are well established (e.g. Hammani et al. Reference Hammani, Marrone, Colagrossi, Oger and Le Touzé2020), and capable of simulating multiple bubbles (Zhang, Sun & Ming Reference Zhang, Sun and Ming2015) and bubble-free-surface interactions (Sun et al. Reference Sun, Le Touzé, Oger and Zhang2021), we seek to avoid the cost of resolving both phases explicitly. We observe that the terms ‘Eulerian–Eulerian’ and ‘Eulerian–Lagrangian’ used above are misnomers (and somewhat ambiguous) in the context of SPH-based methods, while appropriate for Eulerian mesh-based numerical methods. In this work, we refer to the two approaches as ‘continuous–continuous’ and ‘continuous–discrete’, descriptions that remain clear even when used to describe schemes with non-Eulerian methods for the continuous phase. Of the bubble modelling approaches described above, developments in SPH lag behind mesh-based methods. A model based on the continuous–continuous approach has only recently been introduced to SPH (Fonty et al. Reference Fonty, Ferrand, Leroy, Joly and Violeau2019), but with very promising results for air entrainment in flow over a spillway (Fonty et al. Reference Fonty, Ferrand, Leroy and Violeau2020), though the method is currently limited to simplified closure models for interphase momentum exchange. We are not aware of any continuous–discrete SPH models for bubbly flows, although the SPH implementations closest in philosophy to this approach are perhaps the multi-phase dusty gas formulations, developed originally by Monaghan & Kocharyan (Reference Monaghan and Kocharyan1995) and Maddison & Monaghan (Reference Maddison and Monaghan1996), and extended more recently by Laibe & Price (Reference Laibe and Price2012).
Herein, we present an SPH implementation of the continuous–discrete approach for bubbly free-surface flows. The liquid is resolved via large eddy simulations (LES) using a semi-implicit isothermally compressible SPH framework, whilst each bubble is represented as a discrete Lagrangian particle that interacts with the liquid via exchanges of momentum, volume and sub-resolution turbulence closures. We focus particularly on integral models for bubble entrainment, breakup and free-surface interaction, with application to breaking waves.
The remainder of the paper is set out as follows. In § 2, we introduce the governing equations of our model. Section 3 presents details of the numerical implementation. In § 4, we test our simulation framework against numerical and experimental data for bubble plumes, and in § 5, we use our model to simulate the air entrainment in breaking waves. Section 6 is a summary of our conclusions. Further validation of our LES model is provided in the Appendix. A table of all symbols used in the work is given in table 1.
List of notation used herein. Except where stated explicitly, all properties are dimensionless.
2. Governing equations
The system considered is a continuous liquid phase, containing a dispersed bubble phase, as illustrated in figure 1. The liquid phase is governed by the isothermal compressible Navier–Stokes equations, whilst each bubble is modelled as a discrete particle that obeys Newton's second law. In the following, the subscripts 
$l$ and 
$b$ indicate properties in the liquid (continuous) and dispersed bubble phases, respectively. Where the subscript 
$l,b$ appears, it denotes a liquid property evaluated at a bubble. We consider the liquid to be isothermally compressible, with sound speed 
$c$, density 
$\rho _{l}$ and viscosity 
$\mu _{l}$. The bubbles are assumed to be spherical, comprised of gas with density 
$\rho _{b}$, and with a constant liquid–bubble surface tension 
$\gamma$. The entire system is subject to a gravitational acceleration 
$g\boldsymbol {e}_{g}$, where 
$\boldsymbol {e}_{g}$ is a unit vector. Following non-dimensionalisation by suitable integral length and velocity scales, the problem is parametrised by the Reynolds number 
$Re$, the Mach number 
$Ma$, the Froude number 
$Fr$, the (integral scale) Weber number 
$We$, and the density ratio 
$\beta =\rho _{l}/\rho _{b}$. Although our numerical framework is able to capture acoustic signals, in the present work these are not of interest, and are damped out by our implicit treatment of the pressure, hence 
$Ma$ is treated as a numerical, rather than physical, parameter. In the limit 
$Ma=0$ (
$c\to \infty$), our framework collapses to an incompressible framework. This approach of permitting weak compressibility, with an artificial sound speed, is common in SPH (fully explicit weakly compressible SPH being the most widely used variant in engineering simulations), but a difference in our approach is to use an implicit treatment, as in Khayyer & Gotoh (Reference Khayyer and Gotoh2009), permitting larger time steps, and resulting in more accurate pressure fields.
An illustration of the configuration considered herein. The liquid is treated as a continuum (blue), with dispersed bubbles (red) treated as discrete particles.
2.1. Liquid phase
The liquid phase is modelled as an isothermally compressible continuum, with an LES scheme. Due to the assumption of isothermal flow, we can write
\begin{equation} \frac{{\rm d}\tilde{p}}{{\rm d}\rho_{l}}=c^{2}, \end{equation}
where 
$\tilde {p}$ is the (implicitly) filtered pressure, and 
$c$ is an artificial speed of sound. The filtered continuity equation may be then expressed, in a Lagrangian frame of reference, as an evolution equation for the pressure:
\begin{equation} \frac{{\rm d}}{{\rm d}{t}}(\alpha\tilde{p}) ={-}\frac{\alpha}{Ma^{2}}\,\boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}_{l},\end{equation}
where 
$\alpha$ is the liquid volume fraction, and 
$\tilde {\boldsymbol {u}}_{l}$ is the implicitly filtered liquid velocity. Again we note that in the limit 
$Ma=0$, (2.2) becomes 
$\boldsymbol {\nabla }\boldsymbol {\cdot }\tilde {\boldsymbol {u}}_{l}=0$, and incompressibility is recovered. We next assume that the liquid is weakly compressible, the compressibility 
$1/\rho _{l}c^{2}\ll 1$ allowing us to neglect terms involving the density variation in the (filtered) momentum equation, which is written as
\begin{equation} \frac{{\rm d}}{{\rm d}{t}}(\alpha\tilde{\boldsymbol{u}}_{l}) ={-}\boldsymbol{\nabla}(\alpha\tilde{p}) +\frac{\alpha}{Fr^{2}}\,\boldsymbol{e}_{g} +\frac{1}{Re}\,\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\alpha(1+\nu_{srs})\,\boldsymbol{\nabla}\tilde{\boldsymbol{u}}_{l}\right) +\boldsymbol{M}, \end{equation}
where 
$\nu _{srs}$ is a dimensionless sub-resolution viscosity, determined by the LES closure model, and the term 
$\boldsymbol {M}$ represents the momentum exchange between the liquid and bubble phases. In practice, we solve (2.2) and (2.3) in a frame that deviates from perfectly Lagrangian by a small velocity 
$\boldsymbol {u}_{ps}$, referred to in the SPH literature as a shifting velocity, and introduced to add stability to the numerical solution. This results in a small error that scales with the resolution of the discretisation scheme, associated with the advection terms 
$\boldsymbol {u}_{ps}\boldsymbol {\cdot }\boldsymbol {\nabla }$ that are omitted from (2.2) and (2.3). This approach is used widely in the SPH literature (Lind et al. Reference Lind, Xu, Stansby and Rogers2012).
2.1.1. The LES model
We denote the implicit filter width of the SPH framework as 
$\tilde {\varDelta }$. The filtered equations are closed with a model for the sub-resolution viscosity, which is comprised of a shear-induced and bubble-induced eddy viscosity component: 
$\nu _{srs}=\nu _{S}+\nu _{B}$. The bubble-induced viscosity 
$\nu _{B}$ accounts for the production of sub-resolution turbulence by bubbles following the model of Sato & Sekoguchi (Reference Sato and Sekoguchi1975), and is described in § 2.2.2. For the shear-induced turbulence, we use the mixed-scale model of Lubin et al. (Reference Lubin, Vincent, Abadie and Caltagirone2006), with
\begin{equation} \nu_{S}=Re\,{C}_{M}\,\tilde{\varDelta}^{1+\xi}\,\lvert\tilde{\mathcal{S}}\rvert^{\xi/2}(q_{c}^{2})^{(1-\xi)/2}, \end{equation}
where 
$\tilde {\mathcal {S}}=(\boldsymbol {\nabla }\tilde {\boldsymbol {u}}_{l}+(\boldsymbol {\nabla }\tilde {\boldsymbol {u}}_{l})^{\rm T})/2$ is the resolved strain rate tensor, 
$q_{c}^{2}$ is the test-filtered kinetic energy, 
$\xi =0.5$ and 
${C}_{M}=0.06$. We evaluate 
$q_{c}$ by explicitly filtering 
$\tilde {\boldsymbol {u}}_{l}$ with a Shepard filter (see § 3.1.2). We note that multi-phase LES closure models are still an open area of research. The separation of the effective viscosity into shear- and bubble-induced components, following Derakhti & Kirby (Reference Derakhti and Kirby2014), is a modelling simplification that presumes that the effects of bubbles and free surfaces on (2.4) may be neglected. As in Derakhti & Kirby (Reference Derakhti and Kirby2014), the turbulent dissipation rate is proportional to the shear-induced viscosity and the square of the strain rate tensor norm, and is calculated as
\begin{equation} \varepsilon=\frac{1+\nu_{S}}{Re}\,\lvert\tilde{\mathcal{S}}\rvert^{2}. \end{equation}In the development of our method, we explored several additional closure models, including a standard Smagorinsky model and the dynamic Germano model (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Lilly Reference Lilly1992), with both local (via Shepard filtering) and Lagrangian (along streamlines) averaging. We choose the mixed-scale model for our simulations because it is known to yield good results in flows with deforming free surfaces (Lubin et al. Reference Lubin, Vincent, Abadie and Caltagirone2006), and in tests of the decay of single-phase isotropic turbulence (included in the Appendix), it yielded improved results in our numerical framework compared with the other closure models.
2.2. Dispersed bubble phase
The dispersed phase is represented by a discrete set 
$\mathcal {B}$ of 
$N_{b}$ Lagrangian bubbles, each with velocity 
$\boldsymbol {u}_{b}$, position 
$\boldsymbol {r}_{b}$, radius 
$a_{b}$, and volume 
$V_{b}=4{\rm \pi} {a}_{b}^{3}/3$. The system of bubbles is governed by
$$\begin{gather} \frac{{\rm d}\boldsymbol{r}_{b}}{{\rm d}t}=\boldsymbol{u}_{b}, \end{gather}$$
$$\begin{gather}V_{b}\,\frac{{\rm d}\boldsymbol{u}_{b}}{{\rm d}t}=\boldsymbol{F}_{d}+\boldsymbol{F}_{l}+\boldsymbol{F}_{vm}+\boldsymbol{F}_{g}, \end{gather}$$
for each bubble in 
$\mathcal {B}$. Here, 
$\boldsymbol {F}_{d}$, 
$\boldsymbol {F}_{l}$, 
$\boldsymbol {F}_{vm}$ and 
$\boldsymbol {F}_{g}$ are the drag, lift, virtual mass and buoyancy forces acting on the bubble due to the surrounding liquid. These forces are evaluated for each bubble 
$i\in \mathcal {B}$ through closure models as in e.g. Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016) and Derakhti & Kirby (Reference Derakhti and Kirby2014), with
$$\begin{gather} \boldsymbol{F}_{d,i}=\tfrac{1}{2}C_{d}\beta{\rm \pi}{a}_{b,i}^{2}\left\lvert\boldsymbol{u}_{rel,i}\right\rvert\boldsymbol{u}_{rel,i}, \end{gather}$$
$$\begin{gather}\boldsymbol{F}_{l,i}=C_{l}\beta{V}_{b,i}\boldsymbol{u}_{rel,i}\times\boldsymbol{\nabla}\times\boldsymbol{u}_{l,i}, \end{gather}$$
$$\begin{gather}\boldsymbol{F}_{vm,i}=C_{vm}\beta{V}_{b,i}\left(\frac{{\rm d}\boldsymbol{u}_{l}}{{\rm d}t}-\frac{{\rm d}\boldsymbol{u}_{b,i}}{{\rm d}t}\right), \end{gather}$$
$$\begin{gather}\boldsymbol{F}_{g,i}=\frac{(1-\beta)V_{b,i}}{Fr^{2}}\,\boldsymbol{e}_{g}, \end{gather}$$
where 
$\beta$ is the density ratio, the relative velocity between the bubble and liquid phase at bubble 
$i$ is given by 
$\boldsymbol {u}_{rel,i}=\boldsymbol {u}_{l}(\boldsymbol {r}_{b,i})-\boldsymbol {u}_{b,i}$, and 
$C_{d}$, 
$C_{l}$ and 
$C_{vm}$ are drag, lift and virtual mass coefficients, respectively. Note the absence of the tilde in the liquid velocity appearing in the definition of relative velocity, and in (2.7b) and (2.7c). Following Breuer & Hoppe (Reference Breuer and Hoppe2017), the sub-resolution fluctuating part of the liquid velocity is modelled stochastically, with
\begin{equation} \boldsymbol{u}_{l}(\boldsymbol{r}_{b,i}) =\tilde{\boldsymbol{u}}_{l}(\boldsymbol{r}_{b,i}) +\boldsymbol{u}^{\prime}_{l}(\boldsymbol{r}_{b,i}) \end{equation}
comprising the LES filtered velocity and a sub-resolution fluctuating part. The calculation of 
$\boldsymbol {u}^{\prime }_{l}(\boldsymbol {r}_{b,i})$ is described in § 2.2.1. Following Derakhti & Kirby (Reference Derakhti and Kirby2014), the drag coefficient in (2.7a) is modelled using the standard drag curve of Clift, Grace & Weber (Reference Clift, Grace and Weber1978) as
\begin{equation} {C}_{d}= \begin{cases} 0.44, & {Re}_{b,i}>1000,\\ \dfrac{24}{Re_{b,i}}(1+0.15\,{Re}_{b,i}^{0.687}), & Re_{b,i}\le1000, \end{cases} \end{equation}where the relative bubble Reynolds number is related to the integral-scale Reynolds number by
\begin{equation} Re_{b,i}=2a_{b,i}\left\lvert\boldsymbol{u}_{rel,i}\right\rvert{Re}. \end{equation}
Following Derakhti & Kirby (Reference Derakhti and Kirby2014), the coefficient of virtual mass and the lift coefficient are set to 
$C_{vm}=C_{l}=0.5$. The momentum exchange from bubble 
$i$ back to the liquid phase is
\begin{equation} \boldsymbol{M}_{b,i}=(-\boldsymbol{F}_{d,i}-\boldsymbol{F}_{l,i}-\boldsymbol{F}_{vm,i})/\beta, \end{equation}
with 
$\beta$ the density ratio as defined at the start of the section.
2.2.1. Langevin model
The fluctuating velocity component ‘felt’ by the bubbles is calculated through the integration of a stochastic model following work by Pozorski & Apte (Reference Pozorski and Apte2009) and Breuer & Hoppe (Reference Breuer and Hoppe2017). The fluctuation velocity obeys the Langevin equation
\begin{equation} {\rm d}\boldsymbol{u}_{l}^{\prime}={-}{\boldsymbol{\mathsf{G}}}\boldsymbol{u}_{l}^{\prime}\,{\rm d}t +\sqrt{2\sigma_{srs}^{2}}\,{\boldsymbol{\mathsf{B}}}\boldsymbol{dW}, \end{equation}
where 
$\boldsymbol {dW}$ is a Wiener process, and 
$\boldsymbol {\mathsf {G}}$ and 
$\boldsymbol {\mathsf {B}}$ are drift and diffusion matrices, respectively. The quantity 
$\sigma _{srs}$ is the standard deviation of the fluctuating velocity, related to the sub-resolution turbulent kinetic energy 
$k_{srs}$ by
\begin{equation} \sigma_{srs}=\sqrt{\tfrac{2}{3}\,k_{srs}}. \end{equation}The sub-resolution turbulent kinetic energy is estimated from the double filtered velocity as
\begin{equation} k_{srs}=\tfrac{1}{2}\left\lvert\tilde{\boldsymbol{u}}_{l}-\boldsymbol{\hat{\tilde{u}}}_{l}\right\rvert^{2}, \end{equation}
where the 
$\hat {\ }$ indicates explicit filtering with a test filter, described in § 3.1.2.
Following Breuer & Hoppe (Reference Breuer and Hoppe2017), by taking advantage of the fact that a Langevin equation may be integrated analytically, we transform (2.12) into the recursion equation
\begin{equation} \boldsymbol{u}_{l}^{\prime}\left(\boldsymbol{r}_{b,i},t+\delta{t}\right)={\boldsymbol{\mathsf{E}}}\,\boldsymbol{u}_{l}^{\prime}\left(\boldsymbol{r}_{b,i},t\right)+{\boldsymbol{\mathsf{W}}}\boldsymbol{\zeta},\end{equation}
where 
$\delta {t}$ is a time increment, 
$\boldsymbol {\mathsf {E}}$ is the exponential of the drift matrix 
$\boldsymbol {\mathsf {G}}$, 
$\boldsymbol {\mathsf {W}}$ is the square root of the velocity fluctuation covariance matrix, and 
$\boldsymbol {\zeta }$ is a random vector whose components are distributed normally. Denoting the filtered relative velocity (excluding the fluctuations) as 
$\tilde {\boldsymbol {u}}_{rel}$, we define the matrix
\begin{equation} {\boldsymbol{\mathsf{R}}}=\frac{1}{\left\lvert\tilde{\boldsymbol{u}}_{rel}\right\rvert^{2}}\,\tilde{\boldsymbol{u}}_{rel}\otimes\tilde{\boldsymbol{u}}_{rel}, \end{equation}and the exponential of the drift matrix is then given by
\begin{equation} {\boldsymbol{\mathsf{E}}}=\left(E_{{\parallel}}-E_{{\perp}}\right){\boldsymbol{\mathsf{R}}}+E_{{\perp}}{\boldsymbol{\mathsf{I}}}, \end{equation}
where 
$\boldsymbol {\mathsf {I}}$ is the identity matrix, and
$$\begin{gather} E_{{\parallel}}=\exp\left(\frac{-\delta{t}}{\tau_{{\parallel}}}\right), \end{gather}$$
$$\begin{gather}E_{{\perp}}=\exp\left(\frac{-\delta{t}}{\tau_{{\perp}}}\right). \end{gather}$$Here,
$$\begin{gather} \tau_{{\parallel}}=\tau_{srs}\left(1+\frac{\left\lvert\tilde{\boldsymbol{u}}_{rel}\right\rvert}{\sigma_{srs}^{2}}\right)^{-{1}/{2}}, \end{gather}$$
$$\begin{gather}\tau_{{\perp}}=\tau_{srs}\left(1+4\,\frac{\left\lvert\tilde{\boldsymbol{u}}_{rel}\right\rvert}{\sigma_{srs}^{2}}\right)^{-{1}/{2}} \end{gather}$$
are the sub-resolution time scales associated with fluctuations parallel and perpendicular to the filtered velocity, and the sub-resolution time scale 
$\tau _{srs}$ is related to the velocity fluctuations by
\begin{equation} \tau_{srs}=\frac{C\tilde{\varDelta}}{\sigma_{srs}}, \end{equation}
with the constant 
$C=1$. The factors relating 
$\tau _{\parallel }$ and 
$\tau _{\perp }$ to the sub-resolution time scale 
$\tau _{srs}$ account for the crossing trajectory and continuity effects (Pozorski & Apte Reference Pozorski and Apte2009). The square root of the covariance matrix is given by
\begin{equation} {\boldsymbol{\mathsf{W}}}=\left(W_{{\parallel}}-W_{{\perp}}\right){\boldsymbol{\mathsf{R}}}+W_{{\perp}}{\boldsymbol{\mathsf{I}}}, \end{equation}where
$$\begin{gather} W_{{\parallel}}=\sigma_{srs}\sqrt{1-\exp\left(-\frac{2\,\delta{t}}{\tau_{{\parallel}}}\right)}, \end{gather}$$
$$\begin{gather}W_{{\perp}}=\sigma_{srs}\sqrt{1-\exp\left(-\frac{2\,\delta{t}}{\tau_{{\perp}}}\right)}. \end{gather}$$2.2.2. Bubble-induced turbulence model
As mentioned above, we evaluate the bubble-induced turbulence, following Derakhti & Kirby (Reference Derakhti and Kirby2014), based on the model of Sato & Sekoguchi (Reference Sato and Sekoguchi1975), where the contribution of an individual bubble to the turbulent viscosity is proportional to the product of the bubble diameter and the relative velocity. In our discrete bubble framework, the contribution of an individual bubble 
$j$ is
\begin{equation} \nu_{B,b,j}=C_{\nu,B}2a_{b,j}\left\lvert\boldsymbol{u}_{rel,j}\right\rvert, \end{equation}
where the constant 
$C_{\nu,B}=0.6$ as in Derakhti & Kirby (Reference Derakhti and Kirby2014). To obtain the bubble-induced turbulent viscosity in the liquid 
$\nu _{B}$, we interpolate 
$\nu _{B,b}$ from the bubbles to the liquid phase, as described in § 3.1.1.
2.2.3. Bubble entrainment and breakup
Our intention is to simulate flows where bubbles are entrained at the free surface. We use an entrainment model, similar in principle to that of Ma et al. (Reference Ma, Shi and Kirby2011) and Derakhti & Kirby (Reference Derakhti and Kirby2014), in which a fraction of the turbulent kinetic energy of the liquid is assumed to be converted into surface energy as bubbles are created (or entrained) at the free surface. We further include a model for bubble breakup, which is based on the imbalance between the restoring pressure on a bubble due to surface tension, and the deforming stress due to the turbulent motion of the liquid, based on the models of Martínez-Bazán, Montañes & Lasheras (Reference Martínez-Bazán, Montañes and Lasheras1999a,Reference Martínez-Bazán, Montañes and Lasherasb) and Martínez-Bazán et al. (Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montañes and Lasheras2010). These models cannot be explained clearly without reference to our discretisation scheme, and we defer detailed description of them to later, in § 3.4.
3. Numerical implementation
3.1. The SPH discretisation
The liquid phase is represented by set of discrete particles 
$\mathcal {P}$, each of which we label 
$i\in [1,N]$, where 
$N$ is the total number of particles. In the present work, we treat the liquid as weakly compressible, but this weak compressibility is treated implicitly, by solving an elliptic equation (Helmholtz, as opposed to Poisson) as in incompressible SPH frameworks. Hence in the present work, the core of the SPH scheme follows closely King & Lind (Reference King and Lind2021) and Lind et al. (Reference Lind, Xu, Stansby and Rogers2012). All particles carry a (constant) liquid mass 
$m$, and with the assumption of a weakly compressible liquid, the associated liquid volume 
$V_{l}$ is assumed constant. The position of particle 
$i$ is denoted 
$\boldsymbol {r}_{i}$, its smoothing length is 
$h_{i}$, and the total volume associated with the particle is 
$V_{i}$. We denote the difference in the property 
$(\,{\cdot }\,)$ of two particles 
$i$ and 
$j$ as 
$(\,{\cdot}\,)_{ij}=(\,{\cdot}\,)_{i}-(\,{\cdot}\,)_{j}=-(\,{\cdot}\,)_{ji}$. In SPH, values and derivatives of field variables at the location of particle 
$i$ are calculated using a weighted sum of the values of the field variables at the neighbouring particles 
$j\in \mathcal {P}_{i}$, where the weights are obtained from a kernel function 
$W(\lvert \boldsymbol {r}_{ij},h_{i}\rvert )=W_{ij}$ and its derivatives. Here, 
$\mathcal {P}_{i}$ is the set of neighbours of particle 
$i$, and contains all particles 
$j$ with 
$\lvert \boldsymbol {r}_{ij}\rvert \le {r}_{s,i}$, where 
$r_{s,i}$ is the support radius of the kernel of particle 
$i$. Throughout this work, we use the Wendland C2 kernel (Wendland Reference Wendland1995) for which the support radius is 
$r_{s,i}=2h_{i}$. We use initial particle spacing 
$\delta {r}=h_{0}/1.3$, where 
$h_{0}$ is the smoothing length when 
$\alpha =1$. In all cases, we set the implicit filter scale to the local smoothing length: 
$\tilde {\varDelta }_{i}=h_{i}$. For a derivation and analysis of SPH fundamentals, we refer the reader to Price (Reference Price2012), Fatehi & Manzari (Reference Fatehi and Manzari2011) and Monaghan (Reference Monaghan2012). In a perfectly Lagrangian framework, particles follow streamlines, which can result in highly anisotropic particle distributions, particularly around stagnation points, degrading the accuracy of the simulation. To regularise the particle distribution, we use the particle shifting technique of Lind et al. (Reference Lind, Xu, Stansby and Rogers2012) to set 
$\boldsymbol {u}_{ps}$ as in King & Lind (Reference King and Lind2021). The ability of the underlying SPH methodology to simulate wave propagation accurately has been demonstrated previously, for example by Lind et al. (Reference Lind, Xu, Stansby and Rogers2012) and Skillen et al. (Reference Skillen, Lind, Stansby and Rogers2013). Note that we do not include surface tension effects in the single-phase SPH simulation – rather, they appear only in the closure terms governing bubble dynamics.
3.1.1. Interpolation between phases
The interaction between the liquid and bubbles occurs through the modification of the liquid volume fraction 
$\alpha$ due to the presence of the bubbles, the momentum exchange 
$\boldsymbol {M}$ between the phases, and the bubble-induced turbulent viscosity 
$\nu _{B}$. To achieve this, it is necessary to interpolate between the phases, i.e. to calculate the value of a liquid property at a bubble location, and vice versa. Figure 2 provides a diagrammatic overview of the framework, including those properties that are exchanged between phases. Numerically, the interphase interpolation is as follows. The total volume associated with SPH particle 
$i$ is defined as
\begin{equation} V_{i}=V_{l,i}\Bigg(1+\displaystyle\sum_{j\in\mathcal{B}_{i}}W(\boldsymbol{r}_{b,j}-\boldsymbol{r}_{i},h_{i})V_{b,j}\Bigg), \end{equation}
where 
$\mathcal {B}_{i}$ is the set of all bubbles within the support radius of particle 
$i$. The summation in the right-hand side of (3.1) represents the contribution to the volume of SPH particle 
$i$ of all individual bubbles in 
$\mathcal {B}_{i}$. The liquid volume fraction of particle 
$i$ is then
\begin{equation} \alpha_{i}=\frac{V_{l,i}}{V_{i}} . \end{equation}As the volume of the bubbles is accounted for in the liquid phase, the SPH particles effectively expand in the proximity of gas bubbles. To retain an accurate SPH approximation, the smoothing length of each SPH particle must be adjusted accordingly, to maintain
\begin{equation} \displaystyle\sum_{j\in\mathcal{P}_{i}}W(\boldsymbol{r}_{j} -\boldsymbol{r}_{i},h_{i})V_{j}\approx1\quad\forall{i}\in\mathcal{P}. \end{equation}
With both 
$V_{i}$ (through (3.1)) and the partition of unity (3.3) having a nonlinear dependence on 
$h_{i}$, we cannot choose 
$h_{i}$ explicitly to satisfy (3.3). However, by setting
\begin{equation} h_{i}=h_{0}\left(\frac{V_{i}}{V_{l,i}}\right)^{{1}/{d}}, \end{equation}
in which 
$d=3$ is the number of spatial dimensions, we obtain a system in which the SPH particle distribution expands and contracts in response to the bubble volumes. The response is not instantaneous, but occurs over a finite time, as the volume effects propagate through the particle distribution. However, when averaged over time, this system results in a discretisation for which (3.3) is satisfied. This effect is discussed further in § 3.3.
A schematic diagram of the numerical framework, showing the liquid and bubble properties that are interpolated between phases.
The momentum exchange between the bubble and liquid phases is evaluated at each bubble (denoted 
$\boldsymbol {M}_{b}$), and then interpolated back to the liquid phase through
\begin{equation} \boldsymbol{M}_{i}=V_{i}\displaystyle\sum_{j\in\mathcal{B}_{i}} \boldsymbol{M}_{b,j}\,W(\boldsymbol{r}_{b,j}-\boldsymbol{r}_{i},h_{i}). \end{equation}
The bubble-induced turbulence at each bubble, 
$\nu _{B,b}$, is interpolated back to each SPH particle in the same manner to obtain 
$\nu _{B}$. Evaluation of the lift, drag and virtual mass forces for each bubble requires the knowledge of the filtered liquid velocity and liquid velocity gradients, and the turbulent kinetic energy, at each bubble location. Additionally, the bubble entrainment model described in § 3.4 requires the turbulent dissipation rate to be interpolated from the liquid to each bubble location. These properties are interpolated from SPH particles to bubble locations through
\begin{equation} \phi_{b,j}=\displaystyle\sum_{i\in{P}}\phi_{i}\,W(\boldsymbol{r}_{b,j}-\boldsymbol{r}_{i},h_{i})V_{i}\quad\forall{j\in\mathcal{B}},\end{equation}
where 
$\phi _{i}$ is the value at particle 
$i$ of the property to be interpolated, and 
$\phi _{b,j}$ is the interpolated property at bubble 
$j$. In our code, we construct an array for each particle 
$i$ containing the indices 
$\mathcal {P}_{i}$, and a global array containing the indices of the bubble neighbours of all particles: 
$[\mathcal {B}_{1}\ldots \mathcal {B}_{i}\ldots \mathcal {B}_{N}]$.
3.1.2. The SPH operators
For the test filter used to evaluate 
$k_{srs}$, and 
$q_{c}$ in the LES closure model, we use a normalised Shepard filter
\begin{equation} \hat{\tilde{\phi}}_{i}=\frac{\sum_{j\in\mathcal{P}_{i}}\tilde{\phi}_{j}W_{ij}V_{j}}{\sum_{j\in\mathcal{P}_{i}}W_{ij}V_{j}}.\end{equation}First derivatives are discretised according to
\begin{equation}
\langle\boldsymbol{\nabla}\tilde{\phi}\rangle_{i}=\displaystyle\sum_{j\in\mathcal{P}_{i}}(\tilde{\phi}_{j}-\tilde{\phi}_{i})
\boldsymbol{\nabla}{W}^{{\star}}_{ij}{V}_{j},\quad\langle\boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}\rangle_{i}=\displaystyle\sum_{j\in\mathcal{P}_{i}}
(\tilde{\boldsymbol{u}}_{j}-\tilde{\boldsymbol{u}}_{i})\boldsymbol{\cdot}\boldsymbol{\nabla}{W}^{{\star}}_{ij}{V}_{j},
\end{equation}
where the corrected kernel gradient 
$\boldsymbol {\nabla }{W}^{\star }_{ij}$ due to Bonet & Lok (Reference Bonet and Lok1999) is used, as detailed in King & Lind (Reference King and Lind2021). This provides first-order consistency for first derivatives. The angled brackets indicate that the quantity is an SPH approximation to the gradient. The Laplacian is approximated using the formulation of Morris, Fox & Zhu (Reference Morris, Fox and Zhu1997) as
\begin{equation} \langle{\nabla}^{2}\tilde{\phi}\rangle_{i}=\displaystyle\sum_{j\in\mathcal{P}_{i}}\frac{2\tilde{\phi}_{ij}}{\lvert\boldsymbol{r}_{ij}\rvert^{2}}\,\boldsymbol{r}_{ij}\boldsymbol{\cdot}\boldsymbol{\nabla}{W}_{ij}{V}_{j},\end{equation}
and for the inhomogeneous ‘div-grad’ operator with spatially varying coefficient 
$\kappa$, we use
\begin{equation} \langle\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\kappa\boldsymbol{\nabla}\phi\right)\rangle_{i}=\displaystyle\sum_{j\in\mathcal{P}_{i}}\frac{2\bar{\kappa}_{ji}\tilde{\phi}_{ij}}{\lvert\boldsymbol{r}_{ij}\rvert^{2}}\boldsymbol{r}_{ij}\boldsymbol{\cdot}\boldsymbol{\nabla}{W}_{ij}{V}_{j}, \end{equation}
with 
$\bar {\kappa }_{ji}={2\kappa _{i}\kappa _{j}}/(\kappa _{i}+\kappa _{j})$ the harmonic mean. To evaluate the shifting velocity 
$\boldsymbol {u}_{ps}$, we calculate the gradient of the particle number density as
\begin{equation} \boldsymbol{\nabla}{\rho}_{N,i}=\displaystyle\sum_{j\in\mathcal{P}_{i}}\Bigg[1+\frac{1}{4}\left(\frac{W_{ij}}{W_{ii}}\right)^{4}\Bigg]\boldsymbol{\nabla}{W}_{ij}V_{j}, \end{equation}and then set the shifting velocity as
\begin{equation} \boldsymbol{u}_{ps,i}=\frac{h_{i}^{2}}{4\,\delta{t}}\begin{cases}\boldsymbol{\nabla}{\rho}_{N,i} & \forall{i}\in\mathcal{P}_{I},\\ (\boldsymbol{n}_{i}\boldsymbol{\cdot}\boldsymbol{\nabla}{\rho}_{N,i}) \boldsymbol{n}_{i} & \forall{i}\in\mathcal{P}_{FS}, \end{cases} \end{equation}
where 
$\boldsymbol {n}_{i}$ is the unit vector normal to the surface at particle 
$i$, 
$\mathcal {P}_{I}$ is the set of internal particles, and 
$\mathcal {P}_{FS}$ is the set of free-surface particles. This treatment of the shifting velocity at the free surface, and the identification of free-surface particles, follows Lind et al. (Reference Lind, Xu, Stansby and Rogers2012) and King & Lind (Reference King and Lind2021). The surface normal vectors are evaluated as
\begin{equation} \boldsymbol{n}^{{\star}}_{i}=\frac{1}{h_{i}}\displaystyle\sum_{j\in\mathcal{P}_{i}}V_{j}\,\boldsymbol{\nabla}{W}_{ij},\end{equation}
where the term 
$1/h_{i}$ normalises the magnitude of the vector relative to the resolution. Following Chow et al. (Reference Chow, Rogers, Lind and Stansby2018), the surface-normal vectors are then smoothed using the normalised Shepard filter as in (3.7), with 
$\boldsymbol {n}_{i}=\widehat {\boldsymbol {n}_{i}^{\star }}$.
3.2. Fractional step approach for isothermally compressible liquid
Our approach to the time integration of (2.2) and (2.3) combines the methods used in King & Lind (Reference King and Lind2021), and the approach described in Khayyer & Gotoh (Reference Khayyer and Gotoh2009, Reference Khayyer and Gotoh2016). We use a fractional step algorithm to solve (2.2) and (2.3) in our SPH framework. Following the classic projection method of Chorin (Reference Chorin1968), initially introduced to SPH in Cummins & Rudman (Reference Cummins and Rudman1999), the right-hand side of (2.3) is split, with viscous and advective terms being applied in a predictor step, and the pressure gradient and any divergence-free body forces (e.g. gravity) being used in a projection step to obtain a velocity field that satisfies the continuity equation (2.2). Splitting (2.3) as described above, we obtain
$$\begin{gather} \alpha^{n+1}\tilde{\boldsymbol{u}}^{{\star}}_{l}= \alpha^{n}\tilde{\boldsymbol{u}}^{n}_{l}+\delta{t}\left[\frac{1}{Re}\,\boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha^{n}(1+\nu_{srs})\,\boldsymbol{\nabla}\tilde{\boldsymbol{u}}_{l}^{n})+\boldsymbol{M}\right], \end{gather}$$
$$\begin{gather}\alpha^{n+1}\tilde{\boldsymbol{u}}_{l}^{n+1}=\alpha^{n+1}\tilde{\boldsymbol{u}}_{l}^{{\star}}-\delta{t}\left[\boldsymbol{\nabla}(\alpha^{n+1}\tilde{p}^{n+1})-\frac{\alpha^{n+1}}{Fr^{2}}\,\boldsymbol{e}_{g}\right] , \end{gather}$$
where 
$\tilde {\boldsymbol {u}}^{\star }_{l}$ is an intermediate velocity, which is not required to be compatible with the continuity equation (2.2). We require the velocity field at the end of the time step 
$\tilde {\boldsymbol {u}}_{l}^{n+1}$ to satisfy (2.2), and we take the divergence of (3.14b), to get
\begin{align} \alpha^{n+1}\,\boldsymbol{\nabla}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}_{l}^{n+1}&={-}\tilde{\boldsymbol{u}}_{l}^{n+1}\boldsymbol{\cdot}\boldsymbol{\nabla}\alpha^{n+1} +\boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha^{n+1}\tilde{\boldsymbol{u}}_{l}^{{\star}})\nonumber\\ &\quad -\delta{t}\,{\nabla}^{2}(\alpha^{n+1}\tilde{p}^{n+1})+\frac{\delta{t}}{Fr^{2}}\,\boldsymbol{e}_{g}\boldsymbol{\cdot}\boldsymbol{\nabla}\alpha^{n+1}. \end{align}
Substituting the right-hand side of (3.15) into the final term of (2.2) evaluated at time step 
$n+1$, and replacing the time derivative with a backwards Euler difference equation, yields
\begin{align} \frac{\alpha^{n+1}\tilde{p}^{n+1}-\alpha^{n}\tilde{p}^{n}}{\delta{t}}&= \frac{1}{Ma^{2}}\left[\vphantom{\frac{1}{Ma^{2}}}\tilde{\boldsymbol{u}}_{l}^{n+1}\boldsymbol{\cdot}\boldsymbol{\nabla}\alpha^{n+1}-\boldsymbol{\nabla}\boldsymbol{\cdot}(\alpha^{n+1}\tilde{\boldsymbol{u}}_{l}^{{\star}})\right.\nonumber\\ &\left.\quad {}+\delta{t}\,{\nabla}^{2}(\alpha^{n+1}\tilde{p}^{n+1})-\frac{\delta{t}}{Fr^{2}}\,\boldsymbol{e}_{g}\boldsymbol{\cdot}\boldsymbol{\nabla}\alpha^{n+1}\right]. \end{align}
As (3.16) contains 
$\tilde {\boldsymbol {u}}_{l}^{n+1}$ explicitly, we substitute (3.14b) back into (3.16), note that 
$(1/\alpha )\,\boldsymbol {\nabla }\alpha =\boldsymbol {\nabla }\ln \alpha$, and collect terms containing 
$\tilde {p}^{n+1}$, obtaining
$$\begin{gather}
{\nabla}^{2}(\alpha^{n+1}\tilde{p}^{n+1})-\frac{Ma^{2}\,\alpha^{n+1}\tilde{p}^{n+1}}{\delta{t}^{2}}-\boldsymbol{\nabla}(\ln\alpha^{n+1})\boldsymbol{\cdot}\boldsymbol{\nabla}(\alpha^{n+1}\tilde{p}^{n+1})\nonumber\\
=\frac{\alpha^{n+1}}{\delta{t}}\,\boldsymbol{\nabla}\boldsymbol{\cdot}(\tilde{\boldsymbol{u}}_{l}^{{\star}})-\frac{Ma^{2}\,\alpha^{n}\tilde{p}^{n}}{\delta{t}^{2}}.
\end{gather}$$
Note that in the combined limits of single-phase (
$\alpha =1$) incompressible (
$Ma=0$) flow, the standard Poisson equation is recovered.
3.3. Temporal evolution of both phases
We now describe the complete algorithm for the temporal evolution of the complete system. In the following, each operation is applied to every SPH particle 
$i\in \mathcal {P}$ or bubble 
$i\in \mathcal {B}$, and we have dropped the subscripts 
$i$ for clarity. We introduce the superscripts 
$n$, 
$\star$ and 
$n+1$ to represent properties at the current, intermediate and next time steps. The algorithm is as follows.
(i) Advect particles to intermediate positions according to
$\boldsymbol {r}^{\star }=\boldsymbol {r}^{n}+\delta {t}\,\tilde {\boldsymbol {u}}_{l}$, and bubbles to new positions according to $\boldsymbol {r}_{b}^{n+1}=\boldsymbol {r}_{b}^{n}+\delta {t}\,\boldsymbol {u}_{b}^{n}$.(ii) Construct boundary conditions via mirror particles, calculate bubble entrainment and breakup if any (as described in § 3.4), and build neighbour lists
$\mathcal {P}_{i}$ and $\mathcal {B}_{i}$.(iii) Calculate the shifting velocity
$\boldsymbol {u}_{ps}$ from (3.12), based on a modified form of the Fickian shifting introduced by Lind et al. (Reference Lind, Xu, Stansby and Rogers2012), described in detail in King & Lind (Reference King and Lind2021).(iv) Interpolate bubble volumes into SPH particle locations using (3.1) and (3.2) to obtain volume fractions
$\alpha ^{n+1}$. Adjust smoothing lengths via (3.4). Note that this is a first-order approximation to $\alpha ^{n+1}$, but it allows us to retain an explicit scheme for the volume fractions.(v) Evaluate
$\boldsymbol {\nabla }\tilde {\boldsymbol {u}}_{l}$, $\varepsilon$ and $k_{srs}$, and interpolate to bubble positions through (3.6).(vi) Evolve (2.15) to obtain
$\boldsymbol {u}_{l}^{\prime }$ at bubble locations, and evaluate forces in bubbles via (2.7).(vii) Update bubble velocities by integrating (2.6b), evaluate
$\boldsymbol {M}_{b}$ and $\nu _{B,b}$, and interpolate back to SPH particle positions through (3.5).(viii) Evaluate remaining terms in (3.14a), and evaluate
$\tilde {\boldsymbol {u}}_{l}^{\star }$.(ix) Construct and solve (3.17) to obtain
$\alpha ^{n+1}\tilde {p}^{n+1}$, and hence $\tilde {p}^{n+1}$. The system (3.17) is solved using a BiCGStab algorithm with Jacobi preconditioning, with Neumann boundary conditions on solid surfaces, and homogeneous Dirichlet boundary conditions on free surfaces as in King & Lind (Reference King and Lind2021).(x) Evaluate pressure gradient, and project
$\tilde {\boldsymbol {u}}_{l}^{\star }$ onto $\tilde {\boldsymbol {u}}_{l}^{n+1}$ using (3.14b).(xi) Advect particles to final positions
$\boldsymbol {r}^{n+1}=\boldsymbol {r}^{n}+\delta {t}\,(\tilde {\boldsymbol {u}}_{l}^{n}+\tilde {\boldsymbol {u}}_{l}^{n+1}+2\boldsymbol {u}_{ps})/2$.
The value of 
$\delta {t}$ is set adaptively according to criteria for the Courant condition and viscous diffusion, as in Xenakis et al. (Reference Xenakis, Lind, Stansby and Rogers2015):
\begin{equation} \delta{t}=0.2\min\left(\frac{h}{\max_{\mathcal{P}}\left(\lvert\tilde{\boldsymbol{u}}_{0}\rvert\right)},Re\,{h}^{2}\right), \end{equation}
in which 
$\max _{\mathcal {P}}$ is the maximum value over the set of all particles 
$\mathcal {P}$. Although our numerical scheme is capable of capturing acoustic waves in the liquid, this does not impose an additional (stability related) constraint on the time step, as the acoustic physics is treated implicitly in the fractional step approach. The acoustic part of the system is unconditionally stable, although for larger time steps, the acoustic waves are subject to greater numerical dissipation. There is a trade-off here, between computational costs, and the degree to which acoustic information is of interest. Regardless of the application, there is a benefit of the present approach over a perfectly incompressible approach. Whilst both methods result in smooth pressure fields with no spurious noise (as is commonly experienced in explicit, weakly compressible SPH), the additional terms in (3.17) accounting for compressibility increase the diagonal dominance of the linear system that represents the discrete form of this equation. This increased diagonal dominance renders the system more amenable to solution using iterative methods, which will converge more quickly. In the present application, we are not interested in capturing the acoustic waves, and do not impose an additional time step constraint proportional to 
$h\,Ma$. Despite this, for 
$Ma=0.05$, we still obtain a reduction in the number of iterations necessary to solve (3.17) by a factor of typically approximately 
$5$. For our SPH framework, the solution of (3.17) is the most expensive aspect of the simulation, so this presents a significant computational saving. With this in mind, we consider 
$Ma$ in this framework to be a numerical parameter (rather than a physical one), as in weakly compressible SPH schemes. The value 
$Ma=0.05$ is consistent with the artificial sound speeds widely used in weakly compressible SPH, ensures that the density fluctuations are well below 
$1\,\%$ (Monaghan Reference Monaghan1994), and is used throughout this work.
As mentioned in § 3.1.1, SPH particle volumes are adjusted to account for bubble volumes, resulting in a redistribution of SPH particles over a finite time to retain an approximately uniform particle number density. To verify that this effect is small, we perform a simulation over a triply periodic cubic domain with unit side length, and the liquid initially at rest. We set 
$\delta {r}=1/20$, and neglect gravity. After 
$10$ SPH time steps, we introduce a stationary bubble at the centre of the domain. First, we observe that the resulting pressure and velocity fields remain zero everywhere, as does the velocity of the bubble. The response to the addition of the bubble is purely numerical, in that only the particle distribution is affected. We calculate the 
$L_{2}$ norm of the magnitude of 
$\boldsymbol {u}_{ps}$ over the domain, which, as it is proportional to the SPH particle concentration gradient, allows us to quantify the effect. Figure 3 shows 
$\|\boldsymbol {u}_{ps}\|_{2}$ plotted against time step number for several values of bubble radius 
$a_{b}$. For all bubble sizes, there is an initial peak when the bubble is created, where the SPH particle volumes have increased to accommodate the bubble, but no particle redistribution has taken place. The magnitude of this peak is proportional to the cube of the bubble radius, as expected. For bubbles with 
$a_{b}\le 0.1\delta {r}$, the peak in 
$\|\boldsymbol {u}_{ps}\|_{2}$ is below 
$10^{-6}$. After the peak, the traces all decay. The decay is initially exponential, with characteristic time approximately 
$2\delta {t}/3$ (slope lines shown in red in figure 3) for all bubble radii, but at later times, the decay rate slows. We believe that this reduction in decay rate is a consequence of the finite volume of the bubble that must be accommodated within the particle distribution. For larger bubbles, the particle redistribution must be less localised in space, and consequently in time. We have confirmed this by running this test for a range of different values of 
$\delta {t}$ spanning an order of magnitude. We performed the same test with a fixed 
$a_{b}=0.4\,\delta {r}$, for a range of values of 
$h_{0}/\delta {r}\in [1.1,1.5]$ (noting that all other simulations in this work are performed with 
$h/\delta {r}=1.3$), and found a negligible variation in the peak value of 
$\|\boldsymbol {u}_{ps}\|_{2}$ with 
$h_{0}/\delta {r}$, whilst the initial decay time scale varied from 
$0.55\delta {t}$ for 
$h_{0}/\delta {r}=1.1$, to 
$0.8\delta {t}$ for 
$h_{0}/\delta {r}=1.5$. In all these cases, the characteristic time of decay is less than the value of the SPH time step – in other words, the particle redistribution happens quickly. With the results shown in figure 3 in mind, we observe that the mean value of 
$\|\boldsymbol {u}_{ps}\|_{2}$ for the simulation of a breaking wave studied in § 5 is 
$1.4\times {10}^{-2}$. Whilst there are physical processes in our model occurring on shorter time scales than the redistribution process (e.g. bubble breakup, sub-resolution fluctuations), the effect of particle redistribution due to the presence of bubbles is more than an order of magnitude smaller than the particle redistribution due to the fluid motion, which is a standard aspect of SPH.
Variation of the 
$L_{2}$ norm of the shifting velocity magnitude 
$\lvert \boldsymbol {u}_{ps}\rvert$ with time (in units of 
$\delta {t}$) for a triply periodic unit cube of fluid at rest, subject to the instantaneous addition of a single bubble. The different patterned lines indicate different bubble sizes 
$a_{b}$ relative to the SPH particle spacing 
$\delta {r}$. The red lines indicate an exponential decay with characteristic time 
$2\delta {t}/3$.
3.4. Bubble entrainment, breakup and free-surface interaction
3.4.1. Entrainment
Our intention is to simulate flows where bubbles are entrained at the free surface. We use an entrainment model, similar in ethos to that of Ma et al. (Reference Ma, Shi and Kirby2011) and Derakhti & Kirby (Reference Derakhti and Kirby2014), in which a fraction of the turbulent kinetic energy of the liquid is assumed to be converted into surface energy as bubbles are created (or entrained) at the free surface. Based on this energy balance, the energy available for bubble creation at SPH particle 
$i$ in a given time step is
\begin{equation} E_{bc,i}=C_{\varepsilon}\beta\alpha_{i}\varepsilon_{i}V_{i}\,\delta{t}, \end{equation}
where 
$C_{\varepsilon }=0.01$ is a constant set empirically. The surface energy of a bubble of radius 
$a_{b}$ is
\begin{equation} E_{se}=\frac{4{\rm \pi}{a}_{b}^{2}}{We}. \end{equation}
The closure models used to evaluate the forces on each bubble are based on the assumption of spherical non-interacting bubbles (Fraga et al. Reference Fraga, Stoesser, Lai and Socolofsky2016). When the concentration of bubbles is large (and hence 
$\alpha$ is small), these assumptions cease to be valid. Therefore, we impose an additional constraint on bubble entrainment, such that 
$1-\alpha \lesssim 1/3$, by denoting the volume available for bubble entrainment as
\begin{equation} V_{bc,i}=\frac{1}{W(0)}\left(\frac{3}{2}-\frac{1}{\alpha_{i}}\right), \end{equation}
where 
$W(0)$ is the maximum value of the SPH kernel. Here, 
$V_{bc,i}$ is the volume of bubbles that must be entrained at SPH particle 
$i$ to result in 
$\alpha _{i}=2/3$. We note that our algorithm does not impose a hard limit 
$\alpha >2/3$, but rather an approximate limit, as it is possible for bubbles to converge on a particular region of the flow subsequent to entrainment. Although this limit may influence the resulting physics of our simulations, it is a necessary limitation of the model of this form: we can either limit 
$\alpha$, or permit 
$\alpha$ to stray into territory where our underlying assumptions cease to be valid. Given that spherical bubbles in a cubic-packed lattice would yield a volume fraction 
$\alpha \approx 0.52$, this limit does not seem overly stringent. Furthermore, we note that the simulations of breaking waves in Derakhti & Kirby (Reference Derakhti and Kirby2014) yielded values 
$\alpha >2/3$. Following Derakhti & Kirby (Reference Derakhti and Kirby2014), we entrain bubbles only at the free surface, and only when 
$\varepsilon$ exceeds a certain threshold. Where others (Ma et al. Reference Ma, Shi and Kirby2011; Derakhti & Kirby Reference Derakhti and Kirby2014) treated the bubbles as a continuum, we treat them discretely, hence our algorithm differs somewhat, despite the principles of energy balancing being the same. Where Derakhti & Kirby (Reference Derakhti and Kirby2014) modelled the bubbles through a set of bubble groups, each with a characteristic size, we are able to model individual bubbles, and hence obtain a continuous bubble size distribution. To achieve this, our entrainment algorithm is as follows.
At every time step, for every free-surface particle 
$i\in \mathcal {P}_{FS}$ for which 
$\varepsilon _{i}>0.2$ and 
$\alpha >2/3$, proceed through the following steps.
(i) Evaluate the available energy
$E_{bc,i}$ and volume $V_{bc,i}$ for bubble creation, according to (3.19) and (3.21). Initialise counters for new bubbles, and new potential bubbles: $n_{nb}=0$ and $n_{npb}=0$.(ii) Generate a potential bubble radius
$a_{b,p}\in [We^{-3/5}/40,\delta {r}]$ with uniform probability distribution. Increment the new potential bubble counter: $n_{npb}=n_{npb}+1$.(iii) Evaluate the surface energy
$E_{se,p}$ of the potential bubble via (3.20). If $E_{bc,i}\ge {E}_{se,p}$, then there is enough energy. If $E_{bc,i}< E_{se,p}$ and $E_{bc,i}/E_{se,p}>\zeta _{E}$, where $\zeta _{E}\in [0,1]$ is a uniformly distributed random number, then there is not enough energy for this potential bubble.(iv) Evaluate the volume of the potential bubble,
$V_{b,p}=4{\rm \pi} {a}_{b,p}^{3}/3$. If $V_{bc,i}\ge {V}_{b,p}$, then there is enough volume available to create this bubble. If $V_{bc,i}< V_{b,p}$ and $V_{bc,i}/V_{b,p}>\zeta _{V}$, where $\zeta _{V}\in [0,1]$ is a uniformly distributed random number, then there is not enough volume available for this potential bubble.(v) If the checks in steps (iii) and (iv) were passed, then make a new bubble
$j$ with radius $a_{b,j}=a_{b,p}$, position $\boldsymbol {r}_{b,j}=\boldsymbol {r}_{i}+\boldsymbol {\zeta }_{r}\,\delta {r}$, where $\boldsymbol {\zeta }_{r}$ is a random vector with elements in $[0,1]$, $\boldsymbol {u}_{b,j}=\tilde {\boldsymbol {u}}_{l,i}$, and $\delta {r}$ is the initial particle spacing. Denote the time of bubble creation as $T_{b,j}$. Increment the new bubble counter: $n_{nb}=n_{nb}+1$.(vi) If the checks in steps (iii) and (iv) were passed, reduce the available energy and volume by
(3.22a)$$\begin{gather} E_{bc,i}=E_{bc,i}-E_{se,p}, \end{gather}$$(3.22b)$$\begin{gather}V_{bc,i}=V_{bc,i}-V_{b,p}. \end{gather}$$(vii) Check whether to continue entraining bubbles. If all the inequalities
(3.23a)$$\begin{gather} E_{bc,i}>0, \end{gather}$$(3.23b)$$\begin{gather}V_{bc,i}>0, \end{gather}$$(3.23c)$$\begin{gather}n_{nb}<10, \end{gather}$$(3.23d)are satisfied, return to step (ii). Otherwise, move on to the next SPH particle.$$\begin{gather}n_{npb}< n_{nb}+5 \end{gather}$$
The addition of the stochastic processes in steps (iii) and (iv) to allow some bubbles for which 
$E_{se,p}>E_{bc,i}$ and 
$V_{b,p}>V_{bc,i}$, based on the remainders, prevents the creation of large bubbles from being overly suppressed.
The maximum potential bubble size 
$a_{b,p}\le \delta {r}$ is set to ensure that the bubbles are smaller than the implicit LES filter scale 
$\tilde {\varDelta }$. The minimum potential bubble size 
$a_{b,p}\ge {We}^{-3/5}/40$ is equal to 
$0.05a_{H}$, where 
$a_{H}$ is the Hinze scale bubble radius evaluated a priori. This limit is imposed to prevent the generation of a very large number of very small bubbles, which would impose a significant computational cost on the simulation, but which play little role in the overall dynamics of the flow. In future, we consider that the present framework could be extended to treat bubbles with 
$a_{b}\ll {a}_{H}$ as a continuum, as is done for all bubble sizes in Eulerian–Eulerian schemes. The value of 
$C_{\varepsilon }$ influences the total entrained volume, but has little influence on the bubble distribution in time, space and size. For all simulations of breaking waves, we set 
$C_{\varepsilon }=0.01$, which provides an approximate match in terms of the total entrained bubble volume with the direct numerical simulations results of Deike et al. (Reference Deike, Melville and Popinet2016).
3.4.2. Breakup
As in work by Ma et al. (Reference Ma, Shi and Kirby2011), Martínez-Bazán et al. (Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montañes and Lasheras2010, Reference Martínez-Bazán, Montañes and Lasheras1999a,Reference Martínez-Bazán, Montañes and Lasherasb), Chan, Johnson & Moin (Reference Chan, Johnson and Moin2021a) and Chan et al. (Reference Chan, Johnson, Moin and Urzay2021b), we use a stochastic breakup model based on energy balance considerations. Our model is similar to that of Martínez-Bazán et al. (Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montañes and Lasheras2010), and is based on the imbalance between the surface restoring pressure and the stress on the bubble surface due to the motion of the liquid. For a bubble of size 
$a_{b}$, the surface restoring pressure is 
$6/(2a_{b}\,We)$. The stress exerted on the bubble by the turbulent motion of the liquid is 
$\tfrac {1}{2}{c}_{def}(2\varepsilon {a}_{b})^{2/3}$, where 
$c_{def}=8.2$ as given in Batchelor (Reference Batchelor1953). The net deforming stress on the bubble is given by the difference, from which we obtain an expression for a characteristic deformation velocity
\begin{equation} u_{def}=\text{sgn}\left(c_{def}\left(2\varepsilon{a}_{b}\right)^{{2}/{3}} -\frac{6}{{a}_{b}\,{We}}\right)\sqrt{\left\lvert{c}_{def}\left(2\varepsilon{a}_{b}\right)^{{2}/{3}} -\frac{6}{{a}_{b}\,{We}}\right\rvert}. \end{equation}In works such as Ma et al. (Reference Ma, Shi and Kirby2011), Martínez-Bazán et al. (Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montañes and Lasheras2010) and Chan et al. (Reference Chan, Johnson and Moin2021a,Reference Chan, Johnson, Moin and Urzayb), which model the dispersed phase through a continuous bubble number density field, a characteristic time scale of breakup is often given by
\begin{equation} T_{bu}=\frac{2a_{b}}{u_{def}} =\frac{2a_{b}}{\sqrt{c_{def}\left(2\varepsilon{a}_{b}\right)^{{2}/{3}}-{6}/({a}_{b}\,{We})}}, \end{equation}
assuming 
$u_{def}>0$. In this work, we treat each bubble individually, hence we can construct a model that attempts to account for the residence time of the bubble, and the fact that deformation and breakup occur over a finite time (Risso & Fabre Reference Risso and Fabre1998). For every bubble, we integrate numerically the deformation velocity to obtain a ‘deformation distance’
\begin{equation} L_{def}=\max\left\{\displaystyle\int_{t_{bc,i}}^{t}u_{def}(\tau)\,{\rm d}\tau,0\right\}, \end{equation}
which is a measure of how deformed the bubble is. Here, 
$t_{bc,i}$ is the time of creation of bubble 
$i$. Since we allow 
$u_{def}<0$, (3.26) accounts for both deformation of the bubble due to turbulence, and relaxation of the bubble back towards a spherical shape. We assume that the bubble breaks once the deformation distance exceeds the bubble diameter, i.e. when 
$L_{def}>2a_{b}$.
We assume that all breakup events are binary; that is, a parent bubble splits into two child bubbles. Once a bubble has been marked for a breakup event, we set 
$V=V_{b,child}/V_{b,parent}\in [0,1]$ randomly with the probability distribution given by Martínez-Bazán et al. (Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montañes and Lasheras2010), where
\begin{align} P(V)\propto \begin{cases}V^{{-}2/3}(1-V)^{{-}2/3}(V^{2/9}-\varLambda^{5/3}) \left[(1-V)^{2/9}-\varLambda^{5/3}\right], & V\in\left[V_{min},V_{max}\right],\\ 0, & \text{otherwise}.\end{cases} \end{align}
Here, 
$V_{min}$ is set to limit the size distributions to regions where the expression for 
$P(V)$ is positive, with an additional hard limit 
$V_{min}\ge 10^{-3}$, to ensure that breakup events do not occur in a capillary-action-dominated regime. Also, 
$P(V)$ is even about 
$V=1/2$, and 
$V_{max}=1-V_{min}$. To normalise 
$P(V)$, and evaluate the cumulative density function (required for generating a random number with distribution 
$P(V)$), we integrate (3.27) numerically. Here, 
$\varLambda$ represents a critical bubble radius (relative to the parent bubble), below which the confining stresses due to surface tension exceed the turbulent stresses – it is the largest bubble that will not break due to the turbulent flow (at a given instant in time and space). For each breaking event, we evaluate 
$\varLambda$ as
\begin{equation} \varLambda=\left(\frac{12}{c_{def}\,We}\right)^{3/5}\bar{\varepsilon}^{{-}2/5}\left(2a_{b,p}\right)^{{-}1}, \end{equation}where
\begin{equation} \bar{\varepsilon}=\frac{1}{t-t_{bc,i}}\displaystyle\int_{t_{bc,i}}^{t}\varepsilon (\tau)\,{{\rm d}}\tau \end{equation}
is the mean dissipation rate of the bubble lifetime, with 
$\varepsilon$ the instantaneous dissipation rate in the liquid at bubble 
$i$. In the above model, a critical value 
$\varLambda _{crit}=2^{-6/45}\approx 0.912$ exists, and the model is not valid for 
$\varLambda >\varLambda _{crit}$ (
$P(V)<0$ 
$\forall {V}$). Effectively, this imposes an additional limit on the smallest bubble that can break, which is not necessarily consistent with the breakup criteria obtained by integration of 
$L_{def}$. Even if the evaluation of 
$L_{def}$ from (3.26) suggests that a bubble should break, if it has a corresponding 
$\varLambda >\varLambda _{crit}$, then it does not break, as the child bubble sizes are undefined in (3.27). In practice, this situation rarely occurs, and has negligible impact on the overall bubble population dynamics. Figure 4 shows 
$P(V)$ for various values of 
$\varLambda$. We see that for small 
$\varLambda$, the child size distribution is largely flat, but with peaks at very small and very large bubbles (more apparent in the inset). As 
$\varLambda$ increases, these peaks reduce, and for 
$\varLambda =0.4$, the distribution is nearly flat. For larger 
$\varLambda$, the possible range of child sizes decreases, with an increasing dominance of equal-sized breakup events. For 
$\varLambda =\varLambda _{crit}$, all breakup events result in 
$V=1/2$.
The probability density function of child bubble volumes for different values of 
$\varLambda$, given by the breakup model of Martínez-Bazán et al. (Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montañes and Lasheras2010).
Whilst binary breakup models have been utilised in a range of works (Martínez-Bazán et al. Reference Martínez-Bazán, Montañes and Lasheras1999a,Reference Martínez-Bazán, Montañes and Lasherasb, Reference Martínez-Bazán, Rodríguez-Rodríguez, Deane, Montañes and Lasheras2010; Ma et al. Reference Ma, Shi and Kirby2011; Derakhti & Kirby Reference Derakhti and Kirby2014), recent work (Rivière et al. Reference Rivière, Mostert, Perrard and Deike2021; Ruth et al. Reference Ruth, Aiyer, Rivière, Perrard and Deike2022) has suggested that the binary breakup mechanism presents some limitation. Direct numerical simulations of individual bubbles breaking in isotropic turbulence (Rivière et al. Reference Rivière, Mostert, Perrard and Deike2021) show that even for a binary breakup event, a bubble filament joining breaking bubbles forms, and this filament ruptures into one or many small bubbles due to capillary effects. This mechanism provides a non-local element to the bubble size cascade, and has been suggested as being responsible for the majority of sub-Hinze-scale bubble formation (Ruth et al. Reference Ruth, Aiyer, Rivière, Perrard and Deike2022). A breakup model constructed on the framework above, but which includes the formation of multiple small bubbles by capillary action, is an active area of development for us. However, at this stage we prefer to adopt a simple binary breakup model, and focus on the free-surface–bubble interactions.
3.4.3. Free-surface interactions
In the vicinity of free surfaces, the assumption that the fluid around the bubble is uniform and infinite does not hold, and the closure models for the forces in (2.7) are not valid. The dynamics of the interactions between bubbles and free surfaces are complex, even for individual bubbles in stationary reservoirs, and can include a direct rise to surface bursting, oscillatory bouncing behaviour (Sanada, Watanabe & Fukano Reference Sanada, Watanabe and Fukano2005; Suñol & González-Cinca Reference Suñol and González-Cinca2010), and long-term persistence of the bubble at the free surface. Despite this, most mesh-based Eulerian–Lagrangian schemes simply represent the free surface via a rigid free-slip condition (e.g. Fraga et al. Reference Fraga, Stoesser, Lai and Socolofsky2016), or through a volume-of-fluid approach to resolve the gas phase above the free surface (e.g. Zhang & Ahmadi Reference Zhang and Ahmadi2005; Pan et al. Reference Pan, Johansen, Olsen, Reed and Sætran2021). In the latter approach, both the deceleration of bubbles on approach to the free surface and the free-surface deformation are captured, although this occurs solely through modification of the relative densities in the closure models used to evaluate the forces on the bubbles. These closure models still assume a uniform liquid field surrounding the bubbles, and the complex and small-scale physics involved are not included. Detailed modelling of free-surface–bubble interactions are beyond the scope of this work. However, some mechanism to describe the behaviour of the interaction is necessary, to prevent bubbles from simply rising due to buoyancy and leaving the domain. For each SPH particle, a surface normal vector 
$\boldsymbol {n}$ is evaluated (and smoothed), as in King & Lind (Reference King and Lind2021). The surface normal vectors are then interpolated to each bubble position through (3.6) to obtain 
$\boldsymbol {n}_{b}$. For each bubble, we construct a parameter 
$\psi _{fs}$ that identifies when a bubble is in proximity to the free surface:
\begin{equation} \psi_{fs}=\frac{\left\lvert\boldsymbol{n}_{b}\right\rvert}{\left\lvert\boldsymbol{n}_{0}\right\rvert{V}_{lb}}, \end{equation}
where 
$\lvert \boldsymbol {n}_{0}\rvert$ is the magnitude of the surface normal vector at a plane surface, and 
$V_{lb}$ is the SPH volume interpolated to the bubble location (i.e. (3.6) with 
$\phi =1$). This 
$\lvert \boldsymbol {n}_{0}\rvert$ is dependent only on the choice of SPH kernel and ratio of smoothing length to SPH particle spacing, and is hence constant across all our simulations, taking the precomputed value 
$\lvert \boldsymbol {n}_{0}\rvert =0.353$. For bubbles far from the free surface, 
$\psi _{fs}\approx 0$. For bubbles within the support radius of free-surface SPH particles, 
$\psi _{fs}$ increases smoothly, to approximately 
$1$ when bubbles are located at the free surface. Denoting 
$\hat {\boldsymbol {n}}_{b}=\boldsymbol {n}_{b}/\lvert \boldsymbol {n}_{b}\rvert$, the relative normal velocity between a bubble and the free surface is 
$\boldsymbol {u}_{rel}\boldsymbol {\cdot }\hat {\boldsymbol {n}}_{b}$. Bubbles with 
$\psi _{fs}>0.1$ and 
$t-t_{bc}>T_{c}$ are flagged to interact with the free surface. Here, 
$T_{c}=\delta {r}/\lvert \boldsymbol {u}_{rel}\boldsymbol {\cdot }\hat {\boldsymbol {n}}_{b}\rvert$ is a threshold age (which varies as a bubble evolves), designed to prevent bubbles from being destroyed at the moment of entrainment. When a bubble is marked for free-surface interaction, it is given an expected merge time 
$t_{m}=t+T_{c}$, at which it will be located at the free surface. For bubbles interacting with the free surface, (2.6b) is modified as
\begin{equation} V_{b,i}\,\frac{{\rm d}\boldsymbol{u}_{b,i}}{{\rm d}t}=\boldsymbol{F}_{d}+\boldsymbol{F}_{l}+\boldsymbol{F}_{vm}+\boldsymbol{F}_{g}+\boldsymbol{F}_{fs}, \end{equation}
where the free-surface interaction force 
$\boldsymbol {F}_{fs}$ is
$$\begin{gather} \boldsymbol{F}_{fs}=\frac{1}{2}(1+\text{erf}(2\log(5\psi_{fs})))\times \frac{\boldsymbol{n}_{b}}{\left\lvert\boldsymbol{n}_{b}\right\rvert}\left\{\left\lvert\boldsymbol{u}_{rel}\boldsymbol{\cdot}\hat{\boldsymbol{n}}_{b}\right\rvert(1+C_{vm}\beta)\,\frac{V_{b}}{\max\left(t_{m}-t,\delta{t}\right)} \right.\nonumber\\ - \left.\left[\boldsymbol{F}_{b}+\boldsymbol{F}_{d}+\boldsymbol{F}_{l}+\beta{C}_{vm}V_{b}\,\frac{{\rm d}\boldsymbol{u}_{l}}{{\rm d}t}\right]\boldsymbol{\cdot}\frac{\boldsymbol{n}_{b}}{\left\lvert\boldsymbol{n}_{b}\right\rvert}\right\}. \end{gather}$$
The first term, 
$\tfrac {1}{2}(1+\text {erf}(2\log (5\psi _{fs})))$, varies smoothly from 
$0$ when 
$\psi _{fs}$ is small, to 
$1$ for large 
$\psi _{fs}$, and is approximately 
$1$ for 
$\psi _{fs}>0.4$. This term ensures that the free-surface interaction force is switched on smoothly as a bubble approaches the free surface. The free-surface interaction force decelerates the bubble in the direction normal to the free surface, such that the bubble velocity approaches the free-surface velocity over the time scale 
$T_{c}$ (or 
$\delta {t}$, whichever is greater). The momentum exchange is modified to include 
$\boldsymbol {F}_{fs}$ as
\begin{equation} \boldsymbol{M}=(-\boldsymbol{F}_{d}-\boldsymbol{F}_{l}-\boldsymbol{F}_{vm}-\boldsymbol{F}_{fs})/\beta. \end{equation}
Additionally, we impose a special treatment for bubbles interacting with spray. Where an individual SPH particle 
$i$ becomes separated from the bulk of the fluid (such that 
$\mathcal {P}_{i}$ contains only 
$i$), it is subjected only to a gravitational body force, and all other terms in (2.3) are set to zero. If a bubble has fewer than 
$20$ (in three dimensions) SPH particle neighbours, and all of those SPH particle neighbours are identified as free-surface particles, then the bubble is discounted, and assumed to have burst. This provides increased stability in regions of violent spray, such as around breaking waves, and prevents bubbles from becoming completely isolated from the SPH simulation.
The effect of the free-surface interaction model is demonstrated in figure 5. A tank of still water is simulated, with an individual bubble rising due to buoyancy. We vary 
$a_{b}$, and taking the characteristic velocity scale as the terminal velocity of the bubble 
$u_{t}$, the bubble Reynolds number varies as 
$Re=10^{5}u_{t}a_{b}$, and the bubble scale Weber number is 
$We_{b}=1.4\times {10}^{4}a_{b}u_{t}^{2}$. The bubble trajectory is shown in figure 5(a), and rise velocity in figure 5(b), for several values of 
$We_{b}$. In both plots, a time shift has been applied so that the coordinate on the abscissa corresponds to the time since the bubble–surface interaction began. We see that for larger 
$We_{b}$, the bubble is decelerated more quickly. For all 
$We_{b}\in [0.01,4.06]$, the bubble comes to rest on the free surface. For larger bubble Weber numbers 
$We_{b}$ (corresponding to larger 
$a_{b}$), the system becomes less stable after the bubble reaches the surface, and this is due to the settling over a finite time of the SPH particles to accommodate the bubble volume. In all cases, though, the final position of the bubble remains within a distance 
$\delta {r}$ of the free-surface location.
Trajectories of individual bubbles approaching a free surface, for a range of bubble Weber numbers 
$We_{b}$. (a) Variation of distance to the free surface, scaled with particle spacing 
$\delta {r}$, with dimensionless time. (b) Variation of relative normal velocity 
$\lvert \boldsymbol {u}_{rel}\boldsymbol {\cdot }\hat {\boldsymbol {n}}_{b}\rvert$ with dimensionless time, scaled with the terminal velocity of the rising bubble. In (a,b), a time shift has been applied such that 
$t=0$ corresponds to the time at which the free-surface interaction begins.
The final feature of the free-surface interaction model is the persistence time. Once a bubble has reached the free surface, it remains subject to (3.31) for a persistence time 
$T_{p}$. If the motion of the liquid phase is such that the bubble (now moving with the liquid) moves away from the free surface (i.e. if 
$\psi _{fs}<0.1$), then the bubble is assumed to no longer interact with the free surface, and its motion is again governed by (2.6b). In reality, when bubbles reach a free surface, a thin lubrication layer forms, and the liquid in this layer drains until the layer ruptures (see e.g. Modini et al. (Reference Modini, Russell, Deane and Stokes2013) for a description of the mechanism). This process is complex, and is strongly influenced by the local geometry and the salinity (Scott Reference Scott1975) (or other contaminants and surfactants) and gradients thereof. We cannot seek to capture this process accurately in our model, and instead we seek an order of magnitude estimate for 
$T_{p}$ that has a physical basis. We use the (here made non-dimensional) expression of Poulain, Villermaux & Bourouiba (Reference Poulain, Villermaux and Bourouiba2018), where
\begin{equation} T_{p}=\frac{We^{3/4}}{Re}\,Fr^{1/2}\,Sc\,{a}_{b}^{1/2},\end{equation}
in which 
$Sc$ is the Schmidt number, taken as 
$Sc=700$ for sea water at 
$20\,^{\circ }{\rm C}$ (De Bruyn & Saltzman Reference De Bruyn and Saltzman1997). Equation (3.34) is based on a mechanistic model for film drainage, accounting for molecular diffusion and curvature-pressure-induced flow, and provides a rough scaling of the persistence time, with a level of approximation that is consistent with the present framework; a more detailed physical model is beyond the scope of this work. Bubbles that remain in proximity to the free surface beyond 
$t_{m}$ are assumed to burst (and are removed from the simulation) at 
$t=t_{m}+T_{p}$.
4. Bubble plumes
We first use our model to simulate a buoyant bubble plume in a tank. Our simulation is configured to match the experimental set-up of Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016) as follows. The domain is a unit cube of liquid, with no-slip conditions at the lateral and lower boundaries, and a free-surface boundary at the upper surface of the liquid. The origin of our coordinate system is at the centre of the base of the tank, with the unit vector in the 
$z$ direction pointing upwards. A bubble sparger is simulated, centred at 
$(x,y,z)=(0,0,0.09)$, with radius 
$r_{0}=0.02$. The governing parameters are 
$\beta =1000/1.2$, 
$Re=10^{6}$, 
$We=1.4\times {10}^{4}$ and 
$Fr=1/\sqrt {9.81}$. The (dimensionless) volumetric flow rate of the bubble sparger is 
$\dot {Q}=1.67\times {10}^{-5}$. With this dimensionless scaling, the Hinze scale is estimated at 
$a_{H}=We^{-3/5}/2=1.62\times {10}^{-3}$.
Initially, we simulate the release of uniform bubbles with 
$a_{b}=10^{-3}$, as in Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016). Figure 6 shows the velocity and vorticity magnitude fields in the 
$y=0$ plane at 
$t=12$. Figure 7(a) shows the evolution of the mean dissipation rate in the liquid for several resolutions. For all resolutions, there are differences in the dissipation rate, although for 
$\delta {r}<1/100$, the dissipation rates converge approximately for 
$t>6$. We note that for the isotropic turbulence with 
$Re=10^{6}$ investigated in the Appendix, a resolution 
$\delta {r}=1/100$ is sufficient to yield accurate dissipation rates compared with the high-order reference data in Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021). For the remainder of this section, we set 
$\delta {r}=1/100$. The increase in high-frequency noise in the mean dissipation rate at approximately 
$t=6$ visible in figure 7(a) is associated with the bubble plume reaching the free surface. Whilst previous authors (e.g. Fraga et al. Reference Fraga, Stoesser, Lai and Socolofsky2016) have used finite-volume methods with a fixed free surface, our SPH simulation captures the motion of the free surface, and includes a model for the interaction of the bubbles with the free surface. The high-frequency noise is a numerical artefact of the SPH scheme: as the liquid upwells in the centre of the tank and moves radially outwards, SPH particles join and leave the free surface, resulting in a small redistribution of the SPH particle arrangement.
(a) Velocity and (b) vorticity magnitude of the flow around the bubble plume at 
$t=12$. Bubbles are shown as black circles (not to scale) in (a). The images show a slice through the plane 
$y=0$.
(a) Evolution of the mean dissipation rate in the liquid for several resolutions. (b) Variation of dimensionless mean vertical velocity in liquid with dimensionless radial position. Blue symbols indicate the results of our simulations at depths 
$z=0.555$ (circles) and 
$z=0.453$ (triangles), which are compared with the self-similar Gaussian solution (dashed black line), and the experimental (red crosses) and numerical (red circles) results of Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016).
Figure 7(b) shows the variation of the mean vertical velocity in the liquid with radial position, at different depths. Here, the vertical velocity is scaled by the (local) velocity at the plume centre 
$w_{c}(z)$, and the radial position is scaled by the local plume width 
$b_{v}$, taken as 
$w(r=b_{v}(z),z)={\rm e}^{-1}\,w_{c}(z)$. The blue symbols indicate the results of our SPH framework just above (
$z=0.555$, circles) and below (
$z=0.453$, triangles) the centre of the tank. The numerical (red circles) and experimental (red crosses) data from Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016) are also shown, as is a Gaussian distribution (dashed black line). We see agreement at both vertical locations with both the Gaussian profile and the results of Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016). For larger 
$r/b_{v}>1.5$, our results deviate from the Gaussian profile, but match those of Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016), tending towards a non-zero 
$w/w_{c}$. As discussed in Fraga et al. (Reference Fraga, Stoesser, Lai and Socolofsky2016), this deviation from the Gaussian is due to the confined plume in our simulations.
We now simulate a plume in the same configuration, but with a polydisperse bubble size distribution. The initial bubble size distribution is a truncated Gaussian with mean 
$a_{b}=10^{-3}$, standard deviation 
$10^{-3}$, and a minimum radius cut-off at 
$a_{b}=10^{-4}$. Figure 8 plots the magnitude of the relative velocity 
$\lvert \boldsymbol {u}_{b}-\tilde {\boldsymbol {u}}_{l}\rvert$ against the bubble size near the centre of the tank, at 
$z=0.555$. Each data point corresponds to a separate bubble, and the data are collected for 
$t\in [10,20]$. We see a clear, almost linear trend, with larger bubbles rising faster relative to the liquid than smaller ones. There is also an obvious discrepancy between the results with (black symbols) and without (red symbols) the Langevin model for sub-resolution fluctuations. The Langevin model increases the variation of the relative velocity (without changing the mean). This observation is consistent with figure 9, which shows the traces of a random sample of bubbles as they rise through the plume. The positions of bubbles are projected onto polar coordinates 
$(r/r_{0},z)$, and are coloured by the bubble radius. Figure 9(a) shows bubble traces where we set 
$\boldsymbol {u}_{l}^{\prime }=\boldsymbol {0}$, whilst in figure 9(b), we calculate 
$\boldsymbol {u}_{l}^{\prime }$ using the Langevin model. Without the Langevin model, there is a clear trend for larger bubbles to migrate away from the centre of the plume, and smaller bubbles to rise more vertically. This self-organising phenomenon has been observed previously in numerical simulations (Fraga & Stoesser Reference Fraga and Stoesser2016), and experimentally (Ye et al. Reference Ye, Zhou, Shao and Niu2022). The phenomenon can be explained by the linear dependence of the lift force on the relative velocity: as bubbles rise through the plume, the relative velocity is close to orthogonal to the gradient of the velocity in the liquid, and as larger bubbles have a larger relative velocity, they experience an increased lift force, resulting in greater lateral migration. When the Langevin model is included, the bubbles still self-organise, but there is increased lateral migration of small bubbles due to the sub-resolution fluctuations. This self-organisation is depicted clearly in figure 10, which plots the normalised bubble radius 
$a_{b}/a_{H}$ against radial position 
$r/r_{0}$ at several depths. At 
$z=0.1$ (figure 10a), just above the base of the plume, the bubble sizes are distributed evenly across the radius of the plume. Further up the plume, the bubbles migrate radially outwards, and the plume gets wider. There is a clear trend, both with (black symbols) and without (red symbols) the Langevin sub-resolution model, for larger bubbles to migrate further outwards, as can be seen in figures 10(b–d). When the Langevin model is included, this trend remains, although as observed in figure 9, the lateral migration of smaller bubbles is increased – the greater spread in radial position of small bubbles is visible in figures 10(b–d). This behaviour is expected, as small-scale fluctuations have a greater effect on small bubbles than larger bubbles through the increased relative importance of the drag force.
Variation in the relative velocity magnitude (excluding fluctuations) 
$\lvert \boldsymbol {u}_{b}-\tilde {\boldsymbol {u}}_{l}\rvert$ with bubble size at depth 
$z=0.555$, for the polydisperse bubble plume, both with (black symbols) and without (red symbols) the Langevin model for sub-resolution velocity fluctuations.
Bubble traces for the polydisperse plume, (a) without and (b) with the Langevin model for sub-resolution velocity fluctuations. Bubble positions are projected onto polar coordinates 
$(r/r_{0},z)$, where 
$r_{0}=0.02$ is the radius of the plume source. Traces are coloured by bubble radius 
$a_{b}/a_{H}$.
Plots of bubble radius 
$a_{b}/a_{H}$ against normalised radial position 
$r/r_{0}$ at several depths, for simulations with (black symbols) and without (red symbols) the Langevin model for sub-resolution velocity fluctuations: (a) 
$z=0.1$, (b) 
$z=0.3$, (c) 
$z=0.6$, and (d) 
$z=0.9$.
Returning to figure 9, both with and without the Langevin model, the effect of the free-surface interaction model is clear: as bubbles approach the free surface, they decelerate and move with the liquid, which takes them radially outwards. This behaviour is in qualitative agreement with numerical simulations of Cloete, Olsen & Skjetne (Reference Cloete, Olsen and Skjetne2009) and Pan et al. (Reference Pan, Johansen, Olsen, Reed and Sætran2021), although in both cases, these works relied on resolution of the motion of the gas flow above the free surface to capture the bubble–free-surface interactions.
5. Breaking waves
We now use our numerical framework to simulate bubble entrainment in a breaking wave. We consider a periodic third-order Stokes wave, which has been extensively studied in the literature with multi-phase mesh-based schemes (Chen et al. Reference Chen, Kharif, Zaleski and Li1999; Iafrati Reference Iafrati2009; Deike et al. Reference Deike, Popinet and Melville2015, Reference Deike, Melville and Popinet2016; Chan et al. Reference Chan, Lozano-Durán and Moin2020, Reference Chan, Johnson, Moin and Urzay2021b; Mostert et al. Reference Mostert, Popinet and Deike2022). We take the wavelength as the integral length scale, and the Froude speed as the characteristic velocity scale. The domain is a unit cube, periodic in the lateral directions, with a free-slip condition at the lower boundary. We set 
$Re=4\times {10}^{4}$, 
$Fr=1$, 
$We=1.98\times {10}^{4}$ and 
$\beta =1000/1.2$. These parameters correspond to the configuration of Mostert et al. (Reference Mostert, Popinet and Deike2022) (although we note here that with our scaling, wave overturning occurs shortly after 
$t=1$, as in Chen et al. (Reference Chen, Kharif, Zaleski and Li1999) and Chan et al. (Reference Chan, Johnson, Moin and Urzay2021b)), and give Bond number 
$Bo=We\,(\beta -1)/(4{\rm \pi} ^{2}\beta )=500$. Furthermore, we highlight that the Weber number in our numerical framework influences only the bubbles, as we do not include a surface tension model for the resolved liquid free surface. The location of the initial free surface is given by
\begin{equation} \eta(x)=\frac{1}{2{\rm \pi}}\left\{\chi\cos\left(2{\rm \pi}{x}\right)+\frac{1}{2}\,\chi^{2}\cos\left(4{\rm \pi}{x}\right)+\frac{3}{8}\,\chi^{3}\cos\left(6{\rm \pi}{x}\right)\right\}, \end{equation}
where 
$\chi$ is the initial wave steepness. The velocity within the liquid is given by
$$\begin{gather} u\left(x,y,t=0\right)=\frac{1}{2{\rm \pi}}\,\chi\sqrt{1+\chi^{2}}\cos\left(2{\rm \pi}{x}\right)\exp\left(2{\rm \pi}{y}\right), \end{gather}$$
$$\begin{gather}v\left(x,y,t=0\right)=\frac{1}{2{\rm \pi}}\,\chi\sqrt{1+\chi^{2}}\sin\left(2{\rm \pi}{x}\right)\exp\left(2{\rm \pi}{y}\right), \end{gather}$$
and 
$w(t=0)=0$. In many of the results that follow, we refer to the time since impact, 
$t-t_{im}$, where 
$t_{im}$ is the time at which the wave first breaks, indicated by a sharp increase in the (spatial) mean value of 
$\varepsilon$, and confirmed visually. An estimate of the Hinze scale is 
$a_{H}=We^{-3/5}/2\approx 1.36\times {10}^{-3}$, and throughout the following, we report bubble sizes as relative to the Hinze scale.
First, we perform simulations with 
$\chi =0.55$ in the absence of bubbles, to ensure that our SPH model provides a converged solution. Figure 11(a) shows variation of the mean kinetic energy for the single-phase wave with time. We see that the evolution of the energy is approximately converged for resolutions 
$\delta {r}\le 1/300$. For the purposes of the present study, this degree of convergence in the kinetic energy is sufficient, and except where otherwise specified, we set 
$\delta {r}=1/300$ in the following. Note that in our single-phase SPH model, surface tension is neglected, although surface tension effects are included in the closure models governing bubble dynamics. This assumption is reasonable for the present case, as the integral scale Weber number is large. Accordingly, there is no physical limit imposed on the minimum droplet sizes expected during the breaking process. However, we note that the numerics of the SPH algorithm – specifically, the particle shifting technique – introduces a surface-tension-like effect into the simulation. Although we are unable to quantify this effect accurately, we note that droplets consisting of individual SPH particles are formed during the simulation at all resolutions, implying that the effective numerical surface tension of the single-phase SPH simulation is governed by resolution, for the range of resolutions explored.
(a) Variation of mean kinetic energy (over the entire liquid domain) with time, for a single-phase wave breaking event in the absence of bubbles, for various resolutions, with initial wave steepness 
$\chi =0.55$. (b) Variation of dissipation-rate-based breaking parameter 
$b$ with initial steepness 
$\chi$ for our simulations, and experimental data of Drazen, Melville & Lenain (Reference Drazen, Melville and Lenain2008) and Melville (Reference Melville1994).
We further validate the numerical framework by performing simulations for a range of initial wave steepness 
$\chi$, and evaluating the breaking parameter 
$b$, related to the dissipation rate by
\begin{equation} b=\frac{\varepsilon_{l}g}{\rho{c}^{5}}, \end{equation}
where 
$c$ is the phase velocity given by 
$c=\sqrt {g\lambda }$ as defined for experiments, with 
$g$ the acceleration due to gravity, and 
$\lambda$ the wavelength. The quantity 
$\varepsilon _{l}$ is the dissipation rate per unit length of breaking crest, which can be evaluated following Deike et al. (Reference Deike, Popinet and Melville2015) as the product of the initial wave energy and the decay rate during breaking, calculated by numerical integration over the simulation domain. Figure 11(b) shows the variation of the breaking parameter with initial wave steepness for our simulations, alongside numerical data from Melville (Reference Melville1994) and Drazen et al. (Reference Drazen, Melville and Lenain2008). The semi-empirical fit of 
$b=0.4(\chi -0.08)^{5/2}$ (Romero, Melville & Kleiss Reference Romero, Melville and Kleiss2012) is shown by the solid line. For all the breaking cases studied (
$\chi \in [0.33,0.55]$), our simulations show a close match with the experimental data and empirical fit.
5.1. Bubble size distributions and the effect of resolution
We keep 
$\chi =0.55$ and now include bubbles in our simulation. Figure 12 shows the wave at several instants during the first (dimensionless) time unit after breaking, from the side (figures 12a–e) and beneath (figures 12f–j). An animation showing the same views of the wave during the complete breaking process is available in the supplementary material available at https://doi.org/10.1017/jfm.2023.649. As the plunging breaker impacts on the surface below, bubbles are entrained in the region of impact (figures 12a, f). At the early times after impact (figures 12a–c,f–h) the flow is largely two-dimensional, with little variation in the transverse direction. Later, as the (resolved) topologically induced gyre continues to roll, a three-dimensional structure forms under the wave (figures 12d,e,i,j). In our simulation, air entrainment occurs through several mechanisms. First, bubbles are entrained on impact of the plunging breaker. Second, bubbles are entrained at the surface ahead of the wave due to the impact of spray forward from the plunging breaker; this effect is particularly clear in figure 12(i). At approximately 
$t-t_{im}=0.8$ (figures 12e,j), the topologically induced gyre collapses, and results in further entrainment. We note here a limitation of our approach. As our model for the continuum is single-phase (we do not resolve the gas above the free surface), although we predict the entrainment of a large topologically induced gyre, the mass of this entrained air is partially lost when this bubble collapses. This limitation is common to all single-phase SPH models, and also to the model of Derakhti & Kirby (Reference Derakhti and Kirby2014). Our model does predict bubble entrainment as the topologically induced gyre collapses, due to the large values of 
$\varepsilon$ as the free surfaces come together. Combining the present work with a multi-phase SPH scheme would overcome this limitation, allowing large bubbles to be resolved, whilst bubbles with 
$a_{b}<\delta {r}$ are modelled. This is an active area of research for us, but we reserve such a model for a future work.
Visualisation of the wave at several instants in time after breaking. (a–e) The wave from the side. ( f–j) A view from beneath the wave. The free surface is shown in turquoise. In (a–e), the turbulent dissipation rate 
$\varepsilon$ is shown in a blue colour scale. In all plots, the bubbles are coloured and scaled by 
$a_{b}/a_{H}$. The corresponding times are (a, f) 
$t-t_{im}=0.015$, (b,g) 
$t-t_{im}=0.09$, (c,h) 
$t-t_{im}=0.215$, (d,i) 
$t-t_{im}=0.465$, and (e,j) 
$t-t_{im}=0.840$. The SPH particles have been interpolated to a coarse regular mesh for visualisation of the free surface, and the vertical striation visible in ( f–j) is simply a numerical artefact of this grid.
Figure 13(a) shows the evolution of the total volume of entrained bubbles (normalised by the total liquid volume), for several resolutions. There is approximately linear growth volume of air entrainment (the dashed black trend line has slope 
$1$), as found in the adaptive mesh simulations of Mostert et al. (Reference Mostert, Popinet and Deike2022). We see close agreement in the total entrained volume for 
$\delta {r}\le 1/400$, suggesting that the method is converged in terms of the total entrained volume with respect to the SPH resolution. This agreement is particularly close at early times, prior to the collapse of the topologically induced gyre. At later times, there is moderate divergence in the total entrained volume, which we believe is due to the inability of our method to account for different modes of bubble entrainment (discussed further below), and the inevitable dependence of the bubble size distribution on the resolution in a model of this type.
Evolution of the bubble population: (a) variation of total entrained volume with time after impact 
$t-t_{im}$, for several resolutions; (b) time variation of the number ratio of sub- to super-Hinze-scale bubbles, for several resolutions. For (a,b), 
$\chi =0.55$.
Figure 13(b) shows the evolution of the ratio 
$N_{H-}/N_{H+}$ for several resolutions, where 
$N_{H-}$ is the total number of sub-Hinze-scale bubbles, and 
$N_{H+}$ is the number of bubbles larger than the Hinze scale. For each resolution, 
$N_{H-}/N_{H+}$ is approximately constant, though it varies between resolutions. Our framework admits only bubbles smaller than the SPH resolution, so for smaller 
$\delta {r}$, the maximum bubble size is smaller. This reduction in the available range of super-Hinze-scale bubble sizes reduces the number of super-Hinze-scale bubbles. Given that our entrainment model is based on considerations of energy and volume balance (which are roughly invariant with 
$\delta {r}$), the total volume entrained is approximately the same for all six values of 
$\delta {r}$ in figure 13(b). Hence we see a corresponding increase in sub-Hinze-scale bubbles as 
$\delta {r}$ is decreased. The reduction in entrained volume at approximately 
$t-t_{im}=0.2$ for all resolutions (visible in figure 13a) corresponds to the period shown in figures 12(b,c,g,h), when a portion of the bubbles entrained during the initial impact are ejected in the spray. A similar pattern is visible in the total bubble population evolution in the results of Mostert et al. (Reference Mostert, Popinet and Deike2022).
Figure 14 shows the bubble size distribution 
$P(a_{b}/a_{H})$ for various resolutions at two times after impact. Figure 14(a) corresponds to 
$t-t_{im}=0.09$ and figures 12(b,g), whilst figure 14(b) corresponds to 
$t-t_{im}=0.84$, and figures 12(e,j). In both plots of figure 14, the spectrum of Deane & Stokes (Reference Deane and Stokes2002) is plotted (dashed black lines), with slope 
$-3/2$ for sub-Hinze-scale bubbles, and 
$-10/3$ for super-Hinze-scale bubbles. The experimental data from Deane & Stokes (Reference Deane and Stokes2002) are also plotted (grey triangles), along with the variation of the experimental data (grey shaded region corresponds to 
$\pm$ one standard deviation). Our model yields bubble size distributions that closely match the data of Deane & Stokes (Reference Deane and Stokes2002), including the expected super- and sub-Hinze-scale slopes, over more than an order of magnitude of bubble radii. This match includes predicting accurately the change in slope of the bubble size distribution about the Hinze scale. This result is significant, given the simplicity of our entrainment and breakup models, which are based on simple energy and volume balance arguments. The shift in bubble size distribution as 
$\delta {r}$ is reduced is clear in figure 14. For smaller 
$\delta {r}$, 
$P(a_{b}/a_{H})$ drops below the 
$-10/3$ slope and decays to zero at a smaller 
$\delta {r}$, whilst the extent to which the sub-Hinze 
$-3/2$ slope is predicted increases. The match between the expected and simulated bubble size distributions is consistent across both time instants, although 
$P(a_{b}/a_{H})$ shows more noise at early times, as the bubble population provides a smaller sample then. Again, we mention that this result shows that the energetics of our simulations are converged with respect to the resolution. We also highlight the significant result that for all six resolutions tested, the model predicts the slopes of the bubble size distribution, and where the Hinze scale falls within the range of possible bubble sizes, the model predicts the Hinze scale and both the sub- and super-Hinze-scale bubble size distribution slopes. Figure 14 also highlights the limitation of this approach – that the bubble size distribution, particularly at the largest bubble scales, depends on the resolution of the SPH. This further supports our plans as discussed above to extend the SPH model to a multiphase one, to allow the full range of bubble scales to be captured.
Bubble size distribution for 
$\chi =0.55$ at several resolutions (circles): (a) shortly after impact, at 
$t-t_{im}=0.09$; (b) at 
$t-t_{im}=0.840$. The grey triangles correspond to the experimental data of Deane & Stokes (Reference Deane and Stokes2002), and the light grey shaded area shows 
$\pm$ one standard deviation of the experimental measurements. Note that the magnitude of the experimental data has been scaled to match the non-dimensionalisation of our numerical results.
We again mention a limitation of our framework, and of models of this type. Our entrainment model is relatively simple, and as such the local (in time and space) entrainment bubble size spectrum depends only on 
$\varepsilon$ and 
$\alpha$, and not on the mechanism of bubble entrainment. The three separate entrainment processes discussed earlier – by plunging breaker, spray impact, and topologically induced gyre collapse – in reality involve quite different mechanisms. Our model, however, cannot differentiate between them. This applies similarly to the entrainment models used in Eulerian–Eulerian work (Ma et al. Reference Ma, Shi and Kirby2011; Derakhti & Kirby Reference Derakhti and Kirby2014). It is this limitation that drives us to propose the coupling of a multi-phase SPH scheme, with the gas above the liquid surface resolved, to the discrete bubble model. This would enable the entrainment of large-scale bubbles to be captured more accurately, with the computationally cheaper discrete bubble model used once large bubbles have broken to 
$a_{b}\approx \delta {r}$, and for breakup, entrainment and tracking of smaller bubbles. Additionally, in the present model, the equation of motion for the bubbles becomes increasingly stiff for small bubbles, due to the drag model. Meanwhile, if no limit is imposed on the minimum bubble size, then the number of bubbles increases significantly. The value of the time step used to integrate the equation of motion for the bubbles is specific to each bubble, and for small bubbles can be more than an order of magnitude smaller than the SPH time step. Consequently, in the present model, a large population of small bubbles becomes computationally expensive to simulate. A valuable extension to the present model would be to account for the smallest bubbles with a continuum–continuum approach (as in Derakhti & Kirby Reference Derakhti and Kirby2014), with a simplified drag model. Thus we would have a scheme that resolves large bubbles, treats intermediate sized bubbles discretely (as in this work), and treats small bubbles as a continuum, allowing the complete range of bubble scales to be modelled.
In the present work, we have focused on unidirectional periodic waves. An important future development of this work is to include a piston-type wavemaker functionality. There is an extensive body of experimental work on air entrainment in wavemaker-generated waves (e.g. Lamarre & Melville Reference Lamarre and Melville1991; Drazen et al. Reference Drazen, Melville and Lenain2008), which may be used to provide further validation of our numerical framework. This would allow us to model a broader range of representative sea states, increasing the value of the model.
Figure 15 shows the wave at several times after 
$t-t_{im}=1$. At 
$t-t_{im}\approx {2}$, a large cloud of bubbles is entrained through a reverse breaking wave, visible in figure 15(a). Our simulations predict the obliquely descending eddies observed by Nadaoka, Hino & Koyano (Reference Nadaoka, Hino and Koyano1989) and Bonmarin (Reference Bonmarin1989), and modelled using a single-phase SPH scheme by Dalrymple & Rogers (Reference Dalrymple and Rogers2006), Landrini et al. (Reference Landrini, Colagrossi, Greco and Tulin2007) and Farahani & Dalrymple (Reference Farahani and Dalrymple2014). The obliquely descending eddies drag a cloud of bubbles downwards away from the surface, which persists for several time units, and is clearly visible in figures 15(b,c). The observations of bubble distribution between breaking waves predicted by our model, at both early (figure 12) and late (figure 15) times, is in good qualitative agreement with the Eulerian–Eulerian model of Derakhti & Kirby (Reference Derakhti and Kirby2014), the detailed simulations of Chan et al. (Reference Chan, Johnson, Moin and Urzay2021b), and the experimental observations of Rapp et al. (Reference Rapp, Melville and Longuet-Higgins1990).
Visualisation of the wave at several time units after breaking. The free surface is shown in turquoise, the turbulent dissipation rate 
$\varepsilon$ is shown in a blue colour scale, and the bubbles are shown coloured and scaled by 
$a_{b}/a_{H}$. The corresponding times are (a) 
$t-t_{im}=1.96$, (b) 
$t-t_{im}=3.84$, and (c) 
$t-t_{im}=5.71$.
5.2. Influence of wave steepness $\chi$
Finally, we change the wave steepness 
$\chi$, to cover both breaking and non-breaking waves. Figure 16 shows the wave profiles (coloured by dissipation rate), and bubbles (coloured by 
$a_{b}/a_{H}$), for several wave steepnesses 
$\chi \in [0.3,0.5]$. The wave with 
$\chi =0.3$ (figure 16a) does not break. For 
$\chi =0.33$ and 
$\chi =0.35$ the wave breaks by overspilling, whilst as 
$\chi$ is increased further, there is a transition to the plunging breaker studied in the previous subsection. As a limiting test, we observe that for the case where the wave does not break, no bubbles are entrained, and that for all cases where the wave (visibly) breaks, bubbles are entrained. Figure 17 shows the time evolution of the total entrained volume (normalised by the total liquid volume), for various wave steepnesses 
$\chi \in [0.33,0.55]$. For 
$\chi =0.33$ and 
$\chi =0.35$, the total entrained volume remains small, although the wave breaking and entrainment process has a long duration. As 
$\chi$ increases and the plunging breaker regime is approached, the total volume of air entrained increases significantly, whilst the peak moves earlier, to approximately 
$1.5$ time units after breaking. As in figure 13, the linear growth of the entrained volume is clearly visible in the plunging breaker regime.
Wave profiles showing dissipation rate 
$\varepsilon$ (blue colour scale) and bubbles (coloured by 
$a_{b}/a_{H}$), for various wave steepnesses 
$\chi$: (a) 
$\chi =0.3$, (b) 
$\chi =0.33$, (c) 
$\chi =0.35$, (d) 
$\chi =0.4$, (e) 
$\chi =0.45$, and ( f) 
$\chi =0.5$. For (b–f), the wave is shown at time 
$t-t_{im}=0.2$. For (a), in which the wave does not break, 
$t=2.5$.
Evolution of total entrained volume with time for various wave steepnesses 
$\chi$. In all cases, the SPH resolution is 
$\delta {r}=1/300$.
We note here the work of Lamarre & Melville (Reference Lamarre and Melville1991) and Melville (Reference Melville1994), who performed measurements on controlled deep-water breaking waves to obtain a picture of the void fraction distribution and the time evolution of integral properties of the bubble plume. Whilst the results of our simulations are not comparable with these experiments (due to the different types and scales of waves considered), we again highlight that the inclusion of a wavemaker in our numerical framework is an area of ongoing development, and will enable direct comparison with these experimental data, along with simulation of a broader range of sea states.
6. Conclusions
In recent years, the potential of SPH for simulations of free-surface flows has been demonstrated widely, but limitations in adaptivity remain in comparison to adaptive mesh methods, preventing the use of SPH high-fidelity simulations of bubbly flows. Approaches that resolve a continuous liquid phase, and include bubbles as discrete Lagrangian particles, are established in the mesh-based community, but have not been developed previously in a mesh-free framework, despite the benefits of mesh-free approaches for free-surface flows.
In this work, we have presented a numerical framework for LES of bubbly, free-surface flows. The framework employs a continuum–discrete structure, where we use SPH for the LES of the liquid phase, whilst treating bubbles as discrete Lagrangian particles, which interact with the liquid via exchanges of volume and momentum. We introduce a Langevin model for the sub-resolution velocity fluctuations, and several closure models for bubble breakup, entrainment, and free-surface interaction. The modification of a projection method to admit a small degree of compressibility provides a significant reduction in computational costs, whilst preserving smooth pressure fields obtained with incompressible SPH. Exchanges between bubbles and liquid are implemented by SPH interpolation, and the additional construction of neighbour lists required to enable this is straightforward in an SPH framework. Individual bubbles are able to be tracked over a lifetime, including motion due to turbulent structures below the resolution of the SPH scheme. Hence models for bubble breakup (or in future, deformation and oscillation), which happen over a finite time, may be constructed in integral formulations, rather than as probabilistic events based on an instantaneous flow state.
We have demonstrated the ability of our method to simulate bubble plumes and breaking waves, with quantitative agreement with previous numerical and experimental data in terms of mean flow statistics, bubble size distributions, and bubble population evolution. Despite the inclusion of bubble entrainment and breakup through simplified energy-balance models, the numerical scheme is capable of predicting accurately the Hinze scale, and the multi-slope bubble size distribution present in breaking waves, alongside the bubble population growth rate. The bubble distributions between breaking waves generated by our model compare well qualitatively with experimental data, including the generation of bubble clouds dragged downwards by obliquely descending eddies.
Our investigations have highlighted a limitation of models of this type, which is that bubble entrainment models based on turbulence dissipation rates alone cannot account for the different physical mechanisms of entrainment in different flow types. Further developments that are planned include the extension of the SPH model to multi-phase flows, to yield a framework in which large-scale bubbles are resolved, with small-scale bubbles modelled as in the present scheme. As bubbles influence the forces and loads on structures due to wave impacts, and given the strengths of SPH in computationally affordable free-surface flow simulations, the work herein offers the potential for improvements in the accuracy of predictive modelling for wave–structure interactions.
Supplementary material
A supplementary movie is available at https://doi.org/10.1017/jfm.2023.649.
Acknowledgements
We are grateful to several anonymous reviewers whose insightful comments have helped to improve the work.
Funding
We are grateful for financial support from the Leverhulme Trust, via research project grant RPG-2019-206. J.R.C.K. would like to acknowledge funding from the Royal Society via a University Research Fellowship (URFR1221290).
Declaration of interests
The authors report no conflict of interest.
Appendix. Choice and verification of LES closure model
In the development of our numerical framework, we have investigated the use of several LES closure models in our SPH scheme. In addition to the mixed-scale model of Lubin et al. (Reference Lubin, Vincent, Abadie and Caltagirone2006), we tested a standard Smagorinsky model (referred to here as SS), and the dynamic model of Germano et al. (Reference Germano, Piomelli, Moin and Cabot1991) and Lilly (Reference Lilly1992). To avoid locally large fluctuations in 
$\nu _{srs}$ from the dynamic model, we have implemented both local averaging using a Shepard filter as in (3.7) (referred to herein as DSS), and Lagrangian averaging following Meneveau, Lund & Cabot (Reference Meneveau, Lund and Cabot1996) (referred to as DSL). The SS and DSS models have been used previously for LES studies of breaking waves in a finite-volume framework by Christensen & Deigaard (Reference Christensen and Deigaard2001) and Christensen (Reference Christensen2006). Temporarily neglecting the dispersed bubble phase, we test our single-phase LES SPH scheme on the problem introduced in Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021) for weakly compressible SPH. The problem consists of a triply periodic domain with unit side length. Setting 
$Re=10^{6}$, a ramped body force is applied over the first time unit, with forcing proportional to the generalised Beltrami flow described in Antuono (Reference Antuono2020). After this, the flow evolves, with a transition to turbulence occurring within the first few time units. With 
$Re=10^{6}$, our maximum resolution (
$256^{3}$) is significantly larger than the Kolmogorov scale, so this is a challenging test for the LES closure model. We compare the results of our code with both the SPH results and a reference fifth-order finite-volume scheme, presented in Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021). Figure 18(a) shows the evolution of the kinetic energy in the domain with time, for the various LES closure models. All results were obtained with 
$\delta {r}=1/128$ to match the resolution in Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021). The SS model (solid black line) is overly dissipative immediately after 
$t=1$, resulting in significantly delayed transition (identified by the time at which the gradient steepens) compared to the other models. This observation is anticipated: SS models are known to perform poorly, being overly dissipative, in laminar and transitional flows (Derakhti & Kirby Reference Derakhti and Kirby2014). Compared with the SPH simulations of Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021) (dashed blue line), who used a static model, both the dynamic models give an earlier transition, and a post-transition dissipation rate that matches more closely the reference finite-volume simulation (referred to as FVM; solid blue line). Interestingly, there is negligible difference between the DSL and DSS models; for the present case, the choice of averaging procedure has no effect. The mixed-scale model (MSM) provides the best match with the reference solution: the post-transition decay rates (the slopes of the lines) are a close match, as highlighted by the dotted black lines, which are parallel.
(a) Comparison of variation of kinetic energy with time for our approach and the results of Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021) (blue lines) for various LES closure models: standard Smagorinsky (SS) indicated by solid black; dynamic (Germano) model with Lagrangian averaging (DSL) indicated by solid red; Germano model with Shepard averaging (DSS) indicated by dashed black; and mixed-scale model (MSM) indicated by dashed red. All results were obtained with 
$\delta {r}=1/128$. (b) Comparison of kinetic energy variation with time for our approach at various resolutions, and the results of Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021) (blue lines).
Figure 18(b) shows the variation of kinetic energy with time as the resolution is varied, from 
$64^{3}$ particles up to 
$256^{3}$, when using the MSM. In the early stages of decay, the flow is laminar, and the theoretical decay rate has characteristic time 
$Re/48{\rm \pi} ^{2}\approx 2111$. Although this decay rate is never achieved due to the numerical dissipation in all the schemes, we see that as we increase the resolution, the decay rate does reduce (flatter lines immediately after 
$t=1$). For finer resolutions, there is an increased delay to transition, and the convergence of the post-transition decay rate is apparent. We note that our results converge to a slightly larger post-transition decay rate (steeper lines) than the finite-volume scheme used in Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021), likely due to their use of resolution 
$1/128$. In any case, the post-transition decay rates yielded by our scheme with the MSM are notably closer to the high-order reference solution of Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021) than the SPH results of Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021) for all resolutions 
$\delta {r}\le 1/96$.
We note that the choice of filter length scale 
$\tilde {\varDelta }$ in SPH is not as clear as for finite-volume schemes (where the implicit filter scale is the same as the cell size). In the present work, the implicit filter scale is taken as the smoothing length 
$\tilde {\varDelta }=h$, but this choice is by no means unique (indeed, Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021) set 
$\tilde {\varDelta }$ equal to the SPH kernel standard deviation), and may contribute to discrepancies between the LES SPH and the finite-volume reference data. The appropriate specification of 
$\tilde {\varDelta }$ in SPH is an open problem.
The energy spectra were evaluated for the above simulations by interpolating the velocity field from SPH particle positions to a Cartesian grid using a variant of the local anisotropic basis function method (King, Lind & Nasar Reference King, Lind and Nasar2020), providing fourth-order consistency. Figure 19(a) shows the energy spectra (normalised by the instantaneous total energy) for the different closure models. In all cases, the spectra are evaluated when the total energy decays to 
$0.02$, as in Antuono et al. (Reference Antuono, Marrone, Di Mascio and Colagrossi2021). The spectra for models SS and MSM appear relatively similar, although MSM yields a slope closer to 
$-5/3$ (solid blue line) for a slightly greater range of wavenumbers. These similarities might be expected – the MSM is the geometric mean of the SS model and a model based on filtering the turbulent kinetic energy. The spectra for the dynamic models (DSS and DSL) are almost identical to each other, and markedly different from the SS model and MSM. They exhibit a more pronounced peak at the forcing wavenumber (despite the forcing term being switched off several dimensionless time units prior to evaluation of the spectra), and a distinctly steeper slope than 
$-5/3$. Furthermore, they have a pronounced peak at high wavenumbers. This is consistent with our observation, when using the DSL and DSS models for simulations of breaking waves, that the dissipation rate appeared noisy. Specifically focusing on the MSM, figure 19(b) shows the energy spectra for different resolutions. All resolutions have a peak in the energy spectrum at the forcing wavenumber, as expected, followed by a decay. With increasing resolution, the range of wavenumbers for which a 
$-5/3$ slope (solid blue line) is observed increases.
(a) Energy spectra for the different closure models: standard Smagorinsky (SS) indicated by solid black; dynamic (Germano) model with Lagrangian averaging (DSL) indicated by solid red; Germano model with Shepard averaging (DSS) indicated by dashed black; and mixed-scale model (MSM) indicated by dashed red. (b) Energy spectra for the MSM closure model at different resolutions. In (a,b), slope 
$-5/3$ is indicated by a solid blue line.
Finally, whilst the tests here provide some quantitative information on the relative performance of different LES closure models in our SPH framework, we note that the performance of the models may change for different flows. For simulations of breaking waves, the turbulence is patchy and transitional, non-isotropic and multiphase. Quantitative analysis of the performance of LES models under these conditions is not trivial, and remains an open problem.
Open access
