Modeling the Closure of Thief Zones in Carbonate Reservoirs Using

10080

Ind. Eng. Chem. Res. 2010, 49, 10080–10093

Modeling the Closure of Thief Zones in Carbonate Reservoirs Using Acids Hamidreza Salimi,* Karl-Heinz Wolf, and Johannes Bruining Department of Geotechnology, Delft UniVersity of Technology, SteVinweg 1, 2628 CN Delft, The Netherlands

Carbonate oil reservoirs can contain fractures and dissolution-vug pore space that may cause short-circuit flows between injection and production well. As a result, injected water bypasses the larger part of the reservoir, leaving much oil unproduced. This paper investigates the concept of reducing the permeability in the fracture that is causing the short-circuit flow by injecting a mixture of sulfuric acid and hydrochloric acid. During the process, reservoir rock dissolves and gypsum subsequently precipitates. We formulate a mathematical and numerical model that simulates the coinjection of sulfuric acid (H2SO4) with hydrochloric acid (HCl) into carbonate rock for both one- and two-dimensional settings. The carbonate rock contains a single fracture. A one-dimensional simulation of the process shows that upstream dissolution is the dominant process, whereas precipitation of gypsum occurs downstream. The two-dimensional simulation shows that downstream fractures are clogged, whereas continuous paths are created upstream. We find that the injection flow rate and acid ratio are the two main factors that influence precipitation and dissolution patterns within a fracture. The results show qualitative agreement with experimental results obtained in the Dietz-De Josseling laboratory. We conclude that injection of an acid mixture can be used to improve water-drive recovery performance from oil reservoirs that contain fractures between injection and production well. 1. Introduction In petroleum engineering, acids are routinely injected to improve the permeability around the production well. Hydrochloric acid is used for well treatments in various formations. In these treatments, oil and gas wells are stimulated by a matrix acidizing treatment or an acid-fracturing treatment.1-4 In matrix acidizing treatments, acid is used to improve conductivity near the wellbore. However, acids are not contemplated for closing fractures. Araque-Martinez and Lake5 applied the method of characteristics to radial flow around a well to study the effect of dissolution, supersaturation, and precipitation around injection wells. Chemical precipitation can cause a reduction in permeability during the injection or mixing of ionic solutions in a permeable medium. The permeability reduction has been attributed to the migration and plugging of pore throats by solid precipitates. Wu and Sharma6 proposed a new model to characterize the local permeability reduction in terms of precipitated solids-size and pore-size distributions into one of the two dominant pore plugging (or solid trapping) mechanisms: particle-size exclusion by pore throats. Noh7 also presented a mathematical model to simulate hydrodynamics and fluidmineral reaction in a fracture within permeable media as a simplified representation of natural fracture mineralization. He7 showed the time evolution of fracture aperture shrinkage patterns caused by calcite cementation. A chalk (calcium carbonate or CaCO3) reservoir can contain fractures and vugs that cause short-circuit streams between injection and production well, which, in turn, lead to the bypassing of a considerable amount of oil. Vugs include all centimeter-scale and larger pores cavities, caverns, and conduits that were developed by the dissolution of carbonate rock.8 Schuiling9 proposed the use of sulfuric acid to increase the volume of a carbonate layer via a conversion reaction to, among others, “volume-increasing” gypsum, to compensate for subsidence, thus leading to uplifting of the land surface. Following his idea, Maersk Oil (see acknowledgment) proposed the use * To whom correspondence should be addressed. Tel.: +31-152788065. Fax: +31-15-2781189. E-mail: [email protected].

of this concept to close the fracture system. Hence, injected water is diverted and improves the sweep efficiency from the matrix part (see Figure 1). However, the injection of sulfuric acid rapidly clogs the pore space, precluding that the flow can penetrate deeper into the reservoir. Therefore, we propose here to inject a mixture of acids (sulfuric acid (H2SO4) and hydrochloric acid (HCl)) to avoid this problem. This study is focused on the development of a model that can be used to describe this process in the laboratory and in the field. The proposed method may also work in heterogeneous reservoirs when thief zones lead to short-circuit streams. Hoefner and Fogler10 described the phenomenon of wormhole formation after hydrochloric acid injection. They proposed the following mechanism. Because the medium is heterogeneous, there are regions with higher permeability. These regions will have more H+ ions flowing to them than to their less-permeable surroundings. Therefore, in these regions, a higher rate of dissolution of calcium will take place. This increases the permeability even more and therefore attracts more H+ ions to the highly permeable regions. In the end, this can create a channel that goes through the entire sample. Sherwood11 derived a simple model that demonstrated the effect of the Damko¨hler number and of the acid capacity number on the instability of a plane reaction front. Indeed, for HCl injection, the main destabilizing mechanism is that the permeability upstream is higher than the permeability downstream. For H2SO4 injection, the permeability upstream becomes lower than the permeability downstream, and this has a stabilizing effect. Numerous papers have described a normalmode analysis of dissolution for simplified cases.4,10-13 However, no paper has addressed the combined effects of dissolution and precipitation simultaneously. We also note that, in the situation considered here, the unperturbed states are timedependent. State-of-the-art normal-mode analysis has not been developed yet to deal with such complex cases (i.e., dissolution and precipitation at the same time with variation of the porosity, density, and reaction rate). Dobson et al.14 performed a laboratory experiment and onedimensional numerical simulations to investigate mineral dis-

10.1021/ie100616x  2010 American Chemical Society Published on Web 09/23/2010

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010

10081

Figure 1. Top view of a reservoir with a horizontal production and an injection well: carbonate reservoir between two wells containing a fracture (left), starting injection of the acid mixture (middle), and blocking of the fracture caused by precipitation of gypsum and, hence, water flow diversion to lesspermeable reservoir parts (right).

