Field Determination of 

Vertical Permeability to 

Air in the Unsaturated Zone

GEOLOGICAL SURVEY PROFESSIONAL PAPER 1051
Field Determination of 

Vertical Permeability to 

Air in the Unsaturated Zone

By EDWIN P. WEEKS

GEOLOGICAL SURVEY PROFESSIONAL PAPER 1051

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

CECIL D. ANDRUS, Secretary

GEOLOGICAL SURVEY 

H. William Menard, Director

Library of Congress Cataloging in Publication Data

Weeks, Edwin P., 1936-

Field determination of vertical permeability to air in the unsaturated zone.

(Geological Survey Professional Paper 1051)

Bibliography: p. 41

1. Sediments (Geology) Permeability. 2. Sediment (Geology) United States.

I. Title. II. Series: United States Geological Survey Professional Paper 1051.

QE471.2.W43 624'.1513 77-10571

For sale by the Superintendent of Documents, U.S. Government Printing Office

Washington, D.C. 20402 

Stock Number 024-001-03092-6
CONTENTS

Page

Abstract.............................................................................................. 1

Introduction........................................................................................ 1

Previous work.............................................................................. 1

Terminology................................................................................ 2

Notation...................................................................................... 2

Conversion of English to metric units of measurement............... 3

Theory................................................................................................ 3

Sources of error in converting air permeability to hydraulic

conductivity.............................................................................. 4

Effects of moisture content.......................................................... 5

Klinkenberg effect........................................................................ 5

Absorption effects........................................................................ 5

The method......................................................................................... 6

Site inspection ............................................................................. 6

Instrumentation........................................................................... 6

Data obtained.............................................................................. 7

Data analysis............................................................................... 9

Results................................................................................................ 11

Birch Creek Playa, Idaho............................................................ 11

Lubbock, Texas, airport site........................................................ 13

Comparison of results for different sites and dates.............. 15

Page

Results Continued

Lubbock, Texas, airport site Continued

Comparison of results with recharge data............................ 16

Correction for structural alteration of materials by wetting. 18

Comparison of results with laboratory data......................... 19

Discussion of results ............................................................ 20

Long Island, New York .............................................................. 20

Summary and conclusions.................................................................. 22

AIRK-A Fortran IV program for computing hydraulic conductivity

from pneumatic head readings................................................. 24

Computation of pneumatic head values ...................................... 26

Numerical solution of equation 3 ................................................ 26

The automatic search procedure.................................................. 26

Program verification.................................................................... 27

Input data requirements .............................................................. 28

List of mnemonics and constants used in the AIRK program .... 30

AIRK program listing ................................................................. 31

TRIMAT subroutine listing......................................................... 37

Sample computer output............................................................. 38

Sample output data..................................................................... 38

References cited.................................................................................. 41

ILLUSTRATIONS

Page 

FIGURE 1. Piezometer nest used to determine pneumatic head differences in the unsaturated zone...................................................................... 7

2. Photograph showing piezometers, manifold, and manometer used to determine pneumatic head...................................................... 7

3. Diagram showing manifold construction.............................................................................................................................................. 8

4. Comparison of pneumatic heads measured at four depths at Lubbock, Tex., airport site A to best-fitting analog simulation............ 9

5. Sample logs and construction diagrams of wells tested at Idaho National Engineering Laboratory, Idaho......................................... 11

6. Comparison of pneumatic heads measured at the Birch Creek Playa, Idaho, site to the best-fitting numerical simulation................. 13

7. Locations of instrumentation used in air-permeability determinations and recharge experiments near the Lubbock, Tex., 

Regional Airport............................................................................................................................................................................... 14

8. Sample log of test well at Lubbock, Tex., airport spreading site ..............................................................................................:.......... 14

9. Comparison of pneumatic heads measured at five depths at Lubbock, Tex., airport site B to best-fitting numerical simulation ....... 16

10. Cross-section showing unsaturated materials at Nassau County, Long Island, N. Y., Recharge Basin 306......................................... 20

11. Comparison of pneumatic heads measured at two depths at an installation inside Recharge Basin 306, Nassau County, Long Island,

N. Y., to the best-fitting numerical simulation................................................................................................................................. 22

12. Sketch showing the relationship between an air-permeability installation and the nodal configurations used to solve for air permeability.. 25

in
IV CONTENTS

TABLES

Page

TABLE 1. Summary of vertical air permeability values computed from pneumatic head data obtained from installations at the Birch Creek

playa, Idaho......................................,............................................._ 12

2. Depths of piezometers in three installations used to determine vertical air permeability at the Lubbock, Tex., airport spreading site... 14

3. Summary of hydraulic conductivity values computed using pneumatic head data obtained from installations at the Lubbock, Tex., 

airport spreading site, 1969-72.......................................................................................................................................................... 15

4. Comparison of hydraulic conductivity values estimated from recharge data, from the air-permeability determinations, and from 

laboratory analyses of drive cores and wireline cores....................................................................................................................... 17

5. Laboratory data for 30 cores from the Ogallala Formation in the Southern High Plains of Texas, showing comparisons of 

measured hydraulic conductivity to those computed from air permeability ..................................................................................... 19

6. Summary of analyses of pneumatic head data obtained from installations at the Nassau County, Long Island Recharge Basin 306. 21
FIELD DETERMINATION OF VERTICAL 

PERMEABILITY TO AIR IN THE UN SATURATED ZONE

By EDWIN P. WEEKS

ABSTRACT

The vertical permeability to air of layered materials in the unsaturated 

zone may be determined from air pressure data obtained at depth during a 

period when air pressure is changing at land surface. Such data may be 

obtained by monitoring barometric pressure with a microbarograph or 

surveying altimeter and simultaneously measuring down-hole pneumatic 

head differences in specially constructed piezometers. These data, coupled 

with air-filled porosity data from other sources, may be compared with the 

results of electric-analog or numerical solution of the one-dimensional 

diffusion equation to make a trial-and-error determination of the air 

permeability for each layer. The permeabilities to air may in turn be 

converted to equivalent hydraulic conductivity values if the materials are 

well drained, are permeable enough that the Klinkenberg effect is small, and 

are structurally unaffected by wetting.

The method was tried at several field sites. Results from three of these 

sites are compared with laboratory data and air-injection or infiltration 

data. Permeabilities to air determined for playa sediments from a test in 

Idaho are about two times higher than laboratory data, and about 50 

percent lower than the permeability to air determined from air injection 

tests. Equivalent hydraulic conductivities determined from permeability 

values to air for the bottom four of six layers of unsaturated materials at a 

site in the semiarid Southern High Plains of Texas are in good agreement 

with hydraulic conductivity values determined by analysis of infiltration 

data, indicating that the method is useful for determining the hydraulic 

conductivity of layers at depth under certain conditions. However, 

hydraulic conductivity computed from the permeability to air for the 

surface layer at this site is a few times higher than the infiltration-determined value. The same general results were obtained for the surface layer 

from several other tests, suggesting that the method may not be suitable for 

determining hydraulic conductivity near the surface. Comparison of results 

for determinations of permeability to air with infiltration data for two 

installations near a spreading basin on Long Island, N.Y., were contradictory. Hence, applicability of the method to determine the hydraulic 

conductivity of unsaturated materials in which residual moisture content is 

relatively high needs further testing.

Permeabilities to air determined by this method are representative of 

relatively large volumes of materials, and the hazard of disturbing the 

materials during measurement is small. Moreover, the method measures 

vertical permeability, rather than horizontal permeability or some composite of horizontal and vertical permeability. Hence, the method offers 

potential advantages over present methods to evaluate sites for artificial 

recharge by spreading; to evaluate ground-water pollution hazards from 

feedlots, sanitary landfills, and land irrigated with sewage effluent; and to 

evaluate sites for temporary storage of gas in the unsaturated zone.

INTRODUCTION

This report describes application of a method proposed by 

Stallman (1967) to determine, in situ, the vertical permeability to air of materials in the unsaturated zone. The method is 

based upon the observation that, when the atmospheric 

pressure changes at land surface, air moves to or from the 

unsaturated zone to maintain a pressure balance between air 

in the soil and the atmosphere. The rate of this air movement 

and the resultant rate of pressure change at depth are 

affected by both the permeability and the air-filled porosity 

of the materials in the unsaturated zone. Consequently, the 

media properties may be determined by analysis of the time 

lag of pressure changes at depth relative to those at the 

surface.

For this study, the method proposed by Stallman was used 

to determine the vertical permeability to air of each layer of a 

sequence of layered unsaturated materials at several sites, 

based on data from relatively simple and easily instrumented 

installations. Also, applicability of the air-permeability data 

to compute hydraulic conductivity of layers at depth was 

tested at two sites by comparison with hydraulic conductivities estimated from the behavior of artificially created 

perched ground-water mounds.

Stallman's method differs from most previous field 

methods for determining either the permeability to air or the 

hydraulic conductivity of unsaturated material in that it 

provides measurements in the vertical direction both at the 

surface and at depth. Such knowledge of the vertical 

permeability of materials in the unsaturated zone is potentially useful for planning water-spreading operations; for 

assessing the ground-water pollution hazard of sanitary 

landfills, feedlots, or sites of spray irrigation of sewage 

effluent; and for determining the feasibility of temporarily 

storing gas in the unsaturated zone.

PREVIOUS WORK

Movement of air to and from the unsaturated zone in 

response to barometric pressure changes was first described
VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

and analyzed by Buckingham (1904). In addition to describing the phenomenon, Buckingham presented an equation for 

the attenuation of the amplitude and the phase lag of a 

periodic atmospheric pressure wave at any depth in a 

homogeneous layer bounded below by an impermeable 

boundary. However, use of observed air pressure or 

pneumatic head at depth to determine the permeability of 

materials in the unsaturated zone was first proposed by 

R.W. Stallman (written commun., 1962). Stallman (1967) 

and Stallman and Weeks (1969) later described application 

of the method at Badger Wash, near Cuba, N. Mex. Data 

obtained at this site were analyzed based on the assumption 

that the unsaturated materials comprised a single 

homogeneous layer bounded below by an impermeable 

boundary. Because of instrument problems during data 

acquisition, the results of analysis of data from this site were 

inconclusive.

The method apparently was discovered independently by 

G.A. Morris and D.F. Snoeberger (1971). They (in 

collaboration with R.B. Rozsa, J. Baker, and C.J. Morris) 

have applied the method to determine the permeability to air 

of the material comprising several nuclear chimneys at the 

Nevada Test Site (Rozsa, Snoeberger, and Baker, 1975, and 

references contained therein). An elaborate data acquisition 

system was used (Snoeberger, Morris, and Baker, 1972), 

consisting of pressure transducers connected to pipes open at 

different depths in the chimney, and equipped to transmit 

data at 15-minute intervals by telephone to the Lawrence 

Livermore Laboratory at Livermore, Calif. For data 

analysis, the nuclear chimney rubble was assumed to consist 

of a homogeneous unit extending to infinity below land 

surface, and pressure changes at depth were analyzed using 

an analytical equation and the principal of superposition 

(Rozsa, Snoeberger, and Baker, 1975, p. 7). Although the 

assumption that air movement can occur to infinite depth in 

the chimney rubble poorly represents actual boundary 

conditions, measured and computed pressures at depth 

matched closely. Also, computed permeabilities to air 

compared favorably with those determined by numerical 

analysis of air injection data. The authors concluded that 

their application of the method gave results correct within an 

order of magnitude.

Other field methods for determining permeability to air 

have relied upon injecting air at land surface into the soil at 

the center of a circular surface area sealed by paraffin (Evans 

and Kirkham, 1949), into a tube inserted into the soil 

(Kirkham, 1946; Grover, 1955; Tanner and Wengel, 1957), 

into an auger hole (King, 1968), or into a packed-off section 

of drill hole (Boardman and Skrove, 1966). In general, these 

methods measure a permeability that is some composite of 

the vertical and horizontal permeability, and the 

measurements are representative of only a relatively small 

sample of material. Thus, these methods are not strictly 

comparable to the method proposed by Stallman.

TERMINOLOGY

Terminology involving discussion of the flow of both air 

and water through porous media is made more complicated 

by differences in conventional usage that has developed with 

regard to the different fluids. The permeability of a porous 

medium to air generally is described in terms of intrinsic 

permeability, and the term "air permeability" refers to the 

intrinsic permeability of the air-dry medium as measured 

with air or another gas. On the other hand, the permeability 

of a porous medium to water generally is reported as the 

"hydraulic conductivity," a term that is dependent on the 

kinematic viscosity of water at the prevailing temperature as 

well as the properties of the medium. Thus, the dissimilar 

terms air permeability and hydraulic conductivity are used in 

this report to reflect conventional usage.

Other complications arise from use of the term "intrinsic 

permeability," which is often defined as being a property of 

the medium alone (Hubbert, 1940, p. 818-819; Lohman and 

others, 1972, p. 9). As pointed out by Childs (1969, p. 

166-167), this concept may be misleading when applied to 

natural earth materials because of physical changes in the 

medium that occur upon exposure to different fluids, 

including waters of different chemical quality. Hence, 

intrinsic permeability is independent of the viscosity of the 

fluid used, but often is dependent upon the nature of the 

fluid. Because of the need to specify the fluid used to 

determine the intrinsic permeability, the term "permeability" 

is sometimes used in this report to designate intrinsic 

permeability to the specified fluid.

NOTATION 

A - area of perched mound at the top of a zone

being analyzed, m2 ; 

b - initial saturated thickness of water-table aquifer,

m; 

C - capacitance of a single capacitor used in the

resistor-capacitor grid for simulating gas flow

in the unsaturated zone, farads; 

CMC - quantity of expanding lattice clay minerals of the

montmorillonite (smectite) group, expressed as

a decimal fraction of the total mineral content; 

g = acceleration due to gravity, approximately equal

to 9.81 m/sec2 ;

: pneumatic head at start of test, m; 

: pneumatic head at time / during test, m; 

: pneumatic head at end of test, m; 

: hydraulic gradient, m/m; 

: intrinsic permeability,K 2 ; 

: relative permeability to air, defined as the ratio of

the air permeability at a given moisture content

to that for the air-dry medium, dimensionless; 

K - hydraulic conductivity, m/day; 

Kc = hydraulic conductivity based on field-determined

air permeability and corrected for effects of

wetting upon the structure of the medium,

m/day;

hi

ht

hf

I

k

K,ra
THEORY

Ke - hydraulic conductivity computed from

field-determined air permeability, m/day; 

L = thickness of unsaturated zone above capillary

fringe, m;

n - porosity, dimensionless; 

n^ - air-filled porosity at prevailing moisture content,

dimensionless;

p_ = pressure at a point, (N/m2 ); 

P - mean atmospheric pressure during a test period,

N/m2 ; 

PD = pneumatic diffusivity for an air-dry medium,

m2 /day; 

PDe = effective pneumatic diffusivity, determined at the

prevailing moisture content of the medium,

m2 /day;

Q = average inflow rate to spreading basin, m3 /day; 

R = resistance of a single resistor used in the

resistor-capacitor grid for simulating gas flow

in the unsaturated zone, ohms; 

RK - ratio of hydraulic conductivity measured with

water to that computed from air permeability

measured on the sample at its prevailing

moisture content, dimensionless; 

5 = drawdown in water-table aquifer, m; 

t - time in air-permeability test, seconds; 

z - vertical dimension in unsaturated zone, or height

above an arbitrary datum, m; 

H« = absolute viscosity of air at the prevailing soil

temperature, kg/m-sec; 

•\) w - kinematic viscosity of water at a given

temperature, m2 /sec; 

P a - density of air, kg/m3 ; 

(j) = pneumatic head or potential, height of

constant-density air column above reference

plane at the prevailing atmospheric pressure,

m.

CONVERSION OF ENGLISH TO METRIC UNITS OF 

MEASUREMENT

In this report, measurements are generally given in metric 

units, although those measurements relating to equipment

English units Conversion factor Metric units

Length in inches (in)

in feet (ft)

Area in square

feet (ft2)

inches (in2 )

Area in acres

Volume in cubic feet (ft3 )

Weight in pounds (Ib)

Pressure in

inches of mercury

x

x

x

x

x

x

x

x

25.40

.3048 =

.0929 =

645.2

.4047 =

.02832 =

.4536 =

.03386 =

millimeters (mm)

meters (m)

square meters (m2)

square millimeters (mm

hectares (ha)

cubic meters (m3 )

kilograms (kg)

bars (105 Newtons/

2 )

in atmospheres x 

Hydraulic conductivity

in feet per day (ft/d) x 

Permeability

in darcies x

square meter) 

1.013 = bars

.3048 = meters/day 

.987 = square microns

are given in both metric and English units. However, the 

AIRK program was written using English units, and English 

units are used in the description of that program to avoid 

confusion.

THEORY

The properties of a dry porous medium that affect the rate 

of pressure change at depth may be quantified by the term 

(Katz and others, 1959, p. 408):

PD-kP

where k_ - intrinsic permeability, L2 ;

P - mean pressure during pressure change,

M/LT; 

\ia = absolute viscosity of air at the prevailing

temperature, M/LT; 

n = porosity, dimensionless.

The quantity PD is here termed "pneumatic diffusivity" in 

analogy with the comparable term from ground-water 

hydraulics, "hydraulic diffusivity" (Lohman and others, 

1972, p. 8). However, the unsaturated zone generally 

contains water absorbed to the grains and held by capillarity 

at grain contacts, which reduces the permeability of the 

medium to air and the volume of pore space available for air 

storage. Thus, under field conditions, the medium property 

that affects air pressure changes may be termed "effective 

pneumatic diffusivity," or

KmkP

(1)

where Kra is the ratio of the air permeability at the 

prevailing moisture content to that for the dry medium; and 

n^ is interconnected air-filled porosity at the prevailing 

moisture content.

If it is assumed that gas flow due to a change in 

atmospheric pressure occurs only in the vertical direction, 

that absolute pressure is small enough that the ideal gas laws 

apply, that change in pressure with depth has a negligible 

effect on gas density, and that air permeability of the 

medium is large enough that the Klinkenberg1 (1941) effect 

may be ignored, the equation governing the isothermal flow 

of air in the unsaturated zone may be written (Katz and 

others, 1959, p. 408, eq. 10-17):

(2)

'The Klinkenberg effect occurs during gas flow through a capillary or porous medium when the 

mean free path of the gas molecules approaches in size the diameter of the capillary tube or pore. 

Under these conditions, the Hagan-Poiseuille velocity distribution no longer holds, as some gas 

molecules tend to slip along the capillary wall. Thus, for very fine grained materials, the intrinsic 

permeability as measured by gas flow exceeds that measured by liquid flow.
VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

where (]> = pneumatic head, equal to 

p - pressure at a point, N/ m2 ; 

Pa - density of air, kg/ m3 ; 

g = acceleration due to gravity, 9.81 m/sec2; 

and z = the height above a reference plane, m (positive

upward).

This partial differential equation is linear with respect to 

(|> 2 rather than ()>. However, for problems in which the 

pneumatic head varies only slightly from its mean value, the 

equation may be written:

dt

(3)

Strictly speaking, use of this equation requires that head 

terms in the solution be replaced by their squares. For 

example, if initial head for a problem is designated hi and 

final head hft and head at an intermediate time h^ the

dimensionless head, l~ J } would be replaced by the term 

2 hi-hf

J (Katz and others, 1959, p. 413-414). However, for the

i~ f

applications described below, the total pneumatic head 

change was never more than .01 P. Under these conditions,

the value

_ . , 

