WaterSubject

Numerical simulation is commonly used to study flow and transport in porous media in a wide range of applications in environmental sciences and reservoir modeling. The complexity and size of numerical models has increased with the ever increasing speed and memory of computers, accompanied by advances in the numerical algorithms available in the simulation software. The ability to represent hydrostratigraphic details such as thin confining layers and faults has improved with high resolution models [3]. The capability of assessing and controlling numerical truncation error is also greatly improved over early numerical models on coarse grids [3], making it possible to accurately capture features such as simulated pumping wells. Despite improvements in grid resolution, there are processes that occur at scales too small to be resolved accurately using conventional numerical methods. For example, the adsorption of contaminants on mineral phases, and diffusion of contaminants from clay layers, often occurs at the centimeter scale. All but the largest massively parallel computer codes have little hope in representing properties and processes at these scales.

Various approaches have been taken to represent heterogeneities and transport occurring at small-scales. For example, the field of stochastic hydrology attempts to quantify the impact of small-scale permeability heterogeneity by casting the governing equations for hydrologic flow and solute transport as a function of random variables [15]. The resulting set of stochastic partial differential equations can be solved using several simplifying assumptions to provide insight on the impact of permeability heterogeneity on the variable of interest, such as hydraulic head and concentration [15]. Another common approach is to investigate small-scale heterogeneity using random fields of hydrologic properties in a Monte Carlo computational setting [15]. Stochastic simulation and Monte Carlo analyses have been very useful in elucidating fundamental mechanisms and effects of heterogeneities on flow and transport behavior. However, these approaches are still difficult to implement in large-scale, site-specific numerical model studies, which most often still employ uniform properties within a hydrostratigraphic unit.

The dual porosity formulation is another approach to incorporating small-scale heterogeneity into large-scale flow and transport simulations. The original development of dual porosity models began with the work of Barenblatt et al. [1] and Warren and Root [14], who recognized that in fractured porous media, the fractures provide the primary conduit for mass transport, whereas the rock matrix could be represented as a fluid storage medium. In such a system the normal mass and energy balance equations can be written for the fracture domain, and an additional equation is written for the matrix node, which is connected only to its corresponding fracture node in a numerical model. In this formulation, the matrix is not a continuous medium in which the full mass and energy transport equations are solved. Nevertheless, important fluid and energy storage terms within the matrix blocks can be explicitly included.

A significant limitation of the dual porosity concept is that although a matrix node is introduced to capture storage in the medium surrounding each fracture, there is no ability to capture gradients into the matrix with a single node. The solution to the original dual porosity models are sometimes called “quasi-steady” matrix solutions because of this limitation. Zyvoloski [16] extended the treatment of the matrix material by introducing a second matrix node. A related conceptual model, referred to as the multiple interacting continua (MINC) model by Pruess and Narasimhan [9], consists of an arbitrary number of nested continua connected to each other through transfer terms that capture the flow of mass and energy in response to pressure, temperature, and concentration gradients. Furthermore, grid cells within any continuum may or may not be connected to one another. If they are, mass and energy flow occurs in parallel within each continuum. Flow within each continuum was recognized to be important in some unsaturated modeling studies such as transport in unsaturated, fractured tuff [10] because in such media, parallel flow could occur at vastly different velocities in the fractures and surrounding matrix rock. In practice, computational limitations have limited the number of continua employed when simulating large-scale flow and transport processes. Usually, one fracture and one matrix continuum are assumed. Even if computational resources are not an issue, the notion of several nested continua, each of which allows for large-scale heat and mass transport, is a significant model abstraction that is difficult to visualize. Nevertheless, systems with two continua, sometimes called dual permeability models, have been proposed in fractured rocks and soils [4].

As implied by the preceding discussion, the different potential options introduced by the MINC formulation and its predecessors has resulted in a somewhat confusing nomenclature for such models. For this reason, Lichtner [7] introduced a terminology to distinguish models with large-scale transport within each continuum from those in which only one-dimensional transport occurs locally within the matrix adjacent to the corresponding fracture grid cell. Lichtner calls the latter dual porosity, discontinuous matrix (DPDM) models, and the former dual porosity, continuous matrix (DPCM) models. We do not adopt Lichtner’s terminology here because we propose using dual porosity formulations for both fracture/matrix systems and also porous media with heterogeneous materials of contrasting permeabilities. Instead, we introduce the term “primary porosity” to represent the medium in which large-scale global flow and transport occurs, and “secondary porosity” to represent the storage volume, typically of lower permeability, that is connected locally only to the primary porosity. Nevertheless, we note that the model introduced here, which we call a generalized dual porosity model (GDPM), falls in Lichtner’s DPDM category.

This paper introduces an efficient numerical technique for systems containing more than one matrix node. To achieve computational efficiency, a one-dimensional decomposition of the fracture–matrix system is performed that only slightly modifies the Newton–Raphson solution of a single-continuum transport system, and results in significant computational and memory savings. The formulation relies on unstructured grid connectivity, commonly associated with finite element methods, to add connections from the primary grid to the first level of secondary nodes. The technique is independent of the grid representing the primary node set, thereby allowing for single-continuum or dual-continuum simulations to be easily converted into a generalized dual porosity solution with no modification of the existing grid.

An additional attribute of the method is that the number of secondary nodes associated with each primary node is spatially variable. Therefore, areas in which gradients in the secondary porosity are significant can be assigned numerous nodes, whereas areas in which interactions between the primary nodes and secondary nodes are not important need not have as many (or even any) nodes in the secondary node set. The method is independent of the physics of the problem. The present study demonstrates the method for unsaturated flow, coupled heat and flow, and solute transport problems. This differs from models such as Karimi-Fard et al. [6] in which transfer functions are tailored to specific flow physics to develop upscaled flow fields. To our knowledge, this is the first implementation of a spatially variable, dual porosity method, as well as the first time that unstructured connectivity has been used to enhance the utility and applicability of the method. In the examples section, the advantages of the approach are demonstrated.