solution and precipitation during water flow through a planar fracture under given temperature and pressure conditions in the repository. Singurindy and Berkowitz15 studied the effect of injecting acid mixtures of HCl and H2SO4 in chalks. Their main conclusion is that the precipitation and dissolution mechanisms are dependent on the H+/SO42- ratio and flow rate. Consecutive experimental work of Singurindy and Berkowitz16 studied the effect of injecting acids in carbonate samples with fractures. They used acid mixtures of 0.1 M HCl/0.1 M H2SO4 and 0.1 M HCl/0.3 M H2SO4 for injection into the samples. Their conclusion was that the orientation of the fracture, with respect to the injected-flow direction, has a large influence on precipitation. According to Singurindy and Berkowitz,16 clogging occurs in fractures that are not parallel to the flow direction; however, for parallel fractures, it leads to a dissolution-dominated system. Speck et al.17 also have described an experimental study of simultaneous HCl and H2SO4 injection. They concluded that HCl forms wormholes (highly permeable paths), while H2SO4 reduces the permeability. In another paper, Dijk and Berkowitz18 derived a model of acid injection by considering precipitation and dissolution of a single component in a single fracture and obtained the aperture development. We were not able to find a reference where the reactions with HCl and H2SO4 in carbonate reservoirs were treated simultaneously and independently. Therefore, we formulate the transport equations for the H2SO4/ HCl injection process, dealing with nine dependent variables. Moreover, we describe some preliminary experiments in our laboratory where acid flow was forced parallel to the sample axis to study the development of wormholes (highly permeable paths) and clogging of initially fractured zones. The paper is organized as follows: Section 2 describes the general framework of a physical model for sealing highly permeable thief zones in chalk (CaCO3) reservoirs using acid. Here, we also define the reaction rates and a volume-of-mixture relation. Section 3 gives the full one-dimensional and twodimensional numerical model. Section 4 explains the onedimensional numerical implementation, whereas section 5 describes the two-dimensional implementation. Section 6 discusses the one-dimensional and two-dimensional results and makes a comparison with the preliminary experimental results. We end with some conclusions.

2. Physical Model We consider a single fracture that intersects with the injection and production wells. This configuration is chosen because it mimics experiments that are currently carried out in the laboratory. The laboratory uses a cylindrical core with length L and diameter H as a part of the reservoir. A vertical fracture of width δ, length L, and height H divides this block in two equal half-cylinders. For two-dimensional simulations, we use a portion of this block. For this, we imagine that the block is cut horizontally in the middle, leaving a surface with a fracture width δ and length L in the center of a rectangle of width W and length L, where L > W. This is the geometry of our twodimensional model (see Figure 1). In the one-dimensional model, we only consider behavior inside the fracture and ignore the presence of the surrounding matrix. The fracture is also considered as a porous medium with porosity φ and permeability K. We contemplate a fractured region with initial porosity φi. This means that the fraction of the volume in the fractured region arbitrarily occupied with CaCO3 is equal to (1 - φi). For instance, all CaCO3 can be confined to the fracture walls, but it is included in the fractured region. CaCO3 has a density of FC and gypsum has a density of FS. If CaCO3 reacts with H2SO4, it is converted to gypsum (i.e., CaSO4 · 2H2O), which will occupy a fractional volume of (1 φi) × FC/FS. The maximum fractional volume of the solids is equal to one (i.e., the porosity is zero), meaning that the maximum porosity for which complete clogging occurs satisfies the equation FC

(1 - φi) × F ) 1 S This equation leads to φi ) 0.59 if we use the densities FC and FS that are given in Table 1. We note that this value is obtained considering only H2SO4 reaction; however, in our situation, HCl also consumes some of the CaCO3. Hence, we use a value of φ ) 0.50, to allow for the HCl reaction. This description is used both for one-dimensional and two-dimensional modeling. 2.1. Model Description. When H2SO4 and HCl are injected in a chalk (CaCO3) reservoir, reactions occur. We ignore the presence of any impurities in the chalk (such as, e.g., Ca(OH)2). HCl forms soluble CaCl2, dissolved CO2, and water as reaction products. In the model, we disregard the various carbonate ions.

10082


Table 1. General Input Data Used in the Simulations input

value

reaction rate constant of H2SO4, kS reaction rate constant of HCl, kC concentration of H2SO4, c°SO concentration of HCl, c°Cl concentration of CaCl2, c°Ca concentration of H2O, c°OH concentration of CO2, c°CO density of CaCO3, FC density of gypsum, FS water viscosity, µ average matrix porosity, φ fracture porosity, φ number of grid blocks, Ny × Nx pressure at the outlet, p i initial matrix permeability, K m initial fracture permeability, K fi specific surface area, a

2.8 × 10-5 cm/s 2.8 × 10-3 cm/s 0.0187 mol/cm3 0.0329 mol/cm3 0.0194 mol/cm3 0.0556 mol/cm3 0.0227 mol/cm3 0.0271 mol/cm3 0.011 mol/cm3 0.001 Pa s 0.30 0.50 30 × 40 50 bar 18 mD 600 mD 200 cm-1

