Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Adv Water Resour. Author manuscript; available in PMC 2010 August 1.
Published in final edited form as:
Adv Water Resour. 2009 August 1; 32(8): 1121–1142.
doi:  10.1016/j.advwatres.2009.05.010
PMCID: PMC2860156

Thermodynamically Constrained Averaging Theory Approach for Modeling Flow and Transport Phenomena in Porous Medium Systems: 7. Single-Phase Megascale Flow Models


This work is the seventh in a series that introduces and employs the thermodynamically constrained averaging theory (TCAT) for modeling flow and transport in multiscale porous medium systems. This paper expands the previous analyses in the series by developing models at a scale where spatial variations within the system are not considered. Thus the time variation of variables averaged over the entire system is modeled in relation to fluxes at the boundary of the system. This implementation of TCAT makes use of conservation equations for mass, momentum, and energy as well as an entropy balance. Additionally, classical irreversible thermodynamics is assumed to hold at the microscale and is averaged to the megascale, or system scale. The fact that the local equilibrium assumption does not apply at the megascale points to the importance of obtaining closure relations that account for the large-scale manifestation of small-scale variations. Example applications built on this foundation are suggested to stimulate future work.

Keywords: Averaging theory, TCAT, System scale, Model formulation

1 Introduction

Previous manuscripts in this series have laid out, and given examples of, the thermodynamically constrained averaging theory (TCAT) for modeling flow and transport in porous media. The objective of the TCAT approach is to increase the length scale of a model from the pore scale in a rigorous fashion. The method integrates conservation equations and thermodynamic relations over a region of interest to obtain revised equations that are posed at the scale of interest. Thus far, we have provided an overview of the elements of the TCAT approach [13], presented some mathematical identities that are useful when implementing TCAT [23], applied the TCAT formalism to obtain the equations for single-fluid-phase flow in porous media at a length scale on the order of tens to hundreds of pore diameters [14], laid out the additional considerations of importance for multi-species systems in which dispersive processes are operative [24], derived equations for species transport from the perspective of momentum equations for each phase and interface as well as the perspective of each species in the phase and interface [16], and considered additional complexity that arises when modeling multiple-fluid-phase flow due to the interactions between the phases and the deformable interface between those phases [19]. In the applications to particular porous medium systems, one consistent objective has been to obtain macroscale equations for which spatial as well as temporal variation of the upscaled variables are considered. The equations developed were closed by constitutive relations derived at the scale of the model.

In many cases, it is useful to model systems at reduced dimensionality. For example, flow in groundwater aquifers may often be considered to be a two-dimensional problem modeled after integration through the thickness of the flow domain [3, 26]. Many instances of flow and transport in a laboratory column are modeled as one-dimensional problems. Although one can integrate the full three-dimensional equations developed using TCAT through one, two, or all three spatial dimensions to reduce the dimensionality of the problem, this approach also involves integration of approximate closure relations to produce lower dimensional counterparts. An alternative is to formulate the lower-dimension equations directly at the larger scale and then develop closure forms directly at the scale of interest. Theorems exist that allow the formulation of a model at the microscale (i.e., the pore scale) or macroscale (i.e., a length scale of tens to hundreds of pore diameters) in some dimensions while being fully integrated over space in the remaining dimensions [18], i.e. formulated at the megascale in those dimensions. Such a mixed approach has been employed in one analysis of open channel flow which involved study of only a single phase [11]. Here we will examine full integration to the megascale over all dimensions of conservation and thermodynamic equations for the case of single-fluid-phase flow in porous media. This will be accomplished in the context of the formalism prescribed for the TCAT approach as discussed previously [13, 23]. We will make use of averaging theorems that eliminate all spatial derivatives such that the models are posed in terms of average variables representative of the entire system and in terms of fluxes at the boundary of the region [18].

We note in passing that the region of integration may be some subsection of a full physical region of interest. Then application of the approach to a collection of subsections and matching the fluxes at the boundaries between subsections provides a model of the entire system. Thus a domain may be modeled as a collection of uniform regions that are linked together. Alternatively, a portion of a region modeled at the megascale can be linked to adjacent subregions that are modeled at the macroscale, or even the microscale, as desired. Thus the flexibility of being able to make use of different length scales in different subregions of space presents an opportunity to optimize the study of a region based on computer resources, data availability, and the features of the region to be modeled. Therefore, the approach of megascale modeling within the TCAT framework has applications to a range of porous medium problems as well as problems in which no solid phase is present.

In the present work, the elements of the TCAT approach will be reviewed insofar as they need to be revisited to formulate the megascale flow and mass transport equations. The result of the analysis will be closed equations that can be used to model single-fluid-phase flow in porous media at the megascale.

2 System Description

The system under consideration in this work is a single-phase flow of a fluid, designated as the w phase, in a deforming porous medium, where the solid is designated as the s phase, in which species transport will be considered unimportant. Thus the equations obtained will be used to model the dynamics of an entity, but we will not consider transport of species. Separating the two phases is an interface, denoted as the ws interface, for which conservation and thermodynamic equations will also be developed. This system is the same as that considered in Gray and Miller [14]. However, rather than integrating microscale equations over an averaging volume with a length scale much larger than the pore scale but much smaller than the system scale, integration will be performed over the entire domain. Thus spatial variability is integrated out of the system, and the resulting conservation equations are ordinary differential equations in time. The spatial integration of point equations to obtain global equations is not particularly novel, but the use of a constrained entropy inequality to obtain constitutive forms at the system scale, or megascale, is novel. The TCAT approach has yielded macroscale equations that are rigorously defined and consistent with the microscale precursor conservation equations and thermodynamic principles. These same qualities are sought in the current work for the target models at the macroscale.

For the present study, the domain is denoted as Ω and its volume is V. Within this domain is a region occupied by a solid phase, Ωs with a volume of magnitude Vs; a region occupied by the fluid phase, Ωw with a volume of magnitude Vw; and the interfacial region between the solid and fluid phases, Ωws with an area of magnitude Aws. These regions are referred to generically as entities. Microscale conservation equations for mass, momentum, and energy will be written for each entity. The external boundary of domain Ω is denoted as Γe. The part of this boundary that intersects the w entity is Γwe and the part that intersects the s entity is Γse. Furthermore, the interface between phases will intersect the boundary of Ω along a collection of common curves designated Γwse. Although it is possible to work with species based momentum and energy equations, this will not be considered in the present exposition. Furthermore, microscale entropy balance equations will be employed for each entity along with appropriate thermodynamic relations. All these equations will be integrated over the corresponding entities to obtain megascale equations describing the system of interest. At this larger scale, the entropy equations are combined to eliminate the terms describing exchanges of entropy across the ws interface and to obtain the entropy balance for the system. The closed models obtained will account for the dynamics of the two phases and the interface between the phases.

Some shorthand notation will be used that makes the expression of the lengthier equations that arise more compact. The set of entities is denoted as


where Ωι represents the domain of the ι entity within Ω, x2110 = {w, s, ws} is the index set of entity qualifiers or identifiers, and w, s, and ws are specific qualifiers that indicate the fluid, or wetting, phase, the solid phase, and the interfacial area between the fluid and solid phases. In instances when we wish to work only with phase entities, we will denote the index set for the phases as x2110P = {w, s}. The index set of the interfaces for this system consists of a single entity and is denoted as x2110I = {ws}. With this convention, x2110 = x2110P [union or logical sum] x2110I.

A compact notation is also used to identify the connected set for an entity, that is the entities that are in contact with a particular entity. The connected sets for the w, s, and ws entities are, respectively, εcw = {Ωws}, εcs = {Ωws}, and εcws = {Ωws} with corresponding connected entity index sets x2110cw = {ws}, x2110cs = {ws}, and x2110cws = {w, s}.

In the approach employed here, the length scale of Ω is taken to be much larger than the pore scale. Thus, we will define average quantities for the system. Because these averages apply to the system as a whole, they will not have variability with respect to the length scale of the averages. This distinguishes the current analysis from that in previous TCAT studies. Nevertheless, the integrals of microscale properties that certainly may vary within the system domain will appear in the megascale equations; and accounting for these variations presents a significant challenge in obtaining closed equations at the scale of interest. The final equations will not contain explicit information about how the microscale entities are distributed within Ω, but the constitutive results will necessarily rely on some degree of specification of the distribution.

Though mathematically distinct from the previous TCAT studies, the approach employed here follows the same route as used in earlier derivations [14, 16, 19, 24]. The existence of these earlier works allows the present contribution to concentrate on the new features incorporated into TCAT analysis here. In general, the formalism requires mathematical theorems that facilitate the change in scale; conservation equations, an entropy balance, and thermodynamic relations at the larger scale; combination of these equations into an entropy inequality that is constrained by the conservation laws; and development of closure relations for the conservation equations from this inequality. These will be developed subsequently.

3 Averaging Operator and Integration Theorems

The theorems employed to integrate the conservation equations for the entities in this system are the standard divergence and transport theorems in three and two dimensions. The primary important feature is to designate the boundaries of each entity explicitly as divided between those that are within the system (i.e., the interphase boundaries at the ws interface) and those that occur at the exterior boundary of the system. In the following theorems, we will make use of the definitions of integral operators as introduced previously [23] whereby


where Pi is a property to be averaged to the megascale, and the subscripts on the operator correspond, respectively, to the domain of integration of the numerator, the domain of integration of the denominator, and a weighting function applied to the integrands in the definition of the averaging process. Omission of the third subscript on the averaging operator implies a weighting of unity, W = 1.

Shorthand notation will be employed to indicate various special forms of the averages that arise. The first of these is the average of the property of an entity over that entity, denoted with a superscript such that


For cases when the mass density of the entity is used as a weighting function for averaging, we denote the average with an overbar on the superscript such that


In some instances, we will average the property of an entity over the domain of another entity. For example, the property of a phase entity can be averaged over the interface that bounds the entity. This is designated by using a subscript to identity the entity to which the property belongs and a superscript to identify the region of integration such as




