Publications

Abstracts

• Jichun Li, T. Arbogast, and Yunqing Huang, Mixed methods using standard conforming finite elements, submitted.
  
  We investigate the mixed finite element method (MFEM) for solving a second order elliptic problem with a lowest order term, as might arise in the simulation of single phase flow in porous media. We find that traditional mixed finite element spaces are not necessary when a positive lowest order (i.e., reaction) term is present. Hence we propose to use standard conforming finite elements Q_k x (Q_k)^d on rectangles or P_k x (P_k)^d on simplices to solve for both the pressure and velocity field in d dimensions. The price we pay is that we have only sub-optimal order error estimates. With a delicate superconvergence analysis, we find some improvement for the simplest pair Q_k x (Q_k)^d with any k >= 1, or for P₁ x (P₁)^d, when the mesh is uniform and the solution has one extra order of regularity. We also prove similar results for both parabolic and second order hyperbolic problems. Numerical results using Q₁ x (Q₁)² and P₁ x (P₁)² are presented in support of our analysis. These observations allow us to simplify the implementation of the MFEM, especially for higher order approximations, as might arise in an hp-adaptive procedure.
  
  
• T. Arbogast and M. S. M. Gomez, A discretization and multigrid solver for a Darcy-Stokes system of three dimensional vuggy porous media, submitted.
  
  We develop a finite element discretization and multigrid solver for a Darcy-Stokes system of three dimensional vuggy porous media, i.e., porous media with cavities. The finite ele- ment method uses low order mixed finite elements in the Darcy and Stokes domains, and special transition elements near the Darcy-Stokes interface to allow for tangential discontinuities implied by the Beavers-Joseph boundary condition. We design a multigrid method to solve the resulting sad- dle point linear system. The intertwining of the Darcy and Stokes subdomains makes the resulting matrix highly ill-conditioned. The velocity field is very irregular, and its discontinuous tangential component at the Darcy-Stokes interface makes it difficult to define intergrid transfer operators. Our definition is based on mass conservation and the analysis of the orders of magnitude of the solution. The coarser grid equations are defined using the Galerkin method. A new smoother of Uzawa type is developed based on taking an optimal step in a good search direction. Our algorithm has a measured convergence factor independent of the size of the system, at least when there are no disconnected vugs. We study the macroscopic effective permeability of a vuggy medium, showing that the influ- ence of vug orientation, shape, and, most importantly, interconnectivity determine the macroscopic flow properties of the medium.
  
  
• R. Naimi-Tajdar, C. Han, K. Sepehrnoori, T. J. Arbogast, and M. A. Miller, A Fully Implicit, Compositional, Parallel Simulator for IOR Processes in Fractured Reservoirs, SPE Journal 12:3 (September 2007).
  
  Naturally fractured reservoirs contain a significant amount of the world oil reserves. A number of these reservoirs contain several billion barrels of oil. Accurate and efficient reservoir simulation of naturally fractured reservoirs is one of the most important, challenging, and computationally intensive problems in reservoir engineering. Parallel reservoir simulators developed for naturally fractured reservoirs can effectively address the computational problem.
  A new accurate parallel simulator for large-scale naturally fractured reservoirs, capable of modeling fluid flow in both rock matrix and fractures, has been developed. The simulator is a parallel, 3D, fully implicit, equation-of-state compositional model that solves very large, sparse linear systems arising from discretization of the governing partial differential equations. A generalized dual-porosity model, the multiple-interacting-continua (MINC), has been implemented in this simulator. The matrix blocks are discretized into subgrids in both horizontal and vertical directions to offer a more accurate transient flow description in matrix blocks. We believe this implementation has led to a unique and powerful reservoir simulator that can be used by small and large oil producers to help them in the design and prediction of complex gas and waterflooding processes on their desktops or a cluster of computers. Some features of this simulator, such as modeling both gas and water processes and the ability of 2D matrix subgridding to the best of our knowledge are not available in any commercial simulator. The code was developed on a cluster of processors, which has proven to be a very efficient and convenient resource for developing parallel programs.
  The results were successfully verified against analytical solutions and commercial simulators (ECLIPSE and GEM). Excellent results were achieved for a variety of reservoir case studies. Applications of this model for several IOR processes (including gas and water injection) are demonstrated. Results using the simulator on a cluster of processors are also presented. Excellent speedup ratios were obtained.
  
  
• Reprint: T. Arbogast, G. Pencheva, M. F. Wheeler, and I. Yotov, A multiscale mortar mixed finite element method, Multiscale Modeling and Simulation 6 (2007), pp. 319-346
  
  We develop multiscale mortar mixed finite element discretizations for second order elliptic equations. The continuity of flux is imposed via a mortar finite element space on a coarse grid scale, while the equations in the coarse elements (or subdomains) are discretized on a fine grid scale. The polynomial degree of the mortar and subdomain approximation spaces may differ; in fact, the mortar space achieves approximation comparable to the fine scale on its coarse grid by using higher order polynomials. Our formulation is related to, but more flexible than, existing multiscale finite element and variational multiscale methods. We derive a priori error estimates and show, with appropriate choice of the mortar space, optimal order convergence and some superconvergence on the fine scale for both the solution and its flux. We also derive efficient and reliable a posteriori error estimators, which are used in an adaptive mesh refinement algorithm to obtain appropriate subdomain and mortar grids. Numerical experiments are presented in confirmation of the theory.
  
  
• T. Arbogast and D. S. Brunson, A computatonal method for approximating a Darcy-Stokes system governing a vuggy porous medium, to appear in Computational Geosciences
  
  We develop and analyze a mixed finite element method for the solution of an elliptic system modeling a porous medium with large cavities, called vugs. It consists of a second order elliptic (i.e., Darcy) equation on part of the domain coupled to a Stokes equation on the rest of the domain, and a slip boundary condition (due to Beavers-Joseph-Saffman) on the interface between them. The tangential velocity is not continuous on the interface. We consider a 2-D vuggy porous medium with many small cavities throughout its extent, so the interface is not isolated. We use a certain conforming Stokes element on rectangles, slightly modified near the interface to account for the tangential discontinuity. This gives a mixed finite element method for the entire Darcy-Stokes system with a regular sparsity pattern that is easy to implement, independent of the vug geometry, as long as it aligns with the grid. We prove optimal global first order L² convergence of the velocity and pressure, as well as the velocity gradient in the Stokes domain. Numerical results verify these rates of convergence, and even suggest somewhat better convergence in certain situations. Finally, we present a lower dimensional space that uses Raviart-Thomas elements in the Darcy domain and uses our new modified elements near the interface in transition to the Stokes elements.
  
  
• Reprint: T. Arbogast, Ch.-S. Huang, and S.-M. Yang, Improved accuracy for alternating-direction methods for parabolic equations based on regular and mixed finite elements, Mathematical Models & Methods in Applied Sciences 17 (2007), pp. 1279-1305
  
  An efficient modification by Douglas and Kim of the usual alternating directions method reduces the splitting error from O(k²) to O(k³) in time step k. We prove convergence of this modified alternating directions procedure, for the usual non-mixed Galerkin finite element and finite difference cases, under the restriction that k/h² is sufficiently small, where h is the grid spacing. This improves the results of Douglas and Gunn, who require k/h⁴ to be sufficiently small, and Douglas and Kim, who require that the locally one-dimensional operators commute. We propose a similar and efficient modification of alternating directions for mixed finite element methods that reduces the splitting error to O(k³), and we prove convergence in the noncommuting case, provided that k/h² is sufficiently small. Numerical computations illustrating the mixed finite element results are also presented. They show that our proposed modification can lead to a significant reduction in the alternating direction splitting error.
  
  
