Infiltrative instability near topography with implication for the drainage of soluble rocks

We present here numerical modeling of infiltration instability near a topographic edge of a water-saturated porous slice by analogy with a limestone formation devoid of initial heterogeneities such as fractures faults or joints and limited by a vertical cliff. In our runs a first dissolution finger develops near the cliff edge, and ends to intersect it above 5 its mid height. Additional fingers develop upstream with a decreasing growth rate and an increasing width. This results from the decrease of the infiltration velocity with distance to the cliff in our models. A sensitivity study shows that a larger permeability contrast between the fingers and the initial undissolved porous medium produces a larger number of fingers, while increasing the dispersivity (lower Pe) produces wider 10 fingers. A slower reaction rate (lower Da) produces fingers that follow the initial flow lines, since dissolution occurs simultaneously along the entire finger. These results suggest that alteration by dissolution of limestones or other soluble formations may produce different underground channel structures in the same drainage basin due to local changes of the non dimensional Pe and Da numbers. 15


Introduction
When a reactive fluid percolates into a soluble porous rock, dissolution zones concentrate fluid flow, due to their increased permeability, and are therefore subjected to further preferential dissolution. Infiltration instability arises due to this positive feedback between flow and reaction. As a result, the reaction front is unstable and develops so-called dissolution fingers where flow is channeled. The wavelength of the instability can be defined as the distance between successive fingers; it corresponds to the wavelength developing initially in the reaction front. The development of the instability is limited by diffusion, which tends to smooth concentration gradients. Moreover, slow reactions, corresponding to wide reaction fronts can only develop wide fingers, while fast reaction can develop instabilities at any wavelength, which result in finger of any width. The infiltration instability has been studied theoretically by Ortoleva et al. (1987) and Lasaga (1990, 1992), who proved that in the case of limestone, wavelengths larger than 10 µm were unstable and able to form dissolution channels. Therefore, the infiltration instability is in action during karst formation and might be a clue to understand the organization of karstic networks. It has been studied experimentally at the centimetric scale by Daccord and Lenormand (1987) who injected water radially from the axis of a plaster plug, and by Fredd and Fogler (1998), who injected several aggressive fluids in limestone. These authors compare their results with simplified models, involving Diffusion Limited Aggregation for Daccord and Lenormand, and 2-D or 3-D network models for Fredd and Fogler. Steefel and Lasaga (1990), Bekri et al. (1995), Ormond and Ortoleva (2000) and Golfier et al. (2002) used coupled mass transfer numerical models to define the morphology of dissolution channels that develop inside a porous medium, as a function of the non-dimensional parameters of the system, which are the Peclet and Damköhler numbers (hereafter P e and Da). Large P e and Da correspond to the strong instability regime, where fingering may set up on any heterogeneity of the porous medium. For P e<5, dispersion prevents the development of the instability. Moreover, Ormond and Ortoleva (2000) have modeled the fingering process with a Brinckman equation, which allows coupling of Darcy flow in the porous medium with Navier Stokes flow in voids, and compared it with a Darcy model, where dissolution zones were described by a several order of magnitude permeability increase. They conclude that although the propagation of fingers was faster when the entire medium was subject to dissolution, the general wormholing pattern was rather insensitive to the equation used in dissolved area. This suggests that flow focusing and finger competition are not sensitive to the kind of equation used in the finger, once flow is much easier in dissolution zones. Siemers and Dreybrodt (1998) and Kaufmann and Braun (1999) compute the dissolution pattern occurring in a fracture set with realistic boundary conditions including topography. Further publications improve the description of the porous medium, introducing a dual fissures set or a porous matrix exchanging fluid with fractures (Kaufman and Braun, 2000;Gabrovsek et al., 2004). Gabrovsek and Dreybrodt (2001) allow fluctuation of the water table and Kaufmann (2003) considers flow inside the unsaturated zone. These publications consider topography as a driving force for flow. However, their initial fissure network constitutes a set of large-scale heterogeneities that channel dissolution. We propose here that dissolution should be also modeled in an homogeneous porous medium to understand how the instability works by itself in the presence of topography. Our models consider mainly fast reactions, by reference to calcite dissolution, but owing the evidence of pseudo-karstic structures in silicated rocks altered under tropical climate (Trescases, 1975;Wray, 1997;Willems et al., 2002), the influence of a slower dissolution rate is briefly explored. Our results should be viewed as a first step to bridge the gap between infiltration instability models using an initial 1D circulation scheme such as those of Ormond and Ortoleva (2000) and Golfier et al. (2002) and dissolution events occurring at the scale of a geological formation.
The simplifications introduced in the coupled equations of infiltration instability are discussed first. These simplifications allow to develop a compact, robust and efficient numerical model for unstable dissolution of porous media with fast dissolution rate. Then the growth of dissolution structures and a sensitivity analysis are presented. Finally possible consequences for the porosity distribution in carbonated and other soluble rocks are discussed.

