Bounding of Flow and Transport Analysis in Heterogeneous Saturated Porous Media: A Minimum Energy Dissipation Principle for the Bounding and Scale-Up

Abstract

We apply minimum kinetic energy principles from classic mechanics to heterogeneous porous media flow equations to derive and evaluate rotational flow components to determine bounding homogenous representations. Kelvin characterized irrotational motions in terms of energy dissipation and showed that minimum dynamic energy dissipation occurs if the motion is irrotational; i.e., a homogeneous flow system. For porous media flow, reductions in rotational flow represent heterogeneity reductions. At the limit, a homogeneous system, flow is irrotational. Using these principles, we can find a homogenous system that bounds a more complex heterogeneous system. We present mathematics for using the minimum energy principle to describe flow in heterogeneous porous media along with reduced special cases with the necessary bounding and associated scale-up equations. The first, simple derivation involves no boundary differences and gives results based on direct Kelvin-type minimum energy principles. It provides bounding criteria, but yields only a single ultimate scale-up. We present an extended derivation that considers differing boundaries, which may occur between scale-up elements. This approach enables a piecewise less heterogeneous representation to bound the more heterogeneous system. It provides scale-up flexibility for individual model elements with differing sizes, and shapes and supports a more accurate representation of material properties. We include a case study to illustrate bounding with a single direct scale-up. The case study demonstrates rigorous bounding and provides insight on using bounding flow to help understand heterogeneous systems. This work provides a theoretical basis for developing bounding models of flow systems. This provides a means to justify bounding conditions and results.

Plain Language Summary

Complex materials with heterogeneous spatial distributions govern subsurface flow and contaminant transport [1]. In models, it is difficult to represent accurately these complexities and so modelers use more uniform material representations to represent the actual complex distributions. This paper develops a rigorous mathematical method to select a more uniform distribution that bounds the complex material distribution. Our approach applies principles from classical mechanics, originally published by Lord Kelvin in the mid-1800s. This approach determines a more uniform system that is guaranteed to “bound” the problem. We use the term “bound” to mean that an analysis done using the more uniform model will always predict worse conditions than using the more complex representation, but will approximate the actual flow as closely as possible. All subsurface flow and transport models use material distributions that are less complex that the actual material distributions. Our approach insures that material distributions used in models will produce results that both bound the more complex problem and are as similar as possible; this knowledge provides managers with information to support decisions.

We present two mathematical derivations. The first is a simple derivation that determines a bounding system that is homogenous throughout the entire model area or domain. This distribution is guaranteed to bound the problem, but model predictions may be too conservative. We present a second, more complex derivation, that divides the model area or domain into sub-regions or volumes that are different from each other, but are homogeneous within the sub-region or volume. This second approach allows a more accurate model with results closer to the actual complex system, while still guaranteed to bound the results.

These derivations provide a theoretical basis for models that bound flow systems. This means that models developed using these methods can be defended and verify to bound the actual flow conditions.

1. Introduction

A major difficulty in evaluating groundwater contaminant transport problems is our limited capacity to consider the full effects of subsurface material heterogeneity [2,3,4]. Unfortunately, the methods needed to completely characterize and model heterogeneous systems are impractical or not available [5,6]. Even as characterization methods gradually become more practical, time and costs for complete characterization of heterogeneous systems can still be prohibitive for most field problems. In addition, while computer capabilities are significantly increasing, all models require some simplified representation of subsurface materials. There are two approaches to selecting these simplified material distributions, selecting material distributions that produce a best estimate of the problem with the resulting unknown errors to be either conservative or non-conservative, or selecting material distributions that produce a bounding estimate with the errors guaranteed to be conservative [7]. In environmental contamination problems, selecting a material distribution that provides an appropriate bounding of the problem is useful to decision makers. Such a material distribution must provide an appropriate bounding to the problem, that is it must assure that model results will be estimate higher contamination or earlier contamination arrival than actual estimates of contamination transport (e.g., conservative) while still providing a good estimate of system behavior.

The appropriate scale-up of media properties from small areas measured or sampled in the field to the full model domain while representing the appropriate variability, is a related problem for model development. A significant amount of work has been done on the issue of upscaling [8,9,10,11,12,13,14,15,16,17,18]. Hydrologists need ways to provide modeling properties at larger scales that appropriately represent the numerous smaller scale variations present in natural heterogeneous systems [16,19,20,21,22,23]. While the work presented here is closely related to the scale-up problem, it does not directly extend previous work done in this area. Bounding approaches, such as the one presented in this manuscript, however do offer a method for realistic scale-up of material properties for models. Bounding approaches provide material properties for larger model elements that appropriately represent the more variable and smaller scale heterogeneous elements, but are guaranteed to be conservative. This may or may not meet the needs of the modeler, which in many cases would prefer a best estimate rather than a conservative estimate.

Bounding methods have a common requirement, to provide conservative estimates for the subsurface contamination problem. To be useful, the resulting bounding system should have more uniform material properties; it should be less variable or less heterogeneous than the actual system. We present a method for developing bounding conditions for heterogeneous subsurface flow and transport systems using the minimum energy principle of Kelvin from classical in-viscid fluid mechanics [24,25,26,27].

The Kelvin theorem states that irrotational motions minimize the total squared fluid speed throughout the system [28]. Specifically, among all possible rotational fluid motions, kinetic energy is minimized if, and only if, the motion is irrotational [27,28]. We can use this minimum energy principle by evaluating the presence and role of rotational flow components in heterogeneous porous media flow. Equally important is realizing that a reduction in rotational flow components accompanies any reduction in heterogeneity present in a system until, at the limits, in a completely homogeneous system, the flow is entirely irrotational [29,30]. Accordingly, the reduced rotational flows that result from less heterogeneous bounding system, does in fact appropriately bound the actual (more heterogeneous) system.

We begin with specific definitions for the class of subsurface problems initially addressed, including a direct statement of the assumptions. This provides an involved and careful delimitation of the particular heterogeneous systems we will consider theoretically and provides the required foundations for subsequent analyses. We investigate concepts and features from the minimum energy theorem of Kelvin pertinent to bounding material heterogeneity. We show that application of the minimum energy concepts within the framework our defined class of heterogeneous problems, enables deriving two sets of bounding and scale-up equations.

