Metabolic syndrome as a risk factor for coronary heart disease and stroke: An 11-year prospective cohort in Taiwan community

全文

(1)Journal of Contaminant Hydrology 72 (2004) 47 – 66 www.elsevier.com/locate/jconhyd. Interaction of reactive fronts during transport in a homogeneous porous medium with initial small non-uniformity Jui-Sheng Chen a, Chen-Wuing Liu b,* a. Department of Environmental Engineering and Sanitation, Foo-Yin University, Kaohsiung, Taiwan, ROC Department of Bioenvironmental Systems Engineering, National Taiwan University, 1, Roosevelt Road, Sec. 4, R.O.C. 10617, Taipei, Taiwan, ROC. b. Received 30 August 2002; received in revised form 13 October 2003; accepted 31 October 2003. Abstract A reactive fluid circulating within a porous medium can dissolve minerals with which it is out of equilibrium and modify the porosity and permeability. The positive feedback between fluid transport and mineral dissolution causes complex reaction front morphologies such as fingers or wormholes. This study presents a numerical model to investigate reaction front instability, temporal aquifer porosity, and species concentration evolution during reactive transport in a homogeneous porous medium with two small, initially local non-uniformities. Simulation results indicate that a stable planar front develops for a small upstream pressure gradient while the growth of two nonuniformities becomes unstable for a large upstream pressure gradient. Moreover, the unstable reaction front may be either double- or single-finger in shape. Reaction front shape selection depends on the spacing of the two local non-uniformities and the upstream pressure gradients. A behavior diagram is constructed to identify a planar, single- or double-front morphology. The critical nonuniformities spacing at which a reaction front begins to merge into a single-finger decreases with increasing upstream pressure gradient. D 2003 Elsevier B.V. All rights reserved. Keywords: Reaction front; Non-uniformities; Finger; Transport. 1. Introduction Water – mineral interaction is important to several branches of geochemistry, including subsurface flow and chemical processes such as mineral dissolution and precipitation, * Corresponding author. Tel.: +886-2-23626480; fax: +886-2-23639557. E-mail address: [email protected] (C.-W. Liu). 0169-7722/$ - see front matter D 2003 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2003.10.006.

(2) 48. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. adsorption and desorption, and redox transformations. Change in aquifer porosity caused by the interaction of reactive pore fluid with the solid matrix has received particular attention (Ortoleva, 1994; Liu et al., 1997; Renard et al., 1998; Lake et al., 2002). A fluid that is undersaturated with respect to minerals in the matrix will cause matrix dissolution and an increase in porosity and typically also permeability. If the initial rock texture is perfectly uniform, the resulting reaction front remains planar. However, geologic media are seldom uniform giving rise to complex, fingered reaction front morphologies. An important question is the extent to which reaction front fingering is developed from nonuniformity like facies lens or sedimentary bedding. Even slight initial non-uniformities, present in all natural rocks, can lead to reaction front unevenness and fingering when the flow is sufficiently rapid (Steefel and Lasaga, 1990). Using macroscopic, non-equilibrium models, previous researchers have shown that a variety of reaction front pattern can form through the coupling of transport and chemical reaction in porous media (Ortoleva, 1994). Reaction front instability has been extensively studied (e.g., Chadam et al., 1986; Ortoleva et al., 1987a,b; Hoefner and Fogler, 1988; Chen and Ortoleva, 1990; Steefel and Lasaga, 1990; Daccord et al., 1993a,b; Lichtner, 1996; Liu et al., 1997; Renard et al., 1998; Logan, 1999). Although reaction front morphology in porous media has been studied by numerous researchers, little attempt has yet been made to examine the effect of two or more nonuniformities in a porous media. Expanding on previous work of a single non-uniformity on the reaction front morphology (Chen and Liu, 2002), this study presents a two-dimensional model for investigating the temporal evolution of aquifer porosities, permeabilities and species concentration during reactive transport in a porous medium with two initial small, local non-uniformities. Simulation results are used to construct a behavior diagram to characterize the front morphology. Quantitative analysis is also performed to evaluate the relative importance of advective and diffusive fluxes on the development of reaction fronts.. 2. Mathematical model This study presents a two-dimensional numerical model used to analyze how reactions between aquifer minerals and pore fluid affect the spatial and temporal evolution of aquifer porosity, permeability, concentrations of aqueous species in solution, and mechanisms affecting that shape of reaction front. Chadam et al. (1986) provide a detailed derivation of the nonlinear partial differential equations that govern the above phenomenon, and are briefly described below. The conservation equations for flow, mass transport, and chemical reaction in a porous medium are formulated using a continuum approach, in which a representative elementary volume (REV) is assumed to be small compared to the phenomena being monitored, and the various properties within the REV are assumed to remain constant. 2.1. Porosity change due to dissolution reactions Although this mathematical model is simplified, generic situations exist that involve a single solid component such as calcite in the porous medium and a single solute such as.