In the same way, H2SO4 forms solid CaSO4 · 2H2O, dissolved CO2, and water as reaction products. We will ignore the fact that, under some temperature and pressure conditions, anhydrite (CaSO4 without water) is formed. The reaction rates are assumed to be proportional to the reactant surface area of CaCO3 and the acid concentrations, i.e., of HCl and H2SO4. We disregard minor soluble amounts of CaSO4. The CaSO4 · 2H2O has a volume that is roughly twice the volume of CaCO3 per mole of calcium. However, we disregard any structural changes of the porous skeleton, because of the reactions, and, therefore, the excess volume will decrease the pore space and, because of that decrease, the porosity. This affects the permeability, which is assumed to follow the porosity dependency of the CarmanKozeny relation. This leads to a heterogeneous permeability field. The dissolved components used here are H2SO4, HCl, CaCl2, H2O, and CO2. We assume that no gaseous CO2 is formed, which is validated a posteriori by showing that the CO2 concentration remains below the solubility limit. Moreover, we assume that the total aqueous volume is a sum of the individual volumes of the constituents. We disregard any temperature variations caused by the reactions. We use the subscripts SO, Cl, Ca, HO, and CO to denote H2SO4, HCl, CaCl2, H2O, and CO2, respectively. In addition, we use the subscripts “S” and “C” to denote CaSO4 · 2H2O and CaCO3, respectively. The transport of the soluble components is governed by the convection-reaction equation. We disregard diffusion and dispersion. Because we are solving the equations numerically, we inevitably introduce numerical diffusion. The CaCO3 reacts with H2SO4, which leads to the formation and precipitation of CaSO4 · 2H2O and the disappearance of CaCO3. Single (aqueous)-phase Darcy’s law describes the flow. We disregard the density heterogeneity in Darcy’s equation. We consider the density variations in the accumulation term of the convectionreaction equations. It is assumed that these density variations are only dependent on the concentrations, because the pressure dependency of the density is ignored. We also ignore the effect of dissolved components on the viscosity of the aqueous solution. An initial heterogeneous porosity distribution is included to start the local formation of wormholes. 2.2. Model Equations. The derivation of the one-dimensional or two-dimensional equations can be made identical by using the gradient and divergence operators appropriately. The core contains an initial volume fraction νCi of solid CaCO3. On the left side, a mixture of H2SO4 and HCl is injected with a i i concentration of cSO and cCl , respectively. As a result, some of the CaCO3 is converted to insoluble CaSO4 · 2H2O and some is converted to soluble CaCl2. The volume fraction of CaSO4 ·

2H2O is denoted as νS and the volume fraction of CaCO3 is denoted as νC. The porosity (φ) is given by the expression φ ) 1 - ν S - νC

(2.1)

The initial permeability is defined as Ki. The permeability of the sample at any time during the process can be expressed, using the Carman-Kozeny expression, as K ) Ki

[

φ3 (1 - φ)2

]

(2.2)

The reactions considered are CaCO3 + H2SO4 + H2O f CaSO4 · 2H2O(s) + CO2(aq) (2.3) and CaCO3 + 2HCl f CaCl2 + H2O + CO2(aq)

(2.4)

From the results, we can validate that no gaseous CO2 is formed, because CO2 can be accommodated in the aqueous phase (i.e., it will not exceed the solubility limit and form a gaseous phase). Indeed, following Henry’s law, one would predict that, at 50 bar, 1.7 mol/L is dissolved in water, which is much more than the injected acid concentrations. The water solution containing H2SO4, HCl, CaCl2, and CO2 will be transported according to Darcy’s law of single-phase flow. 2.2.1. Volume of Mixture Relation. We assume that the solution is an ideal solution. The term “ideal solution” means that, in an ideal mixture, the volume is the addition of the volumes of its components. The “molar concentrations” of the components are denoted by c°SO, etc., in obvious notation. For an ideal solution, we can write cCl cCa cCO cOH cSO + + + + )1 c°SO c°Cl c°Ca c°CO c°OH

(2.5)

where ci/c°i is equal to the volume fraction of component i. Equation 2.5 implies that the fluid density is the summation of the concentrations of its components. 2.2.2. Reaction Rates. The reaction rates, expressed in units of mol/(m3 reactor/s), are denoted by rS″″ for the reaction that is described by eq 2.3 and by rC″″ for the reaction that is described by eq 2.4. By convention, these rates are denoted using four primes as superscripts.19 We use the simplest possible expressions for the reaction rates and ignore supersaturation conditions or chemically induced reactive surface changes. We assume that these reaction rates satisfy the following equations:4,20 r″″ S ) -kSaVCcSO

(2.6)

r″″ C ) -kCaVCcCl

(2.7)

and

where the rates are referring to the molar rates of consumption of CaCO3. The reaction rate constants are denoted by kS and kC (expressed in units of cm/s), respectively. The reaction rates are proportional to the specific surface area of CaCO3, i.e., aVC. The transport equations in the fluid phase only take into account convection and reaction and ignore diffusion; they are expressed as ∂ (φcSO) + ∇ · (uwcSO) ) r″″ S ∂t

(2.8)


∂ (φcCl) + ∇ · (uwcCl) ) 2r″″ C ∂t

(2.9)

∂ (φcCa) + ∇ · (uwcCa) ) -r″″ C ∂t

(2.10)

∂ (φcOH) + ∇ · (uwcOH) ) r″″ S - r″″ C ∂t

(2.11)

∂ (φcCO) + ∇ · (uwcCO) ) -r″″ S - r″″ C ∂t

(2.12)

where φ is the porosity of the porous medium and uw is the water-phase (aqueous) Darcy velocity vector. We consider the density variations in the accumulation term of eqs 2.8-2.12. When the concentration (molar density) of each component changes during the simulation, the fluid density (total concentration) would also change accordingly. The rates for the formation of the solid components are expressed as ∂ 1 V ) (r″″ + r″″ S) ∂t C FC C

(2.13)

()

∂ 1 V )r″″ ∂t S FS S

(2.14)

where νC and νS are the CaCO3 and CaSO4 · 2H2O volume fractions, respectively, and FC and FS are the molar densities of CaCO3 and CaSO4 · 2H2O. 2.2.3. Initial and Boundary Conditions. The initial conditions of all dissolved components for the one-dimensional (r ) x) and two-dimensional (r ) x, y) model are given as cSO(r,t ) 0) ) cCl(r,t ) 0) ) cCa(r,t ) 0) ) cCO(r,t ) 0) ) 0 (2.15)

tion gradients would also be zero to describe a zero-component flux condition, because the normal flow velocities at these boundaries are zero. In the one-dimensional model, we solve all partial-differential equations only in the fracture domain; therefore, the initial porosity equals 0.50 in the entire domain. Because there are two media in the two-dimensional model, i.e., the fracture and reservoir matrix, the initial porosity is dependent on the space coordinate. In the two-dimensional model, we introduce the heterogeneity in the reservoir rock domain by assigning different porosity values to different grid cells. The porosity values of the reservoir rock are generated by adding a random number uniformly distributed in the interval [-∆φ, ∆φ] to the mean value of the reservoir matrix porosity (i.e., 0.3). For our calculations, we use a value of ∆φ ) 0.1. The boundary conditions at the inlet boundary are uwxcl (x ) 0, y,t) ) 0