• T. Arbogast and Chieh-Sen Huang, A fully mass and volume conserving implementation of a characteristic method for transport problems, SIAM J. Sci. Comput. 28 (2006), pp. 2001-2022
  
  The characteristics-mixed method considers the transport not of a single point or fluid particle, but rather the mass in an entire region of fluid. This mass is transported along the characteristic curves of the hyperbolic part of the transport equation, and the scheme thereby produces very little numerical dispersion, conserves mass locally, and can use long time steps. However, since the shape of a characteristic trace-back region must be approximated in numerical implementation, its volume may be incorrect, resulting in inaccurate concentration densities and, further, inaccurate reaction dynamics. We present a simple modification to the characteristics-mixed method that conserves both mass and volume of the transported fluid regions. Our algorithm also handles boundary conditions through a space-time change of variables in the trace-back routines, which allows the boundary to be treated as if it were interior to the domain. Nearly point sources, such as wells, present special difficulties, since characteristic trace-back curves converge in their vicinity. We also present techniques that allow one to conservatively implement wells. The techniques are illustrated in five numerical examples.
  
  
• T. Arbogast and H. L. Lehr, Homogenization of a Darcy-Stokes system modeling vuggy porous media, Computational Geosciences, 10 (2006), pp. 291-302.
  
  We derive a macroscopic model for single phase, incompressible, viscous fluid flow in a porous medium with small cavities called vugs. We model the vuggy medium on the microscopic scale using Stokes equations within the vugular inclusions, Darcy's law within the porous rock, and a Beavers-Joseph-Saffman boundary condition on the interface between the two regions. We assume periodicity of the medium, and obtain uniform energy estimates independent of the period. Through a two-scale homogenization limit as the period tends to zero, we obtain a macroscopic Darcy's law governing the medium on larger scales. We also develop some needed generalizations of the two-scale convergence theory needed for our bi-modal medium, including a two-scale convergence result on the Darcy-Stokes interface. The macroscopic Darcy permeability is computable from the solution of a cell problem. An analytic solution to this problem in a simple geometry suggests that: (1) flow along vug channels is primarily Poiseuille with a small perturbation related to the Beavers-Joseph slip, and (2) flow that alternates from vug to matrix behaves as if the vugs have infinite permeability.
  
  
• Reprint: T. Arbogast and K. J. Boyd, Subgrid Upscaling and Mixed Multiscale Finite Elements, SIAM J. Numer. Anal., 44 (2006), pp. 1150-1171.
  
  Second order elliptic problems in divergence form with a highly varying leading order coefficient on the scale e can be approximated on coarse meshes of spacing H > > e only if one uses special techniques. The mixed variational multiscale method, also called subgrid upscaling, can be used, and this method is extended to allow oversampling of the local subgrid problems. The method is shown to be equivalent to the multiscale finite element method when one uses the lowest order Raviart-Thomas spaces and provided that there are no fine scale components in the source function f. In the periodic setting, a multiscale error analysis based on homogenization theory of the more general subgrid upscaling method shows that the error is O(e + H^m + sqrt(e/H)), where m = 1. Moreover, m = 2 if one uses the second order Brezzi-Douglas-Marini or Brezzi-Douglas-Duran-Fortin spaces and no oversampling. The error bounding constant depends only on the H^m-1-norm of f and so is independent of small scales when m = 1. When oversampling is not used, a superconvergence result for the pressure approximation is shown.
  
  
• Reprint: T. Arbogast and M. F. Wheeler, A family of rectangular mixed elements with a continuous flux for second order elliptic problems, SIAM J. Numer. Anal., 42 (2005), pp. 1914-1931.
  
  We present a family of mixed finite element spaces for second order elliptic equations in two and three space dimensions. Our spaces approximate the vector flux by a continuous function. Our spaces generalize certain spaces used for approximation of Stokes problems. The finite element method incorporates projections of the Dirichlet data and certain low order terms. The method is locally conservative on the average. Suboptimal convergence is proven and demonstrated numerically. The key result is to construct a flux pi-projection operator that is bounded in the Sobolev space H¹, preserves a projection of the divergence, and approximates optimally. Moreover, the corresponding Raviart-Thomas flux preserving pi-projection operator is an L²-projection when restricted to this family of spaces.
  
  
• Reprint: T. Arbogast, Analysis of a two-scale, locally conservative subgrid upscaling for elliptic problems, SIAM J. Numer. Anal., 42 (2004), pp. 576-598.
  
  We present a two-scale theoretical framework for approximating the solution of a second order elliptic problem. The elliptic coefficient is assumed to vary on a scale that can be resolved on a fine numerical grid, but limits on computational power require that computations be performed on a coarse grid. We consider the elliptic problem in mixed variational form over W x V subset of L² x H(div). We base our scale expansion on local mass conservation over the coarse grid. It is used to define a direct sum decomposition of W x V into coarse and "subgrid" subspaces W_c x V_c and dW x dV such that (1) div V_c = W_c and div dV = dW, and (2) the space dV is locally supported over the coarse mesh. We then explicitly decompose the variational problem into coarse and subgrid scale problems. The subgrid problem gives a well-defined operator taking W_c x V_c to dW x dV, which is localized in space, and it is used to upscale, that is, to remove the subgrid from the coarse-scale problem. Using standard mixed finite element spaces, two-scale mixed spaces are defined. A mixed approximation is defined, which can be viewed as a type of variational multiscale method or a residual-free bubble technique. A numerical Green's function approach is used to make the approximation to the subgrid operator efficient to compute. A mixed method pi-operator is defined for the two-scale approximation spaces and used to show optimal order error estimates.
  
  
• T. Arbogast, An overview of subgrid upscaling for elliptic problems in mixed form, in Current Trends in Scientific Computing, Z. Chen, R. Glowinski, and K. Li, eds., Contemporary Mathematics, AMS, 2003, pp. 21-32.
  
  We present an overview of recent work dealing with upscaling second order elliptic problems in mixed form. We use a direct sum decomposition of the solution space into coarse and localized "subgrid" spaces. We use these to construct a two-scale variational form. A numerical Greens function approach allows for its efficient approximation. A three dimensional computational result representing flow in a porous medium illustrates the performance of the technique in approximating fine scales on coarse grids.
  
  
• T. Arbogast and S. L. Bryant, A Two-Scale Numerical Subgrid Technique for Waterflood Simulations, SPE J., Dec. 2002, pp. 446-457.
  
  We present a two-scale numerical subgrid technique for simulating waterflooding. Local subgrid computations are combined with a coarse grid computation to provide a fine scale resolution of the solution. We use on the fine scale porosity, relative and absolute permeabilities, the location of wells, and capillary pressure curves. No explicit macroscopic coefficients nor pseudo-functions result. The method is several times faster than solving the fine scale problem directly, generally more robust, and yet achieves good results as it requires no ad hoc assumptions at the coarse scale and retains all the physics of the original multiphase flow equations.
  
  
• T. Arbogast, Implementation of a Locally Conservative Numerical Subgrid Upscaling Scheme for Two-Phase Darcy Flow, Computational Geosciences, 6 (2002), pp. 453-48
  
  We present a locally mass conservative scheme for the approximation of two-phase flow in a porous medium that allows us to obtain detailed fine scale solutions on relatively coarse meshes. The permeability is assumed to be resolvable on a fine numerical grid, but limits on computational power require that computations be performed on a coarse grid. We define a two-scale mixed finite element space and resulting method, and describe in detail the solution algorithm. It involves a coarse scale operator coupled to a subgrid scale operator localized in space to each coarse grid element. An influence function (numerical Greens function) technique allows us to solve these subgrid scale problems independently of the coarse grid approximation. The coarse grid problem is modified to take into account the subgrid scale solution and solved as a large linear system of equations posed over a coarse grid. Finally, the coarse scale solution is corrected on the subgrid scale, providing a fine grid representation of the solution. Numerical examples are presented, which show that near-well behavior and even extremely heterogeneous permeability barriers and streaks are upscaled well by the technique.
  
  