Equations
We use equations able to capture the nature of the infiltrative instability and to produce a simple and fast code. We model a water-saturated medium, where the Darcy flow Eq. (1) links the filtration velocity v to the pressure gradient ∇P . Here µ stands for the fluid dynamic viscosity, ρ for its density, k for the porous medium permeability, and g is the gravity.
This equation accounts well for flow in a carbonate-cemented sandstone, where a silicated matrix remains after dissolution of the cement, but is only an approximation inside a carbonate where the whole matrix is soluble. Here, dissolved zones are characterized by a several orders of magnitude permeability increase.
Our reaction-transport equation considers a single chemical species, say H + dissolving the porous matrix with a first order kinetics. The concentration of the active species at equilibrium with limestone is assumed to be negligible. Then, Quintard and Whitaker (1999) and Golfier et al. (2002) have shown that the reactive transport of this species can be described by Eq. (2) where C is the averaged solute concentration of the species.
Here, ϕ is the porosity, the first term depicts the accumulation of the species, the advection velocity in the second term is the filtration velocity, the third term describes mixing of chemical species by a diffusion process of diffusivity D, and the first order kinetics reaction is characterized by the constant α. The use of only one chemical species, of a single diffusivity D for the mixing of species, and of first order kinetics are clearly approximations. However, more sophisticated dispersivity models and reaction kinetics could be easily included. As the time constant needed for Eq. (2) to reach a steady solution is much lower than the one required for a significant change of porosity, the accumulation term may be canceled in the Eq. (2). This point will be further quantified in the section devoted to the numerical model. Moreover, in our models, the amount of fluid advected during a time step is several orders of magnitude larger than the increase of fluid volume due to the increase of porosity by dissolution. This allows also to adopt the steady state approximation of the mass conservation Eq. (3) As a result of the quasi-steady approximation, our equations do not include explicitly time dependence. Time stepping arises only through permeability update during dissolution.
Equations are used in their non-dimensional form (1')-(3'), with the distance scaled by the topographic jump, h, the permeability by k 0 , the permeability before dissolution, the concentration by C 0 , the concentration of the infiltrating fluid, the filtration velocities by v 0 =ρgk 0 /µ which is the velocity induced by the gravity in a vertical slice of the initial porous medium. The non-dimensional pressure is P =(P −ρgz)/ρgh. These non dimensional parameters are valid if the fluid density and viscosity are constant during dissolution, which is the case for low solubility minerals such as limestone and silicates. Fluid density changes, such as those involved during karst development in evaporites, could be considered by introducing explicitly in the Eq. (1) the concentration dependence of density.
Equation (5) includes the non-dimensional numbers P e=ν 0 h/D and Da=αh/ν 0 that are those of the infiltration instability, to which should be added the permeability ratio between the fully dissolved and the non dissolved media (Rk) and the aspect ratio of the modeled area, which is 2 in our models. It can be seen in their definitions that P e measures the ratio between advection and diffusion, while Da measures the ratio of reaction rate by advection.

