Summarizing an Eulerian–Lagrangian model for subsea gas release and comparing release of CO2 with CH4

ContentslistsavailableatScienceDirect

Applied Mathematical Modelling

journalhomepage:www.elsevier.com/locate/apm

Summarizing an Eulerian–Lagrangian model for subsea gas release and comparing release of CO ₂ with CH ₄ ^∗

Jan Erik Olsen

^∗

, Paal Skjetne

SINTEF Industry, Postboks 4760 NO-7465 Trondheim, Torgarden, Norway

a rt i c l e i n f o

Article history:

Received 12 June 2019 Revised 13 October 2019 Accepted 23 October 2019 Available online 31 October 2019 Keywords:

CFD

Mathematical modeling Subsea gas release Safety

a b s t ra c t

Asubseagasreleaseisaconcernforbothsafety andenvironment.Thiscanbeassessed bymathematicalmodels.ThedevelopmentofanEulerian–Lagrangianmodellingconceptto studysubseagasreleasehastakenplaceovermanyyearsandthepiecewiseenhancements havebeendocumentedintheopenliterature.Themodelinitscurrentstateissummarized inthisarticle.Modelsimulationsareshowntobeconsistentwithdifferentexperiments varyingindepthfrom7to138m.Themodelcanbeappliedtoestimatehowgassurfaces intothe atmosphere fromasubsea source.Thisis vitalinput to risk assessments.Due to recentinterest insubseaCO₂ storageand transport,acomparison ofCO₂- and CH₄- releaseshasbeenperformed.Modelresultsshowthatamuchsmallerfractionofreleased CO2reachestheatmospherethanCH4duetothehighsolubilityofCO2inwater.

© 2019TheAuthors.PublishedbyElsevierInc.

ThisisanopenaccessarticleundertheCCBY-NC-NDlicense.

