Surface tension as the destabiliser of a vortical interface
UUnder consideration for publication in J. Fluid Mech. Surface tension as the destabiliser of avortical interface
Rashmi Ramadugu Prasad Perlekar † , and Rama Govindarajan TIFR Center for Interdisciplinary Sciences, Tata Institute of Fundamental Research, 500046,Gopanpally, Hyderabad, India International Centre for Theoretical Sciences, Tata Institute of Fundamental Research,Shivakote, Bengaluru 560089, India(Received xx; revised xx; accepted xx)
We study the dynamics of an initially flat interface between two immisciblefluids, with a vortex situated on it. We show how surface tension causes vorticitygeneration at a general curved interface. This creates a velocity jump across theinterface which increases quadratically in time, and causes the Kelvin-Helmholtzinstability. Surface tension thus acts as a destabiliser by vorticity creation, win-ning over its own tendency to stabilize by smoothing out interfacial perturbationsto reduce surface energy. We further show that this instability is manifestedwithin the vortex core at times larger than ∼ ( kW e ) / for a Weber number W e and perturbation wavenumber k , destroying the flow structure. The vorticitypeels off into small-scale structures away from the interface. Using energy balancewe provide the growth with time in total interface length. A density differencebetween the fluids produces additional instabilities outside the vortex core dueto centrifugal effects. We demonstrate the importance of this mechanism in two-dimensional turbulence simulations with a prescribed initial interface.
1. Introduction
The interaction between a vortex and an interface may be considered a buildingblock in the turbulent flow of immiscible fluids. We study this building blockand show that it is prone to a Kelvin-Helmholtz (KH) instability created bysurface tension. We will distinguish our flow from KH instabilities at immiscibleinterfaces across which a velocity jump is externally imposed. This latter class ofproblems has been well studied, and we begin by discussing a few studies. Theeffect of surface tension on the primary KH roll-up process was studied using two-dimensional numerical simulations by Fakhari & Lee (2013). The main findingwas that surface tension has a stabilising effect on the flow. Hou et al. (1997)showed that KH roll ups form only if surface tension is low, and the range ofunstable scales diminish with increase in surface tension. Rangel & Sirignano(1988) showed for a KH instability that non-zero surface tension results in anincrease of the stable regime. Tauber et al. (2002) investigated KH instability in † Email address for correspondence: [email protected] a r X i v : . [ phy s i c s . f l u - dyn ] F e b R. Ramadugu, P. Perlekar and R. Govindarajan density matched fluids at large Reynolds numbers. They find that the non-linearroll-up at low surface tension is similar to that at zero surface tension and highsurface tension results in a nearly flat interface with no roll up. Consistently acrossthese studies, surface tension thus acts as a stabiliser, by effecting a reduction ininterface length and thus suppressing KH roll-up. In fact, in a vast variety of flowsituations, surface tension suppresses large wave number perturbations, therebydecreasing interface area.It has been long known that the same quality of bringing about a reductionin surface area can make surface tension a destabilising agent, but in othercontexts. Famously, (Strutt & Rayleigh 1878) surface tension can destabiliseliquid jets and break them up into droplets by the so-called Plateau-Rayleighinstability. Again, this happens because of the propensity of higher surface tensionto effect a reduction in surface area in a circular flow geometry. In a planargeometry, Biancofiore et al. (2017) studied two parallel interfaces separatingthree immiscible density-matched fluids with linear shear profiles in a Taylor-Caulfield configuration. This system, which is stable without surface tension,is shown to display an instability when there is a phase lock between counter-propagating capillary waves. Thus, wave interactions cause surface tension to actas destabiliser. By a similar mechanism, surface tension at the interface in planarjets and wakes at high enough levels of shear can produce global instabilities(Tammisola et al. 2012).Surface tension has occasionally been reported as giving rise to small-scalestructures. Zhang et al. (2001) again studied the effect of surface tension on theKH instability. At high surface tension, they showed the generation of small-scale vortices in the late stages of evolution, giving rise to a positive contributionof surface tension to flow enstrophy despite a negative contribution to kineticenergy. The recent study of (Tavares et al. 2020) shows evidence, in RayleighTaylor turbulence, of a greater preponderance of smaller scales in immiscibleflows with surface tension as compared to miscible flows.We propose here a new mechanism for the destabilising action of surfacetension. The dynamics due to a single vortex placed at an interface between twofluids is studied in the absence of gravity and viscosity. This geometry is similar toDixit & Govindarajan (2010), but that study was at zero surface tension. We showanalytically how vorticity is produced by surface tension at the interface, and howthis makes the flow unstable. Our direct numerical simulations (DNS) confirmour analytical predictions, and show the differences in evolution of vorticity andthe interface.
2. Problem description
Two immiscible fluids of constant densities ρ and ρ lie on either side of aninitially flat interface in a two-dimensional system. The fluids are incompressibleand inviscid, and the continuity and momentum equations they each satisfy are DρDt = 0 , ∇∇∇ · u = 0 , (2.1) ρ D u Dt = −∇∇∇ P + F σ , (2.2) urface tension effects on a vortex at the fluid interface. − − − − a ) − − − − b ) 020406080100 − − − − c ) Figure 1: (a) The initial density field with an equal volume of two fluids shown in black andwhite.(b) The initial vorticity field, (Lamb-Oseen vortex with core radius( r c ) = 0.1). (c) Thefluid field at a later time, when the two fluids have the same density and surface tension isvanishingly small. where P is the pressure, u = { u r , u θ } are the radial and azimuthal velocitiesrespectively, F σ = σκδ ( x − x s ) n (2.3)is the surface tension force density, σ is the surface tension, κ is the curvature, x = { r, θ } , the subscript s stands for a location on the interface, δ ( . ) is the Diracdelta function and n is the unit normal to the interface at x s . The interfacepasses through the origin. A Lamb-Oseen vortex of circulation Γ and core radius r c is placed with its centre at the origin at time t = 0, as shown in Fig. 1(a,b).When surface tension is zero, and the two fluids have identical densities, eachfluid particle moves strictly in a circular path, with an azimuthal velocity givenby U = Γ πr [1 − exp ( − q )] , (2.4)where for ease of algebra we have defined q ≡ ( r/r c ) . The total angle, θ s , sweptout up to time t by the interface at r s , is a linearly increasing function of timegiven by θ s = Γ t πr s [1 − exp {− q s } ] . (2.5)For r s (cid:29) r c , Eq. 2.4 reduces to U = Γ/ (2 πr s ) for a point vortex, and an initiallyflat interface will wind up into an ever-tightening spiral (Dixit & Govindarajan2010), given by r s θ s = Γ t . Thus, away from the vortex core, at every instanceof time, the interface describes a different Lituus spiral, which is one among theArchimedean class of spirals, as seen at a typical time in Fig. 1(c).The relevant non-dimensional numbers are the Weber number (
W e ) which isa ratio of the inertial force to the surface tension force and the Atwood number( At ) which is a measure of the density contrast between the two fluids: W e ≡ ρU c r c σ , At ≡ ∆ρρ + ρ , (2.6)where ∆ρ = ρ − ρ , ρ = ( ρ + ρ ) /
2, and U c = Γ/ (2 πr c )[1 − exp( − r c . It is sometimes convenient to choose theradial distance r from the vortex centre as our length scale, so the Weber number W e r thus defined is a local quantity. At a given radius, the inertial time-scale R. Ramadugu, P. Perlekar and R. Govindarajan πr /Γ , for the wind-up of the spiral, is shorter than the time scale T σ = (cid:112) ρr /σ ,at which surface tension effects will be visible in the basic flow, by a factor W e / r .For large W e r and in the absence of instabilities, we may neglect the effect ofsurface tension on the interface shape up to a non-dimensional time of W e / r and the base velocity given by the Lamb-Oseen vortex will dictate the interfaceshape.
3. Vorticity generation on the interface and the Kelvin-Helmholtzinstability
An important aspect of this dynamics is the creation of vorticity at the interfaceby the surface tension and by the density contrast (baroclinic torque). This maybe seen with Eq. 2.2 rewritten in the vorticity formulation as
DΩDt = D ( Ω σ + Ω b ) Dt = 1 ρ ∇∇∇ × F σ − ρ ∇∇∇ ρ × ∇∇∇ P, (3.1)where Ω = ∇∇∇ × u is the vorticity, here pointing in the out-of-plane direction. Ω has contributions from the surface tension, Ω σ , and buoyancy, Ω b , representedrespectively by the first and second terms on the right hand side. An examinationof Eqns. 2.3 and 2.4 makes it clear that no vorticity will be generated by aperfectly flat interface or a perfectly circular one. But when the shape of theinterface deviates from these geometries, vorticity may be generated by bothdensity gradients and surface tension. Our focus here is on surface tension as agenerator of vorticity, so we discuss the case of At = 0 below.Integrating the first term in Eq. 3.1 at a given radial location, we obtain atime-dependent vorticity at the interface as follows. Consider f ( r, θ ) = r − r s ( θ )which vanishes on the interface, with a normal n ≡ ∇∇∇ f / |∇∇∇ f | , and curvature, κ ≡ −∇∇∇ · n . We have ∇∇∇ × F σ = − r ∂ θ F σ,r = σr ∂ θ (cid:18) r s + 2 r s ( ∂ θ r s ) − r s ∂ θθ r s [ r s + ( ∂ θ r s ) ] (cid:19) , (3.2)which leads, upon considerable simplification, to r c U c Ω σ = (1 − e − ) e q W eβ (cid:26) q + β ) (cid:20) − χ (cid:21) − (cid:20) q + 11 q + 18 q + 10 q + 5 +(2 q − q − q − e q + 5 e q (cid:21) (cid:20) − χ (cid:21) − β log χ (cid:27) δ (cid:18) rr c θ − rr c θ s (cid:19) , ≡ ∆UU c δ (cid:18) rr c θ − rr c θ s (cid:19) , (3.3)where β = [ q + 1 − e q ], χ = 1 + [2 t n β/ ( qe q )] , t n = t/T c , and T c = 2 πr c /Γ isthe inertial time scale at r c . Due to the vorticity created at the interface, thetwo fluids on either side move with different velocities. We denote the velocitycomponents parallel to the interface on either side by U || and U || + ∆U , with thejump ∆U across the interface given by Eq. 3.3. It follows that the interface mustbe subject to the Kelvin-Helmholtz (KH) instability. For the complete problem,an analytical dispersion relation is not possible to write down, but approximateestimates of the instability growth rates may be written down in two limitingregimes. urface tension effects on a vortex at the fluid interface. r (cid:28) r c , so q (cid:28)
1. Taylor expanding in this limit,after some algebra, and retaining the first surviving term in the expansion, Eq.3.3 reduces to ∆UU c = − − e − ) t n W e . (3.4)The vorticity produced, and thus the jump in velocity, are quadratic in time andproportional to the surface tension. In this limit, the interface can be closelyapproximated by a straight line rotating at a constant rate. We may then writethe relevant dispersion relation (Chandrasekhar 1981), in the case where there isno density contrast ( At = 0), in non-dimensional form as ωkU c = − ∆U U c ± (cid:34) kr c W e − (cid:26) ∆UU c (cid:27) (cid:35) / , (3.5)where the real and imaginary parts of ω respectively give the frequency andgrowth rate of a perturbation of wavenumber k . The first term within the squarebracket stands for the standard stabilising action of surface tension, increasingwith wavenumber. The second term on the other hand indicates the destabilisingeffect of surface tension. As seen from Eq. 3.4, it is quadratic in the surfacetension, while the stabilising term is only linear in this quantity. Moreover, weknow from Eq. 3.4 that ∆U increases quadratically in time, so the destabilisingaction of surface tension must win over its stabilising action at some time forany Weber number. In other words, the interface within the core becomes KHunstable when t n > (8 kr c W e ) / − e − ) / . (3.6)The higher the surface tension, the faster the instability grows. But at any giventime, for any finite W e , there is a cutoff wavenumber beyond which flow is stable.On the other hand, well outside the vortex core, q (cid:29)
1, and to leading orderin 1 /q we have ∆UU c = − − e − ) t n q W e , (3.7)so the vorticity on the interface decreases as r − . For estimating the instability,we may, by following a procedure analogous to Dixit & Govindarajan (2010),approximate the interface as a circle at a radius r , and perturbing it at azimuthalwavenumber m , to get ωr mU c = q − / + (cid:18) − At (cid:19) (cid:18) ∆UU c (cid:19) ± q − / (cid:20)(cid:26) (1 − At ) q m − (cid:18) − At (cid:19)(cid:27)(cid:18) ∆UU c (cid:19) + − (1 − At ) q / m ∆UU c − Atm + m W e (cid:21) / . (3.8)On examining this expression in the light of Eq. 3.7, we see that at At = 0the largest destabilising term is O ( q − / ) smaller than the stabilizing term, soinstability is not expected. When At >
R. Ramadugu, P. Perlekar and R. Govindarajan .
00 0 .
25 0 .
50 0 .
75 1 . r s /r c − − − ∆ U ( a ) t n = 1 t n = 2 − . − . . . . − . − . . . . b ) − − − − . − . . . . − . − . . . . c ) − − − − Figure 2: Vorticity generation with
W e = 12. (a) Velocity jump measured across the interface,as a function of interface location within the vortex core. Symbols: DNS, line: Eq. 3.3. Theperturbation vorticity from the DNS is shown for time t n = 1 in (b) and for t n = 15 in (c).The instability is well-developed by t n = 15 . All lengths are scaled by r c .
4. Direct Numerical Simulations (DNS)
Simulation details
We conduct DNS using an open-source Volume of Fluid (VOF) code Basilisk tosolve Eqs. 2.1 and Eq. 2.2 (Popinet 2018). We place a Lamb-Oseen vortex at theinterface in the center of the domain as shown in Fig. 1 (a,b) and allow it to evolvewith time. Since our computational domain, of length L = 5 πr c , is much largerthan our vortex, the far-field boundary conditions do not affect the results. We usefree-slip conditions on all sides in our inviscid flow, and have checked that periodicboundary conditions give practically indistinguishable answers. Note that oursimulations include non-Boussinesq effects when At (cid:54) = 0. We conduct DNS in asquare box of length 2 π , discretize it with 2048 collocation points, and vary theWeber numbers ( W e ) : 10040 , . , . , and 12. We restrict our simulations totimes over which the sum of the interfacial and kinetic energies is constant.4.2. Vorticity generation and the resulting instabilities
Up to time t n ∼ W e / we expect the vorticity generated on the interface in thenumerical simulations to closely follow Eq. 3.3. This is shown to be the case inFig. 2(a), where t n = 1 and t n = 2 for W e = 12. To calculate the velocity jumpacross the interface at the point r s , we numerically integrate the vorticity alongthe normal. We have checked that ∆U is insensitive to further grid resolution.As predicted, ∆U ∼ t n . Interestingly, the velocity jump changes sign, going fromnegative to positive below r c . There is one more sign change at large r (notshown here), and the velocity jump in the distant spiral arms is negative again,though very small. Fig. 2(b) and (c) show the perturbation vorticity at two times,calculated by subtracting the initial vorticity from the instantaneous field. Itis seen that at the larger time, instability has set in, which is consistent withthe prediction of instability when t n > .
08 for kr c = 2 π by the Eq. 3.4. Theinstability becomes visible in the simulations to the naked eye at around t n = 5(see movies in the supplementary material). For the case where surface tension( W e = 10040) is much lower, the time above which instability can occur, fromEq. 3.4, is very similar to that for
W e = 12, but because the growth rate isminuscule as predicted by the Eq. 3.5 , an instability does not become visibleduring the time of our simulation.The perturbation vorticity for
W e = 10040 and 12 is plotted in Fig. 3. InFig. 3(a) where At = 0, we find, in accordance with our expectations from the urface tension effects on a vortex at the fluid interface. − . − . . . . − . − . . . . a ) − − − − − b ) − − − − − − c ) − − Figure 3: Perturbation vorticity for
W e = 10040 and 12 with and without density contrast.(a)
W e = 10040, At = 0 .
0. Perturbation vorticity of small magnitude, in the neighbourhoodof r ∼ r c , is generated. Vorticity generated within the core is small. (b) W e = 12, At = 0at t n = 25. (c) W e = 10040, At = 0 .
05 at t n = 35. Here for both At , there are positive andnegative vorticity. For At = 0, the negative vorticity is mainly from the interior of the vortex. discussion in sections 3 and 4.2, that a minuscule amount of vorticity is generatedon the interface but no instability happens. Further, at r ∼ r c , this vorticity ispositive. However, for W e = 12, with At = 0 .
0, strong negative vorticity isgenerated within the core, and significant positive vorticity occurs in the vicinityof r ∼ r c . As time progresses, vorticity on the interface within the core causesa KH instability, and into a breakdown into small patches of negative vorticitywithin the core. At even later times, as is evident from Fig. 3(b) positive andnegative perturbation vorticities are interspersed, and a chaotic state in the coreregion ensues. This state causes the destruction of the initial Lamb-Oseen vortex.We notice a small level of asymmetry in the vorticity contours shown in 3(b), andbelieve this to be a numerical artefact. In a viscous flow, such dynamics wouldspeed up the dissipation of kinetic energy, and also enstrophy in the system, tobring it to a static state. We now examine, in Fig. 3(c), what happens when thereis a density difference in the two fluids. It is clear that when the surface tensionis low and the density contrast is high, there is a generation of alternating spiralsof positive and negative vorticity in the region outside the core. The blue arms ofthe spiral are unstable to centrifugal Rayleigh Taylor (CRT) modes, while the redarms are stable (Dixit & Govindarajan 2010). At later time we see a nonlinearbreakdown of the flow, as in the figure. When surface tension and density contrastare both high (see movie3 in the supplementary material), just outside the vortexcore, the CRT mode is stabilised by surface tension but displayed at long distancesfrom the core, as is to be expected from our simplified theory. And the vortexgets disrupted due to the rapidly growing KH instability. We observe that theheavy fluid moves outward from the core of the vortex due to centrifugal effects.4.3. Evolution of the interface
The interfaces for the cases in Fig. 3 are shown in Fig. 4. At low surface tensionand zero density difference, the spiralling interface of Fig. 4(a) is indistinguishablefrom that at zero surface tension, seen earlier in Fig. 1(c). Fig. 4(b) in combinationwith Fig. 3(b) shows that at high surface tension, the interface shape is completelydifferent from the vorticity distribution. Although vorticity is generated only atthe interface, it rolls up into small-scale structures independent of the interfacedue to the KH instability, and with time is spread through the vortex core. Theinterface meanwhile adopts a relatively short and straight shape in the centralregion, in deference to the high surface tension. In contrast, due to the low surface
R. Ramadugu, P. Perlekar and R. Govindarajan − − − − a ) − − − − b ) − − − − c ) Figure 4: Regime occupied by the two fluids, one shown in black and the other in white. (a)A spiral interface is seen at
W e = 10040 with At = 0 at t n = 25. (b) W e = 12, At = 0, t n = 25. The fine structure is erased, and a straight interface is seen in the central portion.(c) W e = 10040 with At = 0 .
05 at t n = 35. The CRT instability in the nonlinear regime isvisible. tension in Fig. 4(c), the interface closely mimics vorticity contours of 3(c). Herethe CRT instability is on display due to the density contrast.We thus have a competition between the response to the vortex, which increasesthe length of the interface, and surface tension, which acts to reduce it. Anindependent estimate of the interface length can be obtained from energy balance.As there is no external forcing and viscous dissipation here, the total kineticenergy is given by ∂ t E = (cid:104) u · F σ (cid:105) . (4.1)where the angle brackets refer to an average per unit area taken over the wholedomain, and E = (cid:104) ρu / (cid:105) . We also have (cid:104) u · F σ (cid:105) = − A ∂ t (cid:82) σdl where dl is aninterfacial line element and A is the total area of the domain (Joseph 1976), giving ∂ t (cid:18) E + σSA (cid:19) = 0 , i . e ., − E − E σ = SA , (4.2)where S = (cid:82) dl , and E the initial kinetic energy. We obtain the length S of theinterface numerically as a function of time from all four simulations, by trackingjumps in c from 0 to 1 and obtain excellent agreement with Eq.4.2 (not shown).4.4. Two-dimensional turbulence
To evaluate whether the mechanism we propose, of destabilisation of a flow bysurface tension, is of relevance in a general two-dimensional turbulent flow oftwo immiscible fluids, we conduct DNS at At = 0 on a doubly periodic box oflength L = 2 π and discretize it with 2048 collocation points. The two fluidsare initially separated by a flower-shaped interface (volume fraction of minorityphase is 0 . . × – and find no qualitative difference duringthe time of the simulation, though in the viscous case we have slowly decayingturbulence. The turbulent flow Weber number is defined as W e T = ρu rms L/σ where u rms is the root mean square velocity at time t = 0, and the eddy turnover time is τ e = L/u rms .The flow consists of vortices of several scales, encompassing a range of Webernumbers based on each, with the largest being two orders of magnitude smaller urface tension effects on a vortex at the fluid interface. a ) 0 2 4 60123456 (b) 0 2 4 60123456 (c)0 2 4 60246 (d) − − − − − − Figure 5: Density (a-c) and vorticity (d-f) contours for viscous simulations with Re = 5 . × .(a,b): initial profiles, (c,d): profiles at t/τ e = 5 . W e T = 5 × . (e,f): profiles at t/τ e = 5 . W e T = 5 × . Small-scale vortical structures occur in greater abundance as the surfacetension increases. than W e T . Moreover, the interface is not particularly designed to pass throughthe centre of any vortex. So, we do not expect an exact agreement with theory, butdo expect vorticity generation due to surface tension, and instabilities resultingin small-scale vortices. The regions occupied by the two fluids at t/τ e = 5 . W e T = 5 × , the two fluids, which displayed many spiralling interfaces atshort times (see movies in the supplementary material), are mingled intricatelyby this time, with elongated fine structures of the inner fluid. While the coresof the large vortices are still preserved, the smaller vortices undergo mergersand annihilation due to numerical viscosity. This picture is already qualitativelydifferent from flow at zero surface tension, in that we see vorticity generationdue to surface tension along the interfaces, as predicted. But the instability, andthe vorticity being peeled off from the interfaces is not yet visible. The pictureis starkly different at high surface tension. There is an explosion of small-scalevorticity everywhere in the flow, and the initial vortices have been disruptedcompletely. Because of these small-scale vortices, there is an enhancement ofenergy at large wave numbers with increase in surface tension. This is consistentwith earlier numerical findings (Li & Jaberi 2009; Trontin et al. 2010), which donot give the underlying instability mechanism. We have checked at early timesthat instability develops as expected.
5. Conclusion
A new role for surface tension, as a destabiliser, in the vortical flow of immisciblefluids is shown here. Due to surface tension, vorticity is generated practically0
R. Ramadugu, P. Perlekar and R. Govindarajan everywhere on the interface, and this vorticity increases with time. It followsthat the two fluid layers on either side move with different velocities, makingit conducive for the KH instability to manifest itself beyond a critical timeproportional to
W e / . This mechanism acts alongside the CRT instability whenthere are density differences. Our inviscid simulations on a single vortex confirmour theoretical predictions and also reveal a peeling off of small-scale vorticityfrom the interface. This mechanism is shown to have a significant presence inviscous and inviscid simulations of two-dimensional turbulence at low Webernumber, with an increased proportion of energy in small-scales. Acknowledgements
RG acknowledge support of the Department of AtomicEnergy, Government of India, under project no. RTI4001. PP and RR acknowl-edge support of the Department of Atomic Energy, Government of India, underproject no. RTI4007.
REFERENCESBiancofiore, L, Heifetz, E, Hoepffner, J & Gallaire, F
Physical Review Fluids (10), 103901. Chandrasekhar, Subrahmanyan
Dixit, Harish N & Govindarajan, Rama
Journal of Fluid Mechanics ,415–439.
Fakhari, Abbas & Lee, Taehun
Physical Review E (2), 023304. Hou, Thomas Y, Lowengrub, John S & Shelley, Michael J
Physics of Fluids (7), 1933–1954. Joseph, Daniel D
Li, Zhaorui & Jaberi, Farhad A
Physics of Fluids (9), 095102. Popinet, St´ephane
Annu. Rev. Fluid Mech. ,1–28. Rangel, RH & Sirignano, WA
The Physics of fluids (7), 1845–1855. Strutt, John William & Rayleigh, Lord
Proc. LondonMath. Soc (4). Tammisola, Outi, Lundell, Fredrik & S¨oderberg, L Daniel
Journal of Fluid Mechanics , 632–658.
Tauber, Warren, Unverdi, Salih Ozen & Tryggvason, Gretar
Physics of Fluids (8), 2871–2885. Tavares, Hugo S, Biferale, Luca, Sbragaglia, Mauro & Mailybaev, Alexei A arXiv preprint arXiv:2009.00054 . Trontin, P, Vincent, S, Estivalezes, JL & Caltagirone, JP
International Journal ofMultiphase Flow (11-12), 891–907. Zhang, Raoyang, He, Xiaoyi, Doolen, Gary & Chen, Shiyi
Advances in water resources24