Microsoft Word - 61millet.doc CHEMICAL ENGINEERING TRANSACTIONS VOL. 41, 2014 A publication of The Italian Association of Chemical Engineering www.aidic.it/cet Guest Editors: Simonetta Palmas, Michele Mascia, Annalisa Vacca Copyright © 2014, AIDIC Servizi S.r.l., ISBN 978-88-95608-32-7; ISSN 2283-9216 Three-dimensional Modelling of Reactive Solutes Transport in Porous Media Ouacil Saouli*a, Mosaab Bencheikh Lehocineb a Unité de recherche URSMA, Constantine University, Algeria b Laboratoire LIPE, Constantine University 3, Algeria saouli_ouacil2000@yahoo.fr The simulation of solutes transport in porous media is of importance in many disciplines including CO2 sequestration simulations, the evaluation of underground water quality, nuclear waste simulations, and petroleum reservoir engineering. We present a numerical method for coupling 3D transport with equilibrium reactions. Our method is based on a fixed-point algorithm that enables us to couple different transport and chemistry modules. A splitting technique is used for transport that enables us to use different spatial and temporal approximations for different transport processes. The Newton's method is used to solve the chemical system. The reliability of the method is assessed by comparing the obtained results in a column with the one-dimensional problem. 1. Introduction The simulation of reactive solutes transport in porous media is of importance in different domains: in CO2 sequestration simulations, in petroleum reservoirs, and in groundwater flow and pollution. This paper is a follow up in three-dimensional of the one-dimensional modeling described in (Saouli et al, 2011). The coupled transport and chemistry problem is solved in a three-dimensional column having uniform hydraulic conductivity under LifeV library (C++ finite element code). The flow is described by the Darcy law, the transport of solutes is described by the convection-dispersion equation, and whereas solutes chemistry involves the solution of nonlinear algebraic equations (local equilibrium is assumed). We have used operator splitting method (Samper et al, 2009), where transport and chemistry are solved separately at each time step. The outline of this paper is as follows. In Section 2, we describe the discretization of the flow equations using mixed finite element method. In Section 3, we use splitting technique that enables us to use different time approximation for convection and dispersion. Convection is approximated by finite volume method while dispersion is approximated by mixed finite element method. Section 4, presents the chemistry model in which we consider equilibrium reactions and in Section 5 the full transport-chemistry problem is modeled and solved with fixed point method. Finally, numerical results are shown in section 6. 2. Flow model For incompressible Darcy flow problem with no sink and source terms, the equations are given by = 0 = − (1) where u is the Darcy velocity, P is the pressure or the hydraulic head, and K is the hydraulic conductivity. We have used mixed finite element method. It leads to the linear system 0 = 0 (2) DOI: 10.3303/CET1441026 Please cite this article as: Saouli O., Bencheikh Lehocine M., 2014, Three-dimensional modelling of reactive solutes transport in porous media, Chemical Engineering Transactions, 41, 151-156 DOI: 10.3303/CET1441026 151 3. Transport model The advection-dispersion equation is used to model transport problem + (− ∇ + ) = (3) where c is the concentration of solute, D is dispersion tensor, and r are sources/sinks term. The dispersion tensor is given by = + ( )( . ) + (4) where Dm is molecular diffusion coefficient, αL and αT are longitudinal and transversal dispersivities respectively, τ is the tortuosity of the porous medium, and δ is the delta kroneker. In fact, molecular diffusion can be neglected in comparison with mechanical dispersion. To solve Eq(3), we have used a splitting technique. 3.1 Convection step We define cn,m an approximation of c at time tn,m. Eq(3) will come + ( ) = 0 (5) The convection step can be written as |ℎ| , ,∆ + ∑ , ∗ , = 0 (6) where ∆tc is convection time step. ∆tc must satisfy a CFL condition that in 1D takes the form = . ≤ 1 3.2 Dispersion step Once cn, M is obtained from the convection step after M time substeps. We introduce φ = -D ∇ c then the dispersion step is to find c and φ such that |ℎ| ,∆ + ∑ , = 0 (7) Given a righthandside r a transport step from cn to cn+1 is denoted in compact form as = ( , ) 4. Chemical equations In this study, we assume a local chemical equilibrium with both aqueous and sorption reactions at every point. We use the Morel formalism, where each chemical equation describes how a set of secondary species is formed from a set of primary species (Morel and Henring, 1993). The general Morel's table is Table 1: Morel's table for chemical system c s K x A1 0 Kx y Total A3 T A4 W Ky where c is the vector of concentrations of the Nc primary aqueous species, s is the vector of concentrations of the Ns sorbed species, A1, A3 and A4 are stoichiometric matrices, Kx, Ky are the equilibrium constants, x is the vector of concentrations of the Nx secondary aqueous species, y is the vector of concentrations of 152 the Ny secondary sorbed species, T is the vector of total concentrations of the aqueous species ci and W is the vector of the total concentrations of the sorbed species si. Each chemical reaction satisfies a mass action law (for each secondary aqueous and sorbed specie): log( ) = ∑ log + log ( ) , = 1, . . , log( ) = ∑ log + ∑ log + log ( ) , = 1, . . , (8) Moreover each component (primary specie) satisfies a conservation law: = + + , = 1, . . , = + ∑ , = 1, . . , (9) For a given T and W, Eq(8) and Eq(9) form a nonlinear system whose unknows are the vectors c, s, x, and y. After solving this system, we calculate = + , = (10) We will use the compact notation = ( ). 5. Coupled transport and chemistry To obtain the coupled model between transport and chemistry, we write Eq(3) for each specie. We assume that the transport operator is the same for all species. Standard manipulations lead to the coupled system + + ( ) = 0, = + (11) = ( ). where φC is the solution operator for the chemical problem defined in the previous section. After discretization in time and space, one obtains the final discrete coupled system, expressed in terms of the solution operator φT for transport introduced in section 3. = ( , , ) = + (12) = ( ) This nonlinear can be solved by various methods. The classical way is based on a fixed-point method (Amaziane et al, 2008). Suppose that Cn, Fn and Tn are given. We start by initializing Fn+1,0 = Fn, and we iterate with the step number k. The algorithm involves the following steps: 1. Compute , = ( , , , ) (using the transport solver, once for each aqueous species), 2. Compute T , = C , + F , , 3. Compute , = ( , ) (using the chemical solver, once for each degree of freedom of the mesh), 4. Set k=k+1 and go back to 1 unless the norm , − , is small enough and we stop the iterations. 153 6. Numerical results We represent here a 3D calculation design. The example of transport in presence of cation exchangers is adopted. This example is used in the documentation of PHREEQC (Parkhurst and Appelo, 1999) as Example 11. The flow and transport parameters used are presented in Table 2. Table 2: Physical parameters Length of column 0.5 m Diameter of column Hydraulic conductivity Pressure difference αL αT 0.1 m 1.16 10-7 m/s 0.5 m 0.05 m 0.005 m The chemical reactions for this example are: Na+ + X-  NaX K+ + X-  KX Ca2+ + 2 X -  CaX2 where X- is the exchanger. The initial and injected concentrations are : Table 3: Initial and injected concentrations [mol/l]. Component Cinitial Cin Ca2+ Cl- K+ Na+ 0 0 2.0 10-4 1.0 10-3 6.0 10-4 1.2 10-3 0 0 The actual chemical system is given by Table 4: Morel's table for the actual chemical system Ca2+ Na+ Cl- K+ X- K NaX 0 1 0 0 1 1 KX CaX2 Total 0 1 T1 0 0 T2 0 0 T3 1 1 5.01 0 2 2.86 103 T4 Figure 1 shows the pressure and the velocity curves. It can easily be seen that a linear distribution of pressure is obtained, when imposing the pressure difference of 0.5 m. The velocity obtained here is 0.01 m/d= 1.16 10-7 m/s which is in agreement with Darcy's law. Figure 1: Pressure and velocity distributions. 154 In order to follow the solutes concentrations along the column, Figure 2 shows the calculated concentrations of the four species. Figure 2: Concentration profiles of Ca 2+, Cl-, Na+, and K+ at time t = 25 d. These concentrations can be represented in better form (see Figure 3). Figure 3: Concentrations of solutes vs height (dm) at radius r = 0 and time t=25 d. The sorbed potassium and sodium ions are successively replaced by calcium. Because potassium exchanges more strongly than sodium (as indicated by a larger value of log K in the exchange reaction), sodium is released first, followed by potassium. Finally, when the entire concentration has been released, the concentration of calcium increases to its steady-state value, the potassium is displaced from the exchanger, and the concentration in solution increases to balance the Cl − concentration (Amir and Kern, 2010). These results are in good agreement with those obtained by (Xu et al, 1999) and computed by PhreeqC software. 7. Conclusions In this paper, the reactive solutes transport in a porous column with constant flow velocity field is examined using computational simulations. For this purpose, a calculation code was developed under LifeV environment to solve the corresponding model equations. It enables the determination of various relevant fields such as pressure, velocity and solutes concentrations and hence the dispersion tensor. The calculations were performed on 3D case and the results were compared to those obtained by PhreeqC software. 155 Finally, it can be recommended to take into account precipitation-dissolution and kinetic phenomena in the chemical model. References Amaziane B., Ossmani M.E., Serres C., 2008, Numerical modeling of the flow and transport of radionuclides in heterogeneous porous media, Compu. Geosci., 12, 437-449. Amir L., Kern M., 2010, A global method for couplind transport with chemistry in heterogeneous porous media, Comput. Geosci., 14, 465-481. Morel F.M.M., Hering J.G., 1993, Principles and applications of aquatic chemistry, Wiley, New York. Parkhurst D.L., Appelo C., 1999, User's guide to PHREEQC- a computer program for speciation, batch- reaction, one-dimensional transport, and inverse geochemical calculations, Tech. Rep., 99-4259. Saouli O., Bencheikh Lehocine M., 2011, 1-D reactive transport modeling in heterogeneous porous media, Chem. Eng. Trans., Vol. 24, 415-420. Samper J., Xu T., Yang C., 2009, A sequential partly iterative approach for multicomponent reactive transport with CORE2D, Comput. Geosci., Vol. 13, 301-316. Xu T., Samper J., Ayora C., Manzano M., Custodio E., 1999, Modeling of non-isothermal multicomponent reactive transport in field scale porous media flow system, J. Hydrol., 214, 144-164. 156