(3) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. 49. calcium or carbonate ion in the fluid. The matrix comprises soluble grains with an average volume of L3 and grain number density of n. The product nL3 is the number of soluble grains that occupy a volume fraction. Meanwhile, the matrix also contains insoluble grains with a volume fraction /n. It follows that: nL3 þ / þ /n ¼ 1. ð1Þ. where / denotes the aquifer porosity. If G represents the rate of grain-volume change due to reaction, then the following occurs: B/ ¼ nG Bt. ð2Þ. An implicit expression for G is required to fully describe the phenomenon. For most minerals, dissolution reaction is very close to the surface-controlled limit (Berner, 1978). Assume that G can be written as: G ¼ CV6L2 ðc ceq Þ ¼ CL2 ðc ceq Þ. ð3Þ. where CV denotes the reaction rate constant, 6L2 represents the surface area of the soluble grains, CV = 6C, c is the concentration of the species in solution, and ceq denotes the equilibrium concentration. The dissolution reaction is assumed to follow first-order kinetics (Lasaga, 1984; Liu and Narasimham, 1989). From Eqs. (1) – (3), the aquifer porosity change that occurs as a result of mineral dissolution reactions can be expressed by the following equation (detailed derivation is provided in Appendix A): B/ ¼ Cn1=3 ð/f /Þ2=3 ðc ceq Þ Bt. ð4Þ. where /f = 1 /n is the final porosity after the complete dissolution of the soluble grains. The 2/3-power indicates that only surface reactions are considered (Lasaga, 1984). 2.2. Conservation of mass and momentum The continuity equation, which describes fluid mass conservation in the system, combined with the momentum conservation equation, Darcy’s law, can be written as: j ð/kð/ÞjpÞ ¼. B/ Bt. ð5Þ. where P denotes the fluid pressure, and k(/) is the intrinsic permeability divided by product of the water viscosity and porosity. A modified form of the Fair –Hatch relation.

(4) 50. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. (Bear, 1972) is employed herein to describe the relationship between porosity and permeability (Chadam et al., 1986; Ortoleva et al., 1987a). kð/Þ ¼. 1 /3 2 2=3 J lh ½ð1 /f Þ þ n1=3 ð/f /Þ2=3 2. where J denotes a packing factor (f5), l represents the fluid viscosity, and h is a geometric factor (f6 for spherical grain). The permeability relationship may also adopt other forms. For instance, Steefel and Lasaga (1990) used the Kozeny – Carman equation to describe this dependence. However, experimental data that justify one form over another are limited. 2.3. Conservation of solute mass A differential equation expressing conservation of species in solution must be able to explain changes in concentration in the fluid phase. These factors that affect the concentration of solute in an aqueous pore fluid are advection, diffusion, and chemical reaction, expressed in the following equation j ½/Dð/Þjc þ /c kð/Þjp þ qs. B/ Bð/cÞ ¼ Bt Bt. ð6Þ. where C represents models per unit volume of porous medium, D(/) denotes the porositydependent diffusion coefficient, and qs represents the density of soluble grains. A common phenomenological relation for D(/) is (Bear, 1972; Lerman, 1979) Dð/Þ ¼ Di /. m. . 3 5 <m< 2 2. . where Di is the species diffusion coefficient in water. The dispersion effects are neglected. This assumption is valid for slow fluid flow. The first, second, and third terms on the lefthand side of Eq. (6) represent the diffusion, advection, and source terms due to the mineral dissolution, respectively. Generally speaking and particular for minerals, the matrix density qs greatly exceeds the solute concentration. Thus the porosity changes appreciably move very slowly during minerals dissolution. Nondimensionalizing a differential equation is often useful, because the resulting equation is generally simpler and more applicable. Consequently a dimensionless time, s, is introduced and defined by s ¼ eðCn1=3 ceq Þt where e = ceq/qs.. ð7Þ.

