Existence of convective threshold and its role on temperature in reactive flow through fractured rocks: a simulation study in 2D

A reactive fluid flowing through porous or fractured rocks causes changes in its porosity and permeability. Experimental studies on the effects of reactive flow in porous systems are normally confined to finite size rock samples in laboratories as actual field data is not sufficiently available. Such laboratory data and simulations based on them, may not always correctly predict or explain changes that take place in rock basins in geological time. In this study we simulate two cases- transient and steady. In a transient case, reactions occur as soon as a reactive fluid is injected into the rock, whereas in the steady case, reactions begin only after the injected fluid percolates through the sample length. The transient case is comparable to real rock basins, while the steady case mimics the situation in laboratory experiments. Variation of the hydrodynamic parameters such as porosity, permeability and surface area, and their inter-relationships at different temperatures, are discussed in both cases. We establish that for every temperature, there exists a threshold flux above which the two cases give similar results. Prediction of changes in rock property due to reactive flow in real situations can be made from laboratory experiments provided the flux is kept above this threshold value. When the fluid flux is kept above this threshold value, the porosity-permeability relation follows a power law behaviour. While the threshold flux is independent of temperature and channel width, the exponent of the power law is a slow decreasing function of temperature.


Introduction
Much research has been done in modelling reactive flow in connection to geothermal energy extraction, storage of nuclear waste and CO 2 sequestration [1][2][3][4][5][6][7]. The subsurface rock usually contains minerals that can dissolve in a suitable reactive fluid that flows through its pore space. Under certain conditions, precipitation of solutes may occur on the pore walls or through nucleation in solution. Thus dissolution and precipitation processes cause changes in the porosity and therefore permeability of the rock. These chemical reactions occur only as long as the the fluid is not in chemical equilibrium. In equilibrium, the amount of dissolution is offset by an equal amount of precipitation and porosity reaches a constant value. Normally equilibrium is dependent on temperature as the solubility of the mineral ions often increases with temperature. The factors that govern reactive fluid flow like density, viscosity, diffusion constant and kinetic rate coefficient, are also all temperature dependent. Therefore it can be expected that a rock of a certain mineral composition in the Arctic region will show porosity and permeability changes with time in a manner quite different from a similar rock sample found in the Equatorial Belt, even though the same reactive fluid may be flowing through it.
Several experimental and theoretical studies investigating the influence of dissolution processes of porous media have been reported [8][9][10][11][12][13][14][15][16][17][18][19][20]. In earlier studies [21,22], the authors have reported the role of concentration and pressure head of the fluid in the study of dissolution and precipitation processes in rocks by using Random Walk Methodology. One of the problems that the researchers face, is the absence of sufficient field data on the effects of reactive flow in real rock basins. Laboratory experiments are performed on finite size samples, and the results are used to predict the outcome of situations in fields. This extrapolation of knowledge may not always yield the expected results. One of the ways to solve this problem is to make the defining equations dimensionless, but the choice of the normalising factor is not always standard. In this work, we try to bridge the gap between results performed in the laboratory with the results in real rock basins by studying reactive flow in porous rocks under two different circumstances: transient and steady. In the transient case, chemical dissolution is allowed as soon as the reactive fluid is injected into the rock. In this case therefore, the rock structure as 'observed' by the reactive plume head, may have been changed due to chemical reactions at earlier times, and the reactive fluid train is exposed to the changed structure. This is the situation observed in actual rock basins. In the laboratories, the sample size being small, time for the reactive fluid to percolate through the sample is small too. Thus reactions may be assumed to begin approximately at the same time in the entire connected pore space of the sample. This situation is modelled by the steady case where reactions between fluid and minerals are allowed only after the plume of the injected fluid has percolated through the sample. Our aim here is to investigate under what circumstances, if any, do results of the two cases match.
Our model considers the feedback between flow and chemical reaction. For studying the behaviour of such coupled processes, numerical modelling still remains an indispensable tool. Dissolution can be modelled at various length scales starting from pore scale to core scale. We consider dissolution at pore scale where the local chemical kinetics and transport of the solute, control dissolution.
We study the impact of the variation of temperature in dissolution patterns evolving due to a continuous injection of a reactive fluid through a fractured rock porous media, modelled here as a 2-dimensional channel. The changes in porosity, permeability and interface area are tracked with time for different temperatures ranging from 1°C to 75°C , for both the cases. We have discussed the condition under which there is the possibility of a transition from one case to the other. Only under such conditions can laboratory results be used to predict realistic expectations in actual field situations.
In the sections to follow, we detail the theory and simulation algorithm, present our results for both transient and steady conditions, and finally discuss and conclude our study. We have taken the case of CO 2 sequestration in sedimentary rocks as a model for our studies in the sense that all the parameters of our simulation correspond to the reactions of H + and -CO 3 2 ions with forsterite mineral.

