Published online 25 February 2008
Published in Vadose Zone J 7:340-349 (2008)
DOI: 10.2136/vzj2006.0076
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: TOUGH2
Efficient Schemes for Reducing Numerical Dispersion in Modeling Multiphase Transport through Heterogeneous Geological Media
Yu-Shu Wua,* and
P. A. Forsythb
a Lawrence Berkeley National Lab., Berkeley, CA 94720
b School of Computer Science, Univ. of Waterloo, Waterloo, ON, Canada
* Corresponding author (YSWu{at}lbl.gov).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 30 May 2006.
 |
ABSTRACT
|
|---|
When modeling transport of chemicals or solute in realistic large-scale subsurface systems, numerical issues are a serious concern, even with the continual progress made over the past few decades in both simulation algorithms and computer hardware. The problem becomes even more difficult when dealing with chemical transport in a vadose zone or multiphase flow system using coarse, multidimensional regular or irregular grids because of the known effects of numerical dispersion associated with moving plume fronts. We have investigated several total variation diminishing (TVD) or flux-limiter schemes by implementing and testing them in the T2R3D code, one of the TOUGH2 family of codes. The objectives of this paper are (i) to investigate the possibility of applying these TVD schemes, using multidimensional irregular unstructured grids, and (ii) to help select more accurate spatial averaging methods for simulating chemical transport, given a numerical grid or spatial discretization. We present an application example to show that such TVD schemes can effectively reduce numerical dispersion.
Abbreviations: REV, representative elementary volume • TVD, total variation diminishing.
|
INTRODUCTION
|
|---|
Numerical modeling approaches for multiphase flow and tracer or chemical transport in porous media are generally based on methodologies developed for reservoir simulation and groundwater modeling. They involve solving coupled mass-conservation equations that govern the transport processes of all chemical components in isothermal and nonisothermal subsurface systems using finite-difference or finite-element schemes. Since the 1960s, in parallel with rapid advances in multiphase flow simulation and groundwater modeling, significant progress has been made in understanding and modeling solute transport through porous and fractured media (Scheidegger, 1961; Bear, 1972; Huyakorn et al., 1983; Istok, 1989; Falta et al., 1992; Unger et al., 1996; Forsyth et al., 1998; Wu and Pruess, 2000).
Since the 1970s, transport problems involving solute and contaminant migration in porous and fractured formations have received increasing attention in the groundwater and soil science literature. As demanded by site characterization, remediation, and other environmental concerns, many quantitative modeling approaches have been developed and applied (van Genuchten and Alves, 1982; Abriola and Pinder, 1985; Corapcioglu and Baehr, 1987; Adenekan et al., 1993; Forsyth, 1994). More recently, suitability evaluation of storing high-level radioactive wastes in underground, unsaturated fractured rock has generated renewed interest in investigating tracer or radionuclide transport in a nonisothermal, multiphase fractured geological system (e.g., Viswanathan et al., 1998; Robinson et al., 2003; Moridis et al., 2003). In addition, the application of tracer tests, including environmental and manmade tracers, has become an important technique for characterizing subsurface porous-medium systems.
Even with the continual progress in both computational algorithms and computer hardware made in the past few decades, modeling the coupled processes involved in multiphase fluid flow and chemical migration in porous and fractured media remains a mathematical challenge. Many unresolved issues and limitations with current numerical approaches still exist. One concern is that severe numerical dispersion often occurs when using a multidimensional control-volume–type numerical grid in field-scale modeling studies. The problem becomes even greater when dealing with tracer transport, in which a general three-dimensional, coarse, irregular grid is used to solve advection–dispersion–type governing equations for handling tracer transport. To overcome these numerical difficulties, scientists have investigated a number of TVD or flux-limiter schemes and applied them in transport modeling, with varying success (Sweby, 1984; Liu et al., 1994; Unger et al., 1996; Forsyth et al., 1998; Oldenburg and Pruess, 1997, 2000). Many of these investigations, however, were demonstrated using one- or two-dimensional regular grids. Our work continues the effort to reduce numerical dispersion in simulating tracer or chemical plumes as they travel spatially through porous or fractured media. In this study, we examine the effectiveness of these TVD schemes as they are used in two- or three-dimensional, irregular, and unstructured grids. Our objectives are (i) to develop a general scheme for implementing different TVD schemes into multidimensional, irregular, unstructured grids of porous or fractured media; (ii) to investigate the applicability of these TVD schemes to such irregular unstructured grids; and (iii) to help select more accurate spatial averaging methods for simulating chemical transport, given a numerical grid or spatial discretization.
The TVD schemes were implemented using the T2R3D code (Wu et al., 1996), one of the TOUGH2 family of codes (Pruess, 1991). In the TOUGH2 numerical approach, a subsurface domain is discretized using an unstructured integrated finite-difference grid, followed by time discretization using a backward, first-order, finite-difference method. The resulting discrete linear or nonlinear equations are handled fully implicitly, using Newtonian iteration, with the fractured medium handled using a general multicontinuum modeling approach. Also, we present an application example to demonstrate that TVD schemes can in general effectively reduce numerical dispersion when modeling transport through heterogeneous geological fractured reservoirs.
|
Model Formulation
|
|---|
The physical processes associated with fluid flow and chemical transport in porous media are governed by fundamental conservation laws, represented mathematically (on the macroscopic level) by a set of governing equations. In addition, the movement of dissolved mass components or chemical species within a fluid in a multiphase porous-medium system is governed by advective, diffusive, and dispersive processes and is further influenced by other processes, such as radioactive decay, adsorption, dissolution and precipitation, mass exchange or partition between phases, and other chemical reactions. The formulation discussed in this section and the next is, for completeness, to model a fully coupled fluid and heat flow and tracer transport process. For a simple case of steady-state flow, the model formulation could be significantly simplified.
Governing Equation
Consider a multiphase, nonisothermal system consisting of several fluid phases, such as gas, water, and oil (nonaqueous phase liquid), with each fluid phase in turn consisting of a number of mass components. To derive a set of generalized governing equations for multiphase fluid flow, multicomponent transport, and heat transfer, we assume that these processes can be described using a continuum approach within a representative elementary volume (REV) in a porous or fractured medium (Bear, 1972). In addition, a condition of local thermodynamic equilibrium is assumed so that at any time, temperatures, phase pressures, densities, viscosities, enthalpies, internal energies, and component concentrations (or mass fractions) are the same locally within each REV of the porous medium.
According to mass- and energy-conservation principles, a generalized conservation equation of mass components and energy in the porous continuum can be written as follows:
 | [1] |
