## Abstract

We propose a new composite similarity variable, based on which a similarity solution is derived for reaction front propagation in fracture–matrix systems. The similarity solution neglects diffusion/dispersion within the fracture and assumes the existence of a sharp reaction front in the rock matrix. The reaction front location in the rock matrix is shown to follow a linear decrease with distance along the fracture. The reaction front propagation along the fracture is shown to scale like diffusion (i.e. as the square root of time). The similarity solution using the composite similarity variable appears to be applicable to a broad class of reactive transport problems involving mineral reactions in fracture–matrix systems. It also reproduces the solutions for non-reactive solute and heat transport when diffusion/dispersion/conduction are neglected in the fracture. We compared our similarity solution against numerical simulations for nonlinear reactive transport of an aqueous species with a mineral in the rock matrix. The similarity solutions agree very well with the numerical solutions, especially at later times when diffusion limitations are more pronounced.

This article is part of the themed issue ‘Energy and the subsurface’.

## 1. Introduction

Transport phenomena in fractured rock are relevant in the context of several energy-related subsurface activities including nuclear waste disposal, geothermal energy extraction and carbon sequestration. There is a significant body of literature on the transport of mass and energy in fractured rock, which highlights the role of advection in highly permeable fractures coupled to diffusion-controlled transport in the rock matrix [1–4]. Computer simulation of reactive transport and coupled processes in fractured rock has also advanced significantly in recent years [5–8]. Some of these models also incorporate the alteration of fracture aperture and hence permeability, resulting from dissolution and precipitation reactions [9–12]. An alternative class of alteration problems involves reaction with minerals in the rock matrix, which are accessed by matrix diffusion. In these problems, the fracture aperture/permeability remain largely unaltered while the mineral in the rock matrix is depleted, producing a reaction front (moving boundary) that advances into the rock matrix [4]. Similar problems also arise in the context of the weathering of fractured bedrock [13].

Our objective in this paper is to propose a simple similarity solution for reaction front propagation in fracture–matrix systems, applicable to situations where mineral reactions occur in the rock matrix. Our solution is valid for coupling of a pure advection equation in the fracture with a nonlinear diffusion–reaction equation in the rock matrix. In the limit of relatively fast reaction kinetics (i.e. diffusion-limited reactions), the nonlinearity in the reaction terms may be replaced by the nonlinearity associated with a moving boundary as in a Stefan-type problem, solutions to which are well established in the literature [14,15]. Our solution relies on the occurrence of such moving boundary reaction fronts. Our similarity solution predicts the evolution of these reaction fronts in two dimensions (perpendicular and parallel to the fracture). We present the coupled transport equations in §2, and the similarity solution in §3. We compare our similarity solutions to numerical simulations in §4.

## 2. Transport equations for reaction front propagation

The transport equations describing the propagation of a reaction front in a fracture–matrix system are summarized in this section. A schematic representation of a fracture–matrix system is shown in figure 1. The problem under consideration involves continuous flushing of a reacting species (whose molar concentration is denoted by *C*) through a fracture, which diffuses into and reacts with a solid mineral (whose molar density, i.e. moles-mineral per unit volume, is denoted by *W*) in the rock matrix.

Laplace-transform solutions to non-reactive solute and heat transport equations in fracture–matrix systems have been obtained including the influence of diffusion/dispersion [1–4], and rely on numerical inversion of Laplace transforms. When diffusion/dispersion in the fracture is neglected, the Laplace transforms can be inverted analytically. The similarity solution to be developed in §3 below also relies on the neglect of diffusion/dispersion in the fracture, and leads to simple analytical expressions. A similar simplification was employed in [4] to interpret results from numerical simulations of reaction front propagation. Also, in [1–4] and other previous analyses of transport in fracture–matrix systems, concentration variations across the (small) fracture aperture are neglected, and an aperture-averaged fracture transport equation is employed. The aperture-integrated advective transport equation in the fracture with diffusion/dispersion neglected is
2.1where *C*_{f}(*x*,*t*) denotes the aperture-averaged concentration of the reacting species in the fracture, *b* is the fracture aperture (figure 1), *q* is the aperture-integrated flux per unit width in the fracture (with dimensions of L^{2}/T), *ϕ* and *D*_{e}, respectively, denote the porosity and effective diffusivity in the rock matrix, and *C*_{m}(*x*,*z*,*t*) denotes the concentration of the reacting species in the rock matrix. Note that in (2.1), the porosity within the fracture is assumed to be one, and the fracture–matrix interface flux on the right-hand side is assumed to be the same on both sides (*z*′=±*b*/2). As indicated in figure 1, *x* denotes a coordinate along the fracture, *z*′ a coordinate perpendicular to the fracture, with (*z*′=±*b*/2) denoting the fracture–matrix interfaces, and *t* denotes time. The term on the right-hand side of (2.1) arises from the condition of continuity of normal flux at the fracture–matrix interfaces, thus arising naturally in the boundary flux terms resulting from aperture averaging.