The simulation algorithm
The algorithm can be broadly divided into three modules. Firstly, we identify the connected pore space and solve for steady flow through it. Secondly, we solve for transport of the reactive solute along with the fluid. Thirdly, we allow dissolution of solid matrix on reaction with the injected solute in accordance with rules set up from chemical equations. The detailed procedure is given below.

Fluid flow
We discretize the pore space into square grids. On this mesh we solve for steady flow under a constant pressure difference. We implement no slip at the rock-pore interfaces by superimposing imaginary grids on the rock layer at the interface [23,24]. An imaginary layer of grids is also added to the inlet and outlet. The fluid is assumed to enter the sample normally from one side, and after traversing the porous pathway, ejects out normally through the opposite side through these cells [25]. The tangential components of velocity are set zero in these cells. Incorporating this condition in the equation of continuity enables one to obtain suitable expressions for pressure and velocity in the pore pixels of the system.
At the pore scale, flow is governed by the Navier-Stokes equation where ρ is the fluid density, v the flow velocity, p the pressure, η the constant dynamic viscosity, and f e represents the body forces. The fluid considered is incompressible, viscous and Newtonian in nature. The pressure gradient is low enough to yield a slow moving fluid. Hence we disregard the inertial term in the above equation. Furthermore, we assume that there are no external forces, and flow is quasi-steady. With these assumptions, equation (1) simplifies to the Stokes equation, Fluid mass conservation for the incompressible fluid is expressed by Equation (2) is discretized in an implicit scheme such that Using equations (3) and (4), we can write down an equation for pressure at step (n+1) Equations (5) and (6) are then spatially discretized before we finally employ them in our simulations [24,25]. Starting with zero velocity as the initial condition, equations (5) and (6) are iterated until the steady state is achieved. The steady state is reached when the change of the average velocities over several iterations is less than some pre-defined value, ò=10 −6 . Needless to say that lower the value of ò, better is the mass-balance. However it must be borne in mind that beyond a certain limit depending on the configuration of the system, the change over the iterations is insignificant and time consuming. Therefore the limit has to be judiciously set so that while the mass-balance error is small, the solver is economical in time. In all our simulations we have assured that the relative mass balance error is below 10 −2 . Permeability κ is determined from the Darcy's law as where Q is the volumetric fluid flux per unit time, under pressure gradient p, and A is the cross-sectional area of the rectangular flow domain.

