Validity of Dmirey's Law Undf r

Transient Ooiidi

ur,s.
Validity of Darey's Law Under 

Transient Conditions

By CHARLES E. MONGAN

U.S. GEOLOGICAL SURVEY PROFESSIONAL PAPER 1331

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON : 1985
DEPARTMENT OF THE INTERIOR 

WILLIAM P. CLARK, Secretary

U.S. GEOLOGICAL SURVEY 

Dallas L. Peck, Director

For sale by the Distribution Branch, U.S. Geological Survey, 

604 South Pickett Street, Alexandria, VA 22304
CONTENTS

Page

Abstract ..................................................................... j

Introduction.................................................................. j

Application of Darcy's law...................................................... j

Description of analytical model.................................................. 2

Mathematical analysis ......................................................... 3

Interpretation of transient-flow equation.......................................... Q

Significance of first bracketed term ........................................... 7

Significance of second bracketed term ......................................... 9

Significance of third bracketed term .......................................... 9

Comparison with Darcy's flow equation .......................................... JQ

Summary and conclusions ...................................................... n

Selected references ............................................................ 12

Appendix: Supplemental derivation of equations ................................... 14

List of symbols and abbreviations ............................................... 15

ILLUSTRATIONS

Page

FIGURE 1. Conceptual model of a porous body ..................................... 3

2. Apparatus for producing selected hydraulic pressure-gradient functions ...... 4

3. Relation between time constant and radius of flow tube ................... g

4. Comparison between Darcy and transient-flow equations for indicated step increase

in pressure gradient .............................................. 10

TABLE

Page 

TABLE 1. Short table of finite Hankel transforms ................................. 5

ill

VALIDITY OF DARCY'S LAW 

UNDER TRANSIENT CONDITIONS

By CHARLES E. MONO AN

ABSTRACT

Darcy's law, which describes fluid flow through porous materials, 

was developed for steady-flow conditions. The validity of applying this 

law to transient flows has been mathematically verified for most 

ground-water flow conditions. The verification was accomplished 

through application of Hankel transforms to linearized Navier-Stokes 

equations which described flow in a small-diameter cylindrical tube 

chosen to represent a single pore in a porous medum.

INTRODUCTION

When fluid flows are studied, the assumption is frequently made that the flow is in the steady state. 

However, most flows in nature are in transient states. 

Flows in aquifers and oil-bearing strata are examples 

of such transient flows occurring in porous media. Flows 

in porous media are governed by Darcy's law under 

steady-state conditions. This same law is assumed to be 

valid with transient flows. The purpose of this report 

is to develop a transient-flow equation from fundamental physical laws and to compare the results from such 

an equation with the results obtained from an empirical 

Darcy's law. An equality of results will verify that Darcy's law may be applied with assurance to most 

transient-flow conditions.

To treat the problem mathematically, certain idealizations are necessary. These will be discussed in detail in 

the following sections.

APPLICATION OF DARCY'S LAW

Darcy's law (1856) states that the volumetric discharge 

of fluid through a porous body is proportional to the 

saturated cross-sectional area perpendicular to the direction of flow and is inversely proportional to the length 

of the flow path. It also is proportional to the difference 

in hydraulic head between the input and output 

surfaces.

Darcy's law may be written (Hubbert, 1969, p. 31):

Q=-KA (1)

where

Description Dimensions

Q=volumetric flow rate

A=cross-sectional area of porous

medium perpendicular to

direction of flow 

K= hydraulic conductivity of

porous medium 

h=hydraulic head 

L=length of flow path

LT-1

L

L

Darcy's formulation (eq. 1) followed from his numerous 

experiments. Many subsequent experiments, both in the 

laboratory and in the field, have shown its validity for 

steady-state conditions. Because of the law's importance, 

both for the management of ground-water resources and 

for its relevance to petroleum production, vigorous efforts have been made to give it a substantial theoretical 

basis (Gray and O'Neill, 1976).

Following these efforts, it is now customary to write 

the equation in a slightly different form; Hubbert (1969, 

p. 59) writes:

q=-a

dL

(2)

where, in addition to terms previously defined,

Description Dimensions

q=specific rate of fluid flow LT~l 

o=specific fluid conductivity LT1" 1 

<}»=fluid potential; $=gh L2 T~2 

g= acceleration (rate and direction) due to gravity LT~2

Scheidegger (1964) writes as follows:

q= grad p-Qg » (3)
VALIDITY OF DARCY'S LAW UNDER TRANSIENT CONDITIONS

where, in addition to the terms previously defined, 

Description Dimensions

H= dynamic viscosity of fluid

Q= density of fluid

k=intrinsic permeability of

porous medium (see Lohman,

1972) 

p=hydraulic pressure at a given

point

ML-3

L2 

ML-1 ?1-2

As Hubbert points out, the form of Darcy's law is 

similar to that of the other transport laws namely, 

Ohm's law, Fourier's law, and Picks' law, which pertain, 

respectively, to the flow of electrical current, the flow 

of heat, and the diffusion of chemical concentration. In 

all these cases, the laws apply to the steady state.

In the steady state, a flow field is characterized by the 

constancy of the velocity components at each given point 

in the field. In the unsteady state, the velocities at a 

given point change with time. A change with time implies a change of momentum, which, of course, involves 

the density of the fluid.

The absence in the Darcy equation of any explicit 

variable for time and the absence of density in a 

dynamical context indicate that the Darcy equation is 

an approximation for the unsteady state. Although it 

has been clear for some time that the law is an approximation, the real question continues to be: "How good 

an approximation is it?" The intent in this paper is to 

set the stage for further numerical estimates of the error made in using Darcy's law in its customary form 

with respect to current field applications.

The present problem is to obtain a better estimate of 

the error resulting from neglecting the inertial aspects 

of the transient-flow process. This paper seeks to furnish 

