Process-based upscaling of reactive flow in geological formations

Delft University of Technology

Process-based upscaling of reactive flow in geological formations

Meulenbroek, Bernard; Farajzadeh, Rouhi; Bruining, Hans

DOI

10.1016/j.ijheatmasstransfer.2020.119969

Publication date

2020

Document Version

Final published version

Published in

International Journal of Heat and Mass Transfer

Citation (APA)

Meulenbroek, B., Farajzadeh, R., & Bruining, H. (2020). Process-based upscaling of reactive flow in

geological formations. International Journal of Heat and Mass Transfer, 157, 1-12. [119969].

https://doi.org/10.1016/j.ijheatmasstransfer.2020.119969

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

International

Journal

of

Heat

and

Mass

Transfer

journalhomepage:www.elsevier.com/locate/hmt

Process-based

upscaling

of

reactive

flow

in

geological

formations

Bernard

Meulenbroek

a,∗

_,

_Rouhi

_Farajzadeh

b

_,

_Hans

_Bruining

b

a Delft University of Technology, Delft Institute of Applied Mathematics, Mekelweg 4, 2628 CD, Delft, The Netherlands b Delft University of Technology, Department of GeoScience and Engineering, Stevinweg 1, 2628 CN, Delft, The Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 17 July 2019 Revised 30 April 2020 Accepted 16 May 2020 Available online 18 June 2020 Keywords:

Reactive flow modeling

Peclet/Damkohler number phasespace Reservoir conditions

Temperature/pH/salinity dependence

a

b

s

t

r

a

c

t

Recently,thereisanincreasedinterestinreactiveflowinporousmedia,ingroundwater,agriculturaland fuelrecoveryapplications.Reactiveflowmodelinginvolvesvastlydifferentreactionrates,i.e.,differing bymanyordersofmagnitude.Solving theensuingmodel equationscan becomputationallyintensive. Categorizingreactionsaccordingtotheirspeedsmakesitpossibletogreatlysimplifytherelevantmodel equations.Indeedsomereactionsproceedsoslowthattheycanbedisregarded.Otherreactionsoccurso fastthattheyarewelldescribedbythermodynamicequilibriuminthetimeandspatialregionof inter-est.Atintermediaterateskineticsneedstobetakenintoaccount.Inthispaper,wecategorizeselected reactionsasslow,fastorintermediate.Wemodel2Dradiallysymmetricreactiveflowwitha reaction-convection-diffusionequation.We showthat wecan subdividethePeDaIIphasespaceinthreeregions.

RegionI(slowreaction);reactioncanbeignored,regionII(intermediatereaction);initiallykineticsneed tobetakenintoaccount,regionIII(fastreaction);allreactiontakesplacesinaverynarrowregionaround theinjectionpoint.Weinvestigatetheseaspectsforafewspecificexamples. Wecomputethelocation inphasespaceofafewselectedmineralsdependingonsalinityandtemperature.Wenotethatthe con-ditions,e.g.,salinityand temperaturemaybeessentialforassigning thereactiontothe correctregion inphasespace.Themethodologydescribed canbeappliedtoanymineral precipitation/decomposition problemandconsequentlygreatlysimplifiesreactiveflowmodelinginporousmedia.

