On the sea spray aerosol originated from bubble bursting jets

Depósito de Investigación de la Universidad de Sevilla https://idus.us.es/ This is an Accepted Manuscript of an article published by Cambridge University Press Journal of Fluids Mechanics 886 (2020), available at: https://doi.org/10.1017/jfm.2019.1061 Copyright 2019. Cambridge University Press. En idUS Licencia Creative Commons CC BY-NC-ND This draft was prepared using the LaTeX style file belonging to the Journal of Fluid Mechanics 1 On the sea spray aerosol originated from bubble bursting jets Francisco J. Blanco–Rodr´ıguez and J. M. Gordillo† ´ Area de Mec´anica de Fluidos, Departamento de Ingenier´ıa Aeroespacial y Mec´anica de Fluidos, Universidad de Sevilla, Avenida de los Descubrimientos s/n 41092, Sevilla, Spain (Received xx; revised xx; accepted xx) Here we provide a theoretical framework revealing that the radius Rdof the top droplet ejected from a bursting bubble of radius Rband Bo ⩽0.05 can be expressed as Rd/Rb=Kb1−(Oh/Oh0 c)1/2for Oh .Oh0 cor as Rd≈18 µ2 l/(ρlσ) for Oh &Oh0 c with the numerically fitted constants Kb≈0.2, Oh0 c≈0.03, Oh =µl/√ρlRbσ1 the Ohnesorge number, Bo =ρlg R2 b/σ the Bond number, and ρl,µland σindicating the liquid density, dynamic viscosity and interfacial tension coefficient, respectively. These predictions, which do not only have solid theoretical roots but are also much more accurate than the usual 10% rule used in the context of marine spray generation via whitecaps for Rb.1 mm, agrees very well with both experimental data and numerical simulations for the values of Oh and Bo investigated. Moreover, making use of a criterion which reveals the mechanism that controls the growth rate of capillary instabilities, we also explain here why no droplets are ejected from the tip of the fast Worthington jet for Oh &0.04. In addition, our results predict the generation of submicron-sized aerosol particles with diameters below 100 nm and velocities ∼σ/µlfor bubble radii 10 µ m .Rb.20µ m, within the range found in natural conditions and in good agreement with experiments, a fact suggesting that our study could be applied in the modeling of sea spray aerosol production. 1. Introduction The breakup of ocean waves entraps air bubbles that disintegrate within the liquid bulk and burst when they rise back to the surface, ejecting into the atmosphere droplets which, after evaporation, constitute the so called sea spray aerosol (SSA), which influences climate by scattering and absorbing radiation and also because part of the solid particles composing the SSA serve as cloud condensation nuclei and ice nuclei (Veron 2015; de Leeuw et al. 2011;Erinin et al. 2019). Hence, the understanding, description and quantification of a fundamental proccess such as the bursting of bubbles at interfaces, which produce drops with different chemical compositions (Wang et al. 2017) depending on whether drops are formed through the disintegration of the cap film (Lhuissier & Villermaux 2012) –film drops– or via the breakup of the jet ejected upwards from the base of the collapsing cavity (MacIntyre 1972) –jet drops–, will help improve the predictive capability of climate models. It is a common belief that most of the submicrometric particles composing SSA originate from film drops (Veron 2015) but, very recently, Wang et al. (2017) demonstrated that the fraction of submicrometric SSA which can be attributed to jet drops varies from 20% to 40% depending on the seawater chemical composition. The flux of SSA †Email address for correspondence: [email protected] 2F. J. Blanco–Rodr´ıguez and J. M. Gordillo from jet drops is usually calculated using, as an intermediate step, the well known 10% rule, which approximates the drop radii as Rd= 0.1Rb, with Rbthe bubble radius. This approximation is commonly used and has been employed very recently in the paper by Wang et al. (2017). However, the numerical simulations in Duchemin et al. (2002) demonstrated, almost twenty years ago, that the function Rd/Rbvaries non monotonically with Rband also that min(Rd/Rb)∼O(0.01), a fact further confirmed by the more accurate numerical results reported in Brasz et al. (2018). Then, motivated by its atmospheric and climate implications and by the remarkable findings in Wang et al. (2017), this contribution aims at improving the quantification of the SSA produced from jet drops. Here we make use of our theory in Gordillo & Rodr´ıguez-Rodr´ıguez (2019) and provide equations for Rdand for the initial velocity of the droplets ejected, Vd, as a function of Rband of the physical properties of the liquid. In doing so, we will also present a new model describing the spatio-temporal evolution of the jet that, among other things, predicts the behaviour of function Rd/Rbaround its minimum. Let us point out that our study is limited to the range of bubble diameters of interest for SSA production, for which gravitational effects can be neglected in a first approximation. The paper is structured as follows: in §2we briefly describe the numerical simulations and present some of the results obtained. Section §3is devoted to develop a model describing the spatio-temporal evolution of the jet from its inception to its breakup and to compare our predictions with experimental and numerical results. Finally, section §4 summarizes the main results obtained. 2. Numerical setup and revision of previous results With the main purpose of analyzing the ejection of droplets from the very fast Worthington jets ejected after the bursting of a bubble of radius Rb, here we have performed an exhaustive numerical study by varying both the Bond and Ohnesorge numbers, respectively defined as Bo =ρlg R2 b/σ and Oh =µl/√ρlRbσ, with ρl,µland σreferring to the liquid density, dynamic viscosity and interfacial tension coefficient, respectively [see figure 1(a)]. Here, lengths, times and pressures are made non-dimensional using as characteristic values Rb,pρlR3 b/σ and σ/Rband hence, velocities will be made nondimensional using as the characteristic value the capillary velocity, Vc=pσ/ (ρlRb). Dimensionless variables will be written throughout the text in lower-case letters to differentiate them from their dimensional counterparts -in capitals- and the subscript dwill be used to denote quantities referring to the instant the droplet detaches from the tip of the jet. Axisymmetric numerical simulations have been performed using the open-source code GERRIS (Popinet 2003;Deike et al. 2018;Brasz et al. 2018;Gordillo & Rodr´ıguez- Rodr´ıguez 2019) using as values for the density and viscosity ratios those corresponding to air bubbles in water namely, 1.2×10−3and 1.8×10−2. The computations are started by suddenly removing the cap of a static bubble at a free surface. The initial shape of the bubble, [see figure 1(a)] at T= 0, is calculated solving the Young–Laplace equation for different values of the Bond number following the procedure detailed in Lhuissier & Villermaux (2012). Simulations with different grid refinements (maximum level of refinement varying from 10 to 14) have been carried out to check the grid independence of the results. We have further validated our numerics by comparing with Brasz et al. (2018) and the GERRIS script developed here is provided as part of the Supplementary Material. Gordillo & Rodr´ıguez-Rodr´ıguez (2019) explained that bubble bursting jets are pro- On the sea spray aerosol originated from bubble bursting jets 3 ( g !  , µ  ! l , µ l R b Z R T =  V d 2R d V j t 2R j  t 2R d0 TT j  t T  T j  t T  T d 10-3 10-2 10-1 O  10-3 10-2 10-1 r   r  N O D R O P  1 "     c   /      u z ) ) v  ) v ) s)    )    z ) z  )  Figure 1. (a) Sketch showing the sequence of events following the bursting of a bubble at a free interface. The retraction of the rim causes a convergent flow towards the base of the cavity that favours the formation of a fast jet of initial velocity Vjet and initial radius Rjet. For sufficiently small values of the Ohnesorge numbers, a drop of radius Rdand velocity Vdis issued from the tip of the jet. (b) Numerical results for the dimensionless droplet radius at the ejection time (rd0, open markers) and at the break-up time (rd, filled markers) corresponding to Bo = 0.01 and different values of the Ohnesorge number. (c) Sketch defining the main variables used in the model developed here to predict rdand vd. 0 1 2 3 z -1 0 1 -1 0 1 -1 0 1 -1 0 1 -1 0 1 0.4 0.5 0.6 0.7 0.8 z -0.1 0 0.1 r -0.1 0 0.1 r -0.1 0 0.1 r -0.1 0 0.1 r -0.1 0 0.1 r -1 0 1 -0.1 0 0.1 r O  = .0O = . O = .O = . !  !  !  !  !  !  O = . O = . Figure 2. Shapes of the collapsing cavity for different values of the Ohnesorge number at two different instants of time. Observe that a tiny bubble is entrapped below the jet for Ohc⩽Oh ⩽Ohd, with Ohc≃0.02 and Ohd≃0.05. Here, Bo = 0.01. 4F. J. Blanco–Rodr´ıguez and J. M. Gordillo duced by a purely inertial mechanism in which capillarity transforms the initial spherical bubble into a truncated cone, with the radius of the base of such a cone decreasing linearly with λ∗, namely the wavelength of the capillary wave excited during the rim retraction process (see figure 2). Gordillo & Rodr´ıguez-Rodr´ıguez (2019) found, both theoretically and numerically, that λ∗∝Oh1/2and also that the amplitude of this wave is ∝λ∗for Oh .0.02. Then, since the radius of the truncated cone varies linearly with λ∗, it was also predicted in Gordillo & Rodr´ıguez-Rodr´ıguez (2019) that the radius of the incipient jet, rjet can be expressed as rjet ∝1−(Oh/Ohc)1/2if Oh < Ohc(Bo)≃0.02 ,(2.1) with Ohc(Bo) the value of the Ohnesorge number for which the crest of the wave reaches the vertex of the cone. Notice that Ohcdepends on the Bond number because the vertex position depends on the initial shape of the static bubble, which clearly depends on the Bond number. Gordillo & Rodr´ıguez-Rodr´ıguez (2019) also explained that, for Oh .Ohc, the initial velocity of the extremely fast and thin jet, vjet, can be calculated as the flow generated by a continuous line of sinks extending upwards from the base of the truncated cone along the axis of symmetry a distance proportional to λ∗∝Oh1/2. Moreover, the intensity of these sinks also depends on Oh1/2through the capillary pressure, which depends on the cone geometry. All these facts explain why vjet depend on Oh1/2, with the precise functional form of vjet calculated in Gordillo & Rodr´ıguez-Rodr´ıguez (2019) as: vjet ∝Oh1/2for Oh 1 and vjet ∝1−(Oh/Ohc)1/2−1/2 for Oh ≈Ohc. (2.2) The shapes of the collapsing bubbles for different values of the Ohnesorge number is represented in figure 2for two different instants of time. Observe that there is a remarkable difference between the two smaller values of Oh and Oh = 0.02: in this latter case, a tiny bubble is entrapped below the jet. This bubble is created because the crest of the capillary wave reaches the vertex of the cone for Oh =Ohcand a point below the vertex for Oh > Ohc. From equation (2.2), notice that the predicted velocity would tend to infinite for Oh =Ohc, but this singularity is clearly not realistic. Indeed, it was shown in Gordillo & Rodr´ıguez-Rodr´ıguez (2019) that, when a bubble is entrapped, the jet velocity is limited by viscous stresses and, thus, the dimensional jet velocity and jet radius are respectively given by Vjet ∝σ/µland Rjet ∝µ2 l/(ρlσ). These are constant values for fixed physical properties and represent the maximum achievable velocity and the minimum achievable radius (with slight deviations near Oh = 0.03 as it will be commented below) which, in dimensionless terms, can be expressed as rjet ∝Oh2and vjet ∝Oh−1for Oh > Ohc.(2.3) The case Oh = 0.032 in figure 2illustrate a case for which equations (2.3) accurately predicts the numerical results. However, for the largest value of the Ohnesorge number illustrated in figure 2, the amplitude of the capillary wave is very much attenuated by viscosity. Consequently, the crests of the capillary waves do not reach the vertex of the cone for Oh > Ohd, with Ohd(Bo ⩽0.01) ≃0.05. Therefore, bubbles are entrapped below the jet for values of the Ohnesorge numbers within the range Ohc(Bo)⩽Oh ⩽Ohd(Bo) with Ohc(Bo 1) ≈0.02 and Ohd(Bo 1) ≈0.05 . (2.4) Clearly, the predictions in equations (2.1) and (2.2) will deteriorate as Oh approaches Ohcand, also, there will be deviations from the predictions given in equations (2.3) for On the sea spray aerosol originated from bubble bursting jets 5 values of Oh near the boundaries limiting the bubble entrapment regime namely, Ohc and Ohd. In spite of this, the numerical results depicted in figure 1(b) verify the scalings in equations (2.1) and (2.3), confirming the non-monotonic dependence of the initial radius of the jet predicted in Gordillo & Rodr´ıguez-Rodr´ıguez (2019). Figure 1(b) also reveals that the radii of the ejected droplets, rd, which is the truly important value for applications, are far larger than their initial value, rd0. Our goal in this contribution is to provide closed expressions for rdand vdat the instant the top droplet is emitted from the tip of the jet. For that purpose, we will derive a system of differential equations that will be integrated making use of the initial values for the jet radius and the jet velocity, which can be either calculated numerically or approximated using equations (2.1), (2.2) and (2.3). However, to do so, it will be first necessary to describe the spatio-temporal evolution of the jet width and of the jet velocity upstream the jet tip. 3. Modeling the spatio-temporal evolution of the jet The time evolution of the jet tip velocity, v(t), and of the radius of curvature of the tip of the jet, b(t), can be calculated using the following mass and momentum balances (Taylor 1959), ds dt =v, 4 3πdb3 dt =π r2 j(vj−v), 4 3π b3dv dt =−π rj+π r2 j(v−vj)2−6π Oh µg µlb v, (3.1) with rj(t) = r(s(t), t) and vj(t) = u(s(t), t) the values of the radius r(z, t) and of the vertical component of the liquid velocity u(z, t) at the vertical position z=s(t) where a drop of radius b(t) is located [see figure 1(c)]. In equation (3.1), notice that the dimensionless term −π rj, results from the addition of two contributions: the interfacial tension force along the perimeter of a circle of radius rj,−2πrj, plus the force πrj associated with the capillary pressure 1/rjexerted on a circle of area π r2 j. Equations (3.1) express that the jet tip is mainly decelerated by the action of interfacial tension forces, but we have also added the small force exerted by the gas because we checked that the inclusion of the Stokes drag over a sphere of radius b(t), with µgin (3.1) the gas viscosity, slightly improves the agreement between the numerical results and the model. Clearly, the system (3.1) can only be solved once rj(t) and vj(t) are known. For that purpose, here we make use of the ballistic equations for r(z, t) and u(z, t) in the form presented in Gekle & Gordillo (2010) ∂r2 ∂t +∂ ∂z r2u= 0 ⇒D Dt ln r2=−Swith S=∂u ∂z and ∂u ∂t +u∂u ∂z = 0 ⇒Du Dt = 0, (3.2) with D/Dt ≡∂/∂t +u ∂/∂z indicating the material derivative and Sthe strain rate. The parabolic system of equations (3.2) can be solved using the method of characteristics in the spatio–temporal region s(t)⩾z⩾zb(t), with zb(t) the origin of the ballistic region. The solution can be found once the liquid velocity and the jet radius u(τ, zb(τ)) = ub(τ) and r(τ, zb(τ)) = rb(τ) are known functions of time at the spatiotemporal boundary (τ, zb(τ)). Indeed, the integration of the momentum equation in (3.2) yields that, along rays z=zb(τ)+(t−τ)ub(τ) the liquid velocity is u=ub(τ). Moreover, 6F. J. Blanco–Rodr´ıguez and J. M. Gordillo 0.5 1 1.5 2 2.5 3 0 r 0.5 1 1.5 2 2.5 3 z 0 r 0.5 1 1.5 2 2.5 3 0 r 0.5 1 1.5 2 2.5 3 0 r 0.5 1 1.5 2 2.5 3 0 r ( 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6  0 20 40 60 80 100 120 u = ,  ) 0 0.01 0.02 0.03 b 0 0.02 0.04 t ! t j   0 50 100 v 0 0.05 0.1 0.15 0.2      0.8 0.9 1 1.1 1.2 1.3 1.4  "    /  O   .0 O   .0 O   .0 O   . O   . O   .0 O   . O   .   # 2  d Figure 3. (a) Time sequence showing the comparison between the spatio-temporal evolution of the jet radius and the model prediction (blue line) calculated for Bo = 0.01 and Oh = 0.024 at the following instants of time: t= 0.463,0.476,0.489 and 0.502. The last frame in the sequence corresponds to the instant of time, t= 0.506, the jet breaks into drops and includes the shapes of the jet for instants of time close to its inception to stress the fact that the final drop diameter is much larger than the initial one. (b) Comparison between the spatio-temporal evolution of the jet velocity calculated numerically and that predicted by the model (blue line). The predicted values could be closer to the calculated ones if the pressure gradient term was added into the momentum equation (3.2), as it is shown in the Supplementary Material. (c) The different continuous lines indicate the time evolution of the dimensionless wavenumber calculated solving equations (3.1) using the analytical expressions in (3.7) for values of the Ohnesorge number around the one the minimum drop radius is attained. The breakup instants calculated using GERRIS, which are indicated using filled squares, reveal that a drop is ejected when the value of the ratio π rj/b calculated using the system (3.1), is such that π rj/b = 0.85. to integrate the continuity equation in (3.2) it is first convenient to notice that ∂ ∂z ∂u ∂t +u∂u ∂z = 0 ⇒DS Dt +S2= 0 ⇒ −DS S2=Dt ⇒S=Sb(τ) 1+(t−τ)Sb(τ),(3.3) On the sea spray aerosol originated from bubble bursting jets 7 with (see Gekle & Gordillo 2010) Sb(τ) = ∂u ∂z (τ, zb(τ)) = −dub/dτ ub(τ)−dzb/dτ .(3.4) Therefore, the substitution of equation (3.4) into the continuity equation in (3.2) yields D Dt ln r2=−Sb(τ) 1+(t−τ)Sb(τ)⇒r=rb(τ) p1+(t−τ)Sb(τ).(3.5) The inertio–capillary balance existing for values of the Ohnesorge number Oh ∼Ohc (Gordillo & Rodr´ıguez-Rodr´ıguez 2019) during a short interval of time after the jet is ejected causes the time evolution of the jet to be self–similar (Zeff et al. 2000;Brasz et al. 2018;Lai et al. 2018), a fact yielding that the liquid velocity and the jet radius at zb(τ) can be respectively expressed as ub(τ) = u0τ−1/3and rb(τ) = r0τ2/3. Defining the origin of the jet region at the vertical position zb(τ) = u0τ2/3, making use of equation (3.4) and of the fact that ub(τ) = u0τ−1/3and that z=zb(τ) + ub(τ)(t−τ), yields Sb(τ) = 1/τ and τ= (u0t/z)3.(3.6) The substitution of equations (3.6) into u=u0τ−1/3and into equation (3.5) yields u(z, t) = z+z0 t+t0 , r(z, t) = δ(t+t0)3 (z+z0)7/2,(3.7) with δ=r0u7/2 0and where t0and z0account for the fact that both the time and the vertical position in the ballistic equations are shifted with respect to the respective origins in the numerical simulations. The values of t0,z0and δin equations (3.7) are determined making use of the numerical values of vj,rjand of the strain rate, S, calculated from the numerical simulations right after the jet is ejected, [see figure 1]. The parameter δ, which measures the ratio between inertial and interfacial tension forces always satisfies the condition δ1, this being the reason why the capillary pressure gradient could be neglected in the momentum equation (3.2) (see the Supplementary Material for details of the solution including the capillary pressure gradient term). Indeed, in spite of its simplicity, the model expressed by equations (3.1), solved using the analytical expressions in (3.7), correctly reproduces the time evolutions of the radius and vertical liquid velocity within the jet, r(z, t) and u(z, t), in a finite spatial region located upstream the drop, which suffices to correctly evaluate the mass and momentum fluxes entering the drop, see figures 3(a)-(b). Let us point out that the very good agreement depicted in figures 3(a)-(b) between the jet shapes calculated numerically using GERRIS and our analytical prediction in (3.7), could have been extended to a larger spatial region. This could have been done providing correct boundary conditions to the parabolic set of equations (3.2) during a longer time period. However, it is not our goal here to predict the full shapes of the Worthington jets, but the diameters and the velocities of the droplets ejected. In this regard, the insets in figure 3(b) show that the time evolutions of the radius and the velocity of the tip of the jet predicted solving equations (3.1) and (3.7) (blue lines) closely follow the time evolutions calculated numerically using GERRIS (black lines). We checked that the good agreement shown in figure 3(a)-(b) between predictions and numerical results extends to the whole range of values of Oh investigated here. Since the theoretical model presented here is able to correctly capture the time evolutions of both v(t) and b(t) (see the inset in figure 3(b)) the only missing ingredient in order to predict the drop radius and drop velocity is to determine the instant of time the jet tip breaks into drops. 8F. J. Blanco–Rodr´ıguez and J. M. Gordillo 3.1. Jet breakup condition The radius and the velocity of the first droplet ejected is calculated here particularizing the functions b(t) and v(t) (see the inset in figure 3(b)) at the instant the jet tip breaks as a consequence of the growth of a capillary instability. Since the wavelength of the capillary perturbation can be approximated by the diameter of the droplet, as it is suggested by the inset of figure 3(c), our criterion for capillary breakup can be expressed as the instant of time for which the dimensionless wavenumber, κ=πrj(t)/b(t), verifies the condition πrj/b =κc, with κc<1 (Rayleigh 1878) a constant. We do not predict here the value of this free constant κc, but fix it here by comparing the model predictions with the numerical results. Indeed, figure 3(c) represents, with continuous lines, the time evolutions of the ratio πrj(t)/b(t) obtained solving the model equations (3.1) using the analytical expressions in (3.7) for several values of Oh. In each of these curves, we include a filled square corresponding to the jet breakup time calculated using GERRIS and, in all cases considered, we find that the breakup instant predicted using the GERRIS values takes place when the ratio πrj/b calculated using our model verifies the condition πrj/b =κc, with κc= 0.85. This fact confirms our hypothesis that the breakup instant can be characterized through the value of the ratio πrj/b. Notice that the value for κc found here is slightly larger than the dimensionless wavenumber of fastest growth rate of capillary perturbation calculated by Rayleigh (1878) for the case of cylindrical jets, κc≃0.7. However, this is not surprising in our case because Gordillo & Gekle (2010) reported that the wavelengths of capillary perturbations leading to the end-pinching of stretched Worthington jets is close to the jet perimeter, which implies values of the dimensionless wavenumber of fastest growth rate larger than 0.7. The insets in figure 4show the comparison between the values of rdand vdcalculated numerically using GERRIS (full symbols) and those predicted by our model (open symbols), the latter obtained solving equations (3.1) and applying the jet breakup condition πrj/b = 0.85. The inset in figure 4(a) shows that most of the predicted values overlap the ones calculated using GERRIS. The good agreement between predictions and numerical results in the complex region where both rdand vdbehave non monotonically with Oh, validates our approach. 3.2. Scaling the size and velocity of the primary droplet While the model developed in the previous sections is capable predict, with great accuracy, the diameters and velocities of the first drop ejected, it is of little practical use since it requires to solve a system of ordinary differential equations. The purpose here is to derive approximate algebraic equations for both rdand vdmaking use of our recent result in Gordillo & Rodr´ıguez-Rodr´ıguez (2019), where we derived expressions for both the radius and the velocity of the incipient jet, rjet and vjet respectively (see figure 1a), as a function of the Ohnesorge number. The model represented by equations (3.1) and (3.7) and, most importantly, the breakup criterion deduced here, πrj/b = 0.85, are free of parameters, exceptuating the small air drag term in equation (3.1) -which, in fact, we could have neglected. Since the system of differential equations only receives information of the bulk flow through the initial conditions, which depend on Oh1/2for Oh < Ohc, we assume that rd=rjet +f(Oh1/2)≃rjet +C1+C2Oh1/2+O(Oh),(3.8) where we have made use of the fact that Oh 1. Hence, making use of equations (2.1)