On the physics of transient ejection from bubble bursting

J. Fluid Mech. (2021), vol.929, A12, doi:10.1017/jfm.2021.791 On the physics of transient ejection from bubble bursting Alfonso M. Gañán-Calvo1,2,†and José M. López-Herrera1 1Dept. Ingeniería Aerospacial y Mecánica de Fluidos, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092, Spain 2Laboratory of Engineering for Energy and Environmental Sustainability, Universidad de Sevilla, E-41092 Sevilla, Spain (Received 28 March 2021; revised 26 July 2021; accepted 6 September 2021) Using a dynamical scaling analysis of the flow variables and their evolution due to bubble bursting, here we predict the size and speed of ejected droplets for the whole range of experimental Ohnesorge and Bond numbers where ejection occurs. The transient ejection, which requires the backfire of a vortex ring inside the liquid to preserve physical symmetry, shows a delicate balance between inertia, surface tension and viscous forces around a critical Ohnesorge number, akin to an apparent singularity. Like in other natural phenomena, this balance makes the process extremely sensitive to initial conditions. Our model generalizes or displaces other recently proposed ones, impacting on, for instance, the statistical description of sea spray. Key words: bubble dynamics, aerosols/atomization, capillary flows 1. Introduction Everyday experience teaches that radially convergent flows close to a liquid surface produce vigorous transient liquid ejections in the form of a jet perpendicular to the surface, as those seen after bubble bursting (Kientzler et al. 1954), droplet impact on a liquid pool (Yarin 2006)or cavity collapse (Ismail et al. 2018). The mechanical energy of the flow comes from the free surface energy of the initial cavity. Among these processes, bubble bursting can be considered the parametrically simplest and most ubiquitous one: just two non-dimensional parameters, namely the Ohnesorge (Oh)and Bond (Bo) numbers, determine the physics and the outcome (Oh and Bo, weighting the viscous and gravity forces with surface tension, respectively). It has received special attention in the scientific literature due to its impact at global planetary scales (bubble bursting at the sea †Email address for correspondence: [email protected] © The Author(s), 2021. Published by Cambridge University Press. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons. org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited. 929 A12-1 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 A.M. Gañán-Calvo and J.M. López-Herrera surface, e.g. Kientzler et al. (1954), Blanchard & Woodcock (1957),Veron(2015)and Sampath et al. (2019); disease and pandemic transmission, e.g. Bourouiba (2021)) and in more ‘indulging’ applications (e.g. Ghabache et al. 2016; Séon & Liger-Belair 2017). Current works state that the subject is already deeply understood (Berny et al. 2020; Sanjay, Lohse & Jalaal 2021),yet despite optimistic statements, even in the simplest case where the liquid properties are constant and the gas-to-liquid density and viscosity ratios are very small, unsettling but crucial questions remain open. Within an ample range of Oh and Bo, bubble bursting exhibits a strong focusing effect due to the nearly cylindrical collapse of a main wave at the bottom of the parent bubble. In these cases, a tiny bubble gets trapped and an extremely thin and rapid spout emerges. That initial spout grows into a much taller, thicker and slower transient jet that may eventually expel droplets. Their eventual size and speed is determined by Oh and Bo. However, for an apparently critical Oh,(Oh∗), the initial ballistic spout keeps its large speed and extremely thin radius until it ejects a nearly invisible droplet, such that the spout momentum (proportional to its velocity times its volume) seems to vanish. Strikingly, a nearly symmetric behaviour is observed around that Oh∗(Séon & Liger-Belair (2017), figure 16). Around that special parametrical region, experiments and numerical simulations exhibit a strong dispersion attributed to local viscous effects and interaction with the gas environment (Brasz et al. 2018; Dasouqi, Yeom & Murphy 2021). In general, the global streamlines of the flow resemble a dipole with its axis perpendicular to the interface (figure 1), which is responsible for the vigorous liquid–air transfer process characteristic of this phenomenon (Lee et al. 2011). The energy excess of the convergent flow leads to both a fast transient capillary jet and avortex ring in opposite directions but equivalent effective momenta by virtue of Newton’s third law of motion. Their radically different kinematics, due to a large density disparity across the interface, is not an obstacle to their profound symmetry, as we will show here. We claim that this symmetry and its scalings are the keys to determine the eventual ejected droplet size and speed, and the appearance of singular dynamics. In this regard, some of the above cited authors and many others (e.g. MacIntyre 1972; Duchemin et al. 2002; Thoroddsen et al. 2018) have noticed the trapping of tiny bubbles after bursting. In particular, the one trapped at the bottom of the formed cavity is the relic of a singularity produced by the collapse of a pilot capillary wave from the rim formed after the upper film rupture (Duchemin et al. 2002; Thoroddsen et al. 2018). This singularity would explain why the range of Oh numbers where this occurs exhibits the largest ejection speeds and smaller drops. In this work, we present: (i) a universal non-trivial scaling of the evolution of the flow kinematic variables (characteristic lengths and velocity); (ii) a prediction of the eventual ejected droplet size and speed, for the complete range of Oh and Bo for which ejection occurs (which generalizes Gañán-Calvo (2017)); and (iii) an explanation why an apparently singular value of Oh∗appears. The existing data from prior works (e.g. those collected in Gañán-Calvo (2017); Brasz et al. (2018),Deikeet al. (2018), Gañán-Calvo (2018)and Berny et al. (2020)) and novel numerical simulations are compared with our theoretical results. When the droplet size and ejection velocity are properly scaled according to our model, the data show a remarkable collapse around our model except in the vicinity of Oh∗. A careful inspection reveals that all data from numerical simulations exhibit an enormous degree of dispersion that experimental data do not display. In fact, while numerical series should be independently fitted, the available experimental series appear to obey a single universal fitting to our model. Given that it relies on just two non-dimensional parameters alone, Oh and Bo, this has a twofold consequence: (i) it highlights the critical importance of the artificial initial conditions after film bursting in the numerical simulations among the different authors, which precludes universality; and (ii) it provides a proof that the natural 929 A12-2 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 Transient ejection from bubble bursting t < t0 t = ~ t0 t0 – t = 3.54 × 10–3 tc t0 – t = 8.4 × 10–6 tc t0 – t = 4×10 –4 tc t < t0 (t0 – t = 3.54 × 10–3 tc) t < t0 (t0 – t = 8.4 × 10–6 tc) t > t0 (t – t0 = 4×10 –4 tc) t = ~ t0 (t0 – t = 10–6 tc) t > t0 –2.4 –0.6–0.4–0.2 0 0.2 0.4 0.6 –2.2 –2.0 –1.8 –1.6 –1.4 –1.2 –1.0 –0.8 –1.80 –0.08 –0.04 0 0.04 0.08 Trapped microbubble –1.75 –1.70 –1.65 –1.60 –1.55 –1.50 –10 –5 0 5 10 –2 –1 0 1 2 –2 –1 0 1 2 –2 –1 0 1 2 –1.80 –0.08 –0.04 0 0.04 0.08 –1.75 –0.004 –0.003 –0.002 –0.001 –0.001–0.002 0.002 0.0020.001 0.001 –1.70 –1.65 –1.60 –1.55 –1.50 W B B V R L W W L A B (a) (b) Figure 1. General overview of the flow development around the critical time t0of collapse of the main pilot wave at the bottom of the cavity, for Oh =0.032, Bo =0. The three instants here illustrated are t0−t=3.54 × 10−3tc,t0−t=8.4×10−6tcand t−t0=4×10−4tc(with tc=(ρR3 o/σ)1/2).(a) Global flow streamlines (similar to dipole contours at a liquid–gas surface) showing a particular one (line A–B) ending at the point (B) just above where the collapse eventually occurs. The blue sphere of the three-dimensional rendering is the initial bubble immediately after bursting. (b) Local details of the same instants. The stream function levels shown are closer around the tiny trapped bubble to exhibit the vortex ring. Note that the upper axial point of the ellipsoid defining the vortex ring coincides with the point where the surface collapsed at t=t0and remains at that position: observe the horizontal line connecting the three subpanels of panel (b). The main flow velocities W(radial) and V(axial, ejection) are indicated. The dashed lines represent the instantaneous stream tube which feeds the collapsing point or the issuing jet. In these simulations using the volume of fluid method (VOF)(Basilisk, Popinet 2015), the density and viscosity of the liquid is 1000 and 100 times that of the gas, respectively, to reflect similar relations to the air–salt water ones. process should be fundamentally universal in terms of Oh and Bo alone, as expected from the universality of initial film drainage and rupture (Eggers & Fontelos (2015)and references therein), which supports the universal validity of our model once the three proposed physically relevant fitting parameters are found. 2. Formulation and dynamical scaling analysis Consider a gas bubble of initial radius Rovery close or tangent to the free surface of a liquid with density, surface tension and viscosity ρ,σand μ, respectively, in static 929 A12-3 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 A.M. Gañán-Calvo and J.M. López-Herrera equilibrium under a relatively small action of gravitational acceleration gnormal to the free surface far from the bubble. The flow properties are considered constant. At a certain instant, the thin film at the point of tangency breaks and the process of bubble bursting starts (figure 1a). The problem is thus determined by the Ohnesorge (Oh = μ/(ρσ Ro)1/2) and Bond (Bo =ρgRo2/σ) numbers alone (Montanero & Gañán-Calvo 2020). Alternatively, instead of Oh, some authors prefer to use the Laplace number (La), which is simply La =Oh−2(e.g. Duchemin et al. 2002; Lai, Eggers & Deike 2018). Also, the Morton number, Mo =Bo Oh4, has been used as a relative measure of the gravitational forces (e.g. Berny et al. 2020). The liquid rim that is initially formed pilots two main capillary wavefronts, one that advances along the bubble surface towards its bottom, and the other propagating away from the cavity. In some parametrical regions of the domain (Oh,Bo) of interest, a significant fraction of the energy content in the nonlinear wave spectrum of the initial surface configuration (after the bubble film burst) is in the small-wavelength domain (wavelets). For small Oh, these wavelets may arrive first to the bottom, but they do not produce any significant effect compared with the collapsing main wave (Gañán-Calvo 2018) in the parametrical realm considered here. When the main wavefront approaches the bottom, it becomes steep. If the wavefront leads to a nearly cylindrical collapsing neck in a certain region, a tiny bubble gets trapped below the surface once the ejection process commences. The process of radial collapse and bubble trapping is locally described by the theory of Eggers et al. (2007), Fontelos, Snoeijer & Eggers (2011)and Eggers & Fontelos (2015) when a strong asymmetry in the axial direction occurs. However, the asymmetry of bubble trapping here observed leads to an extremely rapid and thin initial ejection in the opposite direction from the trapped bubble. Figure 1(b) shows three illustrating instants of the flow development around the critical instant t0at which the axial speed of the liquid surface on the axis reaches its maximum (immediately after the tiny gas neck collapses). When t<t0, the flow is predominantly radial (spherical) with a characteristic speed Wat the interface. After collapse (t>t0), while the main flow keeps running radially with speed W, the axial asymmetry of the flow and the kinetic energy excess of the liquid is axially diverted in the two opposite directions, with radically different results: (i) towards the open gas volume, producing the liquid spout with speed Vand characteristic radial length R; and (ii) towards the liquid bulk, as a backward reaction vortex ring (figure 1b,incolours) with characteristic length L, around the trapped microbubble. Eventually, the advancing front of the resulting capillary jet expels a droplet or droplet train with characteristic radius R. Obviously, all scales (W, V, R, L) are time-dependent along their evolution. Some relevant prior works (e.g. Zeff et al. 2000;Laiet al. 2018)haveanalysed the time self-similarity of the flow variables. Assuming that those time self-similarities exist, our focus is here to obtain closed relationships among those scales and to predict the size and speed of ejections reflected by the eventual values of Rand Vat the end of the process. This approach, though, demands the identification of the key symmetries appearing in the problem that may lead to an effective problem closure (Gañán-Calvo, Rebollo-Muñoz & Montanero 2013). To this end, a set of relations among the radial and axial characteristic lengths and velocities was formulated in Gañán-Calvo (2017,2018) on the basis that inertia, surface tension and viscous forces should be comparable very close to the instant of collapse of the free surface. Under the light of detailed numerical simulations, the previously proposed simplified models are put into perspective after a rigorous derivation of the new model. Following a similar rationale to that in Gañán-Calvo (2017), this derivation is made in two steps: (i) we obtain explicit relationships among the four time-dependent variables along the process, verifying their validity numerically; 929 A12-4 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 Transient ejection from bubble bursting and (ii) we close the problem identifying the critical key factors on the energy balance of the process. The final model is compared with the ample set of experimental and numerical available data, and its validity is extended to the whole range of Oh numbers where ejection takes place. 2.1. Instantaneous scalings of kinematic variables along the process The general momentum equation of the liquid can be written as ρvt+∇(ρv2/2+p−Pa+ρgz)=ρv∧∇∧v+μ∇2v,(2.1) where vis the velocity vector, subscript tdenotes a partial derivative with time, pis the liquid pressure, zthe axial coordinate and Pathe static gas pressure. Equation (2.1) can be multiplied by the unit vector ltangent to any instantaneous streamline, in particular the streamline flowing through a point Awhere it meets the free surface to a point Bat the vicinity of the point of collapse (see figure 1a). Integrating with respect to the streamline coordinate sfrom Ato Byields ρB A l·vtds+ρv2/2B+σ∇·n|B+ρgz|B A=μB A l·∇2vds,(2.2) since the velocity is negligible and pressure is Paat the point A. Here, nis the unit normal on the liquid surface, and zis the depth of point Brespect to A.Asageneral consideration, the liquid velocities are very small everywhere compared with those at distances Lto the collapsing region, which may exhibit a self-similar flow structure (Zeff et al. 2000; Duchemin et al. 2002;Laiet al. 2018). The length scale Lalso characterizes the inverse of the mean local curvature of the liquid surface around the region of collapse, for any given time t. Thus, Lobviously changes with time around the instant of collapse. To look for symmetries around t0, let us consider two situations, one for t<t0and the other for t>t0such that their characteristic length scales Lare the same (see figure 1). Then, one may estimate the characteristic values of each term of (2.2) for both t<t0and t>t0. In any case, the first two terms of (2.2) are always of the same order for this unsteady problem. For t<t0,boththeleft-hand integral and the kinetic energy term at Bin (2.2) should scale as ρW2. The surface tension term at Bshould be proportional to σ/L, and the gravity term to ρgRo. Finally, the curvature of streamlines suggests that the viscous stresses should scale as μW/L, which fixes the scale of the right-hand side term of (2.2). This is so because boundary layers cannot be thinner than the scale Lof the collapsing main wave, as shown by Gañán-Calvo (2018): due to the near-zero tangential viscous stress at the free surface, boundary layers cannot be maintained by other means than the internal inertia. In effect, at the vicinity of the surface, there are no other sources of momentum and no other concentrated sources for viscous diffusion than the internal flow, whose fate is dictated by the free surface geometry and,consequently, its energy content. Since that characteristic length Lis fixed by the collapsing main wave, this should be the minimum (if any) characteristic length of viscous diffusion as well. Thus, the overall scaling balance of (2.2) can be expressed as ρ(W2+β1gR o)∼σL−1+α1μWL −1,(2.3) where prefactors α1and β1should be universal constants reflecting the dominance or subdominance of the corresponding terms. Using the natural length and velocity, 929 A12-5 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 A.M. Gañán-Calvo and J.M. López-Herrera μ=μ2/(ρσ ) and vμ=σ/μ, respectively, and defining ζ=L/μ,ω=W/vμ=and 1=β1Oh2Bo, this equation can be written as ω2+1∼(1+α1ω)ζ−1.(2.4) When t>t0, the vigorous ejection in the axial direction generates the new characteristic velocity Vand radial size Rscales of the liquid jet. The physical symmetry around to should keep the scale of the radial velocity Win the bulk and the spatial scale L(including the axial scale of the jet) the same: see the central subpanel in figure 1(b)forOh =0.032, Bo =0. Even before collapse, the flow develops a stagnant region in the liquid just below the collapsing point. From this point, a vortex ring of characteristic (growing) size Land speed Wlike those beautifully described by Maxworthy (1972) develops;see figure 4 of his work (dissecting absolute and relative speeds is not trivial in these flows, as Maxworthy showed). By virtue of total energy preservation and assuming that the trapped microbubble is much smaller than the vortex, the mechanical energy of the vortex, ρW2L3, should mirror that of the vigorous microjet, ρV2R2L. Furthermore, the streamlines below the jet form a feeding stream tube coming from deep in the liquid domain (surrounding the vortex ring, see figure 1b,right-hand subpanel), which should ultimately originate far away from the bubble at the free surface. In contrast, Gordillo & Rodriguez-Rodriguez (2019) ignores such a vortex and proposes a simplified radial (cylindrical) flow model, which hampers the correct identification of the local flow scales. In summary, one can write ρV2R2L∼ρW2L3=⇒ VR∼WL=⇒ υχ ∼ωζ, (2.5) where χ=R/μand υ=V/vμ. As a secondary consequence, the incoming conical flow raises the conical surface at characteristic speed W. One can now estimate the scaling of the different terms of (2.2) for a streamline ending at a point Bat the surface of the jet (see figure 1a), and write the following balance: ρ(V2+β2gR o)∼σR−1+α2μWL −1,(2.6) where, again, prefactors α2and β2should be universal constants at the end of the process. Using natural scales and defining 2=β2Oh2Bo,(2.6)reads υ2+2∼χ−1+α2ωζ−1.(2.7) As in Gañán-Calvo (2017), (2.4), (2.5)and(2.7) define three relationships between the time dependent scales (χ,υ,ζ,ω). In effect, an explicit relationship between ωand ζis easily derived after rewriting (2.4): (ωζ )2−α1(ωζ) +ζ(1ζ−1)∼0⇒ωζ ∼ϕ≡α1 2+α1 22 +ζ(1−1ζ)1/2 . (2.8) Also, introducing (2.5)and(2.8)in(2.7), the following relationship between υand ζ results: (ϕυ)2−(ϕυ) +ϕ2(2−α2ϕζ−2)∼0⇒ϕυ ∼φ≡1 2+1 4+ϕ2α2ϕ ζ2−21/2 . (2.9) 929 A12-6 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 Transient ejection from bubble bursting 5 10 50 100 0.1 12 2 5 10 20 0.034 0.033 0.032 0.031 0.030 0.029 0.028 0.026 0.035 0.036 0.037 0.038 51020 50 0.1 0.5 1.0 5.0 10.0 0.1 0.5 1.0 5.0 10.0 0.2 0.5 1.0 2.0 5.0 ζvv χωχ (a)(b)(c) Figure 2. Numerical data (momentum-conserving VOF, Basilisk, minimum cell size 3 ×10−4Rofor this figure) compared with prediction (2.10a–c). The Oh numbers of the different series are given at the right-hand side of the figure. (a)Plot of χversus ζ. Black dashed line (χ∼ζ) shows the essential proportionality between both spatial variables along the whole process as predicted, consistently with prior analyses (Zeff et al. 2000). (b)Plot of υversus ω. Black dashed line (υ∼ω) also shows the proportionality between both variables for υbelow 20; the deviation above that value reflects the very high initial velocity acquired by the spout front at its inception for Oh around Oh∗(focusing effect). (c)Plot of χversus υ. The black dashed line is the theoretical prediction for t−to>tμ(α1,2→0), which yields χ∼υ−2. The blue dashed line would correspond to 0 <t−to<tμ,withα1,21, yielding χ∼υ−1. Observe the visible change of regime in the case Oh =0.038 around χ≃0.65, i.e. around R≃0.65μ. The scaling law proposed by Gañán-Calvo (2017), χ∼υ−5/3, is also plotted as a red dashed line. In summary, the following explicit relationships hold: χ∼ϕ2φ−1,υ∼ϕ−1φ, ω ∼ϕ/ζ, (2.10a–c) which can be tested against numerical simulations. For α1→0, α21and1,2=0, (2.10a–c) yields χ∼υ−5/3,ζ∼υ−4/3and ω∼υ2/3, as predicted by Gañán-Calvo (2017). However, given that the ejection is fundamentally a ballistic (inertial) process, one may expect two possible regimes in the evolution of the ejection. (i) When |t−t0|<tμ=μ3/(ρσ2), the inertia and viscous forces should be dominant (α1,21). In this case, for Bo =0, χ∼ζ∼υ−1∼ω−1. (ii) For |t−to|>tμ, inertia and surface tension forces should prevail (α1,21). Here, if Bo =0, one obtains χ∼ζ∼υ−2∼ω−2, consistently with Zeff et al. (2000). Observe that χ∼ζand υ∼ωalways hold. Testing the validity of (2.10a–c) involves showing that the two limits (i) and (ii) and their corresponding scaling relationships are indeed reached as the time-dependent characteristic scales (χ,ζ,υ,ω) evolve along the process. To do so, we analyse their evolution from numerical simulation (momentum-conserving VOF, Basilisk code). We have recorded the time-dependent values of the maximum radius of curvature in the meridional plane at the point of maximum radial velocity of the interface, and the radius of curvature and velocity of the jet front on the jet axis to represent the time-dependent characteristic values of L,Rand V, respectively, for t>t0(no scales of Rand Vappear before ejection, i.e. for t<t0). Since (2.10a–c) does not provide the time dependencies of the variables, in figure 2(a–c) we plot the non-dimensional evolving values of χversus ζ, ωversus υand χversus υ, respectively, for a detailed range of Oh numbers. Firstly, figure 2(a) demonstrates the fundamental proportionality between ζand χalong the whole process, as anticipated. In this regard, the small constant of proportionality between both scales, χ≃0.03ζ,should be noticed. The importance of this observation will be apparent after the global assessment of the energy balance of the process in §2.2. 929 A12-7 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 A.M. Gañán-Calvo and J.M. López-Herrera D2 D1 Ω(t) Figure 3. The control volume Ω(t)bounded by the free surface D1and the fixed hemispherical surface D2, considered for the integral representation of the mechanical energy balance, close to t0. Secondly, the remarkable collapse of the curves of χversus υin figure 2(c) reveals the two clear regimes as anticipated: viscosity-dominated (α1,21) when χ1(R μ) and surface tension-dominated (α1,2→0) when χ1(Rμ), corresponding to χ∼υ−1and χ∼υ−2, respectively. The robustness of the results is thus confirmed: the formation and pinch-off of a droplet is expected to take place once surface tension becomes important, resulting in a dependency as observed. The relatively small deviations from the theoretical prediction (2.10a–c) for the smallest and largest values of υin figure 2(c) reflect the vibrations of the droplet about pinch-off and ejection (not considered here) and the very initial evolution just after t=t0, respectively. 2.2. Energy balance ForagivenBo that fixes the initial static geometry of the bubble, droplet ejection is prevented by the loss of mechanical energy above a certain Oh threshold value. In the absence of gravitational forces (Bo ≃0), a reported value of that threshold is Ohth = 0.052 (Lee et al. 2011). Below that threshold, the finest and fastest jets and droplets are observed for a certain value Oh∗<Ohth. Therefore, to close the problem, a key question which is not resolved by (2.10a–c) is the mechanical energy available to perform the ballistic ejection close to to. One can express the mechanical energy equation in a sufficiently ample fluid volume Ω(t)around the initial bubble from the instant of bubble bursting (t=0) up to the point of collapse to. The fluid volume Ω(t)can be formed by the free surface and a hemisphere around the cavity with a radius approximately twice or three times larger than Ro(see figure 3): Ω(t) ρv2/2+gzdΩ t0 t=0 =t0 t=0S(t) v·τ−pI·ndAdt,(2.11) where τis the viscous stress and Ithe identity matrix. The surface of the integral domain is S(t)=D1(t)+D2,where D2is a hemisphere with a fixed radius. We define the capillary velocity as Vo=(σ/ρRo)1/2, which is the global expected value of the average velocities along the process. Observing the geometric configuration of streamlines in figure 1, the most relevant feature is a stream tube which ends up at the collapse point, coming from below (figure 1b– 929 A12-8 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 Transient ejection from bubble bursting subsequently also shown in figure 6). This stream tube, with velocity Wand cross-sectional area proportional to L2, would ultimately feed the ejection. The global view of the streamlines’ configuration in figure 1(a) reveals that the characteristic length of that stream tube is indeed Ro, since its configuration remains practically unaltered along the process sufficiently far from the collapse region. Thus, the kinetic energy term should be proportional to Ω(t) (ρv2/2)dΩ t t=0 ∼ρW2×L2Ro.(2.12) Note that the initial kinetic energy is zero. If one neglects viscous and gravity works, expression (2.11) suggests the scaling ρW2L2Ro∼σR2 o∼ρV2 oR3 o=const. ⇒WL ∼VoRo,(2.13) which together with the scaling previously obtained as WL ∼VR after collapse, it would suggest VR ∼VoRo. This relationship was already observed by Duchemin et al. (2002) and Ghabache et al. (2014), among others, for negligible viscosity and gravity. However, it does not faithfully describe the nonlinear behaviour observed for Oh around Oh∗:note that in the vicinity of Oh∗the product VR becomes much smaller than VoRo(Duchemin et al. 2002; Walls, Henaux & Bird 2015; Ghabache & Séon 2016; Gañán-Calvo 2017;Séon & Liger-Belair 2017;Deikeet al. 2018). There is not a general consensus around the value of Oh∗and the origin of this behaviour. Naturally, Oh∗depends on the Bond number through the initial shape of the bubble before bursting, since that shape determines the detailed dynamic path of the process. Reducing that geometry to a sphere, i.e. Bo ≃0, the simulations presented here and those of Berny et al. (2020) and Deike et al. (2018) agree on Oh∗≃0.03; in contrast, Gañán-Calvo (2017) suggests Oh∗≃0.043, in the range that would roughly point to Duchemin’s data. Therefore, a detailed analysis of the process of collapse is needed, in particular the effects of microbubble trapping and the onset of ejection at the time singularity represented by t0(if such a singularity occurs). 2.2.1. The nonlinear behaviour around Oh∗: focusing, bubble trapping and the start of ejection Previous works (Duchemin et al. (2002), Ghabache et al. (2014)and Deike et al. (2018), among others) have postulated that the origin of these ballistic ultrafast jets is a finite-time singularity event at Oh∗. In reality, a true singularity followed by an instantaneous ultrafast jet occurs in all instances where a capillary wave collapses nearly cylindrically at the bottom of the cavity, locally resembling the bubble collapse phenomenon described by Eggers & Fontelos (2015). However, given the geometrical asymmetry of that collapse in the axial direction (see figure 4a), its consequences are twofold. (i) A microbubble gets trapped. The dependence of its size on Oh is summarized in figure 4 foranamplerangeofOh numbers. (ii) More importantly, an extremely thin ligament of vanishing radial size and very large speed is ballistically emitted in the opposite axial direction from the tiny trapped bubble, as illustrated in figure 5. As discussed in § 2, a vortex ring (usually much larger than the tiny trapped bubble) with a global mechanical energy equivalent to that of the ejected liquid spout is emitted towards the liquid bulk (see figure 1)to keep the overall neutrality of the initial axial momentum (which can be termed a bazooka effect). 929 A12-9 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 A.M. Gañán-Calvo and J.M. López-Herrera 0.001 0.001 ξ–1 R/R0 5×10–4 0.005 0.010 0.050 0.100 0.500 0.005 0.010 0.050 Oh Oh∗ = 0.034, α3 = 10–3.5 Duchemin 2002(WR) Sakai 1989 Seon 2017 Brasz sim. 2018 Brasz sim. (Bo =0) Brasz (mp 195–215 µm) Brasz (mp 20–40 µm) This work Berny 2020 Duchemin 2002(TR) Ghabache 2016 Spiel 1995 Blanchard 1989 Garner 1954 Hayami 1958 Tedesco 1979 Oh∗ = 0.03, α3 = 10–2.9 Oh∗ = 0.03, α3 = 10–2.4 Figure 8. Experimental and numerical measurements of the radius of the first ejected droplet from different literature sources (see Gañán-Calvo (2017) for additional information). The WR and TR in Duchemin’s data denote ‘wide rim’ and ‘thin rim’ initial conditions, respectively. Dashed lines are the curves 0.18((Oh/Oh∗− 1)2+α3Oh/Oh∗2).TheOh range covers seawater bubbles in the range from 8 μm to 2 mm. (Note that while the model curves use different Oh∗and α3, data are represented using fixed values Oh∗=0.033 and α3= 10−3.55.) we can define a variable ξ−1R/Roand represent it as a function of Oh, where ξis a function derived from (2.26) which tends to 1 for very small Oh and Bo,givenby ξ=Ψ((Oh∗−Oh)2+α3Oh)−1≡1−β3Bo((Oh∗−Oh)2+α3Oh)−1.(2.27) That plot is given in figure 8, where we represent the variable ξ−1R/Ro(calculated with the fixed values (Oh∗,α 3,β 3)of the best fitting) as a function of Oh. This first representation shows a strong variability of data in the vicinity of Oh∗. However, this variability does not seem entirely stochastic: one may observe that despite that the data have been made non-dimensional with fixed values of (Oh∗,α 3,β 3),the closeness of some data sets to curves corresponding to different values of these parameters suggests that those data sets could be independently fitted to our model, as the dashed curves reveal. Intriguingly, the data sets showing largest variability are the ones corresponding to numerical simulations. 2.4.2. Experimental data sets: universality in terms of Oh and Bo The above observations lead to the imperative need to analyse the origin of the deviations in the numerical data sets. In fact, a careful examination of the nature of the different data sources reveal an apparently paradoxical finding: if one separately fits all the experimental data, i.e. excluding all data from numerical simulations, one obtains a robust best fit around our model for Oh∗=0.0296, α3=0.00405 and β3=−0.000158. This yields a minimum normalized standard deviation of 9 % (significantly smaller than the 24 % obtained including numerical simulations) and a prefactor kr=0.186 Oh∗2.Figure 9 shows the resulting optimum fitting, and the distribution of data around the model compared with a normal distribution. Observe that the data of Brasz et al. (2018),and of Ghabache et al. (2016) and Séon & Liger-Belair (2017), provide very valuable resources for comparison around Oh∗despite the experimental errors and the particular difficulty of that region due to the droplet smallness. 929 A12-16 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 Transient ejection from bubble bursting 0.001 0.01 Oh∗2 Ψ–1 R/R0 0.05 0.10 0.50 1.00 0.005 0.010 0.050 Oh 0.001 0.01 Oh∗2 Ψ–1 R/R0 Oh∗2 Ψ–1 R/R0Oh∗2 Ψ–1 R/R0 0.05 0.10 0.50 1.00 0.005 0.010 0.050 Oh Ghabache 2016 Spiel 1995 Blanchard 1989 Garner 1954 Hayami 1958 Tedesco 1979 Sakai 1989 Seon 2017 Brasz (mp 195–215 µm) Brasz (mp 20–40 µm) Duchemin 2002 (WR) Sakai 1989 Seon 2017 Brasz 2018 (simul.) Brasz (Bo = 0, simul.) Brasz (mp 195–215 µm) Brasz (mp 20–40 µm) This work (simul.) Berny 2020 (simul.) Duchemin 2002 (TR) Ghabache 2016 Spiel 1995 Blanchard 1989 Garner 1954 Hayami 1958 Tedesco 1979 0 5 10 15 p.d.f. 20 25 30 0.05 0.10 0.15 0.20 0.25 0.30 0 5 10 15 20 25 30 0.05 0.10 0.15 0.20 0.25 0.30 (a) (b) (c)(d) Figure 9. The radius of the first emitted droplet, expressed as Oh∗2Ψ−1R/Ro:(a)experimental data;(b)the whole set of experimental and numerical data;(c)probability density function (black line) of data in panel (a) around the fitting, compared with a normal distribution with average 0.186 and standard deviation of 9% (blue dashed line), this distribution is nearly invariant with Oh;(d)probability density function of data in panel (b) around the fitting. The deviation from an approximate normal (average 0.185, standard deviation 11.1 %) shows a very large deviation on the tails of the distribution. In contrast, when the whole set of available data is compared with our model, the dispersion around Oh∗explodes as anticipated in figure 8:this is illustrated in figure 9(b). The statistics of both fittings are also given, noting that the distributions are nearly invariant along the Oh range except at the vicinity of Oh∗. This finding highlights once 929 A12-17 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 A.M. Gañán-Calvo and J.M. López-Herrera more the strong sensitivity of the focusing effect on the geometry and the initial conditions: (i) the artificial shape of the initial rim after bursting, which critically determines the initial surface energy content as the pioneering work of Duchemin et al. (2002) revealed; (ii) initial bubble shape; and (iii) the absence or presence of gravity along the evolution. In addition, the numerical accuracy as well as the models for free surface and surface tension reconstruction, etc., also affect the results. Note that the precision of our simulations does not avoid their maximum deviation: the fact that we used an initial edge shape thickness (initial minimum film thickness) equal to 0.1Roand initial hole radius 0.15Ro,leadstothe largest initial meridional radius of curvature among the published initial geometries (see figure 7), which shifts Oh∗from 0.03 to 0.034. In contrast, the actual initial surface energy content after film bursting in real experiments is expected to be fixed by a self-similar process (Eggers & Fontelos 2015, and references therein) with a universal result only in terms of Oh and Bo alone. In effect, despite the relative paucity of data around Oh∗and the fact that all the experiments available in this work have been performed with air at atmospheric pressure, naturally including interaction with the gaseous environment, gas compressibility, surface effects (presence of contaminants, surfactants, particles, etc.) and nonlinear breakdown and local detachment of the first droplet, the scatter compared with that of the numerical data is minimal. This supports the general assumption that this problem is fundamentally determined in nature by two non-dimensional parameters only: Oh and Bo. On the other hand, the considerably less numerous data available for the ejection velocity Vcould be similarly represented in the form ξ1/2V/Vo, which tends to V/Vo for very small Oh and Bo. However, since our model predicts the ejection velocity at the droplet formation point, while the available data correspond to the ejection velocity measured at a fixed point in space (usually, the initial level of the free surface), a direct comparison with our model predictions is not possible. To compensate for the effect of the different measurement point, which necessarily implies the additional action of viscous dissipation and the work of gravity, we propose a simple velocity correction of the form V/Vo=kvΨ−1/2(1+k1Bo +k2Oh)−1,(2.28) which, under the same approach as the droplet size approach, provides an optimal fit to the experimental data of (Ghabache et al. 2016)fork1=2.27 and k2=16, with a minimum normalized standard deviation below 10 % and a prefactor kv=0.39 (see figure 10a). Including the entire data set available for ejection velocity (numerical data from Berny et al. (2020) and our own), we observe the same drastic increase of the normalized standard deviation for droplet size. 3. Concluding statements Figure 9(a) provides the most robust basis for quantitatively and continuously representing the size of the ejected droplet for the whole range of Oh and Bo where droplet emission occurs. For the convenience of the reader, that universal size can be ultimately written, with typical errors below 10 %, as R=0.18RoOh Oh∗−12 +4.44Oh +0.19Bo;Oh∗=0.03.(3.1) The statistical distribution of experimental data around our model, with an approximately normal distribution and a standard deviation below 10 %, provides solid statistical grounds (i.e. the null hypothesis is virtually excluded) for its consistency and 929 A12-18 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 Transient ejection from bubble bursting 0.001 0.10 0.05 Ψ1/2 (1 + k1 Bo + k2 Oh) V/V0 0.50 1.00 5.00 0.005 0.010 0.050 0.001 0.10 0.05 Ψ1/2 (1 + k1 Bo + k2 Oh) V/V0 0.50 1.00 5.00 0.005 0.010 0.050 Oh Ghabache 2016 Duchemin 2002 (WR) Duchemin 2002 (TR) This work (simul.) Berny 2020 (simul.) Ghabache 2016 (a) (b) Figure 10. The velocity of the first emitted droplet, expressed as Ψ1/2(1+k1Bo +k2Oh)V/Vo,withk1= 2.27, with k2=16, for (a)the experimental data of Ghabache et al. (2016)andDucheminet al. (2002), after optimal data collapse according to (2.28)and (b)the whole set of available experimental and numerical data, using the same fitting as in panel (a). correctness: our model agrees with the available experimental data for the whole (Oh,Bo) range. We emphasize that (3.1) is not a fitting polynomial expression: it reflexes the physics of the phenomenon as described in this work using just two non-dimensional parameters. The main advantage of (3.1) compared with other proposed models which provide good results but are not valid beyond Oh∗(Gañán-Calvo 2017,2018; Gordillo & Rodriguez-Rodriguez 2019) is its uniform validity and continuity along the whole (Oh,Bo)range. Among an ample variety of applications, this provides a fundamental tool to describe the statistical distribution of the sea spray as a function of the observed bubble size distribution near the sea surface. A similarly accurate expression for the ejection velocity is not currently possible from available data, due to its strong dependency on the point where it is measured. Based on data at the level of the unperturbed free surface, we propose the following expression from an optimum fitting: V=0.39VoOh −Oh∗2+0.004Oh +0.00017Bo−1/2 (1+2.27Bo +16Oh)−1. (3.2) 929 A12-19 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 A.M. Gañán-Calvo and J.M. López-Herrera Of course, the complexity of the phenomenon can be augmented including surface viscosity (Ponce-Torres et al. 2017), Marangoni and non-Newtonian effects (Sanjay et al. 2021), the presence of immiscible liquids (Yang, Tian & Thoroddsen 2020)or the statistics of the emitted droplets following the first one. In this regard, Berny et al. (2020)presented a detailed analysis of all issued droplets from the scaling analysis of the first. Despite these are topics of subsequent studies, our results entirely capture the physics needed to describe this phenomenon. Acknowledgements. B. Gañán-Riesco is thanked for figure 7. Funding. This research has been supported by the Spanish Ministry of Economy, Industry and Competitiveness (grant numbers DPI2016-78887 and PID2019-108278RB) and by J. de Andalucía (grant number P18-FR-3623). Declaration of interests. The authors report no conflict of interest. Author ORCIDs. Alfonso M. Gañán-Calvo https://orcid.org/0000-0002-7552-6184; José M. López-Herrera https://orcid.org/0000-0002-3063-8676. Author contributions. A.M.G.C. derived the theory, analysed data, reached conclusions and wrote the paper. J.M.L.H. performed the simulations and helped in the writing. REFERENCES BERNY,A.,DEIKE, L., SÉON,T.&POPINET, S. 2020 Role of all jet drops in mass transfer from bursting bubbles. Phys. Rev. Fluids 5(3), 033605. BLANCHARD, D.C. 1989 The size and height to which jet drops are ejected from bursting bubbles in seawater. J. Geophys. Res. 94 (C8), 10999–11002. BLANCHARD,D.C.&WOODCOCK, A.H. 1957 Bubble formation and modification in the sea and its meteorological significance. Tellus 9(2), 145–158. BOUROUIBA, L. 2021 The fluid dynamics of disease transmission. Annu. Rev. Fluid Mech. 53, 473–508. BRASZ, C.F., BARTLETT, C.T., WALLS, P.L.L., FLYNN, E.G., YU, Y.E. & BIRD, J.C. 2018 Minimum size for the top jet drop from a bursting bubble. Phys. Rev. Fluids 3, 074001. DASOUQI, A.A., YEOM,G.-S.&MURPHY, D.W. 2021 Bursting bubbles and the formation of gas jets and vortex rings. Exp. Fluids 62 (1), 1. DEIKE, L., GHABACHE, E., LIGER-BELAIR,G.,DAS, A.K., ZALESKI,S.,POPINET,S.&SEON,T.2018 Dynamics of jets produced by bursting bubbles. Phys. Rev. Fluids 3, 013603. DUCHEMIN, L., POPINET,S.,JOSSERAND,C.&ZALESKI, S. 2002 Jet formation in bubbles bursting at a free surface. Phys. Fluids 14 (9), 3000–3008. EGGERS, J. 1993 Universal pinching of 3D axisymmetric free-surface flow. Phys. Rev. Lett. 71, 3458–3460. EGGERS,J.&FONTELOS,M.A.2015Singularities: Formation, Structure and Propagation.Cambridge University Press. EGGERS,J.,FONTELOS, M.A., LEPPINEN,D.&SNOEIJER, J.H. 2007 Theory of the collapsing axisymmetric cavity. Phys. Rev. Lett. 98, 094502. FONTELOS, M.A., SNOEIJER,J.H.&EGGERS, J. 2011 The spacial structure of bubble pinch-off. J. Appl. Maths 71, 1696–1716. GAÑÁN-CALVO, A.M. 2017 Revision of bubble bursting: universal scaling laws of top jet drop size and speed. Phys. Rev. Lett. 119, 204502. GAÑÁN-CALVO, A.M. 2018 Scaling laws of top jet drop size and speed from bubble bursting including gravity and inviscid limit. Phys. Rev. Fluids 3, 091601. GAÑÁN-CALVO, A.M., REBOLLO-MUÑOZ,N.&MONTANERO, J.M. 2013 Physical symmetries and scaling laws for the minimum or natural rate of flow and droplet size ejected by Taylor cone-jets. New J. Phys. 15, 033035. GARNER, F.H., ELLIS, S.R.M. & LACEY, J.A. 1954 The size distribution and entrainment of droplets. Trans. Inst. Chem. Engrs 32, 222–235. GHABACHE, E., ANTKOWIAK,A.,JOSSERAND,C.&SÉON, T. 2014 On the physics of fizziness: how bubble bursting controls droplets ejection. Phys. Fluids 26, 121701. 929 A12-20 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791 Transient ejection from bubble bursting GHABACHE, E., LIGER-BELAIR,G.,ANTKOWIAK,A.&SÉON, T. 2016 Evaporation of droplets in a champagne wine aerosol. Sci. Rep. 6, 25148. GHABACHE,E.&SÉON, T. 2016 Size of the top jet drop produced by bubble bursting. Phys. Rev. Fluids 1, 051901. GORDILLO,J.M.&RODRIGUEZ-RODRIGUEZ, J. 2019 Capillary waves control the ejection of bubble bursting jets. J. Fluid Mech. 867, 556–571. HAYAMI,S.&TOBA, Y. 1958 Drop production by bursting of air bubbles on the sea surface (1) experiments at still sea water surface. J. Oceanogr. Soc. Japan 14, 145–150. ISMAIL, A.S., GAÑÁN-CALVO, A.M., CASTREJÓN-PITA, J.R., HERRADA,M.A.&CASTREJÓN-PITA, A.A. 2018 Controlled cavity collapse: scaling laws of drop formation. Soft Matt. 14 (37), 7671–7679. KIENTZLER, C.F., ARONS, A.B., BLANCHARD,D.C.&WOODCOCK, A.H. 1954 Photographic investigation of the projection of droplets by bubbles bursting at a water surface. Tellus 6(1), 1–7. LAI,C.-Y.,EGGERS,J.&DEIKE, L. 2018 Bubble bursting: universal cavity and jet profiles. Phys. Rev. Lett. 121, 144501. LEE, J.S., WEON, B.M., PARK, S.J., JE, J.H., FEZZAA,K.&LEE, W.-K. 2011 Size limits the formation of liquid jets during bubble bursting. Nat. Commun. 2,367. MACINTYRE, F. 1972 Flow patterns in breaking bubbles. J. Geophys. Res. 77 (27), 5211–5228. MAXWORTHY, T. 1972 The structure and stability of vortex rings. J. Fluid Mech. 51 (1), 15–32. MONTANERO,J.M.&GAÑÁN-CALVO, A.M. 2020 Dripping, jetting and tip streamin. Rep. Prog. Phys. 83, 097001. PONCE-TORRES,A.,MONTANERO, J.M., HERRADA, M.A., VEGA, E.J. & VEGA, J.M. 2017 Influence of the surface viscosity on the breakup of a surfactant-laden drop. Phys. Rev. Lett. 118 (2), 024501. POPINET, S. 2015 Basilisk flow solver and PDE library. Available at: http://basilisk.fr/. SAKAI, M. 1989 Ion distribution at a nonequilibrium gas/liquid interface. J. Colloid Interface Sci. 127 (1), 156–166. SAMPATH,K.,AFSHAR-MOHAJER,N.,CHANDRALA, L.D., HEO,W.-S.,GILBERT,J.,AUSTIN,D., KOEHLER,K.&KATZ, J. 2019 Aerosolization of crude oil-dispersant slicks due to bubble bursting. J. Geophys. Res. 124 (10), 5555–5578. SANJAY,V.,LOHSE,D.&JALAAL, M. 2021 Bursting bubble in a viscoplastic medium. J. Fluid Mech. 922, A2. SÉON,T.&LIGER-BELAIR, G. 2017 Effervescence in champagne and sparkling wines: from bubble bursting to droplet evaporation. Eur. Phys. J. Spec. Top. 226, 117–156. SPIEL, D.E. 1995 On the births of jet drops from bubbles bursting on water surfaces. J. Geophys. Res. 100 (C3), 4995–5006. TEDESCO,R.&BLANCHARD, D.C. 1979 Dynamics of small bubble motion and bursting in freshwater. J. Rech. Atmos. 13, 215–226. THORODDSEN, S.T., TAKEHARA,K.,NGUYEN,H.D.&ETOH, T.G. 2018 Singular jets during the collapse of drop-impact craters. J. Fluid Mech. 848,R3. VERON, F. 2015 Ocean spray. Annu. Rev. Fluid Mech. 47, 507–538. WALLS, P.L.L., HENAUX,L.&BIRD, J.C. 2015 Jet drops from bursting bubbles: how gravity and viscosity couple to inhibit droplet production. Phys. Rev. E92, 021002(R). YANG, Z.Q., TIAN,Y.S.&THORODDSEN, S.T. 2020 Multitude of dimple shapes can produce singular jets during the collapse of immiscible drop-impact craters. J. Fluid Mech. 904,A19. YARIN, A.L. 2006 Drop impact dynamics: splashing, spreading, receding, bouncing...Annu. Rev. Fluid Mech. 38, 159–192. ZEFF, B.W., KLEBER,B.,FINEBERG,J.&LATHROP, D.P. 2000 Singularity dynamics in curvature collapse and jet eruption on a fluid surface. Nature 403, 401–404. 929 A12-21 Downloaded from https://www.cambridge.org/core. Universidad de Sevilla, on 24 Feb 2022 at 14:20:51, subject to the Cambridge Core terms of use, available at https://www.cambridge.org/core/terms. https://doi.org/10.1017/jfm.2021.791