• T. Arbogast and S.L. Bryant, Numerical subgrid upscaling for waterflood simulations, In Proceedings of the 16th SPE Symposium on Reservoir Simulation held in Houston, Texas, February 11-14, 2001, SPE 66375.
  
  We present a subgrid-scale numerical technique for upscaling waterflood simulations. We scale up the usual parameters porosity and relative and absolute permeabilities, and also the location of wells and capillary pressure curves. Some of these are critical nonlinear terms that need to be resolved on the fine scale, or serious errors will result. Upscaling is achieved by explicitly decomposing the differential system into a coarse-grid-scale operator coupled to a subgrid-scale operator. The subgrid-scale operator is approximated as an operator localized in space to a coarse-grid element. An influence function (numerical Greens function) technique allows us to solve these subgrid-scale problems independently of the coarse-grid approximation. The coarse-grid problem is modified to take into account the subgrid-scale solution and solved as a large linear system of equations. Finally, the coarse scale solution is corrected on the subgrid-scale, providing a fine-grid scale representation of the solution. In this approach, no explicit macroscopic coefficients nor pseudo-functions result. The method is easily seen to be optimally convergent in the case of a single linear parabolic equation. The method is several times faster than solving the fine-scale problem directly, generally more robust, and yet achieves good results as it requires no ad hoc assumptions at the coarse scale and retains all the physics of the original multiphase flow equations.
  
  
• T. Arbogast, Numerical subgrid upscaling of two-phase flow in porous media, In Z. Chen, R. E. Ewing, and Z.-C. Shi, editors, Numerical treatment of multiphase flows in porous media, volume 552 of Lecture Notes in Physics, pages 35-49. Springer, Berlin, 2000.
  
  We present an approach and numerical results for scaling up fine grid information to coarse scales in an approximation to a nonlinear parabolic system governing two-phase flow in porous media. The technique allows upscaling of the usual parameters porosity and relative and absolute permeabilities, and also the location of wells and capillary pressure. Some of these are critical nonlinear terms that need to be resolved on the fine scale, or serious errors will result. Upscaling is achieved by explicitly decomposing the differential system into a coarse-grid-scale operator coupled to a subgrid-scale operator, which we localize by imposing a closure assumption. We approximate the coarse-grid-scale operator with a mixed finite element method that has a second order accurate velocity coupled implicitly to the subgrid scale. The subgrid-scale operator is approximated locally by a first order accurate mixed method. A numerical Greens influence function technique allows us to solve these subgrid problems independently of the coarse-grid approximation. No explicit macroscopic coefficients nor pseudo-functions result. The method is easily seen to be optimally convergent in the case of a single linear parabolic equation.
  
  
• T. Arbogast, L. C. Cowsar, M. F. Wheeler, and I. Yotov, Mixed finite element methods on non-matching multiblock grids, SIAM J. Numer. Anal., 37:1295-1315, 2000.
  
  We consider mixed finite element methods for second order elliptic equations on non-matching multiblock grids. A mortar finite element space is introduced on the non-matching interfaces. We approximate in this mortar space the trace of the solution, and we impose weakly a continuity of flux condition. A standard mixed finite element method is used within the blocks. Optimal order convergence is shown for both the solution and its flux. Moreover, at certain discrete points, superconvergence is obtained for the solution, and also for the flux in special cases. Computational results using an efficient parallel domain decomposition algorithm are presented in confirmation of the theory.
  
  
• T. Arbogast and S. Bryant, Efficient forward modeling for DNAPL site evaluation and remediation, In Bentley et al., editors, Computational Methods in Water Resources XIII, pages 161-166, Rotterdam, 2000. Balkema.
  
  Although the general characteristics of DNAPL flow and transport in the subsurface are reasonably well understood, it is often difficult and expensive to pinpoint sources of DNAPL contamination. Inversion techniques to improve site characterization rely on a forward model of multiphase flow. Ideally the forward model would be very fast, so that many realizations can be carried out in order to quantify and reduce uncertainty, yet capable of handling large numbers of grid elements, so that more accurate (small scale) determinations of soil properties and DNAPL content can be made. To meet these conflicting requirements of speed and detail in the forward modeling of contamination events, we present a subgrid-scale numerical technique for upscaling multiphase flow. Upscaling is achieved by explicitly decomposing the differential system into a coarse-grid-scale operator coupled to a subgrid-scale operator. The subgrid-scale operator is approximated as an operator localized in space to a coarse-grid element. An influence function (numerical Greens function) technique allows us to solve these subgrid-scale problems independently of the coarse-grid approximation. The coarse-grid problem is modified to take into account the subgrid-scale solution and solved as a large linear system of equations. Finally, the coarse scale solution is corrected on the subgrid-scale, providing a fine-grid scale representation of the solution. In this approach, no explicit macroscopic coefficients nor pseudo-functions result. The method is easily seen to be optimally convergent in the case of a single linear parabolic equation. The method is fast, robust, and achieves good results.
  
  
• M. Wheeler, T. Arbogast, S. Bryant, J. Eaton, Qin Lu, M. Peszynska, and I. Yotov, A parallel multiblock/multidomain approach for reservoir simulation, In Proceedings of the 15th SPE Symposium on Reservoir Simulation held in Houston, Texas, February 14-17, 1999, SPE 51884.
  
  Our approach for parallel multiphysics and multiscale simulation uses two levels of domain decomposition: physical and computational. First, the physical domain is decomposed into subdomains or blocks according to the geometry, geology, and physics/chemistry/biology. These subdomains represent a single physical system, on a reasonable range of scales, such as a black oil region, a compositional region, a region to one side of a fault, or a near-wellbore region. Second, the computations are further decomposed on a parallel machine for efficiency. That is, we use a multiblock or macro-hybrid approach, in which we describe a domain as a union of regions or blocks, and employ an appropriate hierarchical model on each block. This approach allows one to define grids and computations independently on each block. This local grid structure has many advantages. It allows the most efficient and accurate discretization techniques to be employed in each block. The multiblock structure of the algebraic systems allows for the design and use of efficient domain decomposition solvers and preconditioners. Decomposition into independent blocks offers great flexibility in accommodating the shape of the external boundary, the presence of internal features such as faults and wells, the need to refine a region of the domain in space or time (by treating it as a distinct block); interfacing structured and unstructured grids; and accommodating various models of multiscale and multiphysical phenomena. The resulting grid is not suited to direct application of discretization methods. We use mortar space techniques to impose physically meaningful, mass conservative, flux-matching conditions on the interfaces between blocks. We present numerical simulations to illustrate several of these decomposition strategies, including the coupling of IMPES and fully implicit models and upscaling by varying the number of degrees of freedom on the block interfaces.
  
  
• T. Arbogast, C. N. Dawson, P. T. Keenan, M. F. Wheeler, and I. Yotov, Enhanced cell-centered finite differences for elliptic equations on general geometry, SIAM J. Sci. Comput., 19:404-425, 1998.
  
  We consider the efficient implementation of mixed finite elements for solving second order elliptic partial differential equations on geometrically general domains, concentrating on the lowest-order Raviart-Thomas approximating spaces. We consider the standard mixed method and its hybrid form, and the recently introduced expanded mixed method. The standard method yields a saddle-point linear system, and while the hybrid method yields a positive definite linear system, it has many more unknowns, one per element edge or face. The expanded mixed method is similar in its structure; however, we give a generalization of the method combined with a global mapping technique that makes it suitable for general meshes. Moreover, two quadrature rules are given which reduce the method to a cell-centered finite difference method on meshes of quadrilaterals or triangles in 2 dimensions and hexahedra or tetrahedra in 3 dimensions. This approach substantially reduces the complexity of the mixed finite element matrix, leaving a symmetric, positive definite system for only as many unknowns as elements. On smooth meshes that are either logically rectangular or triangular with six triangles per internal vertex, this finite difference method is as accurate as the standard mixed method; on non-smooth meshes it can lose accuracy. An enhancement of the method is defined that combines numerical quadrature with Lagrange multiplier pressures on certain element edges or faces. The enhanced method regains the accuracy of the solution on non-smooth meshes, with little additional cost if the mesh geometry is piece-wise smooth, as in hierarchical meshes. Theoretical error estimates and numerical examples are given comparing the accuracy and efficiency of the methods.
  
  