Solute transport by time domain random walk (TDRW)
We solve for solute transport in the system by using the advection-dispersion equation [26]. Random walkers (RW) mimic the role of solute particles in our simulation. We maintain constant concentration at the inlet. For the movement of the RW, we implement the time domain random walk (TDRW) technique [27,28]. The technique involves a RW to jump from one site to another such that the distance traversed is fixed whereas the transition time is variable depending on the local velocity. We broadly employ the approach discussed by [29] where the probability of transition to the next neighbour pixel and the corresponding time for jump, is formulated from the advection-dispersion equation, Here C is the concentration, D the dispersion coefficient and V is the advection velocity. Finite difference discretization of equation (8) in 2-dimensions gives ¶ ¶ = ---- where, C i is the concentration at the i th pixel. Terms C up x , C dn x denote the concentration values at the next neighbour pixels flanking pixel i along x axis, the magnitude of concentration being higher (up=upstream) and lower (dn=downstream) than C i respectively. The notations C up y and C dn y carry similar meanings corresponding to the y axis, u and v denote the components of V along x and y axes respectively, h is the distance between two pixels. Rearranging the terms in equation (9) yields ¶ ¶ = -+ + + where, The above equation can be generalized to For any transition between neighbours, say from j to i, the transition probability W ij and the corresponding time are given as [29]

Chemical dissolution
We have considered first-order kinetic dissolution mechanism at the fluid-solid interface. The rate law for the first order is given by k where k is the kinetic coefficient and σ the specific reactive surface area. In the particle based formulation outlined in the previous section, this reaction is modelled on the Gillespie algorithm [30]. In this framework, chemical reactions at all the rock-pore interfaces, if allowed, occur simultaneously in a single reaction step that corresponds to a single time step of simulation. In any one reaction step, exactly one reaction event-the one with the shortest reaction waiting time, can occur at only those pixels that are located at the solid-fluid boundary. A reaction at any rock pixel stops with the dissolution of that pixel.
The RW at the interface of a rock pixel, representing a reactive solute ion, has a waiting time before it may react. The reaction waiting time for a single RW r as per the probability distribution function is, Thus, the reaction waiting time of a solute particle r out of an ensemble of n solute particles is given by, The probability that any solute particle r reacts at time τ r residence time τ, is obtained by integrating equation (16), and is given by We have assumed τ r =τ in our simulations, and the RW is transported as per equation (12). To determine if the RW present in the pore pixels adjacent to the interface succeeds in dissolving the neighbouring rock pixel, we draw a random number from the uniform distribution 0<=r<1 and compare it with Prob. If r<Prob, dissolution of the rock pixel is allowed. The rock pixel disappears and is assumed to be flushed away with the flowing fluid.