(

VC(r,t ) 0) ) ViC ) 1 - φi VS(r,t ) 0) ) 0

(uwycl)(x, y ) 0 ∨ H, t) ) 0

(2.21)

l ) Ca, CO, SO, Cl, OH (2.22)

To obtain the two-dimensional (horizontal) velocity field, we will use the single-phase Darcy’s law, i.e.,

( Kµ )∇P

(2.23)

where P is the pressure and µ is the viscosity. The pressure is prescribed at the production side. At the inlet, we specify the velocity and injection fluxes. The velocities at the inlet boundary are not uniform, because of different porosities and, hence, permeabilities. However, the total volume flow rate is constant at the inlet boundary. We calculate the velocity of each inlet boundary grid cell, using Darcy’s law, as follows:

(2.17) uwx,i ) Qwinj

The boundary conditions are cCa(x ) 0, y, t) ) cCO(x ) 0, y, t) ) cSO(x ) 0, y, t) ) cCl(x ) 0, y, t) ) 0 (2.18) and

(

)

(2.16)

This initial concentration of water is c°H2O ) Fw/Mw ) (1000 g/L)/(18 g/mol) ) 55.6 mol/L. The initial solid concentrations are

inj cSO(x ) 0, y, t) ) cSO , cCl(x ) 0, y, t) ) cinj Cl , inj cinj cSO Cl cOH(x ) 0, y, t) ) c°OH 1 c°SO c°Cl

inj cinj cSO Cl c°SO c°Cl

(2.20)

where, again, the y-dependence can be ignored in the onedimensional case. Clearly, the injection fluxes of HCl and H2SO4 inj inj are uwxcSO and uwxcCl , respectively. The boundary conditions of the two-dimensional model at the top and bottom boundaries read as follows:

uw ) -

) c°H2O

l ) Ca, CO

inj uwxcHO(x ) 0, y, t) ) uwx c°HO 1 -

whereas the initial concentration of water is i cOH

10083

)

(2.19)

where the y-dependence has been inserted to extend the validity to the two-dimensional case. We do not need a boundary condition for the solid components, because their transport equation does not contain a convection term. At the injection side, we prescribe the inlet flux of the relevant components. All inlet fluxes are considered constant. The gradients of all concentrations at the outlet boundary are zero. If diffusion were considered, also at the top and bottom boundary, the concentra-

( ) Ki

Ny

∑AK

i ) 1, ..., Ny

(2.24)

j j

j)1

where Qwinj is the injected-volumetric flow rate, and Aj and Kj are the cross section and permeability of each inlet boundary grid cell, respectively. Finally, Ny is the total number of the grid cells in the y (vertical)-direction. 3. Analysis of Equations In a model simulating reactive fluid flow, it is essential to have fast and accurate speciation calculations. For this reason, the species concentration should be considered as primary unknowns.21 In the numerical model derived below, the reaction rates are treated as the secondary variables. Therefore, they are dependent on the primary variables (concentrations) and are not independent. It is convenient to obtain equations for the divergence of the Darcy velocity (i.e., div uw) and the rate of change of the porosity (∂φ/∂t). Therefore, we divide eqs

10084


2.8-2.12 by the pure concentrations c°SO, etc., and then add the ensuing equations. Upon substituting eq 2.5, we obtain r″″ r″″ r″″ r″″ r″″ ∂ S C C S - r″″ C S + r″″ C +2 + φ + ∇ · uw ) ∂t c°SO c°Cl c°Ca c°OH c°CO (3.1) Furthermore, we use eq 2.1 and add eqs 2.13 and 2.14 to obtain an equation for the rate of porosity change, i.e., 1 ∂ 1 φ ) - (r″″ + r″″ r″″: ) Υ(VC,cSO,cCl) S) + ∂t FC C FS S

(3.2)

φ(t + ∆t) - Υ(t + ∆t)∆t ) φ(t)

and for the Darcy velocity, we again apply the implicit upstream finite volume approach on eq 3.3 and obtain the following equation: uwm(t + ∆t) - R(t + ∆t)∆x ) uwm-1(t + ∆t)

r″″ r″″ r″″ r″″ r″″ S C C S - r″″ C S + r″″ C +2 + + c°SO c°Cl c°Ca c°OH c°CO 1 1 (r″″ + r″″ r″″: ) R(VC,cSO,cCl) S) FC C FS S

VC(t + ∆t) -

(3.3)

As mentioned above, in the one-dimensional model, we only consider behavior inside the fracture and ignore the presence of the surrounding matrix. The 1-D results show the sequence of processes that occur along streamlines. However, the 1-D model is poor representation of the behavior in more dimensions because it hides the unstable nature of the displacement process. We include the 1-D model to understand the sequence of waves that occur along streamlines without cross talk between streamlines. All of our flow equations (eqs 2.8-2.12) have the following form:

(4.1)

The fully implicit finite-volume approach and first-order upstream weighting lead to the following equation: ∆t m m (u c )(t + ∆t) - r″″l (VC,cSO,cCl) (φcl)(t + ∆t) + ∆x w l ∆t m-1 m-1 (u c )(t + ∆t) (t + ∆t)∆t ) (φcl)(t) + ∆x w l

Gl (V m)(t + ∆t) +

where m denotes the grid cell number. The advantage of this formulation is that the Darcy velocity is dependent only on the upstream velocity and not on the downstream velocity. Therefore, it is possible to solve for all unknown dependent variables in each cell separately, without having to solve all equations simultaneously in the full domain. The disadvantage is that the numerical diffusion of these first-order schemes is relatively large. However, one can compensate for this phenomenon partially by using a larger number of grid cells. For the porosity, we can write eq 3.2 using the implicit upstream finite-volume approach as follows:

(4.6)

∆t m u F (V m)(t + ∆t) ∆x w l Ωl (V m)(t + ∆t) ) Ll m,m-1