• M. F. Wheeler, T. Arbogast, S. Bryant, and J. Eaton, Efficient parallel computation of spatially heterogeneous geochemical reactive transport, In V. N. Burganos et al., editors, Computational Methods in Water Resources XII, Vol. 1: Computational Methods in Contamination and Remediation of Water Resources, pages 453-460, Southampton, U.K., 1998. Computational Mechanics Publications.
  
  In flow and non-reactive transport problems, an efficient division of labor for a parallel computer is given by a decomposition of the domain that assigns roughly an equal number of grid cells to each processor. Typical geochemically reactive transport problems exhibit reaction fronts or zones that travel through the domain and may be concentrated near heterogeneities in the domain, such as wells or sources of contaminants. Within these fronts, kinetic reactions may be governed by relatively stiff differential equations, and equilibrium reactions may be near poorly conditioned minima in the Gibbs free energy. Thus the geochemical computational workload can vary greatly with position and time. A domain decomposition tailored for flow and transport problems will result in computationally intensive reaction zones that fall heavily on certain parallel processors but lightly on the rest. We present an algorithm to balance the geochemical load while preserving the efficiency of the flow and transport calculations. Significant performance improvements are observed on a variety of simple test problems.
  
  
• T. Arbogast, S. E. Minkoff, and P. T. Keenan, An operator-based approach to upscaling the pressure equation, In V. N. Burganos et al., editors, Computational Methods in Water Resources XII, Vol. 1: Computational Methods in Contamination and Remediation of Water Resources, pages 405-412, Southampton, U.K., 1998. Computational Mechanics Publications.
  
  Permeability and porosity parameters of a porous medium are known only in a statistical sense. For risk assessment, one must perform multiple flow simulations of a single site, varying these input parameters. Because multiple simulations of large sites are computationally prohibitive, upscaling from fine to coarse scales is necessary. Traditional upscaling techniques determine a new effective or upscaled permeability field defined on a coarser scale, which is then used in a standard coarse grid discretization operator. We develop here a method of determining a new coarse grid discretization operator that provides an upscaled solution but bypasses the determination of effective permeability and porosity fields. The method has two steps. We first solve for fine scale flow information internal to each coarse grid cell. Because the problems are small, this step is relatively fast. Then we determine a modified coarse grid operator for solving the upscaled problem that includes the fine scale flow information from the first step. The method is developed for single-phase flow in the context of the mixed finite element method; therefore, the method is locally mass conservative. Unlike traditional upscaling methods (such as homogenization) we do not impose arbitrary boundary conditions on the coarse grid. We present comparisons of our method with the harmonic average permeability upscaling technique.
  
  
• T. Arbogast and I. Yotov, A non-mortar mixed finite element method for elliptic problems on non-matching multiblock grids, Comp. Meth. in Appl. Mech. and Engng., 149:225-265, 1997.
  
  We consider the approximation of second order elliptic equations on domains that can be described as a union of sub-domains or blocks. We assume that a grid is defined on each block independently, so that the resulting grid over the entire domain need not be conforming (i.e., match) across the block boundaries. Several techniques have been developed to approximate elliptic equations on multiblock grids that utilize a mortar finite element space defined on the block boundary interface itself. We define a mixed finite element method that does not use such a mortar space. The method has an advantage in the case where adaptive local refinement techniques will be used, in that there is no mortar grid to refine. As is typical of mixed methods, our method is locally conservative element-by-element; it is also globally conservative across the block boundaries. Theoretical results show that the approximate solution converges at the optimal rate to the true solution. We present computational results to illustrate and confirm the theory.
  
  
• T. Arbogast, M. F. Wheeler, and I. Yotov, Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences, SIAM J. Numer. Anal., 34:828-852, 1997.
  
  We present an expanded mixed finite element approximation of second order elliptic problems containing a tensor coefficient. The mixed method is expanded in the sense that three variables are explicitly approximated, namely, the scalar unknown, the negative of its gradient, and its flux (the tensor coefficient times the negative gradient). The resulting linear system is a saddle point problem. In the case of the lowest order Raviart-Thomas elements on rectangular parallelepipeds, we approximate this expanded mixed method by incorporating certain quadrature rules. This enables us to write the system as a simple, cell-centered finite difference method, requiring the solution of a sparse, positive semidefinite linear system for the scalar unknown. For a general tensor coefficient, the sparsity pattern for the scalar unknown is a nine point stencil in two dimensions, and 19 points in three dimensions. Existing theory shows that the expanded mixed method gives optimal order approximations in the L² and H^-s-norms (and superconvergence is obtained between the L²-projection of the scalar variable and its approximation). We show that these rates of convergence are retained for the finite difference method. If h denotes the maximal mesh spacing, then the optimal rate is O(h). The superconvergence rate O(h²) is obtained for the scalar unknown and rate O(h^3/2) for its gradient and flux in certain discrete norms; moreover, the full O(h²) is obtained in the strict interior of the domain. Computational results illustrate these theoretical results.
  
  
• T. Arbogast, Computational aspects of dual-porosity models, In U. Hornung, editor, Homogenization and Porous Media, Interdisciplinary Applied Math. Series, pages 203-223. Springer, New York, 1997.
  
  Macroscopic microstructure models of dual-porosity type model, among other things, the flow of fluids in highly fractured porous media; that is, media comprised of porous matrix rock divided into relatively small blocks by thin fractures. Mathematically, a dual-porosity model is derived by homogenization of the mesoscopic equations (i.e., the Darcy-scale description that explicitly models flow within the fractures and matrix). This results in a relatively complex system of partial differential equations in seven variables, (t,x,y). At first glance, it may not be apparent that there is any advantage to the macroscopic description versus the mesoscopic. However, the dual-porosity model explicitly captures the length scales of the physical problem, and is thus much easier to approximate computationally. The homogenized, dual-porosity model requires far less computational effort to solve numerically, and it approximates well the mesoscopic model.
  
  
• Peng Wang, I. Yotov, M. Wheeler, T. Arbogast, C. Dawson, M. Parashar, and K. Sepehrnoori, A new generation EOS compositional reservoir simulator: Part I-Formulation and discretization, In Proceedings of the 14th SPE Symposium on Reservoir Simulation held in Dallas, Texas, June 8-11, 1997, SPE 37979.
  
  A fully implicit equation-of-state (EOS) computational simulator for large scale reservoir simulation is presented. The simulator uses a multiblock, domain decomposition approach; that is, the reservoir is divided into non-overlapping subdomains that are solved locally in parallel (inner iteration). The subdomain grids are defined independently of each other and their connections are attained through a global interface problem (outer iteration) formulated in terms of appropriate equations that guarantee continuity of total component fluxes. Parallel, iterative techniques are employed to solve the resulting nonlinear equations. The model formulation has been successfully tested for a dry gas cycling process on a single fault block. The numerical results show that the simulator and fluid-related calculations can be conducted efficiently and robustly. Promising results have been obtained using the proposed multiblock approach for nonmatching grids between fault blocks for two-phase flow problems.
  
  