Our first derivation assumes the same boundary conditions for both the bounding and original heterogeneous systems. This first derivation includes the specific Kelvin minimum energy principle and provides the simplest direct bounding criteria, but yields only a single ultimate scale-up. The second expanded derivation treats differing conditions along flow boundaries between elements or regions. Though involving a more complicated analysis, allowing different conditions along boundaries enables using overall piecewise, less heterogeneous representations for the more heterogeneous system. These piecewise boundary conditions give effectively an unlimited scale-up flexibility, depending only on the discretization of the model domain. Perhaps equally important, unlimited scale-up is a useful evaluation capability for analytic element piece-wise modeling of heterogeneous porous systems [31,32,33].

We close with a case study that illustrates bounding and limited scale-up for a flow domain where we assume the flow boundary conditions are equal, which means that classic Kelvin-type limits apply. This example, though perhaps of less practical use in scaling, provides important insights into the bounding representation derived using the minimum energy theorem. This example introduces the arrival time of subsurface contaminant as a concept. Contaminant arrival time can be employed as an efficient and simple summary for flow and transport effects in field problems.

2. Definitions and Assumptions

We restrict this examination of the bounding of flow and scale-up in heterogeneous porous media to saturated flow of a homogeneous fluid in a heterogeneous permeable media. While these principles are applicable for bounding and scale-up of partially saturated and multi-phase flow systems, we limit our analysis to saturated flow. We also limited this work to homogeneous incompressible fluids, which eliminates any possible confusion between rotational flow components arising from fluid non-homogeneity and those arising strictly from material heterogeneity [34,35].

2.1. Heterogeneous Porous Material

A heterogeneous porous medium has different quantities or characteristics at different locations. If a porous medium is heterogeneous with respect to a property, then that property is functionally dependent on the spatial location. Considering the hydraulic conductivity, K, as heterogeneous, then:

$$\begin{array}{c}K={\mathrm{F}}_1\left(x,y,z\right),\end{array}$$

(1)

where K is the scalar saturated conductivity, expressed in any consistent set of force (F), length (L), and time (T) units, with units of ${\mathrm{L}}^4/\mathrm{FT}$. F₁ is a function of the spatial variables in parenthesis, and that function is at least once differentiable with units of L⁴/FT. The variables x, y, and z are Cartesian spatial coordinates with z always oriented verticality upward, all with units of L.

Similarly, the porosity, p, in a heterogeneous porous material is functionally dependent on the location in space:

$$p={\mathrm{F}}_2\left(x,y,z\right),$$

(2)

where p is the scalar effective porosity of the material, with units of L³/L³ and F₂ is a function of the spatial variables in parenthesis and is once differentiable with units of L³/L³.

Spatial variations of several other properties, including, medium compressibility, thermal conductivity, exchange capacity, etc. are implied by the broad title of heterogeneous porous media. Here, we use the phrase only to imply heterogeneity of hydraulic conductivity and porosity, the two properties that govern flow in porous media.

2.2. Homogeneous Porous Material Defined

If the porous material is homogeneous, then the saturated hydraulic conductivity is constant in space and time, i.e.,

$$K=K_0=constant$$

(3)

and

$$p=p_0=constant$$

(4)

where the subscript 0 indicates the value at the initial time. Equations (3) and (4) are the reduced special homogeneous cases for Equations (1) and (2), respectively.

2.3. Dynamic Equations in Heterogeneous and Homogeneous Systems

We assume that Darcy’s Law is applicable.

Consider the movement of a single homogeneous fluid through a saturated heterogeneous media, we can describe the macroscopic impelling force vector field as the gradient of a scalar energy potential, which defined on a unit volume basis is:

$$\varphi =P+\rho gz,$$

(5)

and using the definition of the heterogeneous hydraulic conductivity (Equation 1) gives as the Darcy dynamic flux vector:

$$\stackrel{\rightharpoonup }{q}=-K\nabla \mathsf{\Phi },$$

(6)

where $\stackrel{\rightharpoonup }{q}$ is the Darcian macroscopic volumetric fluid flux vector in the heterogeneous porous system (${\mathrm{L}}^3/{\mathrm{TL}}^2$); K is the scalar saturated hydraulic conductivity (Equation (1)), (L⁴/FT); $\mathsf{\Phi }$ is the scalar potential energy per unit volume of fluid defined (Equation (5)) in the heterogeneous system and is defined to always be positive in both space and time, ($\mathrm{FL}/{\mathrm{L}}^3$); P is the spatial and time varying fluid pressure in units of (F/L³); ρ is the constant mass density of a homogeneous incompressible fluid, (FT²/L³); g is the gravitational scalar constant, (L/T²); z is the vertical Cartesian Coordinate, (L), $\nabla$ is the mathematical gradient operator, (1/L); and all other terms are as previously defined.

The fluid flux vector defined by Equation (6) for a heterogeneous, isotropic porous media is a classical complex lamellar vector field [28] and so contains additional rotational flow components. Those rotational components are everywhere orthogonal to, and independent of, the strictly translational flux vector. By contrast, the flux vector field, when the media is homogeneous, is always irrotational because no rotational components occur under the reduced special case of constant hydraulic conductivity [29,30]. Specifically, in the homogeneous case, using Equation (3);

$${\stackrel{\rightharpoonup }{q}}_0=K_0\nabla \mathsf{\phi },$$

(7)

where ${\stackrel{\rightharpoonup }{q}}_0$ is the Darcian macroscopic volumetric fluid flux vector in the homogeneous porous system, (L³/TL²); $K_0$ is the homogeneous saturated hydraulic conductivity of the porous material (Equation (3)) (L⁴/FT); φ = P + ρgz is the scalar potential energy per unit volume in the homogeneous porous system, which is always defined to be positive both in space and time (FL/L³); and all other terms are as previously defined.

Since the flux vector field for the homogeneous case of Equation (7) is always irrotational, that flux field warrants careful examination. It provides the potential bounding possibility for the generally more dissipative rotational flows involved in heterogeneous systems as considered in the next section.

3. Homogeneous Representation of Heterogeneous Systems

3.1. Bounding Homogenous Systems

In this section, we show that the total flux squared is minimized for an appropriate homogeneous representation of a heterogeneous permeable system

