- Email: [email protected]
TOUGH+CO2: A multiphase fluid-flow simulator for CO2 geologic sequestration in saline aquifers
Computers & Geosciences 37 (2011) 714–723
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
TOUGH + CO2: A multiphase fluid-flow simulator for CO2 geologic sequestration in saline aquifers Keni Zhang n, George Moridis, Karsten Pruess Earth Sciences Division, Lawrence Berkeley National Laboratory, Earth Sciences Division, 1 Cyclotron Road, Berkeley, CA 94720, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 25 January 2010 Received in revised form 9 September 2010 Accepted 15 September 2010 Available online 10 November 2010
TOUGH+CO2 is a new simulator for modeling of CO2 geologic sequestration in saline aquifers. It is a member of TOUGH+, the successor to the TOUGH2 family of codes for multicomponent, multiphase fluid and heat flow simulation. The code accounts for heat and up to 3 mass components, which are partitioned into three possible phases. In the code, the thermodynamics and thermophysical properties of H2O–NaCl–CO2 mixtures are determined based on system status and subdivided into six different phase combinations. By solving coupled mass and heat balance equations, TOUGH+CO2 can model non-isothermal or isothermal CO2 injection, phase behavior and flow of fluids and heat under typical conditions of temperature, pressure and salinity in CO2 geologic storage projects. The code takes into account effects of salt precipitation on porosity and permeability changes, and the wettability phenomena. The new simulator inherits all capabilities of TOUGH2 in handling fractured media and using unstructured meshes for complex simulation domains. The code adds additional relative permeability and capillary pressure functions. The FORTRAN 95 OOP architecture and other new language features have been extensively used to enhance memory use and computing efficiency. In addition, a domain decomposition approach has been implemented for parallel simulation. All these features lead to increased computational efficiency, and allow applicability of the code to multi-core/processor parallel computing platforms with excellent scalability. Published by Elsevier Ltd.
Keywords: CO2 geologic sequestration Saline aquifer Modeling TOUGH+ TOUGH2 Parallel computing Multiphase flow
1. Introduction CO2 geologic sequestration in saline aquifers involves complex multiphase flow processes, such as advection and diffusion, convective mixing, phase appearance/disappearance, dissolution and precipitation of minerals and other chemical reactions. Mathematical models are effective tools for understanding the behavior of CO2 in geological formations. Numerical modeling can play an important role in evaluating the feasibility and reliability of CO2 disposal. However, modeling of such a complex process as CO2 geologic sequestration in saline represents a large computational challenge. Compositional oil reservoir simulators (such as STAR-CMG, ECLIPSE) can be applied for modeling sequestration of supercritical CO2 into deep saline aquifers. Other widely used codes with CO2 simulation capability include STOMP, FEHM, DuMux, COORES, MUFTE-UG, ROCKFLOW, etc. One of the most popular numerical simulators for this type of process is the general-purpose multiphase flow simulator TOUGH2 in combination with the ECO2N fluid property module (Pruess et al., 1999; Pruess, 2005). TOUGH2/ECO2N
n
Corresponding author. E-mail address: [email protected] (K. Zhang).
0098-3004/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.cageo.2010.09.011
models non-isothermal multiphase flows of mixtures of water, salt (NaCl) and CO2. Because of the complexity of these types of subsurface flow processes, most simulations are limited to simplified systems, and mostly in two-dimensions. However, for field-scale or refined-mesh applications, a simulator with more sophisticated modeling functionalities and efficiency to represent geologic heterogeneities and the complicated multiphase, multicomponent flow processes is needed. In this study, a parallel simulator, named TOUGH+CO2, for modeling CO2 geologic sequestration in saline aquifers was developed. The simulator is a member of TOUGH+ (Moridis et al., 2006), the successor to the TOUGH2 family of codes. TOUGH+CO2 is developed based on and retains all the process-modeling capabilities of TOUGH2/ECO2N. More functions for relative permeability and capillary pressure calculation are included in the new code. Additional approaches for determining the effects of salt precipitation on porosity and permeability changes are taken into account. The code can be used to model non-isothermal multiphase flow in systems with H2O–NaCl–CO2 mixtures. It represents fluids as consisting of two mobile phases: a water-rich aqueous phase, a CO2-rich gas phase and an immobile solid halite phase. The code partitions water and CO2 between gas and aqueous phase under different temperature, pressure, and salinity conditions using the correlations developed by Spycher and Pruess (2005). Dissolution and precipitation of salt is treated by means of local equilibrium
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
solubility. Associated changes in medium porosity and permeability and the wettability phenomena may also be modeled. As special cases, the simulator may be used for investigation of sea water intrusion (no CO2) and saturated/unsaturated groundwater (no CO2 and NaCl). In this new development, the FORTRAN 95 object oriented programming architecture and new language features, such as pointers, lists and trees, data encapsulation, defined operators and assignments, operator extension and overloading, use of generic procedures and maximum use of the powerful intrinsic vector and matrix processing operations, have been extensively used. A completely different program architecture and data structure are used (compared to the TOUGH2). In addition, a domain decomposition approach has been implemented for parallel simulation. TOUGH+ CO2 adopts a similar parallelization scheme as developed for the TOUGH2-MP code by Zhang et al. (2001, 2007), and Wu et al. (2002), which has been successfully applied to many large-scale simulations (e.g. Zhang et al., 2003, 2004; Yamamoto et al., 2009). The domain decomposition approach and parallel computation enhance model simulation capabilities in terms of problem size and complexity to a level that cannot be reached by single-CPU codes. By using the TOUGH+ CO2 simulator, multi-million gridblock problems can be run on a typical Linux cluster with several tens to hundreds of processors to achieve 10–100 times improvement in computational time or problem size. The code was benchmarked through comparison of its simulation results with the corresponding results by TOUGH2/ECO2N. In this paper, the concepts, physics and governing equations of the new simulator are summarized. Design, implementation and parallelization of the code are discussed, and specific features and functionalities are introduced. The simulator is used to simulate a site-scale CO2 storage in deep saline aquifers.
2. Theories and approaches 2.1. Modeled components, phases, flow processes and equations Under typical geologic conditions for CO2 sequestration, fluids may be present in three different phases: an aqueous phase that is mostly water but may contain dissolved NaCl and CO2, a liquid CO2 phase that may contain some water and gaseous CO2 phase that also may contain some water. Liquid and gaseous CO2 may coexist along the saturated vapor pressure curve of CO2, which ends at the critical points (31.04 1C, 73.82 bar; Vargftik, 1975). If the geologic system is a saline aquifer, a third component NaCl (salt) must also be considered. Salt can be dissolved in the aqueous phase, or precipitate as solid phase. In TOUGH+ CO2, the thermophysical properties are accurately calculated for gaseous as well as for liquid CO2, but no distinction between them is made in the treatment of flow. This assumption simplifies the phase combinations in the water–NaCl–CO2 system, but it limits the code applications to only processes that do not involve phase change between liquid and gaseous sub-critical CO2. Table 1 shows the phases of a system that can be simulated.
715
The simulator was designed for simulation of the following processes and phenomena in saline aquifers: (1) the partitioning of the mass components among the possible phases; (2) the flow of gases and liquids in the geologic system; (3) precipitation and dissolution of halite, with associated changes in porosity and permeability and (4) the corresponding heat flow and transport. One of primary assumptions of the model is that Darcy’s law is assumed valid in the simulated domain. All phases may appear or disappear in any location during the course of a simulation. Thermodynamic conditions covered include a temperature range from ambient to 100 1C, pressures up to 60 MPa, and salinity from zero to full halite saturation. TOUGH+CO2 follows the mass and heat balance equations developed by Pruess et al. (1999). Integral finite difference (IFD) method (Edwards, 1972; Narasimhan and Witherspoon, 1976) is used for space discretization, while time is discretized fully implicitly as a first-order backward finite difference. One of the most important advantages of IFD is that the discretized gridblocks can have irregular shape. This feature provides flexibility for discretization of complex simulation domains, such as discrete fractures. Detailed time and space discretization schemes can be found at Pruess et al. (1999). Discretized equations are solved by the Newton/Raphson method (Pruess, 1991). In evaluation of the terms for mass balance equation, mass accumulation, flux, source and sink must be calculated at each Newton iteration step. The general form of the mass accumulation term for CO2 sequestration in saline aquifer is X Mk ¼ fSb rb Xbk ð1Þ where k is the 1,y,3 (representing water, NaCl and CO2 components), b ¼1, 2 3 (representing aqueous, gas and solid phases). Here f is the porosity, rb is the density of phase b, Sb is the saturation of phase b and Xbk is the mass fraction of component k in phase b. Before the calculation of mass accumulation, the parameters on the right hand side of Eq. (1) are calculated as a function of primary variables, which are directly solved from equations for previous Newton iteration step, or from initial conditions at first iteration. The heat accumulation term includes contributions from the rock matrix and fluid and halite phases and is given by equation X Mk ¼ ð1fÞrR CR T þ f Sb rb ub ð2Þ b
where k is the 4 (representing heat component), b the 1, 2, 3 (representing aqueous, gas phases and solid salt). Here rR and CR are the grain density and the specific heat of the host rock, respectively, T is the temperature and ub is the specific internal energy in phase b. The mass fluxes of aqueous and gaseous phases are determined by a multiphase version of Darcy’s law, written in the form b krb rb ðrPb rb gÞ Fb ¼ k0 1 þ Pb mb
ð3Þ
b ¼ 1, 2 (representing aqueous and gas phases).
Table 1 Fluid phase conditions modeled by the simulator.a Phase condition
Major constituent
Minor constituents
Single aqueous phase Single gaseous phaseb Two phases
Liquid water (1) Liquid CO2 and (2) gaseous CO2 Liquid water and CO2 (can be in liquid or gaseous phase)
Dissolved CO2 and NaCl H2O Dissolved CO2 in water; dissolved water in CO2-rich phase
a b
In all cases solid salt (NaCl) may be present or absent. All sub- and super-critical CO2 is considered as a single non-wetting phase and referred as ‘‘gaseous’’ phase.
716
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
Here k0 is absolute permeability, b is the Klinkenberg factor (1941) for gas slippage effect (b ¼0 when b ¼1), krb is the relative permeability to phase b, mb is the viscosity, Pb is the pressure in the b phase and g is the vector of gravitational acceleration. The diffusive fluxes are evaluated by the formulation provided in Pruess et al. (1999). The heat flux term accounts for conduction, advection, and radiative heat transfer, and is given by 2 3 X X k F ¼ 4ð1fÞKR þ f Sb Kb 5rT þfs s0 rT 4 þ hb Fb ð4Þ b ¼ 1,2,3
b ¼ 1,2
where KR is the thermal conductivity of the rock, Kb is the thermal conductivity of phase b, T is the temperature, hb is the specific enthalpy of phase b, fs is the radiant emittance factor and s0 is the Stefan–Boltzmann constant. In Eq. (4), rT4 is defined as T24 T14 , where T2 and T1 are temperatures at gridblocks 2 and 1 for a connection, respectively. The temperature unit in the equation is in Kelvin. TOUGH+ makes the temperature unit conversion internally for the radiation heat calculation. 2.2. Thermophysical properties and constitutive relations Through iteratively solving discretized mass balance equations, a set of primary variables will be obtained at each Newton step. The simulated brine-CO2 system is characterized and its state can be fully defined by the set of primary variables. However, to complete the mathematical description of the multiphase flow and heat transfer in porous or fractured media for CO2 geologic sequestration, the mass balance equations need to be supplemented with a number of constitutive equations. These constitutive correlations express interrelationships and constraints of physical processes, variables and parameters, and allow the evaluation of secondary variables or parameters as functions of the set of primary variables. The constitutive relations include the following: P
Constraint of fluid saturations: b Sb ¼ 1. P Constraint of mass fractions: k Xbk ¼ 1. Relation between Capillary pressure and phase saturation: PC b ¼
PC b ðSb Þ. In the code, 7 different functions have been implemented, including linear, Pickens et al. (1979), Narasimhan et al. (1978), Milly (1982), Udell and Fitch (1985), van Genuchten (1980) and Brooks and Corey (1966). Relation between relative permeability and phase saturation: krb ¼ krb ðSb Þ. The following relative permeability functions have been included: linear, Pickens et al. (1979), Corey (1954), Grant (1977), van Genuchten (1980), Verma et al. (1985) and Aziz and Settari (1979).
Other parameters, such as fluid density, viscosity, specific enthalpy of liquid and gas, thermal conductivity, and mass partitioning between components, are all determined by the primary variables. Many of these correlations for estimating properties and interrelationships are determined by experimental studies from the literature. Mass partitioning and phase change are also very important to the description of the multiphase and multicomponent flow system. The partitioning of water and CO2 among co-existing aqueous and gas phases is calculated from a slightly modified version of the corrections developed by Spycher and Pruess (2005) with Altunin’s correlation (1975) for calculation of CO2 molar volumes. The predicted equilibrium compositions for the two phases are functions of temperature, pressure and salinity. The phase equilibrium between aqueous and gas corresponds to the dissolved CO2 mass fraction in the aqueous phase. When CO2 mass
fraction is less than its equilibrium partitioning mass fraction in liquid, the system will be in single-phase liquid condition; when mass fraction of water does not exceeds its equilibrium partitioning mass fraction in gas, the system will be in single-phase gas condition. Intermediate values of mass fractions correspond to two phase conditions. 2.3. Influence of salt on flow system Precipitation and dissolution of solid salt in a porous medium changes the void space available for fluids. Such change in porosity will cause change in permeability as well. However, permeability changes depend not only on the magnitude of porosity change, but also on geometric properties of the pore channels and on where and how salt deposition in the channels occurs. The uncertainty of pore channel geometries and precipitation processes in porous media may bring difficulty for correlating between porosity and permeability change. One of the most popular conceptual models for converging–diverging pore channels is the tubes-in-series model. In TOUGH+ CO2, several choices for the functional dependence of relative change in permeability on the relative change in porosity are offered, which are developed based on the tubes-in-series model. In general, the relation can be written as k ffc ¼f ð5Þ k0 f0 fc Subscript ‘‘0’’ denotes the original property and fc is ‘‘critical porosity’’. The critical porosity is the porosity at which the throats of flow channel become clogged, so that permeability is reduced to zero. For a simplified model, the fc effect may be neglected and a power law dependence of permeability k on porosity f may be selected as the relational function. In addition, the functions for the tubes-in-series model developed by Verma and Pruess (1988) have also been adopted as an alternative choice. Solid salt deposition may not only involve permeability change, but may also influence capillary pressure. Precipitation of salt will alter the pore size distribution, generally reducing pore sizes, and thereby giving rise to stronger capillary pressures. Scaling factors of intrinsic permeability and capillary pressure are computed with original porous medium (OPM) or evolving porous medium (EPM) models developed by Moridis et al. (2008). The OPM model is based on the treatment of the medium porosity as unaffected by the emergence of salt (although subject to change due to changes in pressure and temperature), (b) of the intrinsic permeability of the porous media as unchanging during the evolution of the solid phases and (c) of the fluid flow as a relative permeability issue controlled by the saturations of the various phases in the pores. The EPM model considers the evolution of salt as tantamount to the creation of a new porous medium with continuously changing porosity and intrinsic permeability, the pore space of which is occupied only by the two fluid phases (aqueous and gas).
3. Program design and implementation 3.1. Primary variables The thermodynamic state and mass partition of different components among phase conditions are determined from the equation of state. Using the same approach as TOUGH2/ECO2N (Pruess, 2005), the multiphase flow system for CO2 sequestration is defined by four primary variables. According to Gibbs’ phase rule, the number of thermodynamic degrees of freedom for a system involving three mass components and a heat component in two (mobile) phases is four, which indicates that the four primary
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
variables can completely define the flow system state of CO2 geologic sequestration. All other (secondary) variables can be computed from the four primary variables. An important consideration in the program design is the choice of the primary variables. When a phase appears or disappears, the set of appropriate thermodynamic variables may change. In this study, a primary variable switching scheme is used. Different primary variables are chosen under different phase conditions (Table 2). Experience has proven variable-switching to be a very robust method for treating such multiphase flow system as CO2 storage in saline aquifers. The four primary variables include pressure, salt mass fraction/solid salt saturation, CO2 mass fraction and temperature for single phase; and pressure, salt mass fraction/ solid salt saturation, gas phase saturation and temperature for two phases. When no solid salt is present, the second primary variable is salt mass fraction dissolved in the two-component water–NaCl system. When solid salt is present, solid saturation is selected as the second primary variable, because the salt mass fraction is no longer an independent variable, as it is determined by the equilibrium solubility of NaCl, which is a function of temperature. In Table 2, we note some primary variables plus 10.0 are used in model calculations, because the second and third primary variables all have a value ranging from 0 to 1.0. By adding a number 10.0 the presence or absence of solid salt, and single or two phase fluid conditions can be recognized simply from the numerical value of the second and third primary variable. For the isothermal cases, heat transport is ignored and temperatures are maintained at their initial values. In TOUGH+ , a ‘‘State Identifier’’ was used for conveniently identifying the state of any location in the modeled system.
3.2. Program implementation In TOUGH+ CO2 simulation, the mass and energy balance equations are discretized in space using the integral finite difference method, and time is discretized fully implicitly as a first-order backward finite difference. The discretization results in a set of strongly coupled nonlinear algebraic equations, with the timedependent primary variables of all gridblocks as unknowns. The nonlinear equations are solved through the Newton–Raphson iteration method. The linear equations arising at each Newton step are solved by a linear solver (sequential or parallel). Primary variables are updated based on the solutions of the linear equations, and secondary variables are computed from the updated set of primary variables at each Newton iteration step. The modeling procedure is implemented with the standard FORTRAN 95/2003, except for the libraries for domain decomposition and parallel linear solver which are written in C by third parties. Compared to TOUGH2, TOUGH+ was completely rewritten and its architecture is drastically different. The code follows the approaches of Object-Oriented Programming, and involves entirely new data structures that describe the objects upon which the code is based. The basic objects are defined through derived data types, and their
717
properties and processes are described in modules and sub-modules, i.e. entities that incorporate the object attributes and parameters in addition to the procedures that describe its behavior and processes. New FORTRAN language features are used to enhance code performance, such as pointers, lists and trees, data encapsulation, defined operators and assignments, operator extension and overloading, use of generic procedures, and maximum use of the powerful intrinsic vector and matrix processing operations. The TOUGH+ code is based on a completely modular structure that is designed for maximum traceability and ease of expansion. In addition, many utility functions are carefully designed for higher efficiency, such as the use of efficient binary index searching scheme. Dynamic memory allocation for all larger arrays together with the domain decomposition approach allows large-scale problems to be solved on distributed memory computers. The code has been designed for maximum portability, scalability and efficiency. In summary, the program was written for simulations in the following steps: (1) Initialize the problem, print program information and declare variables. (2) Read and write input parameters through master processor. (3) Broadcast parameters to all processors from master processor. (4) Perform domain partitioning, and set up local distributed variable block row (DVBR) format arrays for Jacobian matrix at all processors. (5) Read and distribute mesh related data from master to all other processors. (6) Exchange domain border and external set of primary variables with neighboring processors. (7) Set up local Jacobian matrix for iterative solution of the flow equations at all processors. (8) Solve linear equations in parallel. (9) Update thermophysical parameters using current set of primary variables in all processors. (10) Converge for the Newton iteration? If not, go to step 7. (11) If the Newtonian iteration has converged, go to step 6 to proceed to new time step. (12) Output results, can be done at the end of any time steps. The above procedure is for parallel simulation. Sequential simulation is possible using the same procedure, except without multiple processors involved. 3.3. Parallelization Domain decomposition methods (DDM) are used as a divideand-conquer strategy for solving large or time-consuming problems. The idea behind this approach is to divide the computational domain into a series of subdomains. The massively parallel simulation scheme implemented into the TOUGH+CO2 simulator uses the similar approach as discussed in Zhang et al. (2001). An effective method for partitioning unstructured grids is a critical
Table 2 Phases, system states and primary variables.a Phase 1-Phase: 1-Phase: 1-Phase: 1-Phase: 2-Phase: 2-Phase:
gaseous, with salt gaseous, no salt aqueous, with salt aqueous, no salt gas + Aqu, with salt gas + Aqu, no salt
State identifier
Primary variable 1
Primary variable 2
Primary variable 3
Primary variable 4
GsS Gas AqS Aqu AGS AqG
P P P P P P
Ss + 10.0 XNaCl Ss + 10.0 XNaCl Ss + 10.0 XNaCl
XCO2 XCO2 XCO2 XCO2 Sg + 10.0 Sg + 10.0
T T T T T T
a The primary variables: P is the pressure, Ss is the solid salt saturation, XNaCl is the NaCl mass fraction in the two-component system water–NaCl, XCO2 is the CO2 mass fraction in gaseous or aqueous phase, Sg is the gas saturation and T is the temperature.
718
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
first step for a successful parallel scheme using the domaindecomposition approach. To obtain optimal parallel performance, the partitioning algorithm should take into account computational load and communication volume balance. In addition, the partition algorithm should also minimize communication volume in terms of message number and message size. It can be very difficult to achieve an ideal domain partitioning for an unstructured grid. To find an optimal selection and trade-off between these issues, we must take into account computer system characteristics, such as floating-point performance and bandwidth and latency of the communication subsystem. Because of the complexity of the problem, involving many variables, commonly used algorithms and software for partitioning large unstructured grids cannot generally optimize all these issues simultaneously. The current practice typically finds a trade-off between computation load balance and low total communication volume, even though this may not be theoretically optimal. In this study, the code partitions the simulation domain, defined by an unstructured grid, using a partitioning algorithm from the METIS software package (Karypsis and Kumar, 1998). MPI (Message Passing Interface) was used for parallel implementation. In a parallel simulation, each processor will treat one portion of the simulation domain for updating thermophysical properties, assembling mass and energy-balance equations, solving linear equation systems, and performing other local computations. Local linear-equation systems are solved in parallel by multiple processors with the Aztec linear solver package (Tuminaro et al., 1999). This parallel simulator has been built with an efficient communication scheme. Detailed discussion of the prototype of the data-exchange scheme can be found in Elmroth et al. (2001), and improvements in the communication scheme are further discussed in Zhang and Wu (2006).
3.4. Program features TOUGH+ CO2 possesses features from both the TOUGH2/ECO2N and TOUGH+ platform. Some additional useful functions have also been implemented into the code. The most important features include the following:
Dynamic memory allocation: The program allocates arrays
according to the problem size: total number of gridblocks, connections, materials (subdomains) and source/sinks. Memory requirement is shared by all processors involved in a simulation. This guarantees the scalability of the code. In addition, the code does not have limitations on the number of materials. Block-by-block permeability modification: Through this function, heterogeneous flow systems may be specified by providing gridblocks with different hydrogeologic properties through applying permeability modification (PM) coefficients to individual gridblocks. To be consistent with the series-version code, the required random numbers are generated by a single processor and are then distributed to all processors involved in the computation. Boundary conditions: The new code allows time dependent firsttype (Dirichlet) boundary conditions. Internally, it always uses large-volume approach for handling the first type boundary conditions. The large-volume approach uses a huge volume for the gridblocks with first-type boundary condition. Finite amount of mass or heat change will not cause change of thermodynamic properties at these gridblocks. Flow in fractured media: TOUGH+CO2 retains the full functions of TOUGH2 for a flexible description of flow in fractured media through either discrete-fracture, dual-porosity, dual-permeability or MINC approaches (Pruess and Narasimhan, 1982, 1985).
Upper Shale Cap
Y
30m
30m
Shale Layers
30m 184m
Z
30m 30m Horizontal Injection Well
6000m
22m X 6000m Fig. 1. Schematic representation of the 3D geometry for CO2 injection in the Utsira Formation.
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
Efficient computing schemes: In addition to the parallel computing approach, TOUGH+CO2 also adopts efficient algorithms for other computations, such as the use of overshooting technique to avoid flow phase oscillating. The code can take full advantage of the computer hardware development of multiple core/ processor systems with shared or distributed memory.
4. Application example: CO2 injection into a 3D layered saline aquifer Continuous long-term injection of CO2 into a saline aquifer will cause a buildup of groundwater pressures in extensive regions. The injection activities may induce significant hydrological and geochemical impacts on groundwater resources. One of the most common applications of TOUGH+ CO2 is for simulating such a CO2 injection process and its influence on the groundwater flow system. We use a three-dimensional (3D) example to demonstrate application of the code and its performance. The first industrial scale CO2 storage project to become operational is at the Sleipner Vest field in the Norwegian sector of the North Sea, where approximately 106 tons of CO2 per year have been injected since 1996 through a horizontal well into sands of the Utsira formation. A two-dimensional model has been designed based on this project for intercomparison of modeling codes (Pruess et al., 2003; Pruess, 2005). The model was patterned after conditions at Sleipner and was designed to investigate CO2 migration in heterogeneous sand-shale inter-layers. In this paper, a three-dimensional version of the model is developed. Even for the 3D model, significant simplifications have to be made. One of the most important simplifications is the assumption of isothermal conditions with a constant temperature of 37 1C, which is the approximate ambient temperature of the formation. CO2 injection rates, system geometry and system permeabilities correspond approximately to those at Sleipner, although no
attempt was made to accurately represent details of the permeability structure within the host formation. The injected supercritical CO2, which is less dense than the saline formation water, will rise through the formation. Its rate of ascent, however, is limited by the presence of four relatively low permeability shale layers (Fig. 1). The top and bottom of the formation is assumed to be impermeable. The only reactive chemistry considered in this problem is the dissolution of CO2 in the aqueous phase. The sequestration system is idealized as a 3D symmetric domain. The model needs only to cover one quarter of the domain because of symmetry. Fig. 1 shows the 3D model range. The horizontal injection well, which has a screen length of 100 m (50 m for the model) with a 1 m 1 m cross-section area, is along the Y-direction. The well is 30 m below the lowest shale unit, while the bottom of the aquifer is another 22 m below the well. The thickness of the formation at the site is 184 m. The injection point is 940 m below the sea floor, while the ocean depth at the site is 80 m. The formation is assumed to include four lower permeability shale units of 3 m thickness which are distributed within the high permeability sand at 30 m distance. Top and bottom of the simulation domain are impermeable boundaries. The two sides with X¼0 and Y¼0 are also impermeable boundary due to the model symmetry. No heat and mass flux is allowed across the impermeable boundaries. The other two sides with X¼6000 m and Y¼6000 m are chosen as fixed hydrostatic pressure boundaries, which allows flow into and out of the domain. The two constant pressure boundaries are around 6000 m away from the injection well, which is far enough that the CO2 cannot reach the two boundaries after 2 years of injection. The model mesh is designed with refined grid around well and near shale layers where significant gradients may occur. Gridding in the X-direction starts with 1 m increments at the well, and becomes coarser at increasing distance with total 29 gridblocks to reach 6000 m. In the Y-direction, the 50 m length horizontal well is discretized into 25 gridblocks at 2 m each. The Y-direction is discretized into a total of 50 gridblocks. Gridding in the Z-direction also uses a 1 m increment blocks at the well and coarser grid below
CO2 Gas saturation at Y = 1.0m, Time = 2 Yrs
0
Sg
-20
0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
-40 -60 Depth (m)
719
-80 -100 -120 -140 -160 -180 0
1000
2000 X (m)
Fig. 2. Distribution of gas saturation at Y¼ 1.0 m after 2 years of CO2 injection.
720
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
and above with a total 34 of gridblocks. With this discretization scheme, the model consists of 49,300 elements and 143,764 connections. The most important parameters for the model include intrinsic permeability 3 10 12 m2 for sand and 10 14 m2 for shale, corresponding porosity of 0.35 and 0.1025. Other parameters can be found in the work of Pruess et al. (2002). The initial conditions of the model are generated based on the measurements and estimations: (1) isothermal condition with T¼37 1C for whole model domain, (2) hydrostatic pressure condition with a pressure of 110 bars at the injection level, (3) a salinity of 3.2% mass fraction of NaCl, (4) CO2 in the aqueous phase in equilibrium with a PCO2 of 0.5 bars, a typical value for sedimentary formation waters at the temperature of 37 1C, which is about 0.04541% CO2 mass fraction. The initial conditions are obtained by specifying the thermodynamic properties of P¼110 bars, T¼37 1C, salinity Xs ¼0.032, CO2 mass fraction XCO2 ¼ 4:541 104 for all gridblocks with fixed pressure at the horizontal level of the injection well, and running the simulation to steady state for gravity equilibrium. The injection of 106 tons of CO2 per year to the system represents an injection rate of about 31.7 kg/s, corresponding to 0.317 kg/s for each 2 m-length injection gridblock. The total injection applied to the model is 1/4 of the total injection rate. The model is run with the given initial conditions, boundary conditions and injection rate for 2 years of continuous injection. Figs. 2–5 show some simulation results. Fig. 2 shows the distribution of gas saturation at the Y¼ 1.0 m (the vertical cross-section which is perpendicular to the horizontal injection well) after two years of simulation time. The peak gas saturations of approximately 60% occur beneath the shale layers at elevations of 123 and 99 m. When CO2 is injected into the target sand aquifer, it partially displaces the resident brine, and partially dissolves in it. Under the site subsurface temperature and pressure conditions, CO2 is buoyant (less dense) compared to brine, and the injected CO2 moves upward towards the bottom of the
shale layers. Eventually, the carbon dioxide is partially distributed among sand layers beneath the shale, and partially migrates through the low permeable shale layer into the next sand layer. At the time of 2 years, the injected CO2 only goes through one shale layer. When CO2 accumulates beneath the shale layers, it slowly dissolves in brine. Dissolution of CO2 leads to increase of the density of aqueous phase by a small amount. Although small, the density increase may be sufficient to generate convection flow in the formation. The distribution of dissolved CO2 mass fraction (Fig. 3) may confirm this process. The range of dissolved CO2 mass fraction for 4.5% or higher covers most of the two-phase zone. The plume is significantly larger than the gas phase distribution. By comparing the 3D model results with previous 2D results (Pruess et al., 2002), we can find that the 3D model predicts much smaller plumes of CO2 in gas saturation and mass fraction distribution. This suggests that a 3D model may be necessary for field-scale investigation. Fig. 4 shows a contour map of gas saturations at X ¼200 m. The cross-section of X¼200 m is parallel to the horizontal injection well, at a distance of 200 m from the well. The figure shows the front of gas CO2 distribution at the location, which has a similar pattern as the distribution at Y ¼1.0 m, except a smaller range. The figure also shows that CO2 migration along the horizontal well direction is very significant. Fig. 5 shows the three-dimensional view of gas saturation distribution after 2 years of CO2 injection. It is very clear that the injected CO2 moves upward towards the bottom of low-permeability caprock (the shale layers) because of its lower density compared to the saline water in the aquifers. The carbon dioxide is distributed among mobile layers beneath the shale layers. Highest gas saturations of approximately 60% occur beneath the first shale layer from the injection point, and the second peak saturations occur beneath the second shale layers. No gas reaches the third and top shale layers at the time. The range CO2 gas plume under the first shale layer is larger than 500 m, and about 400 m under the second shale layer.
CO2 Mass Fractions at Y = 1.0m, Time = 2 Yrs
0
XL
-20
0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0
-40
Depth (m)
-60 -80 -100 -120 -140 -160 -180 0
1000
2000 X (m)
Fig. 3. Distribution of dissolved CO2 mass fraction at Y ¼1.0 m after 2 years of CO2 injection.
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
721
CO2 Gas saturation at X = 200.0m, Time = 2 Yrs
0
Sg
-20
0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
-40
Depth (m)
-60 -80 -100 -120 -140 -160 -180 0
1000
2000 Y (m)
Fig. 4. Distribution of gas saturation at X ¼ 200.0 m after 2 years of CO2 injection.
Fig. 5. Three-dimensional view of gas saturation distribution after 2 years of CO2 injection.
Fig. 6 shows pressure changes with time at 4 gridblocks. AU111 (0.5, 1.0, 162 m) and AUP11(0.5, 49, 162 m) are gridblocks at the center and end of the horizontal injection well, respectively. Both APQ1C (200, 52, 138 m) and ASQ1C (200, 52, 156 m) are located near the end of the injection well (x¼ 200 m, y¼52 m), but APQ1C is 18 m above ASQ1C. CO2 injection causes pressures to rise initially, most strongly and rapidly in the well center gridblock. The two well gridblocks have very similar pattern of pressure evaluation curve, except the one located at the well end shows less pressure buildup. Pressure changes at APQ1C and ASQ1C are less strong and with some time delay. Once the system establishes a quasi-steady flow condition at the injection well, pressures of the well gridblocks start to decline slowly.
Correctness of the simulation results has been confirmed by comparing with the corresponding 3D model results obtained by TOUGH2/ECO2N code. Simulations for the five examples included in TOUGH2/ECO2N package also show almost identical results by the two codes. The five examples are typical problems which may be seen in solving real CO2 storage problems. The simulations were run on a Linux cluster equipped with 10 nodes using infiniBand switch connection, and each node consists of 2 INTEL Xeon 3.6 GHz CPUs. Series and parallel simulations with different number of CPUs were run for the first 0.2 years of CO2 injection to examine the code performance. The total execution time reduces from 30,504 s by series simulation, to 7112 s by
722
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
1.13E+07
AUP11 ASQ1C APQ1C AU111
1.12E+07
Pressure (Pa)
1.11E+07 1.10E+07 1.09E+07 1.08E+07 1.07E+07 1.06E+07 1.00E-03
1.00E-02
1.00E-01
1.00E+00 1.00E+01 1.00E+02 1.00E+03 Time (days)
Fig. 6. Pressure evolution with time at 4 locations.
parallel simulation with 8 CPUs and 3178 s with 16 CPUs. The code demonstrates very good performance with near linear or even super-linear speedup.
5. Conclusions An efficient numerical simulation program, named TOUGH+ CO2, for modeling of CO2 geologic sequestration in saline aquifers has been developed. The program is a member of TOUGH+ , the successor to the TOUGH2 family of codes for multicomponent, multiphase fluid and heat flow simulation. TOUGH+ CO2 was developed based on TOUGH2/ECO2N. In addition to inheriting all functionalities from the TOUGH+ platform, the code retains all the process-modeling capabilities from TOUGH2/ECO2N. The new code is a three-dimensional, fully implicit model that solves large, sparse linear systems arising from discretization of the partial differential equations for mass and energy balance in porous and fractured media. The code uses the FORTRAN 95 OOP architecture, adopts new language features, and redesigns the program data structure to enhance memory use and computing efficiency. In addition, a domain decomposition approach has been implemented for parallel simulation. All these features lead to increased computational efficiency, and allow application of the code to multi-core/processor parallel computing platforms with excellent scalability. The code was applied to simulate a CO2 sequestration project with a 3D model. The model was patterned after conditions at the Sleipner site and was designed to investigate CO2 migration in heterogeneous sand-shale inter-layers. The exercise demonstrates good performance of the new code. Reasonable simulation results have been reached.
Acknowledgements The authors would like to thank Lehua Pan and two anonymous reviewers for their review of this paper. This work was supported by funding from TOUGH royalty funds at the Earth Sciences Division, Lawrence Berkeley National laboratory. References Altunin, V.V., 1975. Thermophysical Properties of Carbon Dioxide. Publishing House of Standards, Moscow, 551 pp. (in Russian).
Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation. Elsevier, London, New York, 476 pp. Brooks, R.H., Corey, A.T., 1966. Properties of media affecting fluid flow. Journal of the Irrigation and Drainage Division 92 (IR2), 61–88. Corey, A.T., 1954. The interrelation between gas and oil relative permeabilities. Producers Monthly November, 38–41. Edwards, A.L., 1972. TRUMP: A Computer Program for Transient and Steady State Temperature Distributions in Multidimensional Systems. National Technical Information Service, National Bureau of Standards, Springfield, VA. Elmroth, E., Ding, C., Wu, Y.S., 2001. High performance computations for large-scale simulations of subsurface multiphase fluid and heat flow. Journal of Supercomputing 18 (3), 233–256. Grant, M.A.,1977. Permeability reduction factors at Wairakei. In: Proceedings of the AICHE-ASME Heat Transfer Conference, Salt Lake City, Utah, Paper 77-HT-52. Karypsis, G., Kumar, V., 1998. METIS: a software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices, V4.0. Technical Report, Department of Computer Science, University of Minnesota, 44 pp. Klinkenberg, L.J., 1941. The permeability of porous media to liquids and gases. API Drilling and Production Practice, 200–213. Milly, P.C.D., 1982. Moisture and heat transport in hysteretic, inhomogeneous porous media: a matric-head based formulation and a numerical model. Water Resources Research 18 (3), 489–498. Moridis, G., Kowalsky, M., Finsterle, S., Pruess, K., 2006. TOUGH+ : the new generation of object-oriented family of codes for the solution of problems of flow and transport in the subsurface. In: Proceedings of the TOUGH Symposium 2006, Lawrence Berkeley National Laboratory, Berkeley, California. Moridis, G.J., Kowalsky, M.B., Pruess, K., 2008. TOUGH +HYDRATE v1.0 User’s Manual. Report LBNL-0149E, Lawrence Berkeley National Laboratory, Berkeley, CA, 264 pp. Narasimhan, T.N., Witherspoon, P.A., 1976. An integrated finite difference method for analyzing fluid flow in porous media. Water Resources Research 12 (1), 57–64. Narasimhan, T.N., Witherspoon, P.A., Edwards, A.L., 1978. Numerical model for saturated–unsaturated flow in deformable porous media, Part 2: the algorithm. Water Resources Research 14 (2), 255–261. Pickens, J.F., Gillham, R.W., Cameron, D.R., 1979. Finite element analysis of the transport of water and solutes in tile-drained soils. Journal of Hydrology 40, 243–264. Pruess, K., Narasimhan, T.N., 1982. On fluid reserves and the production of superheated steam from fractured, vapor-dominated geothermal reservoirs. Journal of Geophysical Research 87 (B11), 9329–9339. Pruess, K., Narasimhan, T.N., 1985. A practical method for modeling fluid and heat flow in fractured porous media. Society of Petroleum Engineers Journal 25 (1), 14–26. Pruess, K., 1991, TOUGH2—a general-purpose numerical simulator for multiphase fluid and heart flow. Report LBNL-29400, Lawrence Berkeley Laboratory, Berkeley, CA. Pruess, K., Oldenburg, C., Moridis, G., 1999. TOUGH2 User’s Guide, Ver 2, Report LBNL-43134. Lawrence Berkeley National Laboratory, Berkeley, CA, 197 pp. Pruess, K., Garcı´a, J., Kovscek, T., Oldenburg, C., Rutqvist, J., Steefel, C., Xu, T., 2002. Intercomparison of numerical simulation codes for geologic disposal of CO2. Report LBNL-51813. Lawrence Berkeley National Laboratory, Berkeley, CA, 94720, 104 pp. Pruess, K., Bielinski, A., Ennis-King, J., Fabriol, R., Le Gallo, Y., Garcia, J., Jessen, K., Kovscek, T., Law, D.H.S., Lichtner, P., Oldenburg, C., Pawar, R., Rutqvist, J., Steefel, C., Travis, B., Tsang, C.F., White, S., Xu, T., 2003. Code intercomparison builds
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
confidence in numerical models for geologic disposal of CO2. In: Gale, J., Kaya, Y. (Eds.), GHGT-6 Conference Proceedings: Greenhouse Gas Control Technologies, Kyoto, Japan, pp. 463–470. Pruess, K., 2005. ECO2N: A TOUGH2 fluid property module for mixtures of water, NaCl, and CO2. Report LBNL-57952. Lawrence Berkeley National Laboratory, Berkeley, CA, 66 pp. Spycher, N., Pruess, K., 2005. CO2–H2O mixtures in the geological sequestration of CO2: II. Partitioning in chloride brines at 12–100 1C and up to 600 bar. Geochimica Cosmochimica Acta 69 (13), 3309–3320. doi:10.1016/j.gca. 2005.01.015. Tuminaro, R.S., Heroux, M., Hutchinson, S.A., Shadid, J.N., 1999. Official Aztec User’s Guide, Ver 2.1, Massively Parallel Computing Research Laboratory. Sandia National Laboratories, Albuquerque, NM, 63 pp. Udell, K.S., Fitch, J.S., 1985. Heat and mass transfer in capillary porous media considering evaporation, condensation, and non-condensable gas effects. In: Paper presented at the 23rd ASME/AIChE National Heat Transfer Conference, Denver, CO. van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, 892–898. Vargftik, N.B., 1975. Tables on the Thermophysical Properties of Liquids and Gases, second ed. John Wiley & Sons, New York. Verma, A.K., Pruess, K., Tsang, C.F., Witherspoon, P.A., 1985. A study of two-phase concurrent flow of steam and water in an unconsolidated porous medium. In: Proceedings of the 23rd National Heat Transfer Conference, American Society of Mechanical Engineers, Denver, CO, pp. 135–143.
723
Verma, A., Pruess, K., 1988. Thermohydrologic conditions and silica redistribution near high-level nuclear wastes emplaced in saturated geological formations. Journal of Geophysical Research 93 (B2), 1159–1173. Wu, Y.S., Zhang, K., Ding, C., Pruess, K., Elmroth, E., Bodvarsson, G.S., 2002. An efficient parallel-computing scheme for modeling nonisothermal multiphase flow and multicomponent transport in porous and fractured media. Advances in Water Resources 25, 243–261. Yamamoto, H., Zhang, K., Karasaki, K., Marui, A., Uehara, H., Nishikawa, H., 2009. Numerical investigation concerning the impact of CO2 geologic storage on regional groundwater flow. International Journal of Greenhouse Gas Control 3 (5), 586–599. Zhang, K., Wu, Y.S., Ding, C., Pruess, K., Elmroth, E., 2001. Parallel computing techniques for large-scale reservoir simulation of multi-component and multiphase fluid flow, Paper SPE 66343. In: Proceedings of the 2001 SPE Reservoir Simulation Symposium, Houston, Texas, 12 pp. Zhang, K., Wu, Y.S., Bodvarsson, G.S., 2003. Parallel computing simulation of fluid flow in the unsaturated zone of Yucca Mountain, Nevada. Journal of Contaminant Hydrology 62–63, 381–399. Zhang, K., Wu, Y.S., Bodvarsson, G.S., Liu, H.H., 2004. Flow focusing in unsaturated fracture networks: a numerical investigation. Vadose Zone Journal 3, 624–633. Zhang, K., Wu, Y.S., 2006. Enhancing scalability and efficiency of the TOUGH2-MP for linux clusters. In: Proceedings of the TOUGH Symposium 2006, Berkeley, CA. Zhang, K., Doughty, C., Wu, Y.S., Pruess, K., 2007. Efficient parallel simulation of CO2 geologic sequestration in saline aquifers, Paper SPE 106026. In: Proceedings of the 2007 SPE Reservoir Simulation Symposium, Houston, Texas, 46 pp.
Contents lists available at ScienceDirect
Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo
TOUGH + CO2: A multiphase fluid-flow simulator for CO2 geologic sequestration in saline aquifers Keni Zhang n, George Moridis, Karsten Pruess Earth Sciences Division, Lawrence Berkeley National Laboratory, Earth Sciences Division, 1 Cyclotron Road, Berkeley, CA 94720, USA
a r t i c l e i n f o
a b s t r a c t
Article history: Received 25 January 2010 Received in revised form 9 September 2010 Accepted 15 September 2010 Available online 10 November 2010
TOUGH+CO2 is a new simulator for modeling of CO2 geologic sequestration in saline aquifers. It is a member of TOUGH+, the successor to the TOUGH2 family of codes for multicomponent, multiphase fluid and heat flow simulation. The code accounts for heat and up to 3 mass components, which are partitioned into three possible phases. In the code, the thermodynamics and thermophysical properties of H2O–NaCl–CO2 mixtures are determined based on system status and subdivided into six different phase combinations. By solving coupled mass and heat balance equations, TOUGH+CO2 can model non-isothermal or isothermal CO2 injection, phase behavior and flow of fluids and heat under typical conditions of temperature, pressure and salinity in CO2 geologic storage projects. The code takes into account effects of salt precipitation on porosity and permeability changes, and the wettability phenomena. The new simulator inherits all capabilities of TOUGH2 in handling fractured media and using unstructured meshes for complex simulation domains. The code adds additional relative permeability and capillary pressure functions. The FORTRAN 95 OOP architecture and other new language features have been extensively used to enhance memory use and computing efficiency. In addition, a domain decomposition approach has been implemented for parallel simulation. All these features lead to increased computational efficiency, and allow applicability of the code to multi-core/processor parallel computing platforms with excellent scalability. Published by Elsevier Ltd.
Keywords: CO2 geologic sequestration Saline aquifer Modeling TOUGH+ TOUGH2 Parallel computing Multiphase flow
1. Introduction CO2 geologic sequestration in saline aquifers involves complex multiphase flow processes, such as advection and diffusion, convective mixing, phase appearance/disappearance, dissolution and precipitation of minerals and other chemical reactions. Mathematical models are effective tools for understanding the behavior of CO2 in geological formations. Numerical modeling can play an important role in evaluating the feasibility and reliability of CO2 disposal. However, modeling of such a complex process as CO2 geologic sequestration in saline represents a large computational challenge. Compositional oil reservoir simulators (such as STAR-CMG, ECLIPSE) can be applied for modeling sequestration of supercritical CO2 into deep saline aquifers. Other widely used codes with CO2 simulation capability include STOMP, FEHM, DuMux, COORES, MUFTE-UG, ROCKFLOW, etc. One of the most popular numerical simulators for this type of process is the general-purpose multiphase flow simulator TOUGH2 in combination with the ECO2N fluid property module (Pruess et al., 1999; Pruess, 2005). TOUGH2/ECO2N
n
Corresponding author. E-mail address: [email protected] (K. Zhang).
0098-3004/$ - see front matter Published by Elsevier Ltd. doi:10.1016/j.cageo.2010.09.011
models non-isothermal multiphase flows of mixtures of water, salt (NaCl) and CO2. Because of the complexity of these types of subsurface flow processes, most simulations are limited to simplified systems, and mostly in two-dimensions. However, for field-scale or refined-mesh applications, a simulator with more sophisticated modeling functionalities and efficiency to represent geologic heterogeneities and the complicated multiphase, multicomponent flow processes is needed. In this study, a parallel simulator, named TOUGH+CO2, for modeling CO2 geologic sequestration in saline aquifers was developed. The simulator is a member of TOUGH+ (Moridis et al., 2006), the successor to the TOUGH2 family of codes. TOUGH+CO2 is developed based on and retains all the process-modeling capabilities of TOUGH2/ECO2N. More functions for relative permeability and capillary pressure calculation are included in the new code. Additional approaches for determining the effects of salt precipitation on porosity and permeability changes are taken into account. The code can be used to model non-isothermal multiphase flow in systems with H2O–NaCl–CO2 mixtures. It represents fluids as consisting of two mobile phases: a water-rich aqueous phase, a CO2-rich gas phase and an immobile solid halite phase. The code partitions water and CO2 between gas and aqueous phase under different temperature, pressure, and salinity conditions using the correlations developed by Spycher and Pruess (2005). Dissolution and precipitation of salt is treated by means of local equilibrium
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
solubility. Associated changes in medium porosity and permeability and the wettability phenomena may also be modeled. As special cases, the simulator may be used for investigation of sea water intrusion (no CO2) and saturated/unsaturated groundwater (no CO2 and NaCl). In this new development, the FORTRAN 95 object oriented programming architecture and new language features, such as pointers, lists and trees, data encapsulation, defined operators and assignments, operator extension and overloading, use of generic procedures and maximum use of the powerful intrinsic vector and matrix processing operations, have been extensively used. A completely different program architecture and data structure are used (compared to the TOUGH2). In addition, a domain decomposition approach has been implemented for parallel simulation. TOUGH+ CO2 adopts a similar parallelization scheme as developed for the TOUGH2-MP code by Zhang et al. (2001, 2007), and Wu et al. (2002), which has been successfully applied to many large-scale simulations (e.g. Zhang et al., 2003, 2004; Yamamoto et al., 2009). The domain decomposition approach and parallel computation enhance model simulation capabilities in terms of problem size and complexity to a level that cannot be reached by single-CPU codes. By using the TOUGH+ CO2 simulator, multi-million gridblock problems can be run on a typical Linux cluster with several tens to hundreds of processors to achieve 10–100 times improvement in computational time or problem size. The code was benchmarked through comparison of its simulation results with the corresponding results by TOUGH2/ECO2N. In this paper, the concepts, physics and governing equations of the new simulator are summarized. Design, implementation and parallelization of the code are discussed, and specific features and functionalities are introduced. The simulator is used to simulate a site-scale CO2 storage in deep saline aquifers.
2. Theories and approaches 2.1. Modeled components, phases, flow processes and equations Under typical geologic conditions for CO2 sequestration, fluids may be present in three different phases: an aqueous phase that is mostly water but may contain dissolved NaCl and CO2, a liquid CO2 phase that may contain some water and gaseous CO2 phase that also may contain some water. Liquid and gaseous CO2 may coexist along the saturated vapor pressure curve of CO2, which ends at the critical points (31.04 1C, 73.82 bar; Vargftik, 1975). If the geologic system is a saline aquifer, a third component NaCl (salt) must also be considered. Salt can be dissolved in the aqueous phase, or precipitate as solid phase. In TOUGH+ CO2, the thermophysical properties are accurately calculated for gaseous as well as for liquid CO2, but no distinction between them is made in the treatment of flow. This assumption simplifies the phase combinations in the water–NaCl–CO2 system, but it limits the code applications to only processes that do not involve phase change between liquid and gaseous sub-critical CO2. Table 1 shows the phases of a system that can be simulated.
715
The simulator was designed for simulation of the following processes and phenomena in saline aquifers: (1) the partitioning of the mass components among the possible phases; (2) the flow of gases and liquids in the geologic system; (3) precipitation and dissolution of halite, with associated changes in porosity and permeability and (4) the corresponding heat flow and transport. One of primary assumptions of the model is that Darcy’s law is assumed valid in the simulated domain. All phases may appear or disappear in any location during the course of a simulation. Thermodynamic conditions covered include a temperature range from ambient to 100 1C, pressures up to 60 MPa, and salinity from zero to full halite saturation. TOUGH+CO2 follows the mass and heat balance equations developed by Pruess et al. (1999). Integral finite difference (IFD) method (Edwards, 1972; Narasimhan and Witherspoon, 1976) is used for space discretization, while time is discretized fully implicitly as a first-order backward finite difference. One of the most important advantages of IFD is that the discretized gridblocks can have irregular shape. This feature provides flexibility for discretization of complex simulation domains, such as discrete fractures. Detailed time and space discretization schemes can be found at Pruess et al. (1999). Discretized equations are solved by the Newton/Raphson method (Pruess, 1991). In evaluation of the terms for mass balance equation, mass accumulation, flux, source and sink must be calculated at each Newton iteration step. The general form of the mass accumulation term for CO2 sequestration in saline aquifer is X Mk ¼ fSb rb Xbk ð1Þ where k is the 1,y,3 (representing water, NaCl and CO2 components), b ¼1, 2 3 (representing aqueous, gas and solid phases). Here f is the porosity, rb is the density of phase b, Sb is the saturation of phase b and Xbk is the mass fraction of component k in phase b. Before the calculation of mass accumulation, the parameters on the right hand side of Eq. (1) are calculated as a function of primary variables, which are directly solved from equations for previous Newton iteration step, or from initial conditions at first iteration. The heat accumulation term includes contributions from the rock matrix and fluid and halite phases and is given by equation X Mk ¼ ð1fÞrR CR T þ f Sb rb ub ð2Þ b
where k is the 4 (representing heat component), b the 1, 2, 3 (representing aqueous, gas phases and solid salt). Here rR and CR are the grain density and the specific heat of the host rock, respectively, T is the temperature and ub is the specific internal energy in phase b. The mass fluxes of aqueous and gaseous phases are determined by a multiphase version of Darcy’s law, written in the form b krb rb ðrPb rb gÞ Fb ¼ k0 1 þ Pb mb
ð3Þ
b ¼ 1, 2 (representing aqueous and gas phases).
Table 1 Fluid phase conditions modeled by the simulator.a Phase condition
Major constituent
Minor constituents
Single aqueous phase Single gaseous phaseb Two phases
Liquid water (1) Liquid CO2 and (2) gaseous CO2 Liquid water and CO2 (can be in liquid or gaseous phase)
Dissolved CO2 and NaCl H2O Dissolved CO2 in water; dissolved water in CO2-rich phase
a b
In all cases solid salt (NaCl) may be present or absent. All sub- and super-critical CO2 is considered as a single non-wetting phase and referred as ‘‘gaseous’’ phase.
716
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
Here k0 is absolute permeability, b is the Klinkenberg factor (1941) for gas slippage effect (b ¼0 when b ¼1), krb is the relative permeability to phase b, mb is the viscosity, Pb is the pressure in the b phase and g is the vector of gravitational acceleration. The diffusive fluxes are evaluated by the formulation provided in Pruess et al. (1999). The heat flux term accounts for conduction, advection, and radiative heat transfer, and is given by 2 3 X X k F ¼ 4ð1fÞKR þ f Sb Kb 5rT þfs s0 rT 4 þ hb Fb ð4Þ b ¼ 1,2,3
b ¼ 1,2
where KR is the thermal conductivity of the rock, Kb is the thermal conductivity of phase b, T is the temperature, hb is the specific enthalpy of phase b, fs is the radiant emittance factor and s0 is the Stefan–Boltzmann constant. In Eq. (4), rT4 is defined as T24 T14 , where T2 and T1 are temperatures at gridblocks 2 and 1 for a connection, respectively. The temperature unit in the equation is in Kelvin. TOUGH+ makes the temperature unit conversion internally for the radiation heat calculation. 2.2. Thermophysical properties and constitutive relations Through iteratively solving discretized mass balance equations, a set of primary variables will be obtained at each Newton step. The simulated brine-CO2 system is characterized and its state can be fully defined by the set of primary variables. However, to complete the mathematical description of the multiphase flow and heat transfer in porous or fractured media for CO2 geologic sequestration, the mass balance equations need to be supplemented with a number of constitutive equations. These constitutive correlations express interrelationships and constraints of physical processes, variables and parameters, and allow the evaluation of secondary variables or parameters as functions of the set of primary variables. The constitutive relations include the following: P
Constraint of fluid saturations: b Sb ¼ 1. P Constraint of mass fractions: k Xbk ¼ 1. Relation between Capillary pressure and phase saturation: PC b ¼
PC b ðSb Þ. In the code, 7 different functions have been implemented, including linear, Pickens et al. (1979), Narasimhan et al. (1978), Milly (1982), Udell and Fitch (1985), van Genuchten (1980) and Brooks and Corey (1966). Relation between relative permeability and phase saturation: krb ¼ krb ðSb Þ. The following relative permeability functions have been included: linear, Pickens et al. (1979), Corey (1954), Grant (1977), van Genuchten (1980), Verma et al. (1985) and Aziz and Settari (1979).
Other parameters, such as fluid density, viscosity, specific enthalpy of liquid and gas, thermal conductivity, and mass partitioning between components, are all determined by the primary variables. Many of these correlations for estimating properties and interrelationships are determined by experimental studies from the literature. Mass partitioning and phase change are also very important to the description of the multiphase and multicomponent flow system. The partitioning of water and CO2 among co-existing aqueous and gas phases is calculated from a slightly modified version of the corrections developed by Spycher and Pruess (2005) with Altunin’s correlation (1975) for calculation of CO2 molar volumes. The predicted equilibrium compositions for the two phases are functions of temperature, pressure and salinity. The phase equilibrium between aqueous and gas corresponds to the dissolved CO2 mass fraction in the aqueous phase. When CO2 mass
fraction is less than its equilibrium partitioning mass fraction in liquid, the system will be in single-phase liquid condition; when mass fraction of water does not exceeds its equilibrium partitioning mass fraction in gas, the system will be in single-phase gas condition. Intermediate values of mass fractions correspond to two phase conditions. 2.3. Influence of salt on flow system Precipitation and dissolution of solid salt in a porous medium changes the void space available for fluids. Such change in porosity will cause change in permeability as well. However, permeability changes depend not only on the magnitude of porosity change, but also on geometric properties of the pore channels and on where and how salt deposition in the channels occurs. The uncertainty of pore channel geometries and precipitation processes in porous media may bring difficulty for correlating between porosity and permeability change. One of the most popular conceptual models for converging–diverging pore channels is the tubes-in-series model. In TOUGH+ CO2, several choices for the functional dependence of relative change in permeability on the relative change in porosity are offered, which are developed based on the tubes-in-series model. In general, the relation can be written as k ffc ¼f ð5Þ k0 f0 fc Subscript ‘‘0’’ denotes the original property and fc is ‘‘critical porosity’’. The critical porosity is the porosity at which the throats of flow channel become clogged, so that permeability is reduced to zero. For a simplified model, the fc effect may be neglected and a power law dependence of permeability k on porosity f may be selected as the relational function. In addition, the functions for the tubes-in-series model developed by Verma and Pruess (1988) have also been adopted as an alternative choice. Solid salt deposition may not only involve permeability change, but may also influence capillary pressure. Precipitation of salt will alter the pore size distribution, generally reducing pore sizes, and thereby giving rise to stronger capillary pressures. Scaling factors of intrinsic permeability and capillary pressure are computed with original porous medium (OPM) or evolving porous medium (EPM) models developed by Moridis et al. (2008). The OPM model is based on the treatment of the medium porosity as unaffected by the emergence of salt (although subject to change due to changes in pressure and temperature), (b) of the intrinsic permeability of the porous media as unchanging during the evolution of the solid phases and (c) of the fluid flow as a relative permeability issue controlled by the saturations of the various phases in the pores. The EPM model considers the evolution of salt as tantamount to the creation of a new porous medium with continuously changing porosity and intrinsic permeability, the pore space of which is occupied only by the two fluid phases (aqueous and gas).
3. Program design and implementation 3.1. Primary variables The thermodynamic state and mass partition of different components among phase conditions are determined from the equation of state. Using the same approach as TOUGH2/ECO2N (Pruess, 2005), the multiphase flow system for CO2 sequestration is defined by four primary variables. According to Gibbs’ phase rule, the number of thermodynamic degrees of freedom for a system involving three mass components and a heat component in two (mobile) phases is four, which indicates that the four primary
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
variables can completely define the flow system state of CO2 geologic sequestration. All other (secondary) variables can be computed from the four primary variables. An important consideration in the program design is the choice of the primary variables. When a phase appears or disappears, the set of appropriate thermodynamic variables may change. In this study, a primary variable switching scheme is used. Different primary variables are chosen under different phase conditions (Table 2). Experience has proven variable-switching to be a very robust method for treating such multiphase flow system as CO2 storage in saline aquifers. The four primary variables include pressure, salt mass fraction/solid salt saturation, CO2 mass fraction and temperature for single phase; and pressure, salt mass fraction/ solid salt saturation, gas phase saturation and temperature for two phases. When no solid salt is present, the second primary variable is salt mass fraction dissolved in the two-component water–NaCl system. When solid salt is present, solid saturation is selected as the second primary variable, because the salt mass fraction is no longer an independent variable, as it is determined by the equilibrium solubility of NaCl, which is a function of temperature. In Table 2, we note some primary variables plus 10.0 are used in model calculations, because the second and third primary variables all have a value ranging from 0 to 1.0. By adding a number 10.0 the presence or absence of solid salt, and single or two phase fluid conditions can be recognized simply from the numerical value of the second and third primary variable. For the isothermal cases, heat transport is ignored and temperatures are maintained at their initial values. In TOUGH+ , a ‘‘State Identifier’’ was used for conveniently identifying the state of any location in the modeled system.
3.2. Program implementation In TOUGH+ CO2 simulation, the mass and energy balance equations are discretized in space using the integral finite difference method, and time is discretized fully implicitly as a first-order backward finite difference. The discretization results in a set of strongly coupled nonlinear algebraic equations, with the timedependent primary variables of all gridblocks as unknowns. The nonlinear equations are solved through the Newton–Raphson iteration method. The linear equations arising at each Newton step are solved by a linear solver (sequential or parallel). Primary variables are updated based on the solutions of the linear equations, and secondary variables are computed from the updated set of primary variables at each Newton iteration step. The modeling procedure is implemented with the standard FORTRAN 95/2003, except for the libraries for domain decomposition and parallel linear solver which are written in C by third parties. Compared to TOUGH2, TOUGH+ was completely rewritten and its architecture is drastically different. The code follows the approaches of Object-Oriented Programming, and involves entirely new data structures that describe the objects upon which the code is based. The basic objects are defined through derived data types, and their
717
properties and processes are described in modules and sub-modules, i.e. entities that incorporate the object attributes and parameters in addition to the procedures that describe its behavior and processes. New FORTRAN language features are used to enhance code performance, such as pointers, lists and trees, data encapsulation, defined operators and assignments, operator extension and overloading, use of generic procedures, and maximum use of the powerful intrinsic vector and matrix processing operations. The TOUGH+ code is based on a completely modular structure that is designed for maximum traceability and ease of expansion. In addition, many utility functions are carefully designed for higher efficiency, such as the use of efficient binary index searching scheme. Dynamic memory allocation for all larger arrays together with the domain decomposition approach allows large-scale problems to be solved on distributed memory computers. The code has been designed for maximum portability, scalability and efficiency. In summary, the program was written for simulations in the following steps: (1) Initialize the problem, print program information and declare variables. (2) Read and write input parameters through master processor. (3) Broadcast parameters to all processors from master processor. (4) Perform domain partitioning, and set up local distributed variable block row (DVBR) format arrays for Jacobian matrix at all processors. (5) Read and distribute mesh related data from master to all other processors. (6) Exchange domain border and external set of primary variables with neighboring processors. (7) Set up local Jacobian matrix for iterative solution of the flow equations at all processors. (8) Solve linear equations in parallel. (9) Update thermophysical parameters using current set of primary variables in all processors. (10) Converge for the Newton iteration? If not, go to step 7. (11) If the Newtonian iteration has converged, go to step 6 to proceed to new time step. (12) Output results, can be done at the end of any time steps. The above procedure is for parallel simulation. Sequential simulation is possible using the same procedure, except without multiple processors involved. 3.3. Parallelization Domain decomposition methods (DDM) are used as a divideand-conquer strategy for solving large or time-consuming problems. The idea behind this approach is to divide the computational domain into a series of subdomains. The massively parallel simulation scheme implemented into the TOUGH+CO2 simulator uses the similar approach as discussed in Zhang et al. (2001). An effective method for partitioning unstructured grids is a critical
Table 2 Phases, system states and primary variables.a Phase 1-Phase: 1-Phase: 1-Phase: 1-Phase: 2-Phase: 2-Phase:
gaseous, with salt gaseous, no salt aqueous, with salt aqueous, no salt gas + Aqu, with salt gas + Aqu, no salt
State identifier
Primary variable 1
Primary variable 2
Primary variable 3
Primary variable 4
GsS Gas AqS Aqu AGS AqG
P P P P P P
Ss + 10.0 XNaCl Ss + 10.0 XNaCl Ss + 10.0 XNaCl
XCO2 XCO2 XCO2 XCO2 Sg + 10.0 Sg + 10.0
T T T T T T
a The primary variables: P is the pressure, Ss is the solid salt saturation, XNaCl is the NaCl mass fraction in the two-component system water–NaCl, XCO2 is the CO2 mass fraction in gaseous or aqueous phase, Sg is the gas saturation and T is the temperature.
718
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
first step for a successful parallel scheme using the domaindecomposition approach. To obtain optimal parallel performance, the partitioning algorithm should take into account computational load and communication volume balance. In addition, the partition algorithm should also minimize communication volume in terms of message number and message size. It can be very difficult to achieve an ideal domain partitioning for an unstructured grid. To find an optimal selection and trade-off between these issues, we must take into account computer system characteristics, such as floating-point performance and bandwidth and latency of the communication subsystem. Because of the complexity of the problem, involving many variables, commonly used algorithms and software for partitioning large unstructured grids cannot generally optimize all these issues simultaneously. The current practice typically finds a trade-off between computation load balance and low total communication volume, even though this may not be theoretically optimal. In this study, the code partitions the simulation domain, defined by an unstructured grid, using a partitioning algorithm from the METIS software package (Karypsis and Kumar, 1998). MPI (Message Passing Interface) was used for parallel implementation. In a parallel simulation, each processor will treat one portion of the simulation domain for updating thermophysical properties, assembling mass and energy-balance equations, solving linear equation systems, and performing other local computations. Local linear-equation systems are solved in parallel by multiple processors with the Aztec linear solver package (Tuminaro et al., 1999). This parallel simulator has been built with an efficient communication scheme. Detailed discussion of the prototype of the data-exchange scheme can be found in Elmroth et al. (2001), and improvements in the communication scheme are further discussed in Zhang and Wu (2006).
3.4. Program features TOUGH+ CO2 possesses features from both the TOUGH2/ECO2N and TOUGH+ platform. Some additional useful functions have also been implemented into the code. The most important features include the following:
Dynamic memory allocation: The program allocates arrays
according to the problem size: total number of gridblocks, connections, materials (subdomains) and source/sinks. Memory requirement is shared by all processors involved in a simulation. This guarantees the scalability of the code. In addition, the code does not have limitations on the number of materials. Block-by-block permeability modification: Through this function, heterogeneous flow systems may be specified by providing gridblocks with different hydrogeologic properties through applying permeability modification (PM) coefficients to individual gridblocks. To be consistent with the series-version code, the required random numbers are generated by a single processor and are then distributed to all processors involved in the computation. Boundary conditions: The new code allows time dependent firsttype (Dirichlet) boundary conditions. Internally, it always uses large-volume approach for handling the first type boundary conditions. The large-volume approach uses a huge volume for the gridblocks with first-type boundary condition. Finite amount of mass or heat change will not cause change of thermodynamic properties at these gridblocks. Flow in fractured media: TOUGH+CO2 retains the full functions of TOUGH2 for a flexible description of flow in fractured media through either discrete-fracture, dual-porosity, dual-permeability or MINC approaches (Pruess and Narasimhan, 1982, 1985).
Upper Shale Cap
Y
30m
30m
Shale Layers
30m 184m
Z
30m 30m Horizontal Injection Well
6000m
22m X 6000m Fig. 1. Schematic representation of the 3D geometry for CO2 injection in the Utsira Formation.
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
Efficient computing schemes: In addition to the parallel computing approach, TOUGH+CO2 also adopts efficient algorithms for other computations, such as the use of overshooting technique to avoid flow phase oscillating. The code can take full advantage of the computer hardware development of multiple core/ processor systems with shared or distributed memory.
4. Application example: CO2 injection into a 3D layered saline aquifer Continuous long-term injection of CO2 into a saline aquifer will cause a buildup of groundwater pressures in extensive regions. The injection activities may induce significant hydrological and geochemical impacts on groundwater resources. One of the most common applications of TOUGH+ CO2 is for simulating such a CO2 injection process and its influence on the groundwater flow system. We use a three-dimensional (3D) example to demonstrate application of the code and its performance. The first industrial scale CO2 storage project to become operational is at the Sleipner Vest field in the Norwegian sector of the North Sea, where approximately 106 tons of CO2 per year have been injected since 1996 through a horizontal well into sands of the Utsira formation. A two-dimensional model has been designed based on this project for intercomparison of modeling codes (Pruess et al., 2003; Pruess, 2005). The model was patterned after conditions at Sleipner and was designed to investigate CO2 migration in heterogeneous sand-shale inter-layers. In this paper, a three-dimensional version of the model is developed. Even for the 3D model, significant simplifications have to be made. One of the most important simplifications is the assumption of isothermal conditions with a constant temperature of 37 1C, which is the approximate ambient temperature of the formation. CO2 injection rates, system geometry and system permeabilities correspond approximately to those at Sleipner, although no
attempt was made to accurately represent details of the permeability structure within the host formation. The injected supercritical CO2, which is less dense than the saline formation water, will rise through the formation. Its rate of ascent, however, is limited by the presence of four relatively low permeability shale layers (Fig. 1). The top and bottom of the formation is assumed to be impermeable. The only reactive chemistry considered in this problem is the dissolution of CO2 in the aqueous phase. The sequestration system is idealized as a 3D symmetric domain. The model needs only to cover one quarter of the domain because of symmetry. Fig. 1 shows the 3D model range. The horizontal injection well, which has a screen length of 100 m (50 m for the model) with a 1 m 1 m cross-section area, is along the Y-direction. The well is 30 m below the lowest shale unit, while the bottom of the aquifer is another 22 m below the well. The thickness of the formation at the site is 184 m. The injection point is 940 m below the sea floor, while the ocean depth at the site is 80 m. The formation is assumed to include four lower permeability shale units of 3 m thickness which are distributed within the high permeability sand at 30 m distance. Top and bottom of the simulation domain are impermeable boundaries. The two sides with X¼0 and Y¼0 are also impermeable boundary due to the model symmetry. No heat and mass flux is allowed across the impermeable boundaries. The other two sides with X¼6000 m and Y¼6000 m are chosen as fixed hydrostatic pressure boundaries, which allows flow into and out of the domain. The two constant pressure boundaries are around 6000 m away from the injection well, which is far enough that the CO2 cannot reach the two boundaries after 2 years of injection. The model mesh is designed with refined grid around well and near shale layers where significant gradients may occur. Gridding in the X-direction starts with 1 m increments at the well, and becomes coarser at increasing distance with total 29 gridblocks to reach 6000 m. In the Y-direction, the 50 m length horizontal well is discretized into 25 gridblocks at 2 m each. The Y-direction is discretized into a total of 50 gridblocks. Gridding in the Z-direction also uses a 1 m increment blocks at the well and coarser grid below
CO2 Gas saturation at Y = 1.0m, Time = 2 Yrs
0
Sg
-20
0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
-40 -60 Depth (m)
719
-80 -100 -120 -140 -160 -180 0
1000
2000 X (m)
Fig. 2. Distribution of gas saturation at Y¼ 1.0 m after 2 years of CO2 injection.
720
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
and above with a total 34 of gridblocks. With this discretization scheme, the model consists of 49,300 elements and 143,764 connections. The most important parameters for the model include intrinsic permeability 3 10 12 m2 for sand and 10 14 m2 for shale, corresponding porosity of 0.35 and 0.1025. Other parameters can be found in the work of Pruess et al. (2002). The initial conditions of the model are generated based on the measurements and estimations: (1) isothermal condition with T¼37 1C for whole model domain, (2) hydrostatic pressure condition with a pressure of 110 bars at the injection level, (3) a salinity of 3.2% mass fraction of NaCl, (4) CO2 in the aqueous phase in equilibrium with a PCO2 of 0.5 bars, a typical value for sedimentary formation waters at the temperature of 37 1C, which is about 0.04541% CO2 mass fraction. The initial conditions are obtained by specifying the thermodynamic properties of P¼110 bars, T¼37 1C, salinity Xs ¼0.032, CO2 mass fraction XCO2 ¼ 4:541 104 for all gridblocks with fixed pressure at the horizontal level of the injection well, and running the simulation to steady state for gravity equilibrium. The injection of 106 tons of CO2 per year to the system represents an injection rate of about 31.7 kg/s, corresponding to 0.317 kg/s for each 2 m-length injection gridblock. The total injection applied to the model is 1/4 of the total injection rate. The model is run with the given initial conditions, boundary conditions and injection rate for 2 years of continuous injection. Figs. 2–5 show some simulation results. Fig. 2 shows the distribution of gas saturation at the Y¼ 1.0 m (the vertical cross-section which is perpendicular to the horizontal injection well) after two years of simulation time. The peak gas saturations of approximately 60% occur beneath the shale layers at elevations of 123 and 99 m. When CO2 is injected into the target sand aquifer, it partially displaces the resident brine, and partially dissolves in it. Under the site subsurface temperature and pressure conditions, CO2 is buoyant (less dense) compared to brine, and the injected CO2 moves upward towards the bottom of the
shale layers. Eventually, the carbon dioxide is partially distributed among sand layers beneath the shale, and partially migrates through the low permeable shale layer into the next sand layer. At the time of 2 years, the injected CO2 only goes through one shale layer. When CO2 accumulates beneath the shale layers, it slowly dissolves in brine. Dissolution of CO2 leads to increase of the density of aqueous phase by a small amount. Although small, the density increase may be sufficient to generate convection flow in the formation. The distribution of dissolved CO2 mass fraction (Fig. 3) may confirm this process. The range of dissolved CO2 mass fraction for 4.5% or higher covers most of the two-phase zone. The plume is significantly larger than the gas phase distribution. By comparing the 3D model results with previous 2D results (Pruess et al., 2002), we can find that the 3D model predicts much smaller plumes of CO2 in gas saturation and mass fraction distribution. This suggests that a 3D model may be necessary for field-scale investigation. Fig. 4 shows a contour map of gas saturations at X ¼200 m. The cross-section of X¼200 m is parallel to the horizontal injection well, at a distance of 200 m from the well. The figure shows the front of gas CO2 distribution at the location, which has a similar pattern as the distribution at Y ¼1.0 m, except a smaller range. The figure also shows that CO2 migration along the horizontal well direction is very significant. Fig. 5 shows the three-dimensional view of gas saturation distribution after 2 years of CO2 injection. It is very clear that the injected CO2 moves upward towards the bottom of low-permeability caprock (the shale layers) because of its lower density compared to the saline water in the aquifers. The carbon dioxide is distributed among mobile layers beneath the shale layers. Highest gas saturations of approximately 60% occur beneath the first shale layer from the injection point, and the second peak saturations occur beneath the second shale layers. No gas reaches the third and top shale layers at the time. The range CO2 gas plume under the first shale layer is larger than 500 m, and about 400 m under the second shale layer.
CO2 Mass Fractions at Y = 1.0m, Time = 2 Yrs
0
XL
-20
0.05 0.045 0.04 0.035 0.03 0.025 0.02 0.015 0.01 0.005 0
-40
Depth (m)
-60 -80 -100 -120 -140 -160 -180 0
1000
2000 X (m)
Fig. 3. Distribution of dissolved CO2 mass fraction at Y ¼1.0 m after 2 years of CO2 injection.
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
721
CO2 Gas saturation at X = 200.0m, Time = 2 Yrs
0
Sg
-20
0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0
-40
Depth (m)
-60 -80 -100 -120 -140 -160 -180 0
1000
2000 Y (m)
Fig. 4. Distribution of gas saturation at X ¼ 200.0 m after 2 years of CO2 injection.
Fig. 5. Three-dimensional view of gas saturation distribution after 2 years of CO2 injection.
Fig. 6 shows pressure changes with time at 4 gridblocks. AU111 (0.5, 1.0, 162 m) and AUP11(0.5, 49, 162 m) are gridblocks at the center and end of the horizontal injection well, respectively. Both APQ1C (200, 52, 138 m) and ASQ1C (200, 52, 156 m) are located near the end of the injection well (x¼ 200 m, y¼52 m), but APQ1C is 18 m above ASQ1C. CO2 injection causes pressures to rise initially, most strongly and rapidly in the well center gridblock. The two well gridblocks have very similar pattern of pressure evaluation curve, except the one located at the well end shows less pressure buildup. Pressure changes at APQ1C and ASQ1C are less strong and with some time delay. Once the system establishes a quasi-steady flow condition at the injection well, pressures of the well gridblocks start to decline slowly.
Correctness of the simulation results has been confirmed by comparing with the corresponding 3D model results obtained by TOUGH2/ECO2N code. Simulations for the five examples included in TOUGH2/ECO2N package also show almost identical results by the two codes. The five examples are typical problems which may be seen in solving real CO2 storage problems. The simulations were run on a Linux cluster equipped with 10 nodes using infiniBand switch connection, and each node consists of 2 INTEL Xeon 3.6 GHz CPUs. Series and parallel simulations with different number of CPUs were run for the first 0.2 years of CO2 injection to examine the code performance. The total execution time reduces from 30,504 s by series simulation, to 7112 s by
722
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
1.13E+07
AUP11 ASQ1C APQ1C AU111
1.12E+07
Pressure (Pa)
1.11E+07 1.10E+07 1.09E+07 1.08E+07 1.07E+07 1.06E+07 1.00E-03
1.00E-02
1.00E-01
1.00E+00 1.00E+01 1.00E+02 1.00E+03 Time (days)
Fig. 6. Pressure evolution with time at 4 locations.
parallel simulation with 8 CPUs and 3178 s with 16 CPUs. The code demonstrates very good performance with near linear or even super-linear speedup.
5. Conclusions An efficient numerical simulation program, named TOUGH+ CO2, for modeling of CO2 geologic sequestration in saline aquifers has been developed. The program is a member of TOUGH+ , the successor to the TOUGH2 family of codes for multicomponent, multiphase fluid and heat flow simulation. TOUGH+ CO2 was developed based on TOUGH2/ECO2N. In addition to inheriting all functionalities from the TOUGH+ platform, the code retains all the process-modeling capabilities from TOUGH2/ECO2N. The new code is a three-dimensional, fully implicit model that solves large, sparse linear systems arising from discretization of the partial differential equations for mass and energy balance in porous and fractured media. The code uses the FORTRAN 95 OOP architecture, adopts new language features, and redesigns the program data structure to enhance memory use and computing efficiency. In addition, a domain decomposition approach has been implemented for parallel simulation. All these features lead to increased computational efficiency, and allow application of the code to multi-core/processor parallel computing platforms with excellent scalability. The code was applied to simulate a CO2 sequestration project with a 3D model. The model was patterned after conditions at the Sleipner site and was designed to investigate CO2 migration in heterogeneous sand-shale inter-layers. The exercise demonstrates good performance of the new code. Reasonable simulation results have been reached.
Acknowledgements The authors would like to thank Lehua Pan and two anonymous reviewers for their review of this paper. This work was supported by funding from TOUGH royalty funds at the Earth Sciences Division, Lawrence Berkeley National laboratory. References Altunin, V.V., 1975. Thermophysical Properties of Carbon Dioxide. Publishing House of Standards, Moscow, 551 pp. (in Russian).
Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation. Elsevier, London, New York, 476 pp. Brooks, R.H., Corey, A.T., 1966. Properties of media affecting fluid flow. Journal of the Irrigation and Drainage Division 92 (IR2), 61–88. Corey, A.T., 1954. The interrelation between gas and oil relative permeabilities. Producers Monthly November, 38–41. Edwards, A.L., 1972. TRUMP: A Computer Program for Transient and Steady State Temperature Distributions in Multidimensional Systems. National Technical Information Service, National Bureau of Standards, Springfield, VA. Elmroth, E., Ding, C., Wu, Y.S., 2001. High performance computations for large-scale simulations of subsurface multiphase fluid and heat flow. Journal of Supercomputing 18 (3), 233–256. Grant, M.A.,1977. Permeability reduction factors at Wairakei. In: Proceedings of the AICHE-ASME Heat Transfer Conference, Salt Lake City, Utah, Paper 77-HT-52. Karypsis, G., Kumar, V., 1998. METIS: a software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices, V4.0. Technical Report, Department of Computer Science, University of Minnesota, 44 pp. Klinkenberg, L.J., 1941. The permeability of porous media to liquids and gases. API Drilling and Production Practice, 200–213. Milly, P.C.D., 1982. Moisture and heat transport in hysteretic, inhomogeneous porous media: a matric-head based formulation and a numerical model. Water Resources Research 18 (3), 489–498. Moridis, G., Kowalsky, M., Finsterle, S., Pruess, K., 2006. TOUGH+ : the new generation of object-oriented family of codes for the solution of problems of flow and transport in the subsurface. In: Proceedings of the TOUGH Symposium 2006, Lawrence Berkeley National Laboratory, Berkeley, California. Moridis, G.J., Kowalsky, M.B., Pruess, K., 2008. TOUGH +HYDRATE v1.0 User’s Manual. Report LBNL-0149E, Lawrence Berkeley National Laboratory, Berkeley, CA, 264 pp. Narasimhan, T.N., Witherspoon, P.A., 1976. An integrated finite difference method for analyzing fluid flow in porous media. Water Resources Research 12 (1), 57–64. Narasimhan, T.N., Witherspoon, P.A., Edwards, A.L., 1978. Numerical model for saturated–unsaturated flow in deformable porous media, Part 2: the algorithm. Water Resources Research 14 (2), 255–261. Pickens, J.F., Gillham, R.W., Cameron, D.R., 1979. Finite element analysis of the transport of water and solutes in tile-drained soils. Journal of Hydrology 40, 243–264. Pruess, K., Narasimhan, T.N., 1982. On fluid reserves and the production of superheated steam from fractured, vapor-dominated geothermal reservoirs. Journal of Geophysical Research 87 (B11), 9329–9339. Pruess, K., Narasimhan, T.N., 1985. A practical method for modeling fluid and heat flow in fractured porous media. Society of Petroleum Engineers Journal 25 (1), 14–26. Pruess, K., 1991, TOUGH2—a general-purpose numerical simulator for multiphase fluid and heart flow. Report LBNL-29400, Lawrence Berkeley Laboratory, Berkeley, CA. Pruess, K., Oldenburg, C., Moridis, G., 1999. TOUGH2 User’s Guide, Ver 2, Report LBNL-43134. Lawrence Berkeley National Laboratory, Berkeley, CA, 197 pp. Pruess, K., Garcı´a, J., Kovscek, T., Oldenburg, C., Rutqvist, J., Steefel, C., Xu, T., 2002. Intercomparison of numerical simulation codes for geologic disposal of CO2. Report LBNL-51813. Lawrence Berkeley National Laboratory, Berkeley, CA, 94720, 104 pp. Pruess, K., Bielinski, A., Ennis-King, J., Fabriol, R., Le Gallo, Y., Garcia, J., Jessen, K., Kovscek, T., Law, D.H.S., Lichtner, P., Oldenburg, C., Pawar, R., Rutqvist, J., Steefel, C., Travis, B., Tsang, C.F., White, S., Xu, T., 2003. Code intercomparison builds
K. Zhang et al. / Computers & Geosciences 37 (2011) 714–723
confidence in numerical models for geologic disposal of CO2. In: Gale, J., Kaya, Y. (Eds.), GHGT-6 Conference Proceedings: Greenhouse Gas Control Technologies, Kyoto, Japan, pp. 463–470. Pruess, K., 2005. ECO2N: A TOUGH2 fluid property module for mixtures of water, NaCl, and CO2. Report LBNL-57952. Lawrence Berkeley National Laboratory, Berkeley, CA, 66 pp. Spycher, N., Pruess, K., 2005. CO2–H2O mixtures in the geological sequestration of CO2: II. Partitioning in chloride brines at 12–100 1C and up to 600 bar. Geochimica Cosmochimica Acta 69 (13), 3309–3320. doi:10.1016/j.gca. 2005.01.015. Tuminaro, R.S., Heroux, M., Hutchinson, S.A., Shadid, J.N., 1999. Official Aztec User’s Guide, Ver 2.1, Massively Parallel Computing Research Laboratory. Sandia National Laboratories, Albuquerque, NM, 63 pp. Udell, K.S., Fitch, J.S., 1985. Heat and mass transfer in capillary porous media considering evaporation, condensation, and non-condensable gas effects. In: Paper presented at the 23rd ASME/AIChE National Heat Transfer Conference, Denver, CO. van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Society of America Journal 44, 892–898. Vargftik, N.B., 1975. Tables on the Thermophysical Properties of Liquids and Gases, second ed. John Wiley & Sons, New York. Verma, A.K., Pruess, K., Tsang, C.F., Witherspoon, P.A., 1985. A study of two-phase concurrent flow of steam and water in an unconsolidated porous medium. In: Proceedings of the 23rd National Heat Transfer Conference, American Society of Mechanical Engineers, Denver, CO, pp. 135–143.
723
Verma, A., Pruess, K., 1988. Thermohydrologic conditions and silica redistribution near high-level nuclear wastes emplaced in saturated geological formations. Journal of Geophysical Research 93 (B2), 1159–1173. Wu, Y.S., Zhang, K., Ding, C., Pruess, K., Elmroth, E., Bodvarsson, G.S., 2002. An efficient parallel-computing scheme for modeling nonisothermal multiphase flow and multicomponent transport in porous and fractured media. Advances in Water Resources 25, 243–261. Yamamoto, H., Zhang, K., Karasaki, K., Marui, A., Uehara, H., Nishikawa, H., 2009. Numerical investigation concerning the impact of CO2 geologic storage on regional groundwater flow. International Journal of Greenhouse Gas Control 3 (5), 586–599. Zhang, K., Wu, Y.S., Ding, C., Pruess, K., Elmroth, E., 2001. Parallel computing techniques for large-scale reservoir simulation of multi-component and multiphase fluid flow, Paper SPE 66343. In: Proceedings of the 2001 SPE Reservoir Simulation Symposium, Houston, Texas, 12 pp. Zhang, K., Wu, Y.S., Bodvarsson, G.S., 2003. Parallel computing simulation of fluid flow in the unsaturated zone of Yucca Mountain, Nevada. Journal of Contaminant Hydrology 62–63, 381–399. Zhang, K., Wu, Y.S., Bodvarsson, G.S., Liu, H.H., 2004. Flow focusing in unsaturated fracture networks: a numerical investigation. Vadose Zone Journal 3, 624–633. Zhang, K., Wu, Y.S., 2006. Enhancing scalability and efficiency of the TOUGH2-MP for linux clusters. In: Proceedings of the TOUGH Symposium 2006, Berkeley, CA. Zhang, K., Doughty, C., Wu, Y.S., Pruess, K., 2007. Efficient parallel simulation of CO2 geologic sequestration in saline aquifers, Paper SPE 106026. In: Proceedings of the 2007 SPE Reservoir Simulation Symposium, Houston, Texas, 46 pp.