Nutrient dynamics and N-anomaly at the SEATS station

(1)

Numerical simulation of the evolution of aquifer porosity and

species concentrations during reactive transport

$

Jui-Sheng Chen

a

, Chen-Wuing Liu

b,

*

a

Department of Environmental Engineering and Sanitation, Foo-Yin Institute of Technology, Kaohsiung, Taiwan, ROC b

Department of Agricultural Engineering, National Taiwan University, Taipei, Taiwan, ROC Received1 August 2000; receivedin revisedform 5 June 2001; accepted10 June 2001

Abstract

While flowing through a porous medium, a reactive fluid dissolves minerals thereby increasing its porosity and ultimately the permeability. The reactive fluidflows preferentially into highly permeable zones, which are therefore dissolved most rapidly, producing a further preferential permeability enhancement. Thus, the reaction front may be unstable. However, other factors, such as diffusion, suppress the instability of a reaction front. This study presents a numerical model to evaluate the interactions between mechanisms that determine the shape of a reactive front. That is, a methodis developedto solve a set of nonlinear equations coupledwith fluidflow, species transport, androck–fluid reactions and includes the effects of grain dissolution and the alteration of porosity and permeability due to mineral– fluid reactions. The numerical model enables us to evaluate how a dissolution reaction affects the porosity structure and fluid pressure variation, from which local Darcy flux can then be evaluated. In addition, the model is used to examine how upstream pressure gradient affects the morphological instability of the species concentrations and the aquifer porosity. Simulation results indicate that, although stable for small upstream pressure gradients, the growth of a planar front becomes unstable for large upstream ones. Moreover, the diffusive, advective and resultant species fluxes of both these mechanisms are computedandpresentedto further elucidate the behavior of the morphological instability for a planar concentration andporosity front that results from the interactions between diffusion andadvection. r 2002 Elsevier Science Ltd. All rights reserved.

Keywords: Permeability; Evolution; Dissolution; Reactive transport

1. Introduction

In subsurface formations, modeling coupled reactions andﬂow involve critical geochemical functions in terms of interpreting phenomena such as weathering,

diagen-esis andore deposition (Lichtner, 1996). The water–

mineral interaction is of relevant interest in several branches of geochemistry that includes subsurface ﬂow andchemical processes. The change in aquifer porosity owing to the interaction of a penetrating ﬂuidwith a

porous solidhas receivedconsiderable attention. That is, due to complete rock dissolution, these phenomena can produce dissolution fingers and create flow channels (Liu et al., 1997). When flowing through a porous medium, a reactive fluid dissolves the mineral, thereby increasing porosity andultimately permeability. Re-gardless of whether hydraulic or thermal driving forces are maintained, the flow pattern evolves as the perme-ability distribution changes. The reactive fluid flows preferentially into zones of high permeability, which are subsequently dissolved rapidly, ultimately producing a further preferential permeability enhancement. Thus, the reaction front may be unstable. However, other factors such as diffusion consistently compete with the flow focusing mechanism, which prevents the elongation

$

Code available from server at http://www.iamg.org/ CG-Editor/index.htm.

*Corresponding author. Fax:+886-2-23639557.

E-mail address:lcw@gwater.agec.ntu.edu.tw (C.-W. Liu).

(2)

of the finger from proceeding indefinitely. The diffusion effect inhibits nonuniformities of the reaction front, which ultimately produces either a planar or nonplanar reaction front. The instability of reaction fronts has been studied extensively (Chadam et al., 1986). Through an

analytical approach, Sherwood(1987) developeda

simple model that demonstrates how the Damk+ohler

number andthe acidcapacity number affect the instability of a planar reaction front. However, this investigation did not incorporate the porosity variation, which produces variation in microscopic fluid velocity. Hinch andBhatt (1990) examinedthe stability of a moving reaction front using the standard of chemical reaction rate proportional to the product of concentra-tion of fluidandsoluble minerals. Furthermore, Ortoleva et al. investigatedthe stability of a moving reaction front (Chadam et al., 1986; Ortoleva et al., 1987a, b). They usedcomplex reaction equations and included solute diffusion as well as significant changes in the porosity. A series of studies, which considered layeredporous media (Xin et al., 1993), was performed on the shape stability of a moving reaction front via the perturbation method, in which a moving reactive front for viscosity changes (Chadam et al., 1991). According to their results, a planar dissolution front is unstable and small perturbation to the front may develop into long fingers. Although the linear stability analysis of reactive transport system equations demonstrates the morpho-logical instability of planar fronts, the nonlinear problem must be numerically simulatedto understand more thoroughly the dynamics of the reactive fingering

problem. Steefel andLasaga (1990) usedtwo-d

imen-sional simulations to investigate how coupling chemical reactions andfluidflow affect the space–time evolution of rock dissolution patterns. Lichtner (1992) developed a system of reactive transport equations basedon

continuum approach andappliedto study supergene

copper enrichment andmetasomatic zoning (Lichtner andBiino, 1992; Lichtner andBalashov, 1993). Chen and Ortoleva (1990) used a computer code to predict the temporal evolution of the reactive feedback problem. Recently, Liu et al. (1997) proposeda geochemical simulator that coupledreaction andtransport and predicts the formation of dissolution ﬁngers and worm-holes. Renardet al. (1998) conductedan experiment to investigate self-organization phenomenon andcom-paredthe experiment to a numerical model of water– rock interaction. Although there have been numerous researches that attempt to accurately predict the instability of reactive fronts, relatively few have attemptedto thoroughly investigate the competition of relevant mechanism on reaction front shape quantita-tively (Steefel andLasaga, 1994; Lichtner andSeth, 1996; Liu et al., 1997). By expanding upon previous studies, this work presents a two-dimensional numerical model for the reaction-induced feedback problem to

obtain the time evolution of species concentrations and aquifer porosity. Quantitative analysis is also performed to evaluate interactions between relevant mechanisms on reaction front shape.

2. Mathematical model

The proposedmodel analyzes how coupling mineral– ﬂuidreactions aﬀect the space–time evolution of aquifer

porosity andspecies concentrations. Chadam et al.