(http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction

Humanlife, assets andthe environment are atrisk during subsea gasreleases. The main causesof risks are fire and explosions due to surfacingof hazardous gases, capsizing of rigs andships dueto hydrodynamic loads andpollution of oceanwaters.Accidentalsubseagasreleaseisnormallycausedbypipelinefailureordrillingoperations.Numerousincidents happen annuallyandhistoricallymanyofthesehaveresulted inlossoflife [1].Theseincidentsprimarily involvenatural gas or methane. With the recent interest in carbon storage there is a growing risk of subsea release of CO₂ related to asphyxiation.Inordertopreventincidentsandapplyproperinterventionifanincidentoccurs,areliableriskassessmentis aprerequisiteforbothnaturalgasandCO₂releases.

Arisk assessment includesseveral stages.Oneof theseisestimatinghow the gasbubbles aredispersed inthe ocean anddistributed atthe surfaceifthey survive all the wayto thesurface. The predictedsurfaceflux can then be used as input to atmospheric dispersion calculations. Here we only consider the plume in the ocean. Historicallythis has been studiedbyso-calledintegralmodelsassumingcertainprofiles(e.g.Gaussian)forthevelocityandbubblevolumefractionin theplume[2-6].Such estimatescanalsobeperformedbytransientmultidimensionalCFDmodellingoflarge-scalebubble plumes.Withtheemergenceofaffordablecomputers,solving thefullsetoftheNavierStoke’sequationsbecamepossible.

CFD modellingof bubbleplumesdates back to the 1980s when Schwarz andTurner [7]developed an Eulerian-Eulerian

∗ This article belongs to the Special Issue: CFD2018 Melbourne.

∗ Corresponding author.

E-mail addresses: jan.e.olsen@sintef.no , dr.jan.erik.olsen@gmail.com (J.E. Olsen).

https://doi.org/10.1016/j.apm.2019.10.057

0307-904X/© 2019 The Authors. Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license.

( http://creativecommons.org/licenses/by-nc-nd/4.0/ )

modelandJohansen andBoysan [8]developed an axisymmetric Eulerian-Lagrangian modelforbubble plumesin bubbly reactors. Bubble plumes in the ocean differ from reactors by their much larger scale. This wasaddressed by Swan and Moros[9]whoappliedtheaxisymmetricEulerian-Lagrangianconcepttosubseagasreleases.Theywereabletoaccountfor thelargenumberofbubblesinsubseagasreleasesbytrackinggroupsofbubblesinsteadofallindividualbubbles.Buscaglia etal.[10]usedtheEulerian-Eulerianapproachinafull3Dsimulationfora77mdeepreleasewithamoderatereleaserate.

Formoreintense releaserates whichcan be observedin incidentalreleasefrom rupturedpipelines orwell blowouts, CFDmodellingwiththe Eulerian–Eulerianapproach failedto providereasonableresults [11].Forintense releasesthe gas isreleasedasajet.TheEulerian–Eulerian methodrequiresagridresolutionsmallerthanthereleasediametertobreakup theincomingjetintoadispersedflow.Combinedwiththeneedforaverylargecomputationaldomaincoveringtheentire watercolumn, thisisnotfeasibleduetothehighcomputationalcost.Thus,atransientfull 3DEulerian–Lagrangianmodel wasrecommendedby Cloeteetal.[11].Itis alreadyadispersed gasflow by definition(i.e.bubbly flow)asitenters the domain.Thismarked a shiftintechnologywhere full3D transientCFD basedon theNavier-Stokes equationsbecamean alternativeto the integral models historically appliedto thesestudies. The modellingconcept described inthe following sectionsistheresultoffurtherdevelopmentofthe workofCloete etal.[11]partlypublishedbefore[12-14].Thisarticle servesasasummaryofthisresearchwithafulldescriptionofthemostrecentversionofthemodellingconceptincludinga largersetofvalidationcases.AnewstudycomparingreleasebehaviorofCO₂andCH₄ispresented.TheEulerian–Lagrangian modellingconcepthasmorerecentlyalsobeenappliedtosubseagasreleasebyothers,includingworkbyFragaetal[15]. andLietal.[16].

2. CFDmodel

An Eulerian–Lagrangian modelling concept has been chosen to studythe large scale bubble plumes emanating from subseagasreleases.Thetransportofmass,momentumandenergyintheocean(liquid)andatmosphere(gas)issolvedin an Eulerianframe ofreferencebythe Navier–Stoke’sequations.The location ofthesea surfaceistrackedasan interface betweenliquid andgas. The bubble motion is calculated in a Lagrangian frame of referencewith Newton’ssecond law.

Bubbleandwatermotionarestronglycoupledthroughdragandturbulence.Theequationsdefiningthemodellingconcept isdescribedinthefollowingsections.

2.1. Flowdynamics

Theflowofbubbleplumesisdrivenbythe buoyancyprovided bythedensitydifferencebetweenthegasbubblesand thesurroundingwater.As thebubblesrises thedragforce transfersmomentum betweenbubblesandsurroundingwater.

Thisacceleratesthe waterinan upwardmotion.Thewatermoving upwardsisreplacedby waterwhichisentrained into theplumehorizontally.Thisentrainmentandturbulencespreadthebubbleplumeandgivesitashapeofacone.Sincethe bubblesaredrivingtheplumedynamics,thebubblemotionisdescribedfirst.

Groupsof bubbles known asparcels are tracked throughoutthe calculations. All bubbles ina parcel share the same properties(e.g. bubble size,density andviscosity). It is sumof all the bubbles in a parcelwhich defines thesource for interactionswiththecontinuousEulerianfield.Bytrackingthebubblesasparcelsinsteadofindividualbubbles,thecompu- tationalcostofmodellingbillionsofbubblesinaLagrangianreferenceframebecomesacceptable.Aforcebalancebasedon Newton’ssecondlawgivesthebubbleaccelerationas

du_b

dt =g

( ρ

b−

ρ )

ρ

b +FD+FVM (1)

ifwe account forbuoyancy(first termontheright-hand side),drag(secondterm)andvirtualmass (lastterm).Heret is time,u_bisbubblevelocity,gisthegravitationalacceleration,

ρ

bisthedensityofthebubblegasand

ρ

^{isthedensityofthe}

seawater(orcontinuousphase).Weusespecificforces,i.e.forcedividedbymass.Thespecificdragforceis FD= 18

μ

ρ

bd_b² CDRe

24

(

^u−ub

)

⁽²⁾

whereReistheReynoldsnumber,C_D isthedragcoefficient,andd_bisthebubblediameter.Thedragcoefficientisprovided bytheexpressionofTomiyamaetal.[17]forpartlycontaminatedconditions

C_D0=max

min

₂₄

Re

1+0.15Re^0.687

, 72 Re

, 8 3

Eo Eo+4

(3) Athighervolumefractionsofbubbles,thecorrectionofTsujietal.[18]isapplied:

CD=CD0

1− 3

α

b

π

²/3

(4)

where

α

bisthevolumefractionofbubbles.Thiscorrectionprimarilyaccountsforaccelerationoftrailingbubbles.Thedrag forceis proportionalto thevelocitydifference betweenthegasbubbles andthesurroundingliquidu_b−u.Note thatu is theinstantaneousvelocityofthesurroundingliquid

u=U+u (5)

whereUistheaveragevelocityanduistheturbulentfluctuationsofthevelocity.Thus,velocityfluctuationscauseturbu- lentdispersionofbubblesthroughthedragforce.Unlessturbulenceisfullyresolvedinthemodel,asub-modelisrequired to account for turbulent dispersion of the unresolved velocity fluctuations. For Lagrangian tracking of bubbles (or parti- cles/parcels)weapplyarandomwalkmodel[19]wheretheturbulentvelocityfluctuationsiscalculatedby

u=

ξ

k (6)

Here

ξ

^{isa Gaussianrandom numberwithzeromeanandunit variance andk} istheturbulent kineticenergyofthe velocity spectrumnot resolved bythemodel.The time forwhichthisvelocityfluctuation isapplied intheintegrationof thebubbletrajectoryislimitedbytheeddylifetime(orthetimeittakesforabubbletotraversethroughaturbulenteddy).

Theeddylifetimeis[14]

τ

^e=3₂^Cμ k

·MIN

1;

k^3/2

(7)

Virtualmassforce alsoknownasaddedmassforceistheforceperceivedbythebubblesincewhenitisacceleratingit deflectsandacceleratessomeofthesurroundingwater.Thespecificvirtualmassforceisgivenas

FVM=CVM

ρ ρ

b

Du Dt −du_b

dt

(8)

ForthevirtualmasscoefficientCVM=0.5isapplied.Lift forceisnormallyincludedinreactormodelling,butsensitivity studiesshow noeffectoftheliftforceintypical bubbleplumesinopen waters.Thisisduetotheabsenceofwallsclose tothebubbles.Insuchscenariostheshearrateisrelatively smallandtheliftforce canbediscarded[12].Bubbleinduced turbulencecan be importantcloseto thereleasepoint wherethe bubblyflow isdense.Since theflow is notsufficiently resolved in this region andthe mechanism is not properly understood(see discussion by Olsen etal. [14] andSchwarz [20]),bubbleinducedturbulenceisneglected.Thiscancauseanuncertaintyintheestimate offorcesandspreadinginthe plumeclosetotherelease.Fortunately,thisregionisnormallyverysmallformostscenarios.

Anestimateofthebubblesizeisrequiredforcalculatingthedragcoefficient(seeEq.(2))andthemasstransferrate(see below). Indense andturbulentbubbleplumes, the bubblesize isgovernedby turbulent breakup andcoalescence.Mass transferandgasexpansionduetopressuregradientswilldominateindiluteplumes.LauxandJohansen[21]developeda bubblesizemodelforanEulerianframeworkaccountingforbreakupandcoalescence.ThiscanberecastintoaLagrangian framework.Whenalsoincludingtheeffectofmasstransferandgasexpansion,thebubblesizeisexpressedbythefollowing differentialequation

d˙b=d_b^eq−d_b

τ

cb +db

3

m˙b

mb−

ρ

˙b

ρ

b

(9)

wherem_b isthemassofa bubble,m˙_bisthemasstransferratefromabubble,

ρ

˙b istheLagrangiantimederivative ofthe bubbledensity,

τ

cb isthetime scale forcoalescenceor breakup andd_b^eq isthe bubblediameterobtainedby a bubbleif itisexposed togivenflow conditions(turbulentdissipation andvolumeconcentration)foralongtime (i.e.equilibriumis reached).Theequilibriumbubblediameteris

d^eq_b =C1√

α

b

σ

_/

ρ

0.6

^0.4

μ

b

μ

0.25

+C2 (10)

where

α

bisthevolumefractionofbubbles,

σ

^{isthesurfacetensionand}

࢞isthesubgridturbulentenergydissipation.For themodelcoefficientsweassumeC₁=4.0(typicalforbubblesinliquids)andC₂=200

μ

^{m(smallestexpectedbubblesize).}

Forfurtherdetails,includingtimescaleforcoalescenceorbreakup,refertoLauxandJohansen[21].Theinitialbubblesize forbubblesreleasedfromthereleasesourceiscalculatedbyanempiricaljetmodel[22].Asthebubbleschangesize,their massalsochanges.Massconservationisthenassuredbyadjustingthenumberofbubblesintheparcelsuchthattotalmass ofaparcelisunchangedduringtheprocessofcoalescenceandbreakupofbubbles.

The motionof the water(i.e. backgroundfluid) iscoupled to the bubblemotion through dragforcesand turbulence.

FlowofseawaterandtheatmosphereaboveiscalculatedbytheNavier–StokesequationsinanEulerianframeofreference andutilizingaVOFmethod[23]whichalsotrackstheinterfacebetweenthecontinuouswaterandtheatmospherebythe geometric-reconstructscheme[24].Massconservationisensuredbythecontinuityequation:

∂ (

¹⁻

α

b

) ρ

∂

^{t +}

∇

^·

( ρ (

¹⁻

α

b

) v ₎

=S^c_b (11)

whereS_b^cisthesourcetermtothecontinuityequationaccountingformasstransferwiththeLagrangianbubblephase.We assume that theflow isdilute, andtheeffectofthe bubblescan be neglected intheabove equation. Thus weapply the followingcontinuityequation:

∂ρ ∂

^{t +}

∇

^·

( ρ v ₎

=S^c_b (12)

Sincemotioninbothwaterandatmosphereisaccountedfor,thecontinuousphase(i.e.backgroundfluid)isdescribed bymixtureproperties.Thus,thecontinuousphasedensityis

ρ

=

α

w

ρ

w+

(

¹−

α

w

) ρ

a (13)

wheresubindexwsymbolswaterandsubindexasymbolsatmosphere.Viscosityandotherpropertiesarecalculatedequally.

Forconservation ofmomentum we follow the sameprocedure as forconservation ofmass andassume dilute flow. The conservationequationthensimplifiestotheReynoldsaveragedNavier–Stokes(RANS)equationsforsinglephaseflow

ρ

^D_Dt^{U =}

ρ

^g−

∇

^p+

∇

^·

μ

eff

∇

^U+

∇

^UT

+S_b (14)

where

μ

effistheeffectiveviscosity(molecular+turbulent).ThesourceterminEq.(14),S_b,isthesourcetermduetodrag ofbubbles

S_b=¹⁸₂₄^μ_ρ^CDRe

bd2

b

(

^ub−u

) ζ

b

^t

^Vc

(15)

Here

ζ

b isthemassflow rateofbubblesmoving throughthecomputational cell,^{1 tis thetimestep andV}c isthe volumeofthecomputationalcell.ThistermprovidesastrongcouplingtotheLagrangianbubbletracking.Asimilarcoupling also appears for the virtual mass force. This is however almost insignificant andthus the expression is not listed here.

Theseequationsarebasedontheassumptionthat theflowisdilute(i.e.thevolumefractionofthebubblesissmall).This assumptionissometimesviolated intheregionclosetothereleasesource.Aslongastheregionwheretheassumptionis violatedissmallcomparedtothetotaldomain,theoverallresultsarenotaffected[25].Forscenarioswithagreaterextent ofnon-diluteconditions(i.e.shallowreleasewithhighreleaserate),afulltwo-waycoupling[26]couldbeimplemented.

Turbulenceandturbulentviscosity

μ

^{tisaccountedforbyaVLESmodel.Largescaleturbulenceisinherentlycapturedby}

themomentumequationsandonlysubgridscaleturbulenceismodelled.Practicallythisisdonebymodifyingtheturbulent viscosityinthek-

^model,

μ

^t,byintroducingafilterfunction[27]

μ

^t=Cμ

ρ

^k

² ·MIN

1;

k^3/2

(16)

Here^{isafiltersize(e.g.themeshsizecanbeusedforfiltersize)whichprovidesalengthscalebelowwhichturbu-} lenceisnotresolved.C_μ isamodelcoefficient.Theunresolvedturbulence(i.e.sub-gridturbulence)iscalculatedwiththe standardk-

^{model[28].Turbulentkineticenergyandturbulent}dissipationissolvedfrom

∂ ∂

^t

⁽ ρ

^k

)

+

∇

·

( ρ

^uk

)

=

∇

·

μ

+

μ

t

σ

k

∇

^k

+Gk+Gb−

ρ

(17)

∂ ∂

^t

⁽ ρ

)

⁺

∇

^·

( ρ

^u

)

⁼

∇

^·

μ

⁺

μ

^t

σ

∇

+C₁

k

(

^Gk+C₃G_b

)

^−C2

ρ

_k²

+S (18) whereG_kdenotesgenerationofturbulentkineticenergyduetomeanvelocitygradients

G_k=

μ

^t

∇

^u+

∇

^uT

−2 3

δρ

^k

·

∇

^{u (19)}

G_bdenotesgenerationofturbulenceduetobuoyancy G_b=

τ

^e·

μ

^t

∇

^uc+

∇

^ucT

·

∇

^g+

ρ

²₃

τ

^Lg·

∇α

⁽²⁰⁾

C_μ=0.09, C_1ε=1.44,C_2ε=1.92andC_3ε=1.0aremodelcoefficientscalibratedagainstexperimentaldata.Theadvantage oftheVLESmodelcomparedtothek-

^{modelisthatthelargescale turbulenceinherentlycalculatedby themomentum}

equationsdoesnotrelyonthemodelcoefficientswhichiscalibratedforsteadystateconditionsatmuchsmallerscales[29]. Turbulenceisdampedattheinterfacebetweenthecontinuousliquidphaseandthegasphaseabovebecauseturbulent structuresarenotcarriedthroughtheinterface.ThisisnotinherentlyaccountedforbyVOFmodelssincetheinterfacesare nottreatedasboundaries.Thusasourceterminthedissipationequationforturbulenceisaddedtoincreasedissipationand dampenturbulenceattheinterface[30].Itisdesignedsuchthat itforcestheenergydissipation,

ε

^{,closetothesurfaceto}

obtainavalueequivalenttoaturbulentlengthscaleapproachingzeroatthesurface.Thetargetvalueofenergydissipation atthesurfaceisgivenby

^new=C3μ/

κ

^{4kls3/2 (21)}

1The contribution of all parcels in the cell is added with adjustment for residence time in cell relative to time step. Volume fraction is calculated similarly.

l_s isthedistancetothesurfaceand

κ

=0.4isthevonKarmanconstant.Numericallytheenergydissipationismodified toapproachthisvaluebythefollowingsourcetermintheenergydissipationequation:

S_sd=N

(

^new−

)

⁽²²⁾

HereS_sdisthesourcetermduetosurfacedampingandNisalargenumber.Itneedstobesufficientlylargetoforce

ε

toapproach

ε

^{new,butnottoolargewhichwillcausenumericalstabilityissues.}

2.2. Masstransfer

Gas dissolution of gas speciesbetween bubbles andthe surrounding ocean becomes significant if the bubbles reside sufficientlylongintheocean.Gasdissolutionisamasstransferprocessandthemasstransferrateforasinglebubblesis

˙

mi=AbJi=

π

^d2bki·

c^sol_i −c^∞_i

(23)

whered_b is bubblediameter,k_iis masstransfer coefficient ofspeciesi,c^sol_i issolubilityof speciesi andc_i^∞ isthe back- groundconcentrationofspeciesiinthesurroundingwater.Thedrivingforceisthedifferenceinsolubilityandbackground concentration.Itappearsthat masstransferincreaseswithincreasing bubblesize.Thisistrueforsinglebubbles,butfora bubbleplumethegoverningparameteristhemasstransferratepermassofbubbles,i.e.thespecificmasstransferrate

m˙i

m_i= Ab

ρ

bV_iJ_i= 6

ρ

bd_bk_i·

c^sol_i −c^∞_i

(24)

Thetotalinterfacialareapermassorvolumeincreaseswithdecreasingbubblesize.Thus,smallerbubblesgivemoremass transferintotalforthewholeplume.ThebubblesizeisestimatedbyEq.(10).Animportantsourceofpotentialerroneous estimatesofgasdissolutionisthemasstransfercoefficientk_i.Themasstransfercoefficientvariessignificantlybetweenclean andcontaminatedconditions.Thisdependson thetypeandconcentrationofsurfactantsintheocean.Numerous typesof surfactantsexistintheocean.Theyaremainlytheproductsofbiologicalprocessesinvolvingphytoplankton[31]andinclude substancessuchaspolysaccharides,lipidsandmore[32].Thus,theamountofsurfactantsintheoceanvarieswithlocation, depth andseason [33] asdoesthe concentrationofphytoplankton.Largerbubbles tendtorefreshits boundarylayer due to vortexshedding. This increasesthedriving force ofthe masstransfer andhencelarger bubbles aremore likelyto act asbubblesincleanconditions.Althoughtherearelimitedexperimentsonbubblesinseawater,itisbelievedthatthemass transfercoefficientseesashifttocleanconditionswhenincreasinginsizebetween3and4.5mm[34].

Manycorrelationsforthemasstransfercoefficientexistsforbothcleanandcontaminatedconditions[35].Higbie[36]de- rivedanexpressionforthemasstransfercoefficientforcleanconditions

k_i= 2

√

π

ReSc_i D^w_i

d_b (25)

whereD^w_i isthediffusioncoefficient ofspeciesi inseawater.TheSchmidt number,Sc_i, describesthe ratiobetweenmo- mentumandmassdiffusionofspeciesiinthesurroundingwater

Sc_i=

μ ρ

^Dwi

(26)

Frössling[37]derivedanexpressionforrigidspheres whichrepresentsalowvalueformasstransferofbubblesincon- taminatedconditions

ki=

2+0.6Re^1/2Sc¹_i^/3

D^w_i

d_b (27)

Inordertoacknowledgethe observationthat thebubble’sshifttheirbehavior fromcontaminatedtocleanbehavior as thebubblesizeincreases,alinearshiftbetweenthecorrelationsofHigbieandFrösslingisappliedbetween3.5and4.5mm.

Thiscorrelationdescribesthemasstransfercoefficientforbubblesinseawaterandiscategorizedasacorrelationforpartly contaminatedconditions.ThecorrelationisplottedinFig.1.MoredetailsonthisisprovidedbyOlsenetal.[34].

The above explainsthat estimatesofgas dissolutionare sensitive to themechanisms governingthe masstransfer co- efficients.Eqs. (23)and(24)show thatgas dissolutionalsodependson solubility,c^sol_i ,andbackgroundconcentration, c^w_i. Background concentration is inherently calculated by the model via a species conservation equation forthe soluble gas components

∂

^cwi

∂

^{t +}

∇

^·

Uc^w_i

=D^w_i

∇

^2cwi +S^ci_b (28)

HereS^ci_b isthesourcetermaccountingfortheamountofspeciesiwhichisdissolvedintowater.

Fig. 1. Mass transfer coefficient as function of bubble size for clean [36] , contaminated [37] and partly contaminated conditions [34] .

Fig. 2. Gridding concept with 3 levels of grid refinement.

2.3.Gridandnumericalscheme

Theabove equationsaresolved inthecommercial softwareANSYS/Fluent16.2linked toa largeset ofdevelopedcode (user definedfunctions) formaterial properties,mass transfer, turbulencemodel (VLES) andmore.The non-linear set of differentialequationsfortheEulerianphasesissolvedwiththe PISOscheme(pressureimplicitwithsplittingofoperator) forthe pressure-velocity couplingand second orderdiscretization for the conservation equations. Interfacetracking and sharpeningare performedby the geometric-reconstruct scheme[38].Time stepis adjusted to keep theCourant number below0.25.

Thecomputationalmeshisbasedonacrudeuniformmeshwhichisrefinedbytheoct–threemethodaroundtheplume andtheoceansurface.Thisassuresanaffordablecomputationalcost,whichisstillhighcomparedtothetraditionalintegral models. The mesh is typically constructed from a basemesh with5 cellscovering the ocean height. Up to 4 levels of refinementisapplied totheregion aroundtheplume andtheoceansurface. Dueto theextentoftheocean surface,any furtherrefinementarounditisverycostly.Furtherrefinementisdonebelowtheoceansurfacetoassureacertainnumber ofgrid cellsacrossthe plume.The plume normallywidens ina coneshape andthus the finestresolutionis requiredin thelower regionwheretheplume isnarrow.Thisisalso theregionwiththehighestvelocityandthus theconstrainton thetime stepisgivenbythehighvelocityandsmallergridsizeintheregionclosetotherelease.Thegriddingconceptis showninFig.2.

Asensitivitystudyongridrefinementhasbeenperformed.Inordertoproperlyestimatethespreadingoftheplumeit isimportanttoresolvethegoverningphysicsasclosetothereleasepointaspossible.Byinitiatingtheturbulentstructures characteristicoftheVLES-modelasearlyaspossible,goodspreadingpredictionsareobtained.Withacoarsegrid,turbulence willnotbeproperlyestimated,andturbulent dispersionwillbeunderpredicted.Fig.3showsplumeillustrationsbasedon simulations with different grid refinement of a plume with a constant rateof 100kg/sfrom 300 m depth. 6 levels are

Fig. 3. Plume shapes with different levels of grid refinement. Low number is a coarse grid and high number a finer grid.

Fig. 4. Density and solubility of CH 4and CO 2as function of depth at 5 °C.

sometimesenoughforshallowerreleases.7levelsofrefinementissufficientforawiderrangeofdepths.Gridindependence upto1000mhasbeenverified.

2.4. Materialproperties

The propertiesof gas dependon temperature andpressure andthe propertiesof seawateritself is also a functionof salinity.ThedensityandviscosityofseawateristakenfromSharqawyetal.[39].Differentgasspecieshavedifferentprop- erties.Twogasesofconcernaremethane,CH₄,andcarbondioxide,CO₂.CO₂issignificantlyheavierandmoresolublethan CH₄.ThisisseeninFig.4wheredensityandsolubilityareplotted.ThepropertiesfoundinNISTREFPROP[40]areapplied.

Diffusivityisimportantformasstransfer.CO₂ hasroughly2timeshigherdiffusivitythanCH₄,whichenhancesmasstrans- ferforCO₂ comparedtoCH₄.Quantitativenumbersondiffusivityisfoundin[41]andtypicalvaluesare0.0011mm²/sfor CO₂and0.0008mm²/sforCH₄(takenat5°C).

3. Modelvalidation

Themodelhasbeenvalidatedagainstaseriesofexperimentsandobservations.

3.1. Rotvoll– releaseofairfrom7m

Engebretsen etal. [42] performed gas release experimentsat Equinor’s (formerly named Statoil) research facilities at Rotvoll inTrondheim. Theseexperimentswere used asthe mainvalidation setforthemodel initsearly stage whenthe k-epsilonturbulencemodelwasappliedandgasdissolutioncouldbeneglected[11].Aseriesofreleaseswereconductedin arectangularbasinwithadepthof7mandasurfaceareaof6×9m.Thebasinwasfilledwithwaterandairwasreleased atthebottomatgasratesof83,170and750Nl/s(equivalentto0.06,0.12and0.54kg/sreferredtothestateattheinlet).

Theinletwascomprisedofareleasevalvewitharapidlyactingpistoninjectinggasverticallywitharrangementsinfront ofittoreducetheverticalmomentum.Becauseofthismomentumbreaker,thefluctuationsinthegasflowandthelength oftheinletjetwereminimized.

Fig. 5. Comparison of rise time observed in experiments and predicted by native model (k-eps) and enhanced model (SURE III).

SimulationswiththelatestversionofthemodelincludingtheVLESturbulencemodeldescribedabovewasrunwithcase definitionsequivalenttothoseoftheexperimentatRotvoll.Wefocusedontherisetimeofthefirstbubblessurfacing.²This is recorded fromthe simulation results and the experimental observations. Both the old version of the model with the k-epsilonturbulencemodel)andthenewmodelareconsistentwiththeexperimentalresultsasseeninFig.5below.

Notethattherisetime,ortimeoffirstpenetratinggas,isdefinedasthefirstinstancewhenaLagrangiangasbubbleis detectedinacomputational cellwithavolumefractionofair(atmosphere)above40%.Thiscelldefinesthewatersurface.

The choice of 40% is to make sure that the bubbles don’t leave the liquid phase prematurely, but results are not very sensitivetothischoice.

3.2.Trolla– releaseofairfrom30m

AnexperimentpartoftheSUREprojectwithreleaseofairfrom30mdepthoutsideTrondheimatalocationknownas TrollawasconductedbySINTEF[43].Theexperimentwasconductedfromabargewithcompressorsdeliveringairthrough ahosetothespecifieddepth.Thearrangementcausedarampuptimeofthereleaseratewhichwasnotinsignificant.This needstobeaccountedforwhenperformingmodelsimulationsforcomparison.Intheexperimentstherisetimewasbased onthreeobservationtechniques;(1)anEchoscopeimagedtherisingplumebysonarsignals,(2)waveguidesmonitoredthe riseoftheoceansurfacebyconductivitymeasurementsand(3)acamerawasmountedatthetopofcranelookingdownat theoceansurface.Thisprovidedthreesomewhatdifferentdefinitionsofsurfacingoffirstgas.

The resultsare plotted inFig. 6 fordifferent releaserates andfor two differentrelease diameters (1 nozzle and2 nozzleswereapplied).Resultsfrommodelsimulationsandobservationsshow thattherisetimedecreaseswithincreasing gasrate. Thisisasexpectedsinceincreasinggasratescauseshigherinletvelocitiesandahigherbuoyantfluxpervolume.

The modelresults do not matchperfectlywith anyof theobservational techniques,butthe results are not dramatically off.Thesereleaseshaveaveryhighmomentumandtheresultingreleasejetpenetrates quitefarupinthewatercolumn.

Suchscenarios areamongthemostdifficulttopredictwiththemodellingconcept.Thisisparticularlydifficultforthe1 releases,since thereleasediameterisevenlessresolved.Inmostother casesthereleasejetonlyaffectsaverysmallpart ofthecomputationaldomain.

Itshould be notedthatthe 1 releasestake longerto surfacethanthe 2 releases.Dueto largerresistance inthe1 releasenozzle,ramp-uptimewastwicethatofthe2releases(7versus3.5sforthethreelowerrates).Thisexplainswhy the1releasesriseslowerthanthe2releases.

3.3.BuggsSprings– releaseofairfrom50m

Milgram[4] measured velocity profiles resulting from releaseof air in Bugg Springs (Florida) from a depth of 50 m fordifferent release ratesup to 0.71kg/s. Simulations withthe model were carried out withparameters andproperties equivalenttothoseoftheexperimentalvalues.ThecomparisonbetweenmodelandexperimentsisseeninFig.7.Hereradial profilesoftheverticalvelocityatdistancesof26and44mabovethereleaseisplotted.Notethattheexperimentalvelocity profileis basedona velocity measurement averagedover 10min andthencurve-fitted to aGaussian profile.The model resultsareaveragedover2minaftertheplumehasreachedaquasi-steadystateandthusthecurvesarenotassmoothas theexperimentalcurves.Thereisgoodconsistencybetweenexperimentalandcalculatedvelocityprofilesformostprofiles.

Forthevelocity profile44mabovethereleaseforthereleaserateof0.34kg/sthereis,however,anoverpredictionofthe

2The rise time is defined as the time between initial release of gas and surfacing of first gas

Fig. 6. Rise time as function of release rate estimated from different observational techniques and model simulations. Results from release with 1 nozzle is seen to the left and 2 nozzle to the right.

Fig. 7. Comparison of velocity profiles 26 and 44 m above release between model and experimental observations for release rate of 0.34 kg/s (left) and 0.71 kg/s (right).

peakvelocityby12%comparedtotheexperiments.Reasonsforthisdeviationhasnotbeenclarified.Othervelocityprofiles notshownherearemoreconsistent.

AcomparisonbetweentheVLESandk-epsilonturbulencemodelonpredictionofplume angleisshowninFig.8.The VLESmodelismoreconsistentwiththeexperimentalobservationsthanthek-epsilonmodel.Thek-epsilonmodelpredicts aplumeanglealmostindependentofgasrate,whereastheVLESmodelandtheobservationsshowthattheangleincreases withgasrateuptoacertainlevel.

3.4. Flags– releaseofnaturalgasfrom138m

DuringapiggingoperationintheNorthSeareleasesofnaturalgaswasmonitoredbyoceansonar,aerialcameralooking downattheoceansurfaceandotherinstruments[44].Threereleasesfrom138mwere studiedandonefrom380m.The mostintensivereleasesweretwoidenticalreleasesfrom138mwithareleaserateofapproximately17kg/swhichlastedfor 120s.TheascentoftheplumefrontasfunctionoftimeisseeninFig.9wheremodelpredictionandobservationofplume riseisshown.Thereissomesignalnoiseintheobservationsoftheearlypartoftherelease.Towardstheendoftheascent thestochasticbehavioroftheturbulenteddydominatingtheplumefrontissignificantandcancauseadeviationfromthe observation.The modelpredicts arisetimeof84s,while88swasobserved.Thedeviationisacceptable,especiallywhen consideringthestochasticinputfromtheturbulence.Themodelpredictionisconsistentwiththeobservations,butmultiple simulationswithastatisticaltreatmentcanimprovetheconsistencyfurther.

Fig. 8. Comparing plume spreading between Bugg Springs experiments and model predictions with VLES and k-epsilon turbulence models.

Fig. 9. Depth of plume front as function of time for 17 kg/s of natural gas from 138 m.

4. CasestudyonreleaseofCO₂vsCH₄

Releaseofnaturalgas andmethanehasbeenanalyzed fordecades duethesafetyaspects inoffshore operations. Due totherecentinterest incarbonstorageandsubsea transportlinesforCO₂,safetyassessmentsofCO₂ subseareleaseneed moreattention.Asanattempttobringsomeknowledgetothistopic,simulationsonCO₂releaseandCH₄releasehavebeen performed.Themodelpresentedabove wasapplied inthesesimulations.It shouldbenotedthat themodelhasnotbeen validatedagainstplumedataonCO₂releases.Ithasbeenvalidatedagainstplumeobservationsonreleasesofairandnatural gas.Itisthenassumedthatthemodelisvalidforothergasesifpropermaterialpropertiesareapplied.Acomparisontothe singlebubbleexperimentsonCO₂ [45]showedgoodconsistencybetweenmodelandobservations.Similarjustificationwas madebyDissanayakeetal.[46]whostudiedCO₂ plumeswithanintegralmodelaftervalidatingthemodelagainstairand naturalgasplumeobservationsandCO₂ comparisononlytosinglebubbleexperiments.

SinceCO₂ ismuchmoresolubleinwaterthanCH₄,scenarioswithsmallgasrateshasnotbeenchosensincetheseare believedtofullydissolveintheoceanfortheCO₂releases.Releasesfromadepthof100mwasstudiedfirst.Thevolumetric releaseratewaskeptconstantforbothgasesat426MMSCFDand1277MMSCFD.Thisisequivalentto100kg/sand300kg/s ofCH₄and275kg/sand825kg/sofCO₂.Thesearehighreleaserates.TheresultsareseeninFig.10.ForthereleaseofCH₄ of426MMSCFDabout25%ofthegasisdissolvedandfor1277MMSCFDabout16%isdissolved.Hence,themajorityofthe releasedCH₄reachestheatmospherewhilethereleasedCO₂ neverreachestheatmosphereduetosignificantlyhighergas dissolution.

Fig. 10. Images of bubbles plumes from release of CH 4and CO 2from 100 m. Colors indicates distance from the plume axis towards the viewer (equivalent to an image by a sonar).

Fig. 11. Images of bubbles plumes from release of 426 MMSCFD of CH 4and CO 2from 100 m. Colors indicates distance from the plume axis towards the viewer (equivalent to an image by a sonar).

Secondly areleaseof426MMSCFD from30mwasstudied.Theresults areillustrated inFig.11.Atthisdepththe gas reaches thesurface fasterand thereis lesstime forgas dissolution.Still only about1% ofthe released CO₂ reachesthe surfaceand about99% is dissolved. All ofthe CH₄ surfaces.Note that forCH4 this isan extremeblowoutat whichthe resultingvolumefractionsofgasbubblesexceedsthemodelassumptions(typicalforshallowreleaseatveryhighgasrate).

Thus,caremustbetakeninconcludingtoomuchfromtheCH₄resultforthisscenario.Stillitisclearthatgasdissolutionis practicallyinsignificantfortheCH₄release.Theseresultsindicatethatgasdissolutionisadominatingmechanisminsubsea gasreleaseofCO₂ andtheextentofthisreducestheriskforsurfaceoperations.Itshould benotedthatCO₂ dissolved in theoceanaffectsthepHlevelandthustheenvironmentalimpactmaybefarworsethanthesafetyimpact.

5. Conclusions

AnEulerian-LagrangianmathematicalmodellingconceptbasedonCFDforstudyingsubseagasreleaseshasbeenderived anddocumented.Simulationswiththemodelare consistentwith4differentexperimentsvaryingindepthfrom7to138

m.Releasesdominatedbyhighjetmomentumshowslessconsistencybetweenmodelandobservations.Themodelcanbe appliedtoestimatehowgassurfacesintotheatmospherefromasubseasource.Thisisvitalinputtoriskassessments.

AcomparisonbasedonsimulationresultsofCO₂-andCH₄-releaseswasperformed.Duethehighdensityandsolubility ofCO₂comparedtoCH₄,asmallerfractionofreleasedCO₂ willreachtheatmospherethanCH₄.Evenforveryhighrelease ratesandfromrelativeshallowreleases,thereisnotmuchCO₂reachingtheatmosphere.Thisisencouragingfromasafety aspect,buttheCO₂ dissolvedintheoceanaffectsthepHlevelandposesathreattomarinelife.

Acknowledgements

TheworkhasbeensupportedbyTotal,Equinor,WildWellControlandPSANorway.

References

[1] J.E. Olsen , P. Skjetne , Current understanding of subsea gas release - a review, Can. J. Chem. Eng. 94 (2016) 209–219 . [2] T.K. Fanneløp , K. Sjøen , Hydrodynamics of underwater blowouts, Norweg. Marit. Res. 4 (1980) 17–33 .

[3] O. Johansen , H. Rye , C. Cooper , DeepSpill - Field study of a simulated oil and gas blowout in deep water, Spill Sci. Technol. Bull. 8 (2003) 433–443 . [4] J.H. Milgram , Mean flow in round bubble plumes, J. Fluid Mech. 133 (1983) 345–376 .

[5] B.R. Morton , G.I. Taylor , J.S. Turner , Turbulent gravitational convection from maintained and instantaneous sources, Proc. Roy. Soc. A 234 (1956) 171–178 .

[6] L. Zheng , P.D. Yapa , Modeling gas dissolution in deepwater oil/gas spills, J. Marine Syst. 31 (2002) 299–309 .

[7] M.P. Schwarz , W.J. Turner , Applicability of the standard k- turbulence model to gas-stirred baths, Appl. Math. Model. 12 (1988) 273–279 . [8] S.T. Johansen , F. Boysan , Fluid dynamics in bubble stirred ladles: Part II. mathematical modelling, Metallur. Trans. B 19B (1988) 756–764 . [9] C. Swan , A. Moros , The hydrodynamics of a subsea blowout, Appl. Ocean Res. 15 (1993) 269–280 .

[10] G.C. Buscaglia , F.A. Bombardelli , M. García , Numerical modeling of large-scale bubble plumes accounting for mass transfer effects, Int. J. Multiph. Flow 28 (2002) 1763–1785 .

[11] S. Cloete , J.E. Olsen , P. Skjetne , CFD modeling of plume and free surface behavior resulting from a sub-sea gas release, Appl. Ocean Res. 31 (2009) 220–225 .

[12] J.E. Olsen , M. Popescu , On the effect of lift forces in bubble plumes, Progr. Comput. Fluid Dyn. 5 (2014) .

[13] J.E. Olsen , P. Skjetne , Modelling of underwater bubble plumes and gas dissolution with an Eulerian-Lagrangian CFD model, Appl. Ocean Res. 59 (2016) 193–200 .

[14] J.E. Olsen , P. Skjetne , S.T. Johansen , VLES turbulence model for an Eulerian–Lagrangian modelling concept for bubble plumes, Appl. Math. Model. 44 (2017) 61–71 .

[15] B. Fraga , T. Stoesser , C.C.K. Lai , S.A. Socolofsky , A LES-based Eulerian–Lagrangian approach to predict the dynamics of bubble plumes, Ocean Modell.

97 (2016) 27–36 .

[16] X. Li , G. Chen , F. Khan , Analysis of underwater gas release and dispersion behavior to assess subsea safety risk, J. Hazard. Mater. 367 (2019) 676–685 . [17] A. Tomiyama , I. Kataoka , I. Zun , T. Sakaguchi , Drag coeffcients of single bubbles under normal and micro gravity conditions, JSME Int. J. Ser. B 41

(1998) 472–479 .

[18] Y. Tsuji , Y. Morikawa , K. Terashima , Fluid-dynamic interaction between two spheres, Int. J. Multiphase Flow 8 (1982) 71–82 . [19] A.D. Gosman , E. Ioannides , Aspects of computer simulation of liquid-fuelled combustors, J. Energy 7 (1983) 4 82–4 90 .

[20] M.P. Schwarz , Bubble induced turbulence in two-fluid simulation if bubbly flow, in: Proceedings of the 26th International Symposium on Transport PhenomenaLeoben„ Austria, 2015 .

[21] H. Laux , S.T. Johansen , A CFD analysis of the air entrainment rate due to a plunging steel jet combining mathematical models for dispersed and separated multiphase flows, in: Fluid Flow Phenomena in Metals Processing, Symposium, San Diego, CA, 1999, pp. 21–30 .

[22] L. Zhao , M.C. Boufadel , S.A. Socolofsky , E. Adams , T. King , K. Lee , Evolution of droplets in subsea oil and gas blowouts: development and validation of the numerical model VDROP-J, Mar. Pollut. Bull. 83 (2014) 58–69 .

[23] C.W. Hirt , B.D. Nichols , Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225 .

[24] D.L. Youngs, Time-dependent multi-material flow with large fluid distortion, in: K.W. Morton, M.J. Banes (Eds.) Numerical Methods for Fluid Dynamics, Academic Press 1982.

[25] J.E. Olsen , S. Cloete , Eulerian-Lagrangian methods for modeling of gas stirred vessels with a dynamic free surface, in: Proceedings of the 7th Interna- tional Conference on Multiphase Flow - ICMFTampa, Florida, 2010 .

[26] D. Darmana , N.G. Deen , J.A.M. Kuipers , Parallelization of an Euler–Lagrange model using mixed domain decomposition and a mirror domain technique:

application to dispersed gas-liquid two-phase flow, J. Comput. Phys. 220 (2006) 216–248 .

[27] J.E. Olsen , P. Skjetne , S.T. Johansen , VLES Turbulence Model for an Eulerian-Lagrangian Modelling Concept for Bubble Plumes, CSIRO, 2015 . [28] B.E. Launder , D. Spalding , The numerical computation of turbulent flows, Comput. Methods Appl. Mech. Eng. 3 (1974) 269–289 . [29] S.T. Johansen , J. Wu , W. Shyy , Filter-based unsteady RANS computations, Heat Fluid Flow 25 (2004) 10–21 .

[30] Q.Q. Pan , J.E. Olsen , S.T. Johansen , M. Reed , L. Sætran , CFD study of surface flow and gas dispersion from a subsea gas release, in: Proceedings of the ASME 33rd International Conference on Ocean, Offshore and Arctic EngineeringSan Francisco, 2014 .

[31] B. Sabbaghzadeh , R.C. Upstill-Goddard , R. Beale , R. Pereira , P.D. Nightingale , The Atlantic Ocean surface microlayer from 50 °N to 50 °S is ubiquitously enriched in surfactants at wind speeds up to 13ms1, Geophys. Res. Lett. 44 (2017) 1–7 .

[32] I. Leifer , R.K. Patro , The bubble mechanism for methane transport from the shallow sea bed to the surface: a review and sensitivity study, Cont. Shelf Res. 22 (2002) 2409–2428 .

[33] J.D. Pakulski , R. Brenner , Abundance and distribution of carbohydrates in the ocean, Limnol. Oceanogr. 39 (1994) 930–940 .

[34] J.E. Olsen , D. Krause , E.J. Davies , P. Skjetne , Observations of rising methane bubbles in Trondheim Fjord and its implications to gas dissolution, J.

Geophys. Res. Oceans 124 (2019) 1399–1409 .

[35] R. Clift , J.R. Grace , M.E. Weber , Bubbles, Drops, and Particles, Academic Press, New York; London, 1978 .

[36] R. Higbie , The rate of absorption of a pure gas into a still liquid during short periods of exposure, Trans. A.I.Ch.E. 31 (1935) 365–389 . [37] N. Frössling , Über die Verdunstung falender Tropfen, Beitr. Geophys. 52 (1938) 170–216 .

[38] ANSYS, Fluent theory guide, www.ansys.com , 2015.

[39] M.H. Sharqawy , J.H. Lienhard , S.M. Zubair , The thermophysical properties of seawater: a review of existing correlations and data, Desalinat. Water Treat 16 (2010) 354–380 .

[40] E. Lemmon , M. Huber , M. McLinden , NIST Reference Fluid Thermodynamic and Transport Properties Database (REFPROP), 23, NIST Standard Reference Database, 2007 .

[41] W. Hayduk , H. Laudie , Prediction of diffusion coefficients for nonelectrolytes in dilute aqueous solutions, AiChE J. 20 (1974) 611–615 .

[42] T. Engebretsen , T. Northug , K. Sjøen , T.K. Fanneløp , Surface flow and gas dispersion from a subsea release of natural gas, in: Proceedings of the 7th International Offshore and Polar Engineering Conference, Honolulu, International Society of Offshore and Polar Engineers, 1997, pp. 566–573 . [43] P. Skjetne , J.E. Olsen , E.J. Davies , F. Leirvik , D. Krause , G. Eidnes , Comparison of Meso scale subsea gas release with multiphase Eulerian–Lagrangian

CFD model, ESRELSlovenia (2017) .

[44] E.J. Davies, P. Skjetne, D. Krause, B. Welde, D.R. de Miranda, Measurments of large-scale subsea gas plumes, J. Hydraul. Eng., To be submittedBull Trimest Plan Fam(2019).

[45] F. Takemura , A. Yabe , Rising speed and dissolution rate of a carbon dioxide bubble in slightly contaminated water, J. Fluid Mech. 378 (1999) 319–334 . [46] A . Dissanayake , J.A . DeGraff, P.D. Yapa , K. Nakata , Y. Ishihara , I. Yabe , Modeling the impact of CO 2releases in Kagoshima Bay, Japan, Jo. Hydro-Environ.

Res. 6 (2011) 195–208 .