We consider the possibilities for bounding of a heterogeneous permeable system through use of an appropriate more homogeneous representation using a Kelvin-type minimum energy principle [24,25,26,27] or an extension of that principle. Truesdell and Toupin [28] state, “A simple and profound theorem of Kelvin characterizes irrotational motions as minimizing the total squared fluid speed.” In other words, among all possible rotational fluid motions, the minimum possible kinetic energy is attained if, and only if, the motion is irrotational [27]. We seek, using the Kelvin principle and the specific system definitions above, expressions for the minimum squared total flux over the region as the basis for bounding evaluations of large heterogeneous systems.

A starting identity is:

$$\overrightarrow{q}={\stackrel{\rightharpoonup }{q}}_0+\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right),$$

(8)

Forming the dot vector product on both sides of Equation (8) yields:

$$\stackrel{\rightharpoonup }{q}\cdot \stackrel{\rightharpoonup }{q}={\stackrel{\rightharpoonup }{q}}_0\cdot {\stackrel{\rightharpoonup }{q}}_0-2K_0[\nabla \mathsf{\Phi }\cdot \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)]+\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)\cdot \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right),$$

(9)

Using the following identity:

$$-2K_0[\nabla \mathsf{\Phi }·\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)]=-2K_0[\nabla ·[\mathsf{\Phi }\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)]]+2K_0\mathsf{\Phi }[\nabla ·\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)],$$

(10)

Substitution of Equation (10) into (9) provides:

$$\stackrel{\rightharpoonup }{q}\cdot \stackrel{\rightharpoonup }{q}={\stackrel{\rightharpoonup }{q}}_0·{\stackrel{\rightharpoonup }{q}}_0+\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)\cdot \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)-2K_0[\nabla ·[\mathsf{\Phi }\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)]]+2K_0\mathsf{\Phi }[\nabla ·\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)],$$

(11)

However, in Equation (11) the last “del-dot” term$\left(\mathrm{i}.\mathrm{e}.,\nabla ·\right)$, is the conservation of mass expressions for the fluxes in the heterogeneous and homogeneous systems, respectively. Further since both conservation of mass expressions represent the saturated flow of an incompressible fluid: hence are equal to zero, the entire last term vanishes.

Multiplying Equation (11) by the fluid mass density, ρ and a differential volume, dv, then integrating over the sub-region volume, V, of interest yields:

$$\begin{array}{cc}{{\displaystyle \int }}_V^{}\rho [\stackrel{\rightharpoonup }{q}·\stackrel{\rightharpoonup }{q}]dv& ={{\displaystyle \int }}_V^{}\rho [{\stackrel{\rightharpoonup }{q}}_0·{\stackrel{\rightharpoonup }{q}}_0]dv+{{\displaystyle \int }}_V^{}\rho [\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)·\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)]dv\\ & -{{\displaystyle \int }}_V^{}2\rho K_0[\nabla ·[\mathsf{\Phi }\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)]dv\end{array},$$

(12)

Since velocities are involved in the classical work of Kelvin and more specifically in in-viscid hydrodynamics, some may wonder the reasons for initially using fluxes in porous media rather than fluid pore velocities. In porous media flow, the macroscopic fluid flux is the primary dynamic dissipation variable. It is also associated most directly with a measure of bulk volume and surface area. Therefore, we use the direct dissipation variable with consistent surface and volume measures, namely the Darcian fluid flux for this derivation. This establishes a porous media minimum dissipation energy principle for saturated flow in heterogeneous systems.

There are two minimum energy principles involved in porous systems, one minimizing energy dissipation in terms of the Darcy fluid flux and a second minimizing the kinetic energy using the fluid pore velocity. This is in contrast to the single Kelvin minimum kinetic energy principle from classical hydrodynamics [24]. Each of the energy dissipation integral terms in Equation (12) warrants physical interpretation. Specifically, the left side represents:

$${\mathrm{F}}_{Heterog}^2={{\displaystyle \int }}_V^{}\rho \left(\stackrel{\rightharpoonup }{q}\right)·\left(\stackrel{\rightharpoonup }{q}\right)dv,$$

(13)

the instantaneous integrated fluid flux squared, F²_Heterog (units of FL), or total energy dissipation in the heterogeneous flow system inside the selected volume, V.

The first term on the right:

$${\mathrm{F}}_{Homog}^2={{\displaystyle \int }}_V^{}\rho \left({\stackrel{\rightharpoonup }{q}}_0\right)·({\stackrel{\rightharpoonup }{q}}_0)dv,$$

(14)

is the instantaneous integrated fluid flux squared, F²_Homog (units of FL), or total energy dissipation in the homogeneous flow system within the specific volume, V.

It is useful to express the last volume integral term on the right side (Equation (12)) as a surface integral using Green’s Theorem:

$${\mathrm{F}}_{Boundary}^2=2\rho K_0{{\displaystyle \int }}_S^{}\phi \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)·\left(\stackrel{\rightharpoonup }{n}\right)ds,$$

(15)

where F²_Boundary, (units of FL) is the instantaneous weighted difference between the heterogeneous and bounding homogeneous flux components along, n, the unit normal vector to the boundary surface element, ds, which is integrated over the entire surface, S, bounding the original selected volume, V.

The substitution of Equations (13), (14), and (15) into Equation (12) gives the notational simpler form:

$${\mathrm{F}}_{Heterog}^2+{\mathrm{F}}_{Boundary}^2={\mathrm{F}}_{Homog}^2+{{\displaystyle \int }}_V^{}\rho \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)·\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)dv,$$

(16)

which suggests that the heterogeneous total flux energy dissipation tends to exceed the energy dissipation in the homogeneous system by the last terms on each side of the equation. Specifically, the total energy dissipation in the representative homogeneous system of constant hydraulic conductivity, K₀, plus the positive integral over the volume and minus any energy change across the boundary surface, together then equals the total flux energy dissipation over that same volume in the heterogeneous system.

The last term is very important as it represents the total difference in energy dissipation between the heterogeneous and homogeneous systems, and as such, we can drop it from Equation (16) to obtain the basic Kelvin inequality. As a scalar vector product of the same difference vector, the volume integral necessarily must be:

$${{\displaystyle \int }}_V^{}\rho \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)·\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)dv\ge 0,$$

(17)

Since ρ is the constant fluid density, it is always positive, and the dot product must necessarily satisfy the inequality everywhere in V, which from the definition of the vector dot product gives:

$$\Vert \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)\Vert \Vert \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)\Vert \cos 0={\Vert \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)\Vert }^2\ge 0,$$

(18)