Geometry and parameters
The initial condition and boundary conditions are pictured on Fig. 1. Our equation set is solved on a 2-dimensional area consisting of a slice of porous medium with a height h of 100 m, an aspect ratio of 2, an initial permeability k 0 of 3×10 −14 m 2 , and a porosity ϕ 0 equal to 0.3, which is assumed to overlie a horizontal impervious medium. The right boundary of this porous slice is a vertical surface from which water is free to escape and that corresponds to a vertical cliff. The constant atmospheric pressure (P =0) is applied there. This corresponds to P =−z/ h. A symmetry condition is applied at the left boundary. The flow lines are curved and convergent toward the right free boundary. Moreover, the infiltration velocity (i.e. the vertical Darcy velocity) is maximum near the free boundary and decreases toward the upstream region, for example from 0.82 at x=0.2 to 0.12 at x=0.99.
A constant concentration of active species C 0 =10 −5 mol/kg is maintained at the upper surface of the model. The characteristic filtration velocity is ν 0 =3×10 −7 m/s, with µ=10 −3 Pas and g=10 ms −2 , and the characteristic time scale is h/ν 0 =3×10 8 s≈10 yr. In our  models the diffusion term describes hydrodynamic dispersion occurring due to the tortuous paths followed by fluid particles inside the porous medium. It operates mainly near the borders of wormholes, where the concentration gradient is perpendicular to the fluid velocity. Hence, the transverse dispersivity is used, which is estimated to 0.3 m in the scale range 10-100 m, according to Gelhar et al., 1992. With the velocity ν 0 , this correspond to a diffusivity D=10 −7 m 2 s −1 which yields P e=300. The constant α has been estimated for limestone by considering that α=K diss S 0 /C 0 where K diss is the kinetic constant for dissolution of calcite, and S 0 is the specific interfacial area between the solid and the fluid, for which a value of 10 2 m 2 /kg has been adopted. The dissolution of calcite is a fast reaction up to a 60% or to a 90% saturation index, according to Palmer (1991) and to Dreybrodt (1996), respectively. Here, we use the constant α deduced from the results of Dreybrodt (1996) up to complete saturation in order to ensure a first order kinetics. This leads to K diss =2×10 −7 molm −2 s −1 which yields Da=7×10 8 . Considering that our fluid velocity is ν 0 /ϕ and that 1/α is the time constant of the dissolution reaction, fluid covers a distance of 5×10 −7 m during reaction, and the reaction front is much thinner than the grid mesh. This allows us to assume an infinitely thin reaction front and thus infinite Da. With the Da of 10 used for the test of a finite dissolution rate, however the front is 10 m thick, and several grid points are out of equilibrium. P e values in the range of 30-3000 are explored by changing the initial permeability of our layer. R k values between 10 and 10 5 are considered. The characteristic time scale is given by h/v 0 and is equal to 1, 10, and 100 yr for P e equal to 3000, 300, and 30, respectively.