where superscript k is the index for the components, k = 1, 2, 3,..., Nc, with Nc being the total number of mass components and with k = Nc + 1 for the energy "component" (thermal energy is treated as a component for convenience); Mk is the accumulation term of component k; Gk is the decay or internal generation (reaction) term of mass or energy components; qk is an external source–sink term or fracture–matrix exchange term for mass or energy component k and energy; and Fk is the "flow" term, describing mass or energy movement, contributed by multiphase flow, or diffusive and dispersive mass transport, or heat transfer, as discussed and defined below.
Under equilibrium adsorption, the accumulation term Eq. [1] for component k is
 | [2] |
where 
is the porosity of porous media; β is an index for fluid phase (β = g for gas, w for aqueous phase); 
β is the density of phase β; Sβ is the saturation of phase β; Xβ
is the mass fraction of component k in fluid β; s is the density of rock solids; and Kdk is the distribution coefficient of component k between the aqueous phase and rock solids to account for adsorption effects. In the case in which components are subject to a first-order radioactive decay, the decay generation term is
 | [3] |
where 
k is the radioactive decay constant of component k. The generation term may also be subject to other processes such as dissolution and precipitation, mass exchange and partition between phases, or chemical reactions under equilibrium or nonequilibrium condition.
The accumulation term in Eq. [1] for the heat equation is usually defined as
 | [4] |
where Uβ and Us are the internal energies of fluid β and rock solids, respectively.
The mass component transport is governed in general by processes of advection, diffusion, and dispersion. Advective transport of a component or solute is carried by fluid flow, and diffusive and dispersive flux is contributed by molecular diffusion and mechanical dispersion, or hydrodynamic dispersion. These processes are described using a modified Fick's law for the total mass flow term for a component k, by advection and dispersion, written as
 | [5] |
where vβ is a vector of the Darcy's velocity or volumetric flow, defined by Darcy's law to describe the flow of single or multiple immiscible fluids as
 | [6] |
where Pβ, µβ, and g are pressure, viscosity of fluid phase β, and gravitational acceleration constant, respectively; z is the vertical coordinate; k is absolute or intrinsic permeability; krβ is the relative permeability to phase β, and g is gravitational constant (vector). In Eq. [5], Dβk is the hydrodynamic dispersion tensor accounting for both molecular diffusion and mechanical dispersion for component k in phase β, defined by an extended dispersion model (Scheidegger, 1961; Bear, 1972) to include multiphase effects as
 | [7] |