We need to differentiate between the two cases. The first case (Case 1) is where, at a few locations in the volume, V, the magnitude of the vector difference may be zero, but over most of the points in the integral, the volume must be non-zero, yet at least at one point, the magnitude of the difference is positive so Equation (17) is always positive, yielding:

$${\mathrm{F}}_{Diff}^2={{\displaystyle \int }}_V^{}\rho \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)·\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)dv>0,$$

(19)

For the second case (Case 2) everywhere in the integral volume (Equation (18)) the magnitude of the difference vector equals zero, or

$${\Vert \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)\Vert }^2=0 , Everwhere Inside V,$$

(20)

This can only be satisfied when the vectors in the heterogeneous vector field are everywhere inside V, equal to those in the homogeneous vector field; a very special condition. In fact, it represents the mathematically special degenerate case or cases. For any such degenerate cases the integral in Equation (17) would be:

$${{\displaystyle \int }}_V^{}\rho \left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)·\left(\stackrel{\rightharpoonup }{q}-{\stackrel{\rightharpoonup }{q}}_0\right)dv=0,$$

(21)

Equation (21) represents completely degenerate cases for the entire theory derived to this point. We consider the extent to which such cases may be physically significant later.

Except for the degenerate cases, Equation (19) represents the entire general class of saturated flow systems, substituting Equation (19) into Equation (16) and rearranging results in the inequality;

$${\mathrm{F}}_{Homog}^2<{\mathrm{F}}_{Heterog}^2+{\mathrm{F}}_{Boundary}^2,$$

(22)

Physically this expression suggests that the squared integral flux or total energy dissipation in Equation (14) for the homogeneous selected volume, V, is less than the sum of the integral of the squared flux or total in Case 1 within that energy dissipation heterogeneous (Equation (13)) over the same volume, V, plus the boundary in Equation (5) over the integral energy change (Equation (15)) surface, S, which completely encloses the original selected volume, V.

We consider the two possibilities in Equation (22) depending on the term; either (1) the boundary term vanishes or (2) the boundary term is retained. We first treat the vanishing boundary term in the following section that yields the simpler minimum total dynamic energy dissipation principle for porous systems. We present the second derivation in a later section with extended analysis that retains the boundary terms to give more general bounding conditions.

3.2. Minimum Flux Energy Dissipation

In this section, we derive equations defining the minimum flux energy dissipation principle for flow in saturated porous media.

The simplest minimum energy principle for flow and transport in a saturated porous systems results when the boundary surface integral in Equations (22) and (15) is zero, that is:

$${\mathrm{F}}_{Boundary}^2=0,$$

(23)

and Equation (22) reduces to the simplest and specific inequality:

$${\mathrm{F}}_{Homog}^2<{\mathrm{F}}_{Heterog}^2,$$

(24)

which demonstrates the Kelvin-type theorem that states the minimum total energy dissipation in the homogeneous system is always less than the dissipation that occurs for flow in any heterogeneous porous material occupying the same corresponding volume. Such minimum total dynamic energy dissipation ideas may seem more definitely indicated in the expanded inequality form provided by substituting Equations (6) and (7) into Equations (13) and (14) and then combined with Equation (24) to give:

$$Homogeneous<Heterogeneous,$$

$$K_0^2\rho {{\displaystyle \int }}_V^{}\left(\nabla \phi \right)·\left(\nabla \phi \right)dv<\rho {{\displaystyle \int }}_V^{}{K}^2\left(\nabla \mathsf{\Phi }\right)·\left(\nabla \mathsf{\Phi }\right)dv,$$

(25)

In terms of total energy dissipation, Equation (25) states there is less total energy dissipation in a homogeneous flow system than for flow through any heterogeneous porous system occupying that same volume, V. This (Equation (24) or (25)) is the Porous Media Minimum Total Energy Theorem or Principle for the flow of a homogeneous fluid in saturated porous media systems.

Specifically, Equation (25) shows: among all possible rotational fluid motions (i.e., heterogeneous system flows) in porous media, the minimum total energy dissipation is attained if, and only if, the motion is irrotational, that is it is a homogeneous flow system.

A useful specific bounding criteria based upon the minimum total dissipation theorem (Equation (25)) suggests an appropriate bounding homogeneous system should have a total flux energy dissipation that is equivalent to the energy dissipation occurring in the particular heterogeneous porous system. To obtain an expression for the hydraulic conductivity of the bounding homogenous system, we can replace the inequality in Equation (25) with an equality and solve for an effective homogeneous bounding hydraulic conductivity, K_0B, which just balances (i.e., equals) the total energy dissipation of the heterogeneous system. Specifically, Equation (25) provides the bounding homogeneous system with a single hydraulic conductivity, K₀₀, given by:

$$\frac{K_{0B}^2}{K_{00}^2}=\frac{{{\displaystyle \int }}_V^{}{K}^2\left(\nabla \mathsf{\Phi }\right)·\left(\nabla \mathsf{\Phi }\right)dv}{K_{00}^2{{\displaystyle \int }}_V^{}\left(\nabla \phi \right)·\left(\nabla \phi \right)dv},$$

(26)

where K_0B is the hydraulic conductivity for use in the homogeneous bounding system of volume, V, with units (L⁴/FT), and K₀₀ is a convenient constant hydraulic conductivity value used in solving the boundary value problem to determine the potential, Φ, as needed to use in Equation (26), (L⁴/FT).

Equation (26) shows that K_0B, the bounding hydraulic conductivity, is always larger than K₀₀ since the numerator integral for the heterogeneous system is always greater than the homogeneous integral in the denominator (see Equation (25)).

3.3. Special Degenerate Cases

In this section, we present derivations for special degenerate cases where minimum energy dissipation bounding does not apply.

The degenerate special cases in Equation (21), where the minimum flux energy dissipation bounding does not apply, warrant consideration for the insight they provide and more particularly as a potential warning about using degenerate cases. We need to take care about having an undue reliance upon using such very special degenerate mathematical cases as being representative of the entire larger and important class of saturated flow problems.

Though Equation (21) represents the pertinent integral degenerate cases, that integral resulted from Equation (20), in which the magnitude of the difference vector can be zero if, and only if:

$${\stackrel{\rightharpoonup }{q}}_0=\stackrel{\rightharpoonup }{q},$$

(27)

If Equation (21) is identically zero along with Equation (27) combined in Equation (15), then Equation (16) reduces to the equality:

$${\mathrm{F}}_{Homog}^2={\mathrm{F}}_{Heterog}^2,$$