Because the objective of this work is to develop megascale models that are consistent with established microscale conservation and balance equations and thermodynamics, theorems are needed to support rigorously the change of scale that is needed. In previous TCAT papers where equations were developed at the macroscale, the averaging region was fixed. Here, we integrate over an entire region, which may change with time. Therefore, the theorems used to facilitate the change in scale are integration forms rather than averaging forms. Because the microscale entities of concern include phases, which are inherently three-dimensional objects, and an interface, which is a two-dimensional object, theorems are needed to convert both three- and two-dimensional operators to the megascale. Since the megascale models considered in this work are megascale in all spatial dimensions, the class of theorems needed will convert microscale differential operators to entirely megascale forms. The theorems needed have been derived previously and we follow the naming convention established by Gray et al. [18]. The theorems needed to convert three-dimensional differential operators to the three-dimensional megascale form include the divergence theorem for a microscale phase quantity of the form

Theorem 1 (D[3,(0,0),3])


where the microscale spatial vector function fι is defined, continuous and differentiable in Ωι and nι is a unit outward normal vector from the boundary of Ωι.

The gradient theorem for a microscale phase quantity is

Theorem 2 (G[3,(0,0),3])


where the microscale spatial scalar function fι is defined, continuous and differentiable in Ωι, and nι is a unit outward normal vector from the boundary of Ωι.

The transport theorem for a microscale phase quantity is

Theorem 3 (T[3,(0,0),3])


where fι is continuous in time, t, vws is the velocity vector for the ws interface, and vext is the velocity of the exterior portion of the boundary of ι, Γιe.

Eqn (7)Eqn (9) are presented in the specific context of integration over a region in a multiphase porous medium, but they are the same as standard divergence and transport theorems found in standard texts [e.g., 8, 28]. For the single fluid system being considered here, no common curve exists. Therefore, the divergence and transport theorems for the interface between the two phases have boundaries only at the exterior of the system. The divergence theorem for the surface is

Theorem 4 (D[2,(0,0),3])


where ι [sm epsilon] x2110P, fws is a continuous, differentiable spatial vector function defined in Ωws, [nabla]′· is the surface divergence operator, Γwse is the curve on the system boundary where the Ωws surface intersects the boundary, and nws is a unit vector at the boundary of the domain of interest that is tangent to the ws surface and normal to the boundary curve. Note that nws is not necessarily normal to Γe.

The gradient theorem for the surface is

Theorem 5 (G[2,(0,0),3])


where ι [sm epsilon] x2110P, fws is a continuous, differentiable spatial scalar function defined in Ωws, [nabla]is the surface gradient operator, Γwse is the curve on the system boundary where the Ωws surface intersects the boundary, and nws is a unit vector at the boundary of the domain of interest that is tangent to the ws surface and normal to the boundary curve.

The transport theorem for a scalar microscale property of the interface, fws, is

Theorem 6 (T[2,(0,0),3])


where ι [sm epsilon] x2110P, fws is continuous in time and differentiable and [partial differential]′/[partial differential]t is the partial time derivative with surficial coordinates held constant.

These last three theorems are useful in working with surface conservation equations. Although TCAT and earlier averaging methods for porous media have employed such conservations equations [14, 16, 19, 21] and employed surface theorems [12, 18], the derivation of the theorems remains an area of active interest [5, 9]. These theorems will be applied to the microscale species mass, momentum, and energy conservation equations as well as the entropy balance equation and the thermodynamic expressions for the w and s phases and for the ws interface. The result will be the closed megascale relations that will be used in concert to produce closed megascale models.

4 Conservation Equations

The TCAT formulation makes use of conservation equations at the scale of interest for the derivation of closure relations. For the system under consideration here, the equations to be used are conservation equations for mass, momentum, and energy for the fluid, solid, and fluid-solid interface. They are obtained from their standard corresponding microscale forms by integration making use of the theorems in §3. Although a more general system could be examined that includes species transport, we will not consider that case here because the primary objective is to demonstrate the TCAT approach at the megascale. In this section, we will also develop the entropy balance equation that is employed in determining the constitutive forms. The equations needed are derived in the following subsections.

4.1 Mass Conservation

For the fluid and solid phase entities, the microscale equation of mass conservation takes the form


where ρι is the microscale mass density, and vι is the microscale velocity of the ι entity.

Integration of this equation over the domain of the entity and application of the divergence and transport theorems, Eqn (7) and Eqn (9), yields the megascale form





V is the volume of the domain, Mι is the mass of the ι entity in the domain Ω, and Mwsι is the rate of mass exchange from the ws interface to the ι phase per unit volume per unit time.

Note that Eqn (14) does not contain any spatial derivatives. The first term is the rate of change of total mass of entity ι in the domain Ω, the second term accounts for the flux of mass into the entity ι due to mass exchange at the interface interior to the system, and the third term accounts for the flux of mass out of the system at the exterior boundary.

The mass conservation for the ws interface is a generalization of the usual jump condition at a surface of discontinuity in that mass is also allowed to accumulate in the interface. The general microscale mass conservation equation is written as


The first two terms account for the rate of change of mass at a location in the interface surface and the net outward flux from that point in directions tangent to the surface, respectively. If the interface is massless such that ρws = 0, each of these terms is zero. The third term is the net flux from the adjacent phases to the interface. Eqn (17) is integrated over the entire interface domain within the system, and Eqn (10) and Eqn (12) are applied to obtain


The integral in this expression accounts for flow out of the system of surface mass at the system boundary.

4.2 Momentum Conservation

For a phase entity, the microscale momentum conservation equation is


where tι is the stress tensor and gι is the acceleration due to a body force. Often this acceleration will be taken to be gravity, although the derivation of equations will allow for more general contributions, such as Coriolis forces and electrical effects. Integration of Eqn (19) over the phase entity and application of divergence theorem Eqn (7) and transport theorem Eqn (9) yields the megascale equation





and Twsι accounts for the transfer of momentum from the ws interface to the ι phase due to stress and deviations from mean processes per unit volume per unit time.

For the ws interface, the microscale momentum conservation equation is


When the interface is massless, the first two and the fourth terms in Eqn (23) will be zero. The third term accounts for stress in the interface due to its interfacial tension. The last two terms, involving properties of the adjacent phase entities evaluated at the interface, are the standard terms for the jump condition at the interface. Application of divergence theorem Eqn (10) and transport theorem Eqn (12) leads to the megascale momentum equation for the interface


4.3 Energy Conservation

For the phase entities, the microscale equation for conservation of total energy is

t(Eι+12ριvι·vι+ριψι)+·[vι(Eι+12ριvι·vι+ριψι)]     ·qι·(tι·vι)hιριψιt=0ιP,

where Eι is the internal energy density, ψι is the acceleration potential, qι is the non-advective heat flux density vector, and hι is the heat source density.

Eqn (25) may be integrated over the corresponding phase entity. After application of divergence theorem Eqn (7) and transport theorem Eqn (9) and some rearrangement of the resulting expression, we obtain the megascale total energy equation

ει=ddt(𝕍Eι¯¯+12𝕄ιvι¯·vι¯+𝕄ιψι¯+𝕄ιKEι¯¯)   𝕍[Qwsι+vιws¯·Twsι+(𝕍Eιws¯¯𝕄ι+12vιws¯·vιws¯+ψιws¯)Mwsι]   +Γιenι·[(vιvext)(Eι+12ριvι·vι+ριψι)qιtι·vι]  d𝔯   𝕍ιhι¯¯Ωιριψιtd𝔯=0    ιP,

where the averages with the double overbar are particular forms of quantities obtained using the averaging operator such that



the megascale kinetic energy per unit mass due to microscale velocity fluctuations is


and the term accounting for energy exchange at the interface due to heat exchange and deviation terms is given by

𝕍Qwsι=Ωws[qι+tι·(vιvιws¯)]·nιd𝔯     +Ωws(Eιρι𝕍Eιws¯¯𝕄ι)ρι(vwsvι)·nιd𝔯     +Ωws[12(vιvιws¯)·(vιvιws¯)+ψιψιws¯)ρι(vwsvι)·nιd𝔯.

The microscale equation for energy transport at the interface between phases allows for the possibility of accumulation of energy at the interface. Thus, the equation that accounts for interface energy properties along with the usual interface energy jump condition is

t(Ews+12ρwsvws·vws+ρwsψws)+·[vws(Ews+12ρwsvws·vws+ρwsψws)]    ·qws·(tws·vws)hwsρwsψwst    +[qw+tw·vw(Ew+12ρwvw·vw+ρwψws)(vwvws)]·nw|Ωws    +[qs+ts·vs(Es+12ρs|vs·vs+ρsψws)(vsvws)]·ns|Ωws=0.

Integration of Eqn (31) over Ωws and then making use of divergence theorem Eqn (10) and transport theorem Eqn (12) provides the megascale equation, which can be rearranged to the form


where Aws is area of the ws interface in Ω.

4.4 Entropy balance

The entropy balance equation expresses the rate of change of entropy due to transport as well as generation due to dissipative processes. For a phase entity, the microscale entropy balance is

ηιt+·(vιηι)·φιbι=Λι0   ιP,

where ηι is the entropy density, [var phi]ι is the entropy density flux vector, bι is the entropy source density, and Λι is the entropy production rate density.

This equation can be integrated over its corresponding domain. Then, after using divergence theorem Eqn (7) and transport theorem Eqn (9) and rearranging, one obtains

𝒮ι=ddt(𝕍ηι¯¯)𝕍[Φwsι+𝕍ηιws¯¯𝕄ιMwsι]         +Γιe[ηι(vιvext)φι]·nιd𝔯𝕍ιbι=𝕍ιΛι¯¯0   ιP,




and the term accounting for entropy transport at the interior boundary of entity ι is defined as


The microscale entropy balance for the ws interface is

ηwst+·(vwsηws)bws·φws  +[φwηw(vwvws)]·nw|Ωws  +[φsηs(vsvws)]·ns|Ωws=Λws0.

Eqn (38) is integrated over Ωws, and divergence theorem Eqn (10) and transport theorem Eqn (12) are applied, yielding the megascale entropy balance for the interface


In the entity-based entropy equations, Eqn (34) and Eqn (39), the terms for non-advective entropy transfer at the interfaces appear. These terms make analysis of each entity separately infeasible. Also, the averaging process makes each of these entropy equations a relation for only a portion of the system associated with a particular location. Therefore, the entropy balance that is employed is the summation over all entities. This eliminates the non-advective fluxes at the ws interface and provides