where 
Tβ and Lβ are transverse and longitudinal dispersivities, respectively, in fluid β of porous media; 
is tortuosity of the porous medium; dβk is the molecular diffusion coefficient of component k within fluid β; and 
ij is the Kronecker delta function (ij = 1 for i = j, and ij = 0 for i 
j), with i and j being coordinate indices.
Equation [5] indicates that the mass flow consists of two parts. The first part, the first term on the right-hand side of Eq. [5], is contributed by advection in all phases, and the second part, the second term on the right-hand side of Eq. [5], is diffusive flux by hydrodynamic dispersion.
Heat transfer in porous media is in general a result of both convective and conductive processes, although in certain cases, radiation may also be involved. These heat-transfer processes are complicated by interactions between multiphase fluids, multicomponents, and associated changes in phases, internal energy, and enthalpy. Heat convection is contributed by thermal energy carried mainly by bulk flow of all fluids as well as by dispersive mass fluxes. On the other hand, heat conduction or radiation is driven by temperature gradients and may follow Fourier's law or Stefan–Boltzmann's law, respectively. Then, the combined, overall heat flux term, due to convection, conduction, and radiation in a multiphase, multicomponent, porous medium system, may be described as
 | [8] |
where hβ and hβk are specific enthalpies of fluid phase β and of component k in fluid β, respectively; KT is the overall thermal conductivity; T is temperature; 
o is a radiation emissivity factor, and 
o (= 5.6687 x 10–8 J m–2 K–4) is the Stefan–Boltzmann constant.
As shown in Eq. [8], the total heat flow in a multiphase, multicomponent system is determined by heat convection of flow and mass dispersion (the first two terms on the right-hand side of Eq. [8]), heat conduction (the third term on the right-hand side), and thermal radiation (the last term on the right-hand side).
Constitutive Relationships
To complete the mathematical description of multiphase flow, multicomponent transport, and heat transfer in porous media, Eq. [1], a generalized mass- and energy-balance equation, needs 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 and parameters as functions of a set of primary unknowns or variables selected to make the governing equations solvable. Table 1
lists a commonly used set of constitutive relationships for describing multiphase flow, multicomponent mass transport, and heat transfer through porous media. Many of these correlations for estimating properties and interrelationships are determined by experimental studies.
|
Numerical Formulation
|
|---|
The methodology for using numerical approaches to simulate multiphase subsurface flow and transport consists in general of the following three steps: (i) spatial discretization of mass and energy-conservation equations of Eq. [1], (ii) time discretization, and (iii) iterative approaches to solve the resulting nonlinear, discrete algebraic equations. Among various numerical techniques for simulation studies, a mass- and energy-conserving discretization scheme, based on finite volume or integral finite-difference or finite-element methods, is the most commonly used approach, and the one discussed here.
Discrete Equations
The component mass and energy balance, as expressed in Eq. [1], is discretized in space using a control-volume, integrated finite-difference concept (Narasimhan and Witherspoon, 1976; Pruess, 1991). The control-volume approach provides a general spatial discretization scheme that can represent a one-, two-, or three-dimensional domain using a set of discrete meshes. Each mesh has a certain control volume for a proper averaging or interpolation of flow and transport properties or thermodynamic variables. Time discretization is performed using a backward, first-order, fully implicit finite-difference scheme. The discrete nonlinear equations for components in the multiphase system at gridblock or node i can be written in a general form:
 | [9] |
where the superscript k also serves as an equation index for all mass components with k = 1, 2, 3, ..., Nc and k = Nc + 1 denotes the heat equation; superscript n denotes the previous time level, with n + 1 the current time level to be solved; subscript i refers to the index of gridblock or node i, with N being the total number of nodes in the grid; 
t is time step size; Vi is the volume of node i; 
i contains the set of direct neighboring nodes (j) of node i; and Aik and Gik are the accumulation and decay generation terms, respectively, at node i. The "flow" term between nodes i and j, flowijk, and the sink–source term at node i for component k or thermal energy, Qik, are defined below.
Equation [9] has the same form regardless of the dimensionality of the system; that is, it applies to one-, two-, or three-dimensional flow, transport, and heat-transfer analyses. The accumulation and decay generation terms for mass components or thermal energy are evaluated using Eq. [2], [3], and [4], respectively, at each node i. The "flow" terms in Eq. [9] are generic and include mass fluxes by advective and dispersive processes, as described by Eq. [5], as well as heat transfer, described by Eq. [8]. In general, the mass flow term is evaluated as (Wu and Pruess, 2000)
 | [10] |