(28)

or expanding Equations (6) and (7), then using Equations (13) and (14) provides:

$$K_0^2\rho {{\displaystyle \int }}_V^{}\left(\nabla \phi \right)·\left(\nabla \phi \right)dv=\rho {{\displaystyle \int }}_V^{}{K}^2\left(\nabla \mathsf{\Phi }\right)·\left(\nabla \mathsf{\Phi }\right)dv,$$

(29)

The question arises; do special flow situations exist in the real world where Equations (27) and (28) (or (29)) are both satisfied? For example, Equation (27) states that the instantaneous constant fluid flux in the homogeneous media represented by the left side of Equation (27) must everywhere equal the instantaneous flux in the heterogeneous system, a seemingly infeasible condition, or at least occurring only under unusual conditions. Such a condition appears particularly difficult to satisfy and emphasizes the very special nature of such flow systems, if they occur in nature.

Further for any such flow system, Equation (28) or (29) shows the homogeneous system cannot bound the heterogeneous system as was previously proved for the general case involving saturated flow systems. Equation (29) shows the total energy dissipation for the homogeneous system is already equal to the total dissipation in the heterogeneous system. So special are such requirements, that any flow systems satisfying these conditions would appear non-representative of the general class of saturated flow problems.

Unlikely as it may appear, an assumed strictly one dimensional transient saturated flow case may be special enough to satisfy the two conditions in Equations (27) and (29). This one-dimensional case provides a simple example to show how these conditions are satisfied. Assume a one-dimensional flow column of uniform cross-sectional area, A, and of total length, L, with the length along the column being denoted by, λ, where 0 ≤ λ ≤ L. Then Equation (27) becomes, using Equations (6) and (7):

$${\stackrel{\rightharpoonup }{q}}_0=-K_0\frac{d\phi \left(t\right)}{d\lambda }=\stackrel{\rightharpoonup }{q}=-K\left(\lambda \right)\frac{\partial \mathsf{\Phi }\left(\lambda f\right)}{\partial \lambda },$$

(30)

where the variables of time, t, and location along the column, λ, inside parenthesis indicates the functional dependence of the variables, K₀, K, and Φ. Also, in Equation (30):

$$K_0\frac{d\phi \left(t\right)}{s\lambda }=a parameter of time, t, only= \varpi \left(t\right),$$

(31)

Because in a homogeneous column the gradient is constant along λ, upon substituting Equation (31) into (30) gives:

$${\stackrel{\rightharpoonup }{q}}_0=\varpi \left(t\right)=K\left(\lambda \right)\frac{\partial \mathsf{\Phi }\left(\lambda f\right)}{\partial \lambda },$$

(32)

which states that anytime, t, on the right side in the heterogeneous column, the product of the variable hydraulic conductivity, K(λ), and the gradient term of the variable potential, $\frac{\partial \mathsf{\Phi }\left(\lambda f\right)}{\partial \lambda }$ is everywhere constant. Note further under the reduced case of steady flow then the product of the two variables is always a constant rather than being a time varying parameter, $\varpi \left(t\right)$. Such unexpected conditions occur here, where two variables each change everywhere along the heterogeneous column, but with one variable increasing and the other decreasing just enough to completely compensate so the product of the two variables are instantaneously constant everywhere along the column. Just how likely of occurring in nature is this very special case of two completely compensating variables, which always have an instantaneous constant product everywhere along the column length?

If there were some fundamental relationship, as for example in the ideal gas law, then such exact compensation of variables is appropriate. But when any and all possible variations in hydraulic conductivity that may occur in nature are involved, such completely compensating conditions that occur both everywhere and always along every column is certainly unlikely if not an altogether unseen occurrence in nature. In addition, such completely compensating conditions enter implicitly through the assumption that flow is strictly one-dimensional an unlikely case in reality.

Substituting the one-dimensional column results above into Equation (28) yields:

$$K_0^2{{\displaystyle \int }}_0^L{[\frac{d\phi \left(I\right)}{d\lambda }]}^2\mathrm{A} d\lambda ={{\displaystyle \int }}_0^{L}K{\left(\lambda \right)}^2{[\frac{d\mathsf{\Phi }\left(\lambda f\right)}{d\lambda }]}^2\mathrm{A}d\lambda ,$$

(33)

and using Equations (31) and (32) in the left and right sides, respectively then integration gives:

$${\mathrm{F}}_{Homog}^{2}column=A{[\varpi \left(t\right)]}^{2}L\equiv A{[\varpi \left(t\right)]}^{2}L={\mathrm{F}}_{Heterog}^{2}column,$$

(34)

which, since the two integrations are identically equal terms, so also are the outer most terms equal. Therefore, the homogeneous and heterogeneous total energy dissipations are always identically equal for the mathematical degenerate cases of one-dimensional transient flow.

The assumed one-dimensional uniform flow constraint is too restrictive to allow rotational flow in the heterogeneous system, so the dynamic energy dissipation (consistent with Equation (28) or (34)) exactly equals the energy dissipation in the associated homogeneous column. Accordingly, for assumed one-dimensional uniform flow under steady and transient saturated cases the porous media minimum energy dissipation theorem for general homogeneous and heterogeneous flow systems does not hold because such assumed one- dimensional flow cases represent the very special mathematical degenerate flow conditions.

This analysis of the mathematical degenerate cases, satisfied by assuming one-dimensional saturated transient flow, suggest how very special and non-representative that assumed linear uniform flow is within the entire class of groundwater flow and transport problems occurring In nature. Such an over-simplified subsurface representation of natural flow and transport are likely inappropriate and can be ignored in field cases.

3.4. Piecewise Homogeneous Individual Elements

In this section, we derive extended bounding in heterogeneous systems through piecewise scale-up with individual element boundary differences.

We analyze differing conditions along boundaries between individually scaled-up elements, and show that both the bounding and scale-up are more general and flexible. Including boundary differences allows using overall piecewise less heterogeneous representations for bounding the entire more heterogeneous system. We can use essentially any range in element numbers, sizes, and shapes. We can more and more closely represent the heterogeneous system as we increase the resolution of the discretization. The basis of this approach is to represent each smaller individual heterogeneous volume element by an appropriate element that is homogeneous in hydraulic conductivity and bounds the heterogeneous volume. Again, we use the reduced total energy dissipation theory involving flow in the more homogeneous systems as the basis to enable bounding of the small scale and more variable heterogeneous system. The extended derivation begins from the boundary integral in Equation (15) and using Equations (6) and (7) to obtain the expanded boundary expressive:

$$F_{Boundary}^2=-2\rho K_0{{\displaystyle \int }}_S^{}\phi K[\nabla \mathsf{\Phi }·\stackrel{\rightharpoonup }{n}]ds+2p{K}_0^2{{\displaystyle \int }}_s^{}\phi [\nabla \phi ·\stackrel{\rightharpoonup }{n}]ds,$$

(35)

Let the two boundary integrals be further defined as:

$$F_{Boundary}^2=K_0{F}_{Heterog}^{Boundary}-K_0^2{F}_{Homog}^{Boundary},$$

(36)

So using some of Equations (35) and (36) to define first:

$$F_{Heterog}^{Boundary}=-2\rho {{\displaystyle \int }}_S^{}\phi K[\nabla \mathsf{\Phi }·\stackrel{\rightharpoonup }{n}]ds=F_{Heterog}^{Pos}+F_{Heterog}^{Neg},$$

(37)

Then from Equation (37) and the definition of the scalar dot product for two vectors gives:

$$F_{Heterog}^{Pos}=-2\rho {{\displaystyle \int }}_S^{}\phi K\Vert \nabla \mathsf{\Phi }\Vert \Vert \stackrel{\rightharpoonup }{n}\Vert \cos \xi {|}_{\frac{\pi }{2}\le \xi <\frac{3\pi }{2}}ds,$$

(38)

which is always positive in sign because of the evaluation limits on the angle ξ between the vector $\nabla \mathsf{\Phi }$ and the normal vector to the surface, S. Similarly, the negative term is:

$$F_{Heterog}^{Neg}=-2\rho {{\displaystyle \int }}_S^{}\phi K\Vert \nabla \mathsf{\Phi }\Vert \Vert \stackrel{\rightharpoonup }{n}\Vert \cos \xi {|}_{-\frac{\pi }{2}\le \xi <\frac{3\pi }{2}}ds,$$

(39)

Starting again with the other terms from Equations (35) and (36) provides the second expression

$$F_{Homog}^{Boundary}=-2\rho {{\displaystyle \int }}_S^{}\phi [\nabla \mathsf{\Phi }·\stackrel{\rightharpoonup }{n}]ds=F_{Homog}^{Pos}+F_{Homog}^{Neg},$$

(40)

and again using Equation (40) to separate the positive and negative parts of the integral gives:

$$F_{Homog}^{Pos}=-2\rho {{\displaystyle \int }}_S^{}\mathsf{\Phi }\Vert \nabla \phi \Vert \Vert \stackrel{\rightharpoonup }{n}\Vert \cos \theta {|}_{\frac{\pi }{2}\le \xi <\frac{3\pi }{2}}ds,$$

(41)

where θ is the angle between the vector, $\nabla \mathsf{\Phi }$, and the normal vector to the surface S. Also:

$$F_{Homog}^{Neg}=-2\rho {{\displaystyle \int }}_S^{}\phi \Vert \nabla \phi \Vert \Vert \stackrel{\rightharpoonup }{n}\Vert \cos \theta {|}_{-\frac{\pi }{2}\le \xi <\frac{3\pi }{2}}ds,$$

(42)

Equations (38), (39), (41), and (42) provide complete expanded individual positive and negative terms making up the boundary integrals for inclusion in the extended energy dissipation inequality. Accordingly, substituting Equation (36) into (22) provides on rearranging:

$$F_{Homog}^2+K_0^2{F}_{Homog}^{Boundary}<F_{Heterog}^2+K_0{F}_{Heterog}^{Boundary},$$

(43)

As compared to the original entire system bounding expression in Equation (24), the new Equation (43) has both additional homogeneous and heterogeneous boundary terms, which through Equations (38), (39), (41) and (42), introduces additional positive and negative effects on both sides of the inequality. So substituting into Equations (43), both (37) and (40), provides:

$$F_{Homog}^2+[F_{Homog}^{Pos}+F_{Homog}^{Neg}]K_0^2<F_{Heterog}^2+[F_{Heterog}^{Pos}+F_{Heterog}^{Neg}]K_0,$$

(44)

which is the complete inequality for the total energy, but Equation (44) does not always assure a minimum energy dissipation principle. Equation (44) contains homogeneous and weighted heterogeneous terms with each potentially containing both positive and negative values. This is in contrast to the totally proved minimum energy dissipation principal in Equations (24) and (25) and the associated unconditional bounding with Equation (26), which follows directly whenever the boundary integral (Equation (22)) vanishes or is zero.

This extended derivation replaces the earlier vanishing boundary constraint of Equation (22) with a different requirement arising from the added positive and potentially negative boundary terms in Equation (44). Specificity, the positive and negative terms in Equation (44) require each situation or particular case be individually tested to find if a minimum energy dissipation principle applies which must apply before any possibilities for bounding are relevant. The requirement that needs to be satisfied to assure a minimum energy dissipation principle for bounding with Equation (44) is:

$${\psi }_{Bound}=[F_{Homog}^{Pos}+F_{Homog}^{Neg}]K_0^2-[F_{Heterog}^{Pos}+F_{Heterog}^{Neg}]K_0>0,$$

(45)

where Equations (19), (41), (42), (38), and (39), together enable the direct specific evaluation of the inequality. Satisfying the inequality in Equation (45) is needed to assure the minimum energy dissipation principal in Equation (44) is valid. We can conclusively show the specific minimum principle by substituting (45) into Equation (44) to obtain:

$$F_{Homog}^2<F_{Heterog}^2,$$

(46)

demonstrating, that for any particular case where ${\psi }_{Bound}$ in Equation (45) is greater than zero, then the minimum energy dissipation principal shown in Equation (44) is in fact always valid.

To accomplish bounding and the scale-up we need to determine the appropriate homogeneous hydraulic conductivity for the element involved. The bounding element conductivity should have the same (equal) total dynamic energy dissipation that would occur within the heterogeneous flow system within that same larger element. Accordingly, either Equation (44) or the shorter form in Equation (43) can be written to show the weighting with K₀:

$$K_0^2{F}_{Homog}+[F_{Homog}^{Bondary}]K_0^2-[F_{Heterog}^{Boundary}]K_0<F_{Heterog}^2,$$

(47)