Results
To validate the code for solute transport mechanism we compared the analytical results with the simulation results for the concentration distribution at various times in a 2-dimensional rectangular homogeneous domain of unit length. A schematic of the system is shown in figure 1(a). The fluid velocity v is constant from right to left. The initial condition is that the concentration . The solution of equation (8) can be expressed in Laplace space [31] as, Figure 1(b) shows the comparison of the analytical and simulation curves for the concentration profile with distance from the inlet, at different times after injection. The coincidence of the graphs validate the transport model in our simulation.
Dissolution studies are done for two cases -transient case where the experiment starts with the injection of the fluid at the inlet of the sample, and the steady case where the experiment starts once the plume has reached the outlet of the sample. These two cases represent the situation observed in real rock basins in nature, and in laboratory set-ups respectively. From the experiments carried out in situ, scientists try to predict the scenario on a large scale. In this paper we try to analyse to what extent the two cases are comparable. In the following section we put forth studies done on the transient and the steady case, and then present a comparison between the two.
Fracture channels that serve as conduits for reactive fluid flow in sedimentary rocks are complex having a hydraulic radius that can vary widely along the channel throats. To simply this complexity of channel throat width, we restrict our initial study to a relatively simple porous sample that consists of channels of various pore throat sizes: 15 and 20 units. The dimensions of each structure are 311×142 pixels. Size of a pixel is taken to be 5 microns. Simulations were done for temperatures lying within the range of 1°C to 75°C. Parameters such as viscosity and density which are temperature dependent, have been assigned appropriate values using the data from the webpage [32], and are displayed in table 1.
Using the Stokes-Einstein equation, and known diffusivity value at 25°C, the diffusivity values for the other temperatures have been calculated.  In equation (20), k B is the Boltzmann's constant and R is the radius of the diffusing particle which is assumed to be a constant here. The relative values for kinetic rate constants have been calculated using the Arrhenius equation and known value of the activation energy for forsterite at 25°C [33].
The respectively. Increased inlet dissolution is observed with increase in temperature for the transient case whereas in the steady case, there is more homogeneous dissolution throughout the channel length. However in the steady case, the extent of dissolution at the inlet remains almost the same as in the transient case for each temperature. As the diffusion coefficient increases and viscosity of the solvent decreases with increase in temperature, equation (20), the RW representing the reactive ions in the solvent, can access the rock interface quicker. The increase in kinetic rate constant with temperature, further aids the completion of dissolution at higher temperature. This results in a conical shape of the channel.
In the steady case, dissolution starts only after the plume arrives at the outlet. The evolving pore structure shows distinctive features that are characteristic of the temperature. At lower temperatures of 1°C-25°C, the channel exhibits homogeneous dissolution along the entire length. At higher temperatures, there is greater inlet dissolution and the channel develops a conical appearance, figures 2(g)-(h) and 2(o)-(p). Thus the structure evolves quite differently for the two cases studied. A comparison of the images for the different channel widths show that though there is greater net dissolution in the wider channel, the dissolution pattern in both the cases is almost independent of the variation in channel width.
Advection is the main driving force propelling the reactive species forward in the fluid through the structure. When constant flux is maintained, advection is the same at all temperatures (equation (7)). Diffusion aides in distribution of the reactive species brought forward by advection to the pore-rock interface lying on either side of the channel. Once the particle is within the system, the direction in which the particle moves and its time of residence at the current pixel is decided by the local advective velocities and the diffusivity. The effect of , is determined through the competition between the kinetic rate k and the local velocity v, as both increase with increase in temperature. c 0 and n represent the initial concentration and the reaction order respectively. If the kinetic rate constant is higher than the local flow velocity, the reaction is mass-transferlimited, while a higher v compared to k, makes the reaction reaction-rate-limited.
At a lower temperature, if the Pe is low to moderate, a reactive ion lying adjacent to a rock pixel, has a greater time of residence in the cell. This in turn results in an increase in the probability of dissolution (equation (17)). Hence though inlet dissolution is observed for lower temperature in the transient case, there is not significant penetration of the reactive plume in the channel. Figure 3 shows the porosity development with time, for both transient and steady cases, in the 2 channels for the different temperatures studied. For the transient case, figures 3 (c)-(d), show an increase in the rate of porosity changes with time. The increase in the net advective velocity with temperature, causes the reactive fluid plume to penetrate deeper into the channel at higher temperatures. The dissolution process is aided because of an increase in the D and k values. As opposed to the transient case, the porosity versus time curves for the steady case show increase in absolute porosity but decline in the rate of increase with time, figures 3(a)-(b). In the steady case, the reactive species is already uniformly distributed in the channel before the reactions begin. If the time of residence of a reactive species in a cell is smaller than the rate of reaction, i.e. dissolution is reaction-rate-limited, much of the reactive species can get flushed out of the system without completing a reaction. With increase in temperature, the kinetic rate constant increases too. The sharp increase in porosity at 75°C indicates that despite the increase in the net advective velocity, there is a net increase in the Damkholer number D. In the steady case, the porosity rate increases initially, and then tends to saturate. It is clear that the transient and steady cases show different behaviour in changes in porosity with time. Figure 4 shows permeability gradually increases with time at all temperatures for the transient (c) and (d) and thesteady (a) and (b) cases, for the two different channel widths. In both cases, at the lower temperatures, 1°C-25°C, there is very slow increase in permeability with time. The rate of increase in permeability is high at 75°C in both cases, and may be attributed to the increase in the kinetic rate constant k and diffusion coefficient D, with temperature. The increased diffusion coefficient transports the reactive ion to the rock interface where the dissolution reaction is completed quickly. However the increase in the rate of permeability changes is smooth in the transient case and jerky in the steady case. This may be indicative of the different manners of dissolution in the system. Except for a greater value of the permeability in the wider channel, there is not much difference due to variation in channel widths.
The permeability versus porosity curves, figure 5, exhibit a quadratic relationship for both the transient and the steady cases. The temperature effect on the nature of the relationship is not very pronounced as the curves for the different temperatures almost seem to coincide. The increase in permeability with porosity increases drastically with temperature. The substantial broadening of the channel width, figure 2, is responsible for the sharp increase in the curves at 75°C.  The surface area versus permeability graphs, (figure 6), show that at low temperatures of 1°C-25°C, there is an increase in surface area with practically very little increase in permeability. We can infer the mode of dissolution that can give rise to such a situation.
As dissolution starts, the porosity increases, and so does the surface area with increased serration of the smooth walls, figure 7(a). It is observed that the surface area shoots up with the onset of dissolution while the permeability remains almost unchanged. During this phase random dissolution occurs along the sides which increases surface area substantially without any effective increase in pore radius. Hence the permeability remains unchanged. This mode of dissolution causes dead ends in the flow path and does not affect the permeability though the total surface area increases. At low temperatures, 1°C-25°C, the net advective velocity is low, reaction is mass-transfer limited, and dissolution is basically concentrated around the inlet of the channel. With increase in temperature and the net advective velocity, the reactive plume advances further in the channel. New pores are created by dissolution of rock pixels lying adjacent to pores created by dissolution in previous timesteps. From the schematic diagram shown in figures 7(b) and (c) it is easily perceived that there is no new augmentation of surface area in most incidents. Due to the growing inlets, the pores fuse in the later stages drastically bringing down the surface area, figure 7(d). Although the dissolution in the neighbourhood of the newly created pores persists which may bring down the surface area or keep it unchanged, with the quick advancement of the plume in the advective direction, new pores are created along the smooth channel.

