Theory Guide

The Theory Guide describes the governing equations and constitutive relations and their numerical solution as formulated in the STOMP simulator. This guide traces the solution of subsurface flow and transport problems from theory to numerical implementation. The theory discussions begin with descriptions of the basal governing equations, which are conservation equations for component mass (i.e., water, air, oil, and salt), thermal energy, and solute mass. These equations that describe flow and transport through multiple-phase subsurface environments are presented in their partial differential form. The constitutive relations describe the dependence of secondary variables to the primary unknowns and complete the specification of the thermodynamic and hydrological states of the subsurface system. Primary variables (primary unknowns) for the governing equations of the STOMP simulator are selected dynamically and are dependent on the hydrologic phase condition of the local subsurface environment. Numerical solution includes the discretization of the governing conservation equations and constitutive relations to a system of nonlinear algebraic equations, and the transformation of these equations into a linear system using Newton-Raphson iteration. Also included are the details of implementing the various types of boundary conditions and an overview of two linear system solvers incorporated into the STOMP simulator. This guide concludes with a synopsis of the source code architecture, thus providing the final link between numerical theory and practical implementation of the solution scheme.

Governing Equations

In broad terms, the STOMP simulator solves coupled conservation equations for component mass and energy that describe subsurface flow over multiple phases through variably saturated geologic media. The resulting flow fields are used to sequentially solve conservation equations for solute transport with radioactive chain decay over multiple phases through variably saturated geologic media. These conservation equations for component mass, energy, and solute mass are partial differential equations that mathematically describe flow and transport through porous media and are collectively referred to as the governing equations. The STOMP simulator has capabilities for modeling subsurface flow and transport over three distinct phases: aqueous, gas, and NAPL. The aqueous phase primarily comprises liquid water with small quantities of dissolved air and oils. The gas phase comprises variable amounts of air, water vapor, and oil vapors. The NAPL phase comprises liquid oil components with the amounts of dissolved air and water assumed to be negligible. An additional conservation equation, termed the salt mass conservation equation, can be solved coupled with the component mass and energy equations to simulate brines and/or surfactants. Because the salt mass conservation equation is solved coupled with the flow and energy transport equations, phase properties can be dependent on salt concentrations. This approach allows modeling of systems with brines and/or surfactants, where the physical properties (e.g., density, viscosity, relative permeability) are dependent on the salt concentration. In this respect, the term salt refers to a dissolved and transported solute whose concentration influences the physical properties of the solvents. Salts differ from solutes because the primary assumption for solution of the solute transport equation is that solute concentrations are infinitely dilute.

Constitutive Relations

Constitutive relations are functions that relate the primary unknowns of the governing equations to the associated secondary variables. Each governing equation is solved for a single variable that is referred to as the primary unknown. The remaining variables in the governing equation are referred to as secondary variables. A closed system requires that all of the secondary variables be computable from the set of primary unknowns. In the STOMP simulator primary variables (unknowns) for a particular governing equation can vary between operational modes and phase conditions; however, the set of primary unknowns is fixed for each operational mode and phase condition combination. Primary unknowns are macroscopic properties that fix the physical state of the system. The thermodynamic state of the system is independent of the path by which the system arrived at the given state; however, the saturation state of the system may be dependent on the prior saturation history. Secondary variables include intensive and extensive properties (e.g., phase saturation, phase relative permeability, porosity, tortuosity, viscosity, density, enthalpy, saturated vapor pressures, vapor mass fractions, and diffusion coefficients). Because secondary variables are often computed in terms of other secondary variables, the calculation order for secondary variables is critical to closing the thermodynamic and hydrologic system. Each operational mode and phase condition combination follows a specified sequence for computing secondary variables from the primary variables.

Primary Variables

Each governing conservation equation in the STOMP simulator is solved for one independent primary variable. For systems with multiple phases, the primary variable for a particular governing equation can vary between nodes within the computational domain, where primary variables are chosen by the simulator according to the local phase condition. This numerical solution scheme is frequently referred to as primary variable switching. The Gibbs phase rule states that the number of independent intensive properties required to fix the intensive state of a system equals the number of components plus two minus the number of phases. The number of independent intensive properties is frequently referred to as the degrees of freedom for the system. The thermodynamic and hydrologic state of a porous media system is specified when the number of independent intensive properties equals the number of degrees of freedom according to the phase rule. The number of solved conservation equations in the STOMP simulator equals the number of independent intensive properties less the number of intensive properties fixed through assumptions.