where the first term in Equation (47) is defined using Equation (14) with Equations (6) and (7) as:

$$K_0^2{F}_{Homog}=\frac{1}{K_0^2}{F}_{Homog}^2=p{{\displaystyle \int }}_S^{}\nabla \phi ·\nabla \phi dv,$$

(48)

Equation (47) provides the bounding conductivity by replacing the inequality with an equality and then solving the resulting quadratic for the bounding element homogeneous hydraulic conductivity, K_0B, that is:

$$K_{0B}=\frac{[F_{Heterog}^{Boundary}]\pm \sqrt{{[F_{Heterog}^{Boundary}]}^2+4[F_{Homog}+F_{Homog}^{Boundary}]F_{Heterog}^2}}{2[F_{Homog}+F_{Homog}^{Boundary}]},$$

(49)

where each of the terms are available in Equations (36) through (42) to allow detailed evaluation of Equation (49). Only a real positive root is an appropriate solution. After K_0B is obtained from Equation (49), then we use that value in Equation (45) by setting K₀ equal to K_0B to assure that ${\psi }_{Bound}$ is in fact greater than zero thereby assuring the minimum energy dissipation principal is appropriate. If an Equation (49) solution provides a real positive root, then evaluating ${\psi }_{Bound}$ after calculating K_0B should usually be quite satisfactory. The importance lies in appropriately verifying, during the analysis, that the inequality in Equation (45) is satisfied after $K_0=K_{0B}$ is placed into ${\psi }_{Bound}$. With that inequality positive, then Equation (43) or (44) provide a completely valid minimum energy dissipation principle for that particular system.

Some discussion and perspective is useful for expanding the derivations of Equations (43) through (46) and (49) in terms of an extended theory for minimum energy dissipation. One condition is that the inequality in Equation (45) must be satisfied in order to assure existence of the minimum energy principle in Equation (43) or (44).

Two somewhat interrelated and yet important extensions can be made beyond the original simpler theory of Equation (24) or (25). The simpler minimum energy dissipation principle involves primarily one final limit or ultimate state, but the extended theory is far wider ranging with a much greater continuing scope for possible application. The extended theory may use the benefits of any one of many possible reduced energy dissipation states made possible with each of the successively smaller rotational flow effects accompanying any less heterogeneity represented within the system. However, the simple minimum energy principle involves only one ultimate homogeneous case having no rotational flow components present in the system. In other words, the first derived simpler energy dissipation principle allows bounding with just a single homogeneous representation for the entire heterogeneous system. Only one such bounding representation is possible because irrotational flow can only occur in the resulting homogeneous system.

Many more possibilities are available as the extended theory opens entire ranges of potentially less rotational flows with many reduced heterogeneity combinations available through piecewise combinations of possible elements used to cover the entire variable system. With that piecewise representation, and the extended energy dissipation theory, comes almost unlimited flexibility represented by an entire range of element numbers, sizes, and shapes for to represent and bound large systems.

We present these statements to encourage further work into these bounding elements through more investigation and specific applications. We do not suggest the derived results have been broadly enough applied to completely describe all potential effects in using the theory. To the contrary, only a few actual applications using the simple original energy dissipation principle for bounding have been implemented.

We present one such example for the simpler bounding approach using the minimum energy dissipation principle in the following section.

4. Case Study

We present a case study of groundwater transport and show how the minimum energy representation provides conservative estimates of travel time for contaminant transport. This case study is relatively simple, allowing the results to be more easily presented and discussed.

4.1. Description

This case study provides insight and effectively illustrates the practical significance of the simpler minimum energy dissipation principle derived above. This example demonstrates homogeneous bounding (Equation (24)) and a simple single system scaling (Equation (26)). It illustrates the minimum energy dissipation principle in use for subsurface flow, which is a first step toward using formal contaminant transport bounding in saturated heterogeneous systems. We present the initial transport results to provide insight and a more confident appreciation for bounding of flow as a vital first step toward bounding transport. This example illustrates that appropriate bounding of flow may occasionally adequately bound advective transport, however this bounding cannot always be assured without both combined fluid flux and follow-up transport bounding.

We consider groundwater flow and advective transport from a covered water pond or seepage crib located some 750 m (2460′) from a river. Groundwater flow is primarily in a shallow unconfined aquifer with a mean saturated thickness of 10.7 m (35.1′) above an underlying clay, which is the effective bottom of the subsurface flow system. Pump test data and previous saturated flow model calibrations provided the partially characterized lateral representation of spatially heterogeneous hydraulic conductivity distribution (Figure 1). The regions of highest hydraulic conductivities of 120 m/day (393.7 ft/day) are located midway along the trench extension from the crib, and diminishing hydraulic conductivity values occur elsewhere with lowest values of 40 to 50 m/day (131–164 ft/day) near the river.

The long-term steady flow rate entering the crib is 3293 liters/min or 870 gal/min. This is smaller than the seepage intake rate over the rectangular crib area of 73.2 m by 76.2 m (240 × 250 ft) which indicates there was no water ponding to raise the water level in the crib. Thus, the water level in the crib was never high enough to over-top the weir crest, so there was never any flow into the trench. All seepage into the groundwater system was through the bottom of the crib.

4.2. Heterogeneous Model and Results

We used the Variable Thickness Transient (VTT) flow model [36], a finite difference code, to develop a calibrated saturated flow model of the site. We selected the VTT model as it is a relatively simple model that provides a good platform for comparison with our bounding solutions. For the model, we used a 1 m (3.048′) grid spacing. We used the calibration process, along with geologic knowledge, to estimate the heterogeneous spatial distribution of hydraulic conductivity shown in Figure 1. Figure 1 shows the high conductivities near the center of the model of 120 m/day, with conductivity decreasing towards to river and towards the right edge of the figure. A zone of relatively high conductivity, 90 m/day, extends from the crib to the top of the figure and below the crib; it extends to the lower right of the figure.

Based on this spatial conductivity distribution, we used the VTT model to compute the steady state flow and the associated advective transport results for the heterogeneous system.

Figure 2 shows the results which include both the contaminant travel time (red lines) and groundwater head contours (blue lines). Head ranges from 121.3 m (398′) near the crib and gradually diminishes to the 118.3 m (388′) near the river. The head contours show a groundwater mound at the crib, with decreasing head in all directions. The regional head gradient is from the lower left of the figure to the upper right, with head contours approximately parallel to the river. This regional pattern can be seen in regions further from the crib that are not shown here.

