Journal of Hydrology 450–451 (2012) 134–139 Contents lists available at SciVerse ScienceDirect Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol Porosity and permeability changes in sedimentary rocks induced by injection of reactive fluid: A simulation model Supti Sadhukhan a, Philippe Gouze b, Tapati Dutta c,⇑ a Physics Department, Jogesh Chandra Choudhury College, India Geosciences, Universite de Montpellier 2, CNRS, Montpellier, France c Physics Department, St. Xavier’s College, Kolkata 700 016, India b a r t i c l e i n f o Article history: Received 26 November 2011 Received in revised form 26 February 2012 Accepted 9 May 2012 Available online 23 May 2012 This manuscript was handled by Laurent Charlet, Editor-in-Chief, with the assistance of Eric C. Gaucher, Associate Editor Keywords: Reactive transport Porosity Permeability s u m m a r y Numerical programs for simulating flow and reactive transport in porous media is essential for predicting reservoir properties related to CO2 sequestration performance, subsurface storage and risk assessment. In this paper we solve the Navier Stokes’ equation using finite difference method, on a simulated porous rock structure, to study the velocity distribution of fluid flowing through it under a constant pressure gradient. A reactive solute carried through the fluid is allowed to interact with the minerals in the rock. This chemical reaction dissolves the mineral which changes the rock structure thus affecting its flow properties. These changes of flow properties are studied with variation in reactive solute concentration and pressure field. The different mechanisms of dissolution responsible for the variation of flow properties for the different parameters is predicted. Before the onset of homogeneous dissolution, variation in porosity follows a power-law behaviour with change in permeability when the latter is scaled by the concentration of the reactive species. The simulation results are compared with available experimental data and found to give a reasonable match. Ó 2012 Elsevier B.V. All rights reserved. 1. Introduction Transport of reactive fluids through porous rocks can lead to both dissolution of some of the minerals, as also precipitation of newer ones. These processes can create new conducting channels or block existing paths of flow. Both these processes will affect subsurface storage, CO2 sequestration, underground hydrocarbon production and study of risk assessment in underground flow. Accounting for both physical and chemical heterogeneity of the rocks at different scales is a challenge to be addressed by any simulation model in such studies. The definition of a pertinent support volume in the frame of the continuum approach has been studied by (Whitaker, 1999; Hornung, 1997). Li et al. (2006) and Li et al. (2007) have used the pore-scale model to investigate the spatial distribution of reaction rates, while (Kang et al., 2002, 2003, 2005, 2010) have used Lattice Boltzmann model to study reactive transport in porous media. Several studies have considered the up scaling of mass transfer processes in heterogeneous media (Edwards et al., 1993; Litchner and Tartakovsky, 2003; Meile and Tuncay, 2006) using different analytical methods. However physical and chemical mechanisms have often been considered separately. Experimental data on the effects of reactive fluid flow ⇑ Corresponding author. Tel.: +91 9330802208; fax: + 91 033 2287 9966. E-mail address: tapati_mithu@yahoo.com (T. Dutta). 0022-1694/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jhydrol.2012.05.024 through porous rocks are few, Luquot and Gouze (2009) have studied porosity and permeability changes by injection of CO2 through limestone reservoir samples. Real rocks are 3-D structures with highly tortuous and often fractal pore spaces. Before attempting simulation of this daunting geometry, we present here, as a preliminary study a simpler 2-D version. In the following sections, we first briefly discuss our basic algorithm for the generation of the porous rock and determination of velocity field of the fluid flow. Next the algorithm for mimicking reactive fluid flow following the rules of advective–diffusive mass transport is discussed. Changes in porosity and permeability due to dissolution is affected by many parameters such as flow rate, reaction rate constants, Peclet and Damkholer numbers. In this study we focus only on the effect of concentration of the reactive species on dissolution and its effect on the physical and flow properties of the sample. In the next section our results are presented and compared with available experimental data. Finally a discussion of the results follow. 2. Theory Fluid transport in a porous structure under a suitable pressure gradient is described by the Navier Stokes’ equation q dV þ ðV  rÞV þ rP  lr2 V ¼ fe dt ð1Þ S. Sadhukhan et al. / Journal of Hydrology 450–451 (2012) 134–139 where V, P, and fe represent the velocity, pressure and external force per unit volume respectively, q and l are respectively the density and dynamic viscosity of the fluid. Neglecting the inertial term and assuming no external forces acting on the fluid, Eq. (1) simplifies to dV 1 ¼  rP þ gr2 V q dt ð2Þ where g ¼ lq is the kinematic viscosity. For an incompressible fluid, the equation of continuity is rV ¼0 ð3Þ Eqs. (1) and (3) when solved together, give the steady state condition of flow in the structure. If the fluid carries a reactive solute that can interact with the mineral content of the rock matrix, chemical interaction can take place resulting in either dissolution, precipitation or both processes along the interface. This is the principle on which the possibility of CO2 sequestration in sedimentary rocks is being investigated. At the CO2 injection zone, the aquifer fluid is saturated with CO2 under an appropriate high pressure. The rock may contain minerals like magnesium calcite. The dissolution reaction is then described by CaCO3 þ Hþ ! Ca2þ þ HCO3 ð4Þ Precipitation reactions occurs with the formation of magnesium– calcite crystals. In situ experiments have been carried out and discussed in detail in Luquot and Gouze (2009). In this paper we shall be considering only the dissolution process. The dissolution reaction is described in the simplest form by: AþB C ð5Þ where A is the reactive solute carried by the fluid, B the mineral and C the product. After the reaction, C is released in the fluid and is carried along it with appropriate velocity. A pore is left behind at the site of B. If we assume that one reactive solute, i.e. A is responsible for dissolution at one site of B, then the dissolution rate can be described in terms of the rate of change in the concentration of A, i.e., dcA ¼ RcA dt ð6Þ where cA is the concentration of A at any instant of time t, and R is the reaction rate constant. Solving Eq. (6) we get the dissolution probability Pdis ¼ 1  expRtres 135 Tarafdar (2003) and Sadhukhan et al. (2007). The generated sample has a connected rock structure and above the percolation threshold, a connected pore structure too. We refer to this model as the Relaxed Bidisperse Ballistic Deposition Model (RBBDM). The basic algorithm is to deposit particles of two different sizes ballistically. In 2-d (1 + 1 model), we drop square 1  1 and elongated 2  1 grains on a linear substrate. The square grains are chosen with a probability p and elongated grains with probability (1  p). The presence of the longer grains leads to gaps in the structure. The porosity /, defined as the vacant fraction of the total volume (area in 2-d) depends on the value of p. The longer grains topple over when they form unstable overhangs thereby leading to compaction. It is well known that natural sand grains are angular and elongated (Pettijohn, 1984), so the aspect ratio 2 is realistic. It has been shown that in BBDM (Dutta and Tarafdar, 2003) the sample attains a constant porosity only after a sufficient number of grains (depending on sample size) have been deposited to overcome substrate effects. Here, a 128  1000 size sample was generated. Fig. 1 shows the 2-d section of the sample generated with p = 0.5. From this sample, a 32  32 sample was selected after the porosity had stabilized. The 32  32 sample was magnified to 128  128, and simulation was carried out on it. This magnification is necessary to make all channels wide enough, so that a realistic velocity profile with no-slip at the walls can be implemented (Sadhukhan et al., 2007). Reaction on the pore-rock interface is then observed at a finer scale than the scale of pore geometry. In the simulation, the larger grain is identified as mineral D while the smaller grain of the rock composition is identified as mineral B, which only can dissolve to release cation C in the fluid. Mineral D has no role in the dissolution reaction process. Though the composition of the rock can be changed by varying the proportion of D–B, in this paper we report our simulation results with rock sample containing only 50% of B. The porous rock is generated stochastically. All results reported are an average over 150 such configurations. The Hoshen and Kopelman algorithm (Hoshen and Kopelman, 1976) is used to identify the channels spanning the sample. Fig. 2a shows a spanning cluster before reactive transport can begin through it. The pressure and velocity fields are solved by the procedure described by (Sarkar et al., 2004) with some necessary departures appropriate to our problem. The space and time discretized versions of Eqs. (1) and (3) are used iteratively to obtain the steady state flow (Sadhukhan et al., 2007). This was identified ð7Þ tres is the time of residence of the reactive solute in a given grid cell. Also the reaction rate is given by: R ¼ krx ð8Þ where k is the kinetic coefficient, r is the specific surface area, x is the relative concentration of the reactive solute with respect to its equilibrium value. The kinetic coefficient k is given for a particular reaction, and x ¼ ccA0A ; cA0 being the concentration of A at equilibrium. As dissolution occurs, the pore space and geometry keeps changing. This alters the pressure distribution and the velocity field which is updated according to Eq. (1). 3. Simulation We generate a porous stochastic structure in 2-d and simulate flow of a single fluid through it, using a numerical finite difference solution of the steady state Navier Stokes’ equation. The details of the generation of the porous rock sample are given in Dutta and Fig. 1. A 2-d section of the simulated structure for p = 0.5. The elongated pore clusters are along the growth direction. 136 S. Sadhukhan et al. / Journal of Hydrology 450–451 (2012) 134–139 Fig. 2. The different colour shades correspond to different values of velocity with the velocity increasing from darker to lighter shade. (a) Shows a typical spanning cluster before reactive fluid flow begins. The isolated pore clusters and rock matrix have been deliberately masked. (b) Shows the same spanning cluster after 100 time steps of fluid flow for a walker concentration of 1000. The general broadening of the channel is due to homogeneous dissolution that becomes prominent at later times. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) when the difference of velocity between successive time steps of iteration was 109 or less. The reactive solute carried by the fluid is represented by a random walker. The initial concentration of the active solute is determined by the fixed number of walkers that are released simultaneously from one end of the percolating flow channel. Every walker moves from one grid to another according to the local velocity of the grid as determined by the velocity field. A walker moves in equal ’jump steps’, the step size being equal to a grid size. The local velocity determines the time of travel s, of every walker from grid to grid. We assume that a walker divides its time of travel equally between the two grid cells. Therefore the time of residence in any grid cell is t res ¼ s 2 ð9Þ .A walker is capable of either advection or diffusion, which it decides stochastically. The velocity component having higher magnitude at every cell, determines the direction of advection in that cell. Diffusion however can occur in all directions with equal probability. The probability of advection and diffusion is given by Pdiff 1 ¼ 4ðPe þ 1Þ Padv ¼ 1  Pdiff ð10Þ ð11Þ where Pe ¼ UL , U is the local velocity, L is the grid length and D being D the diffusivity. The time of residence of every walker is calculated according to Eq. (9) where s is the inverse of local velocity at the grid. Should a walker undergo diffusion, s is taken to be the inverse of the minimum velocity value as diffusion time is the longest time of travel. cA at any instant is calculated as the number of active walkers per unit volume of the fluid. cA0 of Eq. (8) is a parameters taken to be of unit value. With the help of Eqs. (8) and (7), a walker A dissolves a site of B creating a new pore while releasing a cation C in the fluid. The random walker ’dies’ after it completes dissolution. As there is no further replenishment of the active species in the fluid, as the walkers proceed along the flow channel, their concentration weakens with time. In experimental situations, there is a constant flux of the reactive fluid from one end of the channel. As the reactive fluid percolates, chemical reaction with the minerals present in the rock, occurs. The concentration of the reactive species in the fluid, decreases with distance from the inlet. Our simulation here is a study on a one time injection of reactive fluid. The time development of the changes that the fluid causes as it flows through the sample is tracked. In the simulation, all the walkers that were originally released, are tracked parallely. Dissolution at all possible sites occur simultaneously at any time step. The cations that are released continue to move through the channel in the same way except that they are now not reactive. We have done the simulation assuming that every walker gives rise to one cation upon causing dissolution. After dissolution, we begin a new ’time step’ by updating the pressure and velocity at every point according to Eq. (1). Now the remaining active walkers as also the cations move following the same rules as before. The permeability is calculated at every time-step according to Darcy’s law q¼ j rP l ð12Þ where q is the flux, j is the permeability, l is the viscosity, and rP is the pressure gradient across the sample. Flow modulated properties of the sample such as porosity, permeability, specific surface area are tracked over 150 time-steps. Fig. 2b shows the structure that has evolved from an initial structure as shown in Fig. 2a after 100 time steps. The results are discussed in the following section. 4. Results The 2-dimensional rock structure was generated with p = 0.5. The average value of initial porosity / obtained as the volume fraction of the number of voids in the sample was calculated over 150 configurations. Before the transport of the reactive fluid began, its value was 0.42. This falls in the range of porosity values reported for sedimentary rocks. For a given pressure gradient and value of kinematic viscosity, permeability can be calculated easily using Eq. (12). The reactive fluid is allowed to flow through a constant pressure gradient of 1000 units across the sample. The different number of walkers mimic the concentration of the reactive species in the fluid. The initial average permeability j prior to transport, was 0.00543 cm2. Fig. 3 shows a variation of porosity with time for different walker numbers. Except for very low partial pressure of the reactive species, here represented by 100 walkers, the other curves show definite changes in slope at two different time values which are approximately the same in each case. These are demarcated as t1 and t2 in the graph. t1 is the limit of inhomogeneous dissolution which increases the porosity without causing much change in the permeability as the new pores may end up as dead ends. This is also the region where worm-holing, another important manifestation of reactive transport, is likely to begin. Between t1 and t2 there is a region where porosity starts increasing and beyond t2, S. Sadhukhan et al. / Journal of Hydrology 450–451 (2012) 134–139 137 Fig. 3. Variation of porosity with time for different walker number. Inset shows the experimental results of Luquot and Gouze (2009). this increase becomes rapid. In the case of moderate partial pressure (1000 walkers) of the reactive species, the increase in porosity is less than in the cases of 5000 and 10,000 walkers as obviously the concentration of the reactants decrease as they travel down the sample. This decreases the dissolution probability according to Eqs. (7) and (8). However in the cases of high partial pressure of the reactants, the concentration of the reactants is sufficiently strong to continue dissolution which is primarily homogeneous in nature. This rapidly increases the porosity. The nature of the graphs is qualitatively similar to the experimental observations of (Luquot and Gouze, 2009). Their results are shown as an inset in Fig. 3 where they had noted changes in the porosity with time by passing CO2 through carbonate rocks for different partial pressures of the CO2. The porosity changes rapidly for increased partial pressure of the dissolved gas, shown as D1. The variation of porosity is very little, graph D3 as shown in the inset. This matches the nature of the graphs obtained from simulation. In Fig. 4 the permeability change with time is shown for the different walkers. For low partial pressure of the reactants, the permeability does not change much beyond t1 as the concentration becomes too low for dissolution. For moderate to high partial pressure of the reactants, up to t1, change in the permeability is low as most of the dissolution is inhomogeneous in nature creating dead Fig. 4. Variation of permeability with time for different walker numbers. ends to the channel. Between t1 and t2, homogeneous dissolution sets in gradually increasing the permeability of the sample. Beyond t2, homogeneous dissolution occurs with further smoothing of channel walls. This leads to an average increase in channel width which dramatically increases the permeability. This increase of average channel width is evident in a comparison of Figs. 2a and b. Fig. 5 shows the variation of relative permeability change with time. The relative permeability change is rapid from beyond t2 for different concentrations. The inset shows the variation of permeability change as obtained experimentally by Luquot and Gouze (2009). The qualitative match for the different concentrations of reactants is evident. To investigate the relationship between porosity and permeability, we plot porosity versus permeability for different walker numbers, as shown in Fig. 6. Here too, each graph corresponding to different walker number, shows a linear behaviour up to the value of porosity corresponding to the time t2, i.e. the point before the start of homogeneous dissolution. Beyond t2, permeability shows an abrupt increase with porosity. Clearly this can be attributed only to a different dissolution mechanism. This becomes more clear as we investigate the variation of interface area with time as discussed later. On a log–log scale, shown as an inset in the figure, the graphs for different walkers collapse on almost a single straight line when the permeability is scaled by concentration, i.e., the walker number. Since the change in permeability with porosity beyond t2 is within a range of 0.00015, on a log–log scale this is barely discernible. Thus we may say that the function jC shows a power law scaling with the porosity /, the exponent being 7.6. The dissolution probability Eqs. (7) and (8) is also determined by the interface area exposed to the reactants. We have studied the variation of interface surface area with time, shown in Fig. 7, to try and understand the variation of porosity and permeability due to dissolution. In the case of very low partial pressure of the reactants, the surface area of the interface keeps increasing with time. This is the region where worm-holing is more predominant. Since concentration of walkers is low, once any site is dissolved, the concentration drops in the immediate vicinity preventing further dissolution at the site. The dominant mechanism is shown in Fig. 8, region A, where the change of porosity with specific surface area is noted. This increases the specific surface area and results in an increase in the porosity. However the permeability does not change much as this only increases dead ends. With moderate concentration of walkers, the specific surface area shows very 138 S. Sadhukhan et al. / Journal of Hydrology 450–451 (2012) 134–139 Fig. 5. Variation of relative permeability change with time. Inset shows the same study by Luquot and Gouze (2009). 0.0001 695 y=0.0697*x(7.6) 0.0002 porosity (cm2) permb. t1 690 (cm2) surf. area 0.00025 permb./W(-0.031) 700 t2 685 walkers=100 walkers=1K walkers=5K walkers=10K 680 0.00015 675 0.0001 walkers=0.1K walkers=1K walkers=5K walkers=10K 670 0 0.42 0.44 0.46 40 60 Fig. 7. Variation of interface surface area with time. 0.48 porosity walkers=100 walkers=1K walkers=5K walkers=10K 0.48 Fig. 6. Variation of permeability with porosity for different walkers. Inset shows the power-law behaviour on a log–log scale. The solid line shows the fit. C porosity 0.46 slight change and the dominant mechanism is shown in Fig. 8, region B for 1000 walkers. Once again though the porosity increases rapidly, the increase in permeability is only due to smoothing of rough corners in the channel space. When the concentration of walkers is high, 5000 and 10,000 walkers, homogeneous dissolution is the most prominent process. The mechanism is illustrated in Fig. 8, region C. This increases the permeability rapidly as the channel throats are broadened significantly. The specific surface area decreases as the proportion of number of new interfaces is less than before. The relation between specific surface area and permeability for different walker numbers is shown in Fig. 9. The permeability remains almost constant in the beginning where worm holing is more dominant. At the onset of homogeneous 80 t (sec) 5e-05 0 20 0.44 B 0.42 17 A 17.5 18 18.5 19 specific surface area (1/cm) Fig. 8. Variation of porosity with interface surface area. The different mechanisms of dissolution is marked as discussed in text. S. Sadhukhan et al. / Journal of Hydrology 450–451 (2012) 134–139 700 walker=100 walker=1K walker=5k walker=10k A 695 B surf. area 690 (cm2) 685 C 680 675 670 0 5e-05 0.0001 0.00015 0.0002 0.00025 permb.(cm2) Fig. 9. Variation of permeability with interface surface area. dissolution, permeability starts increasing and becomes more rapid as large areas of the channel wall start dissolving. 5. Conclusion The concentration of reactive species in the fluid modulates the permeability of porous rocks – this is expected. However depending on this concentration, the mode of rock dissolution changes. As the concentration increases, the net porosity of the sample increases but permeability does not change in the same proportion. For very low concentration, dissolution occurs at the cost of more dead ends to the existing channel, which hardly causes any change in the permeability. This kind of dissolution results in worm holes. The specific surface area shows an expected increase. With moderate values of concentration, more homogeneous dissolution sets in. The rough edges of the rock in the channels are smoothed out. The specific surface area starts decreasing. This leads to an increase in the permeability as the pore throats increase. With high concentration of the species in the fluid, there is wearing away of large chunks of the layers of the channels. This increases the channel throats considerably and rapidly. This is manifested in a further decrease of the specific surface area and a rapid increase in the permeability. In the regions before large scale homogeneous dissolution can set in, permeability when scaled by the concentration of the reactive species (walkers), shows a power-law variation with porosity. The exponent is 7.6. The reactive constants and the diffusion coefficient are parameters of this study. These may also be varied to fit different experimental observations. Whether the exponent can be related to these parameters or to the rock composition, shall be investigated in a future study. In this work, we have shown how the porosity and permeability of a sedimentary rock changes with time and the role of the concentration of the reactive species. We have attempted to relate these changes to the changes in specific surface area, thereby proposing that the mechanism of 139 dissolution is different for different concentration as the fluid continues its flow. We have compared our results with available experimental results and found favourable comparisons. Since we use the finite difference method to calculate properties, any change in the grid size will cause a scale change of the values though the nature of variation will remain the same. The study is particularly relevant in the problem of CO2 sequestration in sedimentary rocks where CO2 is dissolved in water under high pressure and this is pumped into reservoirs of sedimentary rocks. The partial pressure of CO2 in water will affect the trapping of the gas as a carbonate in limestone reservoirs. This results in the dissolution of the parent rock which changes in porosity and permeability properties. With more and more dissolution, the transported fluid gets saturated with the products of the reaction. This in turn may lead to precipitation on the channel walls. In a future communication we intend to simulate this process of dissolution and precipitation in flow channels using the diffusion reaction process and study their effect on the transport properties of the porous rock. We hope to report these results in a future communication. Acknowledgements The authors are grateful to Marco Dentz and Sujata Tarafdar for their useful suggestions and discussions. T.Dutta is grateful for a grant to IFCPAR CEFIPRA Project No. 4409-1. References Dutta, Tapati., Tarafdar, S., 2003. J. Geophys. Res. 108 (B2), 20–62. http://dx.doi.org/ 10.1029/ 2001JB000523. Edwards, D.A., Schapiro, M., Brenner, M., 1993. Dispersion and reaction in twodimensional model porous media. Phys. Fluids A 5, 837–848. Hornung, U., 1997. Homogenization and Porous Media. Springer-Verlag Inc., New York. Hoshen, J., Kopelman, R., 1976. Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm. Phys. Rev. B 14, 3438. Kang, Q.J. et al., 2002. Lattice Boltzmann simulation of chemical dissolution in porous media. Phys. Rev. E 65 (3). Kang, Q.J. et al., 2005. Numerical modeling of pore-scale phenomena during CO2 sequestration in oceanic sediments. Fuel Proces. Technol. 86 (14–15), 1647– 1665. Kang, Q.J. et al., 2010. Pore scale modeling of reactive transport involved in geologic CO2 sequestration. Trans. Porous Media 82 (1), 197–213. Kang, Q.J., Zhang, D.X., Chen, S.Y., 2003. Simulation of dissolution and precipitation in porous media. J. Geophys. Res. – Solid Earth 108 (B10). Li, L., Peters, C.A., Celia, M.A., 2006. Up scaling geochemical reaction rates using pore-scale network modelling. Adv. Water Resour., 29. Li, L., Peters, C.A., Celia, M.A., 2007. Effects of mineral spatial distribution on reaction rates in porous media. Water Resour. Res. 43, W01419. Litchner, P.C., Tartakovsky, D.M., 2003. Stochastic analysis of effective rate constant for heterogenous reactions. Stoch. Environ. Res. Risk Ass. 17, 419–429. Luquot, L., Gouze, P., 2009. Experimental determination of porosity and permeability changes induced by injection of CO2 into carbonate rocks. Chem. Geol. 265, 148–159. Meile, C., Tuncay, K., 2006. Scale dependence of reaction rates in porous media. Adv. Water Resour. 29 (1), 62–71. Pettijohn, F.J., 1984. Sedimentary Rocks. Harper & Row Publishers Inc., USA. Sadhukhan, S., Dutta, T., Tarafdar, S., 2007. Geophys. J. Int. 169, 1366–1375. Sarkar, S., Toksoz, M.N., Burn, D.R., 2004. Fluid Flow Modelling in Fractures. MIT Earth Resources Laboratory Industry Consortium Meeting. Whitaker, S., 1999. The method of volume averaging. Kluwer Academic Publishers.