would not vary from that for_{ _ Lby more

than .5 percent. Thus, the use of equation 3 simplifies the 

computation without significantly changing the solution. 

The approach is analogous to that of using 

Dupuit-Forchheimer assumptions to analyze problems involving drawdowns in water-table aquifers. Thus, Jacob's 

(1963, p. 248) correction term, s-s2 /2b, where s=drawdown 

and 6=aquifer thickness, could be used to compute a 

corrected (hf-hj) term if (/z/-/z/) were used for s, and P/ Pag for 

b. Based on this analogy, such corrections, or the use of the 

/zHerms, would be necessary only if (hf-hj) were greater than 

AP.

Assuming that z=0 at the land surface and -L at the top of 

the capillary fringe above the water table, the flow system in 

the field is subject to the boundary conditions:

at z=0;

each contact between layers. For example, at the junction of 

layers n and n- 1, the boundary conditions:

6(t> TV =0 at z--L. Q L

The first of these boundary-condition equations specifies 

that the total pneumatic potential at land surface varies as an 

arbitrary function of time. The second equation specifies 

that no gas flow occurs across the bdundary z~-L 

(presumably the top of the capillary fringe).

The field situation is complicated by the fact that the 

materials in the unsaturated zone generally are layered, with 

the materials in each layer having a different intrinsic 

permeability and air-filled porosity. Under these circumstances, additional boundary conditions are imposed at

and

are imposed. These two equations specify that the flux from 

one layer equals that into the adjoining layer, and that 

pneumatic head is the same in the two layers at their 

common boundary.

Analytical solutions for flow in layered systems are 

complicated and difficult to evaluate (Carslaw and Jaeger, 

1959, p. 326). However, electric-analog or numerical simulation provides a feasible approach to the problem, because 

the internal boundary conditions will be met implicitly in the 

simulation model when the air permeabilities and air-filled 

porosities are correctly simulated.

To analyze data obtained from a layered sequence of 

materials, it is a good strategy to add one layer at a time to 

the simulation model, starting with the bottom layer. By 

following this procedure, and assuming that air-filled 

porosity is known for each layer from other data, the 

permeability to air at the prevailing moisture content of the 

uppermost layer will be the only unknown to be solved for 

during any given simulation. Successive estimates of this 

uppermost layer permeability can be substituted into the 

model until a satisfactory match between simulated and 

observed heads at the bottom of or within the layer is 

obtained, based on measured head variations at the top of 

the layer. Once this permeability has been determined, it may 

be assigned to that layer, and a new layer, with as yet 

unknown permeability, added to the model. The process is 

repeated until all the layers are included in the model, as 

described. in detail in the section on the AIRK (air 

permeability) program.

SOURCES OF ERROR IN CONVERTING

AIR PERMEABILITY TO HYDRAULIC

CONDUCTIVITY

The vertical hydraulic conductivity of materials in the 

unsaturated zone, rather than their air permeability, is 

generally the parameter of interest to hydrologists, and may 

be estimated from the air-permeability data based on the 

well-known viscosity conversion (Muskat, 1946, p. 72). 

However, care must be exercised in using these estimates, as 

they may differ from the true hydraulic conductivity because 

of three factors discussed in detail below.

1. The air-permeability measurements are made at the 

prevailing moisture content, and thus the computed 

product of relative air permeability times intrinsic
SOURCES OF ERROR

permeability (Krak) is lower than the true intrinsic 

permeability.

2. For tight materials, gas permeabilities are greater than 

liquid permeabilities because of the Klinkenberg (1941) 

effect.

3. Water will often react with the medium by absorbing 

either directly to the clay particles or to the cations 

absorbed on the clay surfaces. Such reactions may result 

in structural changes in the medium that drastically alter 

its intrinsic permeability (Johnston and Beeson, 1945).

EFFECTS OF MOISTURE CONTENT

The effect of prevailing soil-moisture content on air 

permeability does not appear to be significant in the semiarid 

regions where the method has been tested. A reasonable 

assumption under these conditions is that moisture content 

is at or below field capacity throughout the unsaturated zone 

down to the capillary fringe. Field capacity is generally 

described as the moisture content at which gravity drainage 

is incipient, and may be related to investigations of 

simultaneous flow of water and air by assuming that field 

capacity is the highest moisture content at which the effective 

hydraulic conductivity is zero. For this moisture content, the 

permeability to air has been found by several investigators 

(Botset, 1940; Osoba and others, 1951; Leverett and Lewis, 

1941) to range from 0.6 to 0.8 that of the air-dry medium. 

Thus, errors in estimating intrinsic permeability from air 

permeability at the prevailing moisture content would be 

relatively small under these conditions.

Errors resulting from the effects of moisture content on 

the permeability of the medium to air may be much greater 

where residual moisture content is greater than field 

capacity, as would commonly be the case in humid regions. 

Consequently, care should be taken to obtain air-pressure 

readings only after long dry periods, during which the 

sediments have had a chance to drain to near field moisture 

capacity.

THE KLINKENBERG EFFECT

This effect, as described by Klinkenberg (1941), occurs at 

low gas pressures in fine-grained materials, and results in gas 

permeabilities measured at low ambient pressures being 

larger than the corresponding liquid permeabilities. Katz 

and others (1959, p. 46) presented a chart for converting gas 

permeability to liquid permeability when the Klinkenberg 

effect is significant. However, because the Klinkenberg effect 

is significant only for materials having an intrinsic 

permeability of less than about 10 millidarcies, it generally 

could be ignored for the conditions encountered during this 

study.

ABSORPTION EFFECTS

Upon wetting, the granular materials comprising the 

unsaturated zone will absorb water, causing, in some cases,

structural and intrinsic permeability changes in the porous 

medium. Such changes generally result from the presence of 

clay minerals in the medium, and their magnitude is 

dependent on the initial moisture content of the mate rial, on 

the mineralogy and absorbed cation chemistry of the 

contained clays, and on the chemical quality of the 

percolating water. The water-quality effects occur because of 

absorption of water in greater or lesser amounts by the 

exchangeable cations associated with the clay minerals. For 

example, exchangeable sodium on clays in the presence of 

water of very low dissolved-solid concentrations will absorb 

relatively large amounts of water. Such absorption will cause 

swelling of montmorillonite and mixed-layer clays, resulting 

in great reductions in permeability (Quirk and Schofield, 

1955; McNeal and Coleman, 1966; McNeal and others, 

1968). Even when only non-swelling clay minerals are 

present, however, the absorption of water by the exchangeable cations results in disaggregation of clay particles 

and the migration into and blockage of pores by clay 

platelets (Johnston and Beescn, 1945, p. 45). This pore 

blockage may result in greatly reduced permeability of the 

medium.

Most comparisons of permeability to air and to water 

have been based on laboratory studies, and indicate that the 

differences in permeability to the two fluids can be extreme. 

For example, Johnston and Beeson (1945, p. 47) report 

reductions in permeability of as much as a few thousandfold 

upon changing the permeant fluid from air to brine, and of as 

much as several thousandfold upon changing from air to 

distilled water. Reeve (1953) also reports differences of as 

much as several thousandfold between permeability to air 

and to distilled water. These differences are undoubtedly 

much greater than those that would occur if the air 

permeability were measured at field moisture capacity rather 

than air dryness. Nonetheless, results of laboratory tests on 

30 cores from the Ogallala Formation in the Southern High 

Plains in Texas, described under "Results," show variations 

of as much as 500 to 1 between hydraulic conductivity 

measured directly and that computed from air permeability 

at the prevailing moisture content.

Differences between the hydraulic conductivities computed from air permeability test data and the hydraulic 

conductivities determined from recharge experiments are 

much smaller than those determined in the laboratory, as 

described under "Results." Sufficient data have been 

collected to indicate that the air permeability test gives a 

value for intrinsic permeability of the surface layer that is 

generally a few times higher than that determined from 

recharge experiments, probably as a result of structural 

changes in the medium due to wetting. On the other hand, 

determinations by the two procedures of intrinsic 

permeability of layers at depth were approximately equal for 

one extensive test, also described under "Results," despite 

the presence of swelling clays in some of the layers. These
VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

results suggest that the permeability of materials at depth 

may not be greatly affected by structural changes due to 

wetting, and that the air-permeability method might be 

useful for estimating their hydraulic conductivity. However, 

any investigator planning to use or test the method should be 

aware of the potential error in hydraulic conductivity 

estimates that might arise because of the effects of wetting on 

the structure of the medium.

THE METHOD

The first requirement for application of the described 

method to determine vertical air permeability at a selected 

site is an investigation to ensure that geologic and hydrologic 

conditions at the site meet, at least approximately, the 

theoretical assumptions on which the method is based. Once 

the investigator is satisfied that the site conditions are 

suitable, the site is instrumented. Wide latitude exists in 

selecting suitable instrumentation and installation 

procedures. The selections described here are relatively 

simple and inexpensive, and have proven successful in 

several field applications. Data collection can generally be 

completed during an afternoon barometric decline if normal 

weather prevails. Data analysis should be relatively straight 

forward for anyone with access to a large digital computer, 

using the FORTRAN IV computer program listed and 

described in the section on the AIRK program.

SITE INSPECTION

Geologic and hydrologic conditions will always be more 

complex than those assumed for the theory of this method. 

Hence, it is important that the conditions be investigated to 

ensure that they are satisfactorily approximated by the 

theoretical assumptions. To apply the method, it is assumed 

that the unsaturated zone consists of a sequence of horizontal layers, each being of large areal extent and uniform 

thickness. Each layer is also assumed to have areally and 

vertically uniform permeability and air-filled porosity. It is 

further assumed that both land surface and the water table 

are horizontal, and that the barometric changes are areally 

uniform. Finally, it is assumed either that no open boreholes, 

cracks, crevices, animal burrows, or other openings exist 

that would provide preferred flow paths for air movement, 

or else that the openings are so widely and uniformly 

distributed that the unsaturated zone still behaves as a 

porous medium. The validity of these assumptions as related 

to actual conditions at any field site should be carefully 

considered before applying the method at the site, and in 

selecting screen depths when the site is instrumented.

Subsurface geologic data, including sample and 

geophysical logs, are needed from at least one test hole in the 

site area. These data are used to determine the depth to 

water, to select the layers into which the unsaturated zone is 

to be divided, and to estimate the air-filled porosity for each 

chosen layer. If data are available for additional test holes, 

the continuity and dip of the beds can also be determined.

Selection of layers within the unsaturated zone will be 

clearcut if the unsaturated zone is strongly stratified into 

only a few units. However, if lithologic changes in the 

unsaturated zone are gradational, or the zone comprises 

thinly interbedded materials, the selection criteria must be 

more subjective, as the number of layers must be kept within 

practical limits. For this study, six was the maximum 

number of layers instrumented, but more could be instrumented if called for by the site geology.

As part of the site investigation, it should be ascertained 

that all wells and test holes in the vicinity of the site are 

grouted or backfilled through the unsaturated zone to ensure 

against air circulation in open well bores.

INSTRUMENTATION

Instrumentation for this method must be emplaced in a 

hole augered or drilled to the desired depth in the unsaturated zone. Drilling procedure is important, as the hole 

must be drilled without adding water or mud that could plug 

the wall of the hole and stop or restrict air movement into or 

out of the adjacent unsaturated material. In addition, 

drilling should be halted above the water table, as saturated 

auger or drill cuttings from below the water table tend to 

plaster the well bore above the water table with 

low-permeability material. Generally, augering is the best 

method for installing the hole, although for deep or 

consolidated sequences of unsaturated materials, the hole 

could be installed by air-rotary methods. All the holes used 

in this study were augered except those at the NRTS 

(National Reactor Testing Station, now Idaho National 

Engineering Laboratory), which were drilled by air-rotary 

methods.

Considerable latitude exists for variation in the construction of the air piezometers, which are installed above the 

capillary fringe in the lowermost layer and at each layer 

boundary, as selected during the site inspection. The screens 

should be relatively short to avoid air circulation within the 

screened section, and the pipes leading to the surface should 

be of small diameter to minimize storage effects within the 

pipe. During this study, the sites usually were instrumented 

using standard 32 mm X610 mm (1 I/4-inch by 24-inch) well 

screens connected to the surface by 6.4 mm (1/4-inch) steel 

pipe. At some sites, however, where the chosen layers were 

thin, 150 mm (6-inch) sections of 32 mm (1 1/4-inch) well 

screen or 25 mm (1-inch) diameter, 100 mm (4-inch) long 

brass foot valves were used to screen the piezometers.

The screens are emplaced in the auger hole to form a 

piezometer nest, as shown schematically in figure 1. The 

bottom screen is emplaced above and as near the top of the 

capillary fringe as feasible. Gravel or coarse sand is poured 

into the annulus around the screen and pipe to about 0.1 

meter above the top of the screen. A layer of dry dirt 

approximately 0.1-0.2 meter thick is dropped on the gravel 

to seal it from the grouting material. Expanding cement 

(Peck, 1964) is poured through a drop pipe or tremie (to
THE METHOD

avoid plugging the side of the hole with cement) to the level 

at which the next screen is to be installed. The process is 

repeated until the uppermost screen is installed, and the hole 

is cemented to land surface, thus completing the piezometer 

nest.

Piezometer 4

3

Land surface

Layer 4

Layer 3

Layer 2

Layer 1

Capillary fringe

r

\

§

-

>

•

i

-

•TlHM

ffiffl

FIGURE 1, Sketch showing typical piezometer nest used to determine 

pneumatic heads at selected depths in the unsaturated zone.

Pressure readings, in terms of liquid head, are obtained by 

connecting the piezometers through a manifold to an 

inclined manometer, as shown photographically in figure 2 

and schematically in figure 3. Differences in pneumatic head 

between each piezometer and the atmosphere are made by 

opening the manifold valve connecting that piezometer to 

the manometer, all the other manifold valves being closed. 

Zero manometer readings are obtained by opening the 

manifold to atmosphere. The change in atmospheric 

pressure during the test is measured with a microbarograph 

and (or) a surveying altimeter.

Under field conditions, a manometer inclined at a 10:1 

slope and filled with fluid of specific gravity 1.20 (commercial manometer fluid for outdoor use) has been found to give 

readings accurate to about 0.25 mm (0.01 inch), which is 

equivalent to about 30 mm (0.1 foot) of pneumatic head. A 

760 mm (30-inch) manometer, having a range of 76 mm (3 

inches) is adequate to measure pneumatic head differences 

under most conditions. Because downhole pneumatic head 

will sometimes be less than that at land surface, it is desirable 

to establish, by slightly overfilling the manometer-fluid 

reservoir, a reading of about 13 mm (0.5 inch) at zero

pressure difference. Thus, for a 760 mm (30-inch) 

manometer, a net range of about 64 mm (2.5 inches) is 

available for measuring the major head differences occurring 

during the test. Such a range generally is adequate even when 

unusually large atmospheric pressure changes occur.

FIGURE 2. Photograph showing piezometers, manifold, and manometer 

used to determine pneumatic head versus depth and time. Photograph 

by Richmond F. Brown.

DATA OBTAINED

Data required for the method include curves or tables of 

pneumatic head versus time at each depth during a period 

when atmospheric pressure is changing. During the investigation, it was found to be most convenient to obtain 

these data during a normal afternoon barometric decline. In 

most of the areas where tests were run, such a pressure 

decline usually is equivalent to 3 to 4 millibars or about 30 to 

36 m of pneumatic head in 4 to 5 hours. Generally, when 

determining the downhole response to diurnal barometric 

changes, differences in pneumatic head between each depth 

and land surface are measured in terms of liquid head at 

about 15-minute intervals. Readings generally are started at 

about 10 or 11 a.m., when the barometer is normally quite 

stable, as such readings provide good initial conditions for 

the analysis. The readings generally are continued through 

the afternoon decline until 6 or 7 p.m., when a diurnal 

barometric rise normally occurs.

For electric-analog analysis, the measured pneumatic 

head differences are added algebraically to the absolute 

pneumatic head, as determined from a microbarograph or 

an altimeter. The resultant pneumatic head values can be 

plotted versus time to produce a family of curves (fig. 4), one 

curve for each piezometer and one for atmospheric pressure,
VERTICAL PERMEABIL11 Y TO AIR IN THE UNSATURATED ZONE

Elbow

Board mount for manometer

Manometer

reservoir

Meniscus 

Hose connecting manometer reservoir and tube

Hose clamp

Vented plug 

(Atmospheric 

pressure)

Petcock valve

Pipe-to-hose 

connection

— Atmosphere 

vent

Globe valve

NOTE! All pipe filtings, valves, 

and hoses are standard 

1/4-inch items

Land surface

Cement grout

Piezometer 

number (pipe)

FIGURE 3. Schematic diagram of a manifold, showing its construction and details of its connection to the piezometers and manifold. The 

connection shown is that used when the downhole pressure is greater than that at land surface, and the valve settings are those for reading 

piezometer 1.

that may be used for the electric-analog analysis. Alternatively, for the digital computer analysis, the 

microbarograph data at selected times are punched on one 

set of cards. The differential manometer readings and their 

elapsed times are punched on a second set of cards, and both 

sets are read as data in the program. These data are 

converted by interpolation and algebraic manipulation 

within the program into a table of pneumatic head values, 

one for each screen and for land surface, for each short 

uniform time step used in the analysis. Details concerning 

data input to the program are described in the section on the 

AIRK program.

Pneumatic head differences occurring during the normal 

afternoon diurnal barometric change were found to be large 

enough to be read and analyzed at sites where the unsaturated zone was more than about 20 meters thick and 

there was at least one layer with a permeability of no more 

than 2 or 3 darcies. At other sites, where the unsaturated 

zone was thinner and the layers more permeable, head 

differences during such a diurnal change were too small to 

measure accurately using the described instrumentation. 

Moreover, short-term atmospheric pressure changes that 

could be characterized as noise tended to mask the 

differences. This noise would interfere with the pressure
THE METHOD

26.475

ANALOG PIEZOMETER 

SIMULATION NO. DATA

0.40

0.30

0.20

0.125

— 887

—885

1600 1800 2000 

MARCH 8,1969

2200 2400 0200 0300 

MARCH 9, 1969

FIGURE 4. Pneumatic heads measured March 8-9,1969, at four depths at Lubbock, Tex., airport site A, and best-fitting analog simulation.

difference readings even if more sensitive instrumentation 

were used. However, usable readings can sometimes be 

obtained for a relatively thin, highly permeable section 

during the passage of a "front" in which the atmospheric 

pressure changes a few millibars in a few minutes to an hour. 

Such events are so rare and unpredictable, however, that the 

author was only able to make use of them twice during 

several attempts. These experiences in "chasing fronts" 

suggest that general use of such events is impractical.

DATA ANALYSIS

Data obtained during this investigation were initially 

analyzed by use of an electric analog model, based on the 

analogy between the flow of a compressible fluid and the 

flow of electricity (Karplus, 1958, p. 356-360). Required 

equipment for the analog analysis includes a device for 

electrically exciting or stressing the system, analogous to the 

atmospheric pressure change at land surface; a system of 

variable resistor-capacitor (R-C) networks or grids to 

simulate the layered unsaturated materials; and a device for 

recording voltage changes at selected points within the R-C

grid due to the electrical excitation. The recorded voltage 

changes are analogous to the pneumatic head changes 

measured at depth. For this study, an electronic line follower 

fitted to an x-y recorder (Rogers and Connolly, 1960, p. 

37-38) was used to vary voltage, through high-gain 

direct-current amplifiers, in the R-C network. Voltage at 

selected junctions within the R-C network was read out 

through a large-value resistor and high-gain amplifiers to a 

second x-y recorder. The time scale for both x-y recorders 

was controlled by an external time-base device.

The electric-analog analyses, as performed, were made 

using a series of R-C grids, each grid consisting of 15 variable 

resistors individually connected to ground through fixed 

capacitors. The capacitors in each grid were of the same size.

Details of the electric-analog analysis are not given, 

because the system has been replaced by a digital computer 

program. However, the basic analytical procedure is similar 

to that described for the section on the AIRK program. As 

for the computer program, the layered system is analyzed by 

first treating the bottom layer as a complete system, and 

when that layer is adequately simulated, adding an R-C grid
10 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

to the system to represent the overlying layer. Air 

permeability of the most recently added layer is determined 

from the analog analysis by trial-and-error manipulation of 

the resistors in that grid until the model results, in terms of 

the voltage variation with time within or at the bottom of 

that grid, match the measured pneumatic head at the 

appropriate screen. Once a match is obtained, air 

permeability for that layer may be computed from the 

relationship among air permeability, air-filled porosity, 

thickness of unsaturated material represented by each 

resistor (equal to I/ 15th the layer thickness), the resistance 

and capacitance of each resistor and capacitor in the grid, 

and the ratio of analog to prototype time.

Because the electric analog system used for the initial 

analyses would not be available to most investigators, a 

digital computer program was formulated to analyze the 

pneumatic head data numerically while this report was in 

review. Comparison of the numerical and analog results 

indicated that the electric-analog analyses were incorrect. In 

retrospect, it was recognized that the capacitors in each grid 

representing a given layer should have been varied proportionally to the air-filled porosity of that layer, rather than 

being fixed at the same value for all the grids. Actually, the 

effect of variation in air-filled porosity in layers above the 

bottom one was rather minor compared to effects of 

variations in air permeability, and the analog system as used 

would have provided reasonably good results if air-filled 

porosity had been assumed constant for all the layers in the 

prototype system. Instead, the analyses were made by 

computing an effective pneumatic diffusivity for each layer 

based on the resistances and capacitances for the grid 

representing the layer. An air permeability was computed by 

multiplying this diffusivity by the appropriate air-filled 

porosity term for the layer. This procedure had no effect on 

determination of air permeability of the bottom layer, but 

for overlying layers with larger or smaller air-filled 

porosities, the computed air permeability was almost 

proportionately too large or too small.

Because of these problems with the electric-analog 

analysis, all of the described analyses were redone using the 

digital-computer program described and listed in the section 

on the AIRK program. The program, which is written in 

FORTRAN IV, uses data provided on the site, including 

screen depths and depth to the water table, estimates of 

air-filled porosity, preliminary estimates of PZ^-based 

hydraulic conductivity, and periodic microbarograph and 

differential manometer readings to compute final estimates 

of air-permeability based hydraulic conductivity.

The basic procedure used in the AIRK program to arrive 

at these final hydraulic conductivity estimates is that of 

numerically solving the finite-difference form of equation 3 

for systematically varied estimates of the air permeability of 

a given layer until the difference between measured and 

computed heads within or at the base of the layer is

minimized. Details of the methods by which the supplied 

data are prepared for use in the finite-difference form of 

equation 3, on the sequential analysis of data for each layer, 

on the numerical solution of equation 3, and on the 

automatic search procedure to select the final hydraulic 

conductivity estimates are all described in the section on the 

AIRK program.

As mentioned above, air permeability at the prevailing 

moisture content, rather than hydraulic conductivity, is 

measured by the method. However, hydraulic conductivity 

is frequently the parameter of interest to hydrologists. 

Consequently, the computer program is set up to convert 

initial user-supplied estimates of hydraulic conductivity to 

air permeability, to search for the best-fitting air-permeability value, and to reconvert that value to hydraulic conductivity before printing data out. Because English units are used in 

the program, the conversion described in the section on the 

AIRK program is confusing to follow. In metric units, 

however, the conversion is straightforward. If the medium 

were air-dry, the conversion would be made by the equation:

_*£_

"Vw (4)

where

and

A^=hydraulic conductivity, m/sec; 

/c=intrinsic permeability to air of the air-dry

medium, m2 ;

g=acceleration due to gravity, 9.81 m/sec2 ; 

v vv =kinematic viscosity of water, m2 /sec, at 

the prevailing temperature.

Because the materials in the unsaturated zone are not 

air-dry, the air permeability values determined by the above 

analyses are equal only to a fraction of the value at 

air-dryness, expressed by the formula

where Kra =ratio of air permeability of the porous 

medium at the prevailing moisture content 

to that of the air-dry medium.

Hence, the hydraulic conductivity determined by the 

analysis of pneumatic head data is approximate, and is given 

by the equation

(5)

where Ke =fie\d air-permeability based hydraulic 

conductivity; in m/sec.

As noted above, Kra generally ranges from 0.6 to 0.8 for 

materials drained to field moisture capacity.

The determined value of Ke is dependent on the prevailing 

soil temperature during the test, because the dynamic 

viscosity of air, which is temperature-dependent, occurs in 

the pneumatic diffusivity term. Moreover, Ke must be 

determined for some specified water temperature,
RESULTS 11

because v^ is also temperature-dependent. For the 

program, both vw and nfl are read in, unless specified 

otherwise, for a temperature of 15.6°C. The value for v^ 

was chosen because the U.S. Geological Survey has 

historically used 15.6° C (60° F) as the reference temperature 

for hydraulic conductivity measurements (Wenzel, 1942, p. 

7). The value for ^a was chosen because the mean annual air 

temperature at Lubbock, Tex., where much of the work was 

done, is about 16°C. Moreover, the viscosity of air varies 

only gradually with temperature, so the value generally 

should be adequate even in more northerly climes. The 

investigator does have the option of supplying values for 

viscosity at other temperatures if conditions warrant, 

however.

RESULTS

Field tests of the method were made at five sites in the 

Southern High Plains of Texas and at one site each at the 

Birch Creek Playa, located at the Idaho National Engineering Laboratory near Idaho Falls, Idaho; at the Central Great 

Plains Experiment Station near Akron, Colo.; in Scott 

County, Kans.; and at Nassau County Recharge Basin 306 

on Long Island, N. Y. Tests were described in detail below for 

the sites in Idaho and New York and a site in the Southern 

High Plains near Lubbock, Tex. A wide range of 

hydrogeologic conditions exists at these sites. Moreover, 

other data are available from these sites for comparison with 

the air permeability data.

BIRCH CREEK PLAYA, IDAHO

The method was applied at this site to determine the air 

permeability of the playa sediments, basalt flows, and a 

sedimentary layer between basalt flows, termed an interflow 

bed. These data were needed to support tests of the feasibility 

of emplacing radioactive waste gases in unsaturated basalt 

underlying the sediments (Schmalz, 1969, p. 1; Robertson, 

1969). Because air permeability, rather than hydraulic 

conductivity, is the parameter of interest for this test, 

complications arising from residual moisture content and 

structural changes of the medium due to wetting is not of 

concern. Thus the test represents a quite straightforward 

application of the air-permeability method.

The general stratigraphic sequence at Birch Creek playa, 

as described by Schmalz (1969), includes playa sediments 

overlying a thick sequence of basalt flows of the Snake River 

Group separated by one interflow sedimentary bed (fig. 5). 

The playa sediments range from 2 to 17 m thick, and consist 

of clay and silt with a small amount of sand. The basalts 

consist of interlayered individual flows ranging from several 

to 15 m thick. Permeability in the basalt is due to joints, 

fractures, and voids, and is very unevenly distributed. The 

sediments in the interflow bed consist of sand, silt, and clay, 

and probably are alluvial in origin. The top of the interflow 

bed lies from 23 to 37.5 m below land surface and ranges 

from 1.5 to 4.6 m thick.

Playa 

sediments

Upper 

basalt

Interflow 

sedimentary bed———'

Lower 

basalt

-100

H50

200

rr 

en 

o

2 

<

UJ 

GO

-0

-15

-30

-45

L 60

UJ

o

LJ 

GO

cn

CE 

UJ

No. 19

EXPLANATION

. DEPTH OF CENTER OF 

SAMPLING PORT, FEET

GRAVEL PACK

EXPANDING CEMENT

I PIEZOMETER

FIGURE 5. Sample logs and construction of wells tested at Idaho National 

Engineering Laboratory (modified from Schmalz, ed., 1969, fig. 9).

All permeability results for this site refer to intrinsic 

permeability in darcies2 , because the hydraulic conductivity 

of the materials is not of concern, and other published 

permeability values for the site are also in darcies.

Data of pneumatic head versus time were obtained from 

two piezometer nests (wells 10 and 19, Schmalz, 1969, p. 47) 

located about 200 m apart during a normal afternoon 

barometric decline of about 4 millibars in a 5-hour period on 

October 24, 1969. Screens for piezometers in each nest 

consist of specially designed cages that were originally

2 A one-centimeter cube of a porous medium has a permeability of one darcy if water at 68° F (20° C) 

(viscosity=l centipoise) will flow through at a rate of one cm3 /sec under a pressure drop of one 

atmosphere. At 60° F (15.6°C), a porous medium with a permeability of one darcy would have a 

hydraulic conductivity of 0.74 m/day.
12 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

installed to provide ports for sampling injected radioactive 

gas (Schmalz, 1969, fig. 10). These screens were gravel-packed and grouted with expanding cement as shown in figure 5. 

Piezometers in well 19 tapped materials below the interflow 

sedimentary zone, and piezometers in well 10 tapped 

materials above the zone. Consequently, data from the two 

nests were combined to make a single analysis for the air 

permeability of layers in the profile.

Analysis of the pneumatic-head data from the Birch Creek 

playa site was complicated by three factors. First, during the 

test, it was observed that pneumatic head differences were 

greater in well 19 at the 32-meter depth than at the 41-meter 

depth, indicating that the basalt surrounding the deeper cage 

is better connected to the surface by fissures, joints, and 

voids than that surrounding the shallower one. Such 

conditions differ from those assumed for the one-dimensional analysis, and consequently accurate simulation of 

flow in the lower basalt is not possible with the available 

data. Second, the gravel pack surrounding the screen at the 

24-meter depth (fig. 5) extends from 19 to 30 m, possibly 

allowing air circulation along the well bore. Third, the well 

screens at the 24- and 32-meter depths are not exactly at the 

upper and lower boundaries of the interflow zone.

Despite these complications, it was necessary to establish 

the effects of the underlying materials on the pneumatic 

heads at the base of the playa sediments before their 

permeability could be determined. Consequently, it was 

assumed that the lower basalt could be analyzed as a unit 

extending from 32 to 62 m, based on matching computed 

and measured pneumatic heads at the 58-meter depth. It was 

further assumed that the pneumatic head in the 24-meter 

piezometer was representative of that in the formation at 

that level. Finally, it was assumed that the material from 24 

to 32 m, which includes the interflow bed, could be analyzed 

as a single layer. The basalt layers were assumed to have an 

air-filled porosity of 0.05, based on data given by Robertson 

(1969, p. 13). The layer from 24 to 32 m was assumed to have 

an air-filled porosity of 0.15, which is the approximate 

thickness-weighted average of 0.05 for basalt and 0.2 

(assumed) for the interflow sedimentary zone.

Results of the analyses are shown in table 1. Although the 

results for the basalt and interflow sedimentary zones are

TABLE I. Summary of vertical air permeability values computed from 

pneumatic head data obtained from installations at the Birch Creek 

playa, Idaho

Nature of 

materials in 

layer

Near-surface playa 

sediments. 

Playa sediments ...............

Upper basalt ....................

Interflow sedimentary 

zone. 

Lower basalt ....................

Depth 

interval 

(meters)

0 - 3.0 

3.0-13.4

13.4-24.4

24.4-31.7 

31.7-60.7

Air-filled 

porosity

0.25 

.25

.05

.15 

.05

Air permeability 

(darcies)

6

2

5

2 

10

plausible, the complications mentioned above suggest that 

they be used with caution. Nonetheless, the analyses of data 

for piezometers at the base of and within the playa sediments 

described below result in a reasonably good match of 

measured and computed values (fig. 6). Consequently, it is 

assumed that the simulation of the basalt layer and interflow 

layers provided a good approximation of the boundary 

conditions prevailing at the base of the playa sediments, and 

that the results for the playa sediments are approximately 

correct.

Data from the 13-meter piezometer, at the base of the 

playa sediments, and from the 3-meter piezometer were 

analyzed to determine the air permeability of both a lower 

and an upper layer of playa sediments. The determined air 

permeabilities were 2 and 6 darcies (table 1) for the two 

layers, based on a laboratory-determined average air-filled 

porosity of the playa sediments of 0.25 for 16 samples. The 

higher permeability of the uppermost section of playa lake 

sediments probably results from the greater prevalence of 

desiccation cracks, root holes, and animal burrows at the 

shallow depth.

The field-determined air-permeability values compare 

favorably with laboratory measurements of air permeability 

at the prevailing moisture content on 11 drive cores, which 

ranged from 0.1 to 8.0 darcies. The geometric mean of 

permeability values, rather than the arithmetic or harmonic 

mean, has been suggested by Bouwer (1969) as providing a 

more valid estimate of the permeability of laterally and 

vertically heterogeneous materials, such as those comprising 

the Birch Creek playa. The geometric mean of the eleven 

sample permeability values is 0.96 darcies. Thus, the field 

value for the 3-13-meter layer of 2 darcies is about twice the 

geometric mean of the laboratory results, possibly because 

of compaction of the drive-cores, or because of the lack of 

secondary permeability features in the cores.

King (1968) also made estimates of air permeability of the 

playa sediments, based on injection pressure measurements 

during air injection at a constant rate into three auger holes. 

Of these three auger holes, two were open from a depth of 1.5 

to 3 m and one from a depth of 1.5 to 4.6m. King's (1968, p. 

39-40) values for several test replications in each of the three 

holes ranged from 5 to 18 darcies, and averaged 9 darcies. 

These values are slightly larger than that of 6 darcies 

determined for the 0-3-meter depth during this study. As 

King (1968, p. 33) noted, his values may be erroneously large 

because the effective radius of his wells may have been larger 

than assumed, which may account for the fact that the results 

of this method are lower than King's.

In conclusion, this method has yielded values for the 

vertical permeability of the Birch Creek playa sediments that 

are larger than those determined by laboratory analysis, but 

less than those determined by injection tests. These 

variations can be reasonably explained, and it is likely that 

the values determined by the present method are the most 

accurate of those determined for the vertical permeability of
RESULTS 13

cc •3 

o 

o: LU

S ,- 

o

U- LU 

O h-

(/) <

LU CC

X hO OQ

LU

- X

3 i

LU —

3= ID 

O £^i

0.12

0.10

0.08

0.06

0.04

0.02

NUMERICAL PIEZOMETER 

SIMULATION DEPTH.IN FEET

189

104

80

44

10

DATA

•

A 

O 

D 

x 

A

K884

 883 en

cc

CO

 882

LU

X

— 881

880

50 100 150 200 250 300 350 

TIME, IN MINUTES, SINCE 12=00 NOON OCTOBER 24, 1969

400

FIGURE 6. - Comparison of pneumatic heads measured on October 24, 1969, at five depths in the unsaturated zone at the Birch Creek playa,

Idaho, site to the best-fitting numerical simulation.

the playa sediments. The described method should be highly 

useful for evaluating the feasibility of temporary storage of 

gas in the unsaturated zone at other sites where the materials 

are well-drained.

LUBBOCK, TEXAS, AIRPORT SITE

A site adjacent to the Lubbock, Tex., Regional Airport 

was selected to test the utility of estimates of hydraulic 

conductivity determined by the described method to predict 

infiltration rates and the effects of perching layers during 

water-spreading operations. The Lubbock site is in a 

semiarid area, and problems arising from high residual 

moisture content were not anticipated, nor did they occur. 

However, the geologic section within the unsaturated zone at 

this site does contain significant quantities of clay minerals 

that were expected to swell or disperse upon wetting.

A great many field and laboratory data were obtained on 

the pneumatic and hydraulic properties and of correlative 

geophysical and lithologic properties of the unsaturated 

zone at this site. Air permeability values and equivalent 

hydraulic conductivities were determined from data obtained at three installations. In addition, laboratory determinations of volumetric moisture content, porosity, 

hydraulic conductivity, air permeability (both at the prevailing moisture content and at oven-dryness), grain-size distribution, pore-size distribution, carbonate content, and clay 

mineralogy were obtained on numerous drive cores from 

two wells, and on wireline cores from three wells. Extensive 

suites of geophysical logs were also obtained from several 

observation wells. Details of the laboratory tests and of 

geophysical logging at this site are described by Keys and 

Brown (1971). Finally, two recharge-basin experiments 

were conducted, and numerous data on infiltration rates, 

piezometric head, and neutron-log-based moisture content 

were obtained. Locations of instrumentation relative to the 

recharge ponds are shown in figure 7.

The site is underlain by near-surface deposits of 

Pleistocene and Holocene age and by the Ogallala Formation, which extends to about 55 m in depth. Before spreading 

operations were begun, the section was unsaturated to a 

depth of about 38 m below land surface. A lithologic log 

based on examination of drive cores obtained at one site 

(modified from Keys and Brown, 1971) is shown in figure 8.
14 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

/^\^2^—/o

w

V North Basin /

4N^ I?TM 13 y

v°^ ^• df,0 | 8tr 14 X *b. I9 0 3C M l8 .

V *° ^s" 

^^^°^^———^

®°<^^ 7G& >

^^

<^^^ South Basin

X*-^

^^-50 ————

o 5

>

/

"T

O22 20CH

a loou_ 

50-

0- 

ensiometers

15 

O

-50 

to

IT 

UJ

1- 

LU

S

-20 

-10 

-0

BELOW LAND SURFACE

O Oi1

co 20- 

o:

UJ iUl

S 

z 25-

x"

i0.

UJ Q 30- 

35- 

40-

Land surface o— - — • ——

20 —

. UJ

o 

<

u. 

cr

" 40- 

o

2 

<

_l

5 

o

ul 60—

CD

1- 

LU 

LU 

U.

^ 80 —

X

ia.

UJ

o 

100—

120 —

'•±:±e

&*

••• ,o~

;o-.£'

7TM" O,

*r^ m

¥&

V:ffi 

5££

• 'O'-'

6: '6. '• • o'.- '

•o • ' ' • 

• • • o

• o' . —

-.•'• .<?'.

Sand, silt, and elay, with 

nodular concretion! of calighe

Caliche and calcium carbonatecemented sand (massive calich 

zone of Gile and others, 1966

Solution openings

Clay

Sand with numerous caliche 

nodules

Sand with few caliche nodules

Clay, red, with caliche streaks 

_ Sand with caliche grain coatings

Clay 

Sand with few caliche nodules

__- Clay ^ 

HV^-Sand /

""v-Sand i Few caliche nodules 

-^Clay 1 

^-Sand J

EXPLANATION 

O OBSERVATION WELL

• A AIR-PERMEABILITY 

INSTALLATION

A PIEZOMETER NEST (South 

basin nest is unnumbered)

£C, POROUS-CUP INSTALLATIONS

FOR WATER-QUALITY SAMPLING

—— 40—— WATER-LEVEL CONTOUR — 

Interval 10 feet (3.05 meters)

FIGURE 1. — Locations of air-permeability installations, observation wells, 

water-level piezometer nests, and spreading basins used in recharge 

experiments near the Lubbock, Tex., Regional Airport. Also shown 

are contours, as of July 1975, of the depth, in feet, from land surface to 

the top of the perched ground-water mound created by infiltration 

from the north basin.

Also shown is the occurrence of caliche based on the 

classification of Gile and others (1966).

Three installations for determining effective pneumatic 

diffusivity, herein labeled sites A, B, and C, (fig. 7) were each 

completed somewhat differently. Site A, the first installed, 

consisted of five piezometers (installed at depths listed in 

table 2) in three holes about 2.5 m apart on a line. These holes

FIGURE 8. Sample log of test well at Lubbock, Tex., airport spreading 

site, showing lithology and nature of caliche occurrence (modified 

from Keys and Brown, 1971).

were grouted with chemical grout. Readings of pressure 

versus time were first obtained at site A on March 7-8,1969, 

during passage of a cold front that raised atmospheric 

pressure by about 17 millibars in a 12-hour period. Data for 

the period 16:00 March 7 to 03:00 March 8, 1969, were 

initially analyzed by use of the electric analog network to 

give the curve match shown in figure 3.

Hydraulic conductivity estimates based on this analysis 

seemed unreasonably high, based on known values of 

hydraulic conductivity for the types of materials present. It

TABLE 2. Depths of piezometers in three installations used to 

determine vertical air permeability at the Lubbock, Tex. 

airport spreading site

Depth of center of screen below land surface (meters)

Piezometer No. Site A Site B SiteC

1

2

3

4

5

6

32

25

16

10

5

35

31

26

22

16

10

35

30

26

22

15

9
RESULTS 15

was originally hypothesized that the high results were 

obtained because a gap in the backfill in an observation well 

6 m east of the easternmost piezometers, detected from 

geophysical logs (Keys and Brown, 1971, fig. 7), allowed 

rapid circulation of air between levels. However, results of 

subsequent numerical analysis of the data, listed in table 3, 

indicate that the problem instead resulted from the fact that 

the lowermost layer had been assumed to have a much lower 

air-filled porosity than the overlying layers. Because of the 

problems associated with the electric analog analysis 

described above, the air permeability and equivalent 

hydraulic conductivity values £or the overlying layers were 

computed to be much too high. In fact, the data for this 

installation are in reasonable agreement with those from the 

other installations.

TABLE 3. Summary of hydraulic conductivity values computed using 

pneumatic head data obtained from installations at the Lubbock, Tex., 

airport spreading site, 1969-72

Depth 1

interval

(meters)

0-10

10-16

16-22

22-26

26-30

30-38

"d

0.20

.20

.15

.21

.03

.22

Site A

March 7-8,

1969

(m/day)

12

3

2

.3

Site B

January 14,

1970

(m/day)

7

5

3

3

.2

.2

May 17,

1972

(m/day)

7

3

3

1

.1

.1

SiteC

January 15,

1970

(m/day)

6

7

4

3

3.2

May 17,

1972

(m/day)

6

5

4

2

.1

'The depth intervals are approximate, because screens were placed at somewhat different depths at 

each site. Exact depth intervals are shown in table 2.

!"rf=air"filled porosity, dimensionless, as determined from core data and neutron logs.

'At Site C, the screen immediately beneath the clay layer at the 30-meter depth was plugged, so the 

material in the depth interval from 26-38 m was analyzed as a single uniform layer.

Because the data for site A were originally considered to 

be in error, site B was installed about 76 m northwest of site 

A (fig. 7). For this installation, six piezometers were installed 

in a single hole to depths shown in table 2 and the intervals 

between screens were grouted with chemical grout. Before 

pressure-time observations were obtained at site B, it was 

found that expanding cement was more effective in 

providing a seal in the well bore than chemical grout. 

Consequently, site C was installed about 30 m east of site B, 

with six piezometers (depths listed in table 2) installed in a 

single hole and grouted with expanding cement.

Observations of pressure versus time were obtained 

during normal afternoon diurnal barometric declines at site 

B on January 14, 1970, and May 17, 1972, and at site C on 

January 15, 1970, and May 17, 1972. Results of the 

determinations are shown in table 3. In addition, a comparison of measured versus computed pneumatic heads for 

the data obtained May 17,1972, at site C is shown in figure 9. 

This comparison is typical of those for the airport spreading 

site data, and shows that the computed data match the 

measurements well.

The results of the analyses show that the upper layer is 

quite permeable to air, possibly because of root holes and 

desiccation cracks. However, the layer from 9 to 15 m, which 

includes a massive caliche zone, is nearly as permeable. In 

fact, the air-permeability-determined hydraulic conductivity 

results for the massive caliche zone at all six sites where the 

Miocene Ogallala Formation was tested are similarly high. 

These results suggest that the commonly held belief that the 

massive caliche layer, or"caprock," is so nearly impermeable 

as to impede recharge (Cronin, 1964, p. 38-39, for example) 

is incorrect. Permeability of the massive caliche zone is 

almost undoubtedly the result of cracks, joints, and solution 

openings, as core examination indicates that the primary or 

intergranular permeability of the caliche is indeed low.

Air-permeability-based hydraulic conductivities of the 

materials in the two layers from 15 to 21 and 21 to 26 m are 

somewhat lower than are those of the caliche bed, but are 

reasonable for the materials tested.

The air-permeability-based hydraulic conductivity values 

for the intervals from 26 to 30 m and from 30 to 38 m are 

considerably higher than was anticipated. Low values were 

expected for these intervals because the upper interval 

contains a 2.4-meter-thick clay bed and the lower interval 

contains interbedded sand and clay beds. Cores show the 

clay beds to have very low hydraulic conductivities and 

geophysical logs indicate that the 2.4-meter-thick clay bed 

and some of the thinner underlying clay beds are continuous 

over the spreading basin area. Initially, it was assumed that 

the values determined by the air-permeability technique 

reflected a permeability to air that would be substantially 

reduced when the clay minerals in the clay beds swelled or 

disaggregated upon becoming water saturated. However, 

the results of subsequent recharge experiments indicate that 

the higher values may instead reflect the effects of cracks or 

fractures in the clay beds on their permeability.

COMPARISON OF RESULTS FOR DIFFERENT 

SITES AND DATES

Field air-permeability based hydraulic conductivities 

(table 3) are in reasonable agreement for all three sites. This 

agreement exists despite the desiccation and cracking of the 

chemical grout at sites A and B, as discovered by backhoe 

excavation. Nonetheless, although the deterioration of the 

chemical grout apparently did not significantly affect the 

results of the tests, expanding cement would be preferred for 

this application as it provides greater ensurance that air 

circulation is not occurring through the bore hole.

Analysis of data for sites B and C obtained in January 

1970 resulted in generally higher computed values of 

hydraulic conductivity than did the data obtained in May 

1972. These differences may result from changes in moisture 

content in the profile, from air circulation in a nearby 

temporarily open test hole, or from experimental error. 

Although repeated neutron logs showed no measurable 

changes in moisture content in the immediate vicinity of sites
16 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

0.96

rr

r> 

o 

rr

LlJ

^ _

Q 

U. LiJ

O Hcn <

LU rr

x HO CD

•z. z> "~ en

? en

9 °

< z

UJ —

X to

o £L

UJ •z. 

CL

— 913

 912

0.90

0.88 —

0.86 —

0.84

tr 

<

00

— 911 ±

LU

X

o

— 910

— 909

50 100 150 200 250 300 350 400

TIME, IN MINUTES AFTER I2M5 PM MAY 17, 1972 

FIGURE 9. Comparison of pneumatic heads measured May 17,1972, at five depths at Lubbock, Tex., airport site B to best-fitting numerical simulation.

B and C, a recharge test was in progress on the south basin 

during the time readings were obtained in 1972. Consequently, a possibility exists that the readings in May were affected 

by moisture changes in the general vicinity of the site. On the 

other hand, test drilling was in progress at the spreading 

basin site during January 1970. Hence, a possibility exists 

that the higher values determined from the January 1970 

data result from air circulation in an open test hole. Whether 

the differences are indeed due to moisture content changes, 

to circulation in an open bore hole, or to experimental error 

is not known, but possibly could be determined from a third 

set of data.

COMPARISON OF RESULTS WITH RECHARGE DATA

Two recharge experiments have been conducted at the 

Lubbock, Tex., site, the results of which were analyzed to 

evaluate the utility of the air-permeability data for predicting 

recharge basin performance. The first experiment was 

conducted from April 12, 1972, to June 11, 1973, on the 0.4

I hectare pond located south of site A (fig. 7). The second 

experiment was conducted from April 28 to November 17, 

1975, on the 0.4 hectare pond north of site A. Both 

experiments were performed using water from Lake 

Meredith, a reservoir on the Canadian River near Amarillo, 

Tex. Hydraulic data collected during the experiments 

include inflow, pond stage, water table altitudes, and 

periodic neutron moisture logs at several locations, and 

piezometric readings of perched and semiperched 

ground-water mounds at six levels at each of five sites. 

Tensiometer data were obtained in the near-surface 

materials during the first test, but not the second.

The recharge-experiment data would best be analyzed by 

numerical simulation to include the effects of partial 

saturation on the movement of water from the spreading 

pond to the ground-water reservoir. Such an analysis, based 

on an approach modified from Green and others (1970), is 

planned. However, as these results are not yet available, 

saturated flow theory was used to provide estimates of the
RESULTS 17

hydraulic conductivity of materials in various zones above 

the water table. Such an analysis is possible because the 

percolating recharge water perched almost to land surface 

during the test. Because piezometer data at six different 

levels within the perched mound were available at four sites 

for the second (North Basin) test, versus only one site for the 

first test, the analyses described below are for the North 

Basin test data.

Estimates of hydraulic conductivity were made for each 

zone bracketed by piezometers based on infiltration rate, 

neutron-log moisture content, and perched zone water level 

data collected at the end of July 1975. This period was 

selected because recharge had been in progress for about 90 

days, allowing for much of the air that was initially 

entrapped in the unsaturated porous medium to be dissolved 

or dissipated. In addition, the infiltration rate was near its 

maximum at that time, and pressure heads within the 

perched mound were stable. Such stability is important, 

because storage changes within the mound are assumed 

negligible for the analysis.

Estimates of hydraulic conductivity from the recharge 

data were based on an application of Darcy's Law,

K=IA (6)

whereA:=hydraulic conductivity, m/day;

Q ^average inflow rate to the spreading basin, m3 / day;

/^hydraulic gradient, m/m;

and ,4=area of the perched mound at the top of the zone 

being analyzed, m2 .

Because the hydraulic conductivity term includes the viscosity of water at the prevailing temperature in its formulation, 

the above equation is correct only for isothermal flow. To 

partially compensate for this, the infiltration rate was 

adjusted by the ratio of the dynamic viscosity at the 

prevailing temperature to that at 15.6°C. With this correction, the infiltration rate at the end of July 1975 was 

approximately stable at about 0.9 m/day.

The hydraulic gradient for each zone bracketed by 

piezometers was computed by dividing the head difference 

between piezometers by the thickness of material separating 

them. The head gradients across a given zone showed 

substantial variation among the four sites, apparently 

because of areal heterogeneity of materials in each layer. For 

this analysis, the head gradient was assumed to be the 

average of that for the four sites.

The areal extent of the perched mound at the top of each 

layer was determined from neutron-log data obtained in July 

and early August. (The wells were not all logged at the same 

time). As shown in figure 7, the mound is not symmetrical 

about the center of the pond, and data are too sparse to 

accurately determine the area of mounding. However, rough 

estimates of the area of mounding were made by assuming 

radial symmetry from the center of the pond, and interpolating mound heights along a line through wells 1 and 5. 

The results of the analyses of the piezometer and 

perched-mound area data are summarized and compared 

with the air-permeability data in table 4. However, of the 

seven tabulated values, only the middle 5 were determined by 

the analyses described above. No formal analysis was made 

of the layer from 0 to 2 m, because one of the piezometers 

open at the 2-meter level was dry. However, as the hydraulic 

gradient was near unity in the other three, it is reasonable to 

assume that the hydraulic conductivity of the layer is 

approximately equal to the infiltration rate of about 1 

m/day. Also, no attempt was made to determine the 

hydraulic conductivity of the layer from 32 to 38 m, because 

the perched mound at that level extended beyond the area 

for which neutron-log data were available. However, the 

hydraulic gradient across the zone was slightly steeper than 

that in the 26 to 32 m zone, indicating that hydraulic 

conductivity of the lower zone might be only slightly less 

than that in the 26 to 32 m zone.

TABLE 4. Comparison of hydraulic conductivity values estimated from 

recharge data, from the air-permeability determinations, and from 

laboratory analyses of drive cores and wireline cores

Hydraulic conductivity (m/day)

Depth

interval

(m)

0-6

2-8

8-17

17-21

21-26

26-32

32-38

Recharge

experiment

data

S-*J 1

2

1

4

1

.1

<A

Field air permeability data

(average for three sites)

Uncorrected Corrected 1

6

5

3

2

.2

.2

1.4

.8

.5

.3

.01

.02

Laboratory data

(geometric mean)

0.2

.01

.09

.004

.06

.006

.02

(4)2

(9)

(5)

(2)

(5)

(11)

(10)

'"Corrected" values have been computed using the scheme described in the text to account for the 

structural effects of wetting on the permeability of the materials. 

2 Numbers in parentheses refer to the number of sample values included in the mean.

As a check of the above values, a different analysis was 

made of data obtained May 22-28, 1975 to determine the 

hydraulic conductivity of the 26-32-meter interval. This time 

period was chosen because saturated flow had been established through the depth interval, but the wetting front 

had not yet reached the water table, at least over most of the 

recharge area. Thus, saturated flow theory could be used to 

establish the hydraulic gradient, and changes in storage 

beneath the interval could be computed from neutron logs to 

provide an estimate of vertical flux.

For the analysis, the average hydraulic gradient across the 

26-32-meter zone was determined for the one-week period at 

each of the four piezometer nests, as described above. As for 

the previous analysis, the gradient showed a wide range over 

the area (from 0.24 to 2.1 m/m). The arithmetic average 

hydraulic gradient of 1.3 m/m was used for the analysis.

Flux through the 26-32-meter interval was estimated by 

determining the change in moisture content in the un-
18 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

saturated zone from May 22 to May 28 using calibrated 

neutron logs for observation wells 13, 20, 21, and 22 (fig. 7) 

and by dividing by the time interval between the logs. For the 

analysis, storage changes were computed in the depth 

interval from about 32 m (the approximate location of the 

wetting front on May 22) and the water table, which was at 

the 37-38-meter depth on the same date. These values of flux 

determined for the individual wells also showed considerable 

variation. They ranged from about 0.03 to 0.12 m/day and 

averaged 0.07 m/day.

Use of the average hydraulic gradient and average flux 

rate in Darcy's Law resulted in a computed hydraulic 

conductivity of 0.06 m/day. Assuming that air bubbles were 

trapped during the period in the pores or cracks that provide 

the permeability of the 26-32-meter layer it is not surprising 

that the estimate is somewhat lower than that for the August 

data.

The above recharge-experiment based values of hydraulic 

conductivity are very approximate, both because of assumptions made in the analyses, and because of the apparent wide 

areal variation in hydraulic conductivity of materials at the 

site. Nonetheless, the results are accurate enough to provide 

meaningful comparisons with the air-permeability data.

Comparison of the recharge-based and 

air-permeability-based hydraulic conductivities (table 4) is 

good for the four lower layers, but the air-permeability-based values for the 0-8- and 8-17-meter layers are about 3 to 4 

times higher than those based on recharge data. The 

differences for the uppermost layer may be due to reduction 

in the intrinsic permeability of the material upon wetting. 

Moisture content in at least the upper part of this layer is 

reduced much below field capacity by evapotranspiration, 

providing an opportunity for desiccation cracks that would 

close upon wetting. In addition, root holes in the material in 

this zone might be open to air flow, but might be blocked to 

water flow by transported sediment or by swelling of clay 

minerals. Moreover, results from other air-permeability 

experiments conducted in the High Plains areas of Texas, 

Kansas, and Colorado all indicate Ke values for the surface 

layer of the same general magnitude as those determined at 

the Lubbock site. However, other recharge experiments 

(Prill, 1976; Aronovici and others, 1972; Gillespie and 

Hargedine, written commun., 1976) all suggest that the 

surface layer hydraulic conductivity is about 1 to 3 m/day, in 

general agreement with the recharge experiment data 

described above. Hence, for the surface layer, hydraulic 

conductivities estimated by the air-permeability method 

appear to be consistently higher than those based on ponded 

infiltration tests.

The lack of agreement between the Ke and 

recharge-experiment-based hydraulic conductivity values 

for the 8-17-meter depth is less readily explained. This zone 

is deep enough below land surface that soil moisture should 

be relatively unaffected by evapotranspiration. Moreover,

the zone includes a massive caliche layer, and its permeability probably is due to secondary porosity. Such permeability 

is not likely to be affected by wetting of the medium. A 

possible explanation for the discrepancy between the 

air-permeability and recharge results exists in the extreme 

heterogeneity of the layer with respect to permeability.

CORRECTION FOR STRUCTURAL ALTERATION OF 

MATERIALS BY WETTING

The good agreement between the recharge experiment and 

air-permeability-based hydraulic conductivities for the bottom four layers was a surprise. It had been anticipated that 

the air-permeability-based hydraulic conductivity values 

would all be too high, because they did not include the effects 

on permeability of structural changes of the medium due to 

wetting. Because of these anticipated effects, a correction 

scheme, based on one used successfully by Reeve and 

Tamaddoni (1965) for an analogous problem, was tested. 

Basically, the scheme relies upon a laboratory investigation 

of the effects of wetting upon permeability. These laboratory 

data are used in the equation:

X RK. (7)

where ATc=corrected hydraulic conductivity;

A^=hydraulic conductivity as computed from

field-determined air permeability;

and RK is the ratio of (1) hydraulic conductivity measured 

on a saturated laboratory sample to (2) the hydraulic 

conductivity of the sample computed from its air permeability as measured at the prevailing moisture content.

To test the correction method, both the air permeability at 

the prevailing moisture content and the saturated hydraulic 

conductivity were measured on 30 cores of the Ogallala 

Formation obtained at the Lubbock site and at other sites in 

the Southern High Plains of Texas. The results of these tests 

(table 5) show that the intrinsic permeability of these samples 

as measured with water collected from the underlying 

saturated part of the Ogallala Formation ranged from about 

1 / 500th to slightly more than their intrinsic permeability as 

measured with air at the prevailing moisture content.

In an attempt to further generalize the correction factor, 

the ratio of water permeability to air permeability, RK, was 

correlated with the clay mineral content of the samples, as 

listed in table 5. Had the overall correction scheme been 

successful, the correlation equation would have provided a 

means to correct field-determined air permeability for the 

effects of wetting, based on x-ray mineralogical data from 

samples. The analysis resulted in the correlation equation:

/MT=0.3exp(-5.7CMQ (8)

where CMC is the quantity of expanding-lattice clay 

minerals of the montmorillonite group, expressed as a 

decimal fraction of the total mineral content.
RESULTS 19

TABLE 5.—Laboratory data for 30 cores from the Ogallala Formation in the 

Southern High Plains of Texas, showing comparisons of measured 

hydraulic conductivity to those computed from air permeability

Sample No. 1 Depth Air-permeability- Water-based

(m)2 based hydraulic hydraulic

conductivity conductivity

(m/day)] (m/day)4

Clay mineral 

content6

70 Tex 48

52a

52b

58

64

25.5

27.3

27.3

30.6

33.7

0.7

.5

.2

1.4

.4

0.9

.04

.006

.6

.3

1.3

.08

.03

.4

.8

0.10

.15

.15

.4

.1

67a 

67b 

71 

87a

87b

Unnumbered7 

Unnumbered7

85

88

90a

90b

71 Tex 98 

100 

102 

107

155

156

157

158

159

170

171

172

173

174

35.5

35.5

36.3

30.0

30.0

19.7

31.7

34.4

34.4

1.7

6.7

9.1

26.7

19.2

22.7

24.4

30.0

32.0

20.0

22.6

24.2

28.5

29.9

.01 

2.9 

7,9

.6

1.0

.0026

.0012

.024

.33

.002

.03

.034 

2.1

.5 

7.6

.07

4.0 

1.7

.002

.02

.1

6.4 

2.4 

1.2 

.16

.00002

1.2

3.2

.2

.6

.00026

.00004

.0007

.03

.00004

.0003

.0001

.02

.016

5.3

.013

1.7

.09

.0002

.015

.0002

1.3 

.2 

.04 

.0008

.002 

.4 

.4 

.3

.01

.03

.03

.1

.02

.01 

.003 

.01 

.03

.4

.05

.1

.7

.002

.2 

.08 

.03 

.005

.1

.15

.15

.15

.15

.6 

.6

.2 

.2 

.6

.6

.40

.15

.10

.00

.15

.10

.25 

.40 

.15 

.50

.10 

.25 

.20 

.40

'Sample numbers are those of the Hydrologic Laboratory, U.S. Geological Survey, Lakewood, 

Colo. Samples in the 70 Tex series were obtained at the Lubbock Airport Spreading Site, those 

numbered 71 Tex 98-107 are from a site near Caprqck, Tex., and those numbered 71 Tex 155-174 are 

from a site near Stanton, Tex.

2Depth is to the top of the 0.15-meter drive core from which the sample was taken for analysis. A 

lowercase a or b after the sample number indicates that two samples were analyzed from that core.

3 As measured at the prevailing moisture content.

4Samples from the Lubbock and Caprock sites were tested using water from a well tapping the 

Ogallala Formation at the Lubbock Airport site, and those from the Stanton site were tested using 

water from a well near that site.

5 RK is the ratio of hydraulic conductivity as determined by the air-permeability measurments to 

that measured using water, and is the same as the ratio of intrinsic permeability that would be 

determined by use of the two fluids.

6As a decimal fraction of the total mineral content. These fractions, which are approximate, were 

determined by Barbara Andersofi of the Hydrologic Laboratory, U.S. Geological Survey, 

Lakewood, Colo.

'These samples are both from the clay bed at approximately 27 m in depth at the Lubbock Airport 

Spreading Site,

The correlation coefficient for . this analysis is 0.54, 

indicating that, although the correlation is significant at the 

1 percent level (Crow and others, p. 240), only about 30 

percent of the variation in the ratio RKcan be accounted for 

by variations in clay mineral content. This correlation is 

poorer than that determined by Reeve and others (1954, p. 

132) for the correlation of the ratio of air- to distilled-water 

permeability with exchangeable sodium and cation-exchange capacity, probably in part because clay chemistry

was not considered in this analysis, and in part because the 

clay-mineral contents are approximate.

Corrected values for hydraulic conductivity were computed for each of the six intervals for which air-permeability 

values were determined, based on a combination of 

equations 7 and 8:

Kc=Kex0.3exp(-5.1CMQ (9)

To apply equation 9, the arithmetic average clay content was 

determined from data on all of the samples collected and 

tested from a given layer. This average value was used to 

represent CMC in the equation. Use of the arithmetic 

average was arbitrary, and may not be completely appropriate because it assigns too little weight to the effects of 

wetting on the discrete clay layers within each interval. 

However, the worth of the correction for the effects of 

wetting would be questionable no matter what averaging 

scheme was used.

Results of applying equation 9 to compute corrected 

hydraulic conductivities are listed in table 4. The results 

appear reasonable for the top two layers, but not for the 

bottom four layers. Moreover, the apparent agreement 

between the corrected air-permeability results and the 

recharge-based hydraulic conductivity for the surface layer 

may be fortuitous. The surface layer probably is indeed more 

sensitive to the effects of wetting than the deeper layers, as 

described above. However, the laboratory-based correction 

factor was determined on samples mainly collected at depth, 

and thus is not necessarily relevant for the surface layer data. 

Agreement of the corrected air-permeability-based 

hydraulic conductivity with the recharge-based value for the 

8-15-meter layer also seems fortuitous, as the massive caliche 

in this zone should not be structurally affected by wetting.

Failure of the correction scheme to improve the 

air-permeability-based hydraulic conductivity estimates for 

the bottom four layers may result from flow occurring 

through secondary porosity features, such as cracks and 

joints, rather than through the intergranular porosity. Thus 

the laboratory results, which include effects of wetting on 

permeability related to intergranular porosity only, might 

not be relevant.

Additional experience with the air-permeability technique 

will be needed to reconcile whether, or under what conditions, a correction for the effects of wetting is warranted 

for layers at depth. On the other hand, sufficient experience 

has been gained to indicate that Rvalues of the near-surface 

materials generally would need to be corrected for the effects 

of wetting. Such a correction could be practically developed 

if the ratio of air-permeability-based to recharge-based 

hydraulic conductivity could be correlated with some readily 

measurable property of the near-surface materials.

COMPARISON OF RESULTS WITH LABORATORY DATA

As described above, a number of sequential laboratory 

analyses of air permeability at the prevailing moisture
20 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

content and of hydraulic conductivity were made to assess 

the effects of wetting on permeability, as tabulated in table 5. 

In addition, hydraulic conductivity only was measured on 

many more cores obtained from the Lubbock Airport site to 

provide an independent assessment of the recharge potential 

at the site. The results of these laboratory analyses are 

compared with the field air-permeability and recharge-experiment results in table 4.

In general, the field-determined Ke values are somewhat 

larger than the laboratory values of hydraulic conductivity 

based on their air permeability at the prevailing moisture 

content. For example, in the 21-26-meter layer, the 

field-determined air-permeability-based hydraulic conductivity for the May 1972 data is about 4 times higher than the 

geometric mean of the laboratory values determined from 

drive-core samples of this unit. Comparison of hydraulic 

conductivity values for the 30-38-meter layer indicates that 

the geometric mean of the laboratory data is about the same 

as the air-permeability based hydraulic conductivity. 

However, this layer is much more strongly stratified than the 

21-26-meter layer, and consists of alternating thin beds of 

relatively clean sand and of clay. Under these conditions, the 

harmonic mean of the laboratory values may represent the 

hydraulic conductivity of the layer as a whole more 

adequately than the geometric or arithmetic means (Terzaghi, 1943, p. 243-244). The harmonic mean of the sample 

hydraulic conductivities, assuming that each sample 

represents a bed of equal thickness, is about one-fifth to 

one-tenth the air-permeability-based hydraulic conductivity. Laboratory air-permeability data are available for too 

few samples to compute averages for the other layers for 

which field data are available. However, the laboratory 

values for samples within these layers are somewhat less than 

the field values.

These comparisons indicate that the air permeability 

values determined in the field are several times higher than 

those determined in the laboratory. These differences are not 

surprising, because of the inevitable compaction of 

drive-core samples during coring operations; and because of 

the small sample size of the core, which may miss secondary 

permeability features. In fact, if such features had appeared 

in a core, they probably would have been attributed to 

disturbance by coring, and the core would not have been 

analyzed.

The geometric mean of laboratory hydraulic conductivity 

values from drive cores obtained in each zone are also shown 

in table 4 to provide a comparison between the air 

permeability method and laboratory methods for estimating 

hydraulic conductivity. The laboratory values were obtained 

on drive cores from wells 3 and 11 by the U.S. Geological 

Survey Hydrologic Laboratory, Lakewood, Colo. Part of 

the analyses were made with Lakewood tap water, and part 

with Canadian River Aqueduct water. However, subsequent 

experiments by the author indicate that the hydraulic

conductivity values should be the same using water from 

either source.

The laboratory-based values are all substantially smaller 

than those based on either the air-permeability-data or the 

recharge-data values, probably in part because of compaction of the samples during coring, and in part because effects 

of secondary permeability are not represented in the 

laboratory samples. The laboratory values are so much 

lower than those determined from the recharge experiment 

that their use in predicting recharge pond performance 

would be seriously misleading.

DISCUSSION OF RESULTS

The comparisons made among the various data obtained 

at the Lubbock Airport site indicate that the results of the 

air-permeability method are repeatable, and that the method 

can provide useful data for predicting the performance of 

recharge basins. The comparisons also indicate that the 

permeability of materials at depth is not affected by wetting 

nearly to the extent anticipated. Although the lack of 

permeability sensitivity may not be general, the results of 

these tests are encouraging that the method can be a useful 

tool for assessing the hydraulic conductivity in the unsaturated zone, particularly if the method is used in 

conjuction with other methods to evaluate the hydraulic 

conductivity of near-surface materials.

LONG ISLAND, NEW YORK

A test of the method was made on Long Island, N.Y., to 

evaluate its applicability to determine hydraulic conductivities in a humid region, where residual moisture content in 

the unsaturated zone might be high. The test, which included 

two installations, was made at a spreading basin (Nassau 

County Recharge Basin 306) for which much information on 

lithology, infiltration, and deep percolation was available 

(Prill and Aaronson, 1973). Materials in the unsaturated 

zone at this site (fig. 10) are glacial in origin, and include two

Site Site 2

FIGURE 10. Cross-section showing assumed uniform layers of unsaturated materials at Nassau County, Long Island Recharge Basin 

306 (modified from Prill and Aaronson, 1973, fig. 2). Not to horizontal 

scale.
RESULTS 21

layers of sand and gravel that are separated by a gravelly 

sandy loam and are overlain by a surface loam.

The first installation (site 1) for determining air 

permeability, which consisted of four piezometers, was made 

in the basin in June 1971. Readings were obtained on July 1, 

1971, during an afternoon diurnal pressure change of about 

3 millibars. These observations were made at a time when 

moisture content of the unsaturated materials was relatively 

low, as there had been only 50 mm of rain in June and no 

storm runoff had collected in the basin during that month. 

Moreover, very little storm runoff had collected in the basin 

during the previous April and May.

Data from the installation were initially analyzed by 

electric analog to determine air permeability and the 

equivalent hydraulic conductivity for four layers, including 

the lower sand and gravel, two layers differentiated on the 

basis of moisture content in the gravelly sandy loam, and a 

layer consisting of the surface loam and upper sand and 

gravel combined. The results of this initial analysis were 

greatly in error. In addition to the erroneous handling of the 

effects of variations in air-filled porosity between layers, as 

noted previously, initial conditions were inadequately 

represented for these data. Substantial differences in 

pneumatic head between depths occurred at the start of this 

test, but the analog system was designed only to handle 

initial conditions of constant pneumatic head with depth. An 

attempt was made to overcome this problem by manipulation of the input head readings, but the subsequent digital 

similation indicates that the attempt was not successful. 

Because of these problems, it was concluded that data from 

the site had given values for hydraulic conductivity that were 

much too low, and that only a poor match had been 

obtained.

Subsequent numerical analysis of this set of data resulted 

in hydraulic conductivity values (table 6) that are quite 

plausible, and the match of computed to measured 

pneumatic head data (fig. 11) is quite good. Moreover, the 

Ke value for the lower part of the sandy gravelly loam agrees 

favorably with the hydraulic conductivity value of 0.03 

m/day determined for this layer by Prill and Aaronson 

(1973, p. 741) from infiltration data. However,, the Ke value 

of 0.06 m/day for the 0-1.2-meter layer is about five times 

smaller than their infiltration-based value of 0.3 m/day for 

this layer.

Because it had been assumed that erroneously low results 

had been obtained for the installation inside the pit, the 

second installation was made outside and about 3 m from the 

edge of the pit. It was assumed that the materials underlying 

the area outside the pit would be much better drained than 

those beneath the pit, because they would not receive deep 

percolation from ponded infiltration.

Data were obtained at this site (site 2 in table 6) on June 

29, 1972, during an afternoon diurnal barometric change of 

1.8 millibars. The data for this test are poor, because only 

very small pneumatic head differences with depth occurred,

TABLE 6. Summary of analyses of pneumatic head data obtained from 

installations at the Nassau County, Long Island Recharge Basin 306

[Ke, air permeability based hydraulic conductivity, m/day; n ,, effective air-filled porosity,

dimensionless]

Nature of 

materials in 

layer

Site 1

Depth interval Kg r 

(m)

Site 2

i^ Depth interval Kg 

(m)

n d

Surface loam .....

Upper gravelly

sandy loam.

Lower gavelly

sandy loam.

Lower sand and

gravel.

'0-1.4

1.4-3.7

3.7-5.2

5.2-17.7

0.1

.9

.03

14

0.15

.1

.05

.2

0-1.7

2 1.7-4.7

4.7-6.4

6.4-20.7

11

8

1.2

62

0.2

.1

.05

.2

'This layer includes 1.2 m of surface loam overlying 1.2 m of sand and gravel. 

2This layer includes a layer of sand and gravel extending from 1.7 to 3.0 m.

and these differences were of the same magnitude as 

wind-induced "noise" in the manometer readings. Nonetheless, results of analysis of these data probably are of the 

correct magnitude. Hydraulic conductivity values determined from the air-permeability data (table 6) range from 

about 6 to about 200 times larger than those determined 

inside the pit. Moreover, the values for the two less 

permeable layers do not agree at all well with the infiltration-based values of Prill and Aaronson (1973, p. 

740-741).

In addition to the field air permeability values determined 

for the two sites, air-permeability values were also determined in the laboratory on five drive-core samples obtained 

from the 3.7-5.2-meter depth at site 1. These samples, which 

were all from the lower gravelly sandy loam layer, had an 

average air permeability at their prevailing moisture content 

of about 0.3 darcies, equivalent to a hydraulic conductivity 

at 15.6°C of 0.24 m/day. This value is about five times 

greater than that determined from the field data for site 1, 

and about five times less than that for site 2. However, Prill 

(oral communication, 1971) found that laboratory hydraulic 

conductivities obtained on samples from the pit were much 

greater than the values suggested by the infiltration data. He 

hypothesized that many of the samples were disturbed 

during the drive-coring operation. If this were the case for 

the samples analyzed here, the results would tend to verify 

the field air-permeability values for site 1.

The differences in the air-permeability-based hydraulic 

conductivities for the two installations might be explained in 

several different ways. For example, it could be hypothesized that high residual moisture content of materials within the 

pit resulted in their having a low air permeability, whereas 

those outside the pit reflect the air permeability at field 

moisture capacity. For this explanation to be plausible, it 

must also be hypothesized that the surface loam and lower 

gravelly sandy loam are sensitive to the effects of wetting, 

with the permeability of these materials to air at field 

moisture capacity being 30 to 40 times greater than their
22 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

cr

z> 

o 

tr

UJ

O H

O

CO <

uj tr ic ho co

2 ~ Z> CO

? co

UJ

UJ ~Z. 

Q_

1.06

1.02

0.98

0.94

0.90

0.86

0.82

0.78

r^—r-i

~-X-^ X 

^ ^T^ A X

-NOTE:

~i— •*— .

-«-x^

A A *

Data for 

lower son 

gravelly s 

plotted, 

for the a 

plots are

•~ |

• i

X— V _

A A x^>

A ^

piezomete 

d and gra 

•andy loa 

They fall 

djacent Ic 

not mea

- — •

_* x _ x

. A **

jrs in the 

vel and tl 

m are not 

so near ( 

lyers that 

ningful

NUMERICAL 

SIMULATION

— -

• • ~-^-^

-*-x

AX

A

16 *

Jata

the ——

- — —

•

^-^

x

\. V

A ;

A*

i

PIEZOMETER 

AT BASE OF C

Lower gravelly 

sandy loam

Surface loam 

Atmosphere

^*^r>

Vx

k x-x '**

•^r~^

A**

)ATA

x

A

* ^"""""v^..

x-^x-^

« >

'*-,

— 1016

— 1014 

-1012

— 1010

:

i 

—1008

CO

cr

CD

o

UJ

X

o

<

s

UJ•z.

Q.

50 100 150 200 250 300 350 

TIME, IN MINUTES, SINCE 10 = 00 AM JULY I, 1971

400 450

FIGURE 11. Comparison of pneumatic heads measured July 1, 1971, at two depths at an installation inside Recharge Basin 306, Nassau

County, Long Island, N. Y., to the best-fitting numerical simulation.

intrinsic permeability as determined from analysis of infiltration of storm runoff water. Based on these hypotheses, 

the apparent agreement between air-permeability-based and 

infiltration-based intrinsic permeabilities for site 1 is fortuitous. Alternatively, it could be hypothesized that the 

different results for the two sites are due to areal heterogeneity of the materials with respect to air permeability. Still 

another explanation might be that the hydraulic conductivity of the materials in the basin has been reduced by 

suspended sediment and suspended or dissolved organic 

matter in the infiltrating storm runoff. In conclusion, 

however, it is not known which, if any, of these explanations 

is true, and the results of the Long Island test remain 

ambiguous. Additional tests are needed to evaluate the 

worth of the air-permeability method in humid regions, 

where high residual moisture content is difficult to avoid.

SUMMARY AND CONCLUSIONS

A method has been described for determining the vertical 

permeability to air (referred to as "air permeability") of 

layers of material in the unsaturated zone, based on

observations of the rate of change in pneumatic head at 

depth following a change in barometric pressure at land 

surface. These pneumatic head changes arise as the soil air 

compresses or expands in response to the surface pressure 

change, and the compression or expansion is slowed by the 

resistance to flow imposed by the medium.

Data collection and analysis for the method are relatively 

straightforward. Pneumatic head data may be obtained 

from nested piezometers in a single hole, and may be read 

fairly simply with an inclined manometer connected to the 

piezometers during a normal afternoon diurnal barometric 

change. Because the mathematical description of flow in 

layered media is complicated, data analysis is best accomplished by analog or by numerical simulation using the 

included digital computer program.

Hydraulic conductivity, rather than air permeability, is 

generally the parameter of interest to hydrologists. Air 

permeability values determined by the described method 

may be used to estimate the hydraulic conductivity of 

materials in the unsaturated zone under certain conditions. 

For example, if all the materials in the unsaturated zone are
SUMMARY AND CONCLUSIONS 23

at or below field moisture capacity, are permeable enough 

that the Klinkenberg effect is small, and are structurally 

unaffected by wetting, the method should provide reliable 

estimates of vertical hydraulic conductivity.

The method was tested at several field sites, including 

three that are described in detail in this report. The 

complexity of the field sites relative to the information 

desired varied among the three sites, and this complexity 

affected the relative success of the method.

For the first test described, the method was used to 

determine the vertical permeability to air of playa sediments 

overlying basalt in Idaho. This information was needed to 

determine the feasibility of temporarily storing radioactive 

gas generated in a planned "loss of fluids" test of a nuclear 

reactor. The values determined are a few times larger than 

laboratory values and a few times lower than values 

determined by analysis of data obtained by injecting air into 

auger holes. The results of this method probably are more 

accurate than those of the other two because of compaction 

of the samples during the coring process and because 

horizontal as well as vertical flow occurred during the air 

injection tests. The results of the test indicate that the 

method is a straightforward one for determining vertical air 

permeability, and generally could be used with confidence to 

determine that parameter.

For the second test described, the method was used to 

estimate the hydraulic conductivity of materials at an 

experimental artificial recharge site near Lubbock, Texas. 

This site is in a semiarid area, and soil moisture is at or below 

field capacity to a depth just above the water table. However, 

the formation includes several clay beds that contain 

significant quantities of expandable-lattice (swelling) clay 

minerals. Hence, parts of the formation might be sensitive to 

the effects of wetting on its intrinsic permeability.

Data from different installations and from the same 

installation on different dates are in general agreement, 

indicating repeatability of the tests. In addition, hydraulic 

conductivity values determined from the air-permeability 

data are in good agreement with values based on 

recharge-experiment data for the lower four of six levels 

tested, despite the presence of swelling clay minerals. In fact, 

a scheme to correct for the effects of wetting failed to 

improve the values for this site. Air-based values for the top 

two layers were both a few times larger than the hydraulic 

conductivity values based on recharge data for the same 

depth intervals. The results for the top layer, but not the 

second layer, appear to be general for the tests run in areas 

where the residual soil moisture content is low, and may be 

due to the effects of wetting on the structure of the 

near-surface materials.

The third test described was made to estimate the vertical 

hydraulic conductivity of materials at an experimental 

recharge facility on Long Island, N.Y., which has a humid 

climate. The potential exists at this site both for problems 

due to the effects of high residual moisture content in the soil

and to the effects on permeability of swelling clay contained 

in the materials in the unsaturated zone. Tests at two 

installations at this site, one within the spreading basin, and 

one outside but near the basin rim, showed substantially 

different results. Hydraulic conductivity estimates obtained 

from the installation inside the basin showed good agreement with hydraulic conductivity values determined from 

recharge data. However, analysis of data from the installation outside the basin resulted in hydraulic conductivity 

estimates from 6 to 200 times higher than those obtained 

inside the basin. The good agreement between air-based 

hydraulic conductivity values for the basin installation and 

recharge-based values coupled with the lack of agreement 

between air-permeability installations make the results of 

the test inconclusive. Additional tests will be necessary to 

assess the worth of the method to determine hydraulic 

conductivity in regions of high residual moisture content.

An added perspective on the worth of the air-permeability 

method may be obtained by considering the techniques 

currently available for determining hydraulic conductivity in 

the unsaturated zone. For most prediction problems, an 

average vertical hydraulic conductivity value is needed over 

some significant area, such as the areal extent of a recharge 

mound, or the area beneath a feedlot or sanitary landfill. 

However, current techniques, such as laboratory core 

analysis, infiltrometer techniques (Bouwer, 1961), or the dry 

auger hole method of Boersma (1965), all measure the 

hydraulic conductivity of relatively small samples of 

material. Consequently, a great many samples or tests would 

be necessary to approximate the average hydraulic conductivity of a heterogeneous layer or of a layer in which 

secondary permeability is predominant. Moreover, the 

hazard of compacting or disturbing the tested materials by 

current techniques is high, resulting in the possibility that the 

hydraulic conductivity determinations will be biased. Consequently, the air-permeability method offers promise in 

providing useful hydraulic conductivity data at relatively 

low cost, despite potential problems of high residual 

moisture content and of structural alteration of the medium 

upon wetting.

A final note of caution is needed in regard to determining 

the air permeability of relatively thin, highly permeable 

materials. Where the depth to water is not more than a few 

meters, and the materials in the unsaturated zone have a 

permeability of a few darcies or more, pressure differences 

would generally be too small to be measured by equipment 

described here. Although more sensitive measuring equipment could be used, short-term erratic atmospheric pressure 

variations that could be characterized as "noise" would be of 

similar magnitude, and the "signal," or desired pressure 

difference, would be difficult to determine.

In conclusion, a relatively simple and reliable field method 

to determine the vertical air permeability of layered 

materials in the unsaturated zone has been developed and 

tested. The method may also be used to estimate the
24 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

hydraulic conductivity of well-drained materials at depth in 

the unsaturated zone under certain conditions. Thus, the 

method is potentially useful for a number of hydrologic 

applications, such as evaluation of spreading sites for 

artificial recharge, of ground-water pollution hazards from 

feedlots or sanitary landfills, or of sites for temporary 

storage of gas in the unsaturated zone. However, applicability of the method to estimate hydraulic conductivity of 

materials at depth that may not be well drained must be 

determined through additional testing.

AIRK — A FORTRAN IV PROGRAM

FOR COMPUTING HYDRAULIC

CONDUCTIVITY FROM PNEUMATIC

HEAD READINGS

As described above, the AIRK program takes data 

acquired at the site, including screen depths and depth to the 

water table, estimates of air-filled porosity and preliminary 

estimates of PDe -based hydraulic conductivity, and periodic 

microbarograph and differential manometer readings and 

computes final estimates of air-permeability based hydraulic 

conductivity for each layer into which the unsaturated zone 

has been divided. These hydraulic conductivity values are 

determined by converting the preliminary estimates of 

hydraulic conductivity to air permeability, computing 

pneumatic heads by numerical solution of the equation for 

one-dimensional flow of gas through a porous medium 

(equation 3), and comparing the computed heads to those 

measured at the appropriate screen. Based on this comparison, a new trial value for air permeability is determined 

and new pneumatic head values computed. This process is 

repeated until the computed and measured heads are 

suitably close.

An important procedure in making the trial-and-error 

selections of air permeability is to perform the computations 

by adding one layer at a time, starting with the layer at 

greatest depth. By using this procedure, the only unknown to 

be solved for in a given series of simulations is the air 

permeability of the last layer added. This overall step wise 

procedure represents the basic function of the AIRK 

program. This procedure is discussed below, followed by a 

detailed discussion of some of the program elements that 

support its basic function. The following discussion is 

somewhat inexact in regard to the actual program sequence, 

but instead emphasizes the logic used in the program.

The stepwise procedure of numerically solving equation 3 

to provide hydraulic conductivity values for a three-layer 

system is illustrated in figure 12. At the start of the 

simulation, the unsaturated zone is divided into JMAX-1 

equal increments or nodes, and the faces of these nodes are 

numbered from 1 to JMAX, starting at land surface. The 

nodal position of each screen, assumed to represent the 

location of the layer boundary, is noted and stored as

JSD(NL) in statement 96. NL represents the screen number, 

as counted from the bottom, and land surface is designated 

as JSD (number of screens plus 1).

In the first sequence of simulations, a uniform trial air 

permeability, based on a user-supplied hydraulic conductivity estimate, a uniform air-filled porosity value, and an 

interpolated initial head value are assigned to each node 

from JSD(2)+1 to JMAX in statements 138-163 (see 

program listing). These nodes are used to simulate layer 1, as 

shown by nodal configuration A in figure 12. Pneumatic 

head values for each time step, as interpolated from those 

measured at screen 2, are assigned to node JSD(2). Based on 

these head values, the trial air permeability value, and the 

air-filled porosity, pneumatic head is computed for each 

time step at each underlying node. The head values 

computed for node JSD(l) are compared to those interpolated from the measurements for screen 1. Based on this 

comparison, a new larger or smaller trial air permeability 

value, as appropriate, is assigned to each node. The logic 

used in the program to select the new trial permeability value 

is described in the section on the automatic search 

procedure.

Once a new trial permeability value is selected, initial head 

conditions are reestablished, and head values at each node 

are again computed. Computed head at node JSD(l) is again 

compared to the measured head at screen 1, and a new trial 

air permeability selected. This process is repeated until the 

sum, for all the time increments, of the differences between 

computed and measured pneumatic head values has been 

minimized. The air permeability value that produces this 

minimum difference then remains assigned to each of the 

nodes from JSD(1)+1 to JMAX for the remainder of the 

analysis.

Once a satisfactory simulation of air movement in the 

bottom layer has been made, the user-supplied trial air 

permeability value for layer 2 is assigned to each node from 

JSD(3)+1 to JSD(2) (nodal configuration B in fig. 12) along 

with the air-filled porosity for that layer. Initial head values 

are assigned to each node from JSD(3)+1 to JMAX, and 

pneumatic head values interpolated from measurements for 

screen 3 for each time step are assigned to node JSD(3). 

Pneumatic heads are computed at each node for each time 

step, based on the trial air permeability value for layer 2 and 

the established air permeability for layer 1. Computed heads 

at node JSD(2) are compared to head values interpolated 

from those measured in screen 2, and a new trial air 

permeability value for layer 2 is determined from this 

comparison. The new trial air permeability is assigned to 

nodes JSD(3)+1 to JSD(2), and initial head values are 

reestablished for the entire system below node JSD(3). 

Heads are again computed, and those for node JSD(2) are 

compared to those for screen 2 to arrive at a new adjusted air 

permeability value. This whole process is repeated until an 

air permeability value is selected for layer 2 through the 

search procedure. At that point, the selected air permeability
AIRK A FORTRAN IV PROGRAM 25

Land surface

___ __ _ II «

Cement1^ <

Izi^zr1! i

Layers: ;

I^-I^-H : — — — _ « — — i — — — f

i i

I 1 '

I | ;

Cementzp !

1 1 •

LayerJ- 1

, I,':? 1 1 •:

1,1':

* . • . 9

. ; — ~~ * ——— * 

. * •

'Layer jji

••••:• .i

A _ _ _

i< ^^_m_ 

JilScreej) IK— ^—

Ql i i

jfl | |

>' 1 1

5 1 1

* 1 1 -; i i

15 1 1

| -^Screen

r-i • So • .

51 ' . •'. l~+~:

A* .' \ ——

1 : •'— •

• ' • • • ^-Screen

• Top of capil lary

. • . fringe • . • .'

A

A

-o- 

-o-

-o- 

-o-

-o-

_'V_

JSD(2)

JSD(I) 

JMAX

-o- 

-o- 

-o- 

-o-

-0-

y/

JSD(3) 

JSD(2) 

JMAX

"W^-

-O-

-o 

-o- 

-o-

-0-

S^v

/CY Ay%

JSD(4) = I

JSD(3)

EXPLANATION

SPECIFIED HEAD NODE

NODE FOR WHICH COMPUTED HEADS — 

are compared to measured heads

IMAGE NODE FOR--Establishing no-flow 

boundary at node JMAX

NODAL AREA FOR WHICH AIR 

PERMEABILITY —

H is being determined by 

>J successive approximations

^ K| has been determined and specified 

y> during first layer simulation

^i K2 has been determined and specified 

^ during second layer simulation

JMAX

Prototype

ABC 

Node configuration

FIGURE 12. Sketch showing the relationship between a hypothetical air-permeability installation at a site including three layers in the 

unsaturated zone and the nodal configurations used to solve for the air permeability of each layer. Nodal configurations A, B, and C 

represent those used to determine the permeability of layers 1, 2, and 3, respectively.

value remains assigned to each node from JSD(3)+1 to 

JSD(2) for the remainder of the analysis.

Once an air permeability value for the second layer is 

selected and assigned to the appropriate nodes, a new layer is 

added for analysis, as shown by nodal configuration C in 

figure 12. For the example shown, this is the uppermost 

layer. The user-supplied trial air permeability value for layer 

3 is assigned to nodes 2 to JSD(3), initial head values 

interpolated for the entire column are assigned to each node, 

and pneumatic head values interpolated from the 

microbarograph or altimeter record for each time step are 

assigned to node 1. Pneumatic heads are computed for each 

node, and those computed for node JSD(3) are compared to 

the measured values for screen 3. This comparison is, as 

above, used to repeatedly adjust the trial air permeability 

value until a satisfactory value is obtained.

At this point in the computer program, a complete 

simulation is made using the selected air permeability values 

for each layer, and pneumatic heads computed for each 

screen position during the simulation are printed out.

Comparison of the data for the final simulation to the 

measured pneumatic heads provides a means of assessing the 

overall adequacy of the analysis.

The final step of the program is to tabulate the hydraulic 

conductivity values computed from the selected air 

permeability values for each layer, along with appropriate 

information on the conditions of the test. (See "Sample 

output data.")

The above discussion describes the overall function and 

operation of the program. However, three of the procedures 

used in accomplishing this function are described in more 

detail below. These include the methods used to (1) 

synthesize head data from the supplied barometer and 

manometer readings, (2) numerically solve equation 3 for 

given values of air permeability and air-filled porosity, and 

(3) to select the most appropriate hydraulic conductivity 

value for each layer through the search procedure. Input 

data requirements are described in detail in a following 

section, as is the program output. Mnemonics used in the
26 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

program are tabulated, the program is listed, and sample 

output data are presented.

COMPUTATION OF PNEUMATIC HEAD VALUES

For the method as described, pneumatic head values are 

read sequentially at somewhat irregular time intervals at a 

few points in the unsaturated zone profile. However, the 

computer program requires pneumatic head values at 

simultaneous short equal time intervals, and also requires, 

for initial conditions, that the pneumatic head be described 

at each node at the start of the test. Consequently, the 

barometric and differential manometer readings must be 

interpolated in time to produce simultaneous values of 

pneumatic head, and the initial readings must be interpolated in space. The first step in establishing a table 

defining pneumatic head versus time at each screen is that of 

interpolating pneumatic head at land surface from the 

microbarograph or altimeter records. This table is formulated by statements 60-72 by linear interpolation of 

supplied values of microbarograph readings versus time, as 

read in card series 9. If altimeter readings were obtained, 

rather than microbarograph records, they must be converted 

by the user to equivalent microbarograph readings before 

they are read into the program.

Pneumatic head values for each screen are determined in 

program steps 73-94 by converting the differential 

manometer readings, read in inches of manometer fluid in 

card series 10, to inches of mercury; and linearly interpolating between the converted readings to obtain values 

coincident with the end of each time step. These interpolated 

manometer readings are added algebraically to the interpolated barometric reading to give values for pneumatic 

head at the end of each time increment for each screen.

Once the differential-manometer readings have been 

converted to pneumatic head, initial head values for each 

space node (required as an initial condition for the numerical 

solution of equation 3) are established by program 

statements 95-116. For each layer except the bottom one, 

these values are determined by linear interpolation of the 

heads at time zero for the screens that bracket the layer. 

However, with regard to the bottom layer, linear extrapolation of the head gradient between screens 1 and 2 (fig. 1) 

would result in an implied violation of the no-flow boundary 

condition at the top of the capillary fringe. Hence, within the 

bottom layer, head values are computed from a sine curve fit 

through the pneumatic head values for screens 1 and 2. In 

addition to meeting the boundary condition at the water 

table, additional theoretical justification exists for assuming 

that the data for this layer should approximate a sine curve. 

Carslaw and Jaeger (1959, p. 99) show that, for four different 

sets of initial conditions for one-dimensional heat flow in a 

slab, the temperature distribution in the slab converges 

toward a sine-curve configuration with time. Based on the 

analogy between heat and fluid flow and the fact that the 

observed initial pneumatic head differences generally result

from atmospheric pressure changes that either occurred or 

began some time before readings were begun, the sine-curve 

assumption seems well justified.

NUMERICAL SOLUTION OF EQUATION 3

For the computer analysis, flow of air in the unsaturated 

zone, as described by equation 3, is analyzed by a numerical 

finite-difference technique. Use of the technique requires 

that the space and time derivatives of pneumatic head be 

approximated by values of the difference quotients of head 

at separate discrete points both in space and in time. The 

discrete points in space are generally termed nodes, and for 

this application it was convenient to express the 

head-difference quotients for equally spaced nodes solved at 

equal time increments. These head-difference quotients may 

be expressed in one of several ways, as described in 

textbooks on numerical analysis. For the AIRK program, 

the head-difference quotients are formulated in statements 

127-167 as a set of backward difference equations (von 

Rosenberg, 1969, p. 20-22), one for each node at a given time 

increment. These formulations, with the inclusion of the 

appropriate boundary and initial conditions, result in a set of 

equations whose nonzero coefficients form a tridiagonal 

matrix that is efficiently solved by use of the Thomas 

algorithm (von Rosenberg, 1969, p. 113). The Thomas 

algorithm is programed as the TRIMAT subroutine in the 

AIRK program.

Assignment of air permeability and air-filled porosity 

values, used in computing coefficients in the finite-difference 

equations, is straightforward except for the nodes representing the layer boundaries. Because a face-centered node 

configuration is used in this analysis, the interblock air-filled 

porosity at each layer boundary is computed as the 

arithmetic average of the values for the two adjoining layers. 

No need arises to compute interblock air permeability for 

this node configuration, however.

THE AUTOMATIC SEARCH PROCEDURE

The automatic search procedure is based on adjusting air 

permeability for a given layer to minimize the sum, for all the 

time increments, of the difference between measured and 

computed pneumatic heads in either the bottom screen, for 

the first layer simulation; or in the screen at the bottom of the 

uppermost layer, in subsequent simulations. As an example, 

if barometric pressure declined during the test, pneumatic 

heads would increase with depth. Moreover, if the trial air 

permeability is larger than the true value, the computed head 

differences across the layer will be less, and the average 

computed pneumatic heads smaller, than those measured. 

Hence, the sum for all of the time steps of the difference 

between computed and measured heads, symbolized HDIF 

in the program, would be negative. Conversely, if barometric 

pressure were rising, the head at depth would be less than 

that nearer the surface. Under these conditions, HDIF 

would be positive if the permeability estimate was too low.
AIRK A FORTRAN IV PROGRAM 27

Hence, for tests involving rising barometric pressure, HDIF 

is multiplied by an IBSN value of -1 to make its sign 

consistent with that for falling barometric pressure.

As indicated by the above discussion, the need for 

increasing or decreasing the trial air permeability estimate 

may be determined by the sign of HDIF. In the program, 

equation 3 is solved for the user-supplied trial value of air 

permeability by the TRIMAT subroutine, HDIF is computed in statements 176-181, and its sign determined by 

statement 186. If HDIF is negative, the factor representing 

air permeability is doubled, and if positive, halved, based on 

statements 190 or 195. In succeeding computations, the sign 

of HDIF is determined for the new air permeability value. If 

the sign of HDIF is the same as that determined in the 

previous step, the air permeability is again halved or doubled 

as before and the process repeated.

A change in the sign of HDIF indicates that the true air 

permeability should lie between the current and previous 

trial values. Hence, once a change in sign of HDIF occurs, an 

air permeability value is interpolated from the current and 

previous HDIF and air permeability values, based on 

statements 200 or 204. Following the interpolation of the 

head difference data, the various hydraulic conductivities 

equivalent to the trial air permeabilities are listed, and the 

computed heads for the node representing the screen for 

which comparisons were made are tabulated for selected 

time increments in a column beneath the respective 

hydraulic conductivity values.

At this point in the program, the interpolated air 

permeability value is used to again compute pneumatic 

heads, and HDIF and its sign are again determined. This 

time, however, the air permeability estimate is adjusted by a 

factor of only 10 percent. Computations are repeated until 

HDIF changes sign from that determined using the interpolated air permeability value. At this time, a second 

interpolation is made to determine hydraulic conductivity, 

and head values for the appropriate screen are computed a 

final time. Each of the hydraulic conductivity values used for 

the second set of computations and the resulting pneumatic 

heads are then tabulated.

Some situations may occur for which the search procedure 

will fail to find an air permeability value that minimizes 

HDIF. For example, head differences across a given layer 

might be so slight that, due to short term random variability 

in the readings, the measured heads in the lower screen may 

show less head difference than do those measured for the 

overlying screen. If, for this or some other reason, HDIF for 

a given layer fails to change sign after 10 trial values of air 

permeability have been used to solve equation 3, the 

program notes its failure to find a solution. The program 

then lists the 10 values of hydraulic conductivity, tabulates 

the computed head values at the node representing the screen 

at the base of the layer for each permeability estimate, and 

terminates.

PROGRAM VERIFICATION

The AIRK program was verified by computing pneumatic 

heads for selected depths and assumed values of permeability and air-filled porosity based on two different analytical 

equations. These computed heads were entered as data into 

the program, along with wrong "estimates" of the equivalent 

hydraulic conductivity and correct values for the air-filled 

porosity. The program was run to compute values of 

hydraulic conductivity, and these values compared to those 

used in the analytical equations for the pneumatic head 

calculations.

For the first verification procedure, head values were 

computed at four levels in a uniform layer of unsaturated 

material, based on an analogous heat-flow equation 

presented by Carslaw and Jaeger (1959, p. 104, eq. 4). This 

equation is based on the assumptions that pneumatic head at 

land surface increases linearly as a function of time, that no 

flow occurs at the water table, and that pneumatic head is 

initially constant throughout the profile. These computed 

heads were inserted into the program, along with different 

wrong values for the hydraulic conductivity of each of the 

four assumed layers. Computed equivalent hydraulic conductivity values for each layer, computed through the 

automatic search procedure, ranged from 2.99 m/d to 3.20 

m/d, compared to the value of 3.05 m/d used to compute the 

head values. A second run using other wrong estimates 

indicate that the search procedure converged to the same 

estimates whether the initial estimates were too high or too 

low.

A second verification run was made by computing heads 

at four levels in a two-layer unsaturated zone, also based on 

an analytical equation given by Carslaw and Jaeger (1959, p. 

323-324). For this equation, pneumatic head at the land 

surface boundary is assumed to undergo an instantaneous 

change at time /=0, and to maintain that new value 

thereafter. The other boundary, corresponding to the water 

table, is assumed to be maintained at the initial pneumatic 

head. Although this boundary condition is different than 

that existing in the field or in the program, boundary 

conditions specified in the program were modified during 

that run to match those assumed for the analytical equation. 

Head values computed from the equation were used in 

conjunction with wrong values of air permeability in the 

AIRK program to estimate the true air permeability of four 

layers, including one that coincided with the lower layer in 

the analytical solution, and three that divide the upper layer. 

Computed values for the equivalent hydraulic conductivity 

are compared to those used in the analytical equation in the 

table below.

Depth

interval,

meters

0-15

15-35

35-50

50-100

Hydraulic Conductivity, m/day

Initial

"estimate"

50.

3.

1.

10.

Final True

computed value value

7.76 8.

8.12 8.

8.55 8.

.970 1.

Air-filled

porosity

0.2

.2

.2

.1
28 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

As may be seen from the table, the computed values range 

from -3 percent to +7 percent of those used in the analytical 

equation.

Results of the two comparisons show that the program 

yields reasonably good approximations of the true hydraulic 

conductivity. Errors arise in part because of truncation and 

roundoff errors in computing heads from the 

finite-difference approximation of the differential equation, 

so that they do not agree exactly with those computed 

analytically. The significance of these errors is magnified in 

the permeability or hydraulic conductivity estimates, 

because that parameter must be changed by a few percent in 

order to alter computed head by a few tenths of a percent. 

This lack of sensitivity of head to small changes in 

permeability or hydraulic conductivity contributes to the 

general uncertainty of any problem concerned with the 

identification of aquifer properties. (See, for example, Hefez 

and others, 1975.) Moreover, as noted above, many other 

sources of error exist in estimating hydraulic conductivity 

from the air permeability at the prevailing moisture content. 

Hence, it was concluded from the above comparisons that 

the program is adequate for its intended purpose, although 

additional work is desirable to improve the comparison with 

the analytical equations.

INPUT DATA REQUIREMENTS

Required data to execute the program include ten cards or 

card series described in sequence and detail below.

Card 1.—Contains in the first 32 columns, under the label 

SINAM, the name of the site at which data were collected.

Card 2.—Contains in columns 1-4 and 7-10, the beginning 

and ending 24-hour clock times, labeled CLCKB and 

CLCKE, for the data being analyzed. In columns 13-24, the 

date the data were collected is given under the label RD ATE. 

Example:

0845bbl615bbMARbl5,bl972

Card 3.—Contains, right-justified in columns 1-2, a 

number labeled IBSN that is used in the computations to 

account for the main direction of the barometer change. 

Because of limitations of the search procedure, the 

pneumatic head data must follow either a generally upward 

or a generally downward trend, although minor deviations 

are permissible. The period of record to be analyzed should 

be selected accordingly. If the barometer fell during most of 

the test, bl or +1 should be entered; if it rose during most of 

the test, -1 should be entered.

On the same card, an option code to weight the data, 

labeled IWTO, may be entered in column 4. If it is felt that 

the later data are more reliable than the early data, and hence 

worthy of greater weight, a "1" should be entered in column 

4. If, on the other hand, the early data are to be weighted, a 

"2" should be entered in column 4. If column 4 is left blank, 

all data are weighted equally.

Card 4. This card contains five numbers. Rightjustified in columns 1-2 is an integer labeled NLYS that 

designates the number of layers used to represent the 

unsaturated zone. This number must be equal to the number 

of screens for which data are to be read in. As dimensioned, 

this number would be limited to 9.

In columns 3-9, a floating-point number labeled UZD is 

entered to designate the thickness of the unsaturated zone, in 

feet. In general, it is convenient to round this number to the 

nearest foot.

In columns 10-13, a right justified integer labeled JM AX is 

entered to designate the total number of nodes to be used to 

simulate the unsaturated zone. It has been convenient for 

this study to divide the unsaturated zone into 1-foot 

increments. Thus, JMAX generally is set equal to 1 plus the 

depth to the water table in feet. As dimensioned, this number 

has a maximum value of 200.

In columns 14-20, a floating-point number labeled TINC 

is entered to designate the length, in minutes, of each time 

step to be used in the computations. For a typical afternoon 

diurnal record, this number would be set equal to 5.

In columns 21-24, a right-justified integer labeled NTS is 

entered to designate the total number of time steps to be used 

in the analysis. This number is equal to one plus the period of 

record in minutes divided by the time increment, which also 

is in minutes. As dimensioned, this number would be limited 

to 150.

Card 5.—This card contains, right-justified in columns 

1-2, an integer labeled NSTEP designating the number of 

time increments after which head data should be printed out. 

This number allows head data to be printed out less 

frequently than the computations are made, providing more 

succinct tables of head versus time for examination or 

plotting. For the usual observations of afternoon diurnal 

fluctuations made during this study, head values were 

computed at 5-minute intervals and printed out at 15-minute 

intervals. NSTEP would be set equal to 3 under these 

conditions.

Card 6.—This card contains two floating-point numbers. 

The first number, labeled SGMF, is punched in columns 1-7 

and designates the specific gravity of the manometer fluid. 

The second number, labeled PBAR,, is entered in columns 

8-14 to designate the mean station barometric pressure, in 

inches of mercury, during the period of record being 

analyzed. It should be kept in mind that the station pressure 

is the actual atmospheric pressure at the site, rather than the 

pressure corrected to sea level generally quoted by the U.S. 

Weather Service. If such data are not available, the value 

may be estimated by correcting the average sea-level 

pressure of 29.92 inches of mercury for the altitude of the 

site, using, for example, the equation:

where

P=29.9 exp (-altitude in feet/27,000.),

P is in inches of mercury.
A1RK A FORTRAN IV PROGRAM 29

Alternatively, the station pressure could be determined 

from handbook values of atmospheric pressure versus 

altitude. It is adequate to enter this value to the nearest 0.1 

inch.

Card 7.—Card 7 is the first of NLYS cards containing 

data on screen depths and assumed or estimated physical 

properties of each layer used to represent the unsaturated 

zone. It contains, in the first 5 columns, a floating-point 

value, labeled SD, representing the depth to the bottom 

screen, from land surface, in feet. An initial estimate of the 

hydraulic conductivity, HYDC, of the lowermost layer, in 

feet per day, is entered as a floating-point number in columns 

6-13. The following table may be used to make such 

estimates for each layer, assuming the record was obtained 

during a normal afternoon barometric pressure decline.

A/! 1 Hydraulic conductivity,

feet per day

.2:5 

.05:2 

.01:05

< .01

0.01

.1 

1.

10. 

100.

'ah is the approximate maximum pneumatic-head difference, in inches of water, measured between 

screens straddling the layer for which estimates are to be made.

The air-filled porosity estimate, labeled DPOR, based on 

neutron logs, core analyses, and any other available data for 

the lowermost layer is entered as a two-digit decimal fraction 

in floating-point format in columns 14-17. DPOR generally 

would be equivalent to specific yield.

If it is desired for any reason to retain the provided 

hyraulic conductivity estimate, a "1" (labeled KOP) should 

be entered in column 20. As an example, a plot of computed 

and measured head data for a previous computer run may 

have indicated a less-than-ideal curve fit for a given layer. 

Under these conditions, the investigator may choose, from 

the plot, a somewhat different value for hydraulic conductivity than was selected by the programed search routine. 

The investigator should then insert his selected value in the 

appropriate space on the card, enter a 1 in column 20 of the 

card representing that layer, and rerun the program to 

obtain improved estimates of the hydraulic conductivity of 

the layers overlying that for which the hydraulic conductivity estimate was revised.

Each additional card of the NLYS cards in the Card 7 

series must contain the screen depth, hydraulic conductivity, 

and air-filled porosity or specific yield for the overlying 

screens and layers. These cards must be in sequence, it being 

remembered that both the screens and layers are numbered 

upward from the water table, but that all depths are 

measured from land surface.

Card 8.—This card contains either a 1 in column 1 

(labeled IVIS) or a blank in column 1 and two numbers in 

double-precision exponent format in columns 2-11 and 

12-21 designating the dynamic viscosity of air, labeled ADV, 

at the prevailing temperature of the soil gas at the site and the 

kinematic viscosity of water, labeled WKV, at a temperature

selected by the investigator. ADV must be in pound-seconds 

per square foot, and WKV in square feet per second. If 1 is 

entered in column 1, ADV and WKV at 60°F are read in by 

default.

Card 9 (series). This series contains data on barometric 

pressure changes during the period of record to be analyzed. 

For the first card, the barometric pressure, BAR1, as 

measured in inches of mercury, preferably precise to 0.001 

inch, is entered in columns 1-6. In general, it is preferable to 

drop the digits to the left of the decimal point, adjusting the 

readings by 0.5 inch or so if necessary to avoid any change in 

dropped digits during the record. Dropping these digits 

reduces round-off error in the finite difference 'computations. The time after the start of the record in minutes is 

entered as a floating point numeral in columns 7-11. The first 

such value, designated Tl, would be 0.

Each subsequent card in the Card 9 series should contain, 

in sequence, the barometer reading designated BAR2 and 

the time, designated T2, in minutes, after the start of the test 

in columns 1-6 and 7-11, respectively. The time readings 

need not be incremented equally. However, the last time 

value must equal the product TINC X (NTS-1).

A card containing 99999. in columns 1 -6 should be entered 

in the data deck to indicate the end of the barometer-time 

file.

Card 10 (series).—This card series contains the differential manometer readings for each screen and time during the 

period of record. As many as six differential manometer 

readings and their corresponding times are entered on each 

card, with readings for the bottom screen in the leftmost 

columns. The .entered differential manometer readings 

represent the net difference, in inches of manometer fluid, 

between the reading when the manometer is opened to a 

given screen and that obtained when both ports of the 

manometer are open to the atmosphere. Readings are 

considered positive when the downhole pneumatic head is 

greater than that at land surface. For each card, initial 

differential manometer readings, designated DMR1 or 

DMR2, for screens 1, 2, 3, ..., are entered as floating point 

numbers in columns 1-6, 12-17, 23-28,..., and the time value 

for each reading in minutes since the test began and 

designated TRl(NL) or TR2(NL), is entered in 

floating-point format in columns 7-11, 18-22, 29-33,.... For 

the first card, all time readings must be 0., and for the last 

card, all time readings must equal TINC X(NTS-l). For each 

of the other cards, the times for readings at each screen may 

differ from the others, and need not be incremented equally. 

Readings for any given screen must be in sequence, however, 

and each card must contain a reading for each screen. 

Consequently, if, during a sequence of measurements, a 

reading for a given screen is not available, a value should be 

interpolated for inclusion on the card.

The last card in the Card 10 series should be followed by a 

card containing 99999. in the first 6 columns to indicate the
30 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

end of file. This card will complete the data deck for the 

program.

LIST OF MNEMONICS AND CONSTANTS 

USED IN THE AIRK PROGRAM

The following list describes the mnemonics and constants 

used in the AIRK program that are not described in the 

section on input data. Items are listed in the order that they 

appear in the program, and, where appropriate, units are 

provided.

NL is a counter designating screen or layer number. It is 1 

for the bottom screen or layer and NLYS for the top screen 

or layer.

MPT designates whether the initial or final hydraulic 

conductivity estimates are to be printed, and determines 

whether the heading for table 1 or table 5 will be printed, as 

appropriate.

HYDS (NL)=final hydraulic conductivity selected 

through the search procedure for the NLth layer, saved for 

printing in table 5.

NBS=total number of observation points, including that 

at land surface.

NT=counter for time steps.

Ml, M2=counters for interpolating barometer and 

differential manometer readings at equal time steps.

13.56=specific gravity of mercury at 60°F.

HR (NBS, NT)=atmospheric pressure at land surface at 

the end of the NTth time step, in inches of mercury, as 

interpolated from input data.

NTL (NL)=counter for keeping track of last previous 

interpolated head value for the NLth screen.

DV=interpolated differential manometer reading, converted to inches of mercury.

HR (NL, NT)=pneumatic head reading at the NLth screen 

and NTth time step, in inches of mercury, as determined by 

interpolation of input data.

JSD is the node number for the NLth screen, or for land 

surface (NL=NBS), as measured from land surface. Land 

surface is at node 1 (i.e., JSD(NBS)=1).

NNL (NL)=number of nodes in NLth layer.

NBOT=number of nodes between screen 1 and screen 2.

ELB=floating-point designation of number of nodes 

representing bottom layer.

XL=ratio of number of nodes between screens 1 and 2 to 

number of nodes in layer 1.

HMAX=pneumatic head at water table (node JMAX).

Jl=counter to vary the upper boundary node designation 

within nested DO-loops as layers are added to the model.

HNS (J)=initial pneumatic head, in inches of mercury, at 

any node J, used to specify initial conditions in the solution 

of equation 3.

J2=counter to designate the bottom node of the layer for 

which trial values of permeability are being established. For 

this program each layer above the first is assumed to extend

to the bottom of the node immediately above the subjacent 

screen.

DELH= difference in pneumatic head at the start of the 

test across any layer except the bottom one.

TIME=elapsed time from start of test, in minutes.

DELT=time increment in days.

DZSQ=square of the depth increment represented by 

adjacent nodes, in square feet.

P (NL)= a factor representing air permeability, in ft5 per 

inch-pound-day. These awkward units arise from expressing 

hydraulic conductivity in ft per day, pneumatic head in 

inches of mercury, dynamic viscosity in pound-seconds per 

ft2 , kinematic viscosity in ft2 per second, and gas specific 

storage (SS(J)) in ft3 per inch-pound.

386.4= acceleration due to gravity, inches per second 

squared.

SS (J)=specific storage term for soil gas for the material 

represented by the Jth node, in ft3 per inch-pound.

846.=weight density of mercury at 60° F, in pounds per ft3 .

R (J)=coefficient for computing time rate of change of 

storage within the material represented by the Jth node due 

to a change in head during the time period DELT.

JJ2, JM2, JPl=delimiters to control the head difference 

quotients at nodes not on either boundary of the system 

being analyzed.

ICT2=a counter to designate whether air permeability is to 

be halved or doubled, or to be increased or decreased by 10 

percent.

JS=node designation of screen within or at the bottom of 

the layer being simulated, used in comparing measured and 

computed heads at that point.

ICT=counter to designate the number of trial air 

permeability values used for a given layer.

A(J), B(J), and C(J)=the finite-difference equation coefficients by which the unknown heads for the current time 

step at nodes J-l, J, and J+l are multiplied in the 

backward-difference formulation of equation 3.

H (J)=head at node J at time step NT, initialized to 

HNS(J) for time step 1.

D (J)=finite-difference equation coefficient times the head 

at the Jth node computed during the previous time step.

TRIMAT=subroutine for solving the set of implicit 

finite-difference equations for pneumatic head at the current 

time step by use of the Thomas algorithm.

HSCS (ICT, NT)=pneumatic head expressed in inches of 

mercury, computed for the ICTth trial hydraulic conductivity and NTth time step for the NLth screen. These head values 

are saved for printing in table 3.

HSCS (INL, NT)=computed head for the NLth screen at 

time NT, based on the final simulation. These head values 

are printed in table 4.

HDIF=the difference between computed and measured 

heads, summed for each time step, for the screen within or at 

the base of the NLth layer.
A1RK A FORTRAN IV PROGRAM 31

ISIGN=counter to keep track of the sign of HDIF in 

repeated simulation trials.

PS=air permeability for the current time step saved for 

possible interpolation of an air permeability value following 

the next simulation trial.

HDS=HDIF for the current trial saved for possible 

interpolation of an air permeability value following the next 

trial.

ISIGN2=counter for sign of HDIF for the current trial air 

permeability, used for comparison with ISIGN, the sign of 

HDIF for the previous trial.

AIRK PROGRAM LISTING

FORTRAN IV G LEVEL 21 MAIN DATE = 77318 18/44/45

0001

0002

0003

0004

0005

0006

0007

0008

0009

0010

0011

0012

0013

AIRK PROGRAM A 

THIS PROGRAM COMPUTES THE PNEUMATIC OlFFUSIVITY OF N LAYERS OF UNSAT A 

URATED MATERIAL, BASED ON READINGS OF PRESSURE VERSUS DEPTH AT EACH A 

LAYER BOUNDARY, AS MEASURED DURING A CHANGE OF ATMOSPHERIC PRESSURE. A 

REQUIRED DATA ARE BAROMETRIC READINGS VERSUS TIME, PRESSURE CHANGES A 

VERSUS TIME FOR EACH DEPTH, THE NUMBER OF LAYERS, THE THICKNESS OF A 

THE UNSATURATED ZONE, THE DEPTH OF EACH SCREEN, MEAN STATIGM PRESSURE A 

DURING THE TEST, SPECIFIC GRAVITY OF THE MANOMETER FLUID, AND INITIAL A 

ESTIMATES OF HYDRAULIC CONDUCTIVITY AND DRAINED POROSITY FOR EACH A 

LAYER, STARTING FROM THE HATER TABLE. A

IMPLICIT REAL*8(A-H,0-Z» . A

COMMON /COMI/ HI200),A(200),B<200|,C(200),0(200) A

DIMENSION HR(9,150i, SSI200), R(200), SDC9), HYDS(9), DPORC9), NTL A 

1(9), NNL(9), DMRK9I, DMR2(9), TRK9), TR2(9), JSD(9), P(9), KOP(9 

2), XI2), HYDC(9,10)

DIMENSION HSCSdO, 200) , HNS(200)

DIMENSION SINAMC8), RDATEC3), NLINK9)

DATA X/l.O,0.l/

X(ICT2)=FACTOR FOR VARYING PERMEABILITY BY 100 PER CENT FOR THE FIRST 

SWEEP, AND BY 10 PER CENT FOR THE SECOND SWEEP.

READ 82, (SINAM(I),1=1,81

SINAME IS A 32- CHARACTER FIELD FOR THE SITE NAME AND LOCATION.

**************************************:

READ 78, CLCKB,CLCKE,CRDATE(II,I=1,3)

t*************************************:

RDATE IS A 12-CHARACTER FIELD FOR THE DATE READINGS WERE OBTAINED. 

CLCKB AND CLCKE ARE 4-CHARACTER FIELDS FOR BEGINNING AND ENJING 24- 

HOUR CLOCK TIMES.

*************** 

READ 69, IBSN,IWTO

ISSN IS SIGN OF BAROMETER CHANGE, -1 FOR RISING BAROMETER; +1 FOR 

A FALLING BAROMETER.

IWTO IS A CODE FOR WEIGHTING THE DATA WITH RESPECT TO TIME, 0 FOR 

EQUAL WEIGHTING, 1 FOR MORE HEAVILY WEIGHTING EARLY DATA, AND 2 FOR 

MORE HEAVILY WEIGHTING LATER DATA.

1

2

3

4

5

6

7

8

9

10

11

12

13

READ 63, NLYS,UZD, JMAX,TINC,NTS 

****************

READ 69, NSTEP 

N=NSTEP+1

NSTEP DESIGNATES THE INTERVAL AT WHICH HEAD READINGS ARE PRINTED.

[*****************]

NLYS = NUMBER OF LAYERS. A

NTS=TOTAL NUMBER OF TIME STEPS. A

JMAX= TOTAT NUMBER OF POSITION NODES. A 51

TINC=TIME INCREMENT FOR PNEUMATIC HEAD COMPUTATIONS, IN MINUTES. A 52

UZO=UNSATURATED ZONE DEPTH, IN FEET. A 53

********************************************************************* A 54

READ 64, SGMF,PBAR A 55

SGM*SPECIFIC GRAVITY OF MANOMETER FLUID. A 57 

PBAR = MEAN PRESSURE DURING TEST, IN INCHES OF MERCURY. A 58

A 14 

A 15 

A 16 

A 17 

A 18 

A 19 

A 20 

A 21 

A 22 

A 23 

A 24 

A 25 

A 26 

A 27 

A 28 

A 29 

A 30 

A 31 

A 32 

A 33

34

35

36

37

38

39

40

A 41 

A 42 

A 43 

A 44 

A 45 

A 46 

A 47 

A 48

49

50
32 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

0014

0015

0016

0017

0018

0019

0020

0021

0022

0023

0024

0025

0026

0027

0028

0029

0030

0031

0032

0033

0034

0035

0036

0037

0038

0039

0040

0041

0042

0043

0044

0045

0046

0047

0048

0049

0050

0051

0052

0053

0054

0055

0056

0057

0058

0059

DO 1 NL=1,NLYS

1 READ 65, SD(NL),HYDC(NL,1I,DPOR(NL),KOPCNL)

Q **********************************************************************

C SD(NL)=DEPTH OF NLTH SCREEN, NUMBERED FROM THE BOTTOM AND MEASURED

C FROM LAND SURFACE.

C DPOR(NL)=DRAINED POROSITY ESTIMATE FOR NLTH LAYER,

C (HYDCCNL,1)=ESTIMATED.HYDRAULIC CONDUCTIVITY, IN FEET PER DAY, FOR

C LAYER 1. THIS VALUE IS CONVERTED TO AIR PERMEABILITY, PINL), WITHIN

C THE PROGRAM.

C KOP=1 WILL CAUSE PROGRAM TO RETAIN THE HYDC VALUE FOR A GIVEN LAYER,

C RATHER THAN TO CONTINUE SEARCHING FOR A BEST-FIT SOLUTION.

C

c **********************************************************************

READ 68, IVIS,ADV,WKV

IVIS IS A CODE TO SPECIFY WHETHER DEFAULT VALUES OF AIR AND WATER 

VISCOSITY ARE TO BE USED. IF IVIS=1, THE VALUES FOR ADV AfiD WKV AT 

60 DEGREES F ARE READ IN AUTOMATICALLY. IF THE USER WISHES TO SUPPLY 

VALUES FOR DIFFERENT TEMPERATURES, IVIS SHOULD BE A BLANK OR ZERO. 

ADV=DYNAMIC VISCOSITY OF AIR AT THE PREVAILING SOIL TEMPERATURE, 

IN LB*SEC/SQUARE FOOT.

WKV=KINEMATIC VISCOSITY OF WATER AT DESIRED TEMPERATURE, IN SQUARE 

FEET PER SECOND.

IF ( IVIS.NE.l) GO TO 2

ADV=3.74D-07

WKV=l.21D-05

2 CONTINUE

THE FOLLOWING STATEMENTS PRINT OUT INPUT DATA AND INFORMATION 

REGARDING THE TEST. 

MPT=0

3 PRINT 76, (SINAMU),I = 1,8) 

PRINT 72, UZD 

PRINT 73, NLYS,JMAX

PRINT 77, CLCKB,CLCKE, (RDATEU ), 1 = 1,3) 

PRINT 74, NTS.TINC 

PRINT 97, PBAR 

PRINT 98, SGMF 

IF (IBSN.EQ.1) GO TO 4 

PRINT 95 

GO TO 5

4 PRINT 96

5 CONTINUE

IF *IWTO-1) 6,7,8

6 PRINT 92 

GO TO 9

7 PRINT 93 

GO TO 9

8 PRINT 94

9 CONTINUE

PRINT 101, ADV,WKV

IF (IVIS.EQ.l) PRINT 102

PRINT 84

IF IMPT.EQ.U GO TC 10

PRINT 79

GO TO 12

10 PRINT 83

DO 11 NL=1,NLYS

11 HYDCINL,1)=HYDS(NLI

12 PRINT 80 

NBS=NLYS+1 

SD(NBS)=0. 

NL=1

PRINT 75, NL,SO<2),UZD,HYDCINL,1),OPOR(NL) 

IF (NLYS.EQ.1) GO TO 14 

DO 13 NL=2,NLYS

13 PRINT 75, NL,SD(NL+1),SD(NLI,HYDCINL,1),DPORINL)

14 IF IMPT.EQ.l) GO TO 62

PRINT 99, SDC1)

THE FOLLOWING STATEMENTS COMPUTE ATMOSPHERIC PRESSURE AT LAND SURFACE FOR EQUAL TIME INCREMENTS BY INTERPOLATION FROM USER-SUPPLIED

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99 

A 100 

A 101 

A 102 

A 103 

A 104 

A 105 

A 106 

A 107 

A 108 

A 109 

A 110 

A 111 

A 112 

A 113 

A 114 

A 115 

A 116 

A 117 

A 118 

A 119 

A 120 

A 121 

A 122 

A 123 

A 124 

A 125 

A 126 

A 127 

A 128 

A 129 

A 130
AIRK A FORTRAN IV PROGRAM 33

0060

0061

0062

0063

0064

0065

0066

0067

0068

0069

0070

0071

0072

0073

0074

0075

0076

0077

0078

0079

0080

0081

0082

0083

0084

0085

0086

0087

0088

0089

0090

0091

0092

0093

0094

0095

0096

0097

0098

0099

0100

0101

0102

0103

0104

0105

0106

0107

0108

0109

C BAROMETER READINGS. A 131

READ 66, BAR1.T1 A 132

C BAR!, T1=RAW MICROBARCGRAPH READINGS AT TIME Ti,-WITH BAR1 IN A 133

C INCHES OF MERCURY, AND Tl IN MINUTES SINCE START OF TEST. A 134

HRINBS, i)=BARl A 135

M2=0 A 136

15 READ 66, BAR2,T2 A 137

IF (BAR2.GT.99.90I GO TO 17 A 138

Ml=M2+l A 139

M2=(T2/TINC)+1 A 140

DO 16 NT=M1,M2 A 141

16 HR<NBS,NTI = BARl + (BAR2-BARl)*((NT-l)*TINC-Tl)/<T2-TU A 142

Ti=T2 A 143

BAR1=8AR2 A 144

GO TO 15 A 145

17 CONTINUE A 146

C THE FOLLOWING STATEMENTS RELATE TO COMPUTING PNEUMATIC HEAD AT EACH A 147

C SCREEN FOR THE END OF EACH TIME INCREMENT. THESE COMPUTATIONS A 148

C ARE ACCOMPLISHED BY CCNVERTING THE MANOMETER READING FOR A GIVEN A 149

C SCREEN AND TIME TO INCHES OF MERCURY, INTERPOLATING TO OBTAIN A A 150

C A VALUE FOR THE READING AT THE END OF THE TIME STEP, AND ALiiE- A 151

C BRAICALLY ADDING THE READING TO THE CONCURRENT BAROMETER READING. A 152

NT=1 A 153

READ 67, (DMRHNL),TR1(NL),NL=1,NLYS) A 154

C DMR(NL)=DIFFERENTIAL MANOMETER READING AT NLTH SCREEN AT START OF A 155

C TEST. A 156

C DMR2INL)=DIFFERENTIAL MANOMETER READING AT NLTH SCREEN AT TIME TL2, A 157

C IN MINUTES AFTER START OF TEST. A 158

00 18 NL=1,NLYS A 159

NTL(NL)=1 A 160

DMRHNL)=DMR1(NL)*SGMF/13,56 A 161

18 HRCNL,NT)=HR<NBS,l)+DMRl(NL) A 162

19 READ 67, (OMR2(NL),TR2<NL),NL=1,NLYSI A 163

C DIFFERENTIAL MANOMETER READINGS ARE CONSIDERED POSITIVE WHEN THE A 164

C DGWNHOLE PNEUMATIC HEAD IS GREATER THAN THAT AT LAND SURFACtE. A 165

IF (DMR2(l).GT.99,90) GO TC 22 A 166

DO 21 NL=l,NLYS A 167

M1=NTL(NL) A 168

DMR2(NL)=DMR2iNL)*SG«F/l3.56 A 169

M2={TR2(NL)/TINC)+l A 170

DO 20 NT=M1,M2 A 171

DV=DMRl(NL)+(DMR2(NL)-DMRliNLJ)*{(NT-1>*TINC-TR1CNL))/(TK2(NL)-TRl A 172

1CNL1J A 173

HRtNL,NT)=HR(NBS,NT)+DV A 174

20 CONTINUE A 175

DMRUNL)*DMR2(NL) A 176

TR1CNL1*TR2(NLI A 177

NTL(NL)=M2+1 A 178

21 CONTINUE A 179

GO TO 19 A 180

22 CONTINUE A 181 

C THE FOLLOWING STATEMENTS RELATE TO COMPUTING HEAD AT EACH SPACE NODE A 182 

C AT THE START OF THE TEST. THESE COMPUTATIONS ARE BASED ON FITTING A 183 

C A SINE CURVE THROUGH THE DATA FOR THE BOTTOM LAYER, AND BY LINEAR A 184

C INTERPOLATION FOR THE OVERLYING LAYERS. A 185

DO 23 NL=1,NBS A 186

23 JSD<NLI=<SD(NU/UZD)*DFLOAT<JMAX-1) + 1 A 187

C JSDINLI=NODAL POSITION OF NLTH SCREEN. A 188

NNLC1»=JMAX-JSD(2) A 189

C NNL(NU = NUMBER OF NODES INCLUDED IN NLTH LAYER. A 190

NBOT=JSDtn-JSD(2) A 191

J1=JSD«2I A 192

PI=3.1415926 A 193

ELB=DFLQATCNNLC1I) A 194

XL=NBOT/ELB A 195

HMAX=(HR(l,n-HR(2,in/DSIN{PI*XL/2.) A 196

NL=i A 197

DO 24 J=J1,JMAX A 198

24 HNSCJ)=HR(NL+1,1)+HMAX*DSIN(PI*<J-J1)/{2.*ELBJ) A 199

IF (NLYS.EQ. II GO TO 28 A 200

DO 25 NL=2,NLYS A 201

25 NNL(NLI=JSDtNL)-JSDCNL+ll A 202
34 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

0110

0111

0112

0113

0114

0115

0116

0117

0118

0119

0120

0121

0122

0123

0124

0125

0126

0127

0128

0129

0130

0131

0132

0133

0134

0135

0136

0137

0138

0139

0140

0141

0142

0143

0144

0145

0146

0147

0148

0149

0150

0151

0152

0153

0154

0155

0156

0157

0158

0159

0160

0161

0162

0163

0164

0165

0166

0167

0168

J2=Jl-l A 203

00 27 NL=2,NLYS A 204

J1=JSOINL*1I A 205

DELH=HR(NL,U-HR<NL + 1,1I A 206

DO 26 J=J1,J2 A 207 

HNS(J)aHR<NL+l,l)+DELH*(J-J1I/NNHNL) A 208

26 CONTINUE A 209

J2=Jl-l A 210

27 CONTINUE A 211

28 PRINT 81 A 212

PRINT 100 A 213

PRINT 89 A 214

DO 29 NT=1,NTS,NSTEP A 215

TIME=(NT-1.)*TINC A 216

29 PRINT 71, NT,TIME,tHR(NL,NT),NL=l,NBS) A 217

C THE FOLLOWING STATEMENTS SET UP INITIAL VALUES FOR PERMEABILITY, A 218

C SPECIFIC STORAGE, AND THE STORAGE-TIME TERM FOR THE DIFFERENCE A 219

C EQUATIONS. A 220

DELT=TINC/1440. A 221

DZSQ=<UZD/(JMAX-1.))**2. A 222

J2=JKAX A 223

DO 31 NL=1,NLYS A 224

P(NL)*HYDC(NL,1)*WKV/(386.4*ADV) A 225

Jl=JSDINL+l) A 226

DO 30 J=J1,J2 A 227

30 SS<JI=DPOR(NL)/<PBAR*846.I A 228

J2=J1-1 A 229

31 CONTINUE A 230

Rtl)=SS(l)*DZSQ/DELT A 231

DO 32 J=2,JMAX A 232

32 R< J)=(SS( J)+SS(J-1M*DZSQ/(DELT*2.) A 233

C THE FOLLOWING STATEMENTS ESTABLISH THE COEFFICIENTS FOR THF UN- A 234

C KNOWN HEAD TERMS IN THE FINITE-CIFFERENCE EQUATION. THEY ARE RE- A 235

C COMPUTED FOR EACH TRIAL PERMEABILITY VALUE. A 236

JJ2=JMAX A 237

DO 60 NL=1,NLYS A 238

JM2=JJ2-1 A 239

PRINT 85, NL,NL A 240

TCT2=1 A 241

J1 = JSD(NL + U+1 A 242

JP1=JH-1 A 243

JS=JSD(NL) A 244

33 CONTINUE A 245

DO 58 ICT=1,10 A 246

34 HDIF=0. A 247

DO 35 J=J1,JM2 A 248

A{J)=P(NL) A 249

BtJ)=-(2.*P(NLl+R(J)) A 250

35 C(J)=P(NL) A 251

IF (JJ2.NE.JMAX) GC TO 36 A 252

A(JJ2)=2.*PINL) A 253

B{ JJ2)=-(2.*P(NLi+R(JJ2)) A 254

GO TO 37 A 255

36 A( JJ2)*=P(NL) A 256

Bt JJ2)*-(P(NL) + P(NL-L)+R( JJ2N A 257

C(JJ2)=P(NL-1) A 258

37 CONTINUE A 259

HYDC(NL,ICT)=P(NL)*386.4*ADV/WKV A 260

C THE FOLLOWING STATEMENTS REESTABLISH INITIAL HEAD CONDITIONS FOR A 261

C EACH TRIAL PERMEABILITY VALUE. A 262

DO 38 J=J1,JMAX A 263

38 H{J)=HNS(J) . A 264

C THE FOLLOWING STATEMENTS SOLVE THE DIFFERENCE ECUATIONS FOR A GIVEN A 265

C TRIAL PERMEABILITY, ONE TIME STEP AT A TIME, AND DETERMINE THE SUM A 266

C OF THE DIFFERENCE FROM MEASURED VALUES. A 267

DO 46 NT=2,NTS A 268

D(J1)=-(P(NLI*HR(NL*1,NT)+R(J1)*H(JU) A 269

DO 39 J=JP1,JMAX A 270

39 D{J)=-(R<JJ*H<J)I A 271

CALL TRIMAT(Jl.JMAX) A 272

C DETERMINES WHETHER TO SAVE HEAD MATRIX AT SCREEN FOR PRINTING AFTER A 273

C FINAL HYDRAULIC CONDUCTIVITY VALUE FOR A GIVEN LAYER IS CHOSEN, OR A 274
AIRK A FORTRAN IV PROGRAM 35

0169

0170

0171

0172

0173

0174

0175

0176

0177

0178

0179

0180

0181

0182

0183

0184

0185

0136

0187

0188

0189

0190

0191

0192

0193

0194

0195

0196

0197

0198

0199

0200

0201

0202

0203

0204

0205

0206

0207

0208

0209

0210

0211

0212

0213

0214

0215

0216

C TO ESTABLISH TABLE OF FINAL HEAD VALUES FCR COMPLETE SIMULATION. A 275

IF (NL.EQ.NLYS.AND.KOP(NLYS).EQ.l) GO TO 40 A 276

HSCSUCT,NT)=H(JS) A 277

GO TO 42 A 278

40 CONTINUE A 279

DO 41 INL=1,NLYS A 280

41 HSCSUNL,NT)=H(JSDUNLM A 281

GO TO 46 A 282

C WEIGHTS DIFFERENCE BETWEEN COMPUTED AND MEASURED HEADS WITH TIME, A 283

C ACCORDING TO SPECIFIED IWTO VALUE. A 284

42 IF (IWTO-1) 45 t 44,43 A 285

43 HDIF=HDIF*((H(JS)-HRJ NL,NT) ) *I BSN*DFLOAT (NT/NTSM A 286

GO TO 46 A 287

44 HOIF=HOIF+C(HCJS)-HR(NL,NT))*IBSN*(1.-NT/NTS}) A 288

GO TO 46 A 289

45 HDIF=HDIF + «H( JS1-HR<NL,NT))*IBSN) A 290

46 CONTINUE A 291

C TESTS FOR END OF COMPLETE SIMULATION, AND ROUTES PROGRAM TO FINAL A 292

C TABLE PREPARATION. A 293

IF {NL.EQ.NLYS.AND.KCP(NLYS).EQ.l) GO TO 60 A 294 

C TESTS FOR END OF SIMULATION OF A GIVEN LAYER, AND ROUTES PROGRAM TO A 295

C PRINT HEAD VS TIME TABLES FCR TRIAL SIMULATIONS. A 296

IF (KQPlND.EQ.l) GC TO 53 A 297

C TESTS FOR FIRST TRIAL SIMULATION OF' A GIVEN LAYER. A 298

C TESTS WHETHER A PREVICUS SIMULATION CF THE HEAD CHANGE ACROSS THE A 299

C LAYER HAS BEEN MADE, AND IF SO, ROUTES THE PROGRAM TO COMPAKE THE A 300 

C SIGN OF THE CURRENT HEAC DIFFERENCE KITH THAT FGR THE PREVIOUS SIMU- A 301

C LATION. A 302

IF (ICT.GT.l) GO TC 49 A 303

IF (HDIF) 47,52,48 A 304 

C FOR AN INITIAL SIMULATICNt ASSIGNS A CODE FOR THE SIGN OF HJIF, SAVES A 305

C THE OLD VALUES OF PERMEABILITY AND HOIF, AND INCREASES OR DECREASES A 306

C PERMEABILITY ACCORDING TO THE SIGN OF HDIF. A 307

47 ISIGN=l A 308

PS=P(NL) A 309

HDS=HDIF A 310

P(NL)=P(NL)/(1.+X(ICT2)) A 311

GO TO 58 A 312

48 ISIGN=2 A 313

PS=P(NL) A 314

HDS=HDIF A 315

P(NL)=P(NL)*(1.+X( ICT2M A 316

GO TO 58 A 317

C ESTABLISHES CODE FOR SIGN OF HOIF FOR CURRENT TIME STEP, COMPARES A 318

C PRESENT AND PREVIOUS HDIF SIGNS, AND IF THE SIGN HAS CHANGEJ, INTER- A 319

C POLATES A PERMEABILITY VALUE FRCM THE PRESENT AND PREVIOUS VALUES FOR A 320

C PERMEABILITY AND HOIF. IF HDIF DOES NOT CHANGE SIGN, ROUTEo PROGRAM A 321

C TO SAVE THE SIGN CODE, HDIF, AND PERMEABILITY FOR THE PRESENT SIMU- A 322

C LATION, AND TO INCREASE OR DECREASE THE TRIAL PERMEABILITY VALUE, A 323

C AS APPROPRIATE. A 324

49 IF (HDIF) 50,52,51 A 325

50 ISIGK2=1 A 326

IF «ISIGN-ISIGN2).EC.O) GG TO 47 A 327

P(NL) = PS-<CPS-P(NL))*(DABS{HGS))/{DABS(HOS)+DABS(HDIF))) A 328

GO TO 52 A 329

51 ISIGN2=2 A 330

IF ( USIGN-ISIGN2).EQ.O) GO TO 48 A 331

PCNL)=PS+(CP{NL)-PS)*(DABS(HDS))/(DABS(HDS)+DABS(HDIF))) A 332

52 CONTINUE A 333

C PRINTS OUT TABLES FOR HEAD VERSUS TIME FOR THE TRIAL SIMULATIONS FOR A 334

C A GIVEN LAYER, ONCE A CHANGE IN THE SIGN OF HDIF HAS BEEN DETECTED. A 335

IF (ICT2.LT.2) GO TC 54 A 336

KOP(NL)=1 A 337

IF (NL.LT.NLYS) GO TO 58 A 338

53 PRINT 87 A 339

GO TO 55 A 340

54 PRINT 86 A 341

55 PRINT 88 A 342

PRINT 90, (HYDCCNL.IC),IC=1,ICT) A 343

PRINT 89 A 344

DO 56 NT=N,NTS,NSTEP A 345

TIME=(NT-1.)*TINC A 346
36 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

0217

0218

0219

0220

0221

0222

0223

0224

0225

0226

0227

0228

0229

0230

0231

0232

0233

0234

0235

0236

0237

0238

0239

0240

0241

0242

0243

0244

0245

0246

0247

0248

0249

0250

0251

0252

0253

0254

0255

0256

0257

0258

0259

0260

0261

0262

0263

0264

0265

56 PRINT 71, NT,TIME, <HSCS(1C,NT), IO1,ICTJ A 347

IF UCT2.E6.2) GO TC 57 A 348

ICT2=2 A 349

GO TC 33 A 350

57 HYDS(NL)=P<NL)*386.4*ACV/WKV A 351

IF (NL.EQ.NLYS) GO TC 34 A 352

JJ2=J1-1 A 353

GO TO 60 A 354

58 CONTINUE A 355

C PRINTS OUT CONVERGENCE FAILURE MESSAGE AND RESULTS OF THE TEN TRIAL A 356

C SIMULATIONS MADE WHEN NC CHANGE IN THE SIGN OF HDIF OCCURS. A 357

PRINT 70, NL,HYQC(NL,ICT) A 35b

PRINT 90, (HYDC(NL,IC),IC=1,ICT) A 359

DO 59 NT=N,NTS,NSTEP A 360

TIME=(NT-1.)*TINC A 361

59 PRINT 71, NT,TIME,(HSCSl1C ,NT),IC = 1,ICT) A 362

GO TO 62 A 363

60 CONTINUE A 364

C -PRINTS OUT COMPUTED PNEUMATIC HEAD VALUES AT EACH SCREEN FOR THE A 365

C FINAL SIMULATION, A 366

PRINT 91 A 367

PRINT 100 A 368

PRINT 89 A 369

DO 61 NT=N,NTS,NSTEP A 370

TIME=(NT-1«)*TINC A 371

61 PRINT 71, NT,TIME, (HSCS(KL,NT),NL=I,NLYS), HP(NBS,NT) A 372

C ESTABLISHES CODE FOR FINAL TABLE AND ROUTES PROGRAM TO PRINT TEST A 373

C INPUT DATA AND A TABLE SHCWING THE FINAL HY9R4ULIC CONDUCT IVITY CHOSE A 374

C FOR EACH LAYER, FOLLOWED BY PROGRAM TERMINATION. A 375

MPT=1 A 376

GO TO 3 A 377

62 CONTINUE A 378

STOP A 379

C A 382

63 FORMAT (12,F7.0,14,F7.0, 14 ) A 383

64 FORMAT (2F7.2) A 384

65 FORMAT (F5.1,F8.4,F4.?,2X,I 1) A 385

66 FORMAT (F6.3,F5.0) A 386

67 FORMAT {&<F6.3 ,F5.0) ) A 387

68 FORMAT (I1,?D10.3) A 388

69 FORMAT (212) A 389

70 FORMAT (' FAILED TC CONVERGE AT»,13,'LAYER. LAST HYDC WAS',E10.3,« A 390

1 FEET/DAY.') A 391

71 FORMAT (1X,I3,F8.0,10F12.4) A 392

72 FORMAT {/,39X,'THE UNSATURATED ZONE IS',F7.0,' FEET THICK AT THIS A 393

1SITE.M ' A 394

73 FORMAT {/,32X,'FOR THE ANALYSIS, THIS ZONE WAS DIVIDED INTCSI3,* A 395

1LAYERS AND 1 , 15, MODES.») A 396

74 FORMAT (/,34X,'FGR PURPOSES CF COMPUTATION, THIS INTERVAL WAS DIVI A 397

10EO«,/,40X,«INTO',15,' TIME STEPS GF«,F5.0,' MINUTES EACH. 1 ) A 398

75 FORMAT (27X,12,8X,F4 .0,'-»,F4.0,1 IX,010.3, 13X,F4.2) A 399

76 FORMAT (1H1,24X,«ANALYSIS OF PNEUMATIC HEAD DATA AT THES8A4,' SIT A 400

IE TO',/,35X, DETERMINE VERTICAL HYDRAULIC CCNDUCTIVITY uF THE UNS A 401

2ATURATED ZONE.') A 402

77 FORMAT (////,30X,« DATA WERE COLLECTED FRO* ',A4, TO ",A4,« HOURS A 403

1 ON ,3A4,////) A 404

73 FORMAT 12 ( A4 ,2X) ,3 A4) A 405

79 FORMAT (25X,'TABLE 1, INPUT DATA ON LAYER BOUNDARIES, INITIAL EST A 406

1IMATES OF HYDRAULIC',/,41X,'CONDUCTIVITY, AND AIR-FILLED POROSITY. A 407

2' ) A 408

80 FORMAT (57X,'ESTIMATED' f /,25X, LAYER ',3X,'LAYER INTERVAL HYDRAUL A 409

1IC CONDUCTIVITY AIR-FILLED',/,25X,'NUMBER (FEET) A 410

2 (FEET PER DAY) POROSITY',/) A 411

81 FORMAT (1H1,» TABLE 2.--MEASURED OR INTERPOLATED PNEUMATIC HEAD DA A 412

1TA, IN INCHES OF MERCURY, AT EACH SCREEN.') A 413

82 FORMAT <8A4) A 414

83 FORMAT (25X,• TABLE 5—FINAL ESTIMATES OF HYDRAULIC CONDUCT IVITY, A 415

IBASED ON THE ASSUMED',/,40X,• LISTED VALUES OF AIR-FILLED POROSITY A 416

2.') A 417

84 FORMAT (/////) A 418

85 FORMAT (IHltlOX,* ASSUMED VALUES FOR HYDRAULIC CONDUCTIVITY FOR LA A 419

1YER',I4,' AND THE RESULTING COMPUTED',/,30X,'HEAD VALUES AT SCREEN A 420
A1RK-A FORTRAN IV PROGRAM 37

0266

0267

0268

0269

0270

0271

0272

0273

0274

0275

0276

0277

0278

0279

0280

0281

0282

0283

2',13,' ARE LISTED BELOW.') A 421

86 FORMAT <//,l5X,'VALUES FOR THE FIRST SET OF COMPUTATIONS (HYORAULI A 422

1C CONDUCTIVITY HALVED OR DOUBLED EACH TIME.)') A 423

87 FORMAT (////,15X,' VALUES FOR THE SECCND SET OF COMPUTATIONS (HYDR A 424

1AULIC CONDUCTIVITY CHANGED BY TEN PER CENT EACH TIME.)') A 425

88 FORMAT {/,20X,» HYDRAULIC CONDUCTIVITY, FEET PER DAY') A 426

89 FORMAT {' TIME TIMEt',/,' STEP MIN',6X,'PNEUMATIC HEAu VALUES, A 427

UN INCHES OF MERCURY') A 428

90 FORMAT 112X, lOi 2X,C10.3)) A 429

91 FORMAT UH1, 16X, 'TABLE 4—PNEUMATIC HEADS AT EACH SCREEN, BASED ON A 430

1 ATMOSPHERIC PRESSURE CHANGES',/,20X,'AND THE FINAL COMPJTER-SELEC A 431

2TED HYDRAULIC CONDUCTIVITY VALUES LISTED IN TABLE 5 BELO«') A 432

92 FORMAT <30X,'FOR THIS SET CF COMPUTATIONS, ALL READINGS ARE WEIGHT A 433

1ED EQUALLY') A 434

93 FORMAT (15X,'FOR THIS SET OF COMPUTATIONS, READINGS OF PNEUMATIC H A 435

1EAD ARE WEIGHTED INVERSELY WITH TIME') A 436

94 FORMAT (15X,'FCR THIS SET OF COMPUTATIONS, READINGS OF PNEUMATIC H A 437

1FAD ARE WEIGHTED DIRECTLY WITH TIME') A 438

95 FORMAT (30X,' THE 8APOMETER WAS RISING DURING MOST OF THE READING A 439

1PERIOD.M A 440

96 FORMAT <30X,' THE BAROMETER DECLINED DURING MOST OF THc READING A 441

1PERIOD.' ) A 442

97 FORMAT <22X, 'MEAN STATION PRESSURE DUPING THESF OBSERVATIONS WAS • A 443

1,F5.2,' INCHES OF MERCURY.') A 444

98 FORMAT (25X,'THE MANCMETER FLUIC USED FOP THIS TEST HAS A SPECIFIC A 445

1 GRAVITY OF »,F5.2) A 446

99 FORMAT <//,25X,'THE BOTTCM SCREEN IS AT THE «,F5.0,' FOOT DEPTH BE A 447

1LOW LAND SURFACE') A 448

100 FORMAT (IOX,'PNEUMATIC HEAC VALUES ARE LISTED IN SEQUENCE FROM SCR A 449

1EEN 1 ON THE LEFT TO LAND SURFACE ON THE RIGHT.') A 450

101 FORMAT (20X,'THE DYNAMIC VISCOSITY VALUE FOR AIR USED IN THIS RUN A 451

US »,D10,3,' LB*SEC/SQUARE FCOT, ' ,/, 19X, ' AND THE VALUE FUR THE KIN A 452

2EMATIC VISCOSITY OF WATER IS ',010.3,' SQUARE FEET PER SECOND.') A 453

102 FORMAT (30X,'THESE VISCOSITY VALUES ARE FOR A TEMPERATURE OF 60 D£ A 454

1GREES F,',/,43X,'ANC WERE REAC IN BY THE DEFAULT OPTION.') A 455

END A 456-

TRIMAT SUBROUTINE LISTING

FORTRAN IV G LEVEL 21 TPIMAT DATE = 77318 18/44/45

0001

0002

0003

0004

0005

0006

0007

0008

0009

0010

0011

0012

0013

0014

0015

0016

0017

0018

0019

SUBROUTINE TRIMATCJl,JNAX) 8 1

C THE TRIMAT SUBROUTINE SOLVES THE SET GF IMPLICIT EQUATIONS JSED TO 82

C DESCRIBE THE HEAD DISTRIBUTION AT THE END OF A GIVEN TIME STEP, BASED B 3

C CN THE THOMAS ALGORITHM (VCN ROSENBEPG, 1969). B 4

C REFERENCE—VON ROSENBERG, DALE, 1969, METHODS FOR THE NUMERICAL SOLU- B 5

C TION OF PARTIAL DIFFERENTIAL EQUATIONS: NEW YORK, AMERICAN ELSEV- B 6

C 1ER, P. 113. B 7

IMPLICIT REAL*8(A-H,C-Z) B 8

COMMON /COM1/ H(200),A(200),B(200),C{200),0(200) B 9

DIMENSION BETA1200), GAMMAC200) B 10

JPI=JI*I e 11

BETACJIMBUl) B 12

GAMMA(J1)=DCJ1)/B<Jl) E 13

00 1 I=JP1,JMAX B 14

BETA(I)=B(I)-A(I)*C( I-1)/BETAU-1) B 15

GAMMAUMCDU )-A(I )*GAPMA( 1-1) }/6ETA( 1 ) B 16

1 CONTINUE B 17

HtJMAX)=GAMMA{JMAX) B 18

JLIM=JMAX-J1 B 19

DO 2 JJ=1,JLIM B 20

J=JMAX-JJ B 21

H(J)=GAMMA(J)-C(J)*H(J+1)/BETA(J) B 22

2 CONTINUE B 23

RETURN B 24

END B 25-
38 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

SAMPLE COMPUTER OUTPUT

Output from the AIRK program is in five tables. Table 1 

lists the input data, with appropriate headings, from cards 

1-4 and 6-8. These data provide a means of checking the 

input data for accuracy.

Table 2 lists the pneumatic heads at each screen at the end 

of specified time steps used in the analysis, in inches of 

mercury. These data are computed within the program from 

the barometric data read in card series 9, and from the 

differential manometer values read in card series 10. 

Readings for a given time are listed from the bottom screen 

in the first column of head values sequentially up to the 

atmospheric pressure in the rightmost column.

The third table (unnumbered) is actually a series of NLYS 

tables, and lists the equivalent hydraulic conductivity values 

used in the programed search procedure for a given layer, 

and in a column immediately below each value, lists

computed head values for the included or underlying screen 

based on that value. Pneumatic head values in this table may 

be compared to those in table 2 to judge the adequacy of the 

simulation. Moreover, these head values may be used to 

make a manual selection of the most appropriate hydraulic 

conductivity value if the fit chosen by the automatic search 

procedure is inadequate. The first page of twelve pages for 

this table is shown in the sample output below.

Table 4 lists pneumatic heads computed at each screen 

based on the final computer-selected air permeability values 

for each layer and on the observed barometric pressure 

changes at land surface. Values in table 4 may be compared 

to those in table 2 to judge the overall adequacy of the 

simulation.

Table 5 lists the selected hydraulic conductivities in the 

same format as in table 1 to provide a summary of the 

computer analysis.

SAMPLE OUTPUT DATA

ANALYSIS OF PNEUMATIC HEAD DATA AT THELUBBOCK AIRPORT SPREADING SITE B SITE TO 

DETERMINE VERTICAL HYDRAULIC CONDUCTIVITY dF THE UNSATURATED ZONE.

THE UNSATURATED ZONE IS 125. FEET THICK AT THIS SITE. 

FOR THE ANALYSIS, THIS ZONE WAS DIVIDED INTO o LAYERS AND 126 NODES.

DATA WERE COLLECTED FRCM 1215 TO 1925 HOURS ON ,-lAY 17, 1972

FOR PURPOSES CF COMPUTATION, THIS INTERVAL WAS DIVIDED

INTO 79 TIME STEPS OF 5. MINUTES EACH.

MEAN STATION PRESSURE CURING THESE OBSERVATIONS WAS 26.50 INCHES OF MERCURY. 

THE MANOMETER FLUID USED FOR THIS TEST HAS A SPECIFIC GRAVITY OF 1.20

THE BAROMETER DECLINED DURING MOST OF THE REAuING PERIOD. 

FOR THIS SET OF CCMPUTATICNS, ALL READINGS ARE WEIGHTED EQUALLY

THE DYNAMIC VISCOSITY VALUE FOR AIR USEC IN THIS RUN IS 0.374D-06 L8*SEC/SQUARE FOOT, 

AND THE VALUE FOR THE KINEMATIC VISCCSITY OF WATER IS 0.121D-04 SQUARE FEET PER SECOND, 

THESE VISCOSITY VALUES ARE FOR A TEMPERATURE OF 60 DEGREES F t 

AND WERE READ IN BY THE DEFAULT OPTION.

TABLE 1. INPUT DATA ON LAYER BCUNCARIES, INITIAL ESTIMATES OF HYDRAULIC 

CCNCUCTIVITY, AND AIR-FILLED POROSITY.

ESTIMATED

LAYER LAYER INTERVAL HYDRAULIC CONDUCTIVITY AIR-FILLED 

NUMBER (FEET) (FEET PER DAY) POROSITY

98.-125, 

90.- 98, 

72.- 90, 

51.- 72, 

32.- 51 

0.- 32,

0.100D-01 

0.100D 00 

0.100D 02 

0.100D 02 

0.100D 02 

0.100D 02

0.22

0.03

0.21

0.15

0.20

0.20

THE BOTTOM SCREEN IS AT THE 115. FOOT DEPTH BELOW LA.4D SURFACE
A1RK-A FORTRAN IV PROGRAM 39

TABLE 2.— MEASURED OR INTERPOLATED PNEUNATIC HEAC DATA, IN INCHES Of MERCURY, AT EACH SCREEN. 

PNEUMATIC HEAD VALUES ARE LISTED IN SEQUENCE FROM SCREEN 1 ON THE LEFT TO LAND 

SURFACE ON THE RIGHT.

TIME TIPE,

STEP

1

4

7

10

13

16

19

22

25

23

31

34

37

40

43

46

49

52

55

58

61

64

67

70

73

76

79

WIN

0.

15,

30.

45.

60.

75.

90.

105.

120.

135.

150.

165.

180.

195.

210.

225.

240.

255.

270.

285.

300.

315.

330.

345.

360.

375.

390.

PNEUMATIC

0,9695

0.9676

0.9679

0.9675

0.9646

0.9623

0.9632

0.9624

0.9613

0.9618

0.9613

0.9596

0.9590

0.9574

0.9545

0.9519

0.9485

0.9462

0.9456

0.9433

0.9413

0.9345

0.9292

0.9241

0.9206

0.9181

0.9144

HEAD VALUES,

0.9615

0.9599

0.9593

0.9590

0.9526

0.9501

0.9497

0.9478

0.9423

0.9409

0.9402

0.9358

0.9330

0.9262

0.9189

0.9201

0.9186

0.9147

0.9117

0.9088

0.9062

0.8975

0.8900

0.8865

0.8837

0.8810

0.8763

IN INCHES CF MERCURY

0.9580

0.9556

0.9536

0.9519

0.9460

0.9431

0.9416

0.9372

0.9314

0.9286

0.9264

0.9195

0.9175

0.9137

0.9074

0.9054

0.9023

0.8972

0.8933

0,8898

0.8856

0.8768

0.8704

0.8676

0.8654

0.8631

0.8586

0.9562

0.9535

0.9510

0.9491

0.9431

0.9398

0.9384

0.9322

0.9270

0.9244

0.9228

0.9161

0.9132

0.9094

0.9028

0.9012

0.8977

0.8916

0.8865

0.8828

0.8777

0.8697

0.8666

0.8628

0.8597

0.8577

0.8551

0.9549

0.9522

0.9490

0.9469

0.9416

0.9366

0.9350

0.92S7

0.9230

0.9196

0.9181

0.9115

0.9077

0.9050

C.9006

0.8959

0.8922

0.8863

0.8803

0.8753

0.8/07

0.8635

0.8602

0.8583

0.8567

0.8555

0.8524

0.9527

0.9500

0.9468

0.9439

0.9388

0.9341

0.9318

0.9240

0.9184

0^9146

0.9113

0.9073

0.9034

0 .8996

0.8940

0.8905

0.8868

0.3804

0.8741

0.8688

0.8641

0.8571

0.8545

0.8541

0.8538

0.8533

0.8502

0.9500

0.9470

0.9440

0.9410

0.9360

0.9300

0.9260

0.9190

0.9140

0.9100

0.9060

0.9030

O.d990

0.8950

0.8900

0.8850

0.8800

O.d730

0.8670

0.8620

0.8580

0.8520

O.d500

0.8500

0.8500

0.8500

0.8430

ASSUMED VALUES FOR HYDRAULIC CONDUCTIVITY FOR LAYER 1 AND THE RESULTING COMPUTED 

HEAD VALUES AT SCREEN 1 ARE LISTED BELOW.

VALUES FOR THE FIRST SET CF CCMPUTATICNS 

(HYDRAULIC CONDUCTIVITY HALVED OR DOUBLED EACH TIME.)

HYDRAULIC CONDUCTIVITY, FEET PER DAY 

0.100D-01 0.2000-01 0.400D-01 0.800D-01 0.1600 00 0.3200 00

TIME

STEP

4

7

10

13

16

19

22

25

28

31

34

37

40

43

46

49

52

55

58

61

64

67

70

73

76

79

TIME,

MIN

15.

30.

45.

60.

75.

90.

105.

120.

J35.

150.

165.

180.

195.

210.

225.

240.

255.

270.

285.

300.

315.

330.

345.

360.

375.

390.

PNEUMATIC

0.9694

0.9694

0.9694

0.9694

0.9694

0.9693

0.9693

0.9693

0.9693

0.9692

0.9692

0.9692

0.9692

0.9692

0.9691

0.9691

0.9691

0.9691

0.9690

0,9690

0.9690

0.9689

0.9689

0.9688

0.9688

0.9687

HEAD VALUES,

0.9694

0.9694

0.9693

0.9693

0.9692

0.9692

0.9692

0.9691

0.9691

0.9690

0.9690

0,9689

0.9688

0.9688

0.9687

0.9686

0.9685

0.9684

0.9682

0.9681

0.9679

0.9678

0.9676

0.9674

0.9672

0.9669

IN INCHES CF

0.9694

0.9693

0.9692

0.9691

0.9690

C.9689

0.9688

0.9687

0.9686

0.9684

0.9682

0.9680

0.9678

0.9675

0.9672

0.9669

0.9665

0.9661

0.9657

0.9652

0.9647

0.9641

0.9636

0.9630

0.9623

0.9616

MERCURY

0.9693

0.9691

0.9689

0.9687

0.9685

0.9683

0.9680

0.9676

0.9672

0.9667

0.9662

0.9657

0.9651

0.9644

0.9636

0.9628

0.9619

0.9609

0.9599

0.9589

0.9578

0.9567

0.9555

0.9542

0.9528

0.9513

0.9691

0.9688

0.9684

0.9680

0.9674

0.9668

0.9661

0.9653

0.9644

0.9634

0.9624

0.9613

0.9601

0.9588

0.9572

0.9556

0.9539

0.9523

0.9505

0.9488

0.9469

0.9450

0.9428

0.9404

0.9380

0.9355

0.9688

0.9680

0.9673

0.9665

0.9653

0.9641

0.9628

0.9614

0.9598

0.9582

0.9566

0.9548

0.9529

0.9506

0.9479

0.9455

0.9431

0.9407

0.9382

0.9357

0.9331

0.9301

0.9267

0.9232

0.9198

0.9164
40 VERTICAL PERMEABILITY TO AIR IN THE UNSATURATED ZONE

VALUES FOR THE SECOND SET CF COMPUTATIONS 

(HYDRAULIC CONDUCTIVITY CHANGED BY TEN PER CENT EACH TIME.)

HYDRAULIC CCNCUCTIVITY, FEET PER DAY 

0.296D 00 0.2690 00 0.291D 00

TIME

STEP

ft

7

10

13

16

19

22

25

28

31

34

37

40

43

46

49

TIKE,

MIN

15.

30.

45.

60.

75.

90.

105.

120.

135.

150.

165.

180.

195i

210.

225.

240.

PNEUMATIC

0.9688

0.9681

0.9674

0.9667

0.9656

0.9644

0.9632

0.9619

0.9604

0.9588

0.9573

0.9556

0.9538

0.9516

0.9491

0.9467

HEAD VALUES,

0.9689

0,9683

0.9676

0.9669

0.9660

0.9649

0.9637

0.9625

0.9611

0.9596

0.9582

0.9566

0.9549

0.9528

0.9505

0.9482

IN INCHES CF MERCURY

0.9688

0,9682

0.9675

0.9667

0.9657

0.9645

0.9633

0.9620

0.9606

0.9590

0.9575

0.9558

0.9540

0.9518

0.9494

0.9470

TABLE 4—PNEUMATIC HEADS AT EACH SCREEN, BASED ON ATMOSPHERIC PRESSURE CHANGES

AND THE FINAL COMPUTER-SELECTED HYDRAULIC CONDUCTIVITY VALUES LISTED IN TABLE 5 8ELOW

PNEUMATIC HEAD VALUES ARE LISTED IN SEQUENCE FROM SCREEN 1 ON THE LEFT TO

LAND SURFACE ON THE RIGHT.

TIME

STEP

4

7

10

13

16

19

22

25

28

31

34

37

40

43

46

49

52

55

58

61

64

67

70

73

76

79

TIME,

MIN

15.

30.

45.

60.

75.

90.

105.

120.

135.

150.

165.

180.

195.

210.

225.

240.

255.

270.

285.

300.

315,

330.

345.

360.

375.

390.

PNEUMATIC

0.9688

0.9682

0.9675

0.9667

0.9658

0.9648

0.9637

0.9624

0.9610

0.9594

0.9577

0.9558

0.9539

0.9518

0.9496

0.9473

0.9450

0.9425

0.9399

0.9372

0,9343

0.9313

0.9283

0.9252

0.9220

0.9189

HEAD VALUES,

0.9603

0.9589

0.9572

0.9554

0.9531

0.9505

0.9476

0.9443

0.9408

0.9373

0.9339

0.9306

0.9273

0.9241

0.9206

0.9171

0.9133

0.9093

0.9050

0.9007

0.8964

0.8922

0.8882

0.8849

0.8821

0.8797

IN INCHES CF MERCURY

0.9559

0.9537

0.9514

0.9488

0.9456

0.9418

0.9378

0.9334

0.9288

0.9245

0,9205

0.9168

0.9131

0.9093

0.9053

C.9010

0.8965

0.8915

0.8863

0.8814

0.8765

0.8719

0.8680

0.8652

0.8632

0.8616

0.9541

0.9517

0.9491

0.9462

0.9425

0.9384

0.9339

0.9291

0.9244

0,9200

0.9160

0.9123

0.9085

0.9045

0.9002

0.8956

0.8909

0.8855

0.3802

0.8752

0.8703

0.8657

0.8622

0.8600

0.8584

0.8570

0.9521

0.9494

0.9467

0.9433

0.9390

0.9347

0.9298

0.9246

0.9199

0.9157

0.9118

0.9081

0.9042

0.9000

0.8955

0.8909

0.8855

0.8798

0.8744

TJ.8696

0.8645

0.8603

0.8577

0.8561

0.8551

0.8538

0.9497

0.9469

0.9440

0.9400

0.9349

0.9305

0.9248

0.9195

0.9151

0.9109

0.9074

0.9036

0.8997

0.8952

0.8904

0.8856

0.8796

0.8736

0.8683

0.8638

0.8585

0.8550

0.8535

0.8527

0.8522

0.8508

0.9470

0.9440

0.9410

0.9360

0.9300

0.9260

0.9190

0.9140

0.9100

0.9060

0.9030

0.8990

0.8950

0.8900

0.8850

0.8800

0.8730

0.8670

0.8620

0.8580

0.8520

0.8500

0.8500

0.8500

0.8500

0.8480

TABLE 5—FINAL ESTIMATES OF HYDRAULIC CONDUCTIVITY, dASED ON THE ASSUMED 

LISTED VALUES CF AIR-FILLED POROSITY.

ESTIMATED

LAYER LAYER INTERVAL HYDRAULIC CONDUCTIVITY AIR-FILLED 

NUMBER (FEET) (FEET PER DAY) POROSITY

98.-12 5, 

90.- 98, 

72.- 90, 

51.- 72, 

32.- 51, 

0.- 32.

0.2910 00 

0.416D 00 

0.4400 01 

0.890D 01 

0.1050 02 

0.2300 02

0.22

0.03

0.21

0.15

0.20

0.20
REFERENCES CITED 41

REFERENCES CITED

Aronovici, V. S., Schneider, A. D., and Jones, O. R,, 1972, Basin recharge

of the Ogallala aquifer: Am. Soc. Civil Engineers Proc., Jour.

Irrigation and Drainage Div., v. 98, no. 1, p. 65-76. 

Boardman, C. R., and Skrove, J. W., 1966, Distribution in fracture

permeability of a granitic rock mass following a contained nuclear

explosion: Jour. Petroleum Technology, v. 18, no. 5, p. 619-623. 

Boersma, L., 1965, Field measurement of hydraulic conductivity above a

water table: in C.A. Black, ed., Methods of soil analysis: Am. Soc.

Agronomy, Agronomy Mon. 9, p. 242-248, Madison, Wis. 

Botset, H. G., 1940, Flow of gas-liquid mixtures through consolidated sand:

AIME Trans., v. 136, p. 91-96. 

Bouwer, Herman, 1961, A double-tube method for measuring hydraulic

conductivity of soil in sites above a water table: Soil Sci. Soc. America

Proc., v. 25, p. .334-342.

 1969, Planning and interpreting soil permeability measurements: Am.

Soc. Civil Engineers Proc., Jour. Irrigation and Drainage Div., v. 95,

no. IR3, p. 391-402. 

Buckingham, Edgar, 1904, Contributions to our knowledge of the aeration

of soils: U. S. Dept. Agriculture Soils Bur. Bull. 25, 52 p. 

Carslaw, H. S., and Jaeger, J. C., 1959, Conduction of heat in solids:

London, Oxford University Press, 510 p. 

Childs, E. C., 1969, An introduction to the physical basis of soil water i

phenomena: New York, John Wiley and Sons, Inc., 493 p. I 

Cronin, J. G., 1964, A summary of the occurrence and development of j

ground water in the Southern High Plains of Texas: U.S. Geol. Survey

Water-Supply Paper 1693, p. 38-39. 

Crow, E. L., Davis, F. A., and Maxfield, M. W., 1960, Statistics manual:

New York, Dover Publications, 288 p. 

Evans, D. D.,and Kirkham, Don, 1949, Measurement of air permeability of

soil in situ: Soil Sci. Soc. America Proc., v. 14, p. 65-73. 

Gile, L. H., Peterson, F. F., and Grossman, R. B., 1966, Morphological and

genetic sequences of carbonate accumulation in desert soils: Soil Sci.,

v. 101, p. 347-360.

 Green, D. W., Dabiri, H. E., Weinaug, C. F., and Prill, R. C., 1970,

Numerical modeling of unsaturated ground water flow and comparison of the model to a field experiment: Water Resources Research,

v. 6, no. 3, p. 862-874. 

Grover, B. L., 1955, Simplified air permeameters for soil in place: Soil Sci.

Soc. America Proc., v. 19, p. 414-418. 

Hefez, E., Shamir, U., and Bear, J., 1975, Identifying the parameters of an j

aquifer cell model: Water Resources Research, v. 11, no. 6, p. 993-1004. j 

Hubbert, M. K., 1940, The theory of ground-water motion: Jour. Geology,

v. 48, no. 8, p. 785-944. 

Jacob, C. E., 1963, Determining the permeability of water-table aquifers://?

Bentall, Ray, compiler, Methods of determining permeability,

transmissibility, and drawdown: U. S. Geol. Survey Water-Supply

Paper 1536-1, p. 245-271. 

Johnston, Norris, and Beeson, Carrol, 1945, Water permeability of

reservoir sands: Petroleum Tech., v. 8, no. 3, p. 43-55. 

Karplus, W. J., 1958, Analog simulation: New York, McGraw-Hill, p.

356-360. 

Katz, D. L., and others, 1959, Handbook of natural gas engineering: New

York, McGraw-Hill, 802 p. 

Keys, W. S., and Brown, R. F., 1971, The use of well logging in recharge

studies of the Ogallala Formation in West Texas: U. S. Geol. Survey

Prof. Paper 750-B, p. B270-B277. 

King, L. O., 1968, Mathematical models for underground injection of

gaseous wastes into the vadose zone: Batelle Northwest Laboratories

Pub. BNWL-945, 55 p. 

Kirkham, Don, 1946, Field methods for determination of air permeability

of soil in its undisturbed state: Soil Sci. Soc. America Proc., v. 11, p.

93-99. 

Klinkenberg, L. J., 1941, The permeability of porous media to liquids and

gases: Am. Petroleum Inst. Drilling and Production Practice, p. 200.

Leverett, M. C., and Lewis, W. B., 1941, Steady flow of gas-oil-water

mixtures through unconsolidated sands: AIME Trans., v. 142, p. 107. 

Lohman, S. W., and others, 1972, Definitions of selected ground-water

terms revision and conceptual refinements: U. S. Geol. Survey

Water-Supply Paper 1988, p. 9. 

McNeal, B. L., and Coleman, N. T., 1966, Effect of solution composition on

soil hydraulic conductivity: Soil Sci. Soc. America Proc., v. 30, p.

308-312. 

McNeal, B. L., and others, 1968, Factors influencing hydraulic conductivity

of soils in the presence of mixed-salt solutions: Soil Sci. Soc. America

Proc., v. 32, p. 187-190. 

Morris, G. A., and Snoeberger, D. F., 1971, Calculations of pressure change

at depth in a nuclear chimney following atmospheric pressure change:

Lawrence Livermore Lab. Report UCID-15963, 6 p. 

Muskat, Morris, 1946, Flow through porous media: New York,

McGraw-Hill, p. 72. 

Osoba, J. S., and others, 1951, Laboratory determination of relative

permeability: AIME Trans., v. 192, p. 47-51. 

Peck, R. L., 1964, Chem-Comp-new shrinkage-compensating cement for

crack-resistant concrete: Modern Concrete, v. 28, no. 6, p. 54-56. 

Prill, R. C., 1976, Movement of moisture in the unsaturated zone in a

loess-mantled area, southwestern Kansas: U. S. Geol. Survey

Open-File Report 76-611, 73 p.

Prill, R. C., and Aaronson, D. B., 1973, Flow characteristics of a subsurface-controlled recharge basin on Long Island, New York: U. S.

Geol. Survey Jour. Research, v. 1, no. 6, p. 735-744. 

Quirk, J. P., and Schofield, R. K., 1955, The effect of electrolyte

concentration on soil permeability: Jour. Soil Sci., v. 6, no. 2, p.

163-178. 

Reeve, R. C., 1953, A method for determining the stability of soil structure

based upon air and water permeability measurements: Soil Sci. Soc.

America Proc., v. 17, p. 324-329. 

Reeve, R. C., and Tamaddoni, G. H., 1965, Effect of electrolyte

concentration on laboratory permeability and field intake rate of a

sodic soil: Soil Sci., v. 99, no. 4, p. 261-266. 

Reeve, R. C., and others, 1954, A comparison of the effects of exchangeable

sodium and potassium upon the physical condition of soils: Soil Sci.

Soc. America Proc., v. 18, p. 130-132. 

Robertson, J. B., 1969, Behavior of xenon-133 gas after- injection

underground: U. S. Geol. Survey open-file report, 37 p. 

Rogers, A. E., and Connolly, T. W., 1960, Analog computation in

engineering design: New York, McGraw-Hill, 443 p. 

Rosza, R. B., Snoeberger, D. F., and Baker, J., 1975, Permeability of a

nuclear chimney and surface alluvium. Area 2, ERDA NTS: Lawrence

Livermore Lab. Report UCID-16722, 11 p. 

Schmalz, B. L., (ed.), 1969, Injection of gas into the lithosphere at the

National Reactor Testing Station: AEC Research and Development

Report IDO-12069, 166 p. 

Snoeberger, D. F., Morris, C. J., and Baker, J., 1972, Chimney permeability

by atmospheric pressure change-instruments and data handling:

Lawrence Livermore Lab. Report UCID-16154, 14 p. 

Stallman, R. W., 1967, Flow in the zone of aeration, in Yen Te Chow, ed.,

Advances in Hydroscience: New York, Academic Press, v. 4, p.

151-195.. 

Stallman, R. W., and Weeks, E. P., 1969, The use of atmospherically

induced gas-pressure fluctuations for computing hydraulic conductivity of the unsaturated zone (abs.): Geol. Soc. America Abs. with

Programs, pt. 7, p. 213. 

Tanner, C. B., and Wengel, R. W., 1957, An air permeameter for field and

laboratory use: Soil Sci. Soc. America Proc., v. 21, p. 663-664. 

Terzhagi, Karl, 1943, Theoretical soil mechanics: New York, John Wiley

and Sons, Inc., p. 243-244. 

von Rosenberg, D. U., 1969, Methods for the numerical solution of partial

differential equations: New York, American Elsevier, 125 p. 

Wenzel, L. K., 1942, Methods for determining permeability of water-bearing materials: U. S. Geol. Survey Water-Supply Paper 887, p. 7.

*U.S. Government Printing Office: 1978-777-102/32