such an estimate on the basis of an analytical study of 

a simple flow-field element. From such an element, an 

approximation to a porous body may be constructed, at 

least conceptually.

A customary approach to a mathematical analysis 

such as that discussed in this report is to define a 

geometric model of the physical circumstances that are 

to be studied. Of course, it is necessary to examine the 

model in some detail to be sure that it properly reflects 

reality. The model should be simple enough so that it 

is amenable to direct mathematical description and 

calculation.

DESCRIPTION OF ANALYTICAL MODEL

There are certain properties that characterize porous 

bodies, one of which is the presence of continuous

passages through the material of the body that are 

capable of conducting fluid. The passages do not have 

constant cross-sectional areas they change both in area 

and direction. In addition, they may be characterized by 

a large surface area-to-volume ratio. The surface is the 

areal extent of the rock passage that is actually wetted 

by the liquid, and the volume is the volume of liquid contained within that surface.

Over the years, various attempts to design 

mathematical models capable of describing in detail the 

flow of fluid through porous earth materials have been 

made. Brief mention will be made of some of these to 

give the tenor of the analytical developments up to the 

present.

Happel and Brenner (1965, p. 389) make reference to 

an earlier survey of the literature on the efforts to 

describe Darcy's law from fundamental principles. Their 

review can be highlighted by the following three summary paragraphs:

1. Slichter (1899) used an arrangement of spheres to 

idealize a porous earth material and assumed a 

triangular cross section for a flow passage. It was 

assumed that Poiseuille's law could be applied to 

describe fluid flow in the triangular section. 

However, this model was inadequate for general 

use because it was oversimplified.

2. Blake (1922) introduced the concept of a hydraulic 

radius, which he defined as the volume filled with 

the liquid divided by the wetted surface. He 

regarded a flow passage in a porous body as a tube 

having a very complicated cross section. Using the 

hydraulic radius so defined, Blake was able to 

achieve results that were quite useful. It should be 

noted that Blake's hydraulic radius is the inverse 

of the surface area-to-volume ratio that plays a considerable role in catalytic chemistry particularly 

in chemical engineering.

3. Additional mathematical models were designed by 

Kozeny (1953), and a semiempirical model was 

made by Carman in 1956. The Kozeny-Carman 

work led to expressions for the constant in the 

Darcy equation which could be related to the void 

space and the hydraulic radius. Empirical 

developments of these constants are now commonly 

used in chemical engineering to analyze flow 

problems.

Happel and Brenner themselves made an analysis 

based on similitude and came up with constants for the 

flow equation that showed that the Darcy constant is 

dimensionally proportional to L squared divided by a 

constant C. This was a useful advance.

In a later study, Scheidegger (1964) made an analysis 

based on probability considerations. The pore structure
MATHEMATICAL ANALYSIS

was supposed to have a random distribution which could 

be given in probability terms. The motion of a fluid particle through such a flow field was assumed to be a random walk. Scheidegger made some use of the Poiseuille 

concepts.

Scheidegger makes the comments that the earlier lines 

of analytical development were twofold: (1) 

mathematical models based on probability, and (2) 

models based on the geometrical stacking of flowpassage components.

Many ground-water investigators since 1965 have examined a wide variety of flow problems and have 

published the results of their theoretical derivations. 

Because none of the foregoing investigative studies permit the Darcy constant to be calculated for an actual 

assemblage of porous earth particles, the fluid-flow 

analysis continues to rest on a somewhat empirical 

basis.

This paper combines some of the respective probabilistic and geometric analytical features. As can be 

seen in figure 1, it is possible to simulate a porous body 

by assembling a bundle of capillary tubes. If the crosssectional area of each tube is described by a Gaussian 

distribution, then we have achieved a random collection 

of the areas through which the liquid moves. A representative sample of such a collection can then be made in 

the three-dimensional form of a small elemental "pancake" (fig. 1).

Elemental 

pancake and 

representative^ 

collection of ^ 

tubes

Individual 

capillary 

tubes

Random crosssectional areas 

of tubes

Tubes oriented 

at random 

angles

Random spacing 

of elemental 

pancakes

FIGURE 1. Conceptual model of a porous body.

A second kind of pancake can be made by having the 

tubes of the same diameter but inclined or oriented in 

a random fashion, and a third arrangement can be made 

by having the separations between the pancakes reflect

a Gaussian distribution. These concepts presuppose complete vertical continuity of flow from one pancake to 

another.

Thus, starting with a single tube, it is possible to build 

up a structure by a sequence of probability distributions 

that leads to the simulation of a porous body. The final 

description of the entire flow field, which can be 

calculated from the Gaussian distributions, hinges on 

and is determined by the behavior of the fluid moving 

through a single capillary tube.

In his discussion of mathematical models, Garding 

(1977, p. 6) remarks that physicists tend to validate their 

models by measurements in the real world; mathematicians, on the other hand, use rules of logic.

With respect to physical measurements, W. O. Smith 

(oral commun., 1968) noted that he was satisfied with 

the laboratory measurements he made on capillary 

tubes using an apparatus substantially as shown in 

figure 2. Capillary tubes and undisturbed samples of 

porous bodies were the subjects of extensive laboratory 

experiments conducted by Smith. Unfortunately, his 

numerical data were not organized and published prior 

to his death and, therefore, only this author's 

mathematical analysis is described in this paper.

The analytical technique adopted by this author is one 

that features a study of the dynamic responses of fluid 

in a single capillary tube subjected to the specific stimuli 

of selected hydraulic pressure-gradient functions. Thus, 

a simple flow element is selected namely, a small slug 

or cylinder of water; it has already been shown that 

there are various ways of constructing the flow field for 

a porous body from such a simple element. The equations that are developed to describe the behavior of the 

simple flow element are then carried through for application to the whole porous body.

The foregoing procedure has a certain resemblance to 