(4.7)

where Ll m,m-1 ) Gl (V m)(t) +

∆t m-1 u F (V m-1)(t + ∆t) ∆x w l

(4.8) We have introduced the notation Ll m,m-1 for the right-hand side of eqs 4.2-4.6, and Gl, Ωl, and Fl are the auxiliary nonlinear functions. Equations 4.3-4.6 can also be written in the form of eq 4.7 by assigning zero values when appropriate. We assume that Flm-1(t + ∆t) and uwm-1(t + ∆t) have been precomputed when solving the equations for the previous cell (m - 1). We emphasize that Ll m,m-1 is not dependent on the condition of cell m at the new time (t + ∆t). 4.2. Solution of the Nonlinear System. The system (eq 4.7) is solved using the Newton-Raphson (NR) method. Given an approximate solution for the kth iteration for V k and uk of eq 4.7, the NR method finds a better approximation for the (k + 1)th iteration. Equation 4.7 becomes Gl (V k+1) +

(4.2)

1 r″″(t + ∆t)∆t ) VS(t) FS S

4.1. Generalized Notation. Let us rewrite eqs 4.2-4.6 with a concise notation V for the five unknown liquid concentrations (cSO, cCl, cCa, cCO, cOH), two solid volume fractions (νC and νS), and the porosity (φ). We obtain the nonlinear implicit scheme

4. Numerical Procedure of the 1-D Equations

∂ ∂ (φcl) + (uwcl) ) r″″l (VC,cSO,cCl) ∂t ∂x l ) Ca, CO, SO, Cl, OH

1 ( ) (r″″(t + ∆t) + r″″ S t + ∆t )∆t ) VC(t) FC C (4.5)

VS(t + ∆t) +

We can still use the five equations (eqs 2.8-2.12) for the liquid components and two other equations (eqs 2.13 and 2.14) for the solid components. We have a total of nine unknown variables, i.e., five unknown liquid concentrations (cSO, cCl, cCa, cCO, cOH), two unknown solid volume fractions (νC and νS), one unknown porosity, and finally one unknown velocity in the onedimensional model or the unknown pressure field in the twodimensional model.

(4.4)

In the one-dimensional situation, it is not necessary to use Darcy’s law, because all the results follow from the incompressibility of the fluids and the rock. Moreover, the numerical scheme for the solid components (eqs 2.13 and 2.14) can be written as follows:

Substitution of eq 3.2 into eq 3.1 leads to an expression for the divergence of the Darcy velocity: ∇ · uw )

(4.3)

∆t k+1 u F (V k+1) - Ωl (V k+1) - Ll m,m-1 ) 0 ∆x w l (4.9)

Substituting V k+1 ) V k + dV, uk+1 ) uk + du, and neglecting second-order terms, we obtain

(

)

∂Gl k ∆t k ∂Fl k ∆t (V ) + u (V ) dV + F (V k)du ) -Rl k ∂V ∆x w ∂V ∆x l (4.10)

where we have defined R lk as Rl k ) Gl(V k) +

∆t ( ∆x )u F (V ) - L k w l

m,m-1

k

l

(4.11)

This is solvable for (dV, du) if u∆t/∆x is not a characteristic speed, which can be achieved by taking ∆t small enough. After

Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010 k

division of this equation by u ∆t/∆x, we obtain the following linear system to be solved at each Newton iteration:

(

∆x M V, - k u ∆t

( )( ) -Rk1

dV1 dV2 l dV8

)

du/uk

-Rk2 ) l -Rk8

(4.12)

-Rk9

where M is the 9 × 9 Jacobian matrix. 4.3. Numerical Implementation. The quantities in the grid cells are computed in the injection-production direction from the left to the right. To facilitate the computation of the fluxes of all components at the injection boundary, we discretize the space variables in a uniform vertex-centered manner. All calculations in the NR scheme are dependent on the old condition of cell m, as well as on available information from cells to the left of cell m. The method of solution is dependent on the new condition of the cell. We stop the one-dimensional model when the porosity of one cell becomes zero. Simulations use a uniform grid with 180 blocks. This 1-D implicit method is inexpensive, because it only involves the solution of many 9 × 9 matrices, as opposed to a single big matrix. More information on this numerical procedure can be found in Bruining and Marchesin.22 5. Numerical Procedure of the 2-D Model For the two-dimensional (2-D) setting, we cannot use the 1-D numerical procedure, because the information of each grid block is dependent on the four neighboring grid blocks. Moreover, we do not know the upstream direction, because the permeability of each grid block changes with time nonlinearly. Automatic implementation of the upstream approximation of 3(uc) uses the following formulation:23 1 (uxcl)(xi+1/2,j) ) [(ux + |ux |)i+1/2,jci,j + 2 l ) Ca, CO, SO, Cl, OH (ux - |ux |)i+1/2,jci+1,j]

(5.1)

To calculate the Darcy velocity at the boundary between each pair of grid blocks, we use the geometric average permeability, i.e., ui+1/2 )

√KiKi+1 µ

(

pi - pi+1 xi+1 - xi

)

(5.2)

We use nonuniform grid blocks in the y-direction and uniform grid blocks in the x-direction. Two-dimensional simulations use 30 grid blocks in the y-direction (Ny ) 30), i.e., 10 grid blocks in the fracture zone and 20 grid blocks in the matrix part, and 40 grid blocks in the x-direction. Therefore, we use 1200 grid blocks for the 2-D simulation. The system of equations again consists of nine partial differential equations (eqs 2.8-2.14, 3.2, and 3.3). For discretization in space, we use the implicit upwind finite-volume method on the cell-centered scheme. We use the NR approximation to linearize the nine nonlinear equations. The total number of unknown variables is equal to 9 × Nx × Ny, i.e., 10 800 in our case. Hence, the size of the Jacobian matrix in the NR method is 10 800 × 10 800, which is big, compared to the size of the 1-D Jacobian matrix, i.e., 9 × 9. However, current software simplifies the 10 800 × 10 800 matrix by discarding all the zeros. After obtaining the permeability field of the domain, we use the effective-medium approximation24 to calculate the fracture and matrix permeability, using