Discussion
We have considered reactive fluid flow through a simple porous structure consisting of channels of varied pore radii in 2-dimensions. The reactive fluid causes dissolution of the solid matrix. Simulations were conducted on this structure for a wide range of temperatures. Results having different initial conditions were compared as in: (i) the transient case where dissolution starts with the injection of the fluid at the inlet, (ii) the steady case, where dissolution starts once the reactive plume arrives at the outlet. Although we have considered two distinct cases, our conjecture is that in nature an intermediate process occurs. To understand the changes that happen in the geological time scale from the time scale of our laboratory experiments, the dissolution process has been hastened by a suitable reaction velocity in our simulation. The disadvantage added due to this modification is that the reactive fluid may flow through a porous structure that is drastically distorted from its original form as it advances through the channel. In nature, the changes are very gradual and the evolving structure can be traced back to the original one. This very gradual change over geological time scales, is captured in the steady case of the problem where we have assumed that the injected fluid is not subjected to modifications in rock structures due to previous reactions. The two cases of our study may be summed up thus: the transient case shall provide information on what a given structure may look like in large time scales, whereas the steady case shall enable us to understand the gradual change in the structure that leads on to the final changes observed.
The variation of the hydrodynamic parameters with time and their inter-relationships were compared in both the cases. The higher rates of porosity and permeability changes at high temperature (75°C) compared to low temperature (1°C-25°C) is not surprising. This can be understood when we remember that increase in temperature increases the diffusive velocity and kinetic rate constant, and decreases the viscosity. Pronounced inlet dissolution is observed, specially in the narrower channel, in the transient case.
In most desktop experiments small sample sizes are involved. Therefore observation times are quite large compared to the time required by the reactive fluid to fully invade through these samples. This situation is similar to the steady case of our simulation study as the entire structure is affected almost simultaneously by the reactive flow.
In sedimentary rock basins however, the situation is closer to the transient case.
To investigate under what circumstances the results from steady case can be used to make predictions in the transient case, we ran the simulations for both the cases for higher fluxes.
When the time development of porosity and permeability for the two cases look similar the prediction of results from the steady case to transient case is signalled. We ran the simulations with a higher constant flux, ∼5.0 cm/s, by increasing the pressure gradient across the sample. Considering the anomalous expansion of water at 4°C, we ran simulations at 4°C, to check if there was any signature of this behaviour in the properties studied. Contrary to our expectations we have not observed any deviation from the trends observed. A possible reason is that the change in density of water is not much compared to the changes in the other hydrodynamic quantities.
We examine the poro-perm relations for the transient and steady cases for higher flux for all the temperatures to find a power law dependence. Figures 11 and 12 display the power law behaviour of the poro-perm graphs for the temperatures 25 and 75°C, with similar behaviour obtained for 1 and 50°C too. The exponent of the power law is approximately the same for both the channel widths studied and for the temperatures 1°C, 25°C and 50°C, for both the transient and steady cases. For 75°C, the value of the exponent decreases but remains the same for both the transient and steady cases. We may conjecture that as the flux is increased, the transient and steady case results converge.
4.1. The case of irregular fracture channels Fracture channels in sedimentary rock basins generally have complex fracture channels that act as conduits to reactive fluid flow. To examine whether there exists a convective flux threshold for flow in complex rock structure, we used the Relaxed Bidisperse Ballistic Deposition Model (RBBDM)to simulate a 2-dimensional rock structure of size 160×160 pixels. The details of this model is discussed in [25]. The RBBDM involves two sizes of grains to be dropped randomly on a substrate line with probability p and (1−p) respectively. The grains are allowed to relax to their minimum potential energy after deposition. By varying p, a porous structure of appropriate porosity may be generated. In order to study reactive transport, a structure spanning fracture was identified using the Hoshen-Kopelman algorithm. A reactive fluid was injected at different fluxes and at different temperatures as before. The simulation was done for both the transient and steady state conditions as discussed earlier. Figure 13 shows the evolved fracture channel after 1000 time steps for the temperatures 25, 50 and 70°C, for both transient and steady cases, for flux ∼14.5 cm/s. It is apparent from the figure that the channels evolve differently with temperature, and are different for transient and steady cases. In comparison, figure 14 shows the time evolution of the channels for different temperatures for both the cases at the threshold flux of ∼87.0 cm/s. In this case, for every temperature, the final channel has a similar appearance for both the Figure 8. Typical images of pore structure for high advective velocity; grey denote pores pixels, white denote the dissolved rock pixels while black denote the other rock pixels. Transient and steady cases are indicated for the two channel widths indicated. Temperature for each case is indicated at the bottom of the respective column with unit in degree Celsius. transient and steady cases. The corresponding poro-perm graphs for different temperatures, for both the transient and steady states, is displayed in figure 15. As for the straight channels of different widths, we note that the poro-perm relationship has a power law dependence that is characterized by an identical exponent of value m∼3.4 for both steady and transient cases. The relationship is almost independent of temperature.  Our study reveals that there exists a threshold flux Q th above which the steady and transient cases give similar results. Above this threshold the poro-perm relations show a power law behaviour. While this threshold Q th is seen to be independent of channel width and temperature, the exponent of the power law shows a slow decrease with increase in temperature. The results of our simulation are based on temperature dependent parameters like density, viscosity, reaction rate constants and diffusivity values of reaction between H CO 2 3 and forsterite.

Conclusions
To conclude, temperature affects reactive flow quite prominently, whether in the transient or the steady case. The mode of dissolution changes as one goes from low to high temperature. It is however independent of channel width as proved from the initial study on different channels of finite widths. As long as the injection flux is below a critical value, the dissolution pattern develops differently for different temperatures. There exists a threshold  flux Q th beyond which desktop experiments involving reactive flow if performed, may be useful in predicting morphological changes that occur real rock basins as the dissolution pattern shows similar behaviour. Above this threshold flux, the poro-perm relation is a power law with the exponent being a slow decreasing function of increasing temperature. We intend to test this hypothesis further on more complex real rock porous structures, whose reports we hope to communicate in the near future.