where FA,ijk and FD,ijk are the net mass fluxes by advection and hydrodynamic dispersion along the connection, respectively, with
 | [11] |
where Aij is the common interface area between connected blocks or nodes i and j; and the mass flow term, Fβ,ij, of fluid phase β is described by a discrete version of Darcy's law. That is, the mass flux of fluid phase β along the connection is given by
 | [12] |
where 
ij is transmissivity and is defined differently for finite-difference or finite-element discretization, and 
β is the flow potential term at node j or i (Pa). If the integral finite-difference scheme (Pruess, 1991) is used, the transmissivity is calculated as
 | [13] |
where Di is the distance from the center of block i to the interface between blocks i and j. The flow potential term in Eq. [12] is defined as
 | [14] |
where Zi is the depth to the center of block i from a reference datum.
For mass component transport, the flow term or the net mass flux by advection and hydrodynamic dispersion of a component along the connection of nodes i and j is determined by
 | [15] |
where nij is the unit vector along the connection of the two blocks i and j.
The total heat flux along the connection of nodes i and j, including advective, diffusive, conductive, and radiation terms, may be evaluated, when using a finite-difference scheme, by
 | [16] |
In evaluating the "flow" terms in Eq. [11–14] and Eq. [16], subscript ij + 1/2 is used to denote a proper averaging or weighting of fluid flow, component transport, or heat-transfer properties at the interface or along the connection between two blocks or nodes i and j. The proper weighting of mass fraction in Eq. [11] for calculating advective mass flux is the objective of this work, which is discussed in detail below. The convention for the signs of flow terms is that flow from node j into node i is defined as "+" (positive) in calculating the flow terms. Wu and Pruess (2000) presented a general approach to calculating these flow terms associated with advective and dispersive mass transport and heat transfer in a multiphase system, using an irregular and unstructured multidimensional grid.
Spatial and Temporal Weighting and Flux Limiter Schemes
As shown in Eq. [11] and [12], in general, two types of spatial weighting schemes are needed in modeling multiphase tracer transport. The first one, (Xβk)ij+1/2), is used in Eq. [11] for estimating the averaged mass fraction for calculating advective flux. The other, (βkrβ/µβ)ij+1/2, is used in Eq. [12] for mobility weighting of the multiphase flow term. In the literature, flux-limiter schemes have been used not only for the first type of weighting but also for the second type (Blunt and Rubin, 1992; Oldenburg and Pruess, 2000). However, it has been observed in practical simulations (Datta-Gupta et al., 1991) that the numerical smearing caused by saturation fronts is in general much less severe than that with dissolved concentration fronts. Therefore, in this work, we focus our attention on the use of the flux limiters for mass fraction averaging or modeling concentration plume only, whereas the traditional, full upstream weighting is used in mobility or relative permeability averaging for estimating fluid displacement or saturation fronts.
In addition to spatial weighting schemes, temporal weighting also needs to be addressed in numerical formulation. Commonly used temporal weighting schemes include fully implicit and Crank–Nicolson methods. The fully explicit weighting is rarely used because of its strict limitation in time-step size. Among these schemes, the fully implicit method has proven itself to be most effective in handling numerical problems associated with solving highly nonlinear multiphase flow equations. In particular, the theoretical analysis of advective–dispersive transport through a one-dimensional finite-volume grid by Unger et al. (1996) indicates that the fully implicit scheme has no limitations in Courant number under various temporal weighting schemes, including flux limiters. They demonstrate how fully implicit temporal weighting leads to unconditionally stable solutions for linear advection–dispersion equations. Although fully implicit weighting is only a first-order approximation, with numerical errors of the same size as the time step, it is our experience (in conducting hundreds of large, field-scale simulations of coupled multiphase flow and chemical transport) that fully implicit temporal schemes always result in stable solutions and that temporal discretization errors, caused by a fully implicit scheme, are of secondary importance in simulation when compared with many other unknowns. The key is to have a robust numerical scheme that leads to reliable and stable solutions under different spatial discretizations and various physical conditions. Because it is impractical to define a Peclet or Courant number for detailed theoretical analyses in most field applications when using multidimensional, irregular, unstructured grids, fully implicit temporal weighting should be a first choice.
Selecting proper spatial-weighting schemes becomes critical when dealing with coupled processes of multiphase flow, chemical transport, and heat transfer in a fractured medium. This is because fracture and matrix characteristics often greatly differ; for example, there can be a many-orders-of-magnitude contrast in flow and transport properties, such as permeability and dispersivity. It is further complicated in that there are no generally applicable weighting schemes or rules applicable to all problems or processes (Wu and Pruess, 2000). The following weighting schemes were used for flux calculation in this work:
- Upstream weighting for relative permeability and/or mobility
- Harmonic or upstream weighting for absolute permeabilities in global fracture or matrix flow
- Matrix absolute permeability, thermal conductivity, and molecular diffusion coefficients for fracture–matrix interaction
- Phase saturation–based weighting functions for determining diffusion coefficients
- Upstream weighted enthalpies for advective heat flow
- Central weighted scheme for thermal conductivities of global heat conduction
Consider the schematic of Fig. 1
, representing a multidimensional irregular, unstructured grid of porous and/or fractured media. To calculate advective flux between nodes i and j, we also need the information from a secondary upstream node (denoted as i2up), which is an upstream node to the upstream one, ups(i,j), between nodes i and j (Unger et al., 1996; Forsyth et al., 1998). As shown in Fig. 1, the node i2up is determined by the maximum potential method in terms of maximum fluid influx into the node ups(i,j), which is implemented at Newtonian iteration level for every connection. Various weighting schemes for spatially averaged mass fraction or concentration for advective flux calculation between nodes i and j are summarized as upstream and central.
View larger version (26K):
[in this window]
[in a new window]
|
FIG. 1. Schematic for determining the second upstream block (i2up) for flow between block i and block j, using the geometric method and the maximum potential method.
|
|
In the upstream weighting scheme,
 | [17] |