Approximations and numerical methods
Two different problems arise in numerical modeling of carbonates dissolution, both of which originate in the sharp  interface between non dissolved limestone and voids created by dissolution. The first one consists in tracking a sharp interface without allowing it to spread on several grid points. The second one consists in accounting the permeability jump arising at this interface.
Assuming instantaneous dissolution allows key simplifications: firstly the time constant required by the mass transfer equation (Eq. 2) to reach steady state depends on advection and diffusion, only. In our models diffusion operates in the fingers and at their border, tending to smooth the concentration of active species inside them (Fig. 2a, and b). Thus, the characteristic time to reach steady state in Eq. (2) is the one required to establish the concentration profile inside the longest and widest finger, which is nearly equal to the time required for the fluid to travel along the whole length of the finger, i.e. of the order of magnitude of 1. It results that steady state in Eq. (2) is reached in a non-dimensional time of 1.
Secondly, the modeled slice can be divided in a non dissolved volume, where no active species may be present, and a fully dissolved volume, where no limestone is present and therefore no reaction may occur. Thus, Eq. (2) can be simplified to yield: Equation (8) expresses that when the limestone is fully dissolved, the transport-reaction equation does not include any reaction term, while Eq. (9) indicates that infinitely fast dissolution reaction brings the concentration of active species to zero as soon as there is still some limestone to dissolve. Since the dissolution front is infinitely thin, the combination of Eqs. (8) and (9) matches the whole volume of the modeled slice.
Then, dissolution reduces to the advection of the infinitely thin porosity front at a velocity γ v N C with γ = 1 (1−ϕ) C 0 C S , v N being the normal velocity at the interface, and (1−ϕ)C s the carbonate volume concentration inside the solid phase (Fig. 2d). Considering a molar mass and density of 100 g and 2.800 kgm −3 for calcite and our concentration of active species in the infiltrating fluid of 10 −5 mol/kg, it comes γ =5×10 −7 . Finally, as γ 1, the interface moves at a much lower velocity than the filtration velocity and the time needed to record a significant porosity change is much larger than 1. This implies that the steady state mass transfer Eq. (5) may be used.
Multigrid methods, known for fast convergence in the presence of steep gradients (Hackbusch, 1985) are employed to solve alternatively Eq. (7) and the combination of Eqs. (8) and (9) on a 257×129 grid with some checks of consistency on a 513×257 grid. Equation (9) is applied as a constraint inside the non dissolved volume. The advection term of Eq. (8) is discretized with the finite-analytic method of Lowry and Li (2002), which prevents spreading by artificial diffusion of concentration gradients. The non diffusive VOF method, developed by Hirt and Nichols (1981) is adopted to track sharp porosity fronts. The time step is adjusted so that less than 1/10 of the volume of a mesh element is dissolved at each iteration. Coupling the set of Eqs. (8) and (9) with Eq. (7) is achieved via permeability changes. The multigrid solver and the VOF method have been checked against analytical solutions. The whole code has been carefully checked against the examples given by Ormond and Ortoleva (2000) in their Figs. 1 and 2.

Results
In our first models, the upper right corner was first dissolved and attracted further dissolution, that was then concentrated at the surface to produce an analog of a river. To constrain vertical infiltration from the surface, a series of vertical high permeability slots, analog to an epikarst, were imposed randomly down to an 8-m depth. No dissolution was allowed inside these slots as well as on a 15 m length at the upper right corner. This does not imply that no dissolution occurs inside the epikarst, but that here dissolution results in vertically organized structures that allow infiltration of rainwater. For easier comparison, the distribution of the vertical slots is identical in all the simulations. However several initial distributions of vertical slots have been checked and the influence of this parameter is discussed in the forthcoming sections.

The reference case
Our reference case is characterized by Rk=100, P e=300 and infinite Da. A first finger develops at x=0.25 in a direction nearly parallel to the initial flow lines. The finger curves in the immediate surrounding of the free boundary and ends crossing the cliff at a z=0.65 after a time lapse of 2×10 5 (Fig. 3a). The position of this first finger depends slightly on the initial permeability distribution. After the breakthrough, neither the porous flow nor the water saturated medium assumptions are valid and dissolution is no longer computed in the finger. The pressure in this finger is almost equal to the atmospheric pressure due to the large permeability ratio. This induces on subsequent fingers an attraction effect similar to that of the right topographic edge.
The second finger grows at x=0.7 and reaches the cliff at z=0.3 after a non-dimensional time lapse of 2.3×10 6 (Fig. 3b). The development of this second finger is significantly slower than the first one even considering its larger length before breakthrough. It is slightly wider than the first one. This can be viewed as a consequence of the decreasing initial fluid velocity in the upstream direction resulting in a decrease of the local Peclet number and an increasing smoothing of aggressive species inside the finger due to diffusion. Thus, only wider fingers can develop, an effect that has already been noted by Steefel and Lasaga (1990).
Two fingers that set on at x=1.1 and x=1.2 are prevented from developing by hydrodynamic capture of the aggressive fluid by the second finger through a superficial horizontal dissolution drain (Fig. 3c). The third finger is wider than the previous ones and develops after a longer while. It follows firstly the initial flow lines and then is captured and converges toward the second finger after a non-dimensional time lapse of 9.8×10 7 (Fig. 3d).
Simulations with various initial vertical high permeability slots distributions gave slightly different finger distributions all with the same common characteristics: (i) a first finger arrives above mid cliff height after a short time lapse, (ii) 2 or 3 additional fingers reach the cliff near z=0.3, their width and the duration of their growth increasing in the upstream direction, (iii) fingers growing near another fully developed finger are stopped by hydrodynamic capture.