• T. Arbogast, C. N. Dawson, P. T. Keenan, M. F. Wheeler, and I. Yotov, The application of mixed methods to subsurface simulation, In R. Helmig et al., editors, Modeling and Computation in Environmental Sciences, volume 59 of Notes on Numerical Fluid Mechanics, pages 1-13, Braunschweig, 1997. Vieweg Publ.
  
  We consider the application of mixed finite element and finite difference methods to groundwater flow and transport problems. We are concerned with accurate approximation and efficient implementation, especially when the porous medium may have geometric irregularities, heterogeneities, and either a tensor hydraulic conductivity or a tensor dispersion. For single-phase flow, we develop an expanded mixed finite element method defined on a logically rectangular, curvilinear grid. Special quadrature rules are introduced to transform the method into a simple cell-centered finite difference method. The approximation is locally conservative and highly accurate. We also show that the highly nonlinear two-phase flow problem is well approximated by mixed methods. The main difficulty is that the true solution is typically lacking in regularity.
  
  
• M. F. Wheeler, T. Arbogast, S. Bryant, C. N. Dawson, F. Saaf, and Chong Wang, New computational approaches for chemically reactive transport in porous media, In G. Delic and M.F. Wheeler, editors, Next Generation Environmental Models and Computational Methods (NGEMCOM), pages 217-226, Philadelphia, 1997. Proceedings of the U.S. Environmental Protection Agency Workshop (NGEMCOM), SIAM.
  
  Flow accompanied by chemical and nuclear reaction and mass transfer plays a central role in many environmental applications. The complexity of coupled flow, transport, and reaction phenomena imposes stringent performance criteria on numerical simulators. This article describes some recent developments that address these criteria. Particular emphasis is placed on robustness and parallel computation. Two practical examples illustrate briefly the role of numerical simulation in subsurface contaminant remediation and assessment.
  
  
• T. Arbogast, S. Bryant, C. Dawson, F. Saaf, Chong Wang, and M. Wheeler, Computational methods for multiphase flow and reactive transport problems arising in subsurface contaminant remediation, J. Computational Appl. Math., 74:19-32, 1996.
  
  A mathematical formulation and some numerical approximation techniques are described for a system of coupled partial differential and algebraic equations describing multiphase flow, transport and interactions of chemical species in the subsurface. A parallel simulator PARSSIM has been developed based on these approximation techniques and is being used to study contaminant remediation strategies. Numerical results for a highly complex geochemistry problem involving strontium disposal in a pit at Oak Ridge National Laboratory are presented.
  
  
• T. Arbogast, M. F. Wheeler, and Nai-Ying Zhang, A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media, SIAM J. Numer. Anal., 33:1669-1687, 1996.
  
  We study a model nonlinear, degenerate, advection-diffusion equation having application in petroleum reservoir and groundwater aquifer simulation. The main difficulty is that the true solution is typically lacking in regularity; therefore, we consider the problem from the point of view of optimal approximation. Through time integration, we develop a mixed variational form that respects the known minimal regularity, and then we develop and analyze two versions of a mixed finite element approximation, a simpler semidiscrete (time continuous) version and a fully discrete version. Our error bounds are optimal in the sense that all but one of the bounding terms reduce to standard approximation error. The exceptional term is a nonstandard approximation error term. We also consider our new formulation for the nondegenerate problem, showing the usual optimal L₂-error bounds; moreover, superconvergence is obtained under special circumstances.
  
  
• T. Arbogast, M. F. Wheeler, and I. Yotov, Logically rectangular mixed methods for flow in irregular, heterogeneous domains, In A. A. Aldama et al., editors, Computational Methods in Water Resources XI, volume 1, pages 621-628, Southampton, 1996. Computational Mechanics Publications.
  
  We develop and analyze a mixed finite element and a cell-centered finite difference method for groundwater flow in an irregular, heterogeneous, multi-block aquifer domain. The methods are designed to handle full tensor hydraulic conductivity with possible discontinuities. The domain can be divided into a series of smaller, non-overlapping sub-domain blocks of irregular geometry. Each is covered by a logically rectangular grid; these are not required to match on the interface so that we can model faults, local refinements, and other internal boundaries, such as interfaces where the conductivity is discontinuous. After continuously mapping each sub-domain to a rectangular reference sub-domain, all computations can be performed in a simple rectangular context. Standard mixed finite element spaces are used for the local sub-domain discretization. A ``mortar'' finite element space is introduced to accurately approximate the pressure along the sub-domain interfaces. Quadrature rules are employed to transform the mixed finite element method into cell-centered finite differences for the pressures. Theoretical and computational results show that the scheme is highly accurate. Super-convergence for both pressure and velocity is obtained at certain discrete points. Three dimensional numerical examples using an efficient parallel domain decomposition solver are also presented.
  
  
• T. Arbogast, C. N. Dawson, and M. F. Wheeler, A parallel algorithm for two phase multicomponent contaminant transport, Applications of Math., 40:163-174, 1995.
  
  We discuss the formulation of a simulator in three spatial dimensions for a multicomponent, two phase (air, water) system of groundwater flow and transport with biodegradation kinetics and wells with multiple screens. The simulator has been developed for parallel, distributed memory, message passing machines. The numerical procedures employed are a fully implicit expanded mixed finite element method for flow and either a characteristics-mixed method or a Godunov method for transport and reactions of dissolved chemical species in groundwater. Domain decomposition, symmetric and nonsymmetric solvers have been developed for solving the systems of equations resulting from the discretization of the model. Results from applying this simulator to a bioremediation field problem with several injection and production wells each having multiple screens are presented.
  
  
• T. Arbogast and Zhangxin Chen, On the implementation of mixed methods as nonconforming methods for second order elliptic problems, Math. Comp., 64:943-972, 1995.
  
  In this paper we show that mixed finite element methods for a fairly general second order elliptic problem with variable coefficients can be given a nonmixed formulation. (Lower order terms are treated, so our results apply also to parabolic equations.) We define an approximation method by incorporating some projection operators within a standard Galerkin method, which we call a projection finite element method. It is shown that for a given mixed method, if the projection method's finite element space M_h satisfies three conditions, then the two approximation methods are equivalent. These three conditions can be simplified for a single element in the case of mixed spaces possessing the usual vector projection operator. We then construct appropriate nonconforming spaces M_h for the known triangular and rectangular elements. The lowest-order Raviart-Thomas mixed solution on rectangular finite elements in R² and R³, on simplices, or on prisms, is then implemented as a nonconforming method modified in a simple and computationally trivial manner. This new nonconforming solution is actually equivalent to a postprocessed version of the mixed solution. A rearrangement of the computation of the mixed method solution through this equivalence allows us to design simple and optimal order multigrid methods for the solution of the linear system.
  
  
• T. Arbogast and M. F. Wheeler, A characteristics-mixed finite element method for advection dominated transport problems, SIAM J. Numer. Anal., 32:404-424, 1995.
  
  We define a new finite element method, called the characteristics-mixed method, for approximating the solution to an advection dominated transport problem. The method is based on a space-time variational form of the advection-diffusion equation. Our test functions are piecewise constant in space, and in time they approximately follow the characteristics of the advective (i.e., hyperbolic) part of the equation. Thus the scheme uses a characteristic approximation to handle advection in time. This is combined with a low order mixed finite element spatial approximation of the equation. Boundary conditions are incorporated in a natural and mass conservative fashion. The scheme is completely locally conservative; in fact, on the discrete level, fluid is transported along the approximate characteristics. A post-processing step is included in the scheme in which the approximation to the scalar unknown is improved by utilizing the approximate vector flux. This has the effect of improving the rate of convergence of the method. We show that it is optimally convergent to order one in time and at least suboptimally convergent to order 3/2 in space.
  
  
