Microfluidic measurement of the dissolution rate of gypsum in water using the reactive infiltration-instability

. We present an original method for measuring the intrinsic dissolution rate of gypsum. We use a simple microfluidic setup, with a gypsum block inserted between two polycarbonate plates, which is dissolved by water. By changing the flow rate and the distance between the plates, we can scan a wide range of Péclet and Damköhler numbers, characterizing the relative magnitude of advection, diffusion and reaction in the system. We find the dissolution to be unstable, with a formation of a characteristic fingering pattern. The dissolution rate can then be calculated from the initial wavelength of this pattern. Alternatively, it can also be estimated from the time it takes for the gypsum chip to get completely dissolved near the inlet channel. The method presented here is general and can be used to assess the dissolution rates of other minerals.


Introduction
One of the main challenges in geochemical modeling is the determination of the rates of reaction of minerals. When a mineral is subjected to a flow of dissolving fluid, the dissolution reaction is not only kinetically controlled by the rate of detachment of the ions from the reactive surface, but also by the ability of the aqueous solution to transport efficiently the products of the reaction away from the surface. This strong coupling between chemistry and transport complicates the extraction of the intrinsic reaction rate from the measurable data. The measurements of gypsum dissolution rate by Raines and Dewers [1] show well the difficulty associated with the measurement and with the diffusion boundary layer. They used a rotating disk set-up combined with a mixed flow reactor and proposed in opposition to the commonly used linear dissolution law, a quadratic law with an intrinsic reaction rate of 3.5×10 −3 mol/m 2 /sec. But their result has been criticized [2,3] for issues in the methodology and especially in the estimation of the diffusion boundary layer. Indeed, they used the thickness of this layer as a fitting parameter for their reaction law and obtained values, which are much lower than the theoretical ones [2]. The difficulty of the measurement explains why the literature actually exhibits numerous values for the intrinsic reaction rate of gypsum dissolution. For example, the 2004 USGS report [4] cites Jeschke et al. [3] with a value of 1.6×10 −3 mol/m −2 /sec obtained using both the rotating disk and batch experiment, but also cites Raines and Dewers [1], four years after Dreybrodt criticisms [2]. A very successful attempt to reconcile literature was put in place by Colombani [5]. Considering that the discrepancy in literature stems from the hydrodynamics associated with different experimental set-ups and the estimation of the diffusion boundary layer, he reassessed all the literature data and came up with the consolidated value of 7×10 −5 mol/m −2 /sec, which is consistent with his result obtained by interferometry (5×10 −5 mol/m −2 /sec for the (010) cleavage plane of a single crystal of gypsum). This conclusion is also corroborated by the experiments of Mbogoro et al. [6], who found a value of (5.7±1.4)×10 −5 mol/m −2 /sec.
In this article, we present an original method for measuring the intrinsic reaction rate of mineral dissolution. Instead of monitoring the concentration of ions in the fluid and calculate the reaction rate through the diffusion boundary layer -as was done before in the literature -we take a different approach and calculate the dissolution rate from the pattern carved by reactive fluid in a dissolving mineral.

Methodology
Let us begin by shortly recapitulating gypsum reactivity. Gypsum dissolves in water according to: (1) Following the transition state theory, the intrinsic reaction rate i.e. the rate of detachment of atoms from the solid can be written as where k0 is the intrinsic reaction rate, sr is the specific reactive surface area, a is the activity of the solid (taken equal to 1), η is a parameter and Ω is the saturation defined as . Here, ax stands for the activity of species X, and K is the thermodynamic equilibrium constant. Expressing the activities through the molalities and the activity coefficients, we obtain the following expression for the saturation is the mean activity coefficient. Since the solubility of gypsum is quite low, the ionic strengths even at saturation are relatively moderate and we can consider that the solution is following an ideal behavior [5]. This means that the mean activity coefficient can be approximated by unity as well as the activity of water (solvent). As a result, assuming as in [4,5] that the reaction rate is linear with molality, we have the following relation for the reaction rate, assuming that molalities can be expressed as the volumetric concentration (density variation of water is also negligible): . The reactive surface area in our experiment is actually considered to be equal to the geometric surface area as argued below.
The microfluidic set-up used in the experiments has been described in [7]. It consists of two disks of 6.5 cm diameter and 1 cm thickness made of polycarbonate. The bottom plate contains a rectangular indentation (3.3 cm x 3.8 cm x 0.5 mm) engraved using the MSG4025 CNC micro-milling machine. This indentation serves as a cast for the plaster that is used as a soluble material in our system. The top plate engraved by the same machine contains a hierarchical cascade of channels connected to large inlet/outlet reservoirs (4.5 cm x 5 mm x 2 mm). Their large sizes help in keeping the uniform pressure along the edges of the system. In between the two plates a spacer made out of the ultra-thin PET-based double-coated tape of a thickness of h0=100μm (ORABOND 1394TM from ORAFOL GmbH), h0=70μm (ORABOND 1398) or h0=210μm (ORABOND 1397TM) creates the aperture of the cell. Experiments were performed at room temperature. The cast was prepared with a 60 % ratio (w/w) of water to plaster. This yields an average porosity of the block of φ = 50 %. A special care during the preparation has been put on the saturation of the medium with water in chemical equilibrium with the gypsum. Before the injection, the chip is placed in vacuum within a beaker of saturated water in order to remove any air bubbles which might have remained in the pore space. Sealing is finally ensured with a silicone joint on the sides of the system. After sealing, we connected the syringe pump (Harvard Apparatus PHD2000) to the chip to inject pure water. We recorded the experiment with a UI 1550LE-C-HQ CCD camera IOS, Germany), acquiring photographic images of the system every 100 sec. In order to ensure homogeneous intensity of light over the system we used a circular fluorescent illuminator. Since the hydraulic resistance of the gypsum is much larger than the hydraulic resistance of the slot between the top cover and the plaster layer, we can safely assume that the whole flow is actually circulating on top of the soluble layer and does not penetrate in the porous medium. Therefore, the reactive surface area that the fluid sees is then only the top surface of the soluble layer. Since the roughness is relatively small because of the polishing step, we can assume that the reactive surface area can be approximated by the geometrical surface area.
The aperture of the fracture is several orders of magnitude smaller than its lateral dimensions, thus the mathematical description of the flow and transport in the fracture can be obtained by depth-averaged equations [7,8,9]. In particular, the aperture (h) increase due to the dissolution of the soluble layer is governed by the following equation:

x y t x y z t c x y z t dz h
is the flow-averaged concentration field, υ denotes the molar volume of gypsum and φ -its porosity. Next, is an effective reaction rate which accounts for the diffusive slowdown of a reaction as the aperture increases. In the above, D is the diffusion coefficient. and the Sherwood number Sh measures the dimensionless diffusive flux transporting the products of reaction away from the reactive surface. Two different parameters control the behavior of the system: the Péclet number Pe =vh0/D, which characterizes the relative magnitude of convective and diffusive transport, and the Damköhler number Da =k/v relating the surface reaction rate to the rate of convective transport. Note that in the experiment we can control both the flowrate (v) and the initial aperture (h0) and thus vary Pe and Da values independently. The dissolution of fractures proceeds nonuniformly -planar dissolution front is unstable and soon the flow becomes spontaneously localized in pronounced dissolution channels ("wormholes"). This is a manifestation of reactive infiltration instability, caused by a positive feedback between the reaction and fluid flow: a faster reaction locally leads to a permeability increase, which speeds up the local flow allowing the reactant to penetrate deeper inside the fracture. The initial wavelength of the instability is uniquely determined by the Peclet and Damköhler number [9]. Measuring the characteristics of the dissolution pattern allows then to trace back the characteristics of the dissolution and calculate the intrinsic reaction rate, k. In order to get the main wavelength of the dissolution front, the front is extracted from the pictures and a Fast Fourier Transform is applied. The tallest peak in the spectrum is then considered as the main wavelength of instability. Agreement between predicted wavelength and measured one is very good, usually below 10%.
The reaction rate can be estimated in yet another way -by measuring the time needed to erode the whole layer of soluble material in the fracture at a point close to the inlet (inlet dissolution time). Indeed, if we consider the erosion equation (3), we can calculate the inlet dissolution time as: where hmax is the maximum aperture (equal to the sum of the initial aperture and the thickness of the soluble layer).

Results and discussion
We have conducted 16 dissolution experiments, with different flow rates (in the range of 0.25 to 2 ml/hr) and different initial apertures (in the range of 70 to 210 μm). The values of the reaction rate estimated based on the experimental data are given in Fig. 2. Some of the experiments (no. 2, 3, 4 and 8) have been conducted several times, which allowed us to estimate the statistical error associated with a given method. The average value of the intrinsic reaction rate was estimated to be (for the estimate of k0 through inlet dissolution time). The lower value of the latter estimate is due to the assumption adopted in deriving Equation (5) that the incoming concentration is zero. In reality, there is a small but nonzero concentration of the calcium ions at the inlet, due to the upstream diffusion. This results in the decrease of the dissolution speed at the inlet, which is consistent with our results.
Currently in the literature, several values for the dissolution rate of gypsum coexist. Based on our experiments we can conclude that the actual intrinsic reaction rate of the dissolution of gypsum is very close to the one obtained by Colombani [5] and his reassessment of the data in the literature, and to the value obtained by Mbogoro et al. [6]. Other values reported in the literature [1][2][3]5] are almost two orders of magnitude higher. As argued by Colombani, such a significant overestimation of the reaction rate comes from the underestimation of the specific reactive area -he advocates using the BET surface area instead of geometric surface area for the normalization of the dissolution rates.
Summarizing, we have devised a novel method of estimating the reaction rate based on the initial wavelength of the reactive infiltration instability in a dissolving microfluidic chip. Importantly, the method can be relatively easily adapted to other reactions and minerals. Indeed, gypsum was used in these experiments because of the simplicity associated with its chemistry and moulding, but nothing prevents the use of other minerals which can be used as compacted powder and inserted in the indentation of the bottom plate. Two precautions have to be taken though: one should avoid solvents which react to polycarbonate (like acetone) and one should avoid the use of metal in the construction of the cell since it is likely to react through oxidation with the ionic solution. Importantly, the method is not limited to linear kinetic laws, since the instability wavelength can also be estimated in a nonlinear case [8]. Thanks to the above, we are confident that the method proposed here would provide a reliable and efficient alternative method for measuring the intrinsic reaction rate of minerals. We acknowledge the support of the National Science Centre (Poland) under research grant 2012/07/E/ST3/01734. The authors wish to thank Jean Colombani from the University of Lyon for the discussions on the reaction rate of gypsum and ROLLFIX ® for providing the tape used in the experiments.