Sensitivity analysis on Rk and P e
As a common feature, the first finger (the rightmost one) intersects the cliff above mid-height. Then a second finger or a series of fingers arrives slightly above the base of the limestone formation for the highest Peclet values. This implies that the base of the formation is not drained by high permeability channels, a drawback of our modelling technique in which the development of a finger is stopped as soon as it arrives at the cliff. Therefore, further dissolution cannot develop from the base of this finger toward the base of the cliff. Smaller Peclet numbers (Fig. 4a and c) allow only the development of wide fingers. This results from diffusion, that smoothes the concentration of active species inside thin fingers to almost zero. The second finger of Fig. 4a and C develops slowly and is at least 200 m wide. At x=0.7 the infiltration velocity is only 0.5 and the resulting local Peclet number is reduced to 15, near the value of 5 under which no fingers are allowed to develop for the high permeability ratios of this study (Ortoleva, et al., 1987;Steefel and Lasaga, 1990). The local Peclet decreases in the upstream direction and reaches 5 at x=1.1. Figure 4a and C show then that dissolution produces large pockets around the widest initial vertical high permeability slots. Probably, these pockets could not develop without an initial high permeability zone able to focus the aggressive fluid flow. They could be compared to dolines that are common dissolution features of karstic landscape (Ford and Williams, 1994).
The runs of Fig. 4a and b are characterized by the largest permeability contrast (Rk=10 5 ) and thus by enhanced hydrodynamic capture of reactive fluid that prevents the simultaneous development of neighbour fingers, an effect that has been previously noted by Ormond and Ortoleva (2000). This is best evidenced in Fig. 4b where the finger at x=0.95 devel-ops on almost half the cliff height before it is stopped by the progression of an horizontal channel that originates from the first finger. Figure 4a shows that after the development of the first finger, only large dissolution pockets grow around large initial permeability zones. They do not evolve as fingers even until 30% of the initial limestone is dissolved. We suspect however that in real geological examples, due to superficial erosion and run-off, rainwater could be focused at the center of these pockets and that instable dissolution could set on as a result of the increased infiltration velocity and Peclet number. Finally the case of Fig. 3d is characterized by the development of several dissolution fingers, as a consequence of the high Peclet number, that enhances instability, and of the low permeability contrast that limits the competition between fingers by hydrodynamic capture.