© 2020TheAuthors.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense.(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Theefficiencyofimprovedorenhancedoilrecovery (IOR/EOR) processes is often influenced by the composition of the flowing aqueousphaseandisthusaffectedbythemassexchangebetween fluid and solid phases [1]. For example, the rheology of water-soluble polymers or the magnitude of the interfacial-tension re-ductionbysurfactantsstronglydependsontheionicstrengthand hardness(concentrationofdivalentcations)oftheaqueousphase and strongly influence the displacement efficiency [2]. In recent years,theadditionaloilextractedbytuningthecompositionofthe injectedwaterhasledtomoredetailedinvestigationofthenature oftheinteractionsbetweentherockandthefluidsresidinginthe pore[3].

Additionally, dissolution and precipitation of minerals have a considerable effect on the permeability of the reservoir [4]. The permeabilityimpairmentcanleadtohighcirculationcostsandloss ofinjectivity[5].Permeabilityimpairmentoccursbothdueto

par-∗ _{Corresponding author.}

E-mail address: b.j.meulenbroek@tudelft.nl (B. Meulenbroek).

ticlereleaseandsubsequentretention[5] andprecipitationof min-erals[4].

When a fluid is injected into an oil reservoir, the chemical speciesinthe solutionwillinteractwiththespeciesinthe reser-voir fluids contained in the pores and/or the substances on the rockgrains. Asthesespeciesaresubjecttotransportinthe reser-voir,in the numerical modelling of the process it is essential to know the final composition of the solution at the end of each timestep [6–8]. This is determined by how far the reactions are completedwithineachgridcellduringeachtimestep[7,8].Forfast reactions,thefrontdevelopsoverashortdistance,whileforslow reactionsfronts developover alongerdistancefromtheinjection point.

Subsurfaceformationsarecomposedofseveralminerals,which adds another complexity in modelling of geochemical reactions

[9,10].The time requiredtoreach equilibriumcanvary from sec-onds to (a few thousand) years for different minerals [11,12]. Therefore, in the numerical simulations, during a single time step,the concentration ofsome species could be over-or under-estimated.Therearetwoimportantmethodsingeo-chemical mod-elling,viz.,the Local Equilibrium Approach(LEA) andthe Kinetic Approach (KA) [13]. The LEA, which is the most common ap-https://doi.org/10.1016/j.ijheatmasstransfer.2020.119969

proachinreactive transportmodelling,assumesthatthereactions takeplace instantaneously andneglects theduration ofthe reac-tion.This approachis preferredbecause ofits simplicityandthe fact that it is computationally less expensive. However, the lo-calequilibriumapproachmayover- orunderestimatethe concen-tration of the species in the solution, e.g., due to heterogeneity

[14–16].

Indeed, the reactions in nature do not always reach equilib-riumwithin thetime andspaceofinterest. TheKineticApproach describes the time-dependency of the reactions. Some reactions complete faster than others and therefore the assumption of lo-cal equilibriumdoesnot holdin some practical cases.For exam-ple, in IOR/EOR applications, the injected fluid (which might be inequilibrium itself) is usually out ofequilibrium with the rock (this isreferred to aspartial equilibrium)[1]. The reservoir rocks consist ofseveral components whose reactions withthe injected fluid havedifferent rates.Therefore, the kinetics ofthe reactions inthis casedetermines the composition of thefluid (and conse-quentlythe efficiency of the injected chemical).This poses chal-lengeswhentheresultsofthelab-scaleexperimentsaretranslated tolarge-scalefieldapplications,becausedifferenttime andlength scales are involved. For this purpose,when the experiments are modelled,specialattention shouldbe given tocorrectly calculate the equilibrium time (Teq) and the equilibrium distance (Req) of therelevantminerals. Thesequantitiesare used asguidelinesfor selectionof the temporaland spatialresolutionof the numerical simulations.

It is difficult (if not impossible) to develop a general upscal-ingtechniquetoaccountforthegeochemicalreactionsindifferent IOR/EORprocesses.Apartfromthegeologicaluncertaintiesandthe intrinsicheterogeneity of the formations, the informationon the mineralogicalcompositionoftherocksisscarce.Evenifthe miner-alogyoftheformationispartiallyknown,thesurfaceareaavailable forthereactionsishardtomeasure.Notably,thebulk mineralogi-calcompositiondeterminedfromX-RayDiffraction(XRD)datacan significantlydifferfromthesurfacecompositionobtainedbyX-Ray PhotoelectronSpectroscopic(XPS)[17].

Therefore, itis suggestedto followa simplerapproach,which isprocessspecific.Insuchanapproach,onlyrelevantminerals af-fectingthe efficiencyofthe processare listed.Themainminerals consideredare Na (K,Ca) -silicate (carbonate) minerals(see [9– 12] forweathering rates). The particular interest is thento eval-uatethe time anddistancethat is requiredfor theinjected fluid toreachtheequilibriumstate. Theequilibriumtime(Teq)andthe equilibriumdistance(Req)areindicationsofthetemporaland spa-tialresolutionsofthesimulations.

It isthereforeourobjectivetoincludemineralreactionsinthe modeling ofreactive flow in porous media. This problemcan be studied on two characteristic length scales R_char: the well radius andthereservoirsize.Thischoicedefinesalsoacharacteristictime

Tchar,seeSection2.3.Thisgivestheopportunitytostudythe equi-libriumlength(Req) andequilibriumtime (Teq) ofthesereactions withrespecttothecharacteristiclengthandtime.Wedistinguish threepossibilities(I)Req Rchar,thewholereservoirwillbeatthe injectionconcentration(II)Req  Rchar,thereservoirwillbeatthe initialconcentrationand(III)Req~ Rchar,theconcentrationis vary-ingthroughoutthewholereservoir.Notethat incases(I)and(II) thereservoircanbeconsideredatequilibriumandkineticsdonot havetobetakenintoaccount.

The paper is organized asfollows. Section 2 gives the physi-cal model and model equations. In Section 2.6 the model equa-tionsarereducedtoasecondorderODEwithboundaryconditions.

Section3 showsandcomparesanalytical resultsfordifferent val-uesof the physical parameters. Section 4 computes Req, and we determinewhetherweareincase(I,II,III)foragivensetof min-erals.WeendthepaperwithsomeconclusionsinSection5.

2. Physicalandmathematicalmodel

Inthissectionwewillpresentsinglephasetransportand reac-tion (dissolution andprecipitationofthe omni-present minerals), using a reaction convection diffusion model. In the model equa-tionswewillusecapitallettersfordimensionalvariables,variables withsubscript c forcharacteristicscales andsmallletters for di-mensionlessvariables.

2.1. Reactionrates

Themethodologypresentedinthisworkcanbeappliedto gen-eraldissolution/precipitationreactions.Asillustrationweconsider thefollowingexplicitreactions:

K-feldspar:KAlSi3O8 +4H2O+4H+↔Al3+ +3H4SiO4 +K+ Kaolinite:Al2Si2O5(OH)4 +6H+↔H2O+2H4SiO4+2Al3+ Anhydrite:CaSO4↔Ca2++SO24−

Calcite:CaCO3↔Ca2++CO2₃−

ForthereactionratesS1weapplytheoftenadoptedsimplifying assumption (see, however, [18]) of using the law of mass action andobtain(see[10])

S1 =

κ

(

A1/V

)(

1 − IAP

K

)

n_{, (1)}

whereIAP is theionic activityproduct,K is theequilibrium con-stant, A1

V [m2/m3]isthereactivesurfaceareaper cubicmeter wa-ter;

κ

[mol/

(

m2 _s

₎

_{] is a temperaturedependent rateconstant. In} practice thevalue ofn is chosen to mimic experimentally deter-minedreactionrates,andherechosentobeequaltoone.Therate constantcanbeconsideredinitsArrheniusform

κ

= k10e−Ea/(RT), (2)

wherek10 isthefrequencyfactor,Ea istheactivationenergy,Ris thegasconstantandTistemperatureinK.

Asconcentrationofinterest,C,weusethedissolvedions,(e.g.,

Al3+) or dissolved molecules (e.g. H4SiO4). Moreover we assume that the reactions occur independent ofeach other and that we cansingleout onespecificionormolecule.Thisapproach cannot describethedecompositionofthemineral,butdoesallowto clas-sifyreactionsintoslow,intermediateandfastreactions.

WeassumefurthermoreisothermalflowatconstantpH,which means thatwe assume that theion activityproductis relatedto thedissolvedmineralconcentrationC,”independent” ofother con-centrations, which are presentin excess so that they can be as-sumedtobeconstant.VarioustemperaturesandpHvaluescanbe considered.Givenitsapproximatenatureitisnotwarrantedto as-sumenon-lineardependences,andweassumethat IAPis propor-tionaltoC,i.e.,

IAP ∼ C. (3)

Wenote that similarsimplifying assumptionsare madeinthe literature(see[7]).Wethenlumptheproportionalityconstantand

KtogetherintotheconstantC0andobtain S1=

κ

A1 V



1 − C C0



. (4) SettingA= A₁

ρwV wehavereactive surfaceinm2/kg-waterandS= S₁

ρw weobtainthereactionratein[mol/kg-water/s]andwefinally

obtain S = A

κ



1 − C C0



. (5) 2.2. Transportequations

Weassumeincompressibleflow,whichmeansthatwehavethe followingequationfortheinterstitialvelocityV

Transport isdescribed withareaction convectiondiffusion equa-tion

∂

C

∂

T +

∇

·

(

CV

)

= D



C + S, (7)

where the reactionrate S isgiven inSection 2.1. The concentra-tioninitially presentinthereservoirisdenotedby C0 [mol/ (kg-water)]; we assume that the reservoir initially is in equilibrium, meaningthatwedisregardslownaturallyoccurringvariationsthat take placeaftertimeslonger thanourtimeframeofinterest.This meansthat,fromapracticalpointofview,thereservoirisinitially atthe equilibriumconcentration (atthe initiallypresent environ-mentalconditions).Thismeansthatreactionwilltakeplaceifwe perturbtheinitialconcentrationC0toC.

We are looking for radially symmetric solutions V=

V

(

r,t

)

ˆr,C=C

(

r,t

)

, which means that equation (6) reduces to 1

R

∂

∂

R

(

RV

)

= 0 (8)

andequations(7)-(5) yield

∂

C

∂

T + 1 R

∂

∂

R

(

RCV

)

= D 1 R

∂

∂

R



R

∂

C

∂

R



+ A

κ

(

1 − C C0

)

, (9) whichreducesto

∂

C

∂

T + V

∂

C

∂

R = D 1 R

∂

∂

R



R

∂

C

∂

R



+ A

κ

(

1 − C C0

)

(10) because

∂

∂

R

(

RCV

)

= RV

∂

C

∂

R+ C

∂

RV

∂

R = RV

∂

C

∂

R (11) duetoequation(8).

Equations(8) and(10) aremadedimensionlessusingthe char-acteristic length Rc, characteristic time Tc (discussed below) and theinitialequilibriumconcentrationC0asfollows:

r = _RR c, T = _TT c ⇒

v

= V ·Tc Rc (12) c = C − C 0 C0 , (13) i.e., c=0 corresponds to the equilibrium concentration.

Equation(8) readsindimensionlessform 1 r

∂

∂

r

(

r

v

)

= 0 ⇒

∂

∂

r

(

r

v

)

= 0 (14)

andequation(10) reads

∂

c

∂

t +

v

∂

c

∂

r = 1 Pe 1 r

∂

∂

r



r

∂

c

∂

r



− Da Ic, (15)

whereweusedthePecletnumber Pe = R2c

DTc = VcRc

D , (16)

whichdenotestheratiobetweenconvectionanddiffusion.Wealso usedthefirstDamkohlernumber

DaI = A

κ

Tc C0 = A

κ

Rc C0Vc, (17) whichdenotestheratiobetweenreactionandconvection. Further-morewewillalsousetheDamkohlernumberII

DaII= PeDaI= A

κ

R2

c C0D

, (18)

whichdenotestheratiobetweenreactionanddiffusion.

2.3.Discussionofthescales

Twoscaleshavetobechosen:thecharacteristiclengthRcanda characteristicinterstitialvelocityVc(oracharacteristictimeTc).We are mainly interestedin clogging around the inlet, which means thatwewillsetRc=Rw(wellradius)inourcomputationsforthe fieldscale.ForreasonsofflexibilitywewillmaintainageneralRc, toallowotherchoices,likeRc=Rres(reservoirradius).

Solving equations(14) yields

v

=K₁

r (where K1 is a constant), whichisindimensionalform

V = VcRc K1

R;

thismeansthatVvariesthroughoutthereservoirandmultiple rea-sonablechoicesforVccanbemade.However,choosingthe charac-tervelocity asthe velocityatRw,i.e.,Vc=V

(

Rw

)

impliesthat the constantK1=1whichsimplifiesourexpressions;inthisworkwe willmakethischoice.

UsingthecharacteristicvelocityVcandthecharacteristiclength

Rcwehavethecharacteristictime Tc =

Rc Vc.

(19) Notethat the characteristictime dependson whetherwe choose thereservoirradiusorthewellradiusasourcharacteristiclength. InSection4 wewilldiscussthebehaviouroffourminerals (an-hydrite, calcite, kaolinite, K-feldspar) under two conditions (low salinityand alkaline flood) both atlow injection velocity/low Pe

(named ”field”) and higher injection velocity/higher Pe (named ”experiment”).Thefieldandexperimentalparametersaregivenin

Tables 2 and3; the parameters concerning thefour minerals are

giveninTable4.

2.4.Summaryofthemathematicalmodel

Thegoverningequationsare

∂

∂

r

(

r

v

)

= 0 (20) and

∂

c

∂

t +

v

∂

c

∂

r = 1 Pe 1 r

∂

∂

r



r

∂

c

∂

r



− Da Ic. (21)

We inject a solution at concentration C₌C_{in j =}C0 at a constant injectionvelocityatr=rw,i.e.,

V = Vc Rc Rw at r = rw, (22) i.e.,

v

= 1 rw at r = rw, (23)

ThismeansthatwehaveaDanckwertsboundaryconditionforthe concentrationatthewellr=rwR_Rw_c (notethatrvisconstantdueto equation(20))

(

r

v

)

c−_P1 _e·



r

∂

c

∂

r



=

(

r

v

)

cat r=rw, (24) wherec isthedimensionlessinjectionconcentration

c= Cin j− C 0

C0 .

(25) ForlargevaluesofPe,theboundarycondition(24) canbe approx-imatedby

Table 1

Nomenclature; capital letters are used to denote dimensional variables, small letters are used to denote the corresponding dimensionless variables.

variable description unit subscript

A reactive surface area m 2 /(kg-water)

0 initial

C concentration mol/(kg-water) inj injection

D diffusion coefficient m 2 /s

ss steady state DaI Damkohler number I (see Eq. 17 ) - tr transient DaII Damkohler number II (see Eq. 18 ) -

Pe Peclet number (see Eq. 16 ) - w well

R (radial) distance m res reservoir

S reaction rate ( Eq. 5 ) [mol / (kg-water s)]

T time [s] c characteristic

V velocity m/s

κ rate constant mol/(m 2 s)

Table 2

Parameters and scales (field).

well radius Rw 10 −1 [m] reservoir radius Rres 10 3 [m] characteristic length Rc 10 −1 [m] characteristic velocity Vc 10 −6 [m/s] diffusion coefficient D 10 −9 [m 2 /s]

Peclet number Pe 100 Table 3

Parameters and scales (experimental). characteristic length Rc 10 −1 [m] characteristic velocity Vc 5 · 10 −5 [m/s] diffusion coefficient D 10 −9 [m 2 /s]

Peclet number Pe 5000

whichmeansthatinthiscasetheconcentration isapproximately theinjectionconcentrationatr=rw.Atr=rreswehave

∂

c

∂

r = 0 , (27)

andatt=0wehavetheinitialcondition

c = 0 at t = 0 . (28)

ThedimensionlessparametersPeandDaIaregivenby Pe= RcVc

D , DaI = A

κ

Rc

C0Vc.

(29) NotethatPeisindependentofthelengthscalechosen(becausewe choseVc=V₍Rw

)

andV∼ _R1), but that both DaI and DaII=DaIPe dependonthelengthscale.

2.5.Finalnoteonthelengthscale

Besides that theDaIdependson thelength scale,thelocation ofthe boundary conditions depends on thischoice aswell. This changesthenatureofthesolutionsofourproblem.

IfweuseRc=Rw asourcharacteristiclengthscale,our bound-aryconditionsareat

r = R_Rw w =

1 and r= R_Rres w ≈ ∞ ,

(30) whereas ifwetake Rc=Rres asour characteristiclengthscale we haveourboundaryconditionsat

r = Rw Rres≈ 0

and r = Rres Rres =

1 . (31)

Thismeans that because thefunctional formof thesolution will change-theproblemsatthesetwolengthscalesshouldbetreated separately.InthisworkwewillfocusontheRc=Rw choice.

2.6. Analyticalsteady-statesolutionandtransients

Inthissectionwewillsolveourproblem(20)-(28) onthe char-acteristicscaleRc=Rw.

Thesolutionofthevelocityprofileisstraightforward

∂

∂

r

(

r

v

)

= 0 ,

v

(

1

)

= 1 ⇒

v

=

1

r. (32)

The solution of problem(20)-(28) is the superposition of a time independentsteadystateprofileandatransientprofile,i.e., c

(

r,t

)

= css

(

r

)

+ ctr

(

r,t

)

, (33) wherethesteadystatesatisfiesasecondorderODE

dcss dr = 1 Pe d dr



rdcss dr



− Da 1rcss, (34)

whichcanberewrittenas d2_c ss dr2 =

(

Pe− 1

)

r dcss dr + Da1Pecss. (35)

Theboundaryconditions(24)-(27) yield css − r Pe dcss dr = cat r = 1 and dcss dr = 0 as r→ ∞ . (36) Table 4

Mineral parameters low sal. Alaska at 60 C (field and experimental) and Alkali field, see [19] . The Damkohler numbers are calculated using equation (18) and Tables 2 and 3 .

mineral A [m 2 /kgw] _κ[mol/(m 2 s)] C

0 [mol/kgw] DaII Anhydrite (low sal.) 40 1 . 2 · 10 −3 _{1 . 74 · 10} −2 _{2.7 · 10} 7

(alkaline) 40 1 . 2 · 10 −3 _{4 . 7 · 10} −2 _{1.0 · 10} 7

Calcite (low sal.) 40 4 . 2 · 10 −6 _{2 . 63 · 10} −3 _{6 · 10} 5

(alkaline) 40 4 . 2 · 10 −6 _{9 . 4 · 10} −6 _{2 · 10} 8

Kaolinite (low sal.) 40 1 . 9 · 10 −17 _{3 . 8 · 10} −6 _{2 · 10} −3

(alkaline) 40 1 . 9 · 10 −17 _{1 . 2 · 10} −3 _{4 . 7 · 10} −6

K-Feldspar (low sal.) 40 3 . 11 · 10 −20 _{1 . 29 · 10} −4 _{1 . 1 · 10} −7

Fig. 1. Steady state concentration profile c ss as function of r (Regime A) for a fixed value of Pe = 10 5 and two values of Da II ; Da II = 10 −4 and Da II = 10 2 .

Fig. 2. Steady state concentration profile c ss as function of r (Regime A) for a fixed value of Pe = 10 5 and Da II = 10 8 . Note the different scale on the r -axis compared to Fig. 1 .

Forthetransientsweneedtosolvetheinitialvalueproblem r

∂

ctr

∂

t +

∂

ctr

∂

r = 1 Pe

∂

∂

r



r

∂

ctr

∂

r



− Da 1rctr (37)

withboundaryconditions ctr − r Pe dctr dr = 0 at r = 1 and dctr dr = 0 as r → ∞ . (38) andinitialcondition

ctr

(

r,0

)

= −c ss

(

r

)

. (39) The solution of the steady state problem is already sufficient to answerourresearchquestion;computationofthetransientsisnot necessaryinthispaper.

3. Analyticalsteady-statesolutionsofthemodelequations

Equations (35) -(36) determine css(r) forall values of Pe and

DaII.Findingthecompleteexplicitsolutionhoweverisdifficultdue thelargevariationinthevaluesofthePeandDaII.Wewill subdi-videthePeDaII− parameterspaceinthreeregimesandforeachof the threeregimes we willdetermine (anapproximation to) css(r) usingadifferentmethod.

RegimeADaII  Pe2,assumptionA:css(r)≈ css,A(r) RegimeB DaII ~ Pe2,assumptionB:css(r)≈ css,B(r) RegimeC DaII  Pe2,assumptionC:css(r)≈ css,C(r)

Notethat theassumptiondependsonthe regime,i.e., we will have a different functional form of the approximation for each regime.The solutionforRegimeAisgiveninSection 3.1,the so-lution for Regime B is given in Section 3.2 and the solution for RegimeCisgiveninSection3.3.Wethenmergetheseresultsand discussthefull PeDaII-phasespaceinSection 3.4,usingthe appro-priate approximations(A,B,C)forcss in each pointof the phases-pace.

3.1. Regimea:DaII Pe2

In Regime A the diffusion term d2css

dr2 can be neglected and

equations (35)-(36) can now be simplified and solved (see

AppendixA.1 fordetails)andthefinalresultisasfollows css = ce−α(r

2₋₁₎

,

α

= DaII

2

(

Pe− 1

)

. (40)

InFigs. 1 and2 we plotcss(r)for afixed Pe=105 andthree dif-ferentDaII numbers;DaII=10−4,DaII=102andDaII=108.We ob-servethreetypesofbehaviour;forsmallDaIIthewholereservoiris almostattheinjectionconcentration,seetheuppercurveinFig.1. forhighDaII thewhole reservoirremains atinitialconcentration, seeFig.2 andnotethedifferentscalesonther-axis.For intermedi-ateDaII weobservemixedbehaviour,seethelowercurveinFig.1. WewilldiscussthisinmoredetailinSection4.

Fig. 3. Steady state concentration profile c ss as function of r (Regime B) for a fixed value of Pe = 10 1 and Da II = 10 2 .

Fig. 4. Steady state concentration profile c ss as function of r (Regime B) for a fixed value of Pe = 10 4 and Da II = 10 8 . The figure seems similar to Fig. 3 . Note, however, the different scale on the r -axis compared to Fig. 3 .

3.2.Regimeb:DaII~ Pe2

InregimeBweusethefollowingrescalingoftheradial coordi-nate r= 1 + ¯r Pe⇒ d dr = Pe d d¯r , (41)

tosolveequations(35)-(36).Wefindthefollowingequationforcss (seeAppendixA.2 foradetailedderivation)

css = c 1 −1 2

(



1 + 4

β

− 1

)

e −1 2( √ 1+4β−1)(r−1)Pe_,

_β

₌ DaII Pe2. (42) Forsmallvalues ofDaII we observea non-horizontalprofile (see Fig.3);forlargervaluesofDaII theentirereservoirremainsatthe initialconcentration(seeFig.4).

3.3.RegimecDaII Pe2

InRegimeCtheconvectionterm (Pe_r−1)dcss

dr canbeneglected,as wewillshowinAppendixA.3.Solvingequations((35)-(36) yields (seeAppendixA.3 foradetailedderivation)

css = cPe



DaII

e−√ DaII(r−1) ₍₄₃₎

In this case we observe that the whole reservoir is at the ini-tialconcentration.Note thatthe widthofthe concentrationlayer

∼_√1

Da_II;thiswidth maybecomeinthe orderofthe porelength, which means that application of our Darcy scale model is ques-tionableinthispartofregimeC.

3.4. Phasespace

Inordertodiscussthebehaviourofcss(r)weintroducethe fol-lowingtwodistances:r0.99andr0.01 (seeFig.6).

The distance r0.99 is defined as follows: the reservoir is al-mostattheinjectionconcentration (99%of theinjection concen-trationorhigher)untilr=r0.99.Thismeansthatthewhole reser-voiris(eventuallyalmost)atthe injectionconcentration forlarge

r=r0.99.

The distancer0.01 isdefinedin a similar way:thereservoir is almostatthe(zero)initialconcentration (1%oftheinjection con-centration orlower) after r=r0.01. Note that the parameter r0.01 is also a measure of the width of the concentration profile, see

Fig. 6). This means that the reservoir remains atthe initial con-centrationifthiswidthr0.01− 1isverysmall,i.e.,r0.01≈ 1.

The distances r0.99 and r0.01 allow us to subdivide the full

PeDaII-parameter spacein four regions with distinct physical be-haviour.

For some values of (Pe, DaII) the reaction is so slow that the whole reservoiris (eventually)atthe injectionconcentration. We will call this subdomain of the (Pe, DaII) phase space ”RegionI”

Fig. 5. Steady state concentration profile c ss as function of r (Regime C) for a fixed value of Pe = 10 1 and Da II = 10 10 . Note the very fine scale on the r -axis, which runs from 1.0 0 0 0–1.0 0 01.

Fig. 6. The concentration profile c ss ( r ) for Pe = 10 3 and Da II = 10 2 (Regime B) The equilibrium distances r 0.99 = 1 . 4 and r 0.01 = 10 . 6 are shown as well. For 1 ≤ r ≤ r 0.99

we have c ss ≈ c ss (1); this part of the reservoir is at injection concentration. For r ≥ r 0.01 we have c ss ≈ 0; this part of the reservoir is at the initial concentration.

(”no reaction”). Region I is characterized by high values of r0.99; forourfigureswewilluser0.99=100asboundaryforRegionI.

Highreactionratesleadtoa situationwhereallreactiontakes place atthe inlet; almost no reaction takesplace inthe remain-derofthereservoir,whichremainsattheinitialconcentration.We will callthissubdomain ofthe (Pe,DaII) phase space”RegionIII” (”equilibrium”).RegionIII ischaracterized by valuesofr0.01 close toone;forourfigureswewilluser0.01=1.01asboundaryfor Re-gionIII.

The range of values of Pe and DaII in between Regions I and IIIinthephasespaceiscalledRegionII(”kinetics”).Inthispartof the phasespace css(r) is neither atthe initial nor atthe injection concentrationbutvariesbetweenthesevalues.InthisRegionII ki-neticshastobetakenintoaccount.

Veryhighreactionratesmayleadtoa situationwherethe re-action takesplacesover distancesoverthe orderofasinglepore

R_charr0.01∼ 10−6 m. Inthis casetheDarcy-approximation maybe questionable. We will call thissubdomain ofthe (Pe, DaII) phase spaceRegionIV(”Darcy-Brinkman”).

Insummary:

RegionI(”noreaction”)The reaction rate is so low that the reaction can be neglected. The whole reservoiriseventually atthe injectionconcentration.

RegionII(”kinetics”) Intermediate reaction rate; kinetics havetobetakenintoaccountinthe wholereservoir.

RegionIII(”equilibrium” High reaction rate; all reactions take place very close to the in-let, where the concentration varies from injectionto initial concentra-tion.Most ofthereservoir remains attheinitialconcentration. RegionIV(”Darcy-Brinkman”) Thewidthofthe reactionzone

be-comesoftheorderofthe poresize, which means that our Darcy-scale equationsmaynotholdanymore.

4. Quantitativebehaviourofmineralsinthreeregionsofthe phasespace

Bywayofexamplewechoosethreemineralsthatare represen-tativeforthethreeregionsI,II,adIII.Forallregionsweintroduce

teq,beingthetimeneededtoreachthesteadystate.ForregionsII andIII,wealsointroducetheequilibriumdistancereq ≈ r0.01.This equilibriumdistancereqisthedistancebeyondwhichthereservoir is(almost)attheinitialconcentration,seeFig.6.

Foran exactcomputationoftheequilibriumtimeteq the tran-sientproblem(37)-(39) hastobesolved.Wewillnotdosointhis paper;instead we estimate teq usingthe solution of thevelocity profile(32).

v

= 1 r ⇒ t = 1 2 r 2₋1 2 , (44)

which means that we can associate a dimensionless equilibrium timeteqtoanequilibriumdistanceasfollows

teq = 1 2

(

r

2

eq− 1

)

. (45)

Wewilluseequation(45) toestimatetheequilibriumtimesinthis section.

In regions I-III, we observe the following three types of be-haviour

RegionI noreaction; thereservoiris (dependingonrandt) ei-theratinjectionoratinitialconcentration.

RegionII req ≈ 1; kinetics have to be taken into account until

t ≈ teq. However, fort  teq the concentration profile maybeapproximatedbythesteadystateprofilecss(r). RegionIII req 1;thewholereservoirremainsattheinitial

con-centration.

4.1. DiscussionofK-Feldspar(regionsiandII)

ForK-FeldsparwehaveDaII=1.0· 10−7(lowsal)orDaII=2.0· 10−9 (alkaline pH > 7) (see Table 4). Using Fig. 7 (yellow dots)

Fig. 7. The full PeDa II -phase space is subdivided in four regions, viz., region I (re- action can be neglected), region II (kinetics has to be takes into account), region III (fast convergence to chemical equilibrium) and region IV (application beyond Darcy scale). The red dotted lines only indicate the different mathematical approximations used. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

we observethat K-Feldsparis (mostly) in RegionIin the PeDaII -phasespace.This is confirmedby the computation ofr0.99; using equation(40) wefindforPe=102_{, Da}

II=1.0· 10−7 ce−α(r2

0.99−1)= 0 .99 c⇒ r

0.99≈ 4 .5 · 10 3

(fortheother twocases r0.99 iseven bigger).Thiscorrespondsto adimensionaldistanceR0.99=450m,whichisoftheorderofthe sizeofthereservoir:eventuallythewholereservoirisatthe injec-tionconcentration.

Wecanapproximatethefulltimedependentconcentration pro-file(usingequation(44) by

css

(

r,t

)

=



0 t<1 2

(

r2− 1

)

c t ≥1 2

(

r 2_{− 1}

₎

. (46)

Thismeansthatat,e.g.,r₌10(correspondingto1m)the concen-trationis equalto zerountil t= 99

2 (corresponding toT= 992Tc≈ 12days)andequaltocafterwards.

4.2.Discussionofkaolinite(regionII)

For Kaolinite we haveDaII=2.0· 10−3 (lowsal)orDaII=5.0· 10−6(alkaline)(seeTable4).UsingFig.7 (purpledots)weobserve thatKaoliniteisinRegionIIinthePeDaII-phasespace.Thisis con-firmedbythecomputationofr0.99;usingequation(40)wefindfor, e.g.,Pe₌102_,Da II=2.0· 10−3, ce−α(r2 0.01−1)= 0 .01 c⇒ r₀ .01 ≈ 6 .8 · 10 2. and ce−α(r2 0.99−1)= 0 .99 c⇒ r₀ .99 ≈ 33 ,

which means that kinetics have to be taken into account into a partofthereservoir(between3.3and68m)untilthesteadystate

Fig. 8. The position of four different minerals in the PeDa II -phasespace under three different conditions. In some cases (K-Feldspar, Kaolinite) the minerals are in the same region for all circumstances investigated, whereas in other cases (Calcite, An- hydrite) the region has to be determined on a case-by-case basis.

is reached. The time required to reach the (spatially dependent) equilibriumsolutioncanbeestimatedusingequation(44) by teq ≈ 1

2

(

r

2

0.01− 1

)

≈ 2.3 · 105,

whichcorrespondstoTeq=2.5· 105 days.Thismeansthatkinetics hastobetakenintoaccount.

4.3. Discussionofcalcite(boundaryregionsIIandIII)

ForCalcitewe haveDaII=6.0· 105 (lowsal)orDaII=2.0· 108 (alkaline) (see Table 4). Using Fig. 7 (green dots) we observe that Calcite is either in RegionII or in Region III in the PeDaII -phasespace. Thisis confirmed by the computationof r0.01; using equation(42) wefindforPe=5· 103_,Da

II=6.0· 105that r0.01= 1 .04 .

The time required to reachthe (spatially dependent)equilibrium solutioncanbeestimatedusingequation(44)

teq ≈ 1 2

(

r

2

0.01− 1

)

≈ 0.04

whichcorrespondstoTeq≈ 0.05days,i.e.,roughly1h.Thismeans that for 0≤ T ≤ Teq kineticsneed to be taken intoaccount. For

T_{≥ T}eq thesteadystatesaturationprofile(42) maybeused. Forthealkalinecasewefind(Pe=102_andDa

II=2.0· 108) us-ingequation(43)

r0.01 = 1 .0 0 03 ,

which means that Calcite remains at the initial concentration in thiscase.

4.4. Discussionofanhydrite(regionIII)

ForAnhydrite we have DaII=2.7· 107 (low sal)or DaII=1.0· 107 _{(alkaline)(seeTable4).UsingFig.7(reddots)weobservethat} Anhydrite is in Region III in the PeDaII-phasespace. This is con-firmed by the computationof r0.01; using equation (43) we find forPe=102_,Da

II=1.0· 107that r0.01= 1 .001

(andsmallerfortheother cases),whichmeansthat Anhydrite re-mainsattheinitialconcentrationinallcases.

5. Conclusion

Reactiveflow modelinginvolves vastlydifferentreactionrates, i.e.,differingbymanyordersofmagnitude.Consequentlysome re-actions proceed so slow that in the time frame considered they canbedisregarded.Otherreactionsoccursofastthattheyarewell describedby thermodynamic equilibriuminthe time and spatial regionofinterest.Veryfastandveryslowreactionscanbe decou-pled fromthereactions thatoccur atintermediate rates;at inter-mediaterateskineticsneeds tobe takenintoaccount. Itis possi-ble to categorize selected reactions asslow, fast orintermediate. We model 2D radially symmetric reactive flow with a reaction-convection-diffusion equation. Using the solution of the steady state profilewe subdivide thePeDaII phasespacein threeregions. RegionI(slowreaction);reactioncanbeignoredandmineralsare convectedalong withtheflow. Theminerals areeitheratthe in-jection or initial concentration, depending on location and time. RegionIII(fastreaction);allreactiontakesplacesinaverynarrow regionaroundtheinjectionpoint.Thereservoirremainsatthe ini-tialconcentration beyondthisnarrow region.Kineticscan be ne-glectedin(almost)allofthecomputationaldomain.Finally,Region II (intermediate reaction);initially kineticsneed to be takeninto account.Eventually,beyondanequilibriumtimeaspatially depen-dent concentration profileisreached; afterthisequilibriumtime, whichcan be muchlongerthan thetime ofinterestthis concen-tration profilecanbeused; soingeneralkineticshastobe taken intoaccount.Weshow inexplicitexamples,thelocationinphase space of a few selected minerals depending on salinity, pH and temperature. The computation of the Pe and DaII numbers for a (any) specific processmakes it possibletodetermine inwhichof the Regimes I-IIIa process is,i.e., for whichtime window kinet-icshastobetakenintoaccount.Wenotethattheconditions,e.g., pH,salinityandtemperaturemaybeessentialforassigningthe re-actionto thecorrectregioninphase space.The methodology de-scribedcanbeappliedtoanymineralprecipitation/decomposition problemandconsequentlygreatlysimplifiesthemodelequations.

AppendixA. Mathematicalderivationsandproofs

A1. Concentrationprofileinregimea

In Regime A the diffusion term d2css

dr2 can be neglected, as

we will show in equations (A.4)-(A.7) below. This means that

equation(35) reducesto 0 =

(

Pe− 1

)

r dcss

dr + DaIIcss, (A.1)

whichcanbesolvedusingseparationofvariables;weobtain css = Ke−αr

2

,

α

= DaII

2

(

Pe− 1

)

, (A.2)

where the integration constant K can be determined using the boundarycondition(BC)(36).Wecomputecss(1)and dc_drss

|

r=1 css

(

1

)

= Ke−α

dcss

dr

|

r=1 = −2

α

Ke −α

andsubstitutetheseexpressionsintheBC(36) toobtain Ke−α

(

1 −−2

α

Pe

)

= c. Note that 2_Peα = Da_II

Pe(Pe−1) 1 which means that this contribution

can beneglected;solving forK yields K=ceα which meansthat weobtain

css = ce−α(r 2₋₁₎

,

α

= DaII

2

(

Pe− 1

)

. (A.3)

Wewillnowfirstshow thatthediffusiontermcanindeedbe ne-glected.Usingthesolution(A.3) wefind

dcss dr2 = K

(

α

2_r2₋

_α

₎

_e−αr2

, (A.4)

whichmeansthat thistermcan beneglected withrespectto the reactionterm(orequivalentlywithrespecttotheconvectionterm) if

DaIIcss = DaIIKe−αr 2

 K

(

α

2_r2₋

_α

₎

_e−αr2

. (A.5)

ThisholdsifbothDaII 

α

andDaII 

α

2,i.e.,if DaII 

α

=

DaII

2

(

Pe− 1

)

⇒ Pe 1 , (A.6)

whichisalwayssatisfiedinourregimeofinterestandif DaII 

α

2=

Da2 II

4

(

Pe− 1

)

2 ⇒ DaII  Pe

2_{, (A.7)}

whichissatisfiedduetoourrestrictiononRegimeA.

A2.Concentrationprofileinregimeb

InregimeBweusethefollowingrescalingoftheradial coordi-nate r = 1 + ¯r Pe⇒ d dr = Pe d d¯r , (A.8)

whichmeansthatequation(35) becomes Pe2d2css d¯r 2 = Pe

(

Pe− 1

)

1 + ¯r Pe dcss d¯r + DaIIcss, (A.9) whichreducesto d2_c ss d¯r 2 = dcss d¯r +

β

css,

β

= DaII Pe2, (A.10) because 1₊ ¯r

Pe≈ 1 and Pe

(

Pe− 1

)

≈ Pe2 because Pe is large. Eq.A.10 isasecondorderlinearodewithconstantcoefficients.The solutionofEq.A.10 is

css¯r =K1e− 1 2( √ 1+4β−1)¯r_+K 2e 1 2( √ 1+4β+1)¯r_{. (A.11)} Duetothesecondboundarycondition(36) wehaveK2=0which meansthatwefind css

(

¯r

)

= K1e− 1 2( √ 1+4β−1)¯r_{= K} 1e− 1 2( √ 1+4β−1)(r−1)Pe_{. (A.12)} Usingthefirstboundarycondition(36) wecandetermineK1 css

(

1

)

− 1 Pe dcss dr = K1

(

1 − 1 2

(



1 + 4

β

− 1

))

= c⇒ K1 = c 1 −1 2

(



1 + 4

β

− 1

)

(A.13) andwefind css = c 1 −1 2

(



1 + 4

β

− 1

)

e− 1 2( √ 1+4β−1)(r−1)Pe _(A.14)

A3.Concentrationprofileinregimec

In RegimeC theconvection term (Pe−1_r )dcss

dr can be neglected. Thismeansthatequation(35) reducesto

d2_c ss

dr2 = DaIIcss, (A.15)

whichcanbesolved css = K1e− √ DaIIr+ _K 2e √ DaIIr, _(A.16)

whereK2=0duetothe secondboundarycondition(36).The in-tegrationconstantK1 canbe determinedusingthefirstboundary condition(36) asfollows css

(

1

)

− 1 Pe dcss dr

|

r=1 = K1

(

1 +



DaII Pe

)

e −√ DaII= _c. Notethat DaII

Pe2  1inRegimeCwhichmeansthatthecontribution

of1+√DaII Pe ≈ √_Da II Pe ;weobtain K1 = ce √ DaII



Pe DaII , andwefinallyfind css = c _Pe



DaII e−√ DaII(r−1). _(A.17) We will now first show that the convection termindeed can be neglected.Usingsolution(A.17) wefind

dcss

dr = −P ece−

√

DaII(r−1), _(A.18)

whichcanbeneglectedwithrespecttothereactiontermDaIIcssif (orequivalentlywithrespecttothedifussionterm)if

DaIIcss = Pe



DaIIce− √ DaII(r−1) P_e

(

_Pe− 1

)

_c_e− √ DaII(r−1). _(A.19) Thisholds if



DaII Pe− 1, i.e., if DaII  Pe2, which is indeed satisfiedduetoourassumptiononregimeC.

AppendixB. AnalysisofthePeDa_II -phasespace

Using thenotionofthedistancesr0.99andr0.01,wecan subdi-videthe(Pe,DaII)-phasespaceinfourregions.

Thefirstdistancer=r0.99isassociatedwithRegionIandis de-finedasfollows

c ≥ 0 .99 cfor 1 ≤ r ≤ r 0.99, (B.1) i.e., the reservoir is at injection concentration until r0.99. This meansthatifr0.99≈ rres theentirereservoiris(eventually)almost atinjectionconcentration.WedefineourRegionInowasfollows: Region I =

{

(

Pe,DaII

)

|

r0.99 ≥ rres

}

, (B.2) wherethedimensionlessreservoir sizerres=103 withourchoice ofscales.

Theseconddistancer0.01 isdefinedsimilarly

c ≤ 0 .01 cfor r ≥ r 0.01, (B.3) i.e.,thereservoirisatinitialconcentrationbeyondr0.01.Thismeans thatifr0.01 ≈ 1mostofthereservoirisattheinitialconcentration. WedefineourRegionIIInowasfollows:

Region III =

{

(

Pe,DaII

)

|

1 + 10 −5≤ r 0.01 ≤ 1 + 10 −2

}

. (B.4) TheregioninbetweenRegionIandRegionIIIisRegionII,i.e.,

Region II =

{

(

Pe,DaII

)

|

r0.01>1 + 10 −2, r0.99<rres

}

. (B.5) Forvery smallvaluesofr0.01 wehaveto be careful.If r0.01<1+ 10−5,thedimensionlesswidthofthereactiveregionequals 10−5,

whichinfulldimensionsequals10−6mwithourlengthscale.This isoftheorderoftheporesize;thevalidityofourmodelequations becomesquestionable.ThisleadstothedefinitionofregionIV: Region IV =

{

(

Pe,DaII

)

|

r0.01<1 + 10 −5

}

. (B.6) Using the explicit analytical expressions for css derived in Section3 wewillsubdivide thefullPeDaII-phasespaceineach of thethreeregimesin(upto) fourregions inthefollowing subsec-tions.

B1. Regimea

InregimeAwehaveDaII Pe2;forournumericalpictureswe willuseDaII≤ 10−2Pe2.Wewillfirstcomputetheboundaryof Re-gionI.Usingtheexplicitformula(A.3)

css = ce−α(r 2₋₁₎

,

α

= DaII 2

(

Pe− 1

)

and combining the definition of r0.99 and the definition of the boundaryofRegionIasfollows

css

(

r = r0.99 = 10 3

)

= 0 .99 c,

wecanderivearelationbetweenDaII andPeontheboundary be-tweenregionsIandII,i.e.,

DaII =−2 ln

(

0 .99

)

10 6_{− 1}

(

Pe− 1

)

≈ 2· 10−8Pe.

ThismeansthatRegionIisbelowthelineDaII=2· 10−8Pe. Asimilarcomputationyields theboundarybetweenRegionsII andIII,wherewehave

css

(

r = r0.01 = 1 + 10 −2

)

= 0 .01 c, whichyields DaII = −2 ln

(

0 .01

)

(

1 + ·10 −2

₎

2_{− 1}

(

Pe− 1

)

≈ 0 .5 · 10 3_Pe.

WefinallycomputetheboundarybetweenregionsIIIandIV, css

(

r = r0.01 = 1 + 10 −5

)

= 0 .01 c, whichyields DaII = −2 ln

(

0 .01

)

(

1 + 10 −5

₎

2_{− 1}

(

Pe− 1

)

≈ 0 .5 · 10 6_Pe.

Combiningthesethree(blue)boundarylineswe obtainthephase spaceinFig.B.9.TheredlineDa=10−2Pe2 _{denotestheboundary} ofRegimeA;theequationsinthissubsectionareonlyvalidbelow theredline.

B2. Regimeb

InregimeB wehaveDaII ~ Pe2;forournumericalpictures we willuse10−2Pe2_<Da

II<102Pe2.Wewillfirstcomputethe bound-arybetweenRegionsIandII.Usingtheexplicitformula(A.12) : css = K1e− 1 2( √ Pe2_+4Da II−Pe)(r−1),

wecanderivearelationbetweenDaII andPe,i.e., −2 ln

(

0 .99

)

999 =



Pe2_{+ 4 Da} II − Pe, whichyields DaII

Pe <10−6whichholdsnowhereinRegimeB. Asimilarcomputationyields theboundarybetweenRegionsII andIII,



Pe2_{+ 4 Da} II− P e = −2 ln

(

0 .01

)

10 −2 = Nl,

wherethenumericalvalueNl ≈ 103.Thisyieldsthelower purple curve

4 DaII = Nl2+ 2 NlPe

inFig.B.10.TheboundarybetweenRegionsIIIandIVisfound sim-ilarly;thisboundaryisgivenbytheequation

4 DaII = Nu2+ 2 NuPe,

whereNu=106 andis representedby the upperpurple curve in Fig.B.10.

The red lines DaII=10−2Pe2 and DaII=102Pe2 denote the boundary ofRegime B; the equationsin thissubsection are only validbetweenthe redlines.Notethatthe”RegionI” inRegimeB isoutofourrangeofinterestof(Pe,DaII)values.

Fig. B.9. The PeDa II -phase space in Regime A. The dotted red line denotes the boundary of regime A, i.e., Da II ≤ 10 −2 Pe 2 . The lowest blue curve denotes the boundary between region I (no reaction) and II (kinetics), the middle blue curve denotes the boundary between region II and region III (equilibrium), and the up- per blue curve denotes the boundary between region III and region IV (Darcy- Brinkman). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. B.10. The PeDa II -phase space in Regime B. The dotted red lines denote the boundaries of regime B, i.e., 10 −2 Pe 2 ≤ Da

II ≤ 10 2 Pe 2 . The lower purple curve de- notes the boundary between region II (kinetics) and region III (equilibrium), and upper purple curve denotes the boundary between region III and region IV (Darcy- Brinkman). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. B.11. The PeDa II -phase space in Regime C. The dotted red line denotes the boundary of regime C, i.e., Da II ≥ 10 2 Pe 2 . The lower green curve denotes the bound- ary between region II (kinetics) and region III (equilibrium), and upper green curve denotes the boundary between region III and region IV (Darcy-Brinkman). (For in- terpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

B3. Phasespaceinregimec

InregimeCwehaveDaII Pe2;forournumericalpictureswe willuseDaII >102Pe2.Wewillusetheexplicitformula(A.17) css =

cPe



DaII

e−√ DaII(r−1)

tocomputetheboundariesbetweenRegionsIIandIIIandbetween RegionsIII andIV;the boundarybetweenRegions IandII isout ofthe range ofour interest of(Pe,DaII) values. We observethat the boundariesare ata constant DaII number. For the boundary betweenRegionsIIandIIIwehave

ln

(

0 .01

)

= −



DaII · 10 −2⇒ DaII ≈ 2 · 10 5

andfortheboundarybetweenRegionsIIIandIVwefindsimilarly ln

(

0 .01

)

= −



DaII · 10 −5⇒ DaII ≈ 2 · 10 11.

Theseboundariesare denotedbythe greenlinesinFig.B.11. Fur-thermoretheredlineDa=102_Pe2_{denotestheboundaryofRegime} C.Notethat the”RegionI” inRegimeBisoutoftherangeofour interestof(Pe,DaII)values.

Supplementarymaterial

Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ijheatmasstransfer.2020. 119969

References

[1] R. Farajzadeh , T. Matsuura , D. van Batenburg , H.e.a. Dijk ,Detailed modeling of the alkali/surfactant/polymer (asp) process by coupling a multipurpose reser- voir simulator to the chemistry package phreeqc, SPE Reservoir Evaluation & Engineering 15 (04) (2012) 423–435 .

[2] L.W. Lake , R. Johns , B. Rossen , G. Pope , Fundamentals of enhanced oil recovery, Soc. of Pet. Eng., Richardson, Tex, 2014 .

[3] H. Farajzadeh R.; Guo , J. van Winden , J. Bruining , Cation exchange in the pres- ence of oil in porous media, ACS Earth & Space Chemistry 1 (2) (2017) 101–112 . [4] H.A. Ohen , F.e.a. Civan , Formation damage in petroleum reservoirs i: modeling,

SPE-19380-MS, Society of Petroleum Engineers, 1989 .

[5] P. Bedrikovetsky , D. Marchesin , F. Shecaira , A. Souza , P. Milanez , Characteri- sation of deep bed filtration system from laboratory pressure drop measure- ments, Petroleum Science and Engineering 32 (2–4) (2001) 167–177 . [6] P. Lichtner , Scaling properties of time-space kinetics mass tranport equations

and the local equilibrium limit, Am. J. Sci. 293 (1993) 257–296 .

[7] R.B. Knapp , Spatial and temporal scales of local equilibrium in dynamic flu- id-rock systems, Geochim. Cosmochim. Acta 53 (8) (1989) 1955–1964 . [8] F. Golfier , C. Zarcone , B. Bazin , R. Lenormand , D. Lasseux , M. QUINTARD , On

the ability of a darcy-scale model to capture wormhole formation during the dissolution of a porous medium, J. Fluid Mech. 457 (2002) 213–254 . [9] J.I. Drever , Surface and ground water, weathering, and soils: Treatise on geo-

chemistry, 5, Elsevier, 2005 .

[10] C.A.J. Appelo , D. Postma , Geochemistry, groundwater and pollution, CRC press, 2004 .

[11] H. Sverdrup , P. Warfvinge , Calculating field weathering rates using a mechanis- tic geochemical model profile, Appl. Geochem. 8 (3) (1993) 273–283 .

[12] H. Sverdrup , P. Warfvinge , Weathering of primary silicate minerals in the nat- ural soil environment in relation to a chemical weathering model, Water Air Soil Pollut. 38 (3–4) (1988) 387–408 .

[13] S.A .-M.A . Lake L.W.; Bryant , Geochemistry and fluid flow, Gulf Professional Publishing, 2002 .

[14] Y. Al-Khulaifi, Q. Lin , M.J. Blunt , B. Bijeljic , Reservoir-condition pore-scale imaging of dolomite reaction with supercritical co2 acidified brine: effect of pore-structure on reaction rate using velocity distribution analysis, Int. J. Greenhouse Gas Control 68 (2018) 99–111 .

[15] Y. Al-Khulaifi, Q. Lin , M.J. Blunt , B. Bijeljic , Pore-scale dissolution by co2 saturated brine in a multimineral carbonate at reservoir conditions: impact of physical and chemical heterogeneity, Water Resour. Res. 55 (4) (2019) 3171–3193 .

[16] H. Erfani , V. Joekar-Niasar , R. Farajzadeh , Impact of microheterogeneity on up- scaling reactive transport in geothermal energy, ACS Earth and Space Chem- istry 3 (9) (2019) 2045–2057 .

[17] W.F. Ramirez , A.C. Oen , J.F. Strobel , J.L. Falconer , H.E. Evans , Surface composi- tion of berea sandstone, SPE Form. Eval. 2 (1986) 23–30 .

[18] W.J. Moore , Basic physical chemistry, Prentice Hall, 1983 .

[19] M. Kazempour , E. Sundstrom , V. Alvarado , Geochemical modeling and exper- imental evaluation of high-ph floods: impact of water–rock interactions in sandstone, Fuel 92 (1) (2012) 216–230 .