• T. Arbogast, Mixed methods for flow and transport problems on general geometry, In G. F. Carey, editor, Finite Element Modeling of Environmental Problems, pages 275-286, Cichester, England, 1995. Wiley.
  
  A groundwater flow or transport problem requires the solution of a second order elliptic equation with a tensor hydraulic conductivity or dispersion. These problems are posed over an aquifer domain with varying topography due to the geological layering, surface features, and the like. We present an expanded mixed finite element method that can efficiently handle these difficulties. The approximating spaces are defined on a smooth curvilinear grid, obtained by a global mapping of a simple, computational grid to the aquifer domain. Quadrature rules are introduced to transform the mixed method into a cell-centered finite difference method for the pressure or concentration. If rectangular computational elements are used, the stencil is 9 points in 2 dimensions and 19 points in 3 dimensions. Triangular computational elements give a 10 point stencil. The resulting scheme is locally mass conservative. For flow, the linear Galerkin finite element method gives first order accurate velocities, while the rectangular mixed method is second order accurate in the interior of the domain.
  
  
• T. Arbogast, P. T. Keenan, M. F. Wheeler, and I. Yotov, Logically rectangular mixed methods for darcy flow on general geometry, In Proceedings of the 13th SPE Symposium on Reservoir Simulation held in San Antonio, Texas, pages 51-59, February 12-15, 1995, SPE 29099.
  
  We consider an expanded mixed finite element formulation (cell centered finite differences) for Darcy flow with a tensor absolute permeability. The reservoir can be geometrically general with internal features, but the computational domain is rectangular. The method is defined on a curvilinear grid that need not be orthogonal, obtained by mapping the rectangular, computational grid. The original flow problem becomes a similar problem with a modified permeability on the computational grid. Quadrature rules turn the mixed method into a cell-centered finite difference method with a 9 point stencil in 2-D and 19 in 3-D. As shown by theory and experiment, if the modified permeability on the computational domain is smooth, then the convergence rate is optimal and both pressure and velocity are superconvergent at certain points. If not, Lagrange multiplier pressures can be introduced on boundaries of elements so that optimal convergence is retained. This modification presents only small changes in the solution process; in fact, the same parallel domain decomposition algorithms can be applied with little or no change to the code if the modified permeability is smooth over the subdomains. This Lagrange multiplier procedure can be used to extend the difference scheme to multi-block domains, and to give a coupling with unstructured grids. In all cases, the mixed formulation is locally conservative. Computational results illustrate the advantage and convergence of this method.
  
  
• T. Arbogast, M. F. Wheeler, and I. Yotov, Logically rectangular mixed methods for groundwater flow and transport on general geometry, In A. Peters et al., editors, Computational Methods in Water Resources X, Vol. 1, pages 149-156, Dordrecht, The Netherlands, 1994. Kluwer Academic Publishers.
  
  We consider an extended mixed finite element formulation for groundwater flow and transport problems with either a tensor hydraulic conductivity or a tensor dispersion. While the aquifer domain can be geometrically general, in our formulation the computational domain is rectangular. The approximating spaces for the mixed method are defined on a smooth curvilinear grid, obtained by a global mapping of the rectangular, computational grid. The original problem is mapped to the computational domain, giving a similar problem with a modified tensor coefficient. Special quadrature rules are introduced to transform the mixed method into a simple cell-centered finite difference method with a 9 point stencil in 2-D and 19 point stencil in 3-D. The resulting scheme is locally mass conservative. In the case of flow, linear Galerkin procedures give first order accurate velocities, while our method is second order accurate. Both computational and theoretical results are given.
  
  
• T. Arbogast, C. N. Dawson, and M. F. Wheeler, A parallel multiphase numerical model for subsurface contaminant transport with biodegradation kinetics, In A. Peters et al., editors, Computational Methods in Water Resources X, Vol. 2, pages 1499-1506, Dordrecht, The Netherlands, 1994. Kluwer Academic Publishers.
  
  We discuss the formulation of a simulator in three spatial dimensions for two phase groundwater flow and transport with biodegradation kinetics that has been developed at Rice University for massively parallel, distributed memory, message passing machines. The numerical procedures employed are a fully implicit mixed finite element method for flow and a characteristics-mixed method for transport and reactions of dissolved chemical species in groundwater. Domain decomposition solvers have been employed for solving the systems of equations resulting from the discretization of the model. Results from applying this simulator to a bioremediation field problem using a recirculation well in an air-water system are discussed.
  
  
• T. Arbogast, C. N. Dawson, and P. T. Keenan, Efficient mixed methods for groundwater flow on triangular or tetrahedral meshes, In A. Peters et al., editors, Computational Methods in Water Resources X, Vol. 1, pages 3-10, Dordrecht, The Netherlands, 1994. Kluwer Academic Publishers.
  
  Simulating flow in porous media requires the solution of elliptic or parabolic partial differential equations. When the computational domain is irregularly shaped, applying finite element methods with triangular elements offers great flexibility. The mixed finite element method has proven useful for solving flow equations. The difficulty with mixed methods in general is in solving the linear algebraic systems that arise. On rectangular elements, the mixed method with lowest-order approximating spaces can be reduced to a simple finite difference method in the primary variable, thus reducing the linear system to a sparse, symmetric, positive (semi-)definite matrix, for which many solution techniques are known. This type of reduction is not straightforward for triangular elements. In this paper we outline new mixed type methods which generalize finite differences in a manner suitable for use with triangular elements. Numerical examples illustrate the accuracy and efficiency of these new methods.
  
  
• T. Arbogast and M. F. Wheeler, A parallel numerical model for subsurface contaminant transport with biodegradation kinetics, In J. R. Whiteman, editor, The Mathematics of Finite Elements and Applications VIII (MAFELAP 1993), pages 199-213, New York, 1994. Wiley.
  
  In this paper we discuss the formulation of a simulator for groundwater flow and transport with biodegradation kinetics that has been developed at Rice University for massively parallel, distributed memory, message passing machines. The numerical procedures employed are a mixed finite element method for flow and the characteristics-mixed method for transport. Kinetics are treated by time splitting. The linear solvers are based on domain decomposition. Application of this procedure to a bioremediation problem as well as numerical experiments on the INTEL i860 and INTEL Delta are discussed. Results indicate that the procedure is theoretically mass conservative over each grid cell and is approximately so in implementation. Preliminary tests indicate that the procedure is robust and applicable to realistic groundwater problems. Moreover, the numerical model scales almost linearly with the number of processors even for fairly coarse grids.
  
  
• T. Arbogast, Gravitational forces in dual-porosity systems. I. Model derivation by homogenization, Transport in Porous Media, 13:179-203, 1993.
  
  We consider the problem of modeling flow through naturally fractured porous media. In this type of media, various physical phenomena occur on disparate length scales, so it is difficult to properly average their effects. In particular, gravitational forces pose special problems. In this paper we develop a general understanding of how to incorporate gravitational forces into the dual-porosity concept. We accomplish this through the mathematical technique of formal two-scale homogenization. This technique enables us to average the single-porosity, Darcy equations that govern the flow on the finest (fracture thickness) scale. The resulting homogenized equations are of dual-porosity type. We consider three flow situations, the flow of a single component in a single phase, the flow of two fluid components in two completely immiscible phases, and the completely miscible flow of two components.
  
  
• T. Arbogast, Gravitational forces in dual-porosity systems. II. Computational validation of the homogenized model, Transport in Porous Media, 13:205-220, 1993.
  
  Three models are considered for single component, single phase flow in naturally fractured porous media. The microscopic model holds on the Darcy scale, and it is considered to govern the system. The macroscopic, dual-porosity model was derived in Part I of this work from the microscopic model by two-scale mathematical homogenization. In this paper, we show that the dual-porosity model predicts well the behavior of the microscopic model by comparing their computed solutions in certain reasonable test cases. Homogenization gives a complex formula for a key parameter in the dual-porosity model; herein a simple approximation to this formula is presented. The third model considered is a single-porosity model with averaged parameters. It is shown that this type of model cannot predict the behavior of the microscopic flow.
  
  