Phase transitions (i.e., phase appearances and disappearances) are characteristic phenomena of multiple phase systems, which require special numerical techniques to resolve. Historically, phase transitions have been handled by restricting phases from completely disappearing and/or through the application of primary variable switching schemes, where the solved primary variables change with the appearance or disappearance of a phase. The STOMP simulator uses a combination of these two techniques. For all operational modes and phase systems, the aqueous phase is restricted from completely disappearing. The aqueous-phase saturation is restricted to values greater than zero through the application of a vapor-pressure lowering scheme and the soil moisture-retention function. For nonisothermal systems, the aqueous-phase saturation is not artificially restricted to values above the residual or irreducible moisture content, but it is restricted to positive nonzero values. For operational modes involving the gas phase, no restrictions are placed on the appearance or disappearance of the gas phase. Gas-phase transitions are treated by switching the primary variable (unknown) for the air-mass conservation equation. Similarly, for operational modes involving NAPL no restrictions are applied to the appearance or disappearance of NAPL. NAPL transitions are numerically treated with a switch in the primary variable for the oil mass conservation equation.

Numerical Solution

Numerical solution refers to the transformation of the governing conservation equations from partial-differential form to algebraic form, algebraic expression of boundary conditions, linearization of the coupled governing equations and constitutive relations and solution of the linear systems. The governing-conservation equations are discretized following the integrated finite difference, which is locally and globally conserving. This transformation requires that the physical domain be spatially discretized into an orthogonal computational domain which comprises nonoverlapping volumes (nodes). Intrinsic properties are assumed to be uniform over the volume domain and are defined for a node point at the geometric center of the volume. Flux quantities are defined at the geometric center of the surfaces between node volumes and along a direction parallel to the surface normal. Fluxes across node surfaces between neighboring inactive nodes and/or adjacent to the domain boundary are controlled through boundary conditions. Solution of the governing-conservation equations in time requires discretization of the time domain. The method employed in the STOMP simulator is implicit using Euler backward-time differencing.

Mass conservation equations for water, air, and VOC components are nearly identical in form, and therefore result in similar algebraic forms. The conservation equation for energy differs from the mass conservation equations because it has diffusive-dispersive and advective components. Discretization of combined diffusive and advective transport requires donor-cell weighting of the transport components, therefore yielding differ algebraic forms than those for the mass conservation equations. The conservation equation for solute or salt transport is similar in form to that of the energy conservation equation but its discretization uses a different donor-cell weighting scheme, therefore resulting in a separate algebraic form. The expressions that result from discretizing the governing equations are nonlinear algebraic equations.

The system of algebraic equations that include the discretized governing-conservation equations and the constitutive functions are nonlinear. Nonlinearities in the soil-moisture retention functions, relative permeability functions, and physical properties near phase transitions are the primary contributors. Conversion of the algebraic equations from nonlinear to linear form follows the iterative Newton-Raphson technique for multiple variables. The technique typically yields quadratic convergence of the residuals, given sufficiently close estimates of the primary unknowns. Each iteration loop requires the solution of a system of linear equations in terms of the equation residuals. Because only orthogonal grid systems are considered, the system of linear equations will have a block-banded structure.

Code Architecture

The primary design guides for the STOMP simulator have been modularity, computational efficiency, and readability. A modular code architecture is beneficial because of the ease of reading, maintaining, and modifying the algorithms and is essential to the variable configuration source code. Computational efficiency refers to both memory requirements and execution speed. The STOMP simulator has been designed with a variable configuration source code that allows the memory requirements and code algorithms to be partially customized to the computational problem. This approach offers considerable advantages with respect to achieving a computationally efficient code design. Within this source code framework, however, many design choices have been made that affect computational efficiency. Algorithm design often offers options between memory and speed. For example, to lessen memory requirements, a code designer may opt to repeatedly compute commonly used variables. Conversely, execution speed may be increased at the cost of increased memory requirements, by storing commonly used variables after their initial computation. Generally, the approach in the STOMP simulator has been to favor increased memory requirements to gain computational speed. This design approach has been chosen because of current state of computer architecture and capabilities. Because the STOMP simulator has been created as a scientific tool, algorithm readability has been an primary design guide. As a scientific tool, the simulator was never expected to remain unmodified, but rather a constantly changing package of software tools which could be applied to new or more complex problems. This design goal makes readability an essential feature of the code. Code readability has been achieved through an extensive use of comments, a modular design, a large group of common blocks, and minimal subroutine and function arguments.