Dispersion of a fluid plume during radial injection in an aquifer
aa r X i v : . [ phy s i c s . f l u - dyn ] J a n Dispersion of a fluid plume during radial injection in anaquifer
Benjamin W.A. Hyatt Yuri Leonenko Department of Earth and Environmental Sciences, University of WaterlooJanuary 22, 2021
Abstract
This study outlines a model for injected fluid flow in a vertically confined porous aquiferwith mechanical dispersion. Existing studies have investigated the behaviour and geometry ofimmiscible fluid flow in this setting, where the injected fluid displaces the resident fluid, forminga sharp interface between the two. The present study extends analytical solutions to includemechanical dispersion of the interface. The solutions are inverted to solve for time as a functionof position ( 𝑟, 𝑧 ) , giving each position in the aquifer an intersection time corresponding to themoment the travelling interface intersects a point of interest. The set of { 𝑟 𝑜 , 𝑧 𝑜 } positions whichshare an intersection time are treated as dummy variables that represent an “effective surface” andare integrated over to solve for the velocity field within the aquifer. Using this velocity field, theconcentration profile resulting from mechanical dispersion can be found analytically. It is shownthat the concentration of the injected fluid smoothly decays around the position of the interfacefrom immiscible solutions, allowing for the injected fluid to be present in detectable quantitiesbeyond the extent of these interfaces. This concentration spread should be considered in definingouter boundaries on fluids in injection well applications such as carbon capture and storage orgroundwater applications. The flow of fluids injected into porous media has been a research area of interest for decades, havingwide ranging applications including enhanced oil recovery (Buckley and Leverett , 1942), assess-ing drinking water quality (Prommer and Stuyfzand , 2005) and carbon capture and storage (CCS)(Gibbins and Chalmers , 2008). This has inspired the development of many analytical solutions andnumerical simulations modelling the behaviour, geometry and evolution of fluids being injected intothe resident brine of aquifers under a variety of circumstances. A common feature throughout manyof these studies is the assumption that the injected fluid and resident fluid are immiscible; that is, theinjected fluid does not dissolve in the resident fluid and forms a sharp interface separating regionsof 100% injected fluid saturation and 100% brine saturation. Another assumption that can allow forsimple approximate solutions for the fluids’ evolution in the aquifer is vertical equilibrium, wherethe flow velocity in the vertical direction is assumed to be negligible, and the velocity field in theaquifer is strictly radial (in the direction outward from the injection site) (Nordbotten et al , 2005;Nordbotten and Celia , 2006; Juanes et al , 2010).With simplifying assumptions and boundary conditions, analytical solutions for the shape orthickness of this sharp interface as a function of radial distance and time can be obtained. This allows1or the furthest radial extent of the interface to be determined after injecting for a certain period.However, sizable traces of the injected fluid may appear beyond this boundary due to mechanicaldispersion, a mass transfer mechanism where velocity changes on the scale of the medium’s porescause the fluid to spread out (Fetter , 2001): this will cause the concentration of the injected fluidto smoothly decay around the region where this sharp boundary would be. This phenomenon hasbeen examined in studies that assume a radial flow of injected fluid with no variation in the verticaldirection (Hoopes and Harleman , 1967; Hsieh and Yeh , 2014; Tang and Babu , 1979). This purelyradial flow, however, can be interpreted as having an underlying cylindrical, immiscible interfaceradiating outward from the injection well, which itself is a special case of a general interface geometrysubject to vertical equilibrium.The present study proposes a mathematical procedure which allows including the mechanicaldispersion process into consideration and applying it to generic interface geometries, extendingbeyond the cylindrical case. This allows for quantifying a more precise upper boundary cutoff forthe distances at which traces of the injected fluid may appear in the aquifer: for example, findingwhere the relative concentration of the injected fluid drops to some acceptable levels, such as 1%,0.5%, etc. Similarly, a lower boundary distance value for a desired relative concentration cutoff(e.g., 99%) can be defined to quantify regions in the aquifer that are effectively saturated with theinjected fluid. The presence of injected contaminants is of great interest in injection well engineeringapplications (Ahmad et al , 2010; Jankovic and Fiori , 2010; Comba and Braun , 2012; Cahill et al ,2014). As such, computing these boundaries between regions saturated with groundwater, injectionfluid, or a transition zone in between will be of practical in determining injection efficiency andsafety parameters.For the purposes of this investigation, the resident and injected fluids are assumed to have velocityfields which obey Darcy’s law and are incompressible with constant densities and viscosities. Thevertical pressure throughout the aquifer is assumed to be hydrostatic, and the aquifer is assumed tobe isothermal and axisymmetric, with a finite height due to impermeable caprock layers at the topand bottom and infinite laterally. Furthermore, the aquifer is assumed to have constant permeability,porosity, and dispersivity scale.
The mathematical expression for the incompressibility condition of a fluid states that the divergenceof its velocity field is zero (White , 2011): ∇ · 𝒗 = 𝒗 ∝ 𝑟 ˆ 𝒓 (2)An expression for the velocity field can be obtained from the volumetric flow rate of the fluid 𝑄 [L /T] evaluated through some characteristic surface (elaborated on below), which is given by thetotal flux of the velocity field through that surface (White , 2011): 𝑄 = 𝜙 ∬ 𝒗 · 𝒅𝑨 (3)2ere, 𝑄 is assumed to be equal to a constant volumetric injection rate at the site of the well. Thedimensionless factor 𝜙 is the porosity of the aquifer the fluid is flowing through, which is the fractionof volume in the medium that is open and unoccupied of material. This porosity factor is introducedto relate the flow velocity with the volumetric flux from Darcy’s law (Haitjema and Anderson , 2015).For the purposes of this investigation, the porosity is assumed to be constant throughout the aquifer.The nature of the integrated over “characteristic surface” mentioned above is at the heart of thisanalysis. If the injected and resident fluids are imagined to be immiscible, this surface representsthe travelling interface between them the moment it intersects a point of interest ( 𝑟, 𝑧 ) in the aquifer.Note that this surface should depend only on the intersection position and is independent of time orwhere the “true” interface is located. Also note that for any point ( 𝑟, 𝑧 ) there is a corresponding setof points { 𝑟 𝑜 , 𝑧 𝑜 } which represent the solution for the interface at the time it crosses a point ( 𝑟, 𝑧 ) .These 𝑟 𝑜 , 𝑧 𝑜 points represent an effective surface for a point ( 𝑟, 𝑧 ) and are treated as dummy variableswhich are used to evaluate the flux integral in equation 3.The velocity field, however, is interpreted to be a function of the specific position ( 𝑟, 𝑧 ) , andnot these dummy variables representing an effective surface; as such, it will be factored out of theintegral. Additionally, as the velocity field is assumed to be strictly radial, the area element dottedwith it can be written as: 𝒗 · 𝒅𝑨 = 𝑣 ( 𝑟 𝑜 𝑑𝜃𝑑𝑧 𝑜 ) (4)Where the cylindrical element 𝑟 𝑜 𝑑𝜃𝑑𝑧 𝑜 is the projection of the 𝒅𝑨 surface element using thedummy variables mentioned above. Consider the example of a cylindrical, immiscible verticalinterface travelling through a confined aquifer of height 𝐻 . In this case, 𝑟 𝑜 = 𝑟 is independent of thevertical 𝑧 𝑜 , and evaluating the flux integral in 3 yeilds: 𝑄 = 𝜋𝜙𝑟𝑣 ∫ 𝐻 𝑑𝑧 𝑜 = 𝜋𝜙𝑟𝑣𝐻 (5)Solving for the velocity field above gives 𝑣 = 𝑄 / 𝜋𝜙𝐻 , which satisfies the incompress-ibility condition and is the same equation as the velocity given by Tang and Babu (1979) andHoopes and Harleman (1967).In general, the interface formed by the injected fluid may not have a flat, cylindrical interface;due to density and viscosity differences between the fluids, the interface may take on a curved shapewhere the injected and resident fluids are vertically segregated. These interface geometries result in 𝑟 𝑜 values which are a function of 𝑧 𝑜 and time in general. Figure 1 depicts the motion of an arbitrarilyshaped immiscible interface over a very short time interval Δ 𝑡 , travelling with radial velocity 𝑣 = 𝐴 / 𝑟 ,where 𝐴 is constant with respect to 𝑟 . Note that the dummy coordinates 𝑟 𝑜 , 𝑧 𝑜 are not used in thefigure as it is depicting the time dependent, “true” position of the interface rather than the effectiveintersection surface which intersects an arbitrary point. If the overall shape of the interface does notsignificantly change over this interval, then the interface’s height at 𝑡 + Δ 𝑡 is approximately the sameas the height of the interface at time 𝑡 but at a postion further back by Δ 𝑟 = 𝑣 Δ 𝑡 , i.e.: 𝑧 ( 𝑟, 𝑡 + Δ 𝑡 ) ≈ 𝑧 ( 𝑟 − Δ 𝑟, 𝑡 ) (6)This relationship can be used in approximating the time derivative of the interface height 𝑧 : 𝜕𝑧𝜕𝑡 ≈ 𝑧 ( 𝑟 − Δ 𝑟, 𝑡 ) − 𝑧 ( 𝑟, 𝑡 ) Δ 𝑡 = 𝑣 (cid:18) 𝑧 ( 𝑟 − Δ 𝑟, 𝑡 ) − 𝑧 ( 𝑟, 𝑡 ) Δ 𝑟 (cid:19) ≈ − 𝑣 𝜕𝑧𝜕𝑟 (7)Thus, 3 r - r z(r,t) v = A/rz(r - r,t) z(r,t + t)r (cid:1) v t Injected fluid Resident fluid rz Figure 1: An immiscible interface travelling over a short time interval Δ 𝑡𝜕𝑧𝜕𝑡 + 𝑣 𝜕𝑧𝜕𝑟 = 𝑧𝑣 / 𝑟 to equation 8, one can obtain a continuity equation for theheight of the immiscible interface with respect to distance and time: 𝜕𝑧𝜕𝑡 + 𝑟 𝜕𝜕𝑟 ( 𝑟𝑧𝑣 ) = 𝜕𝑧𝜕𝑡 − 𝑟𝜙 (cid:18) 𝑔 Δ 𝜌𝑘𝑧 ( 𝐻 − 𝑧 ) 𝑟𝜇 𝑟 𝑧 + 𝜇 𝑖 ( 𝐻 − 𝑧 ) 𝜕𝑧𝜕𝑟 + 𝑄𝜇 𝑖 ( 𝐻 − 𝑧 ) 𝜋 ( 𝜇 𝑟 𝑧 + 𝜇 𝑖 ( 𝐻 − 𝑧 )) (cid:19) = 𝜇 𝑟 and 𝜇 𝑖 terms are the viscosities of the residual and injected fluids respectively, Δ 𝜌 is theirdensity difference, 𝐻 is the height of the aquifer, 𝑘 is its permeability, and 𝑔 is the acceleration dueto gravity. This form of the equation appears in Nordbotten and Celia (2006) as well as Guo et al(2016), who derive approximate solutions for it given various fluid properties.Given an analytical expression for a solution for this governing equation (typically 𝑧 as a functionof 𝑟, 𝑡 ) and solving for the intersection time 𝑡 𝑜 at a point ( 𝑟, 𝑧 ) , one can apply the effective intersectionsurface integral in equation 3 by substituting 𝑡 𝑜 back into the solution in terms of the dummycoordinates 𝑟 𝑜 and 𝑧 𝑜 , and solving for 𝑟 𝑜 in terms of 𝑟, 𝑧, and 𝑧 𝑜 . Finally, this expression can be usedin the flux integral to obtain an expression for the velocity field, which will be of the incompressibleform: 𝒗 = 𝐴 ( 𝑧 ) 𝑟 ˆ 𝒓 (11)In general, the term 𝐴 is a function of 𝑧 characteristic of the analytical interface solution. In thecase of the cylindrical interface, 𝐴 is a constant equal to 𝑄 / 𝜋𝜙𝐻 .4n reality, the interface between injected and residual fluids is never sharp. There is alwaysdispersion of the fluids due to mass transfer effects taking place at the interface. Such processes aregoverned by the advection-diffusion equation (ADE), which has the general form (Stocker , 2011): 𝜕𝑐𝜕𝑡 = ∇ · ( 𝑫 ∇ 𝑐 ) − ∇ · ( 𝒗 𝑐 ) + 𝑃 (12)Here, 𝑐 is the fluid concentration, 𝑫 is the diffusivity (a tensor, in general; see Bear (1961))and 𝑃 represents sources and sinks. In the domain of the aquifer outside of the injection site, itis assumed that 𝑃 =
0. The diffusivity term is taken to be the sum of the molecular diffusion andmechanical dispersion effects, as done by Neuman et al (1987): 𝐷 𝑇 = 𝐷 𝑚 + 𝑑 𝑇 k 𝒗 k , 𝐷 𝐿 = 𝐷 𝑚 + 𝑑 𝐿 k 𝒗 k (13)Here, 𝐷 𝑚 is the molecular diffusion coefficient (assumed constant), 𝐷 𝑇 and 𝐷 𝐿 are the transverse(normal to the velocity field) and longitudinal (parallel to the velocity field) diffusivity componentsrespectively, and 𝑑 𝑇 and 𝑑 𝐿 are the transverse and longitudinal dispersivity scales. In the presentstudy’s case of a strictly radial velocity field in an isotropic medium, the transverse dispersivity istaken to be 0, and the longitudinal dispersivity 𝑑 𝐿 is simply referred to as the dispersivty 𝑑 . It hasbeen suggested that in large field scale transport situations, longitudinal dispersivity approaches aconstant or asymptotic value at larger distances (Pickens and Grisak , 1981). The asymptotic value(which is constant) can serve as a worst case scenario for evaluation of dispersion. Furthermore, ifthe effects of molecular diffusion are assumed to be far less than those of mechanical dispersion (avery good assumption for all practical applications), the ADE can be simplified as: 𝜕𝑐𝜕𝑡 + 𝐴𝑟 (cid:18) 𝜕𝑐𝜕𝑟 − 𝑑 𝜕 𝑐𝜕 𝑟 (cid:19) = 𝑐 ( 𝑟, 𝑡 ) = 𝑐 𝑜 (cid:16) ln (cid:0) 𝑟𝑑 (cid:1) − ln (cid:16) 𝑑 p 𝐴𝑡 + 𝑟 𝑜 (cid:17) (cid:17) (cid:0) 𝐴𝑡 + 𝑟 𝑜 (cid:1) 𝑑 r 𝑑 (cid:16) ( 𝐴𝑡 + 𝑟 𝑜 ) − 𝑟 𝑜 (cid:17) (15)In this equation, 𝑟 𝑜 is the radius of the injection well, 𝑐 𝑜 is the initial concentration of the injectedfluid, and 𝐴 is the term from the incompressible velocity field. For an aquifer of infinite radial extent,the injection well radius 𝑟 𝑜 can be neglected, giving the so called “line solution” as done by Guo et al(2016). The above equation can also be solved for 𝑟 as a function of time and concentration, whichis of practical interest and allows for comparing the “cut off” for the extent of injected fluid presencewith an immiscible interface. After finding the velocity field for general interface geometries, the 𝐴 ( 𝑧 ) term can be substituted into this equation to allow for 𝑧 dependent solutions: this procedure iscarried out below for some analytical interface solutions. To demonstrate this transformation, three approximate immiscible interface solutions for 𝑧 as afunction of 𝑟, 𝑡 from Guo et al (2016) have been chosen for its application. For all these solutions,5he injected fluid is assumed to be less dense than the resident fluid, while the driving effects ofbuoyancy are assumed to be far less than those of injection, with viscosity differences between theinjected and resident fluids being the key factor in each solution’s geometry. All three solutionsare fixed to be 𝑧 = 𝑧 = 𝐻 outside of two moving boundaries characteristic of each solution.For the intersection surface used in the transformation, only the solution in the region within theseboundaries is used as it represents the travelling interface. Under these circumstances, Guo et al (2016) provides the following “travelling interface” heightsolution: 𝑧 ( 𝑟, 𝑡 ) = 𝐻 ( 𝑀 − ) 𝑀 Γ (cid:18) 𝜋𝜙𝐻𝑟 𝑄𝑡 − (cid:19) + 𝐻 , − 𝑀 − 𝑀 Γ < 𝜋𝜙𝐻𝑟 𝑄𝑡 ≤ + 𝑀 − 𝑀 Γ (16)Here, the dimensionless parameters Γ = 𝜋 Δ 𝜌𝑔𝑘𝐻 / 𝜇 𝑟 𝑄 represents the effect of buoyancycompared to injection and 𝑀 = 𝜇 𝑟 / 𝜇 𝑖 is the ratio of the resident and injected fluids’ viscosities.Mathematically, the solution in 16 assumes that Γ ≪ 𝑀 < 𝑡 = 𝑡 𝑜 be the time a specific point ( 𝑟, 𝑧 ) is intersected by the travelling interface. Rearranging16 for this time yields: 𝑡 𝑜 = 𝜋𝜙𝐻𝑟 𝑄 (cid:16) + 𝑀 Γ 𝑀 − (cid:16) 𝑧𝐻 − (cid:17) (cid:17) (17)Substituting this intersection time back into the solution with respect to the dummy variables: 𝑧 𝑜 = 𝐻 ( 𝑀 − ) 𝑀 Γ (cid:18) 𝑟 𝑜 𝑟 (cid:18) + 𝑀 Γ 𝑀 − (cid:18) 𝑧𝐻 − (cid:19) (cid:19) − (cid:19) + 𝐻 𝑟 𝑜 yields: 𝑟 𝑜 = 𝑟 vuuuut + 𝑀 Γ 𝑀 − (cid:16) 𝑧 𝑜 𝐻 − (cid:17) + 𝑀 Γ 𝑀 − (cid:16) 𝑧𝐻 − (cid:17) (19)Substituting this dummy variable 𝑟 𝑜 as a function of 𝑟, 𝑧 and 𝑧 𝑜 into 4 and then into the fluxintegral (3): 𝑄 = 𝜋𝜙𝑣𝑟 r + 𝑀 Γ 𝑀 − (cid:16) 𝑧𝐻 − (cid:17) ∫ 𝐻 s + 𝑀 Γ 𝑀 − (cid:18) 𝑧 𝑜 𝐻 − (cid:19) 𝑑𝑧 𝑜 (20)Finally, evaluating this integral and solving for the velocity field yields: 𝑣 ( 𝑟, 𝑧 ) = 𝑄𝑀 Γ r + 𝑀 Γ 𝑀 − (cid:16) 𝑧𝐻 − (cid:17) 𝜋𝜙𝑟𝐻 ( 𝑀 − ) (cid:16) (cid:0) + 𝑀 Γ 𝑀 − (cid:1) − (cid:0) − 𝑀 Γ 𝑀 − (cid:1) (cid:17) (21)As expected, this velocity field is of the incompressible form:6 = 𝐴 ( 𝑧 ) 𝑟 ˆ 𝒓 , 𝐴 ( 𝑧 ) = 𝑄𝑀 Γ r + 𝑀 Γ 𝑀 − (cid:16) 𝑧𝐻 − (cid:17) 𝜋𝜙𝐻 ( 𝑀 − ) (cid:16) (cid:0) + 𝑀 Γ 𝑀 − (cid:1) − (cid:0) − 𝑀 Γ 𝑀 − (cid:1) (cid:17) (22)Using this velocity field, the 2-D concentration profile can be given by substituting 𝐴 ( 𝑧 ) above intothe equation from Dagan (1971). Similarly, one can also rearrange the concentration equation for 𝑟 ,and find the radial extent of boundaries that are a function of a specific concentration. Two boundariesof practical interest are where 𝑐 / 𝑐 𝑜 = .
99, and 𝑐 / 𝑐 𝑜 = .
01, as these provide an approximation for a“transition zone” where the relative concentration of the injected fluid experiences the most variationdue to dispersion; outside of these boundaries, the aquifer is essentially saturated with the residentor the injected fluid. The inverted concentration equation (including the line solution approximation 𝑟 𝑜 ≈
0) has the form: 𝑟 ( 𝑐, 𝑧, 𝑡 ) = p 𝐴 ( 𝑧 ) 𝑡 exp ( ( 𝐴 ( 𝑧 ) 𝑡 ) r 𝑑 − (cid:26) 𝑐𝑐 𝑜 (cid:27)) (23)To demonstrate the effect of mechanical dispersion, consider the radial extent of the 𝑐 / 𝑐 𝑜 = . 𝑧 = 𝑟 and 𝑟 respectively. Furthermore, let the function 𝑓 ( 𝑀, Γ ) be the collected “ 𝑀 terms” from 𝐴 ( ) in 22,i.e.: 𝐴 ( ) = 𝑄 Γ 𝜋𝜙𝐻 𝑓 ( 𝑀, Γ ) , 𝑓 ( 𝑀, Γ ) = 𝑀 q − 𝑀 Γ 𝑀 − ( 𝑀 − ) (cid:16) (cid:0) + 𝑀 Γ 𝑀 − (cid:1) − (cid:0) − 𝑀 Γ 𝑀 − (cid:1) (cid:17) (24)This allows the term 𝑟 to be expressed as: 𝑟 = s 𝑄 Γ 𝑡 𝑓 ( 𝑀, Γ ) 𝜋𝜙𝐻 exp (cid:16) 𝑄 Γ 𝑡 𝑓 ( 𝑀, Γ ) 𝜋 𝜙𝐻 (cid:17) r 𝑑 − { . } (25)Furthermore, let the dimensionless parameters for time and radial distance be 𝑇 and 𝑅 respectively,defined as: 𝑇 = 𝑄𝑡𝑑 𝜋𝜙𝐻 , 𝑅 = 𝑟𝑑 (26)Using these parameters, equation 25 can be nondimensionalized as: 𝑅 = p 𝑇 Γ 𝑓 ( 𝑀, Γ ) exp (cid:26) − { . }√ ( 𝑇 Γ 𝑓 ( 𝑀, Γ )) − (cid:27) (27)The term 𝑟 can be found by solving the upper bound on the domain in 16 for 𝑟 : 𝑟 = s 𝑄𝑡𝜋𝜙𝐻 (cid:18) + 𝑀 Γ − 𝑀 (cid:19) (28)Which can be similarly nondimensionalized to obtain the expression:7 = s 𝑇 (cid:18) + 𝑀 Γ − 𝑀 (cid:19) (29)Finally, the dimensionless difference between the 1% concentration boundary and the 𝑀 < 𝑧 = 𝑅 − 𝑅 = √ 𝑇 © «p Γ 𝑓 ( 𝑀, Γ ) exp (cid:26) − { . }√ ( 𝑇 Γ 𝑓 ( 𝑀, Γ )) − (cid:27) − s(cid:18) + 𝑀 Γ − 𝑀 (cid:19)ª®¬ (30)Note that the function 𝑅 − 𝑅 in 30 is not well defined on the whole interval 𝑀𝜖 ( , ) . Thefucntion 𝑓 ( 𝑀, Γ ) = 𝑀 = /( + Γ ) , which causes the exponential term in 𝑅 − 𝑅 to diverge.Thus, the domain of 30 is 𝑀𝜖 ( , /( + Γ )) , the upper boundary of which approaches 1 in the limitof Γ → Γ = .
05 and 𝑇 = ,
100 and 1000. The dimensionless separation functionis always positive, and though its time derivative is initially negative, the separation increases withtime after 𝑇 ≈
2; this increased separation with time is illustrated in the plot.
The approximate solution provided by Guo et al (2016) for a more viscous resident fluid ( 𝑀 > ) is: 𝑧 ( 𝑟, 𝑡 ) = 𝐻𝑀 − s 𝑄𝑀𝑡𝜋𝜙𝐻𝑟 − 𝐻, 𝑄𝜋𝜙𝐻 𝑀 < 𝑟 𝑡 ≤ 𝑄𝑀𝜋𝜙𝐻 (31)This produces the 𝑟 𝑜 expression: 𝑟 𝑜 ( 𝑟, 𝑧, 𝑧 𝑜 ) = 𝑟 + 𝑧 ( 𝑀 − ) 𝐻 + 𝑧 𝑜 ( 𝑀 − ) 𝐻 ! (32)Evaluating the flux integral with this 𝑟 𝑜 expression yields the velocity field: 𝒗 = 𝐴 ( 𝑧 ) 𝑟 ˆ 𝒓 , 𝐴 ( 𝑧 ) = 𝑄 ( 𝑀 − ) 𝜋𝜙 ( 𝑧 ( 𝑀 − ) + 𝐻 ) ln 𝑀 (33)Once again, consider the radial extents (at the top of the aquifer) of the 1% relative concentrationboundary ( 𝑟 ) and the immiscible solution ( 𝑟 ) for the 𝑀 > 𝑧 = 𝑟 : 𝑟 = 𝑑 s 𝑄𝑡 ( 𝑀 − ) 𝜋𝜙𝐻 ln 𝑀 exp ( − ( . )√ (cid:18) 𝑄𝑡 ( 𝑀 − ) 𝜋𝜙𝐻 ln 𝑀 (cid:19) − ) (34)The upper boundary on the domain of the solution 31 can be solved again to find 𝑟 : 𝑟 = 𝑑 s 𝑄𝑀𝑡𝜋𝜙𝐻 (35)8 0 Y L V F R V L W \ U D W L R 5 í 5 5 í 5 0 ī 7 7 7 Figure 2: The dimensionless difference between the 1% relative concentration boundary and theimmiscible interface solution at the top of the aquifer ( 𝑀 < ) 𝑇 and 𝑅 as defined above gives thedimensionless separation for the 1% boundary and the immiscible solution at 𝑧 = 𝑀 > 𝑅 − 𝑅 = √ 𝑇 r 𝑀 − 𝑀 exp ( − ( . )√ (cid:18) 𝑇 ( 𝑀 − ) ln 𝑀 (cid:19) − ) − √ 𝑀 ! (36)This separation function only gives physical solutions on a restricted domain. For certain 𝑀 and 𝑇 values, equation 36 will be ≤
0, which would correspond to a concentration boundary laggingbehind the underlying immiscible solution. Additionally, the derivative of 36 with respect to 𝑇 hasa root of its own at a time dependent 𝑀 value between 1 and the original function’s root. While 𝑀 values in between these roots correspond to a 1% boundary that lies ahead of the original interface,this case results in the interface converging toward and then exceeding the concentration boundary.As such, realistic solutions for the 𝑀 > 𝑀 = 𝑀 =
1, thus solutions that occur before earlier 𝑇 values have larger 𝑀 domains. The 𝑀 value of this root at a certain 𝑇 can be found numerically. Forexample, cases that occur up to 𝑇 =
10 have physical solutions approximately on 𝑀𝜖 ( , . ) , whilecases up to 𝑇 = 𝑀𝜖 ( , . ) ; the latter domain is used to plot the dimensionlessseparation of the 1% boundary and the original interface over time in figure 3. As a final example, consider the case for fluids of equal viscosity. The approximate solution providedby Guo et al (2016) is: 𝑧 ( 𝑟, 𝑡 ) = 𝐻 (cid:18) − √ Γ (cid:18) 𝜋𝜙𝐻𝑟 𝑄𝑡 − (cid:19)(cid:19) , 𝑄𝜋𝜙𝐻 ( − √ Γ ) < 𝑟 𝑡 ≤ 𝑄𝜋𝜙𝐻 ( + √ Γ ) (37)As before, solving for the intersection time 𝑡 𝑜 to get the dummy variable 𝑟 𝑜 yields: 𝑟 𝑜 ( 𝑟, 𝑧, 𝑧 𝑜 ) = 𝑟 vuuuut + √ Γ (cid:16) − 𝑧 𝑜 𝐻 (cid:17) + √ Γ (cid:16) − 𝑧𝐻 (cid:17) (38)Carrying out the same flux integral as before provides the velocity field: 𝒗 = 𝐴 ( 𝑧 ) 𝑟 ˆ 𝒓 , 𝐴 ( 𝑧 ) = 𝑄 r Γ (cid:16) + √ Γ (cid:16) − 𝑧𝐻 (cid:17) (cid:17) 𝜋𝜙𝐻 (cid:18) (cid:16) + √ Γ (cid:17) − (cid:16) − √ Γ (cid:17) (cid:19) (39)Evaluating the function 𝐴 ( 𝑧 ) at 𝑧 = 𝑔 ( Γ ) , written as: 𝐴 ( ) = 𝑄𝜋𝜙𝐻 𝑔 ( Γ ) , 𝑔 ( Γ ) = q Γ ( + √ Γ )( + √ Γ ) − ( − √ Γ ) (40)Finally, defining the dimensionless quantites quantities 𝑅 and 𝑅 the same way as before, the 𝑀 = 0 Y L V F R V L W \ U D W L R 5 í 5 5 í 5 0 ! 7 7 7 Figure 3: The dimensionless difference between the 1% relative concentration boundary and theimmiscible interface solution at the top of the aquifer ( 𝑀 > ) ī E X R \ D Q F \ S D U D P H W H U 5 í 5 5 í 5 0 7 7 7 Figure 4: The dimensionless difference between the 1% relative concentration boundary and theimmiscible interface solution at the top of the aquifer ( 𝑀 = ) 𝑅 − 𝑅 = √ 𝑇 p 𝑔 ( Γ ) exp ( − ( . )√ ( 𝑇𝑔 ( Γ )) ) − q + √ Γ ! (41)This separation function is well defined for Γ 𝜖 ( , ) and is always greater than zero. Like the 𝑀 < 𝑇 ≈ Γ ≪ Γ , being 𝑇 ≈ .
009 for
Γ = .
01, and 𝑇 ≈ .
994 for
Γ = . 𝑀 = leakage. In this aproach, instead of injecting pure CO , it is dissolved in brineproduced from a target aquifer and reinjected back into the formation (Cao et al , 2020, 2021).This problem’s formulation becomes identical to injecting contaminant fluid in porous media, butthe importance of despersion is increased due to much larger scales of injection. Typical aquifer12ariable Quantity [Dimensions] Value 𝑄 Volume injection rate [L T − ] 1,000,000 m /year 𝐻 Aquifer height [L] 100 m 𝑀 Viscosity ratio (res./inj. fluid) 1 Γ Buoyancy parameter 0.05 𝜙 Aquifer Porosity 0.1 𝑡 Time elapsed [T] 50 years 𝑑 Dispersivity scale [L] 35 mTable 1: List of CCS quantities used to plot physical interface and transition boundaries.99% boundary Interface solution 1% boundaryAt 𝑧 = 𝑟 =
976 m 𝑟 = 𝑟 = 𝑧 =
100 m (bottom): 𝑟 =
855 m 𝑟 = 𝑟 = 𝑀 = 𝑐 / 𝑐 𝑜 = .
99 and 𝑐 / 𝑐 𝑜 = .
01) comparedto the 𝑀 = It is clear that mechanical dispersion can allow for an injected fluid concentration transition zone thatis spread considerably on either side of an immiscible interface solution. For the three examples usedthroughout section 3, the 1% relative concentration boundary location was found to increasinglyrecede from the original interface solution over time, subject to domain restrictions for physicalsolutions. This is consistent with the interpretation of immiscible interfaces acting as effectivesurfaces underlying a dispersive concentration profile.The three immiscible interface cases examined were solutions from Guo et al (2016), where theviscosity differences between the injected fluid and resident fluid were the primary factor determiningthe interfaces’ geometries. For all three cases, the dimensionless difference between the 1% boundaryand the original interface at the top of the aquifer 𝑅 − 𝑅 is found to lie approximately in therange 6-13, with larger differences occuring at later dimensionless times 𝑇 . This corresponds to 1%concentration traces of the injected fluid for a given interface geometry lying 6-13 dispersivity lengthsahead of the immiscible interface solution. This observation is also consistent with the “effectivesurface” interpretation of the concentration profile, as the 1% relative concentration boundary ofinjected fluid due to mechanical dispersion is found to propagate further than the furthest radial13 U P ] P , Q W H U I D F H D Q G F R Q F H Q W U D W L R Q E R X Q G D U L H V 0 , Q W H U I D F H L P P L V F V R O X W L R Q &