• T. Arbogast, M. Obeyesekere, and M. F. Wheeler, Numerical methods for the simulation of flow in root-soil systems, SIAM J. Numer. Anal., 30:1677-1702, 1993.
  
  We consider the numerical properties of approximation schemes for a model that simulates water transport in root-soil systems. The model is derived in detail. It is based on a previously proposed model which we reformulate completely in terms of the water potential. The system of equations consists of a parabolic partial differential equation which contains a nonlinear capacity term coupled to two linear ordinary differential equations. A closed form solution is obtained for one of the latter equations. Finite element and finite difference schemes are defined to approximate the solution of the coupled system. Some new techniques which have wide applicability for analyzing the nonlinear capacity term are used, and optimal order error estimates are derived. A postprocessed water mass flux computation is also presented and shown to be superconvergent to the true flux. Computational results which verify the theoretical convergence rates are given.
  
  
• J. Douglas, Jr., T. Arbogast, P. J. Paes Leme, J. L. Hensley, and N. P. Nunes, Immiscible displacement in vertically fractured reservoirs, Transport in Porous Media, 12:73-106, 1993.
  
  A dual-porosity model is defined for saturated, two-phase, compressible, immiscible flow in a vertically fractured reservoir or aquifer. This model allows detailed simulation of the matrix-fracture interaction as well as the matrix flow itself. This is accomplished by directly coupling the matrix and fracture systems along the vertical faces of the matrix blocks, incorporating gravitational effects directly, and simulating flow inside the block. Thus fluid segregation due to gravitational effects and heterogeneities can be simulated. We show that our model can be derived via homogenization techniques. The model (in incompressible form for simplicity of exposition) is then approximated by a computationally efficient finite difference scheme. Calculations are presented to show the convergence of the scheme and to indicate the behavior of the model as a function of several parameters.
  
  
• T. Arbogast, The existence of weak solutions to single-porosity and simple dual-porosity models of two-phase incompressible flow, Journal of Nonlinear Analysis: Theory, Methods, and Applications, 19:1009-1031, 1992.
  
  It is shown that there exists a weak solution to a degenerate and singular elliptic-parabolic partial integro-differential system of equations. These equations model two-phase incompressible flow of immiscible fluids in either an ordinary porous medium or in a naturally fractured porous medium. The full model is of dual-porosity type, though the single porosity case is covered by setting the matrix-to-fracture flow terms to zero. This matrix-to-fracture flow is modeled simply in terms of fracture quantities; that is, no distinct matrix equations arise. The equations are considered in a global pressure formulation that is justified by appealing to a physical relation between the degeneracy of the wetting fluid's mobility and the singularity of the capillary pressure function. In this formulation, the elliptic and parabolic parts of the system are separated; hence, it is natural to consider various boundary conditions, including mixed nonhomogeneous, saturation dependent ones of the first three types. A weak solution is obtained as a limit of solutions to discrete time problems. The proof makes no use of the corresponding regularized system. The hypotheses required for various earlier results on single-porosity systems are weakened so that only physically relevant assumptions are made. In particular, the results cover the cases of a singular capillary pressure function, a pure Neumann boundary condition, and an arbitrary initial condition.
  
  
• T. Arbogast, A simplified dual-porosity model for two-phase flow, In T. F. Russell et al., editors, Computational Methods in Water Resources IX, Vol. 2: Mathematical Modeling in Water Resources, pages 419-426, Southampton, U.K., 1992. Computational Mechanics Publications.
  
  A model for two-phase, incompressible, immiscible fluid flow in a highly fractured porous medium is derived as a simplification of a much more detailed dual-porosity model. This simplified model has a nonlinear matrix-fracture interaction, and it is more general than similar existing ``transfer function'' models. It is computationally less complex than the detailed model, and simulation results are presented which assess any loss in accuracy. It is shown that the new model approximates capillary effects quite well, and better than similar existing models.
  
  
• T. Arbogast, M. Obeyesekere, and M. F. Wheeler, Simulation of flow in root-soil systems, In T. F. Russell et al., editors, Computational Methods in Water Resources IX, Vol. 2: Mathematical Modeling in Water Resources, pages 195-202, Southampton, U.K., 1992. Computational Mechanics Publications.
  
  The goal of this work is to develop an understanding of the properties of root-soil systems, which can be modified by using genetic engineering techniques, in order to improve the performance of plants when water availability is limited. We develop a mathematical model of a root-soil system for plants with taproots. We approximate it with a finite difference algorithm which converges at the optimal rate, including a post-processed water flux approximation. The results of some numerical simulations are presented to show the effects of changes to the root-soil system.
  
  
• T. Arbogast, A. Chilakapati, and M. F. Wheeler, A characteristic-mixed method for contaminant transport and miscible displacement, In T. F. Russell et al., editors, Computational Methods in Water Resources IX, Vol. 1: Numerical Methods in Water Resources, pages 77-84, Southampton, U.K., 1992. Computational Mechanics Publications.
  
  Recently, Arbogast and Wheeler formulated and analyzed a modified method of characteristics-mixed method for approximating solutions to convection-diffusion equations. This scheme is theoretically mass conservative over each grid cell; it is approximately so in implementations. We consider application of this procedure to contaminant transport and to miscible displacement with unfavorable mobility ratio. Results in one, two, and three space dimensions are discussed.
  
  
• J. Douglas, Jr., J. L. Hensley, and T. Arbogast, A dual-porosity model for waterflooding in naturally fractured reservoirs, Comp. Meth. in Appl. Mech. and Engng., 87:157-174, 1991.
  
  In this paper, we are concerned with two-phase, immiscible, incompressible flow in naturally fractured reservoirs. We use a dual-porosity model derived in previous work. The model is presented on a double covering of the reservoir, with one cover, which we call Omega, containing the fracture flow. We attach a matrix block to each point in Omega, with these blocks being topologically disconnected, so that flow between a matrix block and its representative point in the fracture cover is permitted (and carefully modeled), while flow between individual blocks is not. An implicit yet effecient and parallelizable finite difference discretization scheme is defined. The fracture equations retain their nonlinearities, while the matrix equations are linearized. We present some simulations of waterflooding to illustrate various features of the model.
  
  
• T. Arbogast, Gravitational forces in dual-porosity models of single phase flow, In Proceedings, Thirteenth IMACS World Congress on Computation and Applied Mathematics, pages 607-608, Dublin, Ireland, July 22-26, 1991. Trinity College.
  
  A dual porosity model is derived by the formal theory of homogenization. The model properly incorporates gravity in that it respects the equilibrium states of the medium.
  
  
• T. Arbogast, M. Obeyesekere, and M. F. Wheeler, Convergence analysis for simulating flow in root-soil systems, In J. R. Whiteman, editor, The Mathematics of Finite Elements and Applications VII (MAFELAP 1990), pages 361-383, London, 1991. Academic Press.
  
  In this paper we consider the numerical properties of approximation schemes for a model that simulates water transport in root-soil systems. The system of equations consists of a nonlinear parabolic partial differential equation coupled to a linear ordinary differential equation, and they model the potential and intake of water in a root-soil system, respectively. Optimal order error estimates are derived for finite element and finite difference approximations to the solution of the coupled equations. The finite difference scheme is shown to be equivalent to a finite element scheme combined with certain quadrature rules; consequently, all of the error estimates are derived from finite element and quadrature analysis. Computational results which verify the theoretical convergence rates are given for the finite difference scheme.
  
  