(1986) have provided a complete derivation of the nonlinear partial differential equations that model this phenomenon. The governing equations employedherein are described briefly below. The conservation equations for flow, mass transport andreactions in a porous medium are formulated using a continuum approach in which a representative elemental volume (REV) is assumedto be smaller than the length scale of the phenomena that is being monitored. In addition, the various properties within the REV are assumedto remain constant (Lichtner, 1992).

2.1. Porosity change due to dissolution reactions Although this mathematical model is simple, there are generic situations of a single solidcomponent in the porous medium and a single solute in the ﬁeld, respectively. The matrix consists of soluble grains that

have an average volume of L3_{andgrain number density}

of n: The product nL3_{is the soluble grain that occupies a}

volume fraction. It also contains insoluble grains that

occupy a volume fraction f_n: It follows that

nL3þ f þ f_n¼ 1; ð1Þ

where f is the aquifer porosity.

If G is the rate of grain-volume change due to reaction, then the following occurs

qf

qt ¼ nG: ð2Þ

An implicit expression for G is requiredto completely describe the phenomenon description.

Assume that G can be written as

G ¼GSðc ceqÞ; ð3Þ

where G denotes the reaction rate constant, S ¼ L2

represents the surface area of the soluble grains, c refers

to the concentration of the species in solution and ceqis

the equilibrium concentration. The dissolution reaction is assumedto follow ﬁrst-order kinetics. The reaction rate is also assumedto be proportional to the surface

area L2_{of the soluble grains.}

Basedon Eqs. (1)–(3), the porosity change that occurred due to the dissolution reactions takes the following form (detailed derivation is provided in the

(3)

appendix): qf qt ¼ Gn 1=3_ðf f fÞ 2=3_{ðc c} eqÞ; ð4Þ

where f_f ¼ 1 fnis the ﬁnal porosity after a complete

dissolution of soluble grains. The 2=3-power indicates that only surface reactions are considered.

2.2. Continuity equation and momentum conservation The continuity equation, which describes ﬂuid mass conservation in the system, combinedwith the momen-tum conservation, Darcy’s law, can be written as

r ðfkðfÞrpÞ ¼qf

qt; ð5Þ

where p denotes the ﬂuid pressure and kðfÞ represents the porosity dependent permeability. The Fair–Hatch relation (Bear, 1972, p. 134) is employedherein in a modiﬁed form to describe the relationship between

porosity andpermeability (Chadam et al., 1986;

Ortoleva et al., 1987b) kðfÞ ¼ 1 Jm_wy2 f3 ½ð1 f_fÞ2=3þ n1=3_ðf f fÞ 2=32;

where J is a packing factor (B5), m_w is the water

viscosity (B0.001–0.01 P) and y is a geometric factor

(B6 for spherical factor). The permeability relationship

may also adopt other forms. For instance, Steefel and Lasaga (1990) usedthe Kozeny–Carman equation to describe this dependence. However, limited experimental data are available to justify one form over another. 2.3. Conservation of solute mass

A diﬀerential equation expressing conservation of species in solution must account for changes in the concentration in the ﬂuidphase. That is, the movement of species into andout of the REV must be includedin the governing equation to conserve the solute mass. Species conservation, in units of moles per unit volume of porous medium per unit time, can be written as

r ½fDðfÞrc þ cfkðfÞrp þ rs

qf

qt ¼

qðfcÞ

qt ; ð6Þ

where DðfÞ is the porosity dependent diﬀusion

coeﬃ-cient and r_sis the density of soluble grains. A common

phenomenological relation for DðfÞ is (Bear, 1972; Lerman, 1979)

DðfÞ ¼ Difm 3₂omo5₂;

where Diis a constant on the species diﬀusion coeﬃcient

order in water. The dispersion eﬀects are neglected. The

ﬁrst, secondandthirdterms on the left-handside of

Eq. (3) represent the diﬀusion, advection and source terms due to the mineral dissolution, respectively.

It is often useful to nondimensionalize a diﬀeren-tial equation because the resulting equation is generally simpler in form andmore applicable. Therefore a dimensionless time, t; is introduced and deﬁned by

t ¼ eðGn1=3_c

eqÞt: ð7Þ

The space variable r ¼ ðx; yÞ; concentration cðx; y; tÞ and pressure pðx; y; tÞ are also convertedto the following dimensionless forms: %r ¼ rðkceqÞ1=2; ð8Þ g ¼ c ceq ; ð9Þ %p ¼pkðffÞ Dðf_fÞ: ð10Þ