The reactive transport equations in the rock matrix are presented below. As in [1–4] and other previous works on matrix diffusion and fracture–matrix exchange of mass or energy, diffusion/conduction parallel to the fracture (i.e. in the *x* direction) is neglected in the rock matrix. The diffusion–reaction equation for the reacting species is
2.2The reaction term on the right-hand side is a bimolecular reaction between the reacting species and the mineral with rate parameter *k* (dimensions of mole^{−1} L^{3} T^{−1}), *W* is the mineral molar density (moles of mineral per unit porous medium volume), and *m* is a stoichiometric coefficient (*m* moles of *C* react with *n* moles of *W*). The corresponding mineral mass balance equation is
2.3The right-hand side terms in (2.2) and (2.3) constitute a reasonable representation of the mineral reaction rate, which goes to zero in the absence of either the reacting species or the mineral locally (e.g. [16,17]). Under diffusion-controlled conditions, the exact form of the reaction kinetics does not affect the overall behaviour of the reaction front, as will become evident below. In general, the initial conditions for (2.1)–(2.3) may be defined as follows:
2.4In (2.4), *W*_{i} is the initial mineral molar density in the unaltered rock matrix.

The boundary conditions employed in the solution of (2.1)–(2.3) are
2.5*a*and
2.5*b*where *C*_{0} is the constant concentration maintained at the inflow boundary. The condition (2.5a) imposes continuity of the concentration across the fracture–matrix interface. An additional boundary condition of *C*_{m}→0 at may be imposed. However, diffusivities are typically very small in the rock matrix and a relatively sharp (moving boundary) reaction front forms, at which the concentration of the reacting species goes to zero and the mineral volume fraction tends to its initial value *W*_{i}. For the case of simple one-dimensional diffusion and reaction between *C* and *W*, Lichtner [14] presented one of the earliest moving boundary analyses. The location of the moving boundary needs to be determined as part of the solution, employing the Stefan condition [14]. A reformulation of (2.2) incorporating the boundary condition at the reaction front is described further below. It is readily evident by symmetry that the solution for *C*_{m} is identical in the rock matrix on either side of the fracture. Thus, it suffices to determine the solution in the matrix on one side of the fracture. For convenience of notation, we introduce a shifted coordinate *z*=*z*′−*b*/2 (figure 1) on the positive side (*z*′>0), and solve for *C*_{m}(*x*,*z*,*t*) below. Thus (2.1) may be rewritten for *z*′>0 as
2.6In general, *ϕ* and *D*_{e} vary across the reaction front, because the mineral is depleted before the front (*W*=0), and unreacted ahead of the front (*W*=*W*_{i}). However, since the reacting species undergoes diffusion only up to the front (at which it is fully consumed by the reaction), only the values of *ϕ* and *D*_{e} before the front are relevant to the analysis. Thus, we hereinafter use *ϕ* and *D*_{e} to denote the constant matrix properties in the zone where the mineral is depleted. As noted above, we invoke the presence of a sharp moving boundary reaction front *z*_{f}(*x*,*t*), at which *C* goes to zero, and ahead of which the mineral is unreacted, and replace (2.2) by
2.7Correspondingly, the mineral mass balance equation (2.3) may be replaced with
2.8The boundary conditions relevant to (2.6)–(2.8) include (2.5*a*,*b*) above, and the condition of zero concentration at the reaction front:
2.9To solve (2.6)–(2.8), as noted above, an additional equation for the moving boundary reaction front *z*_{f}(*x*,*t*) is needed, which is obtained from a mass balance at the reaction front [14,15] (i.e. the Stefan condition):
2.10We acknowledge that if diffusion along *x* were included in the formulation, a more general form of the Stefan condition would involve the diffusive flux of *C*_{m} normal to *z*_{f}(*x*,*t*). However, when *x*-diffusion is neglected, (2.10) is a consistent representation of the Stefan condition.

## 3. Similarity solution