• T. Arbogast, J. Douglas, Jr., and U. Hornung, Modeling of naturally fractured reservoirs by formal homogenization techniques, In R. Dautray, editor, Frontiers in Pure and Applied Mathematics, pages 1-19. Elsevier, Amsterdam, 1991.
  
  Unique forms of the double porosity model are derived for various types of flow in naturally fractured reservoirs. A single component in a single phase and two component miscible and immiscible flows are treated. These models are derived by homogenizing the appropriate equations describing flow in a highly discontinuous single porosity reservoir. The mathematical theory of homogenization is used only in its formal sense.
  
  
• J. Douglas, Jr. and T. Arbogast, Dual-porosity models for flow in naturally fractured reservoirs, In J. H. Cushman, editor, Dynamics of Fluids in Hierarchical Porous Media, pages 177-221. Academic Press, London, 1990.
  
  This is a survey of work in the area of dual-porosity modeling. It covers modeling of single phase and two-phase flows by physical arguments and by mathematical arguments of homogenization. It also covers computational schemes and presents some numerical results.
  
  
• T. Arbogast, J. Douglas, Jr., and U. Hornung, Derivation of the double porosity model of single phase flow via homogenization theory, SIAM J. Math. Anal., 21:823-836, 1990.
  
  A general form of the double porosity model of single phase flow in a naturally fractured reservoir is derived from homogenization theory. The microscopic model consists of the usual equations describing Darcy flow in a reservoir, except that the porosity and permeability coefficients are highly discontinuous. Over the matrix domain, the coefficients are scaled by a parameter epsilon representing the size of the matrix blocks. This scaling preserves the physics of the flow in the matrix as epsilon tends to zero. An effective macroscopic limit model is obtained which includes the usual Darcy equations in the matrix blocks and a similar equation for the fracture system that contains a term representing a source of fluid from the matrix. The convergence is shown by extracting weak limits in appropriate Hilbert spaces. A dilation operator is utilized to see the otherwise vanishing physics in the matrix blocks as epsilon tends to zero.
  
  
• P. J. Paes Leme, J. Douglas, Jr., T. Arbogast, and N. P. Nunes, A tall block model for immiscible displacement in naturally fractured reservoirs, In Proceedings, Society of Petroleum Engineers Latin American Petroleum Engineering Conference, Rio de Janeiro, Brazil, October 15-19, 1990, SPE 21104.
  
  A model for two-phase immiscible, incompressible flow in a naturally fractured reservoir is derived via the technique of homogenization for a new geometry for the fractures. Two families of parallel planes normal to the bedding of the reservoir describes the fracture geometry. A finite difference scheme is introduced to approximate the solution of the model. The method involves a local-in-time linearization of the water and capilary potentials on the blocks in a fashion that allows an implicit treatment of the fracture potentials simultaneously with the matrix block boundary potentials. The algorithm leads to natural parallelization of the block calculations.
  
  
• T. Arbogast and F. A. Milner, A finite difference method for a two-sex model of population dynamics, SIAM J. Numer. Anal., 26:1474-1486, 1989.
  
  An explicit finite difference scheme is developed to approximate the solution of a nonlinear and nonlocal system of integro-differential equations that models the dynamics of a two-sex population. The scheme is unconditionally stable. The optimal rate of convergence of the scheme is demonstrated for the maximum norm. Results from a numerical simulation of U.S. population growth from 1970 to 1980 are presented; these compare favorably with the actual data.
  
  
• T. Arbogast, On the simulation of incompressible, miscible displacement in a naturally fractured petroleum reservoir, R.A.I.R.O. Modél. Math. Anal. Numéer, 23:5-51, 1989.
  
  A double porosity/permeability model is presented to simulate an incompressible, miscible displacement in a naturally fractured petroleum reservoir. Fluid flow is described in the fracture system and in each matrix block by defining for each an elliptic pressure equation and a parabolic concentration equation. The matrix/fracture fluid transfer is modeled by imposing boundary conditions on the matrix equations and by including a macroscopically distributed source/sink in the fracture equations. A finite element procedure is defined to approximate the solution. It uses mixed methods for the pressure equations, a modified method of characteristics for the fracture concentration equation, and standard Galerkin methods for the matrix concentration equations. It is shown that the procedure converges asymptotically at the optimal rate.
  
  
• T. Arbogast, Analysis of the simulation of single phase flow through a naturally fractured reservoir, SIAM J. Numer. Anal., 26:12-29, 1989.
  
  A general form of the double porosity model for single phase flow through a naturally fractured reservoir is derived by explicitly considering fluid flow in individual matrix blocks. The Warren and Root model is shown to be a crude approximation to this model. The general model consists of a parabolic equation coupled to a series of parabolic equations. It is shown that the coupling term can be viewed as a positive-semidefinite perturbation of the time derivative, and hence it is verified that the model is well posed. A finite element method is presented to approximate the solution, and optimal order L²-error estimates are derived.
  
  
• J. Douglas, Jr., T. Arbogast, and P. J. Paes Leme, Two models for the waterflooding of naturally fractured reservoirs, In Proceedings, Tenth SPE Symposium on Reservoir Simulation, pages 219-225, 1989, Paper SPE 18425.
  
  Two models are described to simulate the waterflooding of a naturally fractured reservoir. The first is the relatively standard double porosity model in which the matrix is a distributed source to the fracture system. We explicitly describe the influence of the fracture system on the matrix by imposing the proper boundary conditions on the matrix equations. The second model is a limit form of this model for small matrix blocks. It assumes a capillary equilibrium between the matrix and fracture systems. Finite difference procedures are given for these models. These procedures completely separate the matrix block calculations from the fracture calculations. Numerical results are presented to illustrate and contrast the two models.
  
  
• T. Arbogast, J. Douglas, Jr., and J. E. Santos, Two-phase immiscible flow in naturally fractured reservoirs, In M. F. Wheeler, editor, Numerical Simulation in Oil Recovery, number 11 in The IMA Volumes in Mathematics and its Applications, pages 47-66. Springer-Verlag, 1988.
  
  A model is defined to simulate a multidimensional, naturally fractured petroleum reservoir. The imbibition process is correctly modelled as a boundary condition on each matrix block. The model is presented in terms of a saturation and the global pressure of Chavent. The numerical method is based on the use of a mixed finite element method for the pressure and standard Galerkin methods for the saturations in the fractures and the blocks. Optimal order asymptotic convergence of the approximate solution to that of the differential system is established under the assumption of nondegeneracy of the relative permeability functions.
  
  
• T. Arbogast, The double porosity model for single phase flow in naturally fractured reservoirs, In M. F. Wheeler, editor, Numerical Simulation in Oil Recovery, number 11 in The IMA Volumes in Mathematics and its Applications, pages 23-45. Springer-Verlag, 1988.
  
  The double porosity model for single phase flow in naturally fractured reservoirs is consists of a highly coupled system of equations. A parabolic equation describing flow in the fracture system is coupled to the matrix through a distributed source. For each matrix block there is a parabolic problem coupled to the fracture system through a boundary condition. In the standard model, this condition is constant in space. This boundary condition is here generalized to take into account the fracture system's pressure drop across the block by considering its linear variation. It is shown that the resulting model is mathematically well posed. A finite element method is presented to approximate the solution, and optimal order L²-error estimates are derived.
  
  
• J. Douglas, Jr., P. J. Paes Leme, T. Arbogast, and T. Schmitt, Simulation of flow in naturally fractured reservoirs, In Proceedings, Ninth SPE Symposium on Reservoir Simulation, pages 271-279, 1987, Paper SPE 16019.
  
  This paper describes two models for simulating flow in naturally fractured petroleum reservoirs, one for single phase flow of a fluid of constant compressibility, and the other for two-phase, incompressible, immiscible flow. Both models are based on the dual porosity concept. In each model the flow in an individual matrix block is simulated using the standard equations describing flow in unfractured media, and the matrix/fracture interaction is based on the imposition of proper boundary conditions on the surface of the block. The models are presented in an easily parallelizable form.