Eqs. (4)–(6) can then be written as (dropping the bars

andwriting dðfÞ ¼ fDðfÞ=Dðf_fÞ; lðfÞ ¼ fkðfÞ=kðf_fÞÞ; eqf qt¼ ðff fÞ 2=3_{ðg 1Þ;} _ð11Þ eqðfgÞ qt ¼ r ðdrg þ lgrpÞ þ qf qt; ð12Þ r ðlrpÞ ¼qf qt: ð13Þ 3. Numerical methods

The numerical model is used to analyze the above reactive transport problem. Eﬀorts to develop numerical reactive transport approaches in recent years have focusedon how to couple the reaction andtransport term (Yeh andTripathi, 1991). To solve the coupledset of equations, several methods have been proposed (Rubin, 1983). The most rigorous approach attempts to solve the governing equations simultaneously, which is commonly referredto as the fully-coupledmethod (Steefel andLasaga, 1994). Alternatively, iteratively-coupledtechniques can be employedto calculate the reaction andtransport. The sequential iteration ap-proach (SIA) is a general iteratively-coupledmethod that solves coupledset of equations in which reactions andtransport are solvedsequentially (for a discussion, see Yeh andTripathi, 1991). Notably, the SIA method was adopted herein. Yeh and Tripathi (1989) discussed explicit andimplicit schemes of the SIA. Owing to its superior convergent rate, the latter is more desirable than the former one. Therefore, a numerical solution technique, which is basedon implicit SIA concept, was developed to solve the coupled set of nonlinear diﬀerential equations. Herein, a subscript ðkÞ is usedto designate quantities at the kth time step. The following steps describe the implicit SIA method that advances the

(4)

model from time interval Dt to a system of algebraic equations for the model variables at ðk þ 1Þth time level:

1. Eq. (11) is usedandsolvedto yieldf if g is provided temporarily. Through g; the temporal evolution of f can be obtainedby applying the methodof separating variable andintegration with respect to t over the time

interval, tkptptkþ1; thus obtaining

Z fkþ1 fk df ðf_f fÞ2=3¼ Z tkþ1 tk gn kþ1=2 1 e dt: ð14Þ

Evaluating these integrals andrearranging the terms yields fn_kþ1¼ f_f 1 3 ðgn kþ1=2 1Þ e Dt þ 1 3ðff fkÞ 2=3 3 ; ð15Þ

where Dt ¼ tkþ1 tk; the subscript k þ 1=2 denotes the

intermediate time level between k and k þ 1; the superscript n denotes nth iteration for solving f at time

level k þ 1: The estimates, gn

kþ1=2; are assumedto be

gn

kþ1=2¼ ðgn1kþ1þ gkÞ=2 if n > 1; or gnkþ1=2¼ gk if n ¼ 1:

Evaluating fn_kþ1 as fn_kþ1 approaches f_f is essential. If

the choice of Dt is inappropriate, then the computed

value can fall in a region where fn_kþ1> f_f; thus

rendering the computation physically nonsensical. To

ensure that the computedvalue of fn_kþ1 is in the

invariant interval ½f₀; f_f ; Dt shouldbe smaller

than. 1 3ðff fkÞ 2=3 e gn1 kþ1=2 1 :

As f_k approaches f_f; a smaller Dtis requiredto satisfy

the numerical stability criterion. Nevertheless, there is an increase in the computing time for small Dt: This can

be handled by simply allowing fn_kþ1 to be f_f as f_f

fn_kþ1 is smaller than a given tolerance error, which

allows much larger time steps to be set in calculating temporal changes in porosity.

2. Eq. (12) is also discretized using implicit, ﬁnite

diﬀerence method. If pn

kþ1is provided temporarily, then

the computedporosity, fn_kþ1; from step 1 can be

substituted directly into the diﬀerence equation to

generate a temporarily known concentration, gn

kþ1: The

estimates of pn

kþ1=2 are assumedto be pnkþ1=2¼ ð pn1kþ1þ

pkÞ=2 if n > 1; or pn_kþ1=2¼ pk: The porosity dependent

permeability anddiﬀusion coeﬃcient are always

eval-uatedat the half-time level using ln_kþ1=2¼ lðfn_kþ1=2Þ

and dn

kþ1=2¼ dðf n

kþ1=2Þ: The spatial derivatives of

porosity dependent permeability and diﬀusion coeﬃ-cient in Eq. (12) are also evaluatedat the half-time level using a similar calculation. However, numerical error may appear when the permeability gradients, ql=qx and ql=qy; which are nonlinearly porosity dependent,

are discretized. To overcome these numerical

diﬃculties, an improved method was developed. The gradient is calculated using a semi-analytical derivative rather than the spatial discretization of permeability. The possible error is only the numerical error from the discretization of the spatial porosity derivative. Then, the permeability gradient can be computed as follows: qlðfÞ qx n kþ1=2 ¼ qf qx qlðfÞ qf n kþ1=2 : ð16Þ

Similarly, the diﬀusion coeﬃcient gradient can be computedas qdðfÞ qx n kþ1=2 ¼ qf qx qdðfÞ qf n kþ1=2 : ð17Þ

By grouping andrearranging terms, the final form of the finite difference equation becomes

Agn i; j;kþ1þ Bgniþ1; j;kþ1þ Cgni1; j;kþ1 þ Dgn i; jþ1;kþ1þ Egni; j1;kþ1 ¼ F gi; j;kþ Ggiþ1; j;kþ Hgi1; j;k þ I g_{i; jþ1;k}þ Jg_{i; j1;k}þ K; ð18Þ

where the subscript ði; jÞ denotes grid center,

A ¼ea1 Dtþ e a2 2þ a4 Dy2þ a6 Dz2 a8a9 2 a10lni; j;kþ1 2 a12a13 2 a14lni; j;kþ1 2 ; B ¼ a3 4Dx a4 2Dx2 a7a9 4Dx; C ¼ a3 4Dx a4 2Dx2þ a7a9 4Dx; D ¼ a5 4Dy a6 2Dy2 a11a13 4Dy ; E ¼ a5 4Dy a6 2Dy2þ a11a13 4Dy ; F ¼ea1 Dt e a2 2 a4 Dx2 a6 Dy2þ a8a9 2 þ a10li; j;k 2 þa12a13 2 þ a14li; j;k 2 ; G ¼ a3 4Dxþ a4 2Dx2þ a7a9 4Dx; H ¼ a3 4Dxþ a4 2Dx2 a7a9 4Dx; I ¼ a5 4Dyþ a6 2Dy2þ a11a13 4Dy ; J ¼ a5 4Dyþ a6 2Dy2 a11a13 4Dy ; K ¼ b1;

(5)

(

0

,

Ly

)

_{No Flow}

(

Lx

,

Ly

)

No Flow

(

0

,-

Ly

)

Upstream Pressure Gradient Flow Direction

(

Lx

,-

Ly

)

Fig. 1. Numerical simulation domain. Inlet ﬂuid is pumped in from left wall and ﬂows out of right wall.

Tim e=0.0 3.20 6.40 9.60 12.8 0 16 .0 0 3.20 6.40 9.60 12 .8 0 16 .0 0 (a) Tim e=0.0 0.70 0.70 1.40 1.40 2.10 2.10 2.80 2.80 3.50 3.50 (b) Tim e=0.0 0.24 0.48 0.72 0.9 6 1.20 0.24 0.48 0.72 0. 9 6 1.20 (c)

Fig. 2. Temporal evolution of aquifer porosity ðf ¼ ðf₀þ f_fÞ=2Þ with (a) upstream pressure gradient=0.5, (b) upstream pressure gradient=2.0 and (c) upstream pressure gradient=5.0. Solid line: Cauchy boundary condition, dashed line: Dirichlet boundary condition.

(6)

where a1¼ fn_{i; j;kþ1}þ f_{i; j;k} 2 ; a2¼ b2¼ fn_{i; j;kþ1} f_{i; j;k} 2 ; a3¼ qd qx n i; j;kþ1=2 ; a4¼ a6¼ dn i; j;kþ1þ di; j;k 2 ; a5¼ qd qy n i; j;kþ1=2 ; a7¼ a11¼ ln_{i; j;kþ1}þ li; j;k 2 ; a8¼ ql qx n i; j;kþ1=2 ; a9¼ pn iþ1; j;kþ1=2 pni1; j;kþ1=2 2Dx ; a10¼ p;n iþ1; j;kþ1=22pnj;kþ1=2þ pni1; j;kþ1=2 Dx2 ; a12¼ ql qy n i; j;kþ1=2 ; a13¼ pn iþ1; j;kþ1=2 p n i1; j;kþ1=2 2Dy ; a14¼ pn i; jþ1;kþ1=2 2pni; j;kþ1=2þ pni; j1;kþ1=2 Dy2 ; b1¼ fn_{i; j;kþ1} fi; j;k Dt :

3. Similarly, Eq. (13) is discretized to evaluate ﬂuid

pressure, pn

i;j;kþ1; if f n

i;j;kþ1that was obtainedpreviously

Time=0.0 3.20 6.40 9.60 12 .80 16 .00 3.20 6.40 9.60 12 .80 16 .00 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (a) Time=0.0 0.70 0.70 1.40 1.40 2.10 2.10 2.80 2.80 3.50 3.50 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (b) Time=0.0 0.24 0.24 0.48 0.48 0.72 0.72 0.96 0.96 1.20 1.20 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (c)

Fig. 3. Temporal evolution of species concentration ðg ¼ 0:5Þ with (a) upstream pressure gradient=0.5, (b) upstream pressure gradient=2.0 and (c) upstream pressure gradient=5.0. Solid line: Cauchy boundary condition, dashed line: Dirichlet boundary condition.

(7)

is substitutedinto finite difference Eq. (13). The final form after grouping andrearranging can be expressed again in the following form:

Apn_{i; j;k}þ Bpn iþ1; j;kþ Cp n i1; j;kþ Dp n i; jþ1;kþ Epnij; j1;k¼ F; ð19Þ A ¼ 2c2 Dx2 2c4 Dy2; B ¼ c1 2Dxþ c2 Dx2; C ¼ c1 2Dxþ c2 Dx2; D ¼ c3 2Dyþ c4 Dy2; E ¼ c3 2Dyþ c4 Dy2; F ¼ c5:

After implementing steps 1–3, the computedsolution fn_{i; j;kþ1}; gn

i; j;kþ1and pni; j;kþ1may be inaccurate andcan be

usedas the new guess value of the next solving iteration.

4. Repeat steps 1–3 successively until the solution has converged. The criterion set for convergence is

ðlnþ1_{i; j;kþ1} ln_{i; j;kþ1}Þ=ln_{i; j;kþ1}

_maxpzl; ð20Þ

where l refers to f; g; or p; zl is a speciﬁedresidue

constant and, the subscript max denotes the maximum value over all gridcenters.

4. Simulation results

The developed two-dimensional numerical model is appliedto simulate the permeability change due to the coupledﬂow andreaction as well as to investigate the

0.02 0.02 0.03 0.03 0.04 0.04 0.04 0.04 0.10 0.10 0.13 0.13 0.12 0.14 0.14 0.14 0.47 0.47 0.58 0.58 0.53 0.53 0.53 0.53 2.07 2.07 2.01 2.01 2.66 2.66 2.85 2.85 4.55 4.55 3.61 3.61 3.64 3.64 3.82 3.82 4.74 4.74 0.02 0.08 0.34 1.33 4.98 5.13 0.0 0 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

Species transport by advection

0.02 0.02 0.03 0.03 0.03 0.03 0.03 0.03 0.08 0.08 0.13 0.13 0.14 0.14 0.14 0.14 0.32 0.32 0.62 0.62 0.58 0.58 0.50 0.50 0.81 0.81 0.09 0.09 0.10 0.10 0.15 0.15 0.07 0.07 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.02 0.05 0.15 0.43 0.20 0.00 0.0 0 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

Species transport by diffusion

0.01 0.04 0.33 1.69 4.48 4.74 0.01 0.04 0.33 1.69 4.48 4.74 0.01 0.05 0.37 1.92 3.61 0.01 0.05 0.37 1.92 3.61 0.01 0.02 0.06 2.56 3.64 0.01 0.02 0.06 2.56 3.64 0.00 0.01 0.03 2.69 3.82 0.00 0.01 0.03 2.69 3.82 0.00 0.03 0.18 0.91 4.78 5.13 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

Resultant species transport by advection and diffusion

(a)

Fig. 4. Quantitative analysis of competition between diﬀusion and advection with (a) upstream pressure gradient=2.0 at time=0.7, (b) upstream pressure gradient=2.0 at time=2.1 and (c) upstream pressure gradient=2.0 at time=3.5.

(8)

unstable growth of the dissolution front. A hypothetical situation is deﬁned for simulation and illustration purposes. The problem considered here is a geometri-cally simple rectangular system, which was designed to simulate isothermal ﬂow andreaction in a carbonated cementedsandstone in which a percentage of the rock is reactive (e.g., a carbonate cement), whereas the remain-der is inert (e.g., a quartz at low temperature). Fig. 1 illustrates the simulation domain and the boundary conditions used herein.

It is assumedthat water is imposedin the positive x direction. Furthermore, a constant pressure gradient,

ðp0

fÞ; is applied, which induces a constant velocity ðvfÞ

ﬂow to enter the sandstone. Therefore, the boundary

condition of pressure at x ¼ 0 is prescribedas qp qx¼ p 0 f ¼ vf kf ðx ¼ 0Þ: ð21aÞ

The imposedﬂow is undersaturatedwith respect to the reactive mineral phase andthe value of the inlet concentration is assumedas constant concentration

strength g₀: Chadam et al. (1986) and Steefel and

Lasaga (1990) use a Dirichlet (the ﬁrst type) boundary for concentration at x ¼ 0 as

g ¼ g₀ ðx ¼ 0Þ: ð21bÞ

However, Van Genuchten andAlves (1982) andChen (1987) have suggestedthat as injectedﬂuids enter the 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.03 0.03 0.03 0.03 0.05 0.05 0.10 0.10 0.13 0.13 0.09 0.09 0.19 0.19 0.38 0.38 0.43 0.43 0.22 0.22 0.55 0.55 1.45 1.45 0.46 0.46 0.48 0.48 1.15 1.15 1.08 1.08 0.57 0.57 0.84 0.84 1.82 1.82 0.72 0.72 1.30 1.30 2.40 2.40 1.28 1.28 1.94 1.94 2.01 2.01 2.02 2.02 3.42 3.42 3.50 3.50 2.77 2.77 5.18 5.18 4.06 4.06 4.66 4.66 0.00 0.00 0.00 0.00 0.01 0.02 0.06 0.15 0.32 0.59 0.95 1.45 3.00 5.58 4.91 0.0 0 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.03 0.03 0.02 0.02 0.05 0.05 0.09 0.09 0.14 0.14 0.05 0.05 0.13 0.13 0.30 0.30 0.49 0.49 0.11 0.11 0.29 0.29 0.67 0.67 0.23 0.23 0.17 0.17 0.43 0.43 0.51 0.51 0.00 0.00 0.24 0.24 0.52 0.52 0.08 0.08 0.29 0.29 0.63 0.63 0.01 0.01 0.39 0.39 0.26 0.26 0.00 0.00 0.45 0.45 0.02 0.02 0.00 0.00 0.02 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.03 0.05 0.07 0.09 0.11 0.14 0.50 0.05 0.00 0.0 0 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.03 0.03 0.06 0.06 0.05 0.05 0.06 0.06 0.15 0.15 0.32 0.32 0.35 0.35 0.18 0.18 0.51 0.51 1.44 1.44 0.33 0.33 0.42 0.42 1.11 1.11 1.06 1.06 0.57 0.57 0.77 0.77 1.79 1.79 0.66 0.66 1.22 1.22 2.34 2.34 1.27 1.27 1.81 1.81 1.78 1.78 2.02 2.02 3.03 3.03 3.48 3.48 2.77 2.77 5.16 5.16 4.05 4.05 4.66 4.66 0.00 0.00 0.00 0.00 0.00 0.01 0.03 0.10 0.25 0.50 0.84 1.31 2.50 5.53 4.91 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

(b)

(9)

porous medium, a Cauchy (third type) boundary is more reasonable than a Dirichlet boundary. The Dirichlet boundary may produce a considerably diﬀerent result than the Cauchy boundary. That is, as injection begins, the Dirichlet boundary, which neglects the concentra-tion gradient across the interface, increases the local concentration level instantaneously. This condition causes more solute mass to enter the medium. Thus, the Cauchy boundary is also adopted in this study and formulatedas

vfg0¼ vfg DðfÞ

qg

qx ðx ¼ 0Þ: ð21b

0_Þ

Since the ﬂuiddensity is assumedto be constant, only the net pressure decrease across the region of interest is

important, therefore, the outlet of the system is ﬁxed. It is appropriate to apply the conditions such that

p ¼ pRðyÞ ¼ 0;

qg

qx¼ 0 ðx ¼ LxÞ: ð22Þ

Both sides have no ﬂow boundary condition such that qp

qy¼ 0;

qg

qy¼ 0 ð y ¼ Ly; LyÞ: ð23Þ

Initially, a small, local nonuniformity was introduced, which hada higher porosity near the left-endside of the system as the porosity distribution at time=0.0 as illustratedin Fig. 2a–c. Then, the initial concentration 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.03 0.03 0.04 0.04 0.03 0.03 0.06 0.06 0.11 0.11 0.18 0.18 0.08 0.08 0.16 0.16 0.36 0.36 0.74 0.74 0.17 0.17 0.38 0.38 0.86 0.86 0.88 0.88 0.33 0.33 0.70 0.70 1.48 1.48 0.30 0.30 0.56 0.56 1.09 1.09 2.05 2.05 0.35 0.35 0.85 0.85 1.55 1.55 0.93 0.93 0.52 0.52 1.21 1.21 2.07 2.07 0.77 0.77 1.66 1.66 2.43 2.43 1.36 1.36 2.30 2.30 2.19 2.19 2.09 2.09 3.94 3.94 3.58 3.58 5.08 5.08 3.98 3.98 4.57 4.57 0.00 0.00 0.00 0.00 0.01 0.02 0.06 0.13 0.25 0.43 0.67 0.97 1.35 1.87 3.80 5.55 4.81 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.03 0.03 0.05 0.05 0.02 0.02 0.04 0.04 0.10 0.10 0.20 0.20 0.04 0.04 0.11 0.11 0.26 0.26 0.60 0.60 0.08 0.08 0.20 0.20 0.43 0.43 0.61 0.61 0.12 0.12 0.29 0.29 0.52 0.52 0.06 0.06 0.16 0.16 0.35 0.35 0.59 0.59 0.00 0.00 0.20 0.20 0.40 0.40 0.33 0.33 0.00 0.00 0.23 0.23 0.44 0.44 0.06 0.06 0.27 0.27 0.53 0.53 0.01 0.01 0.37 0.37 0.18 0.18 0.00 0.00 0.28 0.28 0.01 0.01 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.02 0.04 0.06 0.07 0.08 0.09 0.10 0.15 0.42 0.02 0.00 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.04 0.04 0.07 0.07 0.09 0.09 0.05 0.05 0.13 0.13 0.30 0.30 0.67 0.67 0.14 0.14 0.33 0.33 0.82 0.82 0.97 0.97 0.28 0.28 0.66 0.66 1.46 1.46 0.28 0.28 0.50 0.50 1.06 1.06 2.09 2.09 0.35 0.35 0.79 0.79 1.52 1.52 0.86 0.86 0.52 0.52 1.14 1.14 2.04 2.04 0.72 0.72 1.58 1.58 2.35 2.35 1.35 1.35 2.16 2.16 2.03 2.03 2.09 2.09 3.68 3.68 3.57 3.57 5.07 5.07 3.98 3.98 4.57 4.57 0.00 0.00 0.00 0.00 0.00 0.01 0.03 0.09 0.19 0.35 0.58 0.88 1.25 1.72 3.38 5.53 4.81 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

(c)

(10)

andporosity conditions are chosen as follows:

fðx; y; 0Þ ¼ f₀þ ðf_f f₀Þex_; _ð24Þ

gðx; y; 0Þ ¼ ð1 e5xÞð1 ex_Þ; _ð25Þ

where xðx; yÞ ¼ ðx4_{þ y}4_Þ=ðwL

yÞ4is the initial

perturba-tion parameter.

The initial condition represents a situation in which

the porosity is near its initial post-washout values f_f at

the corner of x ¼ y ¼ 0 andis at its initial value f₀

elsewhere. The initial ﬁnger is wLyand2wLyin length

andwidth, respectively. Also, investigation herein is how the upstream pressure gradient, aﬀects the development

of the reaction front. The input parameters usedto simulate upstream pressure gradient, which inﬂuence morphological instability of reaction front, are as

follows: initial porosity, f₀¼ 0:1; ﬁnal porosity, ff ¼

0:2; inlet concentration, g₀¼ 0; length, Lx¼ 18; width,

2Ly¼ 8; upstream pressure gradients, 0.5, 2 and 5 and

initial perturbation parameter, w ¼ 0:1: Furthermore, the ﬁne-grid(Dx ¼ 0:2 and Dy ¼ 0:2Þ spacing both in the

x and y directions are applied to avoid numerical

dispersion.

The choices of Cauchy or Dirichlet boundary condi-tions for inﬂow boundary concentracondi-tions, which may aﬀect the simulation results, are also examined. Fig. 2a–c

Time=0.0 3.20 6.40 9.60 12 .80 16 .00 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (a) Time=0.0 0.70 1.40 2.10 2.80 3.50 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (b) Time=0.0 0.24 0.48 0.72 0.96 1.20 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (c)

Fig. 5. Temporal evolution of aquifer porosity ðf ¼ ðf₀þ f_fÞ=2Þ with (a) upstream pressure gradient=0.5, (b) upstream pressure gradient=2.0 and (c) upstream pressure gradient=5.0.

(11)

andFig. 3a–c plot the temporal evolution of simulated aquifer porosity andspecies concentrations contours for upstream pressure gradient=0.5, 2 and 5, respectively. The solid and dashed lines denote the porosity contours

with ðf₀þ ffÞ=2 ¼ ð0:1 þ 0:2Þ=2 ¼ 0:15

anddimension-less concentration contours with g ¼ ðg₀þ g_fÞ=2 ¼

ð0:0 þ 1:0Þ=2 ¼ 0:5 for Cauchy andDirichlet bound

-aries, respectively. The Cauchy andDirichlet boundaries produce considerably different reaction front advance-ments. Notably, the reaction front advancement for the former boundaries is slower than that for the latter. This difference is only significant for a small upstream pressure gradient. Dirichlet boundary neglects the local diffusion at the interface, as well as forces more solute

mass to enter the porous medium and accelerates reaction front movement. If upstream pressure gradient is larger, then local diffusion becomes less important. The advection dominates the transport that flows into the formation; hence, the Dirichlet boundary condition is close to that of the Cauchy one. Thus, the Dirichlet boundary condition is only valid for the large upstream gradient condition, which is generally not appropriated to model reactive transport through geological forma-tion. The third type Cauchy boundary condition is adopted at the inflow boundary. Additionally, as the upstream pressure gradient increases (and also inlet velocity), the instability of reaction front (concentration and porosity) also increases for both boundary conditions.

Time=0.0 3.20 6.40 9.60 12 .80 16 .00 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (a) Time=0.0 0.70 1.40 2.10 2.80 3.50 0.0 0 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (b) Time=0.0 0.24 0.48 0.72 0.96 1.20 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00 (c)

Fig. 6. Temporal evolution of species concentration ðg ¼ 0:5Þ with (a) upstream pressure gradient=0.5, (b) upstream pressure gradient=2.0 and (c) upstream pressure gradient=5.0.

(12)

The development of reaction front resembles the results of Chadam et al. (1986), who performed linear stability analysis to demonstrate the morphological instability of a planar front. Notably, the development of the instability of the planar front can be explainedas follows: according to Darcy’s law, if a protrusion of the porosity level line in the reacting zone exists at certain

times, the ﬂow of the most reactive ﬂuidtends to be

focusedin the tip of the protrusion. That is, the permeability at the tip is greater than that in the neighboring regions. As additional reactive fluid arrives in the tip, it tends to advance more rapidly, which subsequently causes fingering. Conversely, diffusion owing to the concentration gradient causes the fluid

within tip to be less reactive andhence decelerates its advancement. The competition between these two mechanisms results either in the shape selection of reacting zone (self-organization), which in turn produces the decay of this protrusion, or in the temporal development of successively more elongated protrusion. While attempting to account for competition between diffusion and advection (the fluid focusing at tip) mechanisms, Fig. 4a–c displays the spatial distribution of the diffusive, advective and resultant species flux of both for times 0.7, 2.1 and3.5 with an upstream pressure gradient of 2.0. The arrows and their adjacent values denote the direction and magnitude of the species flux transport, while the solidline represents the concentration

0.02 0.02 0.02 0.02 0.03 0.03 0.04 0.04 0.05 0.05 0.06 0.06 0.08 0.08 0.13 0.13 0.17 0.17 0.18 0.18 0.25 0.25 0.32 0.32 0.60 0.60 0.62 0.62 0.63 0.63 1.31 1.31 1.43 1.43 2.02 2.02 2.29 2.29 2.77 2.77 3.77 3.77 4.76 4.76 4.14 4.14 3.43 3.43 3.42 3.42 4.45 4.45 4.68 4.68 4.30 4.30 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.01 0.01 0.02 0.02 0.03 0.03 0.04 0.04 0.04 0.04 0.05 0.05 0.06 0.06 0.12 0.12 0.17 0.17 0.18 0.18 0.21 0.21 0.19 0.19 0.45 0.45 0.72 0.72 0.66 0.66 0.85 0.85 0.56 0.56 0.56 0.56 0.04 0.04 0.06 0.06 0.02 0.02 0.11 0.11 0.02 0.02 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.04 0.04 0.07 0.07 0.06 0.06 0.02 0.02 0.04 0.04 0.18 0.18 0.46 0.46 0.32 0.32 0.03 0.03 0.76 0.76 0.88 0.88 1.60 1.60 2.24 2.24 2.71 2.71 3.75 3.75 4.65 4.65 4.13 4.13 3.43 3.43 3.42 3.42 4.45 4.45 4.68 4.68 4.30 4.30 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

(a)

Fig. 7. Quantitative analysis of competition between diﬀusion and advection with (a) upstream pressure gradient=2.0 at time=0.7, (b) upstream pressure gradient=2.0 at time=2.1 and (c) upstream pressure gradient=2.0 at time=3.5.

(13)

contour. Obviously, although diﬀusion tends to

decele-rate reaction front advancement, advection

con-sistently accelerates it. Diffusion produces a smaller species flux than advection does and, therefore, the resultant species flux tends to travel to the tip of the protrusion. To further illustrate the phenomenon of reaction front instability, the temporal evolution of aquifer porosity for two, initially small perturbations are investigated. That is, two small, local nonuniformities with higher porosity were introduced near the left-end side of the system, which the porosity distribution at time=0.0 in Fig. 5a–c demonstrates. The initial condi-tion represents a situacondi-tion in which the porosity is

proximal to its initial post-washout value f_f at the

corner of ðx ¼ 0; y ¼71Þ andis at its initial value f₀

elsewhere. Fig. 5a–c andFig. 6a–c display the temporal evolution of simulatedaquifer porosity andspecies concentration contours for upstream pressure gradi-ent=0.5, 2 and5, respectively. The result is similar to that of one finger, however two fingers have emerged for large upstream gradient. Fig. 7a–b display the spatial distribution of the diffusive, advective and resultant species flux of both the mechanisms at times

0.7, 2.1 and3.5 with an upstream gradient of 2.0.

Obviously, the resultant species flux tends to travel to the tips of two fingers. The above discussion on reaction front stability reveals that flow self-focusing that occurs from interactions between advection and diffusion can drive a front out of a planar state into fingering morphology. 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.03 0.03 0.05 0.05 0.06 0.06 0.03 0.03 0.05 0.05 0.10 0.10 0.19 0.19 0.26 0.26 0.07 0.07 0.13 0.13 0.29 0.29 0.67 0.67 0.59 0.59 0.16 0.16 0.27 0.27 0.64 0.64 1.52 1.52 0.40 0.40 0.36 0.36 0.53 0.53 1.15 1.15 1.13 1.13 0.84 0.84 1.19 1.19 1.24 1.24 2.09 2.09 2.00 2.00 3.07 3.07 4.14 4.14 4.20 4.20 3.19 3.19 4.09 4.09 4.61 4.61 4.37 4.37 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.03 0.03 0.05 0.05 0.07 0.07 0.02 0.02 0.03 0.03 0.08 0.08 0.18 0.18 0.33 0.33 0.03 0.03 0.08 0.08 0.19 0.19 0.46 0.46 0.72 0.72 0.07 0.07 0.13 0.13 0.32 0.32 0.66 0.66 0.02 0.02 0.21 0.21 0.19 0.19 0.42 0.42 0.45 0.45 0.00 0.00 0.84 0.84 0.36 0.36 0.60 0.60 0.05 0.05 0.02 0.02 0.22 0.22 0.09 0.09 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.02 0.02 0.01 0.01 0.01 0.01 0.03 0.03 0.07 0.07 0.12 0.12 0.12 0.12 0.04 0.04 0.09 0.09 0.24 0.24 0.61 0.61 0.69 0.69 0.09 0.09 0.22 0.22 0.59 0.59 1.53 1.53 0.39 0.39 0.15 0.15 0.40 0.40 1.07 1.07 0.94 0.94 0.84 0.84 0.77 0.77 0.88 0.88 1.80 1.80 1.96 1.96 3.05 3.05 3.92 3.92 4.11 4.11 3.18 3.18 4.09 4.09 4.61 4.61 4.37 4.37 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

(b)

(14)

5. Conclusion

A mathematical model was presented to elucidate the reaction front of the water–mineral interaction problem. A finite difference method was developed to solve a set of nonlinear equations, which are fully-coupledfluid flow, species transport androck/fluidreactions. Further-more, numerical algorithms were employedto improve the stability of the calculations. Numerical simulations were performedto examine the morphology of the reaction front. Results indicate that when the driving force exceeds a critical value, the planar dissolution

front becomes unstable andﬁnger-shapedfronts

emerge. Moreover, the unstable fronts become more

pronounced under higher upstream pressure gradients. Quantitative analysis of diﬀusive, advective and resul-tant species ﬂuctuations within these mechanisms was also performedto explain the shape selection of a reaction front.

Acknowledgements

The authors wouldlike to thank the National Science Council of Republic of China for ﬁnancially supporting this work under Contract Nos. NSC 87-2313-B-002-051 andNSC 88-2313-B-002-097. 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.03 0.03 0.04 0.04 0.02 0.02 0.03 0.03 0.07 0.07 0.13 0.13 0.19 0.19 0.05 0.05 0.09 0.09 0.21 0.21 0.47 0.47 0.77 0.77 0.12 0.12 0.21 0.21 0.49 0.49 1.13 1.13 0.31 0.31 0.25 0.25 0.41 0.41 0.90 0.90 1.72 1.72 0.54 0.54 0.64 0.64 0.81 0.81 1.50 1.50 1.32 1.32 1.20 1.20 1.89 1.89 2.64 2.64 3.08 3.08 2.72 2.72 3.77 3.77 4.75 4.75 4.54 4.54 3.50 3.50 4.36 4.36 4.20 4.20 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.00 0.00 0.01 0.01 0.02 0.02 0.03 0.03 0.05 0.05 0.01 0.01 0.03 0.03 0.06 0.06 0.13 0.13 0.22 0.22 0.03 0.03 0.06 0.06 0.15 0.15 0.36 0.36 0.89 0.89 0.05 0.05 0.11 0.11 0.26 0.26 0.56 0.56 0.11 0.11 0.12 0.12 0.16 0.16 0.36 0.36 0.67 0.67 0.00 0.00 0.52 0.52 0.24 0.24 0.46 0.46 0.20 0.20 0.00 0.00 0.18 0.18 0.60 0.60 0.43 0.43 0.01 0.01 0.00 0.00 0.03 0.03 0.01 0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.02 0.04 0.04 0.07 0.07 0.06 0.06 0.03 0.03 0.07 0.07 0.16 0.16 0.39 0.39 0.79 0.79 0.08 0.08 0.17 0.17 0.44 0.44 1.10 1.10 0.28 0.28 0.14 0.14 0.33 0.33 0.84 0.84 1.72 1.72 0.54 0.54 0.26 0.26 0.58 0.58 1.37 1.37 1.15 1.15 1.20 1.20 1.72 1.72 2.07 2.07 2.70 2.70 2.71 2.71 3.77 3.77 4.72 4.72 4.53 4.53 3.50 3.50 4.36 4.36 4.20 4.20 0.00 2.00 4.00 6.00 8.00 10.00 12.00 14.00 16.00 18.00 -4.00 -2.00 0.00 2.00 4.00

(c)

(15)

Appendix

In this appendix, Eq. (4) is derived via Eqs. (1)–(3). First, Eq. (1) can be rearrangedas

L ¼ 1 f fn n 1=3 ¼ ff f n 1=3 : ðA:1Þ

Inserting Eq. (A.1) into Eq. (3) yields

G ¼G ff f

n

2=3

ðc ceqÞ: ðA:2Þ

Eq. (4) can be obtainedby substituting Eq. (A.2) into Eq. (2) andexpressedas qf qt ¼ Gn 1=3_ðf f fÞðc ceqÞ: ðA:3Þ References

Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, Amsterdam, 764pp.

Chadam, J., Ortoleva, P., Sen, A., 1986. Reactive-inﬁltration instabilities. IMA Journal of AppliedMathematics 36, 207–220.

Chadam, J., Peirce, A., Ortoleva, P., 1991. Stability of reactive ﬂows in porous media: coupled porosity and viscosity changes. SIAM Journal of AppliedMathematics 31, 684–692.

Chen, C.S., 1987. Analytical solutions for radial dispersion with Cauchy boundary at injection well. Water Resources Research 23 (7), 1217–1224.

Chen, W., Ortoleva, P., 1990. Reaction front ﬁngering in carbonate-cementedsandstones. Earth Science Reviews 29, 183–198.

Hinch, E.J., Bhatt, B.S., 1990. Stability of an acidfront moving through rock. Journal of FluidMechanics 212, 279–288. Lerman, A., 1979. Geochemical Processes, Wiley, New York,

481pp.

Lichtner, P.C., 1992. Time–space continuum description of ﬂuid/rock interaction in permeable media. Water Resources Research 28, 13135–13155.

Lichtner, P.C., 1996. Continuum formulation of multicompo-nent-multiphase reactive transport. In: Lichtner, P.C., Steefel, C.I., Oelkers, E.H. (Eds.), Reactive Transport in Porous Media, Reviews in Mineralogy, Vol. 34. Miner-alogical Society of America, pp. 1–81.

Lichtner, P.C., Balashov, V.N., 1993. Metasomatic zoning: appearance of ghost zones in limit of pure advective mass transport. Geochimica et Cosmochimica Acta 57, 369–387. Lichtner, P.C., Biino, G.G., 1992. A ﬁrst principles approach to supergence enrichment of a prophyry copper protore. I.

Cu–Fe–S–H2O subsystem. Geochimica et Cosmochimica Acta 56, 3987–4013.

Lichtner, P.C., Seth, M.S., 1996, User’s Manual for MULTI-FLO: Part II MULTIFLO 1.0 andGEM 1.0. Multi-component-Multiphase Reactive Transport Model. CNWRA 96-010.

Liu, X., Ormond, A., Bartko, K., Li, Y., Ortoleva, P., 1997. A geochemical reaction-transport simulator for matrix aciding analysis anddesign. Journal of Petroleum Science & Engineering 17, 181–196.

Ortoleva, P., Chadam, J., Merino, E., Sen, A., 1987b. Self-organization in water–rock interaction systems II: the reactive-inﬁltration instability. American Journal of Science 287, 1008–1040.

Ortoleva, P., Merino, E., Moor, C., Chadam, J., 1987a. Geochemical self-organization I: feedback mechanisms andmodelling approach. American Journal of Science 287, 979–1007.

Renard, F., Gratier, J-P., Ortoleva, P., Brosse, E., Bazin, B., 1998. Self-organization during reactive ﬂuid ﬂow in a porous medium. Geophysical Research Letters 25 (3), 385–388.

Rubin, J., 1983. Transport of reacting solutes in porous media: relation between mathematical nature of problem formula-tion andchemical nature of reacformula-tions. Water Resources Research 19, 1231–1252.

Sherwood, 1987. Stability of a plane reaction front in a porous medium. Chemical Engineering Science 42(7), 1823–1829. Steefel, C.L., Lasaga, A.C., 1990, Evolution of dissolution

patterns: permeability change due to coupled ﬂow and reactions. In: Melchior, D.C., Bassett, R.L. (Eds.), Chemi-cal Modeling in Aqueous Systems II: American ChemiChemi-cal Society Symposium, Vol. 416. pp. 212–125.

Steefel, C.L., Lasaga, A.C., 1994. A coupledmodel for transport of multiple chemical species andkinetic precipita-tion/dissolution reactions with application to reactive ﬂow in single phase hydrothermal systems. American Journal of Science 294, 529–592.

Van Genuchten, M.Th, Alves, W.J., 1982. Analytical solutions of the one-dimensional convective–dispersive solute trans-port equation. UnitedStates Department of Agriculture Technical Bulletin, Vol. 1661. 149pp.

Xin, J., Peirce, A., Chadam, J., Ortoleva, P., 1993. Reactive ﬂows in layeredporous media II. The shape stability of the reaction interface. SIAM Journal of AppliedMathematics 53, 319–339.

Yeh, G.T., Tripathi, V.S., 1989. A critical evaluation of recent developments of hydrogeochemical transport models of reactive multichemical components. Water Resources Research 25 (1), 93–108.

Yeh, G.T., Tripathi, V.S., 1991. A model for simulating transport of reactive multispecies components: model development and demonstration. Water Resources Research 27 (11), 3075–3094.