where subscript ups(i,j) stands for the upstream node for fluid flow between nodes i and j.
For the central weighting scheme,
 | [18] |
The several flux-limiter or TVD schemes tested include the van Leer limiter scheme, the MUSCL method, and the Leonard method. The van Leer flux limiter scheme is expressed as
 | [19] |
where (Xβk)dwn(i,j) is the mass fraction of downstream node of i and j, defined as
 | [20] |
The van Leer weighting factor (rij) is defined as
 | [21] |
with the smoothness sensor,
 | [22] |
where Dups(i,j) and Di2up are the distances from the center of block ups(i,j) or its upstream block i2up to their common interface along the connection between the blocks.
The MUSCL method is described by
 | [23] |
where
 | [24] |
 | [25] |
 | [26] |
and
 | [27] |
In Eq. [24], is a small number, which prevents a zero divide.
The Leonard flux limiter is also described by Eq. [19] with the Leonard weighting factor (rij) defined as
 | [28] |
with rij also defined by Eq. [22].
The numerical implementation of these TVD schemes is made in the T2R3D code (Wu et al., 1996) for simulation of tracer transport through an isothermal or nonisothermal system.
Numerical Solution Scheme
Over the past few decades, a number of numerical solution techniques have been developed in the literature to solve the nonlinear, discrete equations of flow and transport in porous media. When handling multiphase flow, multicomponent transport, and heat transfer in a multiphase flow system, investigators predominantly use a fully implicit scheme. This is because of the extremely high nonlinearity inherent in those discrete equations and the many numerical schemes with different levels of explicitness that may fail to converge in practice. In this section, we discuss a general procedure to solve the discrete nonlinear Eq. [9] fully implicitly, using a Newton iteration method.
We can write the discrete nonlinear Eq. [9] in a residual form as
 | [29] |
Equation [29] defines a set of (Nc + 1) x N coupled nonlinear equations that need to be solved for every balance equation of mass components and heat, respectively. In general, (Nc + 1) primary variables per node are needed to use the Newton iteration for the associated (Nc + 1) equations per node. The primary variables are usually selected among fluid pressures, fluid saturations, mass (mole) fractions of components in fluids, and temperatures. In many applications, however, primary variables cannot be fixed and must be allowed to vary dynamically to account for phase appearance and disappearance (Forsyth et al., 1998). The rest of the dependent variables are treated as secondary variables, including relative permeability, capillary pressures, viscosity and densities, partitioning coefficients, specific enthalpies, thermal conductivities, and dispersion tensor (listed in Table 1), as well as nonselected pressures, saturations, and mass (mole) fractions.
In terms of the primary variables, the residual Eq. [29] at a node i is regarded as a function of the primary variables at not only node i but also at all its direct neighboring nodes j. The Newton iteration scheme gives rise to
 | [30] |
where xm is the primary variable m, with m = 1, 2, 3, ..., Nc + 1, respectively, at node i and all its direct neighbors; p is the iteration level; and i = 1, 2, 3, ..., N. The primary variables in Eq. [30] need to be updated after each iteration:
 | [31] |