Effect of a finite dissolution rate
In order to discuss the case of finite dissolution rate of silicates and carbonates near their equilibrium, a model with Da=10, corresponding to a reaction time constant 1/α=1 yr, has been included in our sensitivity analysis, the other parameters being left unchanged from our reference case. The where no additional finger is allowed to develop, is drawn. The permeability ratio between dissolved and non dissolved zones increases from bottom to top, and the P e increases from left to right. The white rectangle labeled as "reference case" represents the final state of the reference case (Fig. 3d). As in Fig. 2, time is in non-dimensional units.
permeability is modeled with a Kozeni-Karman type (i.e. k proportional to φ 3 /(1-φ) 2 ) law, up to a porosity of 0.6, and beyond that value tends asymptotically to 100. The accumulation term of Eq. (2) can still be neglected according to the results of Ogata (1964Ogata ( , 1970 and Lichtner (1988). The dissolution rate is given at each grid point by the term P e×Da C . Figure 5a shows a diffuse dissolution zone of width close to 0.5 developing nearly along the initial flow lines and reaching the cliff near its mid-height. This can be viewed as a consequence of the reaction front width, that is between 0.1 and 0.5 due to the initial infiltration instability velocity of 1 to 5 in the area where the dissolution develops. Thus the reaction front will include a large part of the porous slice thickness. Our results are consistent with those of Golfier et al. (2002) for the style of instability and those of Ormond and Ortoleva (2000) who note that only fingers wider than the dissolution front are allowed to grow. Moreover, due to the slow dissolution kinetics, a finger develops simultaneously on its whole length and therefore is driven in orientation by the initial flow lines. We have noted that this effect was dependent on the porosity-permeability law and was more pronounced if a steep permeability change was introduced at the beginning of dissolution.
One notes also in Fig. 5a the onset of a second dissolution finger at x=1.05, sharper than the first one, due to the sharper dissolution front resulting from the decrease of the infiltration velocity in the upstream direction. This finger develops also along the initial flow lines and reaches the base of the cliff at t=7×10 5 (Fig. 5b).
Thus, a slow kinetics reaction tends to induce dissolution parallel to the initial flow lines and, as in the previous section, the instability style may change as a consequence of the decreasing infiltration velocity and Damköhler away from the cliff edge.

Discussion and conclusion
Assuming instantaneous dissolution of limestone allows key simplifications in infiltrative instability modeling. They are implemented to explore the effect of a lateral topographic jump on the fingering process. Two main results arise: (i) the first finger develops very near the vertical free boundary and crosses it above its mid-height, while a series of fingers develops progressively in the upstream direction. The drainage of the entire slice length may take up to 100 times the time for breakthrough of the first finger (ii) the dissolution style varies with the infiltration velocity, wider fingers a) T=3.8 E5 b) T=7.0 E5 2  are generally produced away from the jump, or fingers are replaced by large dissolution pockets. These effects are the consequence of the decrease of the infiltration velocity and of the resulting decrease of the Peclet number. They were hidden in previous quantitative models of karst generation since dissolution was channeled in permeability heterogeneities introduced by the initial fracture set of these models.
Considering the approximations introduced in the equations, our results have to be compared with great caution to real geological examples. A major assumption is to consider that the porous medium is saturated with water. Thus, the model assumes the infiltration of 10 myr −1 of rainwater for a non-dimensional vertical velocity of 1, a value that can hardly be achieved, even in most rainy climates. Moreover, our models differ primarily from real karst by the amount of connected porosity, since one of our runs produces a connected network with a near 30% porosity, while a fully mature karst could work with a 1% porosity (Mangin, 1995), including large caves and shafts poorly connected with the main drainage network. Further insight into the geometry of the fingering process might be gained by using 3-dimensional modelling, which is achievable with numerical method similar to those of this study and up-to-date computers.
In the case of silicate rocks altered by dissolution and thus possibly subjected to the infiltration instability, two additional effects should be considered. First, incomplete dissolution will result in less permeable channels than in the karst case, which has been shown here to produce the simultaneous growth of numerous fingers. A second effect will arise from the slow dissolution kinetics that will produce a thick dissolution front possibly including the entire formation thickness and prohibiting the development of the instability. However, a thinner dissolution front may result from the decrease of the fluid velocity in some areas of the infiltration basin, so that the infiltration instability could set on there, under the con-dition that the Peclet number has not decreased below 5. In these formations, the infiltration instability is likely not to be detected, since it might govern the permeability distribution in some places only.
In summary, modeling infiltration instability could yield qualitative tools for better understanding the permeability distribution and its time-evolution inside a formation altered by dissolution. In the light of the present results we propose further modeling efforts accounting for dissolution inside the unsaturated zone, 3d aspects and variable density effects in evaporites.