∫

∞

0

h(K)

K - Keff K + (γ-1 - 1)Keff

dK ) 0

10085

(5.3)

where Keff is the average permeability, h(K) the distribution function (h(K) ) 1/(Nx × Ny)), and γ a geometry factor (which, for our purposes, we assume to be 3/10). In three dimensions, the percolation threshold is between 0.2 and 0.3, depending on the coordination number. Here, we have used a value of 0.3. This number is approximately right in three dimensions. (In two dimensions, the percolation threshold is exactly right in the effective-medium approximation.) Further details can be found in Sahimi.24 6. Results and Discussion This section is organized as follows: First, we discuss the one-dimensional (1-D) results. Subsequently, we discuss the two-dimensional (2-D) results. Then, we describe one laboratory experiment and compare it to the two-dimensional simulation case. For all simulations, the general input data are shown in Table 1, and the specific information for all experimental and simulation cases are summarized in Table 2. 6.1. One-Dimensional Results. The 1-D results show the sequence of processes that occur along streamlines. Our model disregards diffusion and, therefore, the one-dimensional results can be used in a simple streamline simulator without cross talk between streamlines. They are vital to interpret the 2-D results. The main parameters are (1) the flow rate, and (2) the acid ratio (HCl/H2SO4). We illustrate the effect of these parameters with three cases. The initial value of the porosity for all three cases is φ ) 0.5. For incompressible flow, the permeability, heterogeneous or not, does not enter in the equation. For all three cases, we obtain the injection velocity by dividing the injection flow rate by the fracture cross section, which is the diameter of the core multiplied by the width of the fracture. For the interpretation, it is useful to define a dimensionless number that characterizes the reaction rates in porous media. Here, we use the Damko¨hler number (Da). The Da number is the ratio between the acid reaction rate and the convection mass transport rate. It also determines the distance that the acids can penetrate into the rock before it becomes consumed. However, it is not possible to define a meaningful Da number for the entire process in time, because the volume fraction of CaCO3 (νC) in the reaction rate varies from its initial value to zero. Therefore, we calculate the initial Da value for the cases considered below. For our situation, we use the following expression for the Damko¨hler number: Dal )

kl aνCL u

l ) S, C

(6.1)

where the subscripts “S” and “C” denote H2SO4 and HCl, respectively. The first case has a Darcy velocity of 0.042 cm/s, an acid ratio of HCl/H2SO4 ) 0.3/0.9 mol/mol, and, consequently, DaS Table 2. Specific Information of All Experimental and Simulation Cases Injected Concentration of Acid (mol/L) length diameter fracture injection rate case (cm) (cm) width (cm) (mL/min) 1 2 3

19.0 19.0 19.0

9.90 9.90 9.90

0.2 0.1 0.2

5 5 5

H2SO4

HCl

0.9 0.9 1.5

0.3 0.3 0.3

10086


Figure 2. Case 1: (a) outlet concentration of the components, (b) volume fractions of the solid components at the time of clogging, and (c) concentration of the components at the time of clogging.

) 1.27 and DaC ) 126.67. Figure 2a shows the produced concentrations of the aqueous components versus time. As shown in Figure 2a, after 3 min, breakthrough of the CO2, CaCl2, and H2SO4 occurs. However, HCl does not reach the production point; this is attributed to the higher reaction rate constant of HCl, with respect to that of H2SO4. Figure 2a shows that the maximum concentration of CO2 occurs at the time of breakthrough and, after that, this concentration gradually decreases. This is caused by the fact that increasingly less H2SO4 is reacting in the core, because some of the injected amount of H2SO4 is produced unchanged. In addition, also minor amounts of CaCO3 are consumed and can no longer participate in the reactions. Note that CO2 is the product of both the dissolution (reaction with HCl) and precipitation reactions (reaction with H2SO4). After breakthrough of H2SO4, not all of the injected H2SO4 is used to form gypsum; therefore, less CO2 forms. In contrast to CO2, CaCl2 is constant after breakthrough, because no breakthrough of HCl occurs, which means all the injected HCl concentration reacts with CaCO3. Note that the injection flow rate and injected HCl concentration are constant. The maximum porosity in this case (i.e., φ ) 0.9) occurs at the inlet. From there, it monotonically decreases to zero porosity at a distance


9.98 cm from the inlet after 2.33 h of injection. This is shown in Figure 2b. In addition, the volume fraction of CaCO3 and CaSO4 · H2O are also shown in this figure. These volume fractions, added to the porosity, add up to one at each location. All the concentrations are equal to the injection concentration upstream of the point where the CaCO3 concentration becomes zero. Because the porosity is a monotonically decreasing function in the flow direction, the point at which porosity equals the initial porosity remains at a fixed position. The concentration profiles within the core at the time of clogging are shown in Figure 2c. The HCl concentration only penetrates half of the core, but the H2SO4 concentration diminishes only slightly. The carbon dioxide (CO2) concentration and CaCl2 start to differ from zero at the boundary at which all CaCO3 is consumed by HCl and H2SO4. The second case has a Darcy velocity of 0.084 cm/s, which is the same acid ratio as the first case, and DaS ) 0.63 and DaC ) 63.33. Figure 3a shows the production concentration of the aqueous components versus time. In this case, breakthrough of CO2, CaCl2, and H2SO4 occurs after 1.7 min of injection, which