In this section, we present a similarity solution to (2.6)–(2.8) with the boundary conditions (2.5*a*,*b*), (2.9), and the Stefan condition (2.10). We assume the following forms for the solutions to *C*_{f} and *C*_{m}, invoking a composite similarity variable that involves both space variables *x* and *z*:
3.1*a*and
3.1*b*where *α*, *β* and the function *f* are to be determined from (2.6) to (2.8). Note that (3.1*a*,*b*) automatically satisfy the interface boundary condition (2.5a). First, we note that if a function *f* can be found that solves (2.6)–(2.8), (2.9) requires that
3.2*a*If *f* is a monotonic function with a single root, it is single-valued, thus implying that
3.2*b*Next, we use (3.1*a*,*b*) in (2.6) to obtain
3.3Note that the proportionality of *α* to *β* suggests that *β* is merely a scaling factor and will not be present in the final solution. Using (3.1b) in (2.7) along with (3.3), and the boundary conditions (2.5a), (2.9), the solution for *f* is obtained as
3.4Returning to (3.2b), the behaviour of the reaction front *z*_{f}(*x*,*t*) is obtained from
3.5*a*which simplifies to
3.5*b*The form of (3.5b) already suggests that for *x*/*vt*≪1, *z*_{f}(*x*,*t*) decreases linearly with *x*, as observed in the numerical simulations of Steefel & Lichtner [4]. The location at which *z*_{f}(*x*,*t*) goes to zero defines the position of the reaction front in the fracture, *x*_{f}(*t*), for which a relationship is obtained from (3.5b)
3.5*c*The final element needed to complete the solution is *z*_{f}(0,*t*), which is obtained from the solution to a one-dimensional diffusion–reaction problem with a moving boundary (i.e. equations (2.7) and (2.8) at *x*=0). As in typical Stefan problems, a form is assumed for the reaction front location in (2.10) to yield
3.6*a*Considering that *C*_{0}/*W*_{i}≪1 due to the large molar density of mineral in the solid phase compared with maximum aqueous phase concentrations, (3.6a) suggests that *D*_{f0}/*D*_{e}≪1. Using approximations to the error function and exponential for *D*_{f0}/*D*_{e}≪1, (3.6a) simplifies to
3.6*b*Because *D*_{f0}/*D*_{e}≪1, and the diffusive transport in the rock matrix is much slower than advection in the fracture, it is reasonable to expect that the propagation of the reaction front in the fracture is much slower than that of a non-reactive advection front, i.e. *x*_{f}(*t*)≪*vt*. Using this approximation in (3.5c) yields
3.7It is thus clear that the reaction front in the fracture propagates relatively slowly as and that the approximation *x*_{f}(*t*)≪*vt* is well justified. With the advancement of the reaction front in the fracture and matrix given by (3.7), it should be noted that the above similarity solutions are only valid for *x*≤*x*_{f}(*t*). For *x*>*x*_{f}(*t*), the concentrations and mineral density remain equal to their initial values (i.e. (2.4)).

The condition that *x*_{f}(*t*)≪*vt* and the recognition that the similarity solutions are only valid for *x*≤*x*_{f}(*t*) also justify dropping *x*/*vt* in (3.5b), and *x*/*v* in (3.1), (3.2), (3.4) and (3.5a). It can be shown readily that these approximations are equivalent to dropping the transient term in the fracture transport equation (2.1) and using the resulting quasi-steady fracture transport equation to derive the similarity solution. A quasi-steady transport equation is also justifiable considering the slow rate of mineral depletion and using the same arguments employed in going from (3.6a) to (3.6b). Incorporating these simplifications, the final form of the similarity solution is summarized below, with *x*_{f}(*t*) given by (3.7) and :
3.8*a*
3.8*b*
3.8*c*The solution for *W* is given by (2.8) for *x*≤*x*_{f}(*t*), *W*=*W*_{i} for *x*>*x*_{f}(*t*); and *C*_{f} and *C*_{m} are zero outside the ranges specified in (3.8b) and (3.8c), respectively.

It is also interesting to examine the behaviour of the spatially integrated mineral depletion rate in the rock matrix, which is defined as
3.9Note that (3.9) assumes an infinitely long fracture and infinitely wide rock matrix. Because the reaction occurs largely at the reaction front, *R*(*t*) can also be expressed as the rate at which mineral is depleted by reaction front propagation
3.10Evaluating the time derivative in (3.10) using the Leibniz rule and (3.7) for *x*_{f}(*t*), it is readily shown that
3.11The spatially integrated reaction rate *R*(*t*) is time invariant because in a semi-infinite system (in both *x* and *z*), the entire influx of reactive species mass is consumed within the domain at all times (after a short initial time as discussed below in the next section). In fact, the constant value of *R*(*t*) is simply equal to the reactive power of the inflow solution, where the factor *n*/*m* denotes the number of moles of *W* required to react with one mole of *C*. For a fracture with finite length (*L*), it is readily shown that *R*(*t*) is given by (3.11) before the reaction front reaches the end of the fracture, i.e. *t*≤*t*_{L}, where *t*_{L} is the time taken by the reaction front to reach the end of the fracture, given by
3.12For *t*>*t*_{L}, *R*(*t*) decreases as the inverse square root of time, according to
3.13

## 4. Comparison with numerical simulations