ι𝒮ι=ι[ddt(𝕍ηι¯¯)+Γιenι·[(vιvext)ηιφι]d𝔯]    Ωws(·nw)nw·φwsd𝔯𝕍wbw𝕍sbs𝔸wsbws =𝕍wΛw¯¯+𝕍sΛs¯¯+𝔸wsΛws¯¯=Λ0.

This entropy inequality, constrained by the conservation equations and thermodynamic relations, will be employed directly at the megascale to obtain constitutive closure relations. However, before proceeding, it is necessary to obtain the megascale thermodynamic relations.

5 Thermodynamic Relations

The use of thermodynamic relations at the megascale is complicated by the fact that the local equilibrium approximation may not apply for all variables at this scale. Thus, the direct postulation of thermodynamic relations at the megascale runs the risk of being posed in terms of variables that have no relationship to their well-understood microscale counterparts and, perhaps, lack a clear physical meaning. This situation is not unlike what has been noted by Maugin [22] in his discussion of variables used in rational thermodynamics.

Therefore in the TCAT approach, we start with microscale thermodynamic relations and integrate them to the larger scale. This ensures that all large-scale variables are clearly and uniquely defined in terms of their microscale antecedents. The local equilibrium approximation need not be invoked at the megascale. Because the entities in this study are thermodynamically different, being a fluid, a solid, and an interface, their thermodynamic formulations will be presented separately. The derivation follows along the lines of Gray [10] and is based on classical irreversible thermodynamics (CIT) at the microscale [7], but the divergence and transport theorems provided above are employed instead of averaging theorems. Other microscale thermodynamic formulations may be used as a basis for deriving the megascale thermodynamic relations as well, which should be considered if the models produced based upon CIT fail to describe a physical system of interest.

5.1 Fluid-phase thermodynamics

The microscale CIT expression for the energy of the w phase per unit volume is [1, 4]


where θw is the temperature, µw is the chemical potential, and pw is the fluid pressure. Integration of this equation over the w entity to obtain the megascale expression yields




The partial derivative of Eqn (41) with respect to time is


which can be rearranged to


This equation may be integrated using the transport theorem Eqn (9). Subsequent substitution of Eqn (41) then provides

ddt(𝕍Ew¯¯)θw¯¯ddt(𝕍ηw¯¯)μw¯d𝕄wdt+Ωwspwvws·nwd𝔯    +Γwepwvext·nwd𝔯+Ωw[ηw(θwθw¯¯)t+ρw(μwμw¯)t]d𝔯=0.

The gradient of Eqn (41) can be evaluated to obtain a form similar to Eqn (45)


This equation may be integrated over Ωw and simplified using gradient theorem Eqn (8) giving

Ωwspwnwd𝔯Γwepwnwd𝔯  +Ωw[ηw(θwθw¯¯)+ρw(μwμw¯)]d𝔯=0.

The dot product of this equation with vw and addition with Eqn (46) yields

𝒯w=ddt(𝕍Ew¯¯)θw¯¯ddt(𝕍ηw¯¯)μw¯d𝕄wdt  +Ωwspw(vwsvw¯)·nwd𝔯+Γwepw(vextvw¯)·nwd𝔯  +Ωw[ηwDw¯(θwθw¯¯)Dt+ρwDw¯(μwμw¯)Dt]d𝔯=0,

where the material derivative is defined as


By expressing Tw in this form, we obtain the condition that the terms accounting for non-local equilibrium from a megascale perspective are the integral terms that go to zero when temperature, chemical potential and the velocity differences are zero. The material derivative defined in Eqn (50) with respect to entity velocity vῑ may be transformed to a form relative to the s entity by using the identity


where solid-phase referenced velocities are defined as


Application of Eqn (51) with ι replaced by w to the two material derivatives in Eqn (49) yields

𝒯w=ddt(𝕍Ew¯¯)θw¯¯ddt(𝕍ηw¯¯)μw¯d𝕄wdt    +Ωwspw(vwsvw¯)·nwd𝔯+Γwepw(vextvw¯)·nwd𝔯    +Ωw[ηwDs¯(θwθw¯¯)Dt+ρwDs¯(μwμw¯)Dt]d𝔯    +Ωwηwvw¯,s¯·(θwθw¯¯)d𝔯+Ωwρwvw¯,s¯·(μwμw¯)d𝔯,

Because the system under consideration is formulated at the megascale, θw¯¯=0andμw¯=0. The microscale Gibbs-Duhem equation for the w phase is


From Eqn (8), it follows that


Eqn (53)Eqn (55) can be combined to write

𝒯w=ddt(𝕍Ew¯¯)θw¯¯ddt(𝕍ηw¯¯)μw¯d𝕄wdt  +Ωwspw(vwsvs¯)·nwd𝔯+Γwepw(vextvs¯)·nwd𝔯  +Ωw[ηwDs¯(θwθw¯¯)Dt+ρwDs¯(μwμw¯)Dt]d𝔯.

5.2 Solid-phase thermodynamics

For a solid, the CIT expression for the internal energy makes use of the fact that the solid particle movement can be tracked such that the energy of the s phase per unit volume is


where σs is the Lagrangian stress tensor, [mathematical sans-serif bold C]s is the Green’s deformation tensor, js is the solid-phase Jacobian, and : indicates a double dot product between the tensors σs and [mathematical sans-serif bold C]s. Integration of Eqn (57) over Ωs yields the megascale expression for the internal energy




The partial derivative of Eqn (57) with respect to time is


Introduction of some macroscale variables allows this equation be rearranged to


Integration of Eqn (61) over Ωs, application of transport theorem Eqn (9), and substitution of Eqn (57) then provides