The red lines in Figure 2 are the path lines that show individual paths of the water transport from the crib to the river. The labels associated with each path line are the elapsed time for a water particle to reach the river, along that pathway, when starting from the crib at time zero. Travel times to the river range from 94 days near the center of the figure, with flow following the shortest route and through zones of higher conductivity, to 975 days in the upper left of the figure, where flow follows the zone of high conductivity (90 m/day) north before bending toward the river. Paths that start to the west (right side of figure) take the longest with travel times over 1000 days. The labeled travel time contain the cumulative effects of the greater energy dissipation occurring due to the rotational flow along each particular path in the heterogeneous flow system.

4.3. Bounding Homogeneous Model and Results

We used Equation (26) to calculate a homogeneous hydraulic conductivity value guaranteed to bound the heterogeneous solution. This homogenous system dissipates the same total amount of energy as in the original heterogeneous system. While some of the energy in the heterogeneous system is dissipated by irrotational flow, the homogenous system results in higher flux, as there is no irrotational flow, insuring conservative results, which in this case are shorter travel times.

We used the PATHS Model [37] because can use analytic solutions to compute results for the homogenous system. The analytic solutions were particularly useful in computing the solutions and evaluating the specific terms required in Equation (26).

The computed bounding homogenous hydraulic conductivity calculated using Equation (26) was 72 m/day (236.2 ft/day). This value represents the homogeneous hydraulic conductivity the provides an energy dissipation in the homogenous system just equal the energy dissipation in the partially characterized heterogeneous aquifer shown in Figure 1. Since no energy is lost to irrotational flow, the total flow is higher, providing a conservative (bounding) solution.

The PATHS Model used an inflow rate to the crib of 3293 L/min (870 gal/min) with a saturated thickness of 10.7 m (35.1 ft). We represented the crib in the PATHS model as a circular area, rather than rectangular, with a radius of 42.12 m (138.2 ft). This circular size gives an equivalent inflow area to match the actual 73.2 × 76.2 m (240 ft × 250 ft) crib. The results, even close around the crib, were very similar to the VTT model. All other boundary conditions were effectively the same in the two models. The boundaries of the bounding homogeneous and heterogeneous systems were effectively equivalent as required in Equation (23) and defined in Equation (15). Accordingly, the total minimum energy dissipation principle in Equation (24) or (25) is satisfied to assure appropriate flow bounding.

4.4. Comparison and Discussion of Model Results

Table 1. Comparison of results for contaminated water from the 1325 Crib travel time to the river for the VTT (heterogeneous model and representation) and PATHS (bounding homogenous model and representation).

Cumulative Flow to River (q/Q0)	Arrival Time at River from 1325 Crib
	VTT Model (Days)	PATHS Model (Days)
0.000	193.58	184.09
0.050	197.24	186.22
0.100	202.71	189.48
0.150	211.85	197.66
0.200	222.80	206.41
0.250	241.07	220.83
0.300	262.98	237.21
0.350	299.51	255.08
0.400	350.64	286.15
0.500	509.52	360.07
0.600	774.33	469.84
0.700	1110.36	634.22
0.800	1519.44	892.33
0.900	2016.18	1363.94
0.963	2534.84	2046.16
0.988	---	2810.42

presents the flow and advective transport results for the heterogeneous system from the VTT model and for the bounding homogeneous system from the PATHS model. The table includes the instantaneous ratios of cumulative contaminated water flow rate, q, entering the river as a ratio to the steady crib inflow rate, Q₀ = 3293 L/min or 870 gal/min, which were the same in both the original heterogeneous and bounding cases.

The relative cumulative flow rate, q/Q₀, of water from the crib entering the river is the first of the useful summary terms for quantifying potential releases to the river in steady state systems. The second is the contaminant travel times (Column 2), along the various specific path travel paths from the heterogeneous system, Column 3 presents the travel times computed by the homogeneous model for comparison. Together the outflow ratios and associated travel times provide the contaminant arrival time distributions or curves, which is a simple response function that completely summarizes the steady flow and transport effects on contamination problems [34,38,39].

Figure 3 summarizes the river arrival model results (Table 1) as water arrival time distributions. The figure shows the arrival curve of the PATHS bounding homogeneous representation is everywhere to the left (earlier arrives) and always is higher (greater outflow) than is the arrival distribution for the heterogeneous or partially characterized porous system. Accordingly, the homogeneous bounding system results always provide worst case estimates, both by giving earlier contaminant arrivals and by greater contaminated water outflow rates into the river.

Equation (24) always assures a consistent vertical difference between the upper homogeneous bounding arrival curve and the lower dashed partially characterized heterogeneous arrival curve (Figure 3) this represents the incremental q/Q₀ value. This shows that the simple minimum energy dissipation principle applied to a heterogeneous system provides a homogeneous system that appropriately bounds the fluid flux, hence the q/Q₀ entering the river.

For arrival times above 300 to 500 days in Figure 3, there is a comfortable horizontal margin in the travel times for bounding between the dashed heterogeneous partially characterized arrival curve and the solid upper homogeneous bounding curve. For shorter arrival times the bounding difference is small, though even the early arrival times are strictly bounding as may be verified in Table 1.

This comparison shows that the model developed using the minimal energy principles bounds the actual flow problem. This provides a means for modelers to better understand groundwater flow behavior and to use these calculations to understand the limits of the flow before collection the data and performing the effort required for a more complex model. In addition, results such as these can help guide development of more complex and accurate models.

5. Conclusions

As in this example, the minimum energy dissipation bounding solution for flow may in some cases also bound the advective transport, but that need not necessarily always do so. Such minimum energy dissipation flow bounding is always the first step required before proceeding to assure the transport bounding through application of the minimum kinetic energy principle.

This example illustrates the effectiveness as a bounding tool of the simplest minimum energy dissipation theorem or principle of Kelvin as derived and applied in this paper to saturated heterogeneous porous systems. Though not illustrated by examples here, the extended theory derived to enable the piece-wise extension of minimum energy dissipation bounding, warrants very careful consideration since likely it has considerably more merits in practical applications.

This study provides a theoretical basis to justify bounding flow models. As such, it is somewhat limited in application. We expect future work to implement these approaches in useable methods and models. For example, combining this approach with analytical element models could result in parameter selections to bound flow conditions that could be defended based on this theoretical work.