In this section, we compare the analytical solutions presented in §3 with numerical simulations of the transport equations (2.1)–(2.5*a*,*b*). The discretization of the fracture transport equation (2.1) employed an Euler-backward discretization of the time-derivative term, upwind differencing for the advection term, and a staggered temporal coupling with the matrix. The source/sink term on the right-hand side of (2.1) was approximated based on the solution to the matrix transport equations in the previous time step. After solving for *C*_{f} along the fracture, it is used as a boundary condition (2.5a) to solve the matrix transport equation in each time step. An Euler-backward discretization in time and the standard second-order approximation to the diffusion term were used for the matrix transport equation. A Newton–Raphson iteration was employed to solve the resulting nonlinear system of equations. The grid spacing was varied in different simulations, depending on the flow rate. A grid spacing (d*x*) of 1 to 2 cm in the fracture and 0.5 to 2 mm in the matrix (d*z*) were employed, and time steps ranged from 20 to 100 *b*×d*x*/*q*. Simulations were performed for three sets of parameter values, as listed in table 1.

Figure 2 shows the numerical simulation results and analytical solutions for *C*_{m} and *W* at different times and *x*-locations for case 1. The analytical solution (3.8c) for *C*_{m} agrees very well with the corresponding simulation results. The simulation results for *W* are more spread out than the sharp-front analytical solutions (2.8). Better agreement between the analytical and numerical solution for concentrations and mineral molar density would be expected if the reaction rate *k* were increased, because the sharp-front approximation (equations (2.7)–(2.9)) is more accurate for higher *k* values. At any given value of *x*, the agreement between the simulation results and analytical solution improves with time, because the diffusion limitation assumed in deriving the analytical solution is more valid at later times. Figure 3 shows that the analytical solution (3.8b) for *C*_{f} along the fracture agrees very well with the numerical simulation results, except very near the advancing concentration front in the fracture (*x*_{f}(*t*)). The analytical solutions based on using *t*−*x*/*v* as in (3.4) or *t* as in (3.8*b*,*c*) are essentially indistinguishable, because *x*_{f}(*t*)≪*vt*. Figure 4 shows the reaction front location in the matrix (*z*_{f}(*x*,*t*)) for all three cases simulated. The agreement between the analytical expression (3.8a) and the numerical solution is excellent, thus confirming the linear decrease of *z*_{f}(*x*,*t*) and equation (3.7) for *x*_{f}(*t*). Finally, we compare the analytical expression (3.11) for the spatially integrated reaction rate *R*(*t*) with the results of the numerical simulations. In the numerical simulations, *R*(*t*) was calculated by integrating the reaction term on the right-hand side of (2.2) and (2.3) over the entire computational domain. Figure 5 compares (3.11) with the spatially integrated rate obtained from numerical simulations. In the numerical simulations, there is an early time regime where *R*(*t*) first increases, and then decreases to approach the value predicted by (3.11) over a duration of about 100 days. This early time behaviour corresponds to the duration before a well-defined reaction front forms and is further described by Arshadi & Rajaram [18].

## 5. Concluding remarks

We have proposed a similarity solution for reaction front propagation in fracture–matrix systems. Although we developed our solution by considering a bimolecular reaction between a reacting species and a mineral, the exact form of the reaction kinetics does not influence the behaviour, and the analytical solutions ((3.1*a*,*b*), with (3.4); or (3.8*a*–*c*)) are expected to be valid for a broad class of transport problems involving fracture–matrix interactions. In fact, our similarity approach reproduces the simplified forms of the Laplace-transform solutions [1–3,19] for non-reactive solute transport or heat transport, when diffusion/dispersion in the fracture are neglected. This may be readily verified by setting in (3.4) or (3.8*b*,*c*). Our similarity solution also reproduces the behaviour of the reaction front observed in numerical simulations of a multi-component reactive transport problem with mineral interactions [4]. In particular, the quantity *q*/(*ϕD*_{e}) is a type of Peclet number that serves as a scaling factor relating reaction front propagation along the fracture to that in the matrix (equation (3.7)). As noted in [4], the overall nature of reaction front propagation in the fracture is thus diffusive in character. Our composite similarity variable and similarity solution naturally capture this scaling. We believe that the similarity variable employed here is a general tool that can facilitate solutions to multi-component transport and coupled alteration processes in fracture–matrix systems. The main limitation of the similarity solution is in the use of the simplified fracture transport equation (2.6), which is, however, a useful approximation in many practically relevant situations. Similarity solutions of the type proposed here will also be useful for verifying numerical simulators of reactive transport and coupled processes in fractured rock.

## Competing interests

We declare we have no competing interests.

## Funding

This research was funded by Institute of Geophysics, Planetary Physics and Signatures (IGPPS) at Los Alamos National Laboratory (subcontract no. 230859) and the National Science Foundation (EAR0739803, EAR1239281).

## Acknowledgements

We thank two anonymous reviewers for their comments and suggestions.

## Footnotes

One contribution of 12 to a theme issue ‘Energy and the subsurface’.

- Accepted July 6, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.