International Journal of Heat and Mass Transfer 75 (2014) 483–496 Contents lists available at ScienceDirect International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt Pore-scale study of diffusion–reaction processes involving dissolution and precipitation using the lattice Boltzmann method Li Chen a,b, Qinjun Kang b,⇑, Bill Carey b, Wen-Quan Tao a a b Key Laboratory of Thermo-Fluid Science and Engineering of MOE, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM, USA a r t i c l e i n f o Article history: Received 17 September 2013 Received in revised form 13 March 2014 Accepted 25 March 2014 Available online 6 May 2014 Keywords: Mass transport Chemical reaction Pore scale Lattice Boltzmann method Dissolution Precipitation a b s t r a c t A pore-scale model combining the lattice Boltzmann method (LBM) and a fluid–solid interface tracking method is employed to simulate the diffusion–reaction processes involving dissolution and precipitation. Coupled sub-processes including mass transport, chemical reactions, and solid structure evolution are considered. Effects of the precipitation of the secondary solid phase on the dissolution of the primary solid phase are investigated under different dissolution–precipitation reaction kinetics, molar volumes of the primary and secondary solid phases, powder size, surface roughness, and nucleation and crystal growth mechanisms. Different morphologies of the precipitates are predicted by the pore-scale simulations. It is found that the precipitation has opposite effects on the underlying dissolution processes. The favorable effect is that the precipitation reaction consumes the product of the dissolution reaction, thus facilitating the dissolution; while the adverse effect is that the generated precipitates cover the surface of the primary solid phase, thus separating the reactive surface from the reactive components. Based on the extent to which the precipitates affect the dissolution, four types of coupled dissolution–precipitation processes are identified and discussed. Ó 2014 Elsevier Ltd. All rights reserved. 1. Introduction Reactive transport processes involving dissolution– precipitation are pervasive in a variety of scientific, industrial and engineering processes. Typical examples include self-assembled patterns such as Liesegang rings or bands [1], formation of mineral deposits in boilers and heat exchangers, biofilm growth in aqueous environment [2], environmental contaminant transport [3,4], recovery of oil and geologic sequestration of carbon dioxide (CO2) in the subsurface [5–7]. Among these reactive processes, it is commonly encountered that a second solid phase precipitates when a primary solid phase dissolves, and the precipitation and dissolution reactions are closely coupled with each other [8,9]. For example, during mineral trapping of CO2, primary silicate mineral dissolves due to decrease of pH caused by the dissolved CO2; meanwhile the dissolved CO2 can react with cations released by the dissolution reaction to form a secondary precipitate of carbonate mineral [6]. Another typical example is the uptake of Cd2+ by carbonate minerals from polluted water. As the CaCO3 dissolves, ⇑ Corresponding author. Address: Computational Earth Science Group (EES-16), Mail Stop T003, Los Alamos National Laboratory, Los Alamos, NM 87545, USA. Tel.: +1 505 665 9663; fax: +1 505 665 8737. E-mail addresses: lichennht@gmail.com (L. Chen), qkang@lanl.gov (Q. Kang), bcarey@lanl.gov (B. Carey), wqtao@mail.xjtu.edu.cn (W.-Q. Tao). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.03.074 0017-9310/Ó 2014 Elsevier Ltd. All rights reserved. it releases carbonate leading to the aqueous solution to be supersaturated with respect to otavite which then precipitates [4]. These non-linear non-equilibrium transport processes with dynamic evolutions of solid structures pose great challenges to the numerical simulations. The numerical models and methods proposed should be capable of simulating the fluid flow, predicting mass transfer, incorporating homogeneous and heterogeneous reactions, updating the solid structures, and taking into account the interactions between different sub-processes. Usually, for simulating such complex reactive transport processes with the solid structure evolution, a solver for transport processes, such as the finite volume method (FVM), is coupled with a model for tracking the fluid–solid interface, such as the phase-field method (PF). Conventionally, the prediction of transport processes is based on solving the macroscopic density, momentum, energy, and concentration conservation equations using a FV, finite difference (FD), or finite element (FE) method. As these methods are continuum based, it is not easy to properly handle the discontinuity of variables at the complex fluid–solid interfaces, limiting their applications for transport processes in domains with complicated structures such as porous media. Pore-scale methods have been developed which have the capacity of taking into account the complex porous structures, such as pore network model [10], smoothed particle hydrodynamics [11,12] and the lattice Boltzmann method (LBM) [13–17]. Evolutions of fluid–solid interfaces 484 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 are ubiquitous in processes including melting–solidification and dissolution–precipitation. For tracking the fluid–solid interfaces, the PF method [18] and the cellular automaton (CA) methods [19] are widely used. Other methods such as Volume of Fluid (VOF) and Level Set (LS) method, which are commonly used for multiphase flow, are also adopted for this purpose [20]. There have been several pore-scale studies in the literature to study the reactive transport processes involving dissolution and precipitation by coupling different transport process solvers and fluid–solid interface tracing methods since the work of Ref. [21]. Single-phase fluid flow and reactive transport with dissolution– precipitation was studied by Kang et al. [22,23] in which the LBM was used to simulate flow, mass transport and reaction, and the volume of pixel (VOP) was adopted to track the moving fluid–solid interfaces caused by dissolution–precipitation. They predicted the relationship between permeability and porosity under different Pe (Peclet number, representing the relative strength of convection to diffusion) and Da (Damköhler number, representing the relative strength of reaction to diffusion) numbers. Li et al. [20] studied a similar problem using LS to track the fluid–solid interfaces. In addition, Luo et al. [24] implemented a model using a diffuse interface method to track the fluid–solid interface, in which they also considered the natural convection caused by concentration gradients. Smoothed particle hydrodynamics was also adopted for investigating dissolution–precipitation processes in fractures and porous media [11,12]. Later, Kang et al. extended their model to multi-component systems [14] and used the model to study reactive transport processes associated with geological CO2 sequestration [15]. Yoon et al. [16] adopted LBM for fluid flow and FVM for mass transport in mixing – induced calcium carbonate dissolution and precipitation processes. Huber et al. [17] combined LBM and PF method for dissolution–precipitation processes involving single or multiple species. All the above studies are for systems of a single phase fluid with one or multiple chemical species dissolved in it. The numerical studies of multiphase fluid flow coupled with reactive transport with moving fluid–solid interfaces are scarce in the literature [25,26]. Recently, Parmigiani et al. [25] used the LBM to study the process of injection of a non-wetting fluid into a wetting fluid coupled with dynamic evolution of the solid geometries. Very recently, Chen et al. [26] constructed a pore-scale model based on the LBM and the VOP to simulate multiphase reactive transport with phase transition and dissolution–precipitation processes [26]. Their pore-scale model can capture coupled non-linear multiple physicochemical processes including multiphase flow with phase separation, mass transport, chemical reaction, dissolution–precipitation, and dynamic evolution of the pore geometries. The above pore-scale studies provide a detailed insight into the coupled mechanisms between the transport processes and the dissolution–precipitation processes. To the best of our knowledge, however, there are few studies taking the coexistence of dissolution and precipitation processes into account [14,15,26], and there are no studies especially devoted to investigating the coupled dissolution–precipitation processes where the precipitates may cover the surface of the solid phase undergoing dissolution. Existing experiments have shown that the way by which the precipitates affect the dissolution of primary solid phase is very complex and cannot be generalized. Instead, it is affected by several factors including the dissolution and precipitation reaction kinetics, concentrations of reactive components, and the nucleation and crystal mechanisms of the precipitates [3,4,6,7]. The objective of the present study is to numerically investigate the effects of the precipitation of the secondary solid phase on the dissolution of the primary solid phase by performing comprehensive studies including dissolution and precipitation reaction kinetics, molar volume ratio, powder size, surface roughness, and nucleation and crystal growth mechanism. The LBM is adopted for solving the diffusion–reaction processes and the VOP is employed for updating the solid structures. With more than 20 year’s development, the LBM has emerged as a powerful tool for the numerical simulations and investigations of a broad class of complex flows, including porous flow, thermal flow, reactive transport, turbulence flow, and multiphase flow. Since it is based on the discrete kinetic theory, the LBM is a promising tool in dealing with complicated non-linear characteristics as well as complex structures. The VOP developed by Kang et al. [13] is a CA method. The VOP has many advantages such as a clear physical concept, simple and stable arithmetic, easy implementation of various kinds of surface reactions, and flexible coupling with different nucleation and crystal growth mechanisms. For this reason, the method has been used successfully to predict many moving solid–fluid interface phenomena, such as crystal growth [13], rock dissolution due to acid injection [22,23], Liesegang bands or rings [1] and dissolution and precipitation involved in CO2 sequestration[15,16]. The rest of the paper is arranged as follows. In Section 2 the physicochemical model is presented. In Section 3 the numerical method is introduced. In Section 4 two-dimensional (2D) simulation results are presented and discussed. Finally, in Section 5 some conclusions are drawn. 2. Physicochemical model 2.1. Chemical reaction A simplified chemical model is constructed that can be readily used for studying transport and reaction kinetic parameters and for identifying dominant processes and mechanisms. It can be described by the following two principal dissolution and precipitation steps with three aqueous species and two solid phases AðaqÞ þ DðsÞ ¡BðaqÞ BðaqÞ þ C ðaqÞ ¡PðsÞ ð1Þ ð2Þ where aq in the parenthesis stands for aqueous species and s denotes solid phase. Eq. (1) is the dissolution reaction in which aqueous A(aq) reacts with solid D(s) generating aqueous B(aq). The generated B(aq) reacts with another aqueous C(aq) according to precipitation reaction Eq. (2), producing the secondary precipitate P(s). Note that elements comprising the secondary precipitate may either be initially present in the aqueous solution or result from the dissolution of the primary phase. In the present study, the element B(aq) required for precipitation reaction is released from the dissolution, resulting in close coupling between the precipitation and dissolution reactions. Eqs. (1) and (2) can be considered as a simplified form of some realistic reactions in experiments [4,6]. For example for carbonation of silicates, A(aq) can be considered as H+, D(s) as the silicate, B(aq) as the cation, C(aq) as the carbonate or bicarbonate and P(s) as the calcite. For uptake of Cd2+ by carbonate minerals, A(aq) can be considered as H+, D(s) as carbonate mineral, B(aq) as the carbonate, C(aq) as the toxic cation and P(s) as the otavite. Although Eqs. (1) and (2) are not the same as the reactions in [4,6], they retain the basic principles of these realistic reactions. 2.2. Governing equation In the present study, fluid flow is not considered for the purpose to focus on the coupled mechanisms between diffusion and dissolution–precipitation reactions. The general processes can be described as follows. Initially, primary solid phase D(s) in equilibrium with A(aq) and B(aq) is placed in the computational domain, and there is no secondary solid phase P(s) (see Fig. 1). Initial concentration of C(aq) is zero and thus the initial condition is not 485 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 Fig. 1. Schematics of the computational domain. supersaturated with respect to P(s). Then, a solution including A(aq) and C(aq) is introduced into the domain from the left inlet. As the inlet concentration of A(aq) is higher than its initial concentration, solid D(s) dissolves and releases B(aq). Precipitation reaction of B(aq) and C(aq) then occurs. The diffusion–reaction equations for aqueous species A(aq), B(aq) and C(aq) are as follows @C AðaqÞ @t @C BðaqÞ @t @C C ðaqÞ @t ¼ DAðaqÞ DC AðaqÞ ð3aÞ ¼ DBðaqÞ DC BðaqÞ  SHN ð3bÞ ¼ DC ðaqÞ DC C ðaqÞ  SHN ð3cÞ where t is time, C is the concentration and D is the diffusivity. The source term SHN is related to the homogeneous nucleation of the SHN ¼ 8 <0 if qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  : 1 C B þ C C  ðC B þ C C Þ2  4ðC B C C  K n Þ if ðaqÞ ðaqÞ ðaqÞ ðaqÞ ðaqÞ ðaqÞ 2dt precipitates and will be introduced in Section 2.4. 2.3. Dissolution reaction On the surface of D(s), dissolution reaction takes place which consumes A(aq) and releases B(aq) as described by Eq. (1). This reaction is described by the surface reaction in the simulation DAðaqÞ DBðaqÞ @C AðaqÞ @n @C BðaqÞ @n ¼ k1 1  K 1 C BðaqÞ C AðaqÞ ¼ k1 1  K 1 ! C BðaqÞ C AðaqÞ ð4aÞ ! ð4bÞ where k1 and K1 are reaction rate and equilibrium constant, respectively. n is the direction normal to the reactive surface pointing to the void space. The density of the corresponding surface D(s) node is given as follows @dDðsÞ @t ¼ AM DðsÞ k1 1  K 1 C BðaqÞ C AðaqÞ !, dV ð5Þ where A is the reactive area, dV is the volume of each node and M is the molar mass. Here density means the total mass at a solid node with a fixed volume which equals the control volume of a computational node, and is not the physical density of the solid phase. as a kinetically hindered, autocatalytic process [27]. This means that the precipitation reaction between components B(aq) and C(aq) has two barriers. One is related to the kinetics of the nucleation process. If there are no precipitates or other impurities that can catalyze the precipitation, then at first this barrier has to be overcome. The nucleation process is described by a step-function Eq. (6). If the product of concentrations of B(aq) and C(aq) reaches a nucleation threshold Kn, nucleation will start. The other barrier, called thermodynamic barrier, is related to the crystal growth. If there are already some precipitates or impurities present in the system, then crystal growth on the surface of these precipitates or impurities will occur once the product of concentrations of B(aq) and C(aq) exceeds a crystal growth threshold Kc. The nucleation has a high threshold than the crystal growth, and both of them require the development of supersaturated conditions with respect to P(s). 2.4.1. Nucleation It is widely accepted that homogeneous nucleation requires significant degrees of supersaturation for overcoming the high surface tension [28]. In the present study, for homogeneous nucleation at a certain node in the computational domain, the product of concentrations of B(aq) and C(aq) must reach the nucleation threshold Kn. Besides, the homogeneous nucleation at this node is considered as a probabilistic process and only happens with a specified probability kp. The homogeneous nucleation described by the following equation, is considered to be much faster than mass transport, and is included in Eq. (3) as a volume source where dt is the time step. The above step function is widely used in C BðaqÞ C C ðaqÞ < K n the fields of studying precipitation patterns [27], and has been proved to generate various precipitate patterns such as regular bands and rings [27]. Once the homogeneous nucleation occurs at a certain node, this node is marked as a solid precipitate node and is assigned with an initial density d0, which is calculated based on the mass balance d0 ¼ SHN M BðaqÞ dt þ SHN M CðaqÞ dt ð7Þ Subsequent crystal growth is allowed on this solid node. 2.4.2. Crystal growth Crystal growth takes place either on the surface of D(s) or on that of P(s). It is reasonable to distinguish the crystal growth threshold Kc on different surfaces due to different surface energies. Hence, two thresholds Kc,D and Kc,P are introduced to denote the crystal growth barrier on the surface of D(s) and P(s), respectively. The crystal growth process around existing surface is taken into account using the crystal growth model developed in [13,14]. In this model, crystal growth is described by surface reaction. On the surface of a precipitate node, B(aq) and C(aq) react, leading to the consumption of B(aq) and C(aq) in the vicinity and the increment of the node density. The surface reaction is expressed as DC ðaqÞ @C CðaqÞ @n 0 ¼ ( ¼ ( 2.4. Nucleation and crystal growth Precipitation is initiated by nucleation, and subsequently takes place as crystal growth. In this study, the precipitation is regarded ð6Þ C BðaqÞ C C ðaqÞ > K n DBðaqÞ @C BðaqÞ @n if C C ðaqÞ C BðaqÞ < K c k2 ð1  K 2 C C ðaqÞ C BðaqÞ Þ if C C ðaqÞ C BðaqÞ P K c 0 if C C ðaqÞ C BðaqÞ < K c k2 ð1  K 2 C C ðaqÞ C BðaqÞ Þ if C C ðaqÞ C BðaqÞ P K c ð8aÞ ð8bÞ 486 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 where k2 and K2 are reaction rate and equilibrium constant, respectively. Correspondingly, the density of the precipitate node is updated by @dPðsÞ @t  . ¼ AM P k2 1  K 2 C C ðaqÞ C BðaqÞ dV ð9Þ Experiments found that the re-dissolution of the secondary precipitate may happen which would greatly affect the reactive transport processes [16]. In order to focus on the coupling mechanism between precipitation of the secondary mineral and dissolution of the primary mineral, re-dissolution of the secondary precipitate is not considered in the present study, as can be seen from Eq. (8). 3. Numerical method 3.1. LB model for solute transport The governing equations of A(aq), B(aq) and C(aq) are solved using the LB model for solute transport [29,30] g ak ðx þ ea dt; t þ dtÞ ¼ g ak ðx; tÞ  1 sk ðg ak ðx; tÞ  g eq ak ðC k ; uÞÞ þ dtJak Sak ð10Þ where gak is the distribution function of the kth component with velocity ea at the lattice site x and time t, and s is the collision time related to the diffusivity. Here, we use the D2Q5 lattice model for solute transport, as it has been shown to have not only comparable accuracy with the D2Q9 model but also better computational efficiency [30–33]. S denotes the source term due to homogeneous reaction. For A(aq), S = 0. For B(aq) and C(aq), S is given by Eq. (6). The equilibrium distribution function geq takes the following form [30] g eq ak ¼ C k J ak ð11Þ where J is given by [30] Ja ¼  J0 ; a¼0 ð12Þ ð1  J 0 Þ=4; a ¼ 1; 2; 3; 4 The rest fraction J0 can be selected from 0 to 1 to obtain different diffusivity. The equilibrium distribution function given by Eq. (12) can cover a wide range of diffusivity by adjusting J0, which is a prominent advantage of such an equilibrium distribution [1,32–35]. Species concentration C and diffusivity are obtained by Ck ¼ X g ak ð13Þ 1 dx2 Dk ¼ ð1  J 0 Þðsk  0:5Þ 2 dt ð14Þ where dx is the length of one lattice. 3.2. Update of the density of solid node The solid density (volume) described by Eqs. (5) and (9) is updated at each time step based on the following equation dDðsÞ ðt þ DtÞ ¼ dDðsÞ ðtÞ  M DðsÞ k1 1  K 1 C BðaqÞ ! dtA=dV C AðaqÞ   dPðsÞ ðt þ DtÞ ¼ dPðsÞ ðtÞ  M PðsÞ k2 1  K 2 C BðaqÞ C CðaqÞ dtA=dV ð15aÞ ð15bÞ Assuming the depth of a computational cell in the z direction (perpendicular to the paper) equals the lattice length dx, the term A/dV in Eq. (15) is reduced to dx1. For dissolution, if dD reaches zero, this solid node is removed and is converted to a fluid node. The mass of P(s) accumulated so far on this node due to crystal growth is redistributed to neighboring solid nodes, so the precipitation reactions and mass accumulation can continue at these nodes. Such a treatment can guarantee the mass conservation [14]. Here, we give all the mass to a single neighboring solid node. For precipitation, when the density of a precipitate node (parent node) exceeds a certain threshold value (2d0 in the present study), a nearest neighboring fluid node (child node) becomes a solid node. There may be several nearest neighboring fluid nodes around the parent node. Different ways of choosing the child node may lead to different crystal growth patterns. This means our model has an extra degree of freedom, which allows for incorporation of different crystal growth mechanisms. In the base case of the present study, the child node is randomly chosen. In Section 4.6 other rules are also studied to investigate the effects of crystal growth mechanisms on the coupled dissolution–precipitation processes. After the solid node is added, the density of the parent node is set back to d0; meanwhile the child node is assigned with an initial density d0 and is marked as a ‘‘crystal growth node’’ as it results from crystal growth [14]. 3.3. Computational domain, initial and boundary conditions The 2D computational domain is a rectangular channel with the size of a 140  90 lm, as shown in Fig. 1. It is discretized with a 140  90 lattice. Here it is worth mentioning that the dimension of a lattice node is about 1 lm, which is somewhat greater than the size of a typical crystal [36]. Thus, the kp mentioned in Section 2.4.1 serves as the probability for the nucleation process occurring at the nano or micro scale to actually result in a stable nucleus at this lattice node, and its value is set as 0.8. Initially, there is no P(s) in the system. Two rectangles representing powders of D(s) are placed adjacent to the top and bottom surface of the channel. The size of each rectangle is 120  30 lm. Initially, A(aq) and B(aq) are in equilibrium with D(s), and the concentration of C(aq) is zero. So there is no dissolution of D(s) or precipitation of P(s). Then, a solution containing only A(aq) and C(aq) is introduced at the left inlet. Because the concentration of A(aq) at the inlet, CA,in, is higher than the initial concentration of A(aq), the initial equilibrium among A(aq), B(aq) and D(s) is broken and D(s) starts to dissolve according to Eq. (1) at the D(s)-fluid interface. Table 1 Parameters for the base case [14,16]. Variable Value of physical unit Value of lattice unit Diffusivity Molar volume, D Vm 1.0  109 m2 s1 3.0  105 m3 mol1 0.2 0.6 Molar volume ratio of P(s) to D(s) Inlet concentration of A(aq) Inlet concentration of B(aq) Inlet concentration of C(aq) Reaction constant Reaction 1 Reaction constant of Reaction 2 Equilibrium constant of Reaction 1 Equilibrium constant of Reaction 2 Homogeneous nucleation threshold Crystal growth threshold on the surface of D(s) Crystal growth threshold on the surface of P(s) u 1 1 CA,in 2.0 M (2000 mol m3) 0.1 CB,in 0M 0 CC,in 2.0 M 0.1 k1 1.0  10 1 1  103 k2 1.0  103 mol m2 s1 1  105 K1 1.0  102 1  102 K2 1.25  105 m6 mol2 Kn,solution Kc,D Kc,P 5 1 mol m 2 2 4  10 mol m (assumed) 0.8  105 mol2 m6 (assumed) 0.8  105 mol2 m6 (assumed) 6 s 5000 1  103 2  104 2  104 487 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 Fig. 2. Effects of k1. Temporal and spatial evolutions of geometries of D(s) and P(s) for different k1: (a) Case 1 (k1 = 1  101 mol m2 s1) and (b) Case 2 (k1 = 1  103 mol m2 s1).Red node is P(s) and black node denotes D. (c) Time variations of volume of D(s) and P(s) for different k1. B(aq) is thus generated and its concentration increases. It further reacts with C(aq) from the inlet and produces the precipitate P(s). From the above description of the diffusion–reaction processes, it is clear that the precipitation and dissolution reactions are strongly interacted with each other. Besides, the generated precipitates, which may cover the surface of D(s), will greatly affect the mass transport processes. Obviously, the diffusion, dissolution and precipitation processes are closely coupled with each other. At the right outlet, concentration gradient of all solutes is set to zero. Both top and bottom boundaries are non-reactive solid walls with zero diffusion flux, which is achieved by bounce-back scheme with no collision step in the LB framework. At the D(s) (or P(s))-fluid interface, where dissolution or precipitation reaction takes place, the boundary conditions given by Eqs. (4) or (8) are enforced through the method proposed in Ref. [31]. The values of variables for the base case are listed in Table 1. Validation of the codes can be found in our previous studies [26,31,33], and are not repeated here for brevity. 4. Results and discussion In this section, effects of dissolution and precipitation kinetics, molar volume of D(s) and P(s), initial powder size of D(s), surface roughness of D(s), and crystal growth mechanisms on the coupled dissolution and precipitation processes are investigated in detail. Entirely, twenty three cases are designed and simulated for Table 2 Volumes of dissolution and precipitation for different cases. Case Dissolution volume Precipitation volume Case Dissolution volume Precipitation volume Case 1 Case 2 Case 3 Case 4 Case 5 Case 6 Case 7 Case 8 0.314 0.180 Case 9 1.000 0.000 0.621 0.174 1.000 0.233 0.107 0.103 0.023 0.062 0.963 0.063 0.538 0.303 0.920 0.015 0.644 0.330 0.960 0.092 1.000 0.087 0.521 0.247 1.000 0.277 0.122 0.100 Case 10 Case 11 Case 1_A Case 2_A Case 5_A Case 10_A Case 11_A 0.045 0.201 Bold Roman: Type 1; Bold italics: Type 2; Roman: Type 3; Italics: Type 4. different purposes. Case 1 is the base case. Cases 2–8 are for understanding the effects of chemical reaction kinetics on the dissolution–precipitation processes. In Case 2, the dissolution reaction rate constant is altered to study its effects. In Cases 3, 4 and 5, the precipitation reaction rate constant is varied. Cases 6, 7 and 8 are for the purpose of investigating effects of equilibrium constant. 488 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 In addition, different molar volume ratios between the secondary and primary solid phases are considered in Cases 9–11. Further, a different initial structure of the primary solid phase is used to investigate the effects of structures in Case 1_A, Case 2_A, Case 5_A, Case10_A and Case 11_A. Besides, different crystal growth mechanisms are taken into account in Case 1_B, Case 10_B, Case 1_C and Case 10_C. The effects of solid-surface dependent crystal growth barrier are studied in Case 10_D1 and Case 10_D2. Finally, in Case 10_E the homogeneous nucleation process, which is not considered in other cases, is added. For cases whose name has a letter suffix (e.g., A, B, C, D, E), the values of various parameters used in the simulations are the same as those in their corresponding cases whose name does not contain such a letter suffix. Without explicitly stated, the values of variables in different simulation cases are the same as the base case. The results of primary interest are the temporal and spatial evolutions of the geometries of D(s) and P(s), as well as the rate and amount of dissolution and precipitation. 4.1. Effects of reaction rate constant 4.1.1. Effects of k1 Fig. 2(a) shows the temporal and spatial evolutions of geometries of D(s) and P(s) for Case 1 (base case). In all the figures throughout the paper, D(s) is denoted by black nodes and P(s) is marked by red nodes, unless otherwise stated. In case 1, k1 equals 1  101 mol m2 s1 with the dimensionless number Dad of 4.5. Here Dad is defined as k1d/DCA,in, with the characteristic length d as the domain height (90 lm), which reflects the relative strength of the dissolution reaction rate to the diffusion process. The dimensionless number Dap defined as k2d/DCA,in, which is the relative strength of the precipitation reaction rate to the diffusion process, is 4.5  102. The reactive transport process in this case is actually a diffusion-controlled process, meaning that the diffusion process is slow compared to reaction. Therefore, the dissolution of D(s) initially occurs on the left surface of D(s) (t = 20 s), as shown in Fig. 2(a). As the dissolution progresses, more and more B(aq) is Fig. 3. Effects of k2. Temporal and spatial evolutions of geometries of D(s) and P(s) for different k2: (a) Case 3 (k2 = 1  101 mol m2 s1), (b) Case 4 (k2 = 5  105 mol m2 s1) and (c) Case 5 (k2 = 1  105 mol m2 s1). (d) Time variations of volume of D(s) and P(s) for different k2. L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 released, and the increased B(aq) combined with C(aq) from the left inlet results in the supersaturated condition with respect to P(s), leading to precipitation of P(s) (t = 20 s). Initially the structure of P(s) is rather porous and has limited effects on the underlying dissolution of D(s) which thus continuously dissolves. However, P(s) becomes increasingly compact due to continuous crystal growth, and finally completely inhibits the dissolution at the left surface of D(s) (t = 100 s). Thus, the dissolution front gradually moves to the right where precipitates are still absent. Finally, the entire surface of D(s) is coated by the precipitates and the dissolution stops (t = 200 s). The precipitation also ceases, because D(s) is separated from A(aq) by P(s) and no more B(aq) is released. Fig. 2(b) shows the temporal and spatial evolutions of geometries of D(s) and P(s) for Case 2, where k1 is decreased to 1  103 mol m2 s1, only a hundredth of that in Case 1. This time Dad is 4.5  102 and dissolution is a reaction-controlled process, implying that diffusion is very quick compared to reaction. Thus, dissolution of D(s) occurs fairly uniformly over the surface of D(s). Correspondingly, the precipitates also exhibit a relatively uniform distribution (t = 200 s). Note that the concentration of B(aq) near the inlet is very small, even lower than that near the outlet. Therefore, P(s) is undersaturated near the left surface of D(s). This is why a small fraction of the surface of D(s) near the inlet is free of precipitates (t = 200 s). As the precipitates become increasingly compact, local underlying dissolution of D(s) is inhibited and the dissolution front moves back to the left (t = 3000 s). Finally, the precipitates coat the entire surface of D(s) and both the dissolution and precipitation cease (t = 6000 s). Fig. 2(c) shows the time evolution of volumes of D(s) and P(s) for different k1. In the figure the volume is normalized by the initial volume of D(s). The slope of the curves can be considered as an indication of the dissolution (precipitation) rate. As can be seen in 489 Fig. 2(c), both the dissolution and precipitation rates increase as k1 rises. The volume of dissolution increases as k1 reduces, due to slow generation of precipitates which allows more D(s) to be dissolved, as shown in Fig. 2(c) and listed in Table 2. The volume of P(s) is almost the same for different k1. 4.1.2. Effects of k2 Fig. 3(a)–(c) displays the temporal evolutions of geometries of D(s) and P(s) for k2 of 1  101 mol m2 s1 (Case 3), 5  105 mol m2 s1 (Case 4) and 1  105 mol m2 s1 (Case 5), respectively. The dimensionless number Dap is 4.5, 2.25  103 and 4.5  104, respectively. Compared with Case 1, in Case 3 the rate of precipitation is so high that the precipitates quickly cover the entire surface of D(s) and completely stop the dissolution (Fig. 3(a)). Decreasing k2 leads to porous precipitates (Fig. 3(b)). When k2 is very small, the precipitates are so rare that the dissolution can be completed (Fig. 3(c)). Note that the dissolution is still ongoing in Fig. 3(b) and (c) and D(s) can be completely dissolved eventually. Fig. 3(d) shows the time variations of volumes of D(s) and P(s) for different k2. More D(s) is dissolved as k2 reduces, which is expected based on the above discussion. Interestingly, the amount of precipitation firstly goes up and then reduces as k2 decreases (as shown in Fig. 3(d) and listed in Table 2), which is attributed to the coupled dissolution and precipitation process. This implies that a moderate precipitate reaction rate is favorable for both precipitation and dissolution. 4.2. Effects of equilibrium constant (solubility) We change K1 in a wide range (1  101–1  104), and find that the dissolution is far from equilibrium over the entire range due to the fast consumption of B(aq). Thus, changing K1 alone has Fig. 4. Effects of K2. Geometry of D(s) and P(s) for different K2: (a) t = 400 s for Case 6 (K2 = 1.25  106 m6 mol2), (b) t = 350 s for Case 7 (K2 = 2.5  106 m6 mol2) and (c) t = 50 s for Case 8 (K2 = 1.25  103 m6 mol2). Note that the geometries of (c) is final, while that of (a) and (b) are intermediate state (D(s) will finally be completely dissolved in (a) and (b)). (d) Time variations of volume of D(s) and P(s) for different K2. 490 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 Fig. 5. Effects of molar volume ratio. Temporal and spatial evolutions of geometries of D(s) and P(s) with different molar volume ratio of P(s) to D(s): (a) Case 9, u = 0.001, (b) Case 10, u = 0.5 and (c) Case 11, u = 5. (d) Time variations of volume of D(s) and P(s) for different molar volume ratio. slight effects on dissolution and precipitation processes. The simulation results are not presented here for brevity. Fig. 4 shows the geometries of D(s) and P(s) at certain time for K2 of 1.25  106 m6 mol2 (Case 6), 2.5  106 m6 mol2 (Case 7) and 1.25  103 m6 mol2 (Case 8). In the figure, the geometries for Case 8 are the final one, while those for Cases 6 and 7 are in intermediate state. Note that K1 and K2 have different units as can be seen from the Eqs. (4) and (8). Higher K2 means lower solubility of P(s) and the solution is more easily to be supersaturated with respect to P(s); thus, the precipitation is quick and the structures of P(s) are more compact, as can be observed in Fig. 4(c). The coupled mechanism between the dissolution and precipitation can be clearly observed in Fig. 4(d) where the time variations of volumes of D(s) and P(s) are plotted. A moderate K2 leads to the maximum amount of precipitation, as shown in Fig. 4(d) and listed in Table 2. 4.3. Effects of molar volume ratio Molar volume ratio of the secondary precipitate phase to the primary dissolution phase plays an important role on the precipitation and dissolution processes. In this section the molar volume ratio of P(s) to D(s) (u) is set as 0.001 (Case 9), 0.5 (Case 10) and 5 (Case 11), respectively (The molar volume of D(s) is fixed), and the simulation results of temporal evolutions of geometries of D(s) and P(s) are shown in Fig. 5(a)–(c). Combined with Case 1 (u = 1) shown in Fig. 2(a), it can be seen that changing u leads to rich variations of the geometric evolutions. When molar volume of P(s) is extremely small (Case 9), no precipitate appears and the dissolution surface is smooth, as shown in Fig. 5(a). As the molar volume of P(s) increases to half that of D(s) (Case 10), very porous precipitates form (Fig. 5(b)). Further increasing u to 5 (Case 11), the precipitation is so strong that precipitates form armors on the surface of the D(s) in a short period of time and completely suppress the dissolution, as shown in Fig. 5(c). Actually, the precipitation has two opposite effects on the dissolution. On one hand, the precipitation is favorable for dissolution as it consumes B(aq) and keeps the dissolution reaction away from equilibrium, facilitating the dissolution reaction. On the other hand, the precipitates impede the transport of reactant towards D(s) surface, thus slowing down the dissolution. Whether L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 the favorable or detrimental factor dominates depends on the reaction kinetics/parameters as well as the geometrical characteristics. In Cases 9 and 10, the favorable factor is dominant, thus the precipitation facilitates the dissolution, as more clearly shown in Fig. 5(d). Such favorable effect was also found in experiments [37]. Murakami et al. [37] found in their experiment that solution was saturated with respect to the primary mineral, but the precipitation of the second mineral increased the dissolution rate, because the precipitate layer is porous and precipitation reaction removed elements generated by dissolution reaction and thus reduced the saturation rate of the solution with respect to the primary mineral. The precipitate in [37] is Boehmite with a molar volume of 2  105 m3 mol1 and the primary dissolving solid phase is Anorthite with a molar volume of 2.7  104 m3 mol1. The molar volume ratio of Boehmite to Anorthite is about 0.074, leading to a porous precipitate layer generated. For u = 1 (Case 1), there is a transition from the favorable effect to the detrimental effect, and the dissolution rate is very quick initially and then slows down. For u = 5 (Case 11), the detrimental factor dominates and the precipitation inhibits the dissolution. 491 4.4. Identification of different types of coupled dissolution– precipitation process Before further discussing other factors, here we try to classify different types of coupled precipitation and dissolution processes predicted so far, based on the extent to which the precipitation affects the dissolution. Four types are identified based on the simulation results. In Type 1, the precipitate generated is so rare, due to low precipitation reaction rate, high solubility or small molar volume of the precipitate, that the primary solid phase D(s) can be completely dissolved. Cases 5, 6 and 9 belong to such type. In Type 2, more precipitates are generated; but they are too porous to inhibit the transport of reactant so that D(s) can still be completely dissolved. Cases 4, 7 and 10 can be classified into this type. For Type 2, the favorable effect of precipitation on dissolution dominates and the dissolution is accelerated. This type leads to both maximum amount of precipitation and dissolution, as listed in Table 2. Further increasing the precipitate reaction rate, decreasing the solubility of the precipitate or raising the molar volume of the precipitate Fig. 6. Effects of initial size of D(s) power. Temporal and spatial evolutions of geometries of D(s) and P(s): (a) Case1_A, (b) Case2_A, (c) Case5_A, (d) Case10_A and (e) Case11_A. Time evolutions of volume of D(s) and P(s): (f) Case 1, Case 2, Case 1_A and Case 2_A and (g) Case 5, Case 10, Case 11, Case 5_A, Case 10_A and Case 11_A. 492 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 leads to Type 3. In this type, precipitates are sufficiently compact that they finally stop the dissolution completely. The dissolution is only allowed to some extent. Cases 1 and 2 are such type. In this type, the transition from favorable effect of precipitation on dissolution to the detrimental effect can be clearly observed. In Type 4, the precipitation is very quick that the precipitates are extremely compact, thus coating the entire surface of the primary solid phase and completely suppressing the dissolution. Cases 3, 8 and 11 present such characteristics. In this type, the detrimental effect of precipitation on dissolution overwhelms. In both Types 3 and 4, the surface of the primary mineral is eventually covered by the secondary precipitate completely. 4.5. Effects of initial powder size In this section, the effects of the initial size of D(s) powder are explored. The two powder particles in the above simulations are divided into 12 smaller powder particles while their total volume is fixed. Selective simulations are performed based on the four types of dissolution and precipitation processes identified. In the four sets of simulations performed in this section, the parameters are the same as that in Cases 1, 5, 10 and 11, respectively. An additional simulation is performed for the reaction-controlled dissolution and the parameter is the same as in Case 2. Correspondingly, the five simulations are called Case 1_A, Case 2_A, Case 5_A, Case 10_A and Case 11_A. Fig. 6(a)–(e) shows the time evolutions of the geometries of D(s) and P(s) for the five cases, respectively. As can be observed, the characteristics of the four types of dissolution and precipitation are generally well maintained. Obviously, decreasing the powder size increases the reactive surface area, which is favorable for the dissolution. The initial surface area for the powder size in Fig. 6(a)–(e) is about three times that of the cases studied in Sections 4.1–4.3. For dissolution and precipitation of Type 3 (Case 1 and Case 1_A, Case 2 and Case 2_A), both the rate and amount of dissolution and precipitation increase as the powder size decreases as shown in Fig. 6(f) and listed in Table 2, due to the increased reactive surface area. In Fig. 6(a), the narrower void space between smaller powders, however, increases the risk of pore clogging by the generated precipitates. The upstream precipitates clog the transport path, leading to no dissolution for the middle right powder. For dissolution and precipitation of Type 1, the increase of dissolution and precipitation rate (Case 5 and Case 5_A) (Fig. 6(g)) due to decreased powder size is not as remarkable as that for Type 3. For Type 2 (Case 10 and Case 10_A, Fig. 6(d)), however, decreasing the powder size has little effect on the rate of dissolution and precipitation (Fig. 6(g)), suggesting the significant accelerating effects of precipitation on dissolution (the favorable effect discussed in Section 4.3). For type 4 (Case 11 and Case 11_A, Fig. 6(e)), the effect of changing powder size is also slight (Fig. 6(g)) since the precipitates completely covers the entire surface of primary solid phase in a short period of time. Dissolution and precipitation with initially rough surface of D(s) are also studied. The surface roughness is generated by arranging rectangular protuberances on the powder surface. The general evolutions of geometries of D(s) and P(s) are similar to that without surface roughness. Besides, increasing surface roughness has similar effects as decreasing the powder size, because both of them increase the surface area of the powder. 4.6. Effects of crystal growth mechanisms of P(s) In the above simulations, the crystal growth is a random process. In this section, the effects of crystal growth mechanisms are studied. The emphasis is placed on Case 1 and Case 10. In addition, only the geometries of D(s) and P(s) are discussed in detail while little attention is paid to the rate and amount of dissolution and precipitation. 4.6.1. Crystal growth towards the direction of maximum concentration gradient For this kind crystal growth mechanism, when the density of an existing precipitate node is increased to 2d0, the nearest neighboring fluid node towards which the local concentration gradient achieves the maximum is chosen to become a new precipitate node [38]. The concentration gradient is calculated by rðC BðaqÞ C CðaqÞ Þ ¼ Fig. 6 (continued) ðC BðaqÞ C CðaqÞ Þnb  C BðaqÞ C C ðaqÞ jxnb  xj ð16Þ where the subscript ‘nb’ denote neighboring. Thus, the crystal will always grow towards the direction of highest concentrations. Fig. 7(a) and (b) shows the simulation results for Case 1_B and Case 10_B, where the parameters are the same as that in Case 1 and Case 10, respectively. The morphology of P(s) generated significantly differs from that of Case 1 in Fig. 2(a) and of Case 10 in Fig. 5(b). It can be seen that the crystals orient towards the surface of D(s), which is expected because C BðaqÞ C C ðaqÞ is high near the D(s)-fluid interface where B(aq) is generated. Besides, the crystals are inclined to agglomerate, because crystals that once randomly grow now gather together to grow directionally towards D(s)-fluid interface. L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 493 Fig. 7. Temporal evolutions of geometries of D(s) and P(s) with different crystal growth mechanisms: crystal grows towards the maximum concentration gradient (a) Case1_B and (b) Case 10_B, and crystal grows with epitaxial growth mechanism (c) Case 1_C and (d) Case 10_C. 4.6.2. Epitaxial crystal growth Prieto et al. [3] and Cubillas et al. [4] found that the identity and relative structures of the precipitate and primary minerals play an important role on the coupled dissolution and precipitation processes. For primary solid phase and secondary precipitates having similar structures and unit cell parameters, the precipitation is mainly through the mechanism of epitaxial growth, in which the precipitates quickly spread and cover the surface of the primary mineral, thus only a small amount of precipitates can completely stop the dissolution. On the other hand, for primary mineral and secondary precipitates presenting different structures, the precipitates are mainly formed through the mechanism of random growth, and rather porous structures of precipitates are likely to form, allowing the underlying dissolution to continue. The crystal growth mechanism in the above simulations (except that in Section 4.6.1) is a random growth mechanism. In this section, the following scheme is proposed to simulate the epitaxial growth in our simulations. The barrier for crystal growth on the surface of D(s) and P(s) is still set as the same, which is reasonable because the structures and the unit cell parameters of D(s) and P(s) are required to be similar in epitaxial growth mechanism. When the density of an existing precipitate node is increased to 2d0, the nearest neighboring fluid node which is directly adjacent to the D(s) surface is chosen for crystal growth; if there is no such node, then the random growth mechanism is adopted. Through such a scheme, the crystal is expected to grow on the surface of D(s) to the largest extent, which is the distinct feature of the epitaxial growth mechanism. Fig. 7(c)–(d) shows the simulated geometries for Case 1_C and Case 10_C, where the parameters are the same as that in Case 1 and Case 10, respectively. It can be seen that the precipitates tend to grow on the surface of D(s) and quickly covers the surface of D(s). In Case 1 (Fig. 2(a)), several layers of precipitates are required to completely passivate the surface of D(s); while in Case 1_C (Fig. 7(c)), only one layer of precipitate is sufficiently dense to separate the surface of D(s) from the fluid region. For Case 10, the precipitates are so porous that D(s) can be thoroughly dissolved; however in Case 10_C, the undissolved D(s) is besieged by P(s) (Fig. 7(d)). On the whole, the inhibition effects of epitaxial growth mechanism on the underlying dissolution are evident and efficient, consistent with the experimental results in Ref. [3]. 4.6.3. Effects of crystal growth barrier In practice, the crystal growth barrier for different surfaces is different due to different surface energies. Here, it is assumed that the crystal is inclined to grow on the surface of P(s). The crystal growth barrier on D(s) surface is set as 1.5Kc,D and 2Kc,D with the value of Kc,D given in Table 1. The crystal growth barrier on P(s) surface, Kc,P, is fixed as that in Table 1. Only two cases called Case 10_D1 (1.5Kc,D) and Case 10_D2 (2Kc,D) (Values of other variables are the same as Case 10) are considered, and the final geometries of D(s) and P(s) are displayed in Fig. 8. For comparison, the final geometry for Case 10 is also presented in Fig. 8. In the figure, the crystals grown on the D(s) surface are denoted by red squares and these on P(s) surface are denoted by green deltas. In the final geometries, there is no D(s) because for Case 10_D1 and Case 10_D2 P(s) is still discontinuous enough that D(s) is completely dissolved. Note that here we use the term ‘‘discontinuous’’ while in Section 4.3 we use the term ‘‘porous’’ when describing the structures of P(s) for Case 10. The discontinuous pattern may be due to either the 494 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 Fig. 8. Final geometries of D(s) and P(s) for crystal growth with different thermodynamic barrier on the surface of the primary mineral. Red square denotes crystal growing on the surface of D(s) and green delta denotes crystal growing on the surface of P(s). (a) Case 10 with Kc,D, (b) Case 10_D1 with 1.5 Kc,D, and (c) Case 10_D2 with 2.0 Kc,D. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 9. Temporal and spatial evolutions of geometries of D(s) and P(s) with homogeneous nucleation. Red square denotes precipitate crystal growing on the surface of D(s) and P(s), Black square denotes D(s) and blue delta denotes precipitate due to homogeneous nucleation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) porous structures of P(s) (Fig. 8(a)) or patchy distributions of compact precipitates (Fig. 8(c)). As shown in Fig. 8, as crystal growth barrier on D(s) surface increases, fewer crystals directly grow on D(s) surface and thus concentrations of B(aq) and C(aq) around these crystals are higher, leading to subsequently more luxuriant growth of P(s) on its own surface. The precipitates are porous themselves when Kc,D is low (Case 10, Fig. 8(a)), but are bigger, more compact and discontinuous as Kc,D increases (Case 10_D1 and 10_D2, Fig. 8(b) and (c)). These simulations clearly explain similar observations obtained in recent experiments where precipitation on the primary phase becomes rate-limited, leading to large separated agglomerates [39]. 4.7. Effects of homogeneous nucleation In the above simulations, the homogeneous nucleation threshold Kn is set as 4  105 mol2 m6, which is high enough to turn down the homogeneous nucleation of P(s) in the solution. Here, the Kn is reduced to 2.4  105 mol2 m6 which allows the homogeneous nucleation to occur. Case 10_E corresponding to Case 10 is simulated and the temporal evolutions of geometries of D(s) and P(s) are presented in Fig. 9. Obviously, the homogeneous nucleation probability kp will affect the homogeneous precipitation patterns. It is set as 0.8 in our simulation, which is an assumed value. In Fig. 9, black squares denote D(s) and blue deltas represent precipitates resulting from homogeneous nucleation. Note that in Fig. 9 the crystals grown on the surface of D(s) or P(s) are not distinguished, and are both denoted by red squares. Compared with Case 10 in Fig. 5 (b), it can be seen that the homogeneous nucleation results in significantly different dissolution patterns. Since only Kn is changed and inlet concentrations of B(aq) and C(aq) are fixed, the maximum concentration product of B(aq) and C(aq) still locate at the D(s)-fluid interface. Thus, homogeneous nucleation first occurs at the left surface of D(s) (t = 60 s) and completely covers the underlying D(s), therefore inhibiting the dissolution. The dissolution front is forced to move downstream. Due to the decreased concentration product of B(aq) and C(aq), the Kn is no longer exceeded and the precipitates are generated mainly by crystal growth on the surface of either D(s) or P(s) (t = 300 s and t = 600 s). 5. Conclusion There have been extensive experimental studies concerning the effects of precipitation of the secondary solid phase on the dissolution of the primary solid phase. However, such effects are debated as conflicting observations in natural systems and in the laboratory have been reported. Nugent et al. [40] suggested that a thin, patchy and natural coating of aluminosilicate partially inhibited the field dissolution rate of feldspar. However, Lee et al. [41] concluded that such coating was discontinuous and played a limited role in reducing the field dissolution rate of feldspar. In the laboratory, Stockmann et al. [7] found that prisms of carbonate were formed on or adjacent to the surface of basaltic glass, and the precipitates did not affect the dissolution rate. Daval et al. [6] carried out carbonation reactions of wollastonite and found that the compact continuous coatings of small calcite crystals formed in circum-neutral pH greatly inhibited the dissolution of wollatonite. An amorphous silica layer has been widely observed during mineral trapping of CO2 by silicate mineral. While Stockmann et al. [7] and Daval et al. [6] reported that such amorphous silica layer was very porous and had slight effect on the dissolution of the underlying mineral, Bearat et al. [42] and Alfredsson et al. [43] pointed out that such layer completely covered the surface of the primary mineral and inhibited the further dissolution. In addition, Prieto et al. [3] and Cubillas et al. [4] suggested that the precipitate mechanism (epitaxial or random precipitation), which was dependent on the similarity between the structures of the primary and the secondary minerals, played an important role on the morphology of the precipitates and thus on the underlying dissolution. There have been some experimental studies suggesting that the precipitation increased the dissolution rate as the precipitation consumed elements released by the dissolution and thus reduced the saturation rate of the solution with respect to the primary mineral [37,44]. The numerical studies in the present study, together with the above experimental studies, suggest that how the precipitation affects the dissolution of primary mineral cannot be generalized and depends on the morphologies of the precipitates. The pore-scale simulations demonstrate how the coupled diffusion–reaction processes are affected by the dissolution and L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 precipitation reaction kinetics, physical characteristic of the primary and secondary minerals and the nucleation and crystal mechanisms of the precipitates. The simulation results reveal the coupled mechanisms between mass transport, dissolution and precipitation. It is found that the precipitation process has opposite effects on the dissolution. The favorable one is that it facilitates the dissolution reaction by consuming the product of the dissolution reaction. Meanwhile, the precipitates impede the transport of reactant towards the surface of the primary solid phase, thus slowing down the dissolution. Depending on the extent to which the precipitates affect the dissolution, four types of coupled dissolution–precipitation processes are identified under different physicochemical conditions. Particularly, Type 2, in which the dissolution can be completed and the precipitation amount is relatively high, indicating a desirable process for the CO2 sequestration and uptake of Cd+ from polluted water. Besides, it is also found that the crystal growth mechanisms, the crystal growth barrier and the homogeneous nucleation also play important roles on the morphologies of the precipitates, consistent with the experiments [3,39]. If the crystal grows towards the direction of maximum concentration gradient, the precipitates always grow towards the surface of the primary solid phase. If epitaxial crystal growth is the underlying crystal growth mechanism, the precipitates tend to cover the surface of the primary solid phase and hinder the dissolution reaction. In addition, as the crystal growth barrier on the primary solid phase increases, the precipitates become bigger, more compact and more discontinuous. With homogeneous nucleation considered, the precipitates quickly cover the surface of the primary solid phase and force the reaction sites to move towards the outlet. Conflict of interest We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled. Acknowledgments We thank the National Nature Science Foundation of China (No. 51136004) for the support of this work. Q. Kang is grateful for the support from the LANL LDRD Program and Institutional Computing Program. We also appreciate the thorough reviews of Dr. Christian Huber and other two anonymous reviewers. References [1] L. Chen, Q. Kang, Y.-L. He, W.-Q. Tao, Mesoscopic study of the effect of gel concentration and material on the formation of precipitation patterns, Langmuir 28 (32) (2012) 11745–11754. [2] T.R.R. Pintelon, C. Picioreanu, M.C.M. van Loosdrecht, M.L. Johns, The effect of biofilm permeability on bio-clogging of porous media, Biotechnol. Bioeng. 109 (4) (2012) 1031–1042. [3] M. Prieto, P. Cubillas, Á. Fernández-Gonzalez, Uptake of dissolved Cd by biogenic and abiogenic aragonite: a comparison with sorption onto calcite, Geochim. Cosmochim. Acta 67 (20) (2003) 3859–3869. [4] P. Cubillas, S. Köhler, M. Prieto, C. Causserand, E.H. Oelkers, How do mineral coatings affect dissolution rates? An experimental study of coupled CaCO3 dissolution—CdCO3 precipitation, Geochim. Cosmochim. Acta 69 (23) (2005) 5459–5476. [5] M. Hänchen, V. Prigiobbe, R. Baciocchi, M. Mazzotti, Precipitation in the Mg– carbonate system—effects of temperature and CO2 pressure, Chem. Eng. Sci. 63 (4) (2008) 1012–1028. [6] D. Daval, I. Martinez, J. Corvisier, N. Findling, B. Goffé, F. Guyot, Carbonation of Ca-bearing silicates, the case of wollastonite: experimental investigations and kinetic modeling, Chem. Geol. 265 (1–2) (2009) 63–78. 495 [7] G.J. Stockmann, D. Wolff-Boenisch, S.R. Gislason, E.H. Oelkers, Do carbonate precipitates affect dissolution kinetics? 1: basaltic glass, Chem. Geol. 284 (3–4) (2011) 306–316. [8] C. Zhu, A.E. Blum, D.R. Veblen, Feldspar dissolution rates and clay precipitation in the Navajo aquifer at Black Mesa, Arizona, USA, in: R.B. Wanty, R.R.I. Seal (Eds.), The Eleventh International Symposium on Water-rock interaction WRI1, Saratoga Springs, New York, 2004, pp. 895–899. [9] C. Zhu, P. Lu, Z. Zheng, J. Ganor, Coupled alkali feldspar dissolution and secondary mineral precipitation in batch systems: 4. Numerical modeling of kinetic reaction paths, Geochim. Cosmochim. Acta 74 (2010) 3963–3983. [10] L. Li, C.A. Peters, M.A. Celia, Upscaling geochemical reaction rates using porescale network modeling, Adv. Water Resour. 29 (9) (2006) 1351–1370. [11] A.M. Tartakovsky, P. Meakin, T.D. Scheibe, R.M. Eichler West, Simulations of reactive transport and precipitation with smoothed particle hydrodynamics, J. Comput. Phys. 222 (2) (2007) 654–672. [12] A.M. Tartakovsky, G. Redden, P.C. Lichtner, T.D. Scheibe, P. Meakin, Mixinginduced precipitation: experimental study and multiscale numerical analysis, Water Resour. Res. 44 (6) (2008) W06S04. [13] Q. Kang, D. Zhang, P.C. Lichtner, I.N. Tsimpanogiannis, Lattice Boltzmann model for crystal growth from supersaturated solution, Geophys. Res. Lett. 31 (2004) L21604. [14] Q. Kang, P.C. Lichtner, D. Zhang, Lattice Boltzmann pore-scale model for multicomponent reactive transport in porous media, J. Geophys. Res. 111 (B05203) (2006). [15] Q. Kang, P. Lichtner, H. Viswanathan, A. Abdel-Fattah, Pore scale modeling of reactive transport involved in geologic CO2 sequestration, Transp. Porous Media 82 (1) (2010) 197–213. [16] H. Yoon, A.J. Valocchi, C.J. Werth, T. Dewers, Pore-scale simulation of mixinginduced calcium carbonate precipitation and dissolution in a microfluidic pore network, Water Resour. Res. 48 (2) (2012) W02524. [17] C. Huber, B. Shafei, A. Parmigiani, A new pore-scale model for linear and nonlinear heterogeneous dissolution, precipitation and sorption reactions, Geochim. Cosmochim. Acta 124 (2014) 109–130. [18] L.-Q. Chen, Phase-filed models for microstructures evolution, Annu. Rev. Mater. Res. 32 (2002) 113–140. [19] D. Sun, M. Zhu, S. Pan, D. Raabe, Lattice Boltzmann modeling of dendritic growth in a forced melt convection, Acta Mater. 57 (6) (2009) 1755–1767. [20] X. Li, H. Huang, P. Meakin, Level set simulation of coupled advection–diffusion and pore structure evolution due to mineral precipitation in porous media, Water Resour. Res. 44 (12) (2008) W12407. [21] S. Békri, J.F. Thovert, P.M. Adler, Dissolution of porous media, Chem. Eng. Sci. 50 (17) (1995) 2765–2791. [22] Q. Kang, D. Zhang, S. Chen, Simulation of dissolution and precipitation in porous media, J. Geophys. Res. 108 (B10) (2003) 2505. [23] Q. Kang, D. Zhang, S. Chen, X. He, Lattice Boltzmann simulation of chemical dissolution in porous media, Phys. Rev. E 65 (3) (2002) 036318. [24] H. Luo, M. Quintard, G. Debenest, F. Laouafa, Properties of a diffuse interface model based on a porous medium theory for solid–liquid dissolution problems, Comput. Geosci. 16 (4) (2012) 913–932. [25] A. Parmigiani, C. Huber, O. Bachmann, B. Chopard, Pore scale mass and reactant transport in multiphase porous media flows, J. Fluid Mech. 686 (2011) 40–76. [26] L. Chen, Q. Kang, B.A. Robinson, Y.-L. He, W.-Q. Tao, Pore-scale modeling of multiphase reactive transport with phase transitions and dissolution– precipitation processes in closed systems, Phys. Rev. E 87 (4) (2013) 043306. [27] A. Büki, É. Kárpáti-Smidróczki, M. Zrínyi, Computer simulation of regular Liesegang structures, J. Chem. Phys. 103 (1995) 10387. [28] J.W. Mullin, Crystallization, Butterworth-Heinemann, 1993. [29] L. Chen, Q. Kang, Y.-L. He, W.-Q. Tao, Pore-scale simulation of coupled multiple physicochemical thermal processes in micro reactor for hydrogen production using lattice Boltzmann method, Int. J. Hydrogen Energy 37 (19) (2012) 13943–13957. [30] S.P. Sullivan, F.M. Sani, M.L. Johns, L.F. Gladden, Simulation of packed bed reactors using lattice Boltzmann methods, Chem. Eng. Sci. 60 (12) (2005) 3405–3418. [31] Q. Kang, P.C. Lichtner, D. Zhang, An improved lattice Boltzmann model for multicomponent reactive transport in porous media at the pore scale, Water Resour. Res. 43 (2007) W12S14. [32] L. Chen, H.-B. Luan, Y.-L. He, W.-Q. Tao, Pore-scale flow and mass transport in gas diffusion layer of proton exchange membrane fuel cell with interdigitated flow fields, Int. J. Therm. Sci. 51 (2012) 132–144. [33] L. Chen, Y.-L. He, Q. Kang, W.-Q. Tao, Coupled numerical approach combining finite volume and lattice Boltzmann methods for multi-scale multiphysicochemical processes, J. Comput. Phys. 255 (2013) 83–105. [34] L. Chen, H. Luan, Y. Feng, C. Song, Y.-L. He, W.-Q. Tao, Coupling between finite volume method and lattice Boltzmann method and its application to fluid flow and mass transport in proton exchange membrane fuel cell, Int. J. Heat Mass Transfer 55 (13–14) (2012) 3834–3848. [35] L. Chen, Y.-L. Feng, C.-X. Song, L. Chen, Y.-L. He, W.-Q. Tao, Multi-scale modeling of proton exchange membrane fuel cell by coupling finite volume method and lattice Boltzmann method, Int. J. Heat Mass Transfer 63 (2013) 268–283. [36] A.-J. Xie, L. Zhang, J. Zhu, Y.-H. Shen, Z. Xu, J.-M. Zhu, C.-H. Li, L. Chen, L.-B. Yang, Formation of calcium oxalate concentric precipitate rings in twodimensional agar gel systems containing Ca2+–RE3+ (RE = Er, Gd and La)– C2O2 4 , Colloids Surf., A 332 (2–3) (2009) 192–199. 496 L. Chen et al. / International Journal of Heat and Mass Transfer 75 (2014) 483–496 [37] T. Murakami, T. Kogure, H. Kadohara, T. Ohnuki, Formation of secondary minerals and its effect on anorthite dissolution, Am. Mineral. 83 (1998) 1209– 1219. [38] G. Lu, D.J. DePaolo, Qinjun Kang, D. Zhang, Lattice Boltzmann simulation of snow crystal growth in clouds, J. Geophys. Res. 114 (2009) D07305. [39] F. Xia, J. Brugger, G. Chen, Y. Ngothai, B. O’Neill, A. Putnis, A. Pring, Mechanism and kinetics of pseudomorphic mineral replacement reactions: a case study of the replacement of pentlandite by violarite, Geochim. Cosmochim. Acta 73 (7) (2009) 1945–1969. [40] M.A. Nugent, S.L. Brantley, C.G. Pantano, P.A. Maurice, The influence of natural mineral coatings on feldspar weathering, Nature 395 (1998) 588–591. [41] M.R. Lee, D.J. Brown, M.E. Mackensize, M. Smith, Weathering microenviroments on feldspar surfaces: implications for understanding fluid–mineral reactions in soils, Mineral. Mag. 71 (2008) 1319–1328. [42] H. Béarat, M.J. McKelvy, A.V.G. Chizmeshya, D. Gormley, R. Nunez, R.W. Carpenter, K. Squires, G.H. Wolf, Carbon sequestration via aqueous olivine mineral carbonation: role of passivating layer formation, Environ. Sci. Technol. 40 (15) (2006) 4802–4808. [43] H.A. Alfredsson, B.S. Hardarson, H. Franzson, S.R. Gislason, CO2 sequestration in basaltic rock at the Hellisheidi site in SW Iceland: stratigraphy and chemical composition of the rocks at the injection site, Mineral. Mag. 72 (2009) 1–5. [44] J.G. Anderson, M.A. Larson, L.K. Doraiswamy, Microphase-assisted ‘‘autocatalysis’’ in a solid–liquid reaction with a precipitating product-II. Experimental, Chem. Eng. Sci. 53 (13) (1998) 2459–2468.