is about half of the breakthrough time of the first case. This is to be expected, because the injection velocity of the second case is about twice the injection velocity of the first case. The fact that it is not exactly inversely proportional to the velocity is due to variation of the velocity inside the core. The HCl reaches the production point after 2 h of injection. Figure 3a indicates that, after breakthrough of HCl, the CaCl2 concentration starts to decrease. We explain this observation by noting that, after CaCl2 breakthrough, the total chloride (of CaCl2 and HCl) concentration becomes constant and equal to the injected HCl concentration if we ignore the small variations of the Darcy velocity. The maximum CO2 concentration in the second case is smaller than the maximum CO2 concentration in the first case. However, the maximum concentrations of CaCl2 for both cases are the same. This means that the residence time affects the H2SO4 reaction but does not influence the HCl reaction. Hence, we conclude that the characteristic reaction time of H2SO4 is of the same order of magnitude as the residence time. However, the characteristic reaction time of HCl is much smaller than the residence time. The maximum porosity in this case (i.e., φ ) 0.92) occurs at the inlet (see Figure 3b), which is larger than that of the first case. From there, it monotonically decreases to almost zero porosity at the production side after 2.33 h of injection. Figure 3b also shows the volume fraction of CaCO3 and CaSO4 · 2H2O. Figure 3c shows the concentration profiles within the core 2.33 h after the start of the injection. The HCl concentration penetrates the entire core. The H2SO4 concentration is almost constant. Again, as in the first case, the CO2 concentration and CaCl2 start to differ from zero at the boundary at which all CaCO3 is consumed by HCl and H2SO4. The third case has the same Darcy velocity and, consequently, the same Da value as the first case, but with a different acid ratio, i.e., HCl/H2SO4 ) 0.3/1.5 mol/mol. Figure 4a shows the production concentration of the aqueous components versus time. In this case, breakthrough time (3 min) is almost the same as that in the first case, because the injection velocity is the same for both cases. Again, as in the first case, HCl does not reach the production point. The maximum CO2 concentration is larger than that of the first case, because of the higher injected H2SO4 concentration. The maximum porosity in this case (i.e., φ ) 0.85) occurs at the inlet (see Figure 3b), which is smaller than that of the first and second case. From there, it monotonically decreases to zero porosity at a distance of 6.16 cm from the injection point after 1.27 h of injection. Note that this clogging time is smaller than that of the first case. Figure 4c shows the concentration profiles within the core at the time of clogging. The HCl concentration penetrates one-third of the core. Nevertheless, the H2SO4 concentration is present over the entire core. Again, the CO2 concentration and CaCl2 concentration start to differ from zero at the boundary at which all CaCO3 is consumed by HCl and H2SO4 (see Figure 4c and 4b). We investigate the flow-rate effect by comparing the firstcase and second-case results and using the Da number. The maximum porosities in the first and second case are 0.9 and 0.92, respectively. In the first case, the clogging occurs at a distance of 9.98 cm; however, in the second case, complete clogging does not occur. Now, the porosity reaches a finite (minimum) value at the production point. This is also reflected by the fact that the initial Da numbers for the first case (DaS ) 1.27, DaC ) 126.67) are larger than for the second case (DaS ) 0.63, DaC ) 63.33). This means that the residence time of the acids in the first case is twice as long as that in the second case. Consequently, the acid-reaction zone is shorter in the first case. In both cases, the simulations show that the dissolution is

10087


the dominant process upstream, whereas precipitation occurs downstream. We define the dissolution-dominant length as the length at which the porosity is bigger than the initial porosity (0.5). It starts from the injection point and ends at the point at which the summation of the porosity and gypsum volume fraction becomes one, i.e., the intersection point of the porosity and the CaSO4 · 2H2O volume fraction plot in Figures 2b, 3b, and 4b. From the injection point to this point, the CaCO3 concentration is zero. After this, the solution upstream of this point remains invariant. Based on this definition, we obtain dissolution-dominant lengths of 1.9 and 3.7 cm for the first and second case, respectively. Therefore, a higher flow rate leads to more dissolution-dominant reactions upstream. Moreover, the clogging point would be further away from the injection point, because a smaller Da number would lead to a longer reaction zone. If the length of the domain were shorter than the distance between the injection point and clogging point, clogging cannot happen within the core. This situation is not desirable for our purpose, which is reducing the fracture permeability.

10088


Figure 5. Porosity surface contours (horizontal cross-section) for Case 1: (a) at t ) 15 min, (b) at t ) 1 h, (c) at t ) 3 h, and (d) at t ) 6 h. The instability develops in panels (c) and (d) on one side of the fracture; the initial perturbation is due to the random nature of the initial matrix porosity (and permeability) field.

We compare the first-case and third-case results to illustrate the acid-ratio effect. The maximum porosity in the third case is smaller than that in the first case. In the third case, the clogging occurs at a distance of 6.16 cm, whereas in the first case, the clogging happens at a distance of 9.98 cm from the injection point. Again, the third case results show that dissolution is the dominant process upstream. The dissolution-dominant length of this case is (1.1 cm) smaller than that of the first case (1.9 cm). Hence, a lower acid ratio leads to less upstream dissolution and more downstream precipitation. Furthermore, the clogging point would be closer to the injection point. Whether a case with a lower acid ratio is desirable is discussed here below in the 2-D results. 6.2. Two-Dimensional Results. Figure 1 shows the geometric outlay for the 2-D simulations. Initially, the rectangular domain consists of two subdomains, viz, one domain that carries the fracture and one domain for the matrix block. The fracture has a width of a few millimeters and runs over the length of the entire core (see Table 2). Note that, in Figures 5, 7, and 9 the number of grid cells in the y-direction are indicated, seemingly exaggerating the width of the fractures. Figures 5, 7, and 9 show the porosity development and pertain to three different cases, whereas Figures 6, 8, and 10 use the Carman-Kozeny relation and subsequently use the effectivemedium approximation (eq 5.3) to find the average fracture and matrix permeability for all three cases. Initially, the fracture permeability is 300 times as large as the matrix permeability. The displayed porosity-field history in Figure 5 shows that the porosity increases upstream, but it also shows that the porosity decreases downstream. This continues until the fracture becomes clogged and an alternative high-porosity pathway is created along the axes of the fracture. These highly permeable pathways are also called wormholes. Figure 5 (case 1) differs from Figure

Figure 6. Relative effective permeability evolution for Case 1.