ddt(𝕍Es¯¯)θs¯¯ddt(𝕍ηs¯¯)μs¯d𝕄sdtσs¯¯:ddt(𝕍s𝗖sjs)    +Ωs[ηs(θsθs¯¯)t+ρs(μsμs¯t+𝗖sjs:(σsσs¯¯)t]d𝔯=0.

The gradient of Eqn (57) can be evaluated to obtain a form similar to Eqn (61)


Integration of Eqn (63) over Ωs and use of gradient theorem Eqn (8) in light of Eqn (57) gives


Taking the dot product of Eqn (64) with vs and adding the result to Eqn (62) yields

𝒯s=ddt(𝕍Es¯¯)θs¯¯ddt(𝕍ηs¯¯)μs¯d𝕄sdtσs¯¯:ddt(𝕍s𝗖sjs)    +Ωw[ηsDs¯(θsθs¯¯)Dt+ρsDs¯(μsμs¯)Dt+𝗖sjs:Ds¯(σsσs¯¯)Dt]d𝔯=0.

Following the approach in Gray and Schrefler [17], we can show that

σs¯¯:ddt(𝕍s𝗖sjs)+Ωw[𝗖sjs:Ds¯(σsσs¯¯)Dt]d𝔯=  Ωwsns·[2jsσs:(XxXx)·(vsvs¯)]d𝔯  Γsens·[2jsσs:(XxXx)·(vsvs¯)]d𝔯  Ωs{·[2jsσs:(XxXx)]σs:𝗖sjs}·(vsvs¯)d𝔯  Ωwsns·(vwsvs)σs:𝗖sjsd𝔯Γsens·(vextvs)σs:𝗖sjsd𝔯,

where [nabla]Xx is the derivative of a microscale location on the solid phase with respect to its initial location [8].

Substitution of this identity into Eqn (65) and minor rearrangement provides the megascale thermodynamic expression for the solid phase

𝒯s=ddt(𝕍Es¯¯)θs¯¯ddt(𝕍ηs¯¯)μs¯d𝕄sdt  +Ωw[ηsDs¯(θsθs¯¯)Dt+ρsDs¯(μsμs¯)Dt]d𝔯  Ωwsns·[2jsσs:(XxXx)·(vsvs¯)]d𝔯  Γsens·[2jsσs:(XxXx)·(vsvs¯)]d𝔯  +Ωs{·[2jsσs:(XxXx)]σs:𝗖sjs}·(vsvs¯)d𝔯  Ωwsns·(vwsvs)σs:𝗖sjsd𝔯Γsens·(vextvs)σs:𝗖sjsd𝔯.

5.3 Interface thermodynamics

A surface is a two-dimensional entity and the stress in a surface is due to interfacial tension effects. The microscale thermodynamic expression obtained from CIT for a surface is


where γws is the interfacial tension. Integration of this expression over Ωws within the megascale volume yields


The partial time derivative of Eqn (68) is taken while holding the surface coordinates constant such that


or after introduction of macroscale variables


Integration of Eqn (71) over Ωws and use of transport theorem Eqn (12) and identity Eqn (68) yields

ddt(𝕍Ews¯¯)θws¯¯ddt(𝕍ηws¯¯)μws¯d𝕄wsdtΩws(·nw)γwsvws·nwd𝔯    Γwseγwsvext·nwsd𝔯+Ωws[ηws(θwsθws¯¯)t+ρws(μwsμws¯)t]d𝔯=0.

The surface gradient of Eqn (68) can be evaluated to obtain a form similar to Eqn (71)


Integration of Eqn (73) over Ωws accompanied by use of the gradient theorem Eqn (11) and substitution of Eqn (68) gives

Ωws(·nw)γwsnwd𝔯+Γwseγwsnwsd𝔯    +Ωws[ηws(θwsθws¯¯)+ρws(μwsμws¯)]d𝔯=0.

The dot product of Eqn (74) with vws¯ and addition to Eqn (72) yields

ddt(𝕍Ews¯¯)θws¯¯ddt(𝕍ηws¯¯)μws¯d𝕄wsdtΩws(·nw)γws(vwsvws¯)·nwd𝔯    Γwseγws(vextvws¯)·nwsd𝔯    +Ωwsηws[(θwsθws¯¯)t+vws¯·(θwsθws¯¯)]d𝔯    +Ωwsρws[(μwsμws¯)t+vws¯·(μwsμws¯)]d𝔯=0.

Using the definition


Eqn (75) may be restated as

ddt(𝕍Ews¯¯)θws¯¯ddt(𝕍ηws¯¯)μws¯d𝕄wsdtΩws(·nw)γws(vwsvws¯)·nwd𝔯    Γwseγws(vextvws¯)·nwsd𝔯    +Ωwsηws[Ds¯(θwsθws¯¯)Dt+(vws¯vs¯)·(θwsθws¯¯)]d𝔯    +Ωwsρws[Ds¯(μwsμws¯)Dt+(vws¯vs¯)·(μwsμws¯)]d𝔯=0.

However, spatial derivatives of megascale quantities are zero and the Gibbs-Duhem equation provides


Therefore, Eqn (77) may be simplified to

𝒯ws=ddt(𝕍Ews¯¯)θws¯¯ddt(𝕍ηws¯¯)μws¯d𝕄wsdt    Ωws(·nw)γws(vwsvs¯)·nwd𝔯Γwseγws(vextvs¯)·nwsd𝔯    +Ωws[ηwsDs¯(θwsθws¯¯)Dt+ρwsDs¯(μwsμws¯)Dt]d𝔯=0.

6 Constrained Entropy Inequality

The next step in the development of the TCAT model is to employ the conservation equations and the thermodynamic relations as constraints on the entropy inequality Eqn (40). The constraint equations developed previously (i.e., Eqn (14), Eqn (18), Eqn (20), Eqn (24), Eqn (26), Eqn (32), Eqn (49), Eqn (65), and Eqn (79)) may be employed to write the entropy inequality as


where the coefficients λει,λ𝒫ι,λιandλ𝒯ι are Lagrange multipliers whose values will be selected. Note that the quantities being added to the entropy inequality are all zero and thus the expression for entropy generation provided by entropy inequality Eqn (40) is preserved. The Lagrange multipliers are thus free parameters such that the generation term Λ is unaltered. By selecting the multipliers judiciously, most of the time derivatives in the equation can be eliminated so that the generation term takes on a form that is a collection of force-flux products [7, 20]. Because of the integration over the entire domain to the macroscale, the force-flux relations are those for the entire domain plus those at the boundaries of the region. Those at the boundaries must also be evaluated at the megascale. Thus the generation term is of the general form


where each of the independent driving forces, designated as a scalar F, vector F, or tensor [mathematical sans-serif bold F], along with a subscript, is zero at equilibrium. Additionally, all the conjugate fluxes, designated as a scalar by J, a vector by J, or a tensor by [mathematical sans-serif bold J] will also be zero at equilibrium.

With this in mind, the multipliers that cancel the time derivatives arising in Eqn (80) are





Substitution of these values into Eqn (80) followed by algebraic rearrangement of expressions and some cancellations of terms yields the constrained entropy inequality (CEI) in the form

1θw¯¯Γwenw·(vwvext)[12ρw(vwvw¯)·(vwvw¯)         +ρw(μw+ψwμw¯KEw¯¯ψw¯)]d𝔯   1θs¯¯Γsens·(vsvext)[12ρs(vsvs¯)·(vsvs¯)         +ρs(μs+ψsμs¯KEs¯¯ψs¯)]d𝔯   1θws¯¯Γwsenws·(vwsvext)[12ρws(vwsvws¯)·(vwsvws¯)           +ρws(μws+ψwsμws¯KEws¯¯ψws¯)]d𝔯   Γwenw·[φwqwθw¯¯(vwvext)ηwθw(1θw1θw¯¯)]d𝔯   Γsens·[φsqsθs¯¯(vsvext)ηsθs(1θs1θs¯¯)]d𝔯   Γwsenws·[φwsqwsθws¯¯(vwsvext)ηwsθws(1θws1θws¯¯)]d𝔯   +1θws¯¯Ωws(·nw)nw·(qwsφwsθws¯¯)d𝔯   +(hw¯¯θw¯¯bw)𝕍w+1θw¯¯Ωw[ηwDs¯(θwθw¯¯)Dt+ρwDs¯(μw+ψwμw¯KEw¯¯ψw¯)Dt]d𝔯   +(hs¯¯θs¯¯bs)𝕍s+1θs¯¯Ωs[ηsDs¯(θsθs¯¯)Dt+ρsDs¯(μs+ψsμs¯KEs¯¯ψs¯)Dt]d𝔯   +(hws¯¯θws¯¯bws)𝔸ws   +1θws¯¯Ωws[ηwsDs¯(θwsθws¯¯)Dt+ρwsDs¯(μws+ψwsμws¯KEws¯¯ψws¯)Dt]d𝔯   vw¯,s¯θw¯¯·{𝕍[Twsw+12(vwws¯vs¯)Mwsw+12(vwws¯vw¯)Mwsw]         +Γwenw·twd𝔯+𝕄wgw}   +vws¯,s¯θws¯¯·{𝕍[Twsw+12(vwws¯vs¯)Mwsw+12(vwws¯vws¯)Mwsw]         +𝕍[Twss+12(vsws¯vs¯)Mwss+12(vsws¯vws¯)Mwss]Γwsenws·twsd𝔯𝕄wsgws¯}   +1θw¯¯Γwenw·(tw+pw𝗜)·(vwvs¯)d𝔯   +1θws¯¯Γwsenws·(twsγws𝗜)·(vwsvs¯)d𝔯   +1θs¯¯Γsens·[ts2jsσs:(XxXx)]·(vsvs¯)d𝔯   +1θws¯¯Ωws(·nw)nw·tws·(vwsvws¯)d𝔯+1θw¯¯Ωwsnw·(vwsvs¯)pwd𝔯   +1θws¯¯Ωwsρws(vwsvs¯)·nwnw·gwsd𝔯1θws¯¯Ωws(·nw)nw·(vwsvs¯)γwsd𝔯   +(1θw¯¯1θws¯¯)𝕍[Qwsw+(𝕍Ewws¯¯𝕄w+ψwws¯ψw¯μw¯KEw¯¯)Mwsw]   +(1θw¯¯1θws¯¯)𝕍(vwws¯vs¯)·[Twsw+12(vwws¯vs¯)Mwsw]   +1θws¯¯[(μws¯+KEws¯¯+ψws¯)(μw¯+KEw¯¯+ψw¯)]𝕍Mwsw   +(1θs¯¯1θws¯¯)𝕍[Qwss+(𝕍Esws¯¯𝕄sμs¯KEs¯¯+ψsws¯ψs¯)Mwss]   +(1θs¯¯1θws¯¯)𝕍(vsws¯vs¯)·[Twss+12(vsws¯vs¯)Mwss]   +1θws¯¯[(μws¯+Kwws¯¯+ψws¯)(μs¯+KEs¯¯+ψs¯)]𝕍Mwss   1θs¯¯Ωwsns·[2jsσs:(XxXx)]·(vsvs¯)d𝔯   +1θs¯¯Ωs{·[2jsσs:(XxXx)]σs:𝗖sjs}·(vsvs¯)d𝔯   1θs¯¯Ωwsns·(vwsvs)σs:𝗖sjsd𝔯=Λ0.

This general form does not rely on any approximations or restrictions, except the selection of the microscale form of the thermodynamic relation, CIT. This equation therefore can be used as a general starting point for the analysis of a range of problems. We will consider a restricted set here.

7 Simplified Equation

Eqn (85) can be simplified to obtain an expression useful for extracting closure relations for certain cases. Here, we will impose one set of conditions that are not strongly restrictive and thus does not eliminate a large number of important systems.

As the first condition, we will restrict interest to what is commonly called simple systems [8]. The non-advective heat flux appears as a microscale quantity integrated at the boundary. Thus, we will invoke the condition

φι=qιθι    ι.

This condition indicates that the heat flux divided by the temperature is equal to the entropy flux at the microscale, a condition that allows the entropy flux to be eliminated from the CEI. Because the interface between the phases is two-dimensional, the non-advective fluxes qws (and therefore [var phi]ws) and tws only have components that are tangent to the interface. The energy source term appears at the megascale and for a simple megascale system is related to the entropy source according to

(hι¯¯θι¯¯bι)𝕍ι+1θι¯¯Ωι[ηιDs¯(θιθι¯¯)Dt        +ριDs¯(μι+ψιμι¯KEι¯¯ψι¯)Dt]        d𝔯=0    ιP,


(hws¯¯θws¯¯bws)𝔸ws+1θws¯¯Ωws[ηwsDs¯(θwsθws¯¯)Dt      +ρwsDs¯(μws+ψwsμws¯KEws¯¯ψws¯)Dt]  d𝔯=0.

The portions of these terms involving material derivatives do not appear in the microscale definition of a simple system [8]. By definition, at the microscale, the “average” values would be the point value; and thus the differences between microscale and megascale values in the material derivatives would be zero by definition.

Let us also consider the case where there is no mass exchange between entities such that

Mwsι=0   ιP.

When there is no mass exchange at the ws interface, the microscale condition also applies such that


Also, let the interface be massless such that




Microscale non-advective fluxes in the interface normal to the interface are taken to be zero so that




Eqn (86)Eqn (94) can be used to reduce Eqn (85) to

1θw¯¯Γwenw·(vwvext)[12ρw(vwvw¯)·(vwvw¯)    +ρw(μw+ψwμw¯KEw¯¯ψw¯)]d𝔯  1θs¯¯Γsens·(vsvext)[12ρs(vsvs¯)·(vsvs¯)    +ρs(μs+ψsμs¯KEs¯¯ψs¯)]d𝔯  Γwenw·[qw(vwvext)ηwθw](1θw1θw¯¯)d𝔯  Γsens·[qs(vsvext)ηsθs](1θs1θs¯¯)d𝔯  Γwsenws·[qws(vwsvext)ηwsθws](1θws1θws¯¯)d𝔯  vw¯,s¯θw¯¯·(𝕍Twsw+Γwenw·twd𝔯+𝕄wgw¯)  +vws¯,s¯θws¯¯·(𝕍Twsw+𝕍TwssΓwsenws·twsd𝔯)  +1θw¯¯Γwenw·(tw+pw𝗜)·(vwvs¯)d𝔯  +1θws¯¯Γwsenws·(twsγws𝗜)·(vwsvs¯)d𝔯  +1θs¯¯Γsens·[ts2jsσs:(XxXx)]·(vsvs¯)d𝔯  1θws¯¯Ωwsnw·(vwvs¯){(·nw)γwspwns·[2jsσs:(XxXx)]·ns}d𝔯  +(1θw¯¯1θws¯¯)[𝕍Qwsw+Ωwsnw·(vwvs¯)pwd𝔯+𝕍Twsw·(vwws¯vs¯)]  +(1θs¯¯1θws¯¯){𝕍QwssΩwsns·(vwvs¯)ns·[2jsσs:(XxXx)]·nsd𝔯    +𝕍Twss·(vsws¯vs¯)}  1θs¯¯Ωwsns·[2jsσs:(XxXx)]·𝗜·(vsvs¯)d𝔯  +1θs¯¯Ωs{s·[2jsσs:(XxXx)]σs:𝗖sjs}·(vsvs¯)d𝔯=Λ0,

where the surface identity tensor is defined as

𝗜=𝗜nιnι    ιP.

We will make use of microscale constitutive relations for the fluid and the interface stress tensor. For the fluid let


For the interface, assume that the stress is fully accounted for by the interfacial tension such that


Since the system of interest is a porous medium, we will consider the case where the solid velocity is very small and the solid behaves as a quasi-equilibrium system. Thus we make use of the microscale constitutive relation for the solid phase


with the off-diagonal elements of this tensor being zero at the ws interface such that

ns·ts·𝗜=0    onΩws.

Because the magnitude of the solid phase velocity is very small relative to other dynamic motions, we will eliminate all integrals over Ωs containing vsvs in Eqn (95). The velocities of all entities are small so that terms involving the velocity squared can be considered to make a negligible contribution to energy and to the entropy generation relative to other terms. Finally, we will assume that the dynamics of the movement of the Ωws interface are so slow that the equilibrium force balance identity can be applied at all times on the surface (i.e., a semi-equilibrium approximation) with


relative to the other terms in Eqn (95).

Thus, the entropy inequality becomes

1θw¯¯Γwenw·(vwvext)[ρw(μw+ψwμw¯ψw¯)]d𝔯  1θs¯¯Γsens·(vsvext)[ρs(μs+ψsμs¯ψs¯)]d𝔯  Γwenw·[qw(vwvext)ηwθw](1θw1θw¯¯)d𝔯  Γsens·[qs(vsvext)ηsθs](1θs1θs¯¯)d𝔯  Γwsenws·[qws(vwsvext)ηwsθws](1θws1θws¯¯)d𝔯  vw¯,s¯θw¯¯·(𝕍Twsw+Γwenw·twd𝔯+𝕄wgw¯)  +vws¯,s¯θws¯¯·(𝕍Twsw+𝕍TwswΓwsenws·twsd𝔯)  +1θw¯¯Γwenw·τw·(vwvs¯)d𝔯  +(1θw¯¯1θws¯¯)[𝕍Qwsw+Ωwsnw·(vwvs¯)pwd𝔯+𝕍Twsw·(vwws¯vs¯)]  +(1θs¯¯1θws¯¯){𝕍QwssΩwsns·(vwvs¯)ns·ts·nsd𝔯        +𝕍Twss·(vsws¯vs¯)}=Λ0.

The integrals over the various portions of the external boundaries of the system typically involve products of microscale quantities. In fact, these products involve microscale forces and fluxes. It is useful to rearrange these terms so that the integrals can be evaluated. This is accomplished by making two significant, though exact, alterations to the system.

First, although the integrals are over the entire external boundary of an entity, these integrals can conveniently be subdivided into portions of the area that can be more easily evaluated separately. Then, integration over the external boundary of an entity is expressed as the sum of the integrals over each section. With this employed as needed, Eqn (102) becomes

1θw¯¯iΓweΓweinw·(vwvext)[ρw(μw+ψwμw¯ψw¯)]d𝔯  1θs¯¯iΓseΓseins·(vsvext)[ρs(μs+ψsμs¯ψs¯)]d𝔯  iΓweΓweinw·[qw(vwvext)ηwθw](1θw1θw¯¯)d𝔯  iΓseΓseins·[qs(vsvext)ηsθs](1θs1θs¯¯)d𝔯  iΓwseΓwseinws·[qws(vwsvext)ηwsθws](1θws1θws¯¯)d𝔯  vw¯,s¯θw¯¯·(𝕍Twsw+Γwenw·twd𝔯+𝕄wgw¯)  +vws¯,s¯θws¯¯·(𝕍Twsw+𝕍TwssΓwsenws·twsd𝔯)  +1θw¯¯iΓweΓweinw·τw·(vwvs¯)d𝔯  +(1θw¯¯1θws¯¯)[𝕍Qwsw+Ωwsnw·(vwvs¯)pwd𝔯+𝕍Twss·(vwws¯vs¯)]  +(1θs¯¯1θws¯¯){𝕍QwssΩwsns·(vwvs¯)ns·ts·nsd𝔯      +𝕍Twss·(vsws¯vs¯)}=Λ0,



i is a general index, Γιei is a subset of the external boundary Γιe, and x2110Γιe is an index set that includes indexes of all external boundary subsets for the external boundary of entity ι.

Second, we will add and subtract averages over the boundaries to these terms so that the product of forces and fluxes can be written as a macroscale pair plus a pair that is still within the integral

1θw¯¯iΓweΓwei[nw·(vwvext)υoutwei¯][ρw(μw+ψwμw¯ψw¯)nw·τw·nw]d𝔯    1θw¯¯iΓweυoutwei¯[ρwei(μwei¯+ψwei¯μw¯ψw¯)(nw·τw·nw)wei]𝔸wei    1θs¯¯iΓseΓsei[ns·(vsvext)υoutsei¯][ρs(μs+ψsμs¯ψs¯)]d𝔯    1θs¯¯iΓseυoutsei¯ρsei(μsei¯+ψsei¯μs¯ψs¯)𝔸sei    iΓweΓwei{nw·qwqoutwei[nw·(vwvext)υoutwei¯]ηwθw}(1θw1θw¯¯)d𝔯    iΓwe(qoutweiυoutwei¯ηweiθwei¯¯)(1θwei¯¯1θw¯¯)𝔸wei    iΓseΓsei{ns·qsqoutsei[ns·(vsvext)υoutsei¯]ηsθs}(1θs1θs¯¯)d𝔯    iΓse(qoutseiυoutsei¯ηseiθsei¯¯)(1θsei¯¯1θs¯¯)𝔸sei    iΓwseΓwsei{nws·qwsqoutwsei        [nws·(vwsvext)υoutwsei¯]ηwsθws}(1θws1θws¯¯)d𝔯    iΓwse(qoutwseiυoutwsei¯ηwseiθwsei¯¯)(1θwsei¯¯1θws¯¯)𝕃wsei    +1θw¯¯iΓweΓweinw·τw·nwnw·(vextvs¯)d𝔯    +1θw¯¯iΓweΓweinw·τw·𝗜·(vwvs¯)d𝔯    vw¯,s¯θw¯¯·(𝕍Twsw+Γwenw·twd𝔯+𝕄wgw¯)    +vws¯,s¯θws¯¯·(𝕍Twsw+𝕍TwssΓwsenws·twsd𝔯)    +(1θw¯¯1θws¯¯)[𝕍Qwsw+Ωwsnw·(vwvs¯)pwd𝔯+𝕍Twsw·(vwws¯vs¯)]    +(1θs¯¯1θws¯¯){𝕍QwssΩwsns·(vwvs¯)ns·ts·nsd𝔯        +𝕍Twss·(vsws¯vs¯)}=Λ0,







and Lwsei is the length of the ith segment of the external boundary of the ws entity. These quantities are all averages over a portion of the external system boundary. For example, υoutιei¯ is the net (i.e., relative to the system boundary velocity) density weighted average outward normal velocity of entity ι over segment i of the total external boundary that intersects entity ι.

Recall that the development of megascale equations is intended to be useful for systems where smaller scale variations in system properties can be neglected. Consistent with that objective, we can eliminate some of the terms in the preceding equation. We make the assumption that the velocity of the external boundary is determined by the movement of the solid phase. Thus, for purposes of evaluating the integrals, we will neglect the terms [nι·(vιvext)υoutιei¯] as well as terms involving nι·(vextvs¯). The outward flow of solid material at the boundary is also zero so that υoutse¯0. Also, we will neglect the deviation term involving non-advective heat transfer that appears within integrals, nι·qιqoutιe|wei. Finally, two of the terms involving τw are arguably small in comparison to other tems in Eqn (105). The first is the term be integrated over the boundary that multiplies the small velocity difference nw·(vextvs¯). Second, the tangential component of the microscale fluid velocity is considered to be negligible on the boundary so that this term may also be dropped. With all these observations implemented in Eqn (105), the remaining entropy inequality is

1θw¯¯iΓweυoutwei¯[ρwei(μwei¯+ψwei¯μw¯ψw¯)(nw·τw·nw)wei]𝔸wei  iΓwe(qoutweiυoutwei¯ηweiθwei¯¯)(1θwei¯¯1θw¯¯)𝔸wei  iΓseqoutsei(1θsei¯¯1θs¯¯)𝔸sei  iΓwse(qoutwseiυoutwsei¯ηwseiθwsei¯¯)(1θwsei¯¯1θws¯¯)𝕃wsei  1θw¯¯(𝕍Twsw+Γwenw·twd𝔯+𝕄wgw¯)·vw¯,s¯  +1θws¯¯(𝕍Twsw+𝕍TwssΓwsenws·twsd𝔯)·vws¯,s¯  +[𝕍Qwsw+Ωwsnw·(vwvs¯)pwd𝔯+𝕍Twsw·(vwws¯vs¯)](1θw¯¯1θws¯¯)  +[𝕍QwssΩwsns·(vwvs¯)ns·ts·nsd𝔯    +𝕍Twss·(vsws¯vs¯)](1θs¯¯1θws¯¯)=Λ0.

This equation is in the flux-force form suggested by Eqn (81) with eight types of conjugate pairs. The actual number of conjugate pairs depends on the number of sections of the surface of the system over which integrations have been performed. In this implementation of the TCAT approach, there are conjugate scalar and vector pairs, but no tensor conjugate pairs. The set of forces in Eqn (111) that go to zero at equilibrium is


8 Linearized Closure Relations

The TCAT approach leads to the development of an SEI, which in turn constrains the allowable forms of closure relations needed to produce well-posed, solvable models. Such closure relations are not unique in form, or even the order of the approximation. As a result, an infinite number of closure relations exist for a given model. It is both reasonable and common practice to start with a relatively simple set of closure relations and to increase the complexity of these relations as needed to describe physical systems of concern. A natural place to start is with a linear form of the closure relations. The most general linearization of the flux-force pairs is to make each flux a linear function of the full set of forces. Here we will consider a somewhat reduced, but nevertheless rather general, linear coupling. We will consider that there is coupling involving the average velocity and the average frictional drag within the system. Other coupling is assumed to be unimportant, although these assumptions can be examined further in the context of a physical system. Thus the restrictions employed are:

  • The average velocity of the solid phase is essentially equal to the average velocity of the interface between the fluid and solid such that the force vws¯,s¯0 and the conjugate flux of this term is also approximately zero.
  • The heat transferred into the system at the system boundary where the interface intersects the boundary is negligible so that the heat flux at any boundary curve Γwsei is restricted by (qoutwseiυoutwsei¯ηwseiθwsei¯¯)0 and its conjugate force is also negligible.
  • Each member of the set of boundary flow fluxes, υoutwei¯ will depend on the conjugate force at that location as well as the force vwvs. This is expressed as
    where r^υwei is a closure vector.
  • The flux with conjugate force vw,s also depends on the forces conjugate to the boundary flow fluxes, which is of the form
    𝕍Twsw+Γwenw·twd𝔯+𝕄wgw¯=𝗥^w·vw¯,s¯   iΓwer^μwei[ρwei(μwei¯+ψwei¯μw¯ψw¯)(nw·τw·nw)wei],
    where r^μwei is a closure vector.
  • Each of the non-advective heat flux terms for the phases at the boundary will depend only on the conjugate force at that boundary. Therefore
    where kθ is a heat flux closure parameter.
  • The non-advective heat flux terms within the system will depend only on the conjugate force. This is stated as
    𝕍QwssΩwsns·(vsvs¯)ns·[2jsσs:(XxXx)]·nsd𝔯   +𝕍Twss·(vsws¯vs¯)=k^Qs(1θs¯¯1θws¯¯),
    where kQ is an inter-entity heat transfer closure parameter.

Note that if the system is isotropic, a linear dependence of a vector flux on a scalar force or of a scalar flux on a vector force will not exist because the proportionality parameter is a vector. No isotropic coefficient vector exists. For the megascale case, “isotropy” includes interaction with boundary fluxes as well as the actual properties of the materials making up the porous medium system. For a more general form, it is a small matter to include additional coupling between flow induced by temperature differences or energy transport induced by high flow rates. Closed equations may be obtained under the conditions and restrictions identified above.

9 Closed Conservation Equations

Closure relations may be implemented in the megascale conservation equations to provide a set of equations that describes the system. It is important to realize that some approximations will be made in this analysis consistent with the fact that the conservation equations are integrated over the entire system such that variability of properties within the system are treated only through average quantities. Such a description can be useful in many applications or in some regions of interest where the variability of properties is not important to answering a question of interest concerning the system. Here we will provide closed forms of the mass, momentum, and energy conservation equations that can serve as a basis for megascale analysis.

9.1 Closed mass conservation equations

The mass conservation equation for the phases is given in Eqn (14). When there is no mass exchange between phases, this simplifies to

d𝕄ιdt+Γιeρι(vιvext)·nιd𝔯=0   ιP.

Making use of the definition of the outward volumetric flow per unit area given by Eqn (106) and the stipulation of the constraint that the boundary of the system moves with the solid phase, one can express this equation for the w and s phase, respectively, as




Eqn (121) indicates that the solid mass in the system does not change with time while Eqn (120) expresses the change in mass of phase w as a function of the fluxes at the various locations on the boundary where fluxes occur. There are two options for specifying these fluxes. In one case, the flux might be specified directly, for example at a location where fluid is being added or withdrawn at a specified rate. A second option would involve the imposition of the head at the boundary so that the flow out would be due to a head gradient generated by the dynamics of the system. For this type of boundary condition, we can make use of the linear closure relation for υoutwei¯ out defined in Eqn (113). With this expression implemented at some locations, a form of mass conservation Eqn (120) that could be employed is

d𝕄wdt+iΓweυ(ρweiυoutwei¯𝔸wei)   iΓweμρweiK^wei[ρwei(μwei¯+ψwei¯μw¯ψw¯)(nw·τw·nw)wei]𝔸wei   iIΓweμρweir^υwei·vw¯,s¯𝔸wei=0,

where x2110Γweυ is the set of locations where the velocity is specified and x2110Γweµ is the set of locations where the potential is specified with x2110Γwe = x2110Γweυ [union or logical sum] x2110Γweµ. This distinction accounts for the fact that with some systems, at some boundary locations fluxes will be employed while at other locations the flux potential is used. The formulation is capable of handling this mix of conditions. In this equation [K with circumflex]wei is a megascale parameter that may depend on position on the boundary of the system as well as time. This parameter may be interpreted as similar to a hydraulic conductivity, K^Hwei with


where l is a characteristic length such that the difference in potentials divided by this length approximates the normal gradient in potential at the boundary and [sm epsilon]wei is the porosity at the boundary. The parameter r^υwe accounts for coupling at the boundary between an exit, or entrance, velocity and the average velocity within the system. The frictional loss term within the fluid in Eqn (122) may be neglected since most of the losses are due to fluid-solid interactions that account for the difference between chemical plus gravitational potential at the boundary and the average of the potentials in the interior of the system. In that instance, the mass conservation equation for the fluid may be simplified to

d𝕄wdt+iΓweυ(ρweiυoutwei¯𝔸wei)   iΓweμρweiK^wei[ρwei(μwei¯+ψwei¯μw¯ψw¯)]𝔸wei   iΓweμρweir^υwei·vw¯,s¯𝔸wei=0.

Since the interface is considered to be massless and no phase change occurs, mass conservation Eqn (18) for the interface is trivial with all terms being zero.

9.2 Closed momentum conservation equations

The megascale momentum equation for the phases when no mass exchange is taking place reduces from Eqn (20) to

d(𝕄ιvι¯)dt𝕍Twsι𝕄ιgι¯+Γιe[ριvι(vιvext)tι]·nιd𝔯=0    ιP.

The closure relation for this equation for the w phase is given by Eqn (114). Substitution into Eqn (125) when ι = w then gives

d(𝕄wvw¯)dt+𝗥^w·vw¯,s¯+Γweρwvw(vwvext)·nwd𝔯    +iΓwer^μwei[ρwei(μwei¯+ψwei¯μw¯ψw¯)(nw·τw·nw)wei]=0.

The first term in this equation is the rate of change of w momentum. The second accounts for resistance to fluid flow due to frictional interactions of the fluid with the solid. The integral term accounts for momentum added or subtracted by flow across the system boundary, and the summation term accounts for the addition of momentum due to work done at the boundary of the system. Note that the integral term, in effect, involves the velocity squared.

It can be approximated as


where at each section of the boundary


Substitution of Eqn (127) and Eqn (120) into Eqn (126) with the intrafluid frictional stress term neglected then gives

𝕄wdvw¯dt+𝗥^w·vw¯,s¯+iΓweρwei(vwwei¯¯vw¯)υoutwei¯𝔸wei    +iΓwer^μweiρwei(μwei¯+ψwei¯μw¯ψw¯)=0.

When the rate of change of average velocity of the w phase and the advective momentum flux at the boundary are small, the simplest form of the momentum equation is obtained as


To obtain the closed form of the solid phase momentum equation, we start with momentum Eqn (24) for the interface which simplifies for the case of no mass exchange and a massless interface to


Combination of this expression with Eqn (114) yields

𝕍Twss=Γwenw·twd𝔯+Γwsetws·nwsd𝔯+𝕄wgw¯+𝗥^w·vw¯,s¯    +iΓwer^μwei[ρwei(μwei¯+ψwei¯μw¯ψw¯)(nw·τw·nw)wei].

Substitution into Eqn (125) for the solid phase, while making use of the fact that no solid phase enters or leaves the system, then provides the momentum equation

𝕄sdvs¯dt𝕄wgw¯𝕄sgs¯𝗥^w·vw¯,s¯   iΓwer^μwei[ρwei(μwei¯+ψwei¯μw¯ψw¯)(nw·τw·nw)wei]   Γwenw·twd𝔯Γsenw·tsd𝔯Γwsenws·twsd𝔯=0.

This equation provides an expression for the response of the solid phase momentum to gravitational forces, interaction with the fluid in the system, and forces applied at the boundary of the system. For conditions under which the flow is described by Eqn (130) and the acceleration of the solid phase is negligible, this equation may be simplified to a simple balance of forces acting on the solid


9.3 Closed energy conservation equations

The development of the closed energy equation for the w phase begins with energy conservation Eqn (26) simplified to account for the fact that no phase change is being modeled. Additionally, we will neglect the subscale kinetic energy term, KEw¯¯ and assume that the only body force is due to gravity so that the partial time derivative of the microscale gravitational potential is zero. With these terms deleted and the time derivative of the kinetic energy expanded out using the product rule, the energy equation for entity w is

ddt(𝕍Ew¯¯)+vw¯·ddt(𝕄wvw¯)12vw¯·vw¯ddt(𝕄w)+ddt(𝕄wψw¯)    𝕍Qwswvwws¯·𝕍TwswΓwenw·(qw+tw·vw)d𝔯    +Γwenw·(vwvext)(Ew+12ρwvw·vw+ρwψw)d𝔯    𝕍whw¯¯=0.

Next substitute mass conservation Eqn (119) and momentum conservation Eqn (125) for the w entity into this equation to eliminate the time derivatives of mass and momentum. Then rearrange terms to obtain

ddt(𝕍Ew¯¯)+(vw¯vwws¯)·𝕍Twsw    𝕍QwswΓwenw·[qw+tw·(vwvw¯)]d𝔯    +Γwenw·(vwvext)[Ew+12ρw(vwvw¯)·(vwvw¯)]d𝔯    𝕍whw¯¯=0.

Make use of the linearized force-flux relation of Eqn (117) to restate this equation as

ddt(𝕍Ew¯¯)+vw¯,s¯·𝕍Twsw    +Ωwsnw·(vwvs¯)pwd𝔯k^Qw(1θw¯¯1θws¯¯)    Γwenw·[qw+tw·(vwvw¯)]d𝔯    +Γwenw·(vwvext)[Ew+12ρw(vwvw¯)·(vwvw¯)]d𝔯    𝕍whw¯¯=0.

Constitutive Eqn (114) is used to eliminate 𝕍Twsw from Eqn (137) to obtain, under conditions where Eqn (130) holds

ddt(𝕍Ew¯¯)vw¯,s¯·[𝗥^w·vw¯,s¯+𝕄wgw¯]    vw¯,s¯·iΓwer^μwei[ρwei(μwei¯+ψwei¯μw¯ψw¯)]    Γwenw·[qw+tw·(vwvs¯)]d𝔯    +Γwenw·(vwvext)[Ew+12ρw(vwvw¯)·(vwvw¯)]d𝔯    +Ωwsnw·(vwvs¯)pwd𝔯k^Qw(1θw¯¯1θws¯¯)    𝕍whw¯¯=0.

Some additional rearrangement may be made involving the integral terms and using Eqn (97) to obtain

ddt(𝕍Ew¯¯)vw¯,s¯·𝗥^w·vw¯,s¯    vw¯,s¯·iΓwer^μwei[ρwei(μwei¯+ψwei¯μw¯ψw¯)]    +vw¯,s¯·(Γwenwpwd𝔯+Ωwsnwpwd𝔯𝕄wgw¯)    vw¯,s¯·Γwe𝗜·τw·nwd𝔯    Γwenw·[qw+tw·(vwvw¯)]d𝔯    +Γwenw·(vwvext)[Ew+12ρw(vwvw¯)·(vwvw¯)]d𝔯    +Ωwsnw·(vwvw¯)pwd𝔯k^Qw(1θw¯¯1θws¯¯)    𝕍whw¯¯=0.

Terms that are second order in velocity will now be eliminated in light of the porous medium system being studied to obtain

ddt(𝕍Ew¯¯)Γwenw·[qw+tw·(vwvw¯)]d𝔯    +Γwenw·(vwvext)Ewd𝔯+Ωwsnw·(vwvw¯)pwd𝔯    k^Qw(1θw¯¯1θws¯¯)𝕍whw¯¯=0.

The linearized closure relation based on Eqn (115) is then substituted into Eqn (140) at locations where the temperature, rather than the heat flux, is specified at the boundary. This provides

ddt(𝕍Ew¯¯)iΓweqqoutwei𝔸wei+iΓweθk^θwei(1θwei¯¯1θw¯¯)𝔸wei    +Γwenw·(vextvw¯)pwd𝔯+Ωwsnw·(vwvw¯)pwd𝔯    Γwenw·τw·(vwvw¯)d𝔯+Γwenw·(vwvext)pwμwd𝔯    k^Qw(1θw¯¯1θws¯¯)𝕍whw¯¯=0,

where x2110Γweq is the index set of locations on the w external boundary where the flux is specified, x2110Γweθ is the index set of locations on the w external boundary where the temperature of the w phase is specified, and x2110Γwe = x2110Γweq [union or logical sum] x2110Γweθ.

The integrals that remain in Eqn (141) contain some challenges. The first two are related to the change in volume of the system in response to the pressure. The third integral involves frictional effects at the boundary of the system and is typically small. The last term can be managed by breaking the boundary into segments where flow is occurring and then using a closure relation such as


where the macroscale quantities μwwei¯¯ are specified with an equation of state.

For the solid phase, the energy equation analogous to Eqn (136) may be obtained as

ddt(𝕍Es¯¯)+(vs¯vsws¯)·𝕍Twss    𝕍QwssΓsens·[qs+ts·(vsvs¯)]d𝔯    +Γsens·(vsvext)[Es+12ρs(vsvs¯)·(vsvs¯)]d𝔯    𝕍shs¯¯=0.

Substitution of Eqn (118) into Eqn (143) yields

ddt(𝕍Es¯¯)Ωwsns·(vsvs¯)ns·ts·nsd𝔯k^Qs(1θs¯¯1θws¯¯)    Γsens·[qs+ts·(vsvs¯)]d𝔯𝕍shs¯¯=0,

where we have made use of the definition of the system such that the normal component of the solid phase velocity at the boundary defines the velocity of the boundary. Substitution of the closure expression for the non-advective heat flux at the system boundary at locations were the temperature is known and leaving the heat flux explicitly where it is specified gives

ddt(𝕍Es¯¯)iΓseqqoutsei𝔸sei+iΓseθk^θsei(1θsei¯¯1θs¯¯)𝔸sei    Γsens·(vextvs¯)ns·ts·nsd𝔯Ωwsns·(vsvs¯)ns·ts·nsd𝔯    k^Qs(1θs¯¯1θws¯¯)𝕍shs¯¯Γsens·ts·𝗜·(vsvs¯)d𝔯=0,

where x2110Γseq is the index set of locations on the s external boundary where the flux is specified, x2110Γseθ is the index set of locations on the s external boundary where the temperature of the s phase is specified, and x2110Γse = x2110Γseq [union or logical sum] x2110Γseθ.

10 Discussion

The purpose of this series of manuscripts on the TCAT approach is to lay the foundation for rigorous model building for a variety of porous medium systems. The approach that has evolved in this series to accomplish this goal is to present the components of the machinery needed to formulate models and to use this machinery to create frameworks that support a hierarchy of models. Thus our intent with these works is not to present solutions for a single problem, but to establish a series of frameworks that can be applied to entire families of problems, albeit with some additional effort. An important feature of this work is that the model formulation process proceeds in a series of formal, well-defined intermediate steps that provide convenient starting points for additional individual efforts without requiring the substantial foundational efforts needed to derive starting points such as the CEI or the SEI.

Given the primary restrictions on the class of problem of concern in this work and the applicability of CIT as an appropriate thermodynamic theory, the CEI that is derived is an exact expression, which can be broadly applied to single-fluid-phase megascale problems in general. The CEI is not strictly in the desired form of a set of forces and fluxes, which is needed to guide model closure. Because of this, a series of approximations and restrictions are applied to reduce the CEI to the SEI, which is in strictly force-flux form. The simplifications applied are a reasonable starting point, but may not be the final word. Alternative sets of simplifications to achieve force-flux forms can be explored starting from the CEI to produce alternative forms of the SEI if such efforts are deemed necessary to describe a physical system of concern. It should be noted that reducing the CEI to the SEI is much less effort than deriving the CEI, which requires a substantial series of manipulations. The sort of manipulations required to derive the CEI have been detailed in previous papers in this series [15, 16].

The models considered in this work are megascopic in nature, by which we mean that a solution is sought in terms of variables averaged over the entire spatial domain of the system. Because of this, spatial variability within the domain is not resolved at the scale of the solution. It follows then that the theorems applied to derive the conservation equations of concern transform the spatial derivatives to boundary terms, and representing the boundary terms appropriately becomes a central focus of the work.

As an example, and initial starting point for modeling megascale single-fluid-phase systems, a set of linear closure relations are derived from the SEI. These closure relations are in turn used to derive a set of closed conservation of mass, momentum, and energy equations for a megascale system. These closed equations include model parameters that must be determined for application to any specific physical problem, which is typical of any mechanistic model. Because all quantities appearing in the final models are rigorously defined averages of microscale precursors, the means exists to link models across scales and to investigate the dependence of the resultant model parameters on the details of a sub-scale system. Multiscale computation would be one way to approach this problem, while alternative analytical approaches of linking megascale model parameters to microscale system properties are possible under certain, more restrictive conditions. Future works building upon this foundational work along these lines seems potentially fruitful. For example the theoretical underpinnings of Darcy’s [6] foundational experimental work and data correlation can be explored. This can be important in determining problem parameters and also extensions to multiphase systems.

Megascale models are of most applicability when the spatial variations within the system are sufficiently simple in form to result in meaningful megascopic model parameters of a correspondingly simple form. When this is no longer the case, megascale models lose their utility. One way in which megascale models may be used advantageously is to subdivide a system into a set of domains that can each be approximated with a megascopic model and linked together at the boundaries. A second extended way to use megascopic models is to separately treat a simple portion of a larger system with a megascopic model that is coupled to a more highly resolved macroscale model for the more complex portions of the domain that does represent spatial variability. An example of such a system is a highly dynamic unsaturated porous medium zone, or a dynamic system such as a lake, coupled to a more slowly changing saturated porous medium zone. These sorts of extensions are also potentially fruitful ways to employ this foundational work in rigorous study of multiscale systems [e.g., those described in 2, 25, 27].

11 Conclusions

Reduced dimensionality models are important tools for both preliminary analysis and as component parts of larger, more complex models. Derivation of these megascale models poses challenges that are identical to those faced in deriving macroscale models. In both cases, a rigorous foundation requires an exact representation of the averaged variables of concern in terms of the microscale precursors, constraints to ensure that the second law of thermodynamics is enforced, and closure relations that represent sub-scale processes adequately at the larger scale of concern. Corresponding parallels exist between the derivation of TCAT models at the macroscale and at the megascale in concept, although the detailed components of the process must be extended; these extensions were detailed in this work.

Conservation of mass, momentum, and energy and balance of entropy equations were derived for a megascopic single-fluid-phase porous medium system starting from the microscale. These equations do not include spatial derivatives, but they do include a set of boundary integral terms that represent the net input of the conserved or balanced quantity of concern. These conservation and balance equations are defined as precise averages of microscale quantities, which allows for connection across scales.

Classical irreversible thermodynamics is used along with the megascale conservation and balance equations to produce an augmented entropy inequality. A series of manipulations are applied to produce a constrained entropy inequality (CEI) for the megascale single-fluid-phase system of concern, which is an exact expression for the system of concern and the thermodynamic theory chosen. The CEI is an important inequality that can serve as a starting point for a variety of megascale models for single-fluid-phase systems of varying complexity. Because the significant amount of effort required to derive this equation does not need to be repeated, varying applications of this inequality are considered relatively straightforward.

Because the CEI is not strictly in force-flux form, a series of approximations are applied to reduce this inequality to the simplified entropy inequality (SEI), which is necessarily approximate in nature unlike the CEI. The approximations used to derive the SEI are clearly summarized and the resultant SEI formulation is detailed. While perhaps not the final word, the megascale SEI derived is still a rather general expression that may apply to many different systems of varying complexity. A significant number of models can be derived based upon this SEI and only if these models prove inadequate to describe a physical system of concern would the SEI need to be revisited.

The SEI is used to guide the formulation of closure relations. A linear form of these closure relations is detailed, which are deemed a reasonable first approximation but certainly neither unique nor of high complexity. The closure relations are in turn used to produce a complete set of closed equations for conservation of mass, momentum, and energy. These closed equations can be applied to model a range of physical systems. Resolution of the relation between the precise functional dependence of the megascale model parameters appearing in the closed models in relation to specific underlying microscale, or macroscale, systems is an area of additional research deserving of attention using either analytical or numerical means.

In addition to purely megascopic models, the foundational work presented herein can be used to produce megascale models that are coupled to other megascale models, or megascale models that are coupled to macroscale models. The foundation laid in this work should enable such applications in a relatively straightforward manner.


This work was supported by National Science Foundation grant DMS-0327896 and by National Institute of Environmental Health Sciences grant P42 ES05948.


Roman letters
area, l2
entropy source density
[mathematical sans-serif bold C]
Green’s deformation tensor, ––
internal energy density
the set of entities
conservation of energy equation, ml2/t3
connected set of entities
scalar force for entropy generation
vector force for entropy generation
[mathematical sans-serif bold F]
tensor force for entropy generation
set of forces for entropy generation
general scalar function
general vector function
acceleration vector due to an external force, such as gravity, l/t2
scalar magnitude of gravitational acceleration, l/t2
heat source density
[mathematical sans-serif bold I]
identity tensor, – –
[mathematical sans-serif bold I]
surface identity tensor, ––
index set of entities
index set of connected entities
index set of interface entities
index set of phase entities
index set of boundary subsets
scalar flux for entropy generation
vector flux for entropy generation
[mathematical sans-serif bold J]
tensor flux for entropy generation
solid-phase Jacobian, ––
[K with circumflex]
megascale conductivity parameter
megascale kinetic energy per unit mass due to microscale velocity fluctuations, l2/t2
[K with circumflex]H
megascale hydraulic conductivity parameter, l/t
inter-entity heat transfer closure parameter
heat flux closure parameter
length, l
characteristic length, l
mass, m
conservation of mass equation, m/t
transfer of mass from the κ entity to the ι entity per unit volume per unit time, m/(tl3)
outward unit normal vector from the surface of a phase when ι [set membership] x2110P and the unit normal tangent to the surface and outward normal from the bounding curve when ι [set membership] x2110I
general property
conservation of momentum equation, ml/t2
fluid pressure, m/(lt2)
transfer of energy from the κ entity to the ι entity resulting from heat transfer and deviation from mean processes per unit volume per unit time, m/(lt3)
non-advective heat flux density vector
is the outward normal component of the heat flux vector for the ι entity evaluated on the boundary of the domain
[mathematical sans-serif bold R with ^]
resistance tensor
[r with circumflex]υ
closure vector
[r with circumflex]µ
closure vector
entropy balance equation, ml2/(t3T)
CIT-based thermodynamic equation for material derivative of internal energy
transfer of momentum from κ entity to the ι entity due to stress and deviation from mean processes per unit volume per unit time, m/(l2t2)
stress tensor
time, t
volume, l3
initial volume of the solid phase, l3
velocity, l/t
velocity of the exterior portion of the boundary of the ι entity, l/t
is the normal velocity of the ι entity relative to the normal velocity of the external boundary evaluated at the external boundary, l/t
mass averaged velocity of the ι entity relative to the mass averaged velocity of the solid phase, vwvs , l/t
weight function
material coordinate position vector, l
position vector in the solid phase, l
Greek letters
boundary of domain of interest
interfacial tension
measure of quantity of entity ι per macroscale volume
entropy density
entity qualifier
entropy production rate density
vector of Lagrange multipliers
Lagrange multiplier
chemical potential
mass density
Lagrangian stress tensor for the solid phase
stress tensor
transfer of entropy from the κ entity to the ι entity per unit volume per unit time
[var phi]
entropy density flux vector
acceleration potential (e.g., gravitational potential)
spatial domain
Subscripts and superscripts
energy equation qualifier
exterior of domain qualifier
general index
general index
general index
mass equation qualifier
momentum equation qualifier
residual portion of specified equation
index that indicates a solid phase
ith location where the s phase intersects the system boundary
thermodynamic equation qualifier
entity index corresponding to the wetting phase
ith location where the w phase intersects the system boundary
entity index corresponding to the wetting-solid interface
ith common line where the ws interface intersects the system boundary
entity qualifier
Other mathematical symbols
closure of set (overline)
left angle bracket right angle bracket
averaging operator
material derivative
[partial differential]′/[partial differential]t
partial derivative of a point on a potentially moving interface
microscale surficial del operator on an interface
spatial derivative of the thermodynamic function (·) taken with the quantities listed following the vertical bar held constant
averaged classical irreversible thermodynamics
augmented entropy inequality
constrained entropy inequality
classical irreversible thermodynamics
entropy inequality
representative elementary volume
simplified entropy inequality
thermodynamically constrained averaging theory


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


1. Bailyn M. A Survey of Thermodynamics. New York: American Institute of Physics Press; 1994.
2. Barenblatt GE, Zheltov IP, Kochina IN. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. Journal of Applied Mathematics (USSR) 1960;24(5):1286–1303.
3. Bear J. Hydraulics of Groundwater. New York: McGraw-Hill; 1979.
4. Callen HB. Thermodynamics and an Introduction to Thermostatistics. New York: Wiley; 1985.
5. Cermelli P, Fried E, Gurtin ME. Transport relations for surface integrals arising in the formulation of balance laws for evolving fluid interfaces. Journal of Fluid Mechanics. 2005;544:339–351.
6. Darcy H. Determination of the laws of flow of water through sand. In: Freeze RA, Back W, editors. Physical Hydrology. Stroudsburg, PA: Hutchinson Ross; 1983.
7. De Groot SR, Mazur P. Non-equilibrium Thermodynamics. New York: Dover; 1984.
8. Eringen AC. Mechanics of Continua. 2nd ed. Huntington, NY: Krieger; 1980. Notes.
9. Fosdick R, Tang H. Mathematics and Mechanics of Solids. 2008. Surface transport in continuum mechanics. page, OnlineFirst.
10. Gray WG. On the definition of derivatives of macroscale energy for the description of multiphase systems. Advances in Water Resources. 2002;25(8–12):1091–1104.
11. Gray WG, Ghidaoui MS. Thermodynamic analysis of stream flow hydrodynamics. Journal of Hydraulic Research. 2010 in print.
12. Gray WG, Hassanizadeh SM. Averaging theorems and averaged equations for transport of interface properties in multiphase systems. International Journal of Multiphase Flow. 1989;15(1):81–95.
13. Gray WG, Miller CT. Thermodynamically constrained averaging theory approach for modeling flow and transport phenomena in porous medium systems: 1. Motivation and overview. Advances in Water Resources. 2005;28(2):161–180.
14. Gray WG, Miller CT. Thermodynamically constrained averaging theory approach for modeling flow and transport phenomena in porous medium systems: 3. Single-fluid-phase flow. Advances in Water Resources. 2006;29(11):1745–1765. [PMC free article] [PubMed]
15. Gray WG, Miller CT. Consistent thermodynamic formulations for multiscale hydrologic systems: Fluid pressures. Water Resources Research. 2007;43(9)
16. Gray WG, Miller CT. Thermodynamically constrained averaging theory approach for modeling flow and transport phenomena in porous medium systems: 5. Single-fluid-phase transport. Advances in Water Resources. 2009 In press. [PMC free article] [PubMed]
17. Gray WG, Schrefler BA. Analysis of the solid phase stress tensor in multiphase porous media. International Journal for Numerical and Analytical Methods in Geomechanics. 2007;31:541–581.
18. Gray WG, Leijnse A, Kolar RL, Blain CA. Mathematical Tools for Changing Scales in the Analysis of Physical Systems. Boca Raton, FL: CRC Press, Inc.; 1993.
19. Jackson ABS, Miller CT, Gray WG. Thermodynamically constrained averaging theory approach for modeling flow and transport phenomena in porous medium systems: 6. Two-fluid-phase flow. Advances in Water Resources. 2009 In press. [PMC free article] [PubMed]
20. Jou D, Casas-Vazquez J. Extended irreversible thermodynamics and its relation with other continuum approaches. Journal of Non-Newtonian Fluid Mechanics. 2001;96(1–2):77–104.
21. Marle CM. On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media. International Journal of Engineering Science. 1982;20(5):643–662.
22. Maugin GA. The Thermomechanics of Nonlinear Irreversible Behaviors: An introduction. Singapore: World Scientific Press; 1999.
23. Miller CT, Gray WG. Thermodynamically constrained averaging theory approach for modeling flow and transport phenomena in porous medium systems: 2. Foundation. Advances in Water Resources. 2005;28(2):181–202. [PMC free article] [PubMed]
24. Miller CT, Gray WG. Thermodynamically constrained averaging theory approach for modeling flow and transport phenomena in porous medium systems: 4. Species transport fundamentals. Advances in Water Resources. 2008;31(3):577–597. [PMC free article] [PubMed]
25. Murad MA, Cushman JH. Multiscale flow and deformation in hydrophilic swelling porous media. International Journal of Engineering Science. 1996;34(3):313–338.
26. Pinder GF, Gray WG. Essentials of Multiphase Flow and Transport in Porous Media. Hoboken, NJ: John Wiley & Sons, Inc.; 2008.
27. Reggiani P, Rientjes THM. Flux parameterization in the representative elementary watershed approach: Application to a natural basin. Water Resources Research. 2005;41(W04013)
28. Whitaker S. Introduction to Fluid Mechanics. Englewood Cliffs, NJ: Prentice-Hall, Inc.; 1968.