(5) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. 51. The space variable r = (x,y), concentration c(x,y,t), and pressure p(x,y,t) are also converted to the following dimensionless forms: r¯ ¼. r ½Dð/f Þ=kn1=3 ceq 1=2. c¼. c ceq. p¯ ¼. pkð/f Þ Dð/f Þ. ð8Þ. ð9Þ. ð10Þ. Eqs. (4) – (6) can then be written as (dropping the bars and writing, d(/) = /D(/)/D(/f), k(/) = k(/)/k(/f)), e. B/ ¼ ð/f /Þ2=3 ðc 1Þ Bs. ð11Þ. e. Bð/cÞ B/ ¼ r ðdrc þ kcrpÞ þ Bs Bs. ð12Þ. r ðkrpÞ ¼. B/ Bs. ð13Þ. 3. Numerical methods The numerical model is used to analyze the above reactive transport problem. Efforts to develop numerical reactive transport approaches in recent years have focused on how to couple the reaction and transport term (Yeh and Tripathi, 1991). To solve the coupled set of equations, several methods have been proposed (Rubin, 1983). The most rigorous approach attempts to solve the governing equations simultaneously, which is commonly referred to as the fully coupled method (Steefel and Lasaga, 1994). Alternatively, iteratively coupled techniques can be employed to calculate the reaction and transport. A general iteratively coupled method to solve the coupled set of equations is the sequential iteration approach (SIA) where reactions and transport are solved sequentially (for a discussion, see Yeh and Tripathi, 1991). The SIA method was adopted here. Yeh and Tripathi (1989) discussed explicit and implicit SIA schemes and noted that implicit SIA schemes have faster convergence rates. For this reason, an implicit SIA numerical solution was used in the present study to solve the coupled set of nonlinear differential equations (Eqs. (11)– (13)). The subscript (k) is used to designate quantities at the kth time step. The following steps describe an implicit SIA method that advances a system of algebraic equations representing the differential equations (Eqs. (11)– (13)) by a time interval Ds to the (k + 1)th time level..