7 (case 2), in that the dissolution-dominated length is larger for case 2. In Figure 7, the high-porosity flow path is larger than in Figure 5. The low-porosity region in the fracture shows much less variation in Figure 5 with respect to the region in Figure 7. Figure 9 (case 3) shows the effect of the acid ratio. The upstream dissolution-dominant zone is shorter than that in the other two cases. In the third case, the high-porosity region starts at the “bottom”, but it is not as porous as that in the other two cases. In addition, the instability develops in Figures 5, 7, and 9 on one side of the fracture; the initial perturbation is due to the random nature of the initial matrix porosity (and permeability) field. The porosity of the surrounding matrix zone is more affected than that in the other cases. Figure 6 (case 1) clearly indicates that the minimum values of the fracture and matrix relative (Keff /Ki) average permeabilities are 0.04 and 0.79 after 6 h, respectively. On the other hand, Figure 8 (case 2) shows that the minimum values of the fracture and matrix relative


10089

Figure 7. Porosity surface contours (horizontal cross-section) for Case 2: (a) at t ) 15 min, (b) at t ) 1 h, (c) at t ) 3 h, and (d) at t ) 6 h. The instability develops again on one side of the fracture, because of the random nature of the initial matrix porosity (and permeability) field.

Figure 8. Relative effective permeability evolution for Case 2.

average permeability are 0.15 and 0.8, respectively, after 6 h. This means that the matrix relative permeability is not very sensitive to the fracture width. Moreover, it indicates that the acid treatment also reduces the matrix permeability, albeit less than the fracture permeability. However, acid treatment in the second case with half of the fracture width of the first case shows less reduction of the fracture permeability. Therefore, higher flow velocity inside the fracture leads to a more dissolutiondominated zone and a more-pronounced high-flow path. Figure 10 (case 3) shows that the minimum values of the fracture and matrix relative average permeabilities are 0.006 and 0.7, respectively, after 6 h. Comparing the results of the third case with the results of the first and second case reveals that a higher acid ratio HCl/H2SO4 leads to more precipitation for both the fracture part and the matrix part. Reducing the fracture permeability is desirable, but reducing the matrix permeability is not

favorable. Hence, a higher acid ratio has one advantage and one disadvantage, i.e., (1) it blocks the fracture, which avoids short-circuit flows between the injection and production well, but (2) it makes diversion of water flow to the matrix more difficult. 6.3. Some Preliminary Experimental Results. Parallel to the theoretical research, we perform an experimental program with acid injection in chalks under various pressure/temperature (P,T) conditions. Hence, it is possible to compare preliminary experimental results with theory. In the experiments, we used a cylindrical core with a length of ∼20 cm, a diameter of ∼10 cm, and a fracture width of 1 mm. The experiments were run at ∼80 °C, an annular pressure of ∼10 MPa, a pore pressure of ∼5 MPa, and an injection rate of 5 mL/min. The experimental results discussed here are associated with Case 2 (see Table 2). At the start, the core contains brine (with a total NaCl amount of 0.21 mol) instead of fresh water. Thereafter, a mixture of HCl and H2SO4 was injected, resulting into the following generic features. The experiment discussed here (Experiment 1) has high flow rates. The initial pore volume of the core divided by the injection rate is ∼1.5 h, which is used as a measure of the initial residence time. In total, 1.36 mol of H2SO4 and 0.45 mol of HCl were injected. Figure 11 shows the produced ionic concentrations of the components versus time obtained from a high-performance liquid chromatography (HPLC) measurement. Based on the HPLC results (Figure 11) during the experiment, 0.0071 mol of CaCl2 has been produced, at a constant rate for 5 h. A total of 0.37 mol of gypsum (CaSO4 · 2H2O) has been formed, which mainly fills the original fracture, but also the matrix. A total of 0.56 mol of Cl- was produced. Following the mass balance

10090


Figure 9. Porosity surface contours (horizontal cross-section) for Case 3: (a) at t ) 15 min, (b) at t ) 1 h, (c) at t ) 3 h, and (d) at t ) 6 h. The instability now develops on the other side of the fracture, with respect to Figures 7 and 8. We used another realization for the initial matrix porosity (and permeability) field.

Figure 10. Relative permeability evolution for Case 3.

considerations, an estimated 0.1 mol of CaCl2 remains in the core. Figure 12a shows a computer tomography (CT)-scan image of the core before acid injection. A CT-scan image of the core after the experiment (Figure 12b) shows a single remaining wormhole extending from the injection to the production side. The formation surrounding the wormhole has been mostly converted to gypsum. All of this is consistent with the following conceptual model. Most of the fluid is produced through the wormhole. Because the wormhole is surrounded by gypsum, the HCl passes the core without much reaction to CaCl2. Occasionally, new branches of the wormhole are formed and the old branch is clogged with gypsum; however, this is not so evident from the experiment discussed here. Figure 13 shows a perpendicular slice through the matrix after acidizing with cemented-wormhole

Figure 11. HPLC results during acid injection, showing a decrease of Cland Na+ and increase of Ca2+ and SO2-. After ∼1.5 h, acid production starts.

areas. Therefore, the small amount of CaCl2 produced must be due to some reaction of HCl with CaCO3. To give a general idea, we also make some remarks about experiments performed under different conditions. An elaborate discussion on the experimental results can be found in the work by De Iongh.25 Comparing Experiment 1 to an experiment where the acid injection rate has been reduced by a factor of 2.5, while keeping the fracture width and acid ratio constant, reveals that the fracture has almost been fully blocked by precipitation of gypsum. However, since the SO42- concentration remained low in the production fluid, we conclude that dissolution of carbonate and precipitation of gypsum continued until the end of the experiment, as opposed to the early acid breakthrough that was


10091

Figure 12. Three-dimensional (3-D) rendering of CT scans, visualizing the open holes and wormholes from the injection to the production side: (a) CT-scan image of the initial condition (i.e., before acid injection), where the round spots are the pads to keep the fracture open; and (b) Case 2 after acid injection. Pores

Modeling the Closure of Thief Zones in Carbonate Reservoirs Using

Recommend Documents