The Newton iteration process continues until the residuals Rnk,n+1 or changes in the primary variables xm,p+1 over an iteration are reduced below preset convergence tolerances.
A numerical method is used to construct the Jacobian matrix for Eq. [30], as outlined in Forsyth et al. (1995). At each Newton iteration, Eq. [30] represents a system of (Nc + 1) x N linearized algebraic equations with sparse matrices, which are solved by a linear equation solver. When using the flux limiter schemes, as discussed in the last subsection, advective mass flux terms in the discrete Eq. [21] may depend on primary and secondary variables beyond the direct neighboring nodes, such as at the node i2up. In such a situation, the Newton iteration discussed here becomes inexact because the Jacobian matrix does not include the contributions with respect to the primary variables beyond neighboring nodes. Nevertheless, converged solutions should be correct because the residuals are exact. This omission in these Jacobian calculations may make solution convergence more problematic. Many numerical tests, however, have been made for multiphase tracer transport, and no significant numerical problems have been observed.
Fractured Media
The mathematical formulations and flux-limiter schemes discussed above are applicable to both single-continuum and multicontinuum media, as long as the physical processes involved can be described in a continuum sense within either continuum. When handling flow and transport through a fractured rock using the numerical formation of this section, fractured media (including explicit fracture, dual, or multiple continuum models) can be considered special cases of unstructured grids (See Fig. 1). Then, a large portion of the work consists of generating a mesh that represents both the fracture and the matrix system under consideration. Several fracture and matrix subgridding schemes exist for designing different meshes for different fracture–matrix conceptual models (e.g., Pruess, 1983).
Once a proper unstructured grid of a fracture–matrix system is generated, fracture and matrix blocks are identified to represent fracture and matrix domains, separately. Formally they are treated identically for the solution in the model. However, physically consistent fracture and matrix properties, parameter weighting schemes, and modeling conditions must be appropriately specified for both fracture and matrix systems.
|
Application
|
|---|
One example is presented here to demonstrate the application of the TVD schemes discussed above in handling transport through fractured media. The sample problem is based on a two-dimensional site-scale model developed for investigations of the unsaturated zone at Yucca Mountain in Nevada. This example shows transport of one conservative (nonadsorbing) tracer through unsaturated fractured rock using a two-dimensional unstructured grid and a dual-permeability conceptualization for handling fracture and matrix interaction.
The two-dimensional west–east cross-sectional model grid, shown in Fig. 2
, has a total of 30,000 fracture-matrix gridblocks and 74,000 connections between them in a dual-permeability mesh. The potential repository is located in the middle of the model domain, discretized with a locally refined grid (Fig. 2), at an elevation above 1100 m.
View larger version (59K):
[in this window]
[in a new window]
|
FIG. 2. Two-dimensional west–east cross-sectional model domain and grid showing lateral and vertical discretization, hydrogeological layers, repository layout, and several faults incorporated. The cross-section is located between two Nevada coordinates of (169, 372 m, 233, 052 m) and (171, 844 m, 233, 858 m).
|
|
The two-dimensional model uses the ground surface as the top model boundary and the water table as the bottom boundary. Both top and bottom boundaries of the model are treated as Dirichlet-type boundaries; that is, constant (spatially distributed) pressures, liquid saturations, and zero initial tracer concentrations are specified along these boundary surfaces. In addition, on the top boundary, a spatially varying, steady-state, present-day infiltration map (as shown in Fig. 3
), determined by the scientists of the USGS, is used in this study to describe the net water recharge, with an average infiltration rate of about 5 mm yr–1 over the model domain. In addition, an isothermal condition is assumed in this study. The properties used for rock matrix and fractures in the dual-permeability model, including two-phase flow parameters of fractures and matrix as well as faults, were estimated from field tests and model calibration efforts (Wu et al., 2002).
View larger version (18K):
[in this window]
[in a new window]
|
FIG. 3. Net infiltration rate along the west–east cross-sectional model as surface water recharge boundary condition.
|
|
We consider a conservative liquid tracer migrating from the repository downward by advective and dispersive processes, subject to the ambient steady-state unsaturated flow condition. A constant effective molecular diffusion coefficient of 3.2 x 10–11 (m2 s–1) is used for matrix diffusion of the conservative component. Dispersivities used are for matrix: = 5 m and T,M = 0.5 m; for fractures: L,F = 10 m and T,F = 1 m. Transport starts with a finite amount of the tracer initially released into the fracture elements of the repository blocks. After the simulation starts, no more tracer will be introduced into the system, but the steady-state water recharge from the top boundary continues. Eventually, all the tracer will be flushed out from the two-dimensional system through the bottom, the water table boundary, by advective and diffusive processes.
Figures 4
and 5
show normalized tracer concentration contours in the fracture continuum within the two-dimensional model at 10 yr of tracer release, simulated using various weighting schemes of spatially averaged mass fraction for advective flux calculation. Comparisons of simulated concentrations between Fig. 4 (central weighting) and Fig. 5 (TVD–MUCSL) show a large difference at the time of 10 yr. Because all three TVD schemes implemented in this study give similar results for this problem, only the results with MUSCL are shown for the TVD cases in Fig. 5. Figure 6
presents fractional cumulative mass breakthrough curves at the water table, also showing some significant differences between the results using the TVD schemes and the central weighting.
View larger version (29K):
[in this window]
[in a new window]
|
FIG. 4. Concentration (C) distributions within the two-dimensional model at 10 yr, simulated using the central weighting scheme.
|
|
View larger version (28K):
[in this window]
[in a new window]
|
FIG. 5. Concentration (C) distributions within the two-dimensional model at 10 yr, simulated using the total variation diminishing (TVD) scheme (MUSCL).
|
|
View larger version (20K):
[in this window]
[in a new window]
|
FIG. 6. Breakthrough curves of fractional cumulative tracer mass arriving at the water table, since release from the repository, simulated using the different weighting schemes.
|
|
Overall, the simulation results indicate that at early time periods, such as in the first 10 yr (Fig. 4), the central weighting scheme underestimates advective transport, while at later time periods (t > 100 yr) it overestimates advective transport, because of the selection of too high or too low averaged concentration values. In comparison, the TVD schemes provide more accurate, consistent numerical solutions at all the times. In addition, the TVD schemes are tested and found to have much better numerical performance than the central weighting scheme with respect to taking larger time steps or stability.
|
Summary and Conclusions
|
|---|
We have investigated several TVD schemes by implementing them into the TOUGH2 family of codes, using multidimensional irregular unstructured grids. Our test results show that such TVD schemes can reduce numerical dispersion effectively, if used properly. In addition, numerical performance with TVD schemes is significantly better than commonly used central weighting and is comparable to fully upstream weighting. It is encouraging that under multiphase conditions using relatively coarse spatial discretization, these TVD schemes provide more accurate simulation results for modeling large-scale field tracer transport processes through heterogeneous, fractured rock than the traditional modeling approaches.
|
Appendix: Nomenclature
|
|---|
|
ACKNOWLEDGMENTS
|
|---|
The authors thank Lehau Pan for his help with the model grid. Thanks are also due to Guoxiang Zhang, Curt Oldenburg, and Dan Hawkes for their review of the paper. This work was supported in part by the U.S. Department of Energy. The support is provided to Berkeley Laboratory through the U.S. Department of Energy Contract No. DE-AC03-76SF00098.
|
REFERENCES
|
|---|
- Abriola, L.M., and G.F. Pinder. 1985. A multiphase approach to the modeling of porous media contamination by organic compounds: 1. Equation development. Water Resour. Res. 21(1):11–18.[CrossRef]
- Adenekan, A.E., T.W. Patzek, and K. Pruess. 1993. Modeling of multiphase transport of multicomponent organic contaminants and heat in the subsurface: Numerical model formulation. Water Resour. Res. 29(11):3727–3740.[CrossRef]
- Bear, J. 1972. Dynamics of fluids in porous media. American Elsevier, New York.
- Blunt, M., and B. Rubin. 1992. Implicit flux limiting schemes for petroleum reservoir simulation. J. Comp. Phys. 102:194–210.[CrossRef]
- Corapcioglu, M.Y., and A.L. Baehr. 1987. A compositional multiphase model for groundwater contamination by petroleum products: 1. Theoretical considerations. Water Resour. Res. 23(1):191–200.[CrossRef]
- Datta-Gupta, A., L.W. Lake, G.A. Pope, K. Sepehrnoori, and M.J. King. 1991. High resolution monotonic schemes for reservoir fluid flow simulation. In Situ 15(3):289–317.[Web of Science]
- Falta, W.R., K. Pruess, I. Javandel, and P.A. Witherspoon. 1992. Numerical modeling of steam injection for the removal of nonaqueous phase liquids from the subsurface: 1. Numerical formulation. Water Resour. Res. 28(2):433–449.[CrossRef]
- Forsyth, P.A. 1994. Three-dimensional modeling of steam flush for DNAPL site remediation. Int. J. Numer. Methods Fluids 19:1055–1081.[CrossRef]
- Forsyth, P.A., Y.S. Wu, and K. Pruess. 1995. Robust numerical methods for saturated–unsaturated flow with dry initial conditions in heterogeneous media. Adv. Water Resour. 18:25–38.[CrossRef]
- Forsyth, P.A., A.J.A. Unger, and E.A. Sudicky. 1998. Nonlinear iteration methods for nonequilibrium multiphase subsurface flow. Adv. Water Resour. 21:433–499.[CrossRef]
- Huyakorn, P.S., B.H. Lester, and J.W. Mercer. 1983. An efficient finite element technique for modeling transport of fractured porous media: 1. Single species transport. Water Resour. Res. 19(3):841–854.[CrossRef]
- Istok, J. 1989. Groundwater modeling by the finite element method. Water Resources Monogr. 13. American Geophysical Union, Washington, DC.
- Liu, J., M. Delshad, G.A. Pope, and K. Sepehrnoori. 1994. Application of higher-order flux-limited methods in compositional simulation. Transp. Porous Media 16:1–29.[CrossRef]
- Moridis, G.J., Q. Hu, Y.S. Wu, and G.S. Bodvarsson. 2003. Preliminary 3-D site-scale studies of radioactive colloid transport in the unsaturated zone at Yucca Mountain, Nevada. J. Contam. Hydrol. 60:251–286.[CrossRef][Medline]
- Narasimhan, T.N., and P.A. Witherspoon. 1976. An integrated finite difference method for analyzing fluid flow in porous media. Water Resour. Res. 12(1):57–64.[CrossRef]
- Oldenburg, C.M., and K. Pruess. 2000. Simulation of propagating fronts in geothermal reservoirs with the implicit Leonard total variation diminishing scheme. Geothermics 29:1–25.[CrossRef][Web of Science]
- Oldenburg, C.M., and K. Pruess. 1997. Higher-order differencing for geothermal reservoir simulation. In Proc. Workshop on Geothermal Reservoir Engineering, 22nd, Stanford, CA. 27–29 Jan. Stanford Univ., Stanford, CA.
- Pruess, K. 1983. GMINC: A mesh generator for flow simulations in fractured reservoirs. LBL-15227. Lawrence Berkeley National Laboratory, Berkeley, CA.
- Pruess, K. 1991. TOUGH2: A general-purpose numerical simulator for multiphase fluid and heat flow. LBL-29400. Lawrence Berkeley National Laboratory, Berkeley, CA.
- Robinson, B.A., C. Li, and C.K. Ho. 2003. Performance assessment model development and analysis of radionuclide transport in the unsaturated zone, Yucca Mountain, Nevada. J. Contam. Hydrol. 62–63:249–268.[CrossRef][Medline]
- Scheidegger, A.E. 1961. General theory of dispersion in porous media. J. Geophys. Res. 66:3273–3278.[CrossRef]
- Sweby, P.K. 1984. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM J. Num. Anal. 21:995–1011.[CrossRef]
- van Genuchten, M.Th., and W.J. Alves. 1982. Analytical solutions of the one-dimensional convective-dispersive solute transport equation. USDA Technical Bull. no. 1661. USDA, Washington, DC.
- Unger, A.J.A., P.A. Forsyth, and E.A. Sudicky. 1996. Variable spatial and temporal weighting schemes for use in multi-phase compositional problems. Adv. Water Resour. 19(1):1–27.[Medline]
- Viswanathan, H.S., B.J. Robinson, A.J. Valocchi, and I.R. Triay. 1998. A reactive transport model of neptunium migration from the potential repository at Yucca Mountain. J. Hydrol. 209:251–280.[CrossRef]
- Wu, Y.S., C.F. Ahlers, P. Fraser, A. Simmons, and K. Pruess. 1996. Software qualification of selected TOUGH2 modules. LBL-39490, UC-800. Lawrence Berkeley National Laboratory, Berkeley, CA.
- Wu, Y.S., L. Pan, W. Zhang, and G.S. Bodvarsson. 2002. Characterization of flow and transport processes within the unsaturated zone of Yucca Mountain, Nevada. J. Contam. Hydrol. 54:215–247.[CrossRef][Web of Science][Medline]
- Wu, Y.S., and K. Pruess. 2000. Numerical simulation of non-isothermal multiphase tracer transport in heterogeneous fractured porous media. Adv. Water Resour. 23:699–723.[CrossRef]
This article has been cited by other articles:
|
|
|
|
 
H.-H. Liu and T. H. Illangasekare
Preface: Recent Advances in Modeling Multiphase Flow and Transport with the TOUGH Family of Codes
Vadose Zone J.,
February 25, 2008;
7(1):
284 - 286.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