(6) 52. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. (1) Eq. (11) can be used to solve for / if c is provided temporarily. Through c, the temporal evolution of / can be obtained by applying the method of separation of variables and integrating with respect to s over the time interval, sk V s V sk + 1, thus obtaining: Z. /kþ1 /k. d/ ð/f /Þ2=3. ¼. Z. skþ1. sk. cnkþ1=2 1 e. ds. ð14Þ. Evaluating these integrals and rearranging the terms yields: ( /nkþ1. ¼ /f . " n # )3 1 ðckþ1=2 1Þ 1 2=3 Ds þ ð/f /k Þ 3 3 e. ð15Þ. where Ds = sk + 1 sk, the subscript k + 1/2 denotes the intermediate time level between k and k + 1, the superscript n denotes the nth iteration in solving for / at time level k + 1. The estimates cnk + 1/2 are assumed to be cnk + 1/2 =(cnk +11 + ck)/2 if n >1, or cnk + 1/2 = ck if n = 1. Evaluating /nk + 1 as /nk + 1 approaches /f is essential. If the choice of Ds is inappropriate, then the computed value can fall in a region where /nk + 1 > /f, thus rendering the computation to be physically nonsensical. To ensure that the computed value of /kn + 1 is within the interval [/0,/f], Ds should be chosen smaller than 1 1/3(/f /k)2/3[e/(cnk +1/2 1)] (Chen and Liu, 2002). As /k approaches /f, a smaller Ds is required to satisfy the numerical stability criterion, increasing the computing time as the simulation progresses. This problem can be circumvented by simply allowing /kn + 1 to be /f, as /f /n 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 an implicit, finite difference method. If pkn + 1is provided temporarily, the computed porosity, /kn + 1, from step 1 can be substituted directly into the difference equation to generate a temporarily known concentration, ckn + 1. The estimates of pnk + 1/2 are assumed to be pnk + 1/2=( pkn +11 + pk)/2 if n >1, or pnk + 1/2 = pk if n = 1. The porosity-dependent permeability and diffusion coefficient are always evaluated at the half-time level using knk + 1/2 = k(/nk + 1/2) and dkn+ 1/2 = d(/nk + 1/2). The spatial derivatives of porosity-dependent permeability and diffusion coefficient in Eq. (12) are also evaluated at the half-time level using a similar calculation. However, numerical errors may appear when we discretize the permeability gradients, Bk/Bx and Bk/By, which have a nonlinear porosity dependence. To overcome these numerical difficulties, an improved method was developed. The gradient is calculated using a semi-analytical derivative instead of the spatial discretization of permeability. The possible error is only the numerical error from discretization of the spatial porosity derivative. The permeability gradient can then be computed from the following: . n . n Bkð/Þ. B/ Bkð/Þ. ¼ Bx kþ1=2 Bx B/ kþ1=2. ð16Þ.

(7) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. Similarly, the diffusion coefficient gradient can be computed from . n . n Bdð/Þ. B/ Bdð/Þ. ¼ Bx kþ1=2 Bx B/ kþ1=2. 53. ð17Þ. By grouping and rearranging terms, the final form of the finite difference equation is: n n n n n Aci;j;kþ1 þ Bciþ1;j;kþ1 þ Cci1;j;kþ1 þ Dci;jþ1;kþ1 þ Eci;j1;kþ1. ¼ Fci;j;k þ Gciþ1;j;k þ Hci1;j;k þ Ici;jþ1;k þ J ci;j1;k þ K. ð18Þ. where the subscript (i,j) denotes the grid center, n. A¼e. 3 B ¼ a4Dx . C¼. a5 a6 a11 a13 ; 4Dy 2Dy2 4Dy. a5 a6 a11 a13 ; þ 2 4Dy 2Dy 4Dy. F ¼ eaDs1 e 3 G ¼a4Dx þ. H ¼. I¼. a4 a7 a9 ; 2Dx2 4Dx. a3 a4 a7 a9 ; þ 4Dx 2Dx2 4Dx. D¼. E¼. n. a1 a2 a4 a6 a8 a9 a10 ki;j;kþ1 a12 a13 a14 ki;j;kþ1 þe þ 2 þ 2 ; Ds 2 Dy Dz 2 2 2 2. a2 a4 a6 a8 a9 a10 ki;j;k a12 a13 a14 ki;j;k 2 2þ þ þ þ ; 2 Dx Dy 2 2 2 2. a4 a7 a9 ; þ 2Dx2 4Dx. a3 a4 a7 a9 þ ; 4Dx 2Dx2 4Dx. a5 a6 a11 a13 þ ; þ 2 4Dy 2Dy 4Dy. J ¼. a5 a6 a11 a13 þ ; 4Dy 2Dy2 4Dy. K ¼ b1 ;.

(8) 54. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. where a1 ¼. /ni;j;kþ1 þ /i;j;k 2. a 2 ¼ b2 ¼. ;. /ni;j;kþ1 /i;j;k ; 2. Bd. n a3 ¼ ; Bx i;j;kþ1=2 a 4 ¼ a6 ¼. a5 ¼. n di;j;kþ1 þ di;j;k. a9 ¼. kni;j;kþ1 þ ki;j;k 2. ;. n Bk. : Bx i;j;kþ1=2 pniþj;kþ1=2 pnij;kþ1=2 2Dx. a10 ¼. a12. ;. n Bd. ; By i;j;kþ1=2. a7 ¼ a11 ¼. a8 ¼. 2. ;. p;niþ1;j;kþ1=2 2pnj;kþ1=2 þ pni1;j;kþ1=2 Dx2. ;. n Bk. ¼ ; By i;j;kþ1=2. a13 ¼. a14 ¼ b1 ¼. pniþ1;j;kþ1=2 pni1;j;kþ1=2 2Dy. ;. pni;jþ1;kþ1=2 2pni;j;kþ1=2 þ pni;j1;kþ1=2 Dy2. ;. /ni;j;kþ1 /i;j;k : Ds. n (3) Similarly, Eq. (13) can be discretized to evaluate the fluid pressure, pi,j,k + 1 , if the n previously obtained value of /i,j,k + 1 is substituted into the finite difference equation.

(9) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. 55. Table 1 The variables and parameters used in this study and their meaning and dimensions Variable 3. L n / /n G CVor C c ceq p k(/) l J h D(/) qs s c p¯ e. Meaning. Dimension. average volume grain number density aquifer porosity volume fraction of insoluble grain rate of grain-volume change due to reaction reaction rate constant concentration of the species in solution equilibrium concentration pressure porosity-dependent permeability function fluid viscosity packing factor geometric factor diffusion coefficient density of soluble grain dimensionless time dimensionless pressure dimensionless concentration e=ceq/qs. [L3] [L3] [–] [–] [L3 T1] [M1 L4 T1] [M L3] [M L3] [M L1 T2] [M1 L3 T] [M L1 T1] [–] [–] [L2 T1] [M L3] [–] [–] [–] [–]. (Eq. (13)). The final form after grouping and arranging can be expressed in the following form Apni;j;k þ Bpniþ1;j;k þ Cpni1;j;k þ Dpni;jþ1;k þ Epnij;j1;k ¼ F 2c2 2c4 ; Dx2 Dy2 c3 c4 þ E¼ ; 2Dy Dy2 A¼. B¼. c1 c2 þ ; 2Dx Dx2. C¼. c1 c2 þ ; 2Dx Dx2. ð19Þ D¼. c3 c4 þ ; 2Dy Dy2. F ¼ c5. n n After implementing the procedure of steps 1 –3, the computed solution /i,j,k + 1 , ci,j,k + 1 , n and pi,j,k + 1 may be inaccurate and can be used as the new guess value for the next solving iteration. n n (4) The procedure outlined in steps 1 –3 for solving for /ni,j,k + 1, ci,j,k + 1 , and pi,j,k + 1 is repeated iteratively until the following convergence criterion has been met:. nþ1. ð20Þ. ðki;j;kþ1 kni;j;kþ1 Þ=kni;;j;kþ1 Vfk. max. where E refer to f, g, or p; fk is a specified residue constant; and the subscript max denotes the maximum value over all grid centers (Table 1). 4. Simulation result 4.1. Temporal –spatial evolution of the reaction front A numerical model of water –rock interaction, NSPCRT, has been developed by Chen and Liu (2002) and is applied here to investigate the instability morphology of.

(10) 56. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. Fig. 1. Numerical simulation domain. Inlet fluid is injected in from the left wall and flows out of the right wall.. the reaction front during transport. This code has capacity to simulate twodimensional temporal – spatial changes in mass transfer, porosity, permeability, and reaction front morphology induced by infiltration of reactive fluids into porous rocks, and thus is well suited for modeling water – rock interaction. The scenario considered in this study is isothermal flow and reaction in carbonate cemented sandstone within a rectangular problem domain (Fig. 1). In the model, carbonate cement is considered reactive (for example, a carbonate cement) whereas the remaining aquifer matrix is inert (for example, a quartz at low temperature). Fig. 1 illustrates the model boundary conditions. Constant pressure was assigned to the vertical boundaries of the grid, producing a constant pressure gradient ( pfV) that led to a constant fluid velocity in the sandstone. The boundary condition for pressure at x = 0 is thus prescribed as: Bp ¼ pf V Bx. ðx ¼ 0Þ. ð21aÞ. The imposed flow is undersaturated with respect to the reactive mineral phase, and the value of the inlet concentration is assumed a constant value of c0. Chadam et al. (1986) and Steefel and Lasaga (1990) use a Dirichlet (the first type) boundary for concentration at x = 0 as c ¼ c0. ðx ¼ 0Þ. ð21bÞ. Only the net pressure drop across the region of interest is important since the fluid density is assigned to be constant, and thus the outlet of the system is prescribed. It is appropriate to use boundary conditions such that p ¼ pR ðyÞ ¼ 0;. Bc ¼ 0; Ax. ðx ¼ Lx Þ. ð22Þ.

(11) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. 57. No flow boundary condition has been assigned to the horizontal boundary of the grid, i.e. Bp ¼ 0; By. Bc ¼ 0; By. ðy ¼ Ly ; Ly Þ. ð23Þ. To illustrate reaction front instability, the temporal – spatial evolution of aquifer porosity and species concentration in a porous medium with two initially small nonuniformities is investigated. The non-uniformities are zone of higher porosity and permeability located the left boundary of the system, as shown in Fig. 2. The initial conditions for the concentration and porosity were set as follows: /ðx; y; 0Þ ¼ /0 þ ð/f /0 Þen. ð24Þ. cðx; y; 0Þ ¼ ð1 e5x Þð1 en Þ. ð25Þ. where n(x,y)=[x4+( y F a)4]/(wLy)4, w denotes an initial perturbation parameter which is used to perturb the initial concentration and porosity, and 2a represents the spacing of two non-uniformities. The input parameters used to simulate the morphological instability of the reaction front were: initial porosity, /0 = 0.1; final porosity, /f = 0.2; inlet concentration, c0 = 0; length, Lx = 18; width, 2Ly = 8; and initial perturbation parameter, w = 0.1. A fine grid spacing (Dx = 0.2 and Dy = 0.2) in both the x- and y-directions was used to avoid numerical dispersion. The two parameters of particular interest in this study are the upstream pressure gradient (inlet velocity) and spacing of the two small initial non-uniformities. A series of simulations was performed for different upstream gradients and spacing of non-uniformities. Fig. 2a– c displays the temporal – spatial evolution of simulated aquifer porosity and species concentration contours (calcium or carbonate ion) for upstream pressure gradients of 0.5, 2 and 5, respectively, with a fixed non-uniformities spacing of 1.5. The solid and dashed lines denote the aquifer porosity contours with (/0 + /f)/2 and the dimensionless concentration contours with (c0 + cf)/2, respectively. Different reaction front morphologies are observed. The reaction front is planar for an upstream pressure gradient of 0.5, but becomes increasingly unstable for larger upstream pressure gradients of 2 and 5. Interestingly, a double-finger reaction front emerges with an upstream gradient of 5, while a single-finger reaction front emerges for an upstream pressure gradient of 2. The upstream pressure gradient governs the development of different reaction front morphologies. The planar reaction front results from the fact that diffusion dominates over the flow-focusing mechanism, and such that the two initial higher permeability bumps are not amplified. Consequently, the reaction front is planar. Chadam et al. (1986) explained such results as being caused by the inlet velocity being smaller than the primary bifurcation value using bifurcation theory. The unstable single- or double-finger reaction fronts indicate that the inlet velocity exceeds the primary bifurcation value, or that the diffusion effect cannot overcome.

(12) 58. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66.

(13) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. 59. the flow-focusing mechanism. Consequently, the flow tends to be focused through the tips of the bumps, and two bumps simultaneously exist in the reaction zone. The fact that the concentration is smaller inside the bump (concave on the upstream side) than in neighboring regions causes diffusion from the outside bump (convex on the downstream side) into the inside bump. For the single-finger, the diffusion effect cannot inhibit the self-focusing mechanism, but drags the two fingers to merge into a wider single finger. However, the double-finger front means that the strong self-focusing mechanism caused by the fast inlet flow exceeds the finger-merging effect driven by diffusion, thus maintaining the double-finger front. In terms of bifurcation theory, the critical upstream pressure gradient (inlet flow value) at which the fingertip merging begins is termed a secondary bifurcation value. To further investigate the phenomenon, Fig. 3 presents the porosity contours at the same upstream pressure gradient of 0.5, 2 and 5, but with the two initial nonuniformities spacing being increased to 2 instead of 1.5. The simulation results are the same as that for upstream gradient 0.5 and 5, but the reaction front morphology is different for an upstream gradient of 2. Obviously, a double-finger emerges instead of a single-finger. This result is attributed to the fact that the enlarged non-uniformity spacing mitigates the finger-merging mechanism caused by diffusion and a self-focusing mechanism prevails, thus maintaining a double-finger. Fig. 4 presents the porosity contours at the same non-uniformity spacing of 2 but with an upstream pressure gradient of 1. Two non-uniformities gradually develop into a single-finger, confirming that the decreased upstream pressure gradient reduces the self-focusing mechanism and causes the double-finger to merge into a single-finger. Figs. 5 and 6 provide further insight into the development of reaction front shape. Displayed in these figures are the spatial distributions of the diffusive, advective and resultant species fluxes resulting from non-uniformity spacings of 1.5 and 2, respectively, at a fixed upstream pressure gradient of 2.0, and time of 1.4. The arrows and their adjacent values denote the direction and magnitude of the species flux transport, while the dashed line represents the 0.5 concentration contour. Although diffusion tends to pull the fingertip back, advection tends to travel along the fingertip, and thus the resultant of both advective and diffusive species flux determines the shape of the reaction front. The above discussion on reaction front shape reveals that the interaction between a selffocusing flow mechanism and a diffusion-caused finger-merging mechanism can drive two local non-uniformities into different fingering morphologies. Not only can the planar front become unstable, but single- or double-finger fronts can also emerge. Shape selection is. Fig. 2. Temporal evolution of aquifer porosity and species concentration (described by the lines /=(/0 + /f)/ 2 = 0.15 and c=(c0 + cf)/2 = 0.5) with (a) upstream pressure gradient = 0.5, (b) upstream pressure gradient = 2.0 and (c) upstream pressure gradient = 5.0 at a fixed non-uniformities spacing of 1.5. Solid line: aquifer porosity, dashed line: species concentration, the numbers on the curves are the dimensionless time, and the numbers on the x- and y-axes are the dimensionless distance scale..

(14) 60. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. Fig. 3. Temporal evolution of aquifer porosity (/=(/0 + /f)/2 = 0.15) and species concentration c=(c0 + cf)/2 = 0.5 with (a) upstream pressure gradient = 0.5, (b) upstream pressure gradient = 2.0, (c) upstream pressure gradient = 5.0 at a fixed non-uniformities spacing of 2. Solid line: aquifer porosity, dashed line: species concentration..

(15) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. 61. Fig. 4. Temporal evolution of aquifer porosity (/=(/0 + /f)/2 = 0.15) and species concentration c=(c0 + cf)/2 = 0.5 with upstream pressure gradient = 2 at a fixed non-uniformities spacing of 2. Solid line: aquifer porosity, dashed line: species concentration.. determined by the competition between the flow-focusing mechanism and tip-merging effect driven by diffusion. 4.2. Morphological behavior diagram The criteria that determine the development of a planar, single- or double-finger reaction front morphology are of particular interest. These criteria can be conveniently illustrated using a ‘‘behavior diagram’’ (Daccord et al., 1993b). For the present problem, such a diagram delineates the control space domain for which various reaction morphologies exist. The control parameters here include the spacing of the upstream pressure gradient and the two non-uniformities. Fig. 7 displays a two-dimensional behavior diagram plot that uses an upstream pressure gradient and spacing of two non-uniformities as two variables, while the values of the other parameters are fixed. The behavior diagram defines the occurrence conditions of the planar front, single-finger, double-finger and transition zone. The development of a planar front depends only on the upstream pressure gradient, while the single- and double-finger developments are controlled by both the upstream pressure gradient and the non-uniformities spacing. A transition zone lies among the regions of the planar front, and the single- and double-finger fronts. The transition processes are relatively slow. Although the reaction front is unstable in the transition zone, the development of a single- or double-finger front is not clearly pictured and is difficult to distinguish. As illustrated in Fig. 4, a single-finger front slowly developed at the boundary of the transition zone, where the upstream pressure gradient is 1 and the non-uniformities spacing equals 2. Meanwhile, a double-finger front occurred in conditions of medium to large non-uniformity spacing and a high upstream gradient, while a single-finger front developed in conditions of medium non-uniformity spacing and a medium to high upstream pressure gradient. Notably, the behavior diagram can.

(16) 62. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. Fig. 5. Quantitative analysis of competition between advection and diffusion with upstream pressure gradient = 2.0 and two non-uniformities of 1.5 at time = 1.4..

(17) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. 63. Fig. 6. Quantitative analysis of competition between advection and diffusion with upstream pressure gradient = 2.0 and two non-uniformities of 2.0 at time = 1.4..

(18) 64. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. Fig. 7. Behavior diagram showing reaction front morphology is a function of upstream pressure gradient and spacing of non-uniformities.. define the morphology of the reaction fronts for any combination of values of these two parameters.. 5. Conclusion This study investigates the instability of a flow-driven reaction front that occurs via a flow-focusing mechanism during reactive transport in a porous medium with two small, local non-uniformities. A two-dimensional mathematical model is used to simulate this phenomenon, and a series of numerical simulations were conducted to examine the morphology of such reaction fronts. Interestingly, two local non-uniformities in an otherwise uniform system may cause the reaction front to evolve into quite complicated patterns. The planar front occurs when the velocity of the incoming fluid is slow. Meanwhile, the finger-shape front emerges when the upstream velocity exceeds the critical value and destabilizes the reaction front. Moreover, the reaction front initiated by two local non-uniformities may become a double-finger or merge into a single-finger depending on the spacing of two initial local non-uniformities and the upstream pressure gradient. Quantitative analysis of diffusive, advective and resultant species fluxes of both.

(19) J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. 65. mechanisms are computed and presented to further elucidate the behavior of the reaction morphology instability. Furthermore, a behavior diagram is presented that defines the reaction front morphology as a function of upstream pressure gradient and the spacing of the two local non-uniformities. Other control parameters, such as the domain width, can also be used to construct a behavior diagram in delineating the reaction front morphology.. Acknowledgements The authors are grateful to Professors de Marsily and Appold for their detailed helpful suggestions and comments during the review process. We would like to thank the National Science Council of the Republic of China for financially supporting this work under contract nos. NSC 87-2313-B-002-051 and NSC 88-2313-B-002-097.. Appendix A In this appendix, Eq. (4) is derived via Eqs. (1) –(3). First, Eq. (1) can be rearranged as L¼. 1 / /n n. . 1=3 ¼. /f / n. 1=3 :. ðA1Þ. Inserting Eq. (A1) into Eq. (3) yields /f / 2=3 G¼C ðc ceq Þ: n. ðA2Þ. Eq. (4) can be obtained by substituting Eq. (A2) into Eq. (2) and expressed as B/ ¼ Cn1=3 ð/f /Þðc ceq Þ: Bt. ðA3Þ. References Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, Amsterdam. 764 pp. Berner, R.A., 1978. Rate control of mineral dissolution under earth surface conditions. Am. J. Sci. 278, 1235 – 1252. Chadam, J., Ortoleva, P., Sen, A., 1986. Reactive-infiltration instabilities. IMA J. Appl. Math. 36, 207 – 220. Chen, J.S., Liu, C.W., 2002. Numerical simulation of the evolution of aquifer porosity and species concentrations during reactive transport. Comput. Geosci. 28, 485 – 499. Chen, W., Ortoleva, P., 1990. Reaction front fingering in carbonate-cemented sandstones. Earth-Sci. Rev. 29, 183 – 198. Daccord, G., Lenormand, R., Lietard, O., 1993a. Chemical dissolution of a porous medium by a reactive fluid: I. Model for the ‘‘wormholing’’ phenomenon. Chem. Eng. Sci. 48, 169 – 178. Daccord, G., Lietard, O., Lenormand, R., 1993b. Chemical dissolution of a porous medium by a reactive fluid: II. Convection vs. reaction, behavior diagram. Chem. Eng. Sci. 48, 179 – 186..

(20) 66. J.-S. Chen, C.-W. Liu / Journal of Contaminant Hydrology 72 (2004) 47–66. Hoefner, M.L., Fogler, H.S., 1988. Pore evolution and channel formation during flow and reaction in porous media. AIChE J. 34, 45 – 54. Lake, L.W., Bryant, S.L., Araque-Martinez, A.N., 2002. Geochemisty and Fluid Flow. Elsevier, Amsterdam. 226 pp. Lasaga, A.C., 1984. Chemical kinetics of water – rock interactions. J. Geophys. Res. B6, 4009 – 4025. Lerman, A., 1979. Geochemical Processes. Wiley, New York. 481 pp. Lichtner, P.C., 1996. Continuum formulation of multicomponent – multiphase reactive transport. In: Lichtner, P.C., Steefel, C.I., Oelkers, E.H. (Eds.), Reactive Transport in Porous Media, Reviews in Mineralogy, vol. 34. Mineralogical Society of America, Washington D.C., pp. 1 – 81. Liu, C.W., Narasimham, T.N., 1989. Redox-controlled multiple-species reactive chemical transport, I. Model Development. Water Resour. Res. 25, 869 – 882. Liu, X., Ormond, A., Bartko, K., Li, Y., Ortoleva, P., 1997. A geochemical reaction-transport simulator for matrix aciding analysis and design. J. Pet. Sci. Eng. 17, 181 – 196. Logan, J.D., 1999. Reaction fronts in porous media with varying porosity: an exact solution. Comm. Appl. Nonl. Analy. 6 (4), 45 – 50. Ortoleva, P.J., 1994. Geochemical Self-Organization. Clarendon Press, Oxford. 411 pp. Ortoleva, P., Merino, E., Moor, C., Chadam, J., 1987a. Geochemical self-organization: I. Feedback mechanisms and modelling approach. Am. J. Sci. 287, 979 – 1007. Ortoleva, P., Chadam, J., Merino, E., Sen, A., 1987b. Self-organization in water – rock interaction systems: II. The reactive-infiltration instability. Am. J. Sci. 287, 1008 – 1040. Renard, F., Gratier, J.-P., Ortoleva, P., Brosse, E., Bazin, B., 1998. Self-organization during reactive fluid flow in a porous medium. Geophys. Res. Lett. 25 (3), 385 – 388. Rubin, J., 1983. Transport of reacting solutes in porous media: relation between mathematical nature of problem formulation and chemical nature of reactions. Water Resour. Res. 19, 1231 – 1252. Steefel, C.L., Lasaga, A.C., 1990. Evolution of dissolution patterns: permeability change due to coupled flow and reactions. In: Melchior, D.C., Bassett, R.L. (Eds.), Chemical Modeling in Aqueous Systems II: American Chemical Society Symposium, Washington D.C., vol. 416, pp. 125 – 212. Steefel, C.L., Lasaga, A.C., 1994. A coupled model for transport of multiple chemical species and kinetic precipitation/dissolution reactions with application to reactive flow in single phase hydrothermal systems. Am. J. Sci. 294, 529 – 592. Yeh, G.T., Tripathi, V.S., 1989. A critical evaluation of recent developments of hydrogeochemical transport models of reactive multichemical components. Water Resour. Res. 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 Resour. Res. 27 (11), 3075 – 3094..

(21)