differential calculus, whereby complicated volumetric 

structures are obtained by integration from elemental 

volumes. The procedure also is similar to finite-element 

programming, in which the properties of the element are 

known and are combined into larger structures by rules 

that define the interconnections.

MATHEMATICAL ANALYSIS

Consider the experimental flow system consisting of 

a horizontal capillary tube (fig. 2). The advantage of this 

orientation is that the elevation component of the 

hydraulic head is removed from the calculations. The 

outlet end of the tube vents to the atmosphere so that 

a constant pressure (zero reference) prevails at least for 

the duration of a given experiment. The inlet end of the 

tube is in a closed water-supply tank or reservoir provided with an automatic control system that allows the
VALIDITY OF DARCY'S LAW UNDER TRANSIENT CONDITIONS

internal pressure to be changed with respect to a time 

function as may be prescribed for the experiment.

3. instantaneous imposition of a ramp-function change 

(increase or decrease) in the pressure gradient.

Pressure sensor

Water-supply reservoir Open to 

atmosphere

Capillary tube

Pressure 

generator

Controller

Programer

Clock

FIGURE 2. Apparatus for producing selected hydraulic pressure-gradient functions.

The experimental setup, for which the mathematical 

analysis is developed, would consist of (1) initially 

establishing a suitable inlet pressure, (2) allowing the 

system to operate until steady-state equilibrium is 

achieved, and (3) imposing on the reservoir, at an arbitrarily chosen reference time t=Q, any one of the 

following three pressure-gradient conditions or any combination thereof:

1. instantaneous release of reservoir pressure to atmospheric pressure, which allows flow system to 

"decay" exponentially;

2. instantaneous imposition of a step-function change 

(increase or decrease) in the pressure gradient; and

To describe the operation of this experimental flow 

system mathematically, the Navier-Stokes equations are 

used; they are deemed valid because they are derived 

directly from Newton's basic laws.

It is well known that the Navier-Stokes equations are 

analytically intractable. It is necessary to make some 

adaptation to reduce them to a form that can be solved.

According to certain observations by Reynolds (1883), 

in the laminar-flow state, a particle of fluid in a tube 

moves in a straight line. This is a physical observation 

that can be expressed in mathematical terms. If this is 

done, the Navier-Stokes equation (Aris, 1962, p. 182) is 

simplified and can be written as
MATHEMATICAL ANALYSIS

dt Q dx dr

(4)

where, in addition to terms previously defined, 

Vx = fluid velocity, axial direction of tube, 

r = radial distance from tube axis or centerline,

cylindrical coordinates, 

t = elapsed time,

v = ^, kinematic viscosity of fluid,

x = distance along tube axis or centerline, and

p hydraulic pressure.

Inspection of the simplified Navier-Stokes equation 

shows that it is of the parabolic type (Hadamard, 1952, 

p. 38). The boundary conditions require that there be 

a statement of the initial condition and two statements 

regarding the fluid velocity at two different points in 

space. The initial condition used is simply the statement 

that the velocity of a fluid particle is a function of its 

radial distance, r, from the tube axis. This has the 

mathematical form

(5)

where

B0 - preexisting steady pressure gradient at time t=Q 

and

a = radius of flow tube.

The velocity at the wall is zero at all times. This meets 

the requirement for one spatial statement:

Vx (a, 0=0. (6)

The other spatial statement comes from an examination 

of the conditions on the tube axis, which reveals that 

the velocity at r=0 is a singular point (see Powers, 1979, 

p. 21). However, for physical reasons, we know that the 

velocity is finite and we can therefore write

Vx (0, t) (7)

where C is some constant. The desired spatial statement, 

however, is distilled from this by noting that on the tube 

axis, the velocity gradient dVy/dr is zero at all times, 

and we can write

dVx 

dr

(0, t)=Q. (8)

The problem is now well presented.

The equation is similar to the equations describing the 

flow of heat and also diffusion. There are several 

methods available for solving the simplified NavierStokes equation, including the transform method or the 

separation of variables; a numerical solution also may 

be developed. A solution for the case of a step-function 

drive and zero initial conditions was obtained by Holland 

and Brenner (1965) by the application of the separation 

of variables technique. There seems to be an advantage 

in using a transform technique in particular, a finite 

Hankel transformation. The conditions along the axis 

of symmetry of the tube and the no-slip condition at the 

wall are the same for many cases. For this analysis, a 

finite-transform technique was used; special use also was 

made of a Bessel function as a kernel. This technique 

has certain advantages when applied to this particular 

problem, because the initial and boundary conditions 

can be fitted very precisely to the Bessel functions.

The partial differential equation to be solved (see steps 

1-4 in "Appendix") is

(9)

dr3- rdr

jdVx uv x

dT [i

where 

T = t, a selected change in time scale (see step 3 in

"Appendix");

G! = a step-function pressure-gradient change; and 

Gr = a ramp-function pressure-gradient rate of change.

If we now operate on equation 9 and take the finite

Hankel transform, as expressed by Sneddon (1951) and

in table 1, we obtain

«

i+^rTrU dO)

TABLE I. Short table of finite Hankel transforms 

[Sneddon, 1951, p. 83, 531]

Real domain Transformed domain

a'-r2

d2 Vr 1 dV.X

dr2 r dr

aC K-.

-K?V

1C= constant.

where

V = fluid particle velocity in the transformed domain, 

KI=symbol denoting order of roots used in Hankel 

transform, and
VALIDITY OF DARCY'S LAW UNDER TRANSIENT CONDITIONS

«/! =Bessel function, first order, first kind, 

The operation of "taking the transform" causes one 

variable to be "integrated out" (Sneddon, 1972, p. 10). 

Because r=0 and r=a are the limits of integration, the 

spatial boundary conditions are automatically entered. 

Equation 10 reflects the results of these operations.

Equation 10 is an ordinary differential equation, to 

be integrated with respect to the variable, T, and with 

initial conditions available to determine the constant of 

integration. Recall that the initial conditions are a 

parabolic velocity distribution which is transformed by 

applying the results in table 1 to equation 5.

Integration equation 10 with respect to T gives

-bT (11)

where

b=K2i and

L

The initial conditions have been incorporated through 

evaluation of the constant of integration. This equation, 

in the transformed domain, shows the time evolution of 

the velocity for all cylindrical tubes. This has the advantage, therefore, of removing the space variable and 

describing the transient situation with respect to the one 

variable of time. This equation form can reveal the time 

evolution of the flow-velocity profile.

The inverse finite Hankel transformation (Sneddon, 

1951, p. 83) can now be performed on equation 11 to yield 

the following desired solution in the real domain:

(12)

where, in addition to terms previously defined, 

z = r/a, radial distance as a fraction of tube radius;

o TC -

j=- -2 time constant, ith term;

«/o =Bessel function, zero order, first kind; and 

aj =ith root, Bessel function; zero order, «70 (aj)=0.

INTERPRETATION OF TRANSIENT-FLOW 

EQUATION

Equation 12 has a structure that may be interpreted 

from mathematical and physical points of view. The 

most logical approach seems to warrant looking first at 

some of the mathematical features. As opportunities 

arise, however, physical interpretations are interwoven.

The dependent variable, V& is by definition the fluid 

velocity in a direction parallel to the axis of the tube 

at any given point a radial distance, r, from that axis. 

Vx is also a function of time, and thus its complete range 

of definition extends from r=0 to r=a in the space domain and from t=Q to t=°° in the time domain.

In examining equation 12, the velocity, Vx, is seen to 

comprise the sum of three bracketed terms which, in the 

order given, describe (1) a component related solely to 

the preexisting (prior to t=0) steady pressure gradient, 

B0, and to the exponential decay of the flow system after 

that gradient is dropped to zero; (2) a component related 

solely to the imposition of a pressure-gradient step function, GI, and (3) a component related solely to the imposition of a pressure-gradient ramp function, Gr. This 

suggests a very desirable versatility when it comes to 

analyzing practical problems in the field. Almost any 

conceivable pressure-gradient pattern should be 

amenable to analysis by subdividing it into the best 

possible arrangement of equivalent pieces that represent 

B0 , G1? and Gr pressure-gradient terms. Through the 

principle of superposition, a valid analysis can be made. 

Furthermore, this versatility of equation 12 is not contingent on the initial conditions being zero.

It is of interest to note that a "snapshot" of the pointvelocity profile across the tube radius during a transient 

condition is not a parabola. Nor is a series of such snapshots the case of a small parabola growing exponentially 

into a large parabola. Instead, there is throughout the 

evolutionary process a slight warpage away from the 

parabolic form. Although this deviation may be small, 

there is the view expressed by Birkhoff (1960) that, in 

hydrodynamics, there are paradoxes whereby (1) symmetrical causes do not always produce symmetrical effects and (2) small causes sometimes produce large 

effects. Thus, it seems appropriate to carry the fine structure of the point-velocity analysis to the useful limits 

afforded by equation 12 before combining the results 

with any further mathematical operations.

Equation 12 should first be tested to illustrate that 

the initial (at t=0) and boundary conditions have indeed 

been satisfied. If t=Q, each exponential term equals
INTERPRETATION OF TRANSIENT-FLOW EQUATION

unity and, therefore, a simple inspection of equation 12 

shows that the third bracketed term is zero. There now 

remains, in each of the first and second bracketed terms, 

a series summation, the limiting value of which (when 

£=0) was given by Bowman (1958, p. 17) as (I-s?)/4. This 

shows that for £=0, the second bracketed term also is 

zero and, therefore, equation 12 reduces to the simple 

form

Vx (13) (r, 0=-B0- 

M 4

By inspection, equation 13 satisfies the two stated 

spatial boundary conditions if it is noted that the expression dVy/dr is zero when 2=0 (flow-tube center) and that 

the velocity, Vx, is zero when 2=1 (flow-tube wall).

An assessment of the individual significance of the 

three bracketed terms (and their coefficients) in equation 12 is now presented.

SIGNIFICANCE OF 

FIRST BRACKETED TERM

Expansion of the indicated series summation may be 

written as

i t

aJ,(ai)

As the index i increases by unit steps, the successive 

roots a i become larger when evaluating the series summation for a given time, t. (See values given on fig. 3.) 

In other words, successive exponential terms become 

smaller, and if the time factor is now allowed to increase 

the combined effect is that "spectral terms" of higher 

order decrease faster than those of lower order. Thus, 

the "fine structure" in the velocity profile disappears 

early, leaving the first term in the series summation (eq. 

14) as the principal descriptor of the transient-flow 

condition.

An important aspect in the analysis of transient flow 

is the action of the inertial terms as exhibited by the 

time constants. The numerical value of the time constant for a simple tube representing a pore can be 

calculated. The value of the time constants changes enormously with the radius of the tube. The transient 

process that is, the evolution of the velocity patterncan be visualized as starting at the tube wall and proceeding inward to the axis. This mode of analysis provides additional insight into the transient process.

A partial insight into the relative significance of the 

first three terms in the series summation of equation 

14 is gained from the data plotted as figure 3, although 

they were developed primarily to show how the time constant, Tc^, increases with an increase in flow-tube 

radius. By definition, the parameter Tc^ is a function 

of Q, n, a, and aj in the manner repeated in figure 3. A 

dimensional analysis of Tc^ shows that it does indeed 

have the dimensions of time, and numerical values can 

be computed for any fluid if Q and ju have been determined. To prepare figure 3, values of e and ju were taken 

for distilled water at 68 °F (20 °C) and values of the time 

constant, Tc^ and then computed for selected tube 

radii first for the Bessel function root at from the first 

term of equation 14, next for the root a2 from the second 

term of equation 14, and finally for the root a3 from the 

third term of equation 14. The three curves plotted in 

figure 3 graphically illustrate the dominant importance 

of the first term in equation 14 an importance that is 

considerably strengthened by noting that the Bessel 

roots, aj, appear to the cubed power in the term 

denominators in that series-summation equation. Thus, 

the value of af is only about 14 in the first term 

denominator, whereas, it has grown to about 648, or 

nearly fifty fold, in the third term.

The physical significance of just the first bracketed 

term (and its coefficient) in equation 12 can be 

demonstrated by neglecting, for the moment, the second 

and third bracketed terms by assuming that both Gx and 

Gr are equal to zero. By definition, B0 is the preexisting 

steady-pressure gradient maintained on the flow system 

up to time £=0. In the experimental setup shown in 

figure 2, this presumes that steady pressure (and the 

necessary amount of fluid) ha& been maintained in the 

water-supply reservoir, which has a very large volume 

relative to the tube, at some convenient value above atmospheric pressure, for a long enough period so that the 

rate of flow through the capillary tube has stabilized and 

reached an unvarying or steady state. At time £=0, the 

pressure in the reservoir is instantaneously lowered to 

atmospheric. Nothing more is done no energy is externally added to or subtracted from the flow system, and 

the system is simply allowed to "run down." This is a 

classic exponential decay situation, and the exponential 

term in the first bracketed term of equation 12 precisely describes the manner in which the velocity, Vx, 

"bleeds off." As t-+-<*>, the value of the exponential term 

(and hence V^h^O. Prior to £=0, the fluid flowed through 

the tube at a steady rate, and a certain level of kinetic 

energy obviously was being maintained. With removal 

of the constant pressure gradient at £=0, the dissipative 

processes set in and the kinetic energy, acting through 

the viscous forces, converted to heat, which decreased 

fluid velocity.
VALIDITY OF DARCY'S LAW UNDER TRANSIENT CONDITIONS

15

14

13

12

11

10

O Q

111 y 

(/j

H ' 

(/J

O 

O 

LU 6

Roots, Bessel function

JQ (y) =o

a l =2.405 

a 2 =5.540 

a 3 = 8.654 

<* 4 =11.792 

a 5 = 14.931

Note: a// . ,

- r. = =0 -9962

where ^? = 998.2 kg/m 3 (ML" 3 ) *

pt = 1.002 x 10" 3 Pa-s (ML" 1 T~1 ) *

s = tube radius in millimeters 

(* Distilled water at 20° C.)

345678

TUBE RADIUS, a . IN MILLIMETERS 

FIGURE 3. Relation between time constant and radius of flow tube.

10
INTERPRETATION OF TRANSIENT-FLOW EQUATION

SIGNIFICANCE OF 

SECOND BRACKETED TERM

To examine the significance of this term, neglect the 

first and third bracketed terms of equation 12 by assuming that B0 and Gr are equal to zero. Physically, this 

means that, prior to time £=0, the flow system for the 

experimental setup of figure 2 is quiescent at atmospheric pressure that is, there is no pressure gradient on or flow through the capillary tube. At time £=0, 

a step-function pressure gradient, G1? is applied to the 

tube by instantaneously raising the pressure in the 

water-supply reservoir to some convenient value above 

atmospheric pressure and by holding that imposed 

pressure constant thereafter. Holding the pressure constant means that energy (and fluid) must be fed continuously into the flow system. The supplied energy is 

accounted for in two ways: (1) energy that must overcome inertia and accelerate the fluid is stored as momentum, and (2) the remaining energy that must overcome 

the viscous forces in the fluid is dissipated as heat. If 

the rate of energy dissipation has increased to a value 

that equals the constant rate at which energy is being 

fed into the flow system, no further acceleration of the 

fluid can occur and the system has reached a steady 

state.

The second bracketed term of equation 12 contains the 

identical series summation that has already been shown, 

in its expanded form, as equation 14. Discussion of that 

equation showed that, as time elapses, the value of the 

summation decreases and very quickly approaches zero. 

As this happens, the only term left in the brackets is 

(l-2?)/4, which contains no time factor. Thus there is 

no further change with time that is, the steady, or unvarying, state has been reached.

SIGNIFICANCE OF 

THIRD BRACKETED TERM

To examine the unique significance of this term, 

neglect the first and second bracketed terms of equation 

12 by assuming that B0 and Gx are equal to zero. As in 

the preceding discussion, this means that, prior to time 

£=0, the flow system or the experimental setup of figure 

2 is quiescent at atmospheric pressure. At time £=0, a 

ramp-function pressure gradient, Gr, is applied to the 

tube by raising the pressure in the water-supply reservoir at a steady rate. The nature of the flow system 

response can be addressed by noting that the series summation is now a bit different from that exhibited in both 

the first and second bracketed terms of equation 12. 

Although the expansion of the series summation has the 

same form as that shown as equation 14, the multiplier 

of each expanded term is Tc^(l-e~ ci). As time increases, this term in the series quickly approaches zero

and, thus, the summation expansion reverts to the 

simpler form:

2 T

2 Tc J0

+

a\

|

2 7k J0 (a32)

2 Ten J0 (anz)

(15)

The root characteristics of the Bessel function J0 and J: 

resemble a cosine function and a sine function, respectively. Their combined effect, as successive terms in the 

series summation (equation 15) are evaluated, is to approach a finite limit. This occurs after only the first two 

or three terms are used; however, each still needs to be 

multiplied by the appropriate time-constant parameter, 

Tc^ Furthermore, a time parameter, t, also is the 

multiplier for the first term appearing in the third 

bracketed term of equation 12.

One way of visualizing what is described by the third 

bracketed term of equation 12 as time elapses is to say 

that the response, Vx, although trying to catch up with 

the stimulus, Gr, succeeds only in approaching a finite 

lag for any given radial distance, r, from the flow-tube 

axis. This finite lag is different for each radial distance. 

Thus, the velocity profile "warps" with the time. 

Because the flow system has inertia, the response can 

never completely catch up with the stimulus, which is 

increasing continuously at a steady or linear rate.

Interpretation of the transient-flow equation is not 

complete without elaborating on the nature of the time 

constant, Tc^ This is a key parameter in characterizing the transient-flow process.

The time constant, as calculated from its definition in 

equation 12, reflects the nature of the assumptions made 

in devising a single capillary-tube model. Thus, it would 

not be expected to yield numerical values that necessarily agreed with similar constants determined from other 

investigators' works, inasmuch as the latter are based 

on different sets of modeling assumptions.

The time constant may be regarded as the time interval that extends from the moment a linear time invariant system is perturbed by, say, a step function until the system has reached 63 percent of the new steadystate condition. A time invariant system is one whose 

properties, such as dimensions, viscosity, density, and 

temperature, are constants. A linear system is one to 

which the mathematical principle of superposition applies. Time constants appear in the exponential terms 

in the transient-flow equation; their significance may
10 VALIDITY OF DARCY'S LAW UNDER TRANSIENT CONDITIONS

be illustrated by three distinct fluid-flow situations:

1. time constant is very short compared with the duration of the flow process,

2. time constant is of the same order as the duration 

of the process, and

3. time constant is long compared with the process 

duration.

In the first case, which is common in geology, the transient flow effects can be neglected; in the second case, 

those flow effects should be estimated in order to make

a good approximation; and in the third case, transient 

flow effects must be taken into account to provide a 

realistic view of the whole flow process.

COMPARISON WITH DARCY'S FLOW EQUATION

A method of comparing the analytical results obtained 

from the Darcy (steady-state) and transient-flow equations is computing the volume rate of discharge for a 

given flow condition with each of the two equations (fig. 

4). Consider a representative sample of a porous body

SE

A. Imposed pressure-gradient step

B. Discharge computed from Darcy's Law through 

a representative porous sample

(D 

cc.

o

CO C. Discharge computed from transient-flow 

equation through one capillary tube

Q

D. Superposed discharge curves for Darcy and 

and transient flows

2 3 

MULTIPLES OF TIME-CONSTANT UNIT,

1

FIGURE 4. Comparison between Darcy and transient-flow equations for indicated step increase in pressure gradient.
SUMMARY AND CONCLUSIONS 11

initially quiescent and then subjected to a step-function 

pressure gradient, Glf (fig. 4A). From Darcy's law, an 

expression for the volume rate of discharge, Qd, can be 

written directly from equation 1 in the form

(16)

where

Equation 16 contains no time factor and obviously 

describes an unvarying or steady-state discharge condition (fig. 4B).

Now consider a capillary tube of radius a, in a quiescent state, and then apply the same step-function 

pressure gradient (Gx). No other change is imposed, and 

the response is uniquely described by simply the second 

bracketed term (and its coefficient) of equation 12, which 

defines, as a function of time, a velocity profile at right 

angles to the tube axis. This term can be integrated as 

a function of radial distance from the tube centerline, 

between the limits r=0 and r=a. This yields the following expression for the volume rate of discharge, Qt, as 

a function of time:

no*G1

-t/T

£ (17)

Because the series summation in the referenced term 

of equation 12 converges rapidly, only the first term of 

the expansion needed to be integrated. Also, because of 

a-!4 equals 33.2, the coefficient of the exponential term 

in equation 17 can be replaced with the value 1. Equation 17 may, therefore, be rewritten in the form

Qt = l-e (18)

Equation 18 may, therefore, be regarded as a 

mathematical model for the transient-flow state showing the discharge response from a capillary tube subjected to a step-function pressure gradient.

If the value of t becomes sufficiently large compared 

with TCi , the exponential term in the brackets (equation 18) becomes negligible. Indeed, if t=5TCi, the value 

of the exponential term is less than 0.02. Thus, for any 

larger values of t, the bracketed part of equation 18 may 

be considered equal to unity, and, thus, Qt has reached 

an unvarying or steady state. The foregoing observations 

are readily apparent in figure 4C, which is a plot of equation 18.

To compare the difference between the two preceding 

methods for computing rate of discharge under the 

stipulated flow conditions, the discharge curve shown

in figure 4C is superposed on the curve shown in figure 

4B. Note, however, that the curve in figure 4C 

represents the rate of discharge through just one 

capillary tube. Because this single tube is representative 

of all tubes, it is permissible to multiply the rate of 

discharge in the one by some number, n (number of 

tubes), such that the combined or total rate of discharge 

exactly equals the steady rate of discharge (graphed in 

fig. 4B) for the representative sample of the porous body 

when steady state is reached.

The superposed results are as shown in figure 4Z>; the 

hachured area indicates the difference between the two 

forms of analysts during the transient or unsteady part 

of the flow period. Observe that the two curves are indeed coincident when the elapsed time, t, has increased 

to 5TC . Thereafter, the two forms of analysis yield identical results because steady-state conditions prevail in 

the flow system.

For an aquifer of medium gravel, in which the average 

pore-space radius might be 1 mm, the corresponding 

value of Tc , as read from figure 3, is about 0.5 second. 

This suggests that, in response to some step-function 

change imposed on the aquifer, the duration of the transient flow period would be about 2.5 seconds.

The validity of the preceding discussion, which compares the results obtained from the Darcy equation 

rather than the transient-flow equation, is restricted in 

no way by the arbitrary choice of a step-function 

pressure gradient as the stimulus to the flow system. 

The choice was primarily for purposes of illustration.

SUMMARY AND CONCLUSIONS

The mathematical analysis in this paper is based on 

equation 4. If we multiply each term in that equation 

by the fluid density, Q, the left-hand side represents a 

change in the momentum or a force. The first term on 

the right-hand side represents the causative force provided by the pressure gradient. The last term on the 

right-hand side can then be regarded as the dissipative 

or opposing force due to viscosity in effect, a frictional 

force. Hence, this version of equation 4 is a "balance of 

forces" equation.

Consider each term of the foregoing version multiplied 

by a velocity, Vx. The left-hand side of the resulting 

equation may be interpreted to mean the change in 

kinetic energy. The first term on the right-hand side 

represents energy input per unit time to maintain the 

pressure gradient. The second term on the right-hand 

side gives the rate at which frictional energy is 

dissipated as heat. We now have a representation of the 

fluid-flow process in terms of energy. In the steady state, 

the kinetic-energy term is constant that is, energy is 

being introduced into the system at a rate just sufficient
12 VALIDITY OF DARCY'S LAW UNDER TRANSIENT CONDITIONS

to compensate for the frictional losses. The advantage 

in considering the fluid-flow process from this point of 

view a balance in the rates of energy flow is that it 

gives a further insight into the transient as well as the 

steady-state flow conditions.

In a fluid-flow system, the essence of a transient process is its time evolution. Therefore, the time variable 

is normally the dominant parameter in any mathematical description of that system. In this report, the 

measure of time that characterizes the continuously 

changing fluid-flow process is identified as the time constant, and it first appears in equation 12. If the description of the transient process is given in terms of the time 

constant, the time evolution is independent of real time. 

Therefore, it is a simple matter to use equation 12 to 

analyze specific flow-field problems.

A closer look at the structure of the time constant is 

appropriate. The time constant obviously includes two 

kinds of information; the first relates to the properties 

of the fluid (/^ and Q), and the second relates to the 

geometry of the flow tube (radius r, and Bessel root, aj, 

which implies cylindrical shape). The dynamic viscosity, n, of the fluid is associated with the dissipative action in the flow system, and the fluid density, Q, is 

associated with the storage of kinetic energy. Thus, for 

any given flow system, the time constant must reflect 

the ratio between the storage of kinetic energy and the 

dissipative action, inasmuch as it includes the ratio of 

Q to p.

On the basis of mathematical analysis in this paper, 

equation 12 is shown to be a convenient tool for examining specific ground-water situations, regardless of 

whether they are simulated in the laboratory or encountered in the field. Equation 12 includes pressuregradient terms that are measures of energy input. In 

the laboratory, such inputs can be precisely controlled 

as to amplitude and time of application. Then, the 

response is predictable and reflects the design of the apparatus. Similar controlled inputs can be made in the 

field environment. Thus, an aquifer-recharge experiment of a selected magnitude can be started at any given 

moment and with any convenient well. The response will 

be predictable and will reflect the known physical properties and dimensions of the porous-earth material comprising the aquifier. Depending on the nature of the 

experiment laboratory or field the numerical value of 

the time constant can be determined (core sample 

analysis) and a judgment made as to the merits of using the transient-flow analysis.

The time constant increases as the square of the radius 

of the capillary flow-tube, which in turn is intended to 

represent a flow path through some porous earth 

material. If a tube radius of 3 mm, for example, is 

selected (that radius would be exceeded in the field only

by the pore radius of very coarse gravels or cavernous 

limestone), the time constant would be about 1.7 seconds 

(fig. 3). Thus, if a step-function type of pressure gradient 

were applied to this size of flow tube, the steady-state 

condition would be reached in about 8.5 seconds. This 

helps to place a practical upper limit on the value of the 

time constant. In other words, if a detailed description 

of the rate of discharge is needed within the first 8 

seconds of the foregoing flow conditions, equation 12 

should be used.

If a linear type of pressure-gradient change had been 

postulated as the basis for analysis instead of a step function, similar arguments could have been made with only 

some differences in detail, the principal difference being the development of a standing time lag in the 

response when compared with results obtained by applying Darcy's law. The time lag would be of the order 

of two time constants; again, equation 12 would permit 

the calculation for any specific set of conditions.

It may be concluded that, for most of the commonly 

encountered ground-water flow systems, the simple 

Darcy equation remains an adequate mathematical 

description. The interesting point is that any reasonable 

set of flow-field circumstances can now be evaluated 

rapidly and the most appropriate analytical procedure 

followed. Of course, experimental verification would be 

highly desirable.

One final observation seems pertinent to the direction 

that might be taken for any further research. In the 

mathematical analysis given herein, the velocity profile was described by Fourier series. This is a powerful 

analytical technique, and its application to further 

studies of the "structure" of fluid flows would enhance 

our knowledge and understanding of the intricacies of 

unsteady fluid flow processes.

The analysis presented in this report suggests directions for future research. One direction would be to examine the behavior of fluid flow through assemblages 

of tubes and connections organized from simple elements 

by probability rules. Another direction would be to 

measure, in the laboratory, the transient response of 

fluid flow through a core sample of porous earth 

material.

SELECTED REFERENCES

Aris, Rutherford, 1962, Vectors, tensors and the basic equations of fluid 

mechanics: Englewood Cliffs, N.J., Prentice-Hall, 286 p.

Birkhoff, Garrett, 1960, Hydrodynamics: A study in logic, fact and 

similitude: Princeton, N.J., Princeton University Press, 184 p.

Blake, F. C., 1922, The resistance of packing fluid flow: American Institute of Chemical Engineers Transactions, v. 14, 415 p.

Bowman, Frank, 1958, Introduction to Bessel functions: New York, 

Dover Publications, 135 p.
SELECTED REFERENCES 13

Carman, P. C., 1956, Flow of gases through porous media: New York:

Academic Press. 

Darcy, Henri, 1856, Les fontaines publiques de la ville de Dijon: Paris,

Victor Dalmont. 

Garding, Lars, 1977, Encounter with mathematics: New York,

Springer, 270 p. 

Gray, W. G., and O'Neill, Kevin, 1976, On the general equations for

the flow in porous media and their reduction to Darcy's law: Water

Resources Research, v. 12, no. 2, p. 148-154. 

Hadamard, Jaques, 1952, Lectures on Cauchy's problem in linear partial differential equations: New York, Dover Publications, 316 p. 

Happel, John, and Brenner, Howard, 1965, Low Reynolds number

hydrodynamics with special applications to particulate media:

Englewood Cliffs, N.J., Prentice-Hall, 553 p. 

Holland John, and Brenner, Howard, 1965, Low Reynolds number

hydrodynamics with special applications to particulate media:

Englewood Cliffs, N.J., Prentice-Hall, 553 p. 

Hubbert, M. K., 1969, The theory of ground-water motion and related

papers: New York, Hafner Publishing Co., 310 p. 

Kozeny, J., 1953, Hydraulik: Vienna, Springer, 390 p.

Lohman, S. W., 1972, Ground-water hydraulics: U.S. Geological Survey 

Professional Paper 708, 70 p.

Powers, David, 1979, Boundary value problems (2d ed.): New York, 

Academic Press, 351 p.

Reynolds, Osborne, 1883, An experimental investigation of the circumstances which determine whether the motion of water shall 

be direct or sinuous and of the law of resistance in parallel channels: London, Royal Society Philosophical Transactions, v. 174, 

935 p.

Scheidegger, A. E., 1964, Statistical hydrodynamics in porous media, 

in Chow, V. T., Advances in hydroscience: New York, Academic 

Press, v. 1, p. 161-181.

Slichter, C. S., 1899, Theoretical investigation of the motion of ground 

waters: U.S. Geological Survey 19th annual report, pt. 2, 

p. 301-384.

Sneddon, I. N., 1951, Fourier transforms: New York, McGraw-Hill, 

542 p.

 1972, The use of integral transforms: New York, McGrawHill, 539 p.
14 VALIDITY OF DARCY'S LAW UNDER TRANSIENT CONDITIONS

APPENDIX: SUPPLEMENTAL DERIVATION 

OF EQUATIONS

The following numbered steps show the manner in 

which the transient-flow equation was derived, starting 

from the modified Navier-Stokes equation.

1. Navier-Stokes equation simplified:

dVx== _l 3p+ v\d 2 Vx ,13^3 

3* e 3* L ar2 r dr

2. Introduction of the time-dependent gradient, G^ Grt, 

a simple polynomial forcing function that reflects 

a physically attainable state:

3. Insert change of time scale:

T = vt and dT = vdt.

4. Equation to be solved:

5. Boundary conditions:

a. no slip at the flow-tube wall:

Vx (a, t] = 0 

b. finite velocity on the flow-tube axis:

Vx [0, t] < C 

c. symmetry:

-^[0, t] = 0

d. initial conditions:

6. Equation transformed by finite Hankel transform:

7. Integration of the transformed equation:

The primitive is:

\

w

= \--G1 C1 eKZiTdT

r*

C0

vVe

^cl^-JJ/F+c .... 1\ 1^2 IT'4 /e ' ^0

8. Evaluate the constant of integration, using the 

transformed initial conditions.

atT=0, V = --J, 

M Ki

c° = 7T K*- JI (aKd + j

9. Solution in the transformed domain:

+

bv

10. Perform the inverse transformation. Use table of 

transform pairs and also definition of inverse 

transformation.

- ~ T , i aJ0(ar)

Tc - fth term Tc{ =
SYMBOLS AND ABBREVIATIONS 15

LIST OF SYMBOLS AND ABBREVIATIONS

Symbol Dimensions Description

A L2 Cross-sectional area of porous medium perpendicular to

direction of flow.

B0 ML~2T~2 Preexisting steady pressure gradient at time t=Q.

C0 T~l Constant of integration.

r C T a u i * <fcWff> t L Symbol for ^ 

D ML-'T*

G! ML~2 T~2 A step-function pressure-gradient change.

Gr MLr*T~3 A ramp-function pressure-gradient rate of change.

e/0 Bessel function, zero order, first kind.

e/i Bessel function, first order, first kind.

K LT~* Hydraulic conductivity of porous medium.

Ki ................ Symbol denoting roots used in the Hankel transform;

J0(aKi)=Q.

L L Length of flow path. 

Q L3 ?1-1 Volumetric rate of flow.

Qd L3 ?1-1 Volumetric rate of flow calculated according to Darcy's

law.

Q t L3T~l Volumetric rate of flow calculated according to equation

12.

T U Changed time scale; T=vt.

DO2

Tc . L2 Time constant, ith term; T= - M°i

Vx LT~l Fluid velocity in axial direction of tube.

V T~l Fluid particle velocity in transformed domain.

a L Radius of flow tube.

b .......... .... .. Symbol for K*.

c ................ Constant.

g LT~Z Acceleration due to gravity.

h L Hydraulic head.

i ................ Index in integers.

k L2 Intrinsic permeability of porous medium.

p ML~1 T~2 Hydraulic pressure at a given point.

q LT~l Specific rate of fluid flow.
16 VALIDITY OF DARCY'S LAW UNDER TRANSIENT CONDITIONS

Symbol Dimensions Description

r L Radial distance from tube axis.

t T Elapsed time.

x L Distance along tube axis or centerline.

z ................ Radial distance as a fraction of tube radius; z=r/a.

M ML~1 T~1 Dynamic viscosity of fluid.

v UT~^ Kinematic viscosity of fluid.

6 ML~3 Density of fluid.

a T Specific conductivity of fluid.

<f> L2T~2 Fluid potential; $=gh.