PREFERENTIAL FLOW THROUGH EARTHEN LANDFILL COVERS
By
Duraisamy Soundararajan Saravanathiiban

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
Civil Engineering - Doctor of Philosophy
2014

ABSTRACT
PREFERENTIAL FLOW THROUGH EARTHEN LANDFILL COVERS
By
Duraisamy Soundararajan Saravanathiiban
In order to minimize infiltration of precipitation into the waste, final covers are constructed once
a landfill reaches its capacity. In and to sub-humid climates, earthen final covers have been
permitted. However, formation of macropores (or shrinkage cracks) and flow through relatively
large pores can significantly increase percolation through earthen landfill covers during service.
Most commonly used water balance models that are used for predicting percolation through
earthen cover are based on Richards’ equation that simulates only the micropore flow. Hence, a
model that can simulate the micropore and macropore flow is required to simulate the long-term
hydrology of earthen covers. In this study, validation of Root Zone Water Quality Model
(RZWQM) using data collected form an instrumented field-scale test section and development of
a model capable of simulating macropore flow using lattice Boltzmann method were carried out.
An instrumented field-scale test section of an earthen landfill cover, made up of 1.5 m thick
compacted clay overlain by 0.3 m thick topsoil, was constructed at a landfill located in Detroit,
MI and monitored for about four years. Measured annual percolation increased by an order of
magnitude during the second and the third year of service. Controlled irrigation tests conducted
on the test section confirmed macropore dominated flow through the test section. Estimated
effective field hydraulic conductivities of the test section increased by an order of magnitude
during the 4th year of service compare to the 1st year of service. Field methane tracer tests
confirmed the presence and locations of macropores.

Water balance of the field test section was simulated using the model RZWQM and a
commonly used numerical model UNSAT-H. For the first year data, both models simulated
percolation relatively accurately. However, the numerical predictions of percolation were not
accurate for the second and the third year when the effect of macropores was ignored. The
macropore parameters required for RZWQM were calibrated using field irrigation test results.
RZWQM predictions using calibrated macropore parameters yielded relatively accurate
prediction of percolation for the second and the third years of the field data.
Measurement of macropore flow through clay in lab-scale samples is relatively challenging
due to the effect of confining walls and the relatively small size of the sample. A new laboratory
technique was developed to consistently fabricate clay samples containing macropores. High
resolution X-ray CT images of compacted clay specimens were taken to visualize the 3D
structure of macropores. A 3D lattice Boltzmann model (LBM) that can simulate saturated flow
through micropore and macropore flow was developed. Verification of the LBM was carried out
using analytical solutions. The LBM was validated using laboratory measurement of saturated
hydraulic conductivities (ksat) of compacted clay specimens containing macropores. A prediction
equation is formulated to predict the rate of flow of an arbitrary shape and tortuous macropore
using the flow rate of straight vertical cylinder. The predicted ksat using the proposed formulation
and calculated ksat using the LBM matched very well.

Copyright by
DURAISAMY SOUNDARARAJAN SARAVANATHIIBAN
2014

To my Parents
Soundararajan & Rajeswary

v

ACKNOWLEDGEMENTS

I would like to thank my major advisor Dr. Milind V. Khire for his continuous invaluable
guidance and encouragement over the last four years of my study at Michigan State University
which paved the way for the successful completion of this dissertation. Thank you Dr. Khire for
the opportunity to work on a variety of exciting projects, I learned a lot from them.
I would like to thank my co-advisor Dr. M. Emin Kutay for his mentoring and guidance
without which I would not have been able to complete this dissertation. Dr. Kutay, thank you for
your continuous support and patience.
I would like to thank my dissertation committee members Dr. Phanikumar Mantha and
Dr. Wei Zhang for reading my dissertation and their valuable comments. My special thanks to
the members of my PhD qualifying exam committee: Dr. Karim Chatti and Dr. Farhad Jaberi.
My sincere thanks to the National Science Foundation (NSF) for financial support of this
project (Grant No. CMMI-1100020). Also, I would like to thank the Environmental Research
and Education Foundation (EREF) and Waste Management, Inc. for their support for the field
test section and laboratory tests.
I acknowledge the Department of Civil and Environmental Engineering, the College of
Engineering, and the Graduate School of Michigan State University for the financial support
provided to me as teaching assistantships, research fellowships, and travel grants.
I would like to thank Ms. Pat Bartling, Dr. Liwang Ma, and Dr. Laj Ahuja, USDA
Agricultural Research Service for their technical support and suggestions related to RZWQM.

vi

I would like to thank Craig Burck for his technical assistance for the laboratory work.
And, I would like to thank Lori Larner, Laura Taylor, Mary Mroz, and Margaret Conner for their
assistance in processing all the paperwork on time.
I would like to thank my colleagues for their support and friendship.
Finally, I would like to thank my family and friends for all of their love and support.
Their love, encouragement, and understanding have made this dissertation possible.

vii

TABLE OF CONTENTS

LIST OF TABLES .......................................................................................................................... x
LIST OF FIGURES ....................................................................................................................... xi
KEY TO ABBREVIATIONS ...................................................................................................... xvi
CHAPTER 1: INTRODUCTION ................................................................................................... 1
1.1 Background ........................................................................................................................... 1
1.2 Objectives .............................................................................................................................. 4
1.3 Outline of the Dissertation .................................................................................................... 4
CHAPTER 2: LITERATURE REVIEW AND BACKGROUND ................................................. 6
2.1 Earthen Landfill Covers ........................................................................................................ 7
2.2 Preferential Flow ................................................................................................................... 8
2.3 Field Studies on Macropore Flow through Earthen Landfill Covers .................................. 10
2.4 Macropore Flow Models ..................................................................................................... 12
2.5 Root Zone Water Quality Model ......................................................................................... 14
2.6 Lattice Boltzmann Model .................................................................................................... 15
2.7 Synthesis of Previous Work and Motivation for Current Study ......................................... 19
CHAPTER 3: FIELD OBSERVATION OF MACROPORE FLOW .......................................... 22
3.1 Field-Scale Test Section ...................................................................................................... 22
3.2 Temperature Specific Calibration of Water Content Sensors ............................................. 25
3.3 Field versus Laboratory Measured ksat ................................................................................ 32
3.4 Soil Water Storage (SWS)................................................................................................... 34
3.5 Percolation ........................................................................................................................... 36
3.6 Irrigation Tests .................................................................................................................... 41
3.7 Effective Field Hydraulic Conductivity .............................................................................. 46
3.8 Gas Injection Tests and Identification of Macropores ........................................................ 49
CHAPTER 4: WATER BALANCE MODELING ....................................................................... 56
4.1 UNSAT-H Model ................................................................................................................ 56
4.2 Root Zone Water Quality Model (RZWQM) ...................................................................... 57
4.2.1 Macropore Flow Component ........................................................................................ 60
4.3 Input Parameters .................................................................................................................. 62
4.3.1 Field Measured Unsaturated Hydraulic Properties ....................................................... 62
4.3.2 Vegetative Data ............................................................................................................ 72
4.3.3 Meteorological and Hydrologic Data ........................................................................... 76
4.3.4 Initial and Boundary Condition .................................................................................... 77
4.3.5 Macropore Flow Data ................................................................................................... 77
4.4 Evaluation of RZWQM for Water Balance Modeling of Earthen Landfill Cover .............. 78
4.4.1 Simulations of Irrigation Tests ..................................................................................... 78
viii

4.4.2 Simulations of Year 1 ................................................................................................... 83
4.4.3 Simulations of Year 2 ................................................................................................... 87
4.4.4 Simulations of Year 3 ................................................................................................... 91
4.5 Macropore Flow Simulations using RZWQM .................................................................... 92
4.5.1 Calibration of RZWQM for Macropore Flow .............................................................. 95
4.5.2 Year 2 Simulation with Macropore Flow ..................................................................... 98
4.5.3 Year 3 Simulation with Macropore Flow ................................................................... 102
CHAPTER 5: LABORATORY TESTING OF MACROPOROUS COMPACTED CLAY
SAMPLES................................................................................................................................... 105
5.1 Self-healing of Cracks in Compacted Clay Soil ................................................................ 106
5.2 Macroporous Specimen Preparation ................................................................................. 109
5.3 Measurement of Macropore Flow ..................................................................................... 116
5.4 High-Resolution X-Ray CT Images .................................................................................. 119
5.5 Macroporosity and Characterization of Macropores ......................................................... 122
CHAPTER 6: LATTICE BOLTZMANN MODELS ................................................................. 124
6.1 Saturated (single phase) LBM ........................................................................................... 124
6.1.1 Real-numbered Solid Density Micropore and Macropore Flow Model ..................... 127
6.2 Theoretical Verification of the D2Q9 and the D3Q19 Model .......................................... 128
6.3 Laboratory Validation of the LBM ................................................................................... 133
6.4 Effect of Morphology and Tortuosity of Macropore on Macropore Flow ........................ 140
6.4.1 Effect of Macropore Shape and Shape Parameters .................................................... 145
6.4.2 Effect of Tortuosity .................................................................................................... 152
6.5 Development of a Simplified Predictive Equation for Macropore Flow .......................... 156
6.6 Verification of the Simplified Predictive Equation for Macropore Flow ......................... 157
6.7 Multiphase LBM ............................................................................................................... 161
6.7.1 Shan-Chen (SC) Multiphase Model ........................................................................... 161
6.7.2 Color Indices (CI) Model ........................................................................................... 168
CHAPTER 7: SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS FOR FURTHER
STUDIES .................................................................................................................................... 171
7.1 Summary ........................................................................................................................... 171
7.1.1 Field Data ................................................................................................................... 172
7.1.2 Water Balance Modeling ............................................................................................ 173
7.1.3 Formation of macropore in Laboratory Specimen ..................................................... 175
7.1.4 Lattice Boltzmann Model (LBM) ............................................................................... 175
7.2 Conclusions ....................................................................................................................... 177
7.3 Recommendations for Further Studies .............................................................................. 179
REFERENCES ........................................................................................................................... 181

ix

LIST OF TABLES

Table 4.1: Input parameters for water balance modeling using UNSAT-H and RZWQM .......... 63
Table 4.2: Saturated and unsaturated hydraulic input parameters used for modeling. ................. 65
Table 4.3: In-situ saturated and unsaturated hydraulic properties. ............................................... 65
Table 5.1: Measured ksat of macroporous specimens .................................................................. 118
Table 5.2: Characteristics of formed macropore and estimated macroporosity ......................... 123
Table 6.1: Macropore parameters investigated ........................................................................... 144
Table 6.2: Characteristics of macropore, and calculated and predicted ksat ............................... 160

x

LIST OF FIGURES

Figure 1.1: Schematic of flow through an earthen cover containing only micropores (a); versus
micropores and macropores (b). ......................................................................................... 3
Figure 3.1: Cross section of the compacted clay field test section (not to scale) ......................... 24
Figure 3.2: Photograph of WCR ................................................................................................... 27
Figure 3.3: Drying of compacted soil sample in a temperature chamber for WCR calibration. .. 27
Figure 3.4: Variation of WCR period with soil temperature and VWC: Topsoil (a); and Storage
layer (b). ............................................................................................................................ 30
Figure 3.5: VWC versus WCR period and soil temperature for storage layer (a), and topsoil (b).
........................................................................................................................................... 31
Figure 3.6: Comparison of measured VWC with and without temperature specific calibration of
WCR. ................................................................................................................................ 32
Figure 3.7: Compacted clay sample (a); and flexible wall permeameter (b). ............................... 33
Figure 3.8: Soil water storage of storage layer ............................................................................. 35
Figure 3.9: Measured percolation and precipitation recorded at site vs. data collected by NOAA
weather station .................................................................................................................. 38
Figure 3.10: Historical precipitation (a); and recorded solar radiation and air temperature (b). .. 39
Figure 3.11: Irrigation setup and the wetted area. ........................................................................ 42
Figure 3.12: Comparison of measured percolation during irrigation tests. .................................. 43
Figure 3.13: Soil water content of cover, cumulative precipitation and irrigation during irrigation
Test-1 (a); and Test-2 (b). ................................................................................................. 45
Figure 3.14: Estimated keff and degree of saturation (S) of the storage layer. .............................. 48
Figure 3.15: Schematic of the field test section with flux chambers and gas probes. .................. 50
Figure 3.16: Plan view of the field test section showing the flux chamber, gas probes, and
observed macropore locations........................................................................................... 51
Figure 3.17: Measured chamber fluxes from October 12, 2012 to January 4, 2013 .................... 55
xi

Figure 4.1: Schematic showing conceptual model for RZWQM ................................................. 58
Figure 4.2: Co-located water content and water potential sensors ............................................... 66
Figure 4.3: Estimated field-based SWCCs (a); and measured and predicted unsaturated hydraulic
conductivities (b) for storage layer. .................................................................................. 70
Figure 4.4: Sample collected from the test section for root-length density measurement. ........... 73
Figure 4.5: Measurement of root-length density function. ........................................................... 74
Figure 4.6: Measured normalized root-length density function.................................................... 75
Figure 4.7: Measured and simulated cumulative percolation (a); and soil water storage (b) during
the first irrigation test (2010). ........................................................................................... 80
Figure 4.8: Simulated cumulative surface runoff during the first irrigation test (2010)............... 81
Figure 4.9: Measured and simulated cumulative percolation (a); and soil water storage (b) during
the second irrigation test (2013). ...................................................................................... 82
Figure 4.10: Simulated surface runoff during the second irrigation test (2013). .......................... 83
Figure 4.11: Measured and simulated cumulative percolation (a); and soil water storage (b)
during Year 1. ................................................................................................................... 85
Figure 4.12: Simulated cumulative surface runoff (a); and cumulative evaporation (E) and
cumulative potential evaporation (PE) (b) during Year 1. ................................................ 86
Figure 4.13: Measured and simulated cumulative percolation (a); and soil water storage (b)
during Year 2. ................................................................................................................... 89
Figure 4.14: Simulated cumulative evaporation (E) and cumulative potential evaporation (ET)
(a); and cumulative surface runoff (b) during Year 2. ...................................................... 90
Figure 4.15: Measured and simulated cumulative percolation (a); and soil water storage (b)
during Year 3. ................................................................................................................... 93
Figure 4.16: Simulated cumulative surface runoff (a); and cumulative evaporation (E) and
cumulative potential evaporation (ET) (b) during Year 3. ............................................... 94
Figure 4.17: Measured and simulated percolation (a); and soil water storage (b) during the
second irrigation test with micropores and macropores. .................................................. 96
Figure 4.18: Simulated surface runoff during the second irrigation test with and without
macropores. ....................................................................................................................... 98
xii

Figure 4.19: Measured and simulated percolation (a); and soil water storage (b) during Year 2
with and without macropores. ......................................................................................... 100
Figure 4.20: Simulated cumulative surface runoff during Year 2 with and without macropores.
......................................................................................................................................... 101
Figure 4.21: Measured and simulated cumulative percolation (a); and soil water storage (b)
during Year 3 with and without macropores................................................................... 103
Figure 4.22: Simulated cumulative surface runoff during Year 3 with and without macropores.
......................................................................................................................................... 104
Figure 5.1: Particle size distribution (a); and compaction curve (b) of the clay used for
macroporous specimens. ................................................................................................. 107
Figure 5.2: Cracked (a); and self-healed (b) sample-1; and cracked (c); and self-healed (d)
sample-2. ......................................................................................................................... 108
Figure 5.3: Measured water content of Sample 2 during drying and wetting............................. 110
Figure 5.4: Macroporous specimen preparation: at the end of the first lift (a); and at the end of
compaction (b). ............................................................................................................... 111
Figure 5.5: Measured hydraulic conductivity of compacted clay intact samples with time (a); and
with number of pore volume (b). .................................................................................... 113
Figure 5.6: Ratio of outflow to inflow (a); and cumulative volume change of the specimens (b)
during flexible wall permeability test. ............................................................................ 114
Figure 5.7: Measured hydraulic conductivity of macroporous specimen using flexible wall
permeameter. ................................................................................................................... 115
Figure 5.8: Double ring rigid wall permeameter (DRRP) setup: bottom plate (a); specimen inside
the rigid wall cell (b); bottom of specimen (c); and test setup (d). ................................. 117
Figure 5.9: Macroporous specimen before test (a); and after test (b). ........................................ 118
Figure 5.10: X-ray CT scan image: original greyscale image (a); binary image (b); and variation
of threshold along the Line A-A’ (c). ............................................................................. 120
Figure 5.11: 3D structure of macropores: S1 (a); S2 (b); and S3 (c). ......................................... 121
Figure 5.12: Frequency distribution of equivalent radius of the formed macropore. ................. 123
Figure 6.1: Microscopic velocity directions of D2Q9 model (a); and D3Q19 model (b). ......... 125

xiii

Figure 6.2: Verification of D2Q9 permeable solid model using Poiseuille flow simulation
between two parallel plates. ............................................................................................ 129
Figure 6.3: Verification of D3Q19 model using Poiseuille flow simulation through a cylindrical
pipe.................................................................................................................................. 130
Figure 6.4: Velocity profiles between two parallel plates for different ns values ....................... 131
Figure 6.5: Hydraulic conductivity curve under gravity driven flow using analytical solution by
Walsh el al. (2009) and LB model prediction. ................................................................ 133
Figure 6.6: Construction of 3-D structure of the specimen using image slices. Two image slices
of specimen S1 namely S1-035 and S1-035, and interpolation based on three repetitions
of slice S1-035. ............................................................................................................... 134
Figure 6.7: Measured and calculated ksat of macroporous specimens. ....................................... 137
Figure 6.8: Effect of hydraulic gradient and resolution on calculated ksat of macroporous
specimens using LBM..................................................................................................... 137
Figure 6.9: Effect of selected threshold on calculated ksat of macroporous specimens using LBM.
......................................................................................................................................... 138
Figure 6.10: Effect of decrease or increase in macropore width in one pixel increment. .......... 140
Figure 6.11: Sectional shapes of macropores used to generate 3D image: Circle (a); diamond (b);
ellipse (c); and rectangle (d). .......................................................................................... 141
Figure 6.12: Tortuosity parameters and different tortuosity in same size of domain ................. 141
Figure 6.13: Ratio of macropore flow rate per area to micropore flow rate per area of different
shape of macropores: MP1 (a); and MP2 (b). ................................................................. 146
Figure 6.14: Effect of macropore shape on the equivalent hydraulic conductivity ................... 148
Figure 6.15: Effect of perimeter of macropore on macropore flow rate ..................................... 149
Figure 6.16: Effect of aspect ratio of macropore on macropore flow rate .................................. 150
Figure 6.17: 3D visualization of absolute value of velocity distribution for different cross
sections. ........................................................................................................................... 151
Figure 6.18: Velocity vectors and streamlines along the X-Z plane of circular cross section. .. 153
Figure 6.19: Comparison of flow rate through tortuous path and straight inclined path. ........... 154
Figure 6.20: Influence of tortuosity on macropore flow rate. ..................................................... 155
xiv

Figure 6.21: Comparison of simulated and predicted macropore flow rate. .............................. 157
Figure 6.22: Cross section of macropore relevant to minimum equivalent radius (not to scale).
......................................................................................................................................... 159
Figure 6.23: Calculated vs. predicted and measured ksat of macroporous samples .................... 160
Figure 6.24: Phase transition simulated using single component SC multiphase model............ 163
Figure 6.25: Simulated different contact angles (color bar shows the density of fluid) ............. 166
Figure 6.26: Simulation of infiltration of liquid through macropores using SC model .............. 167
Figure 6.27: Simulation of infiltration of water through unsaturated macropores using CI model.
......................................................................................................................................... 170

xv

KEY TO ABBREVIATIONS

1-D = one-dimensional
2-D = two-dimensional
3-D = three-dimensional
ASTM = American Society for Testing and Materials
CT = Computed tomography
GCL = Geosynthetic clay liner
HDMP = Heat dissipation matric potential
LB = Lattice Boltzmann
LBM = Lattice Boltzmann model
NOAA = National Oceanic and Atmospheric Administration
MSW = Municipal solid waste
RZWQM = Root zone water quality model
SWS = Soil water storage
TDR = Time domain reflectometry
VWC = Volumetric water content
WCR = Water content reflectometer

xvi

CHAPTER 1: INTRODUCTION

1.1 Background
In the U.S., landfilling is the most common means to dispose of municipal solid waste (MSW).
Landfill final cover is constructed to cap the landfill and to isolate the MSW from the external
environment. The conventional landfill final covers recommended by Resource Conservation
and Recovery Act (RCRA) consist of a hydraulic barrier layer, typically a geosynthetic material
called geomembrane underlain by compacted clay, to impede percolation (Albright et al. 2004).
Due to economical and sustainability reasons, the use of earthen cover as an alternative to
conventional cover has been extensively tested during the last two decades (Montgomery and
Parsons, 1989; Benson and Khire, 1995; Khire et al. 1997; Melchior, 1997; Benson et al. 2001;
Albright et al. 2003; Scanlon et al. 2005). However, a few studies have suggested that earthen
covers exhibit various degrees of macropore flow (preferential flow) during the several decades
of service life (Montgomery and Parsons, 1989; Khire et al. 1999; Albright et al. 2006a, Benson
et al. 2007). Macropore flow can increases the percolation into the waste significantly. However,
very little data and modeling of macropore flow through earthen cover can be found in literature.
Commonly used numerical models to simulate percolation through earthen landfill covers
are UNSAT-H, HYDRUS, and VADOSE/W (Khire et al. 1999, Mijares and Khire 2012, Benson
et al. 2001). These models are based on classical theory of water flow as formulated by Richards’
equation for flow through saturated-unsaturated porous materials that exhibit capillarity. Several
studies have evaluated the ability of these models to predict field water balance (Khire et al
1

1997, Benson et al 2005). However, these models do not simulate flow through macropores.
Figure 1.1(a) shows a schematic of an as-built compacted clay landfill cover where the flow is
dominated by micropores. The transient liquid flow in micropores is downward or upward
depending on the hydraulic gradient due to the soil capillarity and gravity, respectively. Figure
1.1(b) shows a schematic of the landfill cover after the formation of desiccation cracks
(macropores) where liquid flow is through micropores as well as through the macropores. Flow
through macropores is predominantly downward, similar to in a conduit. Inability to model
macropore flow is one of the key limitations of commonly used water balance models. Hence,
models capable of simulating micropore flow as well as macropore flow need to be developed
and/or validated in order to design earthen landfill covers.
Root Zone Water Quality Model (RZWQM) is a two-domain water balance model
developed by Agricultural Research Service (ARS), United States Department of Agriculture
(USDA) (Ahuja et al. 1995). Two-domain models use gravity gradients for macropore flow and
Richard’s equation for micropore flow. In this 1D vertical multilayer flow model, flow through
cylindrical macropores in surface layer and cylindrical or in plane macropores in the layers
beneath the surface layer can be simulated. The RZWQM model has not been used for predicting
water balance of engineered earthen covers and validation of its macropore flow component has
not been done using field-based data for engineered earthen covers.
For saturated flow, a numerical model based on the Lattice Boltzmann method (Succi 2001,
Wolf-Gladrow 2000) was used. The lattice Boltzmann model (LBM) requires high resolution
pore structure geometry that is generally developed from an X-ray Computed Tomography (CT)
imaging. The LBM was studied and validated to measure the saturated hydraulic conductivity of
porous media by many researchers (Succi et al. 1989; Kutay et al. 2006; Sukop et al. 2013).
2

However, published works are generally limited to impermeable solid medium with micropores
or macropores.

Figure 1.1: Schematic of flow through an earthen cover containing only micropores (a);
versus micropores and macropores (b).

3

1.2 Objectives
The key objectives of this study are:


Quantify micropore flow vs. macropore flow through an engineered earthen landfill final
cover;



Identify and validate numerical model(s) capable of simulating micropore flow and
macropore flow through earthen landfill covers; and



Develop a model based on lattice Boltzmann method that can simulate saturated flow
through micropores and macropores.

1.3 Outline of the Dissertation
This dissertation is organized as seven chapters including this first chapter of introduction.
Chapter 2 presents the literature review and background on earthen landfill covers,
macropore flow through earthen landfill covers, models capable of simulating macropore flow,
and the relatively new numerical approach of LBM.
Chapter 3 presents the observations of macropore flow through a field test section of earthen
landfill cover. Data collected from the field test section is presented. Also, field tests such as
irrigation tests to evaluate macropore flow and gas injection tests used to identify the macropores
are discussed in detail.
Chapter 4 presents the validation of RZWQM for water balance modeling of the earthen
landfill cover using the data presented in Chapter 2. Results of micropore only flow simulations

4

of RZWQM were compared with the commonly used UNSAT-H model simulation results.
Results of micropore and macropore flow using RZWQM are presented
Chapter 5 presents the laboratory tests conducted to evaluate the macropore flow. The
challenges in developing macroporous specimen due to the self-healing nature of cracks in the
lab-scale compacted clay specimens are discussed. A new technique to develop macroporous
specimen in laboratory is presented. Also, calculation of macroporosity and prediction of ksat
using high resolution X-ray CT images are presented.
Chapter 6 presents the model developed based on LBM. Verification of LBM algorithm
using analytical solution(s) is presented. Also, validation of LBM using X-ray CT images of
macroporous specimens and saturated hydraulic conductivities measured in the laboratory are
presented. A systematic analysis of morphology and tortuosity of macropores on macropore flow
rate is presented. Moreover, formulation of an equation to predict flow rate through arbitrary
shape and tortuous macropores using the flow rate of equivalent straight and vertical cylindrical
macropore is presented. In addition, process of development of multiphase LBM is presented.
Chapter 7 presents the summary and conclusions.

5

CHAPTER 2: LITERATURE REVIEW AND BACKGROUND

According to the estimate of the United States Environmental Protection Agency (USEPA), there
are about 1,750 active landfills in the U.S. as of 2006 and over 10,000 closed municipal landfills
(EPA 2007). Landfill final covers are constructed to isolate the waste from the environment
once the landfill reaches its capacity. The final cover recommended by Resource Conservation
and Recovery Act (RCRA) consist of a hydraulic barrier layer to restrict percolation, and are
generally called as conventional covers (Albright et al. 2004). The barrier layer of a conventional
cover may compose of compacted clay having a relatively low saturated hydraulic conductivity
(ksat) or a composite barrier consisting of a geosynthetic material called geomembrane underlain
by compacted clay depending on the landfill liner below the waste (EPA 1992). Due to
economical and sustainability reasons, the use of earthen covers as an alternative to
conventional cover, has been extensively tested during the last two decades (Montgomery and
Parsons, 1989; Benson and Khire, 1995; Khire et al. 1997; Melchior, 1997; Benson et al. 2001;
Albright et al. 2003; Scanlon et al. 2005). Hauser et al. (2001) estimated a potential savings of
50% on landfill cover construction costs by using earthen covers (without geomembrane) instead
of conventional cover. However, according to RCRA, percolation through the designed earthen
landfill covers should not be more than the percolation rate through the conventional cover.
Lysimeters are used to measure percolation in short term and numerical models are used to
predict the long term percolation through earthen landfill cover.

6

2.1 Earthen Landfill Covers
The conventional landfill final cover system required by RCRA is relatively expensive and the
long term performance of the compacted clay barrier especially to erosion protection is unknown
(Albright et al. 2004). Hence, alternative final covers are widely tested. Mostly the alternative
final covers rely on the water-holding capability of soils and the fact that a relatively large
percentage of precipitation returns to the atmosphere through evapotranspiration (Albright et al.
2004). Hence, the term “evapotranspirative (ET) covers” or “water balance covers” are often
used to refer to these alternative final covers. Also, the landfill final covers without the
geosynthetic materials are generally referred to as earthen landfill covers.
The earthen landfill covers are often classified into two types based on how they function:
(1) monolithic covers; and (2) capillary barriers (Khire et al. 2000). Monolithic covers have
relatively simple cross-sectional profile consisting of 150 to 300 mm-thick topsoil layer
underlain by a storage layer of earthen material having adequate water storage capacity. The
function of topsoil layer is to act as a medium for vegetative growth and erosion protection. The
thickness of the storage layer is designed according to required storage for infiltration until the
stored water is removed by evapotranspiration. Finer-grained soil layer is underlain by coarsergrained layer in capillary barriers. The difference in hydraulic properties between the two layers
forms a capillary interface that forces water to be stored in the upper finer-grained layer thus
minimizing infiltration into the underlying coarser-grained layer until the upper layer is almost
saturated.
Earthen landfill covers required to provide performance that is equivalent to conventional
covers with regard to reduction in percolation and for erosion protection (Albright et al. 2004;
7

USEPA 1992). Lysimeters are built in the field to quantify percolation through the landfill cover.
Nests of sensors that measure soil suction and volumetric water content are installed along the
depth of the proposed alternative earthen cover to monitor changes in the soil water storage and
meteorological data is collected to assess and model the long-term percolation from earthen
covers before they are permitted (Mijares 2011). Several field-scale demonstration projects have
been carried out across the United States to assess the performance of earthen covers for a wide
range of geo-climatic conditions.
Albright et al. (2004) carried out a majority of the field demonstration projects under the
Alternative Cover Assessment Program (ACAP) of the USEPA. Field data from such
demonstration projects is collected over a period of only few years. However, if the earthen
cover is permitted for permanent use at the site, it is required that it perform satisfactory beyond
the landfill’s 30-year post-closure monitoring period. Montgomery and Parsons (1989), Benson
and Khire (1995), and Khire et al. (1997) presented field data on macropore flow through earthen
landfill cover. More detailed literature review on macropore flow through earthen covers is
presented in Section 2.3.

2.2 Preferential Flow
Preferential flow includes funnel flow, finger flow, and macropore flow where water flows in
paths bypassing other capillary fraction of the soil matrix (Gerke, 2006). The funnel flow is
lateral redirection of infiltrating water and funneling because of the presents of textural
boundaries in soil where water flows through least resistance and bypasses the less permeable
8

zones (Kung 1990, Gerke 2006). Finger flow results due to spatial heterogeneity in soil moisture
content or air entrapment that causes instability in the wetting front leading to the formation of
thin finger-like flow patterns (Hardie et al. 2011, Glass et al. 1988). Water bypass through
relatively large continuous opening (macropores) in the field such as desiccation cracks due to
drying and wetting, cracks due to freezing and thawing, vegetation root penetrations, inter-clod
voids, rodent holes, and worm holes is called as macropore flow. Among these three kinds of
preferential flow, macropore flow in the storage layer is the key concern related to performance
of engineered earthen landfill covers because macropore flow could significantly increase
percolation into the underlying waste (Song and Yanful 2010).
Beven and Germann (1982) presented an overview of the relevant processes and conditions
that cause the macropore flow in soils. Initiation and maintenance of macropore flow in the
macropore system require a supply of water exceeding all losses to the soil matrix (Beven and
Germann 1982). This most commonly occurs at soil surface when rainfall intensity exceeds
infiltration capacity (Beven and Germann 1982). Macropore flow reduces the surface runoff and
increases total infiltration (Workman and Skaggs, 1990). Interaction between the domains, the
macropore and the matrix, depends on the supply of water within the macropores and the
hydraulic conditions within the soil matrix (Beven and Germann 1982). Voids of “equivalent
cylindrical diameter” larger than about 0.3 - 0.5 mm are considered as macropores (Jarvis 2007).
Daniel (1984) reported 5 to 10,000 times higher ksat in the field than measured in the
laboratory for compacted fine grained soils. Larger clods size and hence the larger inter-clod
voids in field results in higher saturated hydraulic conductivity in the field. Benson and Daniel
(1990) studied the influence of clods on hydraulic conductivity of compacted clay and found that
9

the clod size had large influence on the ksat of compacted clay, especially when compacted dry of
optimum. Lower dry unit weight is associated with soils having large inter-clod voids. Soils that
do not have larger inter-clod voids have higher dry unit weight (Benson and Daniel 1990).
Benson and Daniel (1990) reported compacted samples with relatively large clods have up to six
orders of magnitude higher ksat compared to the compacted samples that contain smaller size
inter-clod voids. Phifer et al. (1994) and Albrecht and Benson (2001) presented laboratory test
results that demonstrate an increase in ksat of compacted clay samples when subjected to wet-dry
cycles. Kim and Daniel (1992), Othman and Benson (1994), and Benson and Othman (1993)
showed changes in hydraulic conductivity of compacted clay due to freeze-thaw cycles.

2.3 Field Studies on Macropore Flow through Earthen Landfill Covers
Montgomery and Parsons (1989) presented three years of data of two compacted clay covers of
field test section located near Milwaukee, WI (US) that had same clay layer thickness (1.22 m,
CL USGS) and different topsoil thicknesses (0.15 m or 0.46 m). The shrinkage limit of the clay
soil was in between 10 to 12. The percolation was relatively high during the 2nd and 3rd year
compared to the 1st year. They noticed formation of cracks after the drought season.
Observations from test pits excavated after two years showed occasional crack size from 6 to 12
mm wide and 1.3 to 1.4 m deep from the surface, and continuous root penetration up to 0.7 m.
Moreover, they found root penetration along the crack plane up to 1.2 m from the surface.
Montgomery and Parsons (1989) concluded that the increase in percolation is due to the
formation of macropores and the thickness of the topsoil layer did not affect the percolation rate.

10

Melchior (1997) monitored two compacted clay covers at a site near Hamburg, Germany,
for 8 years. Both the test sections had identical profiles, 0.6 m thick clay layer overlain by 0.25 m
sand drainage layer and 0.75 cm surface layer, and different slopes, 4% and 20%. A 10 m × 50 m
lysimeter was used to measure percolation. The measured percolation for the first year was
between 1.9 – 8.4 mm and the percolation increased to between 150 – 201 mm during the 8th
year. Observations from test pits excavated 7 years after the construction showed that the clay
layer was cracked and invaded with plant roots through the entire depth. Tracer test carried out
by Melchior (1997) confirmed the preferential flow paths.
Khire et al. (1997) published data for two field test sections of compacted clay covers one in
Atlanta (humid climate) and another one in East Wenatchee, Washington (semi-arid climate).
The test sections were 30×30 m2 in areal extent and 12.2 ×18.3 m2 pan-type lysimeter was used
to measure the water balance. The test sections consists of 0.9 m of compacted clay (ML-CL
USCS) overlain by 0.15 m thick topsoil at site in Atlanta and 0.6 m of compacted clay (MH
USCS) overlain by 0.15 m thick topsoil at the site in Wenatchee. The measured annual
percolation at Wenatchee was about 7.5 mm during the 1st year and it increased to about 22 mm
during the 3rd year. During the 3rd year, the measured cumulative percolation showed step like
increase. Test pit excavated within the test sections showed cracks at the Wenatchee site and no
sign of cracks at Atlanta site. Also, they observed breakthrough while the water content at the
deepest water content sensor remained relatively low and constant. Hence, Khire et al. (1997)
concluded that the percolation measured at Wenatchee was possibly influenced by macropore
flow.

11

Albright et al. (2006a) studied the field performance of a compacted clay landfill cover
located in Albany, Georgia (humid climate) for 4 years. A 10 × 20 m2 pan-type instrumented
lysimeter was used to measure the water balance. The test section consisted of 0.45 m of
compacted clay (SC USCS) overlain by 0.15 m of uncompacted clay (SC USCS). They observed
increase in percolation about eight months after the construction following a 7-week drought.
Measured cumulative percolation during the first seven months of the year was 44 mm and the
measured cumulative percolation increased to 230 mm during the next five months. They noticed
desiccation cracks at the end of drought season. Also, the pattern of percolation changed from
steady and relatively independent of timing of precipitation before the drought to rapid and
closely related to timing of precipitation events after the drought. Moreover the hydraulic
conductivity from laboratory tests, field tests, and from percolation measurements before and
after the drought showed an increase by about three orders of magnitude over the service life of
the cover. Dye tracer test and soil structure analysis confirmed the presence of the cracks. They
concluded that the formation of macropores due to desiccation or freeze-thaw cycling resulted in
the preferential flow.

2.4 Macropore Flow Models
Water balance models which can simulate macropore flow fall under two categories: dualpermeability/porosity models and two-domain models. The dual permeability/porosity models
use the Richards’ equation for the macropore region. HYDRUS-1D is one of the dual
permeability/porosity models (Simunek et al. 1998, 2008; Simunek and van Genuchten 2008).
One of the limitations of the dual permeability models is that the models simulate upward flow in
12

the macropores during evapotranspiration (Kohne et al. 2009). On the other hand, two-domain
models implement gravity flow in macropore flow and Richard’s equation for micropore flow.
RZWQM is a two-domain water balance model developed by Agricultural Research Service
(ARS), United States Department of Agriculture (USDA) (Ahuja et al. 1995). A detailed
literature review of RZWQM is presented in Section 2.5.
A relatively new numerical technique used to model fluid flow through macropores in two or
three dimensional is LBM. The LBM was studied and validated to measure the saturated
hydraulic conductivity of porous media by many researchers (Succi et al. 1989, Kutay et al.
2006, Sukop et al. 2013). Also, LBM is capable of simulating interfaces between fluids without
the need of interface tracking (Sukop and Or 2004). Three approaches are at the research level
for multiphase LBM (Sukop and Or 2004). Gunstensen et al. (1991) presented a multiphase
model where two immiscible fluids were given different colors (red and blue) and the collision
rules were modified to obtain the surface tension between two fluids. Shan and Chen (1993)
proposed a multiphase and multicomponent model which applies interactions similar to van der
Waals attractions between fluid particles. Another approach called “Free energy” was presented
by Swift et al. (1996). However, the free energy approach has been criticized as a “top-down”
approach where macroscopic properties such as the surface tension are supplied directly to the
model (Chin et al. 2002, Sukop and Or 2004). A detailed literature review on LBM is presented
in Section 2.6.

13

2.5 Root Zone Water Quality Model
The USDA-Agricultural Research Services developed RZWQM that includes a macropore flow
model along with a water balance model (DeCoursey and Rojas, 1990). Ahuja et al. (1993)
carried out numerical analysis on the influence of macropore flow and chemical transport by
systematically changing the rainfall sequence, macropore radius, and initial soil water content,
and using three chemicals for a constant macroporosity. They considered a 1.5 m deep soil
profile of a silty clay loam. They found that the macropore size (diameter 2 mm vs. 0.25 mm)
has only a small effect for a fixed macropore volume fraction (0.05%), hence one needs to only
know the overall maximum macropore flow rate than the macropore size distribution in order to
characterize the field macropore flow. However, Ahuja et al. (1993) did not compare the model
results with field or lab data. Ahuja et al. (1995) evaluated and refined the macropore component
using a column experiment (15 cm diameter and 30 cm long) with a 3 mm diameter artificial
vertical hole.
Kumar et al. (1998) carried out field evaluation of macropore flow component RZWQM
using data collected for three years from subsurface drains located at 1.2 m deep and 28.5 m
apart in corn plots near Nashua, Iowa. They used the first year data to calibrate the model and the
following two years data to validate the model. Saturated and unsaturated hydraulic properties
were estimated using one of the approaches given in RZWQM. They interpreted the
macroporosity from infiltration experiments. However, they calibrated the saturated hydraulic
conductivity of surface layer and lateral saturated hydraulic conductivity using the data. Also,
they assumed the diameter of the cylindrical macropore. They concluded that a single value of
macroporosity for all years may not produce accurate results.
14

Malone et al. (2001, 2004) evaluated the macropore component of RZWQM for chemical
transport using undisturbed block samples (30 × 30 × 30 cm3) and pan lysimeters (61 × 61 cm2),
respectively. They calibrated the soil parameters related to hydrology using all the measured data
(eg. saturated hydraulic conductivity) and tested only the chemical sub-component of macropore
flow. Malone et al. (2001) concluded that chemical transport through macropore flow was
sensitive to effective macroporosity, macropore radius, and effective soil radius. Malone et al.
(2004) concluded that the RZWQM can be used to simulate chemical (metribuzin) transport in
the field when macropore flow and hydrology are accurately simulated.

2.6 Lattice Boltzmann Model
Succi et al. (1989) demonstrated the adherence of the LBM to Darcy’s law for a complex 3D
porous media. They used artificial porous media and imposed no-slip boundary condition for
solids (impermeable solids). Cancelliere et al. (1990) compared the numerical (LBM) analysis of
flow through artificial porous media developed using sphere solid shape with permeability
theories. Their LBM implementation was similar to Succi et al. (1989). They concluded that at
solid fractions less than 0.2, the simulated permeabilities match with the predictions of
Brinkman's effective-medium theory, whereas at higher solid fractions, good match is achieved
with Kozeny-Carman equation.
Shan and Chen (1993) presented LBM (herein called SC model), which has the capability of
simulating multiphase and multicomponent immiscible fluids with constant temperature. Shan
and Chen (1994) discussed the application of SC model to single component multiphase (SCMP)
15

fluid simulations. Martys and Chen (1996) studied the flow of a multicomponent fluid subject to
gravitational and surface tension forces in complex 3D porous media using SC model.
Dardis and McCloskey (1998) and Freed (1998) independently proposed an alternative
approach to model Darcy flow in LBM in which the permeability of the media is related to a
momentum sink (damping). They found that assigning solid porosity parameter (ns) at each node
to represent the effect of structure on measured fluid flow is good enough to produce results that
match well with predictions. This approach allows avoiding the requirement of a detailed
micropore structure (geometry) when modeling flow through micropores.
Spaid and Phelan (1997) published a modification to standard LBM equivalent to solving a
Brinkman formulation to simulate fluid flow through macropores and micropores. The Brinkman
equation (Brinkman 1947) is a generalized Darcy’s equation that allows the matching of
boundary conditions at an interface between macropores and micropores (Martys 2001). Martys
(2001) presented an improvement to the Brinkman equation approximation. The Brinkman’s
equation (Brinkman 1947) is as follows,

p   e  2 u 


K

u

Eq. 1

where P is the pressure; u is the velocity; K is the intrinsic permeability; and  e is an effective
viscosity parameter (a parameter that allows for matching of the shear stress boundary condition
across the free-fluid (macropore/porous medium interface). Since the hydraulic conductivity k is
k  K  , where  is the unit weight, the equation Eq. 1 can be rewritten as follows,



p   e  2 u  u
k
16

Eq. 2

It should be noted that the first part in the right side of Eq. 1 is actually the Stokes’ equation,
which is solved by the LBM formulations, whereas the second part is the Darcy’s equation. In
macropore space, since k is infinite, second part goes zero and the Stokes’ equation holds.
Whereas in the permeable medium, since the velocity is small, the first part negligibly small
value and second part dominates (Darcy’s law). The second component can be implemented in
LBM as a body force term ( ak ) applied to the permeable solid regions that causes momentum
sink as suggested by Martys (2001) (Eq. 3)

ak  


g
u u
k
k

Eq. 3

where  is the density of the fluid; and g is the gravitational acceleration. This body force term
is added to body force term in standard LBM. In addition to adding the body term, the relaxation
time to be used in the permeable medium is modified as,

e  3

e
 0.5


Eq. 4

where  e is the modified relaxation time assigned for lattice notes in permeable medium.
Kutay et al. (2006) conducted saturated fluid flow in granular materials using 2D and 3D
LBM. They used 3D geometries of compacted aggregates and asphalt specimens generated from
X-ray Computed Tomography (CT) technique as input for the LBM. The compacted samples
were 10 cm diameter and 13 cm height. They verified the accuracy of the LBM by comparing the
results with analytical solutions of simple geometries and hydraulic conductivity measurements
on the compacted aggregates and hot mix asphalt specimens. They concluded that the results of

17

LBM simulations were in excellent agreement with those obtained from analytical calculations
and laboratory measurements.
Sukop and Or (2004) demonstrated the ability of SC single component multiphase model to
simulate multiphase flow in slit and simple angular pores. Pan et al. (2004) studied the SC
multicomponent model to simulate immiscible fluid flow through sphere packs. They found
encouraging agreement as compared to experimental data for two-phase displacement process,
where one fluid (non-wetting) enters the system and replaces the other fluid (wetting).
Korner et al. (2005) presented a multiphase LBM for free surface modeling. In this model,
the density distribution function evolves only in the liquid region. That is the Boltzmann
algorithm is applied only in the liquid region. Hence, capturing of the interface is necessary for
the free surface model. Also, the free surface model cannot be used to study the liquid-liquid or
liquid-vapor systems where the individual phases effect each other (Attar and Korner, 2009). The
model is suitable only for those liquid-gas systems where the gas phase has negligible effect on
liquid phase (Attar and Korner, 2009). The model requires keeping track of the mass exchange
between cells. The interface advection is identified from the mass fraction variable that is
updated by recording the inflow and outflow of mass using the distribution function. Also, this
particular model requires an additional boundary condition at the interface cells because the
distribution functions coming from gas cells after advection are unknown. Owing to the extreme
complexity in the implementation of this approach, this model is not considered in this study.
Misici and Palpacelli (2005) studied the water percolation through coffee using multiphase
LBM based on an innovative methodology that includes ‘color indices’ for different fluid phases.
More detail on Misici and Palpacelli (2005) model is presented in Section 6.5. Huang and Lu
18

(2009) analyzed the relative permeability of heterogeneous porous media through numerical
study using single component SC model.
Sukop et al. (2013) studied the saturated non-Darcy flow in macroporous limestone aquifer
samples using 3D LBM. They used 3D images of limestone aquifer samples from X-ray CT as
input for the LBM. The domain size of 9×9×9 cm3 was chose such that the irregular boundary of
samples lies inside the domain. The simple no-slip wall boundary condition was used for solid
rock boundaries. They observed reduction in the apparent hydraulic conductivity as the applied
gradient increased to levels representative of those observed in the field due to non-Darcy flow.
Newmen and Yin (2013) presented the effect of geometry on the flow transition (from Darcy
flow to non-Darcy flow) using LBM. They used synthetic 2D porous media for their simulations.
Again they used the no-slip boundary condition for the solids. They found that geometries that
are microscopically heterogeneous experience early change from Darcy flow to non-Darcy flow
systems due to the contrast in the size of pore bodies and the size of pore throats.

2.7 Synthesis of Previous Work and Motivation for Current Study
The field data presented by Montgomery and Parsons (1989), Melchior (1997), Khire et al.
(1997), and Albright et al. (2006) suggest that macropore flow influences percolation in earthen
landfill covers during service. However, a systematic quantitative evaluation of macropore flow
through earthen landfill cover needs to be carried out. Also, there is a lack of published work on
modeling the macropore flow through earthen landfill covers. Khire et al. (1997) and Khire and
Saravanathiiban (2013) used UNSAT-H to predict percolation through earthen landfill covers
19

that had macropores. However, they pointed out that the model does not simulate flow through
macropores and hence underestimates percolation. Hence, there is a need to develop model
capable of simulating flow through micropore as well as macropore to design earthen landfill
cover that meet long-term performance criteria.
The two-domain model RZWQM is capable of simulating gravity flow through cylindrical
macropores. Ahuja et al. (1993) carried out numerical analysis using macropore flow component
of RZWQM. Ahuja et al. (1995), Kumar et al. (1998), Malone et al. (2001), and Malone et al.
(2004) conducted laboratory or field evaluation of the macropore component of RZWQM.
However, previous studies on validation of macropore flow module of RWQM have primarily
concentrated on chemical transport. The hydraulic properties of soils were not measured. The
hydraulic properties of soils were calibrated based on model predictions. There is no published
study that validates RZWQM for water balance of engineered earthen landfill covers. Hence,
validation of RZWQM for water balance modeling of engineered earthen landfill cover and
macropore flow module of the model using field measured data is required.
LBM for fluid flow through porous media has been validated by researchers (Succi et al.
1989, Cancelliere et al. (1990), Kutay et al. 2006, and Sukop et al. 2013). However, published
work on flow through macropores of porous media using LBM is limited to impermeable solids.
Although permeable solids (micropores) algorithms are published, an extended work on porous
media with micropore flow (i.e., flow through the ‘solid’ areas) and macropore flow has not been
done. The major reason for the lack of interest in permeable solids is probably the fact that flow
through micropores is negligibly small during saturated flow as compared to the flow through
macropores. However, micropore flow may dominate during unsaturated flow through soils.
20

Hence a comprehensive physically based model to simulate flow through micropore and
macropore and a systematic analysis is the first step towards modeling unsaturated flow through
soils that has macropores.

21

CHAPTER 3: FIELD OBSERVATION OF MACROPORE FLOW

3.1 Field-Scale Test Section
A compacted clay final cover test section with instrumentation was constructed in late October
2009 at a MSW landfill located in Detroit, Michigan (Mijares 2011). The climate at the site is
classified as sub-humid and the average annual precipitation is approximately 83 cm (32 in). The
annual potential evapotranspiration (PET) to precipitation ratio is about 1.3. A typical cross
section passing through the test section is shown in Figure 3.1. Percolation through the test
section was measured using a pan lysimeter. The lysimeter consisted of a trapezoidal excavation
of the test section having areal dimensions of about 8.5 m x 8.5 m and 0.6 m deep. The lower
boundary of the lysimeter was impermeable (Mijares et al. 2012). The lysimeter was filled with
clean pea gravel. The hydraulic conductivity of the pea gravel was equal to 1 cm/s.
The upper boundary of the gravel layer was vented to the atmosphere using a perforated
HDPE pipe to maintain atmospheric pressure on the lysimeter. The storage layer of the test
section was constructed using a 1.5 m thick native glacial till (clay) which was overlain by a 0.3
m thick topsoil layer. The liquid limit (LL) and the plasticity index (PI) of the clay are 29 and 13,
respectively. The shrinkage limit (SL) of the clay as predicted using the LL and the PI based on a
procedure presented in Holtz and Kovacs (2010) is 13.25. Percentage of fines in the clay soil is
about 78%. The clay soil and topsoil were classified as CL and ML, respectively according to the
Unified Soil Classification System (USCS). The optimum water content and maximum dry unit
weight of the clay for Standard Proctor (SP) effort was 13.4% and 18.5 kN/m3, respectively. The
22

clay was compacted at 4% dry of optimum water content in 30 cm thick lifts using a sheepsfoot
compactor. The average dry unit weight was about 98% of the maximum dry unit weight.
Water content reflectometer (WCR) (Model: CS616, Manufacturer: Campbell Scientific,
Inc.) and heat dissipation-based matric potential (HDMP) sensors (Model: CS431, Manufacturer:
Campbell Scientific, Inc.) were installed at six vertical locations in three nests within the test
section to measure the volumetric water content (VWC) and matric suction, respectively. A
detailed temperature depended calibration was carried out to develop the calibration curve for
WCRs. The ground freeze and thaw was monitored using a thermistor sensor inserted into the
topsoil about 10 cm below the ground surface. Pressure transducers (Model CS431) having
measurement range of 0 to 90 cm with 0.1% accuracy were used to monitor the level of water in
the lysimeter to measure percolation. Calibration curves were developed to estimate the volume
of water collected in the lysimeter and to estimate percolation from the measured water level in
the lysimeter (Mijares 2011). The site had a weather station to record air temperature,
precipitation, solar radiation, barometric pressure, and relative humidity.
The saturated hydraulic conductivity of the storage layer was measured immediately after the
construction of the compacted clay test section, that is before the test section experience any wetdry cycles, using a 1 m by 1 m single-ring infiltrometer (ASTM 2010b, Guyonnet et al. 2000)
and the average value of the saturated hydraulic conductivity was 4.4 × 10-6 cm/s (Mijares et al.
2012). More details pertaining to the field test section and instrumentation can be found in
Mijares (2011).

23

Sprinkler
for irrigation

Top sensor
nest

Bottom
sensor nest
0.3 m thick topsoil
Lysimeter
sensor nest
Subsurface
lateral
drainage
trench

0.3 m

1.5 m thick
compacted
clay

Vent

Subsurface
lateral
drainage
trench

Water content
and
water
potential
sensors

Lysimeter
filled with pea gravel

0.6 m
Liquid level sensor

Geomembrane + GCL

8.5 m
Waste (MSW)
30 m

Figure 3.1: Cross section of the compacted clay field test section (not to scale)

24

3.2 Temperature Specific Calibration of Water Content Sensors
Time domain reflectometry (TDR) is the most common method used for measuring in situ VWC
non-destructively (Benson and Wang, 2006). However, conventional TDR (GHz frequency) can
be expensive due to required specialized electronic interfaces or multiplexers (Udawatta and
Anderson, 2011; Benson and Wang, 2006). Newer TDR system such as water content
reflectometer (WCR) includes all of the electronics in the head of the probe and hence easy to
use compared to conventional GHz frequency TDR. However, the WCR operates at lower
frequency (MHz range) that makes the water content measurement more sensitive to factors like
soil type or clay content, temperature, and soil electrical conductivity (Benson and Wang, 2006).
Also, the dielectric permittivity, which is used as an indirect measurement of soil water content,
varies with the temperature (Udawatta and Anderson, 2011). Hence, a detailed soil, compaction,
and temperature specific calibration was carried out for the WCRs to measure in situ water
content.
Figure 3.2 shows a photograph of a CS616 WCR. The WCR consists of two metallic rods,
each 30 cm long and 0.32 cm diameter with 3.2 cm spacing, connected to a probe head. Voltage
pulse is generated inside the probe head and propagates along the rods. The pulse travels back to
the probe head once it is reflected at the end of the rod and triggers the next pulse. The travel
time varies depending on dielectric permittivity of the soil in which the rods are inserted. The
data logger records the reflection as period after the reflections are divided by a scaling factor
(Udawatta and Anderson, 2011; Benson and Wang, 2006). With a proper calibration equation,
the measured period in microseconds is calibrated to read volumetric water content.

25

Soil samples collected from the field were compacted in a 5 gallon bucket with the
dimension of 30 cm diameter and 36 cm height. In order to saturate the sample from bottom up,
geocomposite (geotextile and geonet) was placed at the bottom of the bucket and connected to
two tygon tubes along two diametrically opposite edges of the geonet. Field compaction criteria
were used to compact the soil in 5 gallons bucket. For the storage layer clay soil, sample was
compacted in five lifts. Each lift was about 6.6 cm thick. Achieved dry density and molding
water content were 18.05 kN/m3 and 9.4%, respectively. Density of the field topsoil was
measured using Shelby tube samples collected from the field test section and the average dry
density was 14.45 kN/m3. Topsoil sample was compacted in two lifts. Sample was gradually
saturated and the weight of the sample was recorded continuously. The topsoil reached about
95% saturation and the clay soil reached about 93% saturation. Once the sample was as saturated
as it can be, a WCR and a CS 650 sensor, about 10 cm apart, were inserted in the soil vertically.
The CS 650 was used to record soil temperature. The TDRs were connected to data logger. The
sample was sealed and placed in an environmental chamber. The sample in the environmental
chamber with TDRs connected to the data logger is shown in Figure 3.3. Initially the chamber
was set to room temperature, and WCR period and temperature of the sample was monitored.
Once the period and temperature reached steady state, the water content of the sample, period,
and temperature were recorded. Then the chamber temperature was changed by ±10oC. The
temperature was varied from 5oC to 35 oC and at each temperature the period was recorded at
steady state. Then the sample was dried to next target water content. Target water contents were
set in an interval of 0.05 VWC. WCR periods were recorded with varied temperature at each
target water contents.

26

WCR rod
WCR head

Figure 3.2: Photograph of WCR

Figure 3.3: Drying of compacted soil sample in a temperature chamber for WCR calibration.

27

Variation of WCR period with temperature for constant water contents are shown in Figure
3.4. The period changes linearly with temperature. Clay (Figure 3.4a) shows slightly more
sensitivity to temperature compared to topsoil (Figure 3.4b). At saturated water content, the
period of topsoil increased from 31.94 μs at 5oC to 36.13 μs at 36oC whereas, the period of clay
increased from 32.07 μs at 5oC to 37.14 μs at 34oC. Clay content of the soil influences the
sensitivity of WCR to temperature. Figure 3.5 shows the surfaces fitted to the calibration data.
Topsoil data was fitted using a 3rd order polynomial equation (Figure 3.5(a)) and the clay was
fitted using a 2nd order polynomial equation (Figure 3.5b). The calibration equations are as
follows:
Topsoil:

VWC T ,   2.56  5.4e 3T  0.26  1.6e 4T 2  6.7e 4T  8.5e 3 2  2.5e 7T 3  5.8e 6T 2
 1.98e 5T 2  1.0e 4 3
Eq. 3.1
Where, T is the soil temperature; and τ is the period in μs.
Clay (storage layer):

VWC T ,   0.333  0.042  0.0013 2  0.01T  5.6e 4T  5.9e 5T 2

Eq. 3.2

Eq. 3.3 presents calibration equation for the clay without the effect of temperature (Mijares
2011),

VWC    0.615  0.0689  0.00255 2  3.86e 5 3

Eq.3.3

Calibration equations were used to measure VWC of the field test section from measured
WCR periods. Figure 3.6 presents the VWC measured using WCR at 45 cm depth of lysimeter
28

nest using soil specific calibration equation (Eq. 3.3), and soil and temperature specific
calibration equation (Eq. 3.2) along with the temperature measured using a co-located HDMP
sensor. During winter when the temperature goes below 10oC, water content is underestimated
by about 40% to 50% and during the summer when temperature goes above 23oC, water content
is overestimated by about 10% if temperature specific calibration is ignored (Eq. 3.3). The water
content measurement is not impacted by temperature when the temperature is around 20oC.

29

40

(a) Topsoil

CS616 Period (s)

35

VWC
0.397
0.326
0.244
0.156
0.084
0.040

30

25

20

15
0

10

20

30

40

50

40

50

o

Soil Temperature ( C)
40

(b) Storage layer

CS616 Period (s)

35

30

VWC

25

0.315
0.257
0.158
0.131
0.105
0.073
0.061

20

15
0

10

20

30
o

Soil Temperature ( C)

Figure 3.4: Variation of WCR period with soil temperature and VWC: Topsoil (a); and Storage
layer (b).

30

0.6

(a)
(a)

0.5

3

Volumetric Water Content (cm /cm )

0.6

3

0.5
0.4

0.4
0.3

0.3

0.2
0.1

0.2

0
36
34

40

32
30

20

28

CS616 (s)

0.1

30
10

26

Temperature (C)

(b)
(b)

0.5

0.45

3

Volumetric Water Content (cm /cm )

0.6
0.4

3

0.5
0.4

0.35

0.3

0.3

0.2

0.25

0.1

0.2

0
0.15
-0.1
0.1

35
45
40

30

35

0.05

30
25

25

20

0

15

CS616 (s)

20

10
5

Temperature (C)

Figure 3.5: VWC versus WCR period and soil temperature for storage layer (a), and topsoil (b).

31

0.45

30

0.35

25

Soil temperature

20

0.30

15

0.25

10
VWC
(soil
specific
calibation)

0.20

o

3

0.40

Soil Temperature ( C)

VWC
(soil and
temperature
specific
calibation)

3

Volumetric Water Content, VWC (cm /cm )

45 cm depth of lysimeter nest

5

0.15
Nov/1/2011

Aug/1/2011

May/1/2011

Feb/1/2011

Nov/1/2010

Aug/1/2010

May/1/2010

Feb/1/2010

Nov/1/2009

0

Figure 3.6: Comparison of measured VWC with and without temperature specific calibration of
WCR.

3.3 Field versus Laboratory Measured ksat
ksat was measured using the falling head and rising tail flexible wall permeameter test (ASTM
2010a). The soil was prepared by passing through size 4.8 mm (No. 4) sieve. The sample was
compacted at Standard Proctor effort in 10.2 cm diameter and 11.6 cm tall compaction mold
according to ASTM (2007). However, an additional 5 blows per lift were required (i.e., 30 blows
32

per lift) to achieve the field dry unit weight. The achieved dry unit weight of the compacted
sample was 17.9 kN/m3. The compacted sample is shown in Figure 3.7a and the sample during
flexible wall test is shown in Figure 3.7b. Confining pressure of 68.8 kPa was applied during the
flexible wall tests. The influent and effluent pressures were maintained at 55.0 kPa and 13.8 kPa,
respectively. The test was continued for about a week until the ratio between inflow and outflow
reached unity and the measured ksat reached a steady-state.

(a)

(b)

Figure 3.7: Compacted clay sample (a); and flexible wall permeameter (b).

The laboratory measured ksat was 6.6 x 10-8 cm/s. The field ksat measured using single-ring
infiltrometer immediately after the test section construction was 4.4 × 10-6 cm/s (Mijares et al.
2012) which is one order in magnitude greater compared to the ksat measured using the flexible
wall permeameter. The effect of clod size was not well captured in the lab-scale test is most
likely the reason for this difference (Benson and Daniel 1990). In addition, the influence of

33

higher effective stress in the lab versus what the soil experienced in the field may be another
reason for lower ksat measured by flexible wall permeameter (Trast and Benson 1995).

3.4 Soil Water Storage (SWS)
SWSs of the storage layer were calculated by integrating the water contents measured at three
sensor nest locations using WCRs (using temperature specific calibration equation). The SWSs
for the top nest, lysimeter nest, and bottom nest are presented in Figure 3.8. Also, degree of
saturation (S) of the storage layer is presented on the right vertical axis of Figure 3.8. The peak S
(field saturation) for the bottom, the lysimeter, and the top nests were about 90%, 88%, and 82%,
respectively. The SWS of the bottom nest is always higher than the other two nests due to lateral
flow. However, lysimeter nest shows more drying than top nest during dry season. The
percolation collected in lysimeter is removed from the system in addition to evapotranspiration.
On the other hand, for the top nest and the bottom nest percolation entered into the waste can be
pulled back into the storage layer during drying. Hence, the top next holds more water than the
lysimeter nest. The influence of bottom boundary on water balance of earthen cover is presented
in Mijares and Khire (2012).

34

Irrigation
125 mm
(on Lysimeter)

Bottom Nest

100
94
88

50
82
76

top nest
40

70
64
Lysimeter Nest
58

30
52

Year 1

Year 3

Year 2

46

Year 4

Figure 3.8: Soil water storage of storage layer

35

May/15/2013

Nov/11/2012

May/11/2012

Nov/8/2011

May/8/2011

Nov/4/2010

May/4/2010

20

Degree of Saturation, S (%)

Irrigation
208 mm
(on Lysimeter)

Nov/1/2009

Soil Water Storage, SWS (cm)

60

Compared to the first year, during the second and the third years, higher SWSs were noticed
in all three nests. Increase in precipitation during the second year could be the reason for the
increase in the SWSs for the second year. Precipitation recorded during the second year was 1,
104 mm and it represents a relatively wet year with 98 percentile precipitation. However, unlike
the SWSs of the lysimeter and the top nests, SWS of the bottom nest showed relatively steady
value for long duration during the first year indicating field saturation. And, during the second
year the SWS of the bottom nest further increased and stayed steady. The sudden increase in
SWS of the bottom nest during the second year could be because of the formation of macropores
due to wet-dry cycles or animal activity. Macropore flow increases infiltration and decreases the
surface runoff which leads to higher SWS (Workman and Skaggs 1990). Because the SWSs of
the top and the lysimeter nests did not show the indication of field saturation during the first
year, the increase in SWSs during the second and the third year cannot be regarded as the
influence of macropore flow based on the SWSs alone. Also, water content measurements are
point measurements which do not capture the flow in between sensors. Unless macropore
penetrate through the WCR rods, the WCRs will not capture the rapid macropore flow. Hence,
SWS alone cannot be used to identify the presence of macropore flow in landfill earthen covers.

3.5 Percolation
Measured percolation from the lysimeter and recorded precipitation are presented in Figure 3.9.
Precipitation recorded at the site did not include snow because the tipping bucket at the site did
not have heater to melt the snow. Hence, precipitation recorded by National Oceanic and
Atmospheric Administration (NOAA) weather station located within 5 km of the site is plotted in
36

Figure 3.9 for comparison. Recorded precipitation at the site underestimates the precipitation by
about 8% to 18%. During the first year, precipitation was about 792 mm, while cumulative
percolation was 7 mm. During the second year, the precipitation was 1, 104 mm, while
cumulative percolation was about 75 mm. And, during the third year, the precipitation was 828
mm, while percolation was about 136 mm. The reason for about an order increase in percolation
during the second year compared to the first year is partly due to record high precipitation during
the second year. Annual precipitation of the NOAA weather station during the last 50 years is
presented in Figure 3.10(a) that showed the second year (2010) was the wettest year in the 50
year history. Precipitation received during the first year (2009) is equal to about 35 percentile.
During the second year (2010), precipitation received at the site corresponds to about 98
percentile. Also, the precipitation received during winter and spring for the second year was 482
mm (44% of total precipitation) whereas the winter and spring precipitation for the 1st year was
only 241 mm (30% of total precipitation). The solar radiation is lower during winter and spring
season compare to summer which results in less evapotranspiration (ET). Recorded solar
radiation and recorded average temperature at the site is presented in Figure 3.10(b). During
winter solar radiation drops below 500 MJ/m2/day and during summer solar radiation is as high
as 18,000 MJ/m2/day. Also, ET is limited by presence of snow (Albright et al. 2010,
Apiwantragoon 2007). Moreover, presence of snow reduces surface runoff hence increases
infiltration which results in higher percolation (Albright et al. 2010). In addition, snowmelts
directly applied to the cover.

37

350

3500
Year 3

Year 4

828 mm (48 percentile)

2800

1104 mm (98 percentile)

210

2100
~ 136 mm

Precipitation:
NOAA
Site

140

1400
Percolation

~ 75 mm

792 mm (35 percentile)

70

700
~ 7 mm

May/1/2013

Oct/30/2012

May/1/2012

Oct/31/2011

May/2/2011

Oct/31/2010

0
May/2/2010

0

Figure 3.9: Measured percolation and precipitation recorded at site vs. data collected by NOAA weather station

38

Cumulative Precipitation (mm)

Year 2

280

Nov/1/2009

Cumulative Percolation (mm)

Year 1

140

(a)

2011 - 121 cm

Precipitation (cm)

120

2010 - 82 cm

Average:
83 cm

100
80

2012 - 69 cm

60
40

0

1960
1963
1966
1969
1972
1975
1978
1981
1984
1987
1990
1993
1996
1999
2002
2005
2008
2011

20

Year

40
Solar
Radiation

Temperature

3000

-10

0

-20

0

0

May/1/13

6000

Oct/30/12

10

May/1/12

9000

Oct/31/11

20

May/2/11

12000

Oct/31/10

30

May/2/10

15000

Daily Average Temperature ( C)

(b)

Nov/1/09

2

Recorded Solar Radiation (MJ/m /day)

18000

Figure 3.10: Historical precipitation (a); and recorded solar radiation and air temperature (b).
39

Measured cumulative percolation during the third year was about 20 times that of the first
year and about two times that of the second year cumulative percolation. However, precipitation
received during the third year was only 42 percentile and winter and spring precipitation was 480
mm (58% of total precipitation) which is about the same as the second year. The peak solar
radiation during the third year was below 15,000 MJ/m2/day whereas peak solar radiation during
the 1st and 3rd year was as high as 18,000 MJ/m2/day. Also, average solar radiation for the 3rd
year was about 6,000 MJ/m2/day which is less than average solar radiation of 7,500 MJ/m2/day
and 6,750 MJ/m2/day for the 1st and the 2nd year, respectively. Lower solar radiation in the 3rd
year resulted in less ET in the 3rd year. Hence, compared to the 1st and the 2nd years, the test
section was at field saturated condition (about 86%) for relatively longer period of time during
winter and spring of the 3rd year (Figure 3.8). Hence, the increase in percolation during the third
year could be partly because test section stayed close to saturation for a relative longer period of
time.
Moreover, macropores such as shrinkage cracks and animal barrows were noticed on the
surface of the test section during the second year and after. The PI and SL of the clay suggest
medium level of susceptibility for desiccation cracking (Holtz and Kovacs 2010). These
macropores might have contributed to higher percolation observed in the 2 nd and the 3rd years in
addition to the reasons discussed previously. Gradual slope of the percolation vs. time plot is
typical of micropore flow. Whereas, sharp increase in the rate of flow is indicative of macropore
flow (Khire et al. 1997).

Hence, cumulative percolation plot indicated the presence of

macropore flow better than the SWS of the test section. However, only the use of this data
(increase in total precipitation, winter and spring precipitation, and the SWSs) makes it harder to
derive conclusions on macropore flow.
40

3.6 Irrigation Tests
Irrigation tests were carried out on the test section to induce breakthrough and estimate fieldscale near saturation hydraulic conductivity. Centrifugal pump was used to pump water from a
storage reservoir into a rotating vane sprinkler jet installed on a 2.2 m tall fence post inserted
within the lysimeter area (Figure 3.1). The pumping rate was about 7 lpm (112 gph). Figure 3.11
shows a photo of the field setup. The white flags in Figure 3 indicate the wetted area which was
circular having a diameter of about 8.4 m. Two irrigation events were carried out on 8 June and
10 June, for 9 hours and 12 hours, respectively (Mijares et al. 2012). The surface application rate
achieved was about 1 cm/h. During the two irrigation days, the site also received precipitation
equal to 0.75 mm and 4.83 mm on 8 June and 9 June, respectively. The SWS of the cover
estimated by integrating the water contents measured by the sensors located within the lysimeter
nest (Figure 3.8) before the irrigation was 48.0 cm on 8 June and 48.1 cm on 10 June. The SWSs
correspond to about 81% of S of the storage layer.
About 3.5 years after the construction of the test section and about 3 years after the first
irrigation test, the second irrigation test was carried out in May 2013. Similar to the first
irrigation test, water was pumped from the storage reservoir. The pumping rate achieved was
about 7.8 lpm (124 gph). Wetted area due to the irrigation was circular having an average
diameter of 10.7 m. Two irrigation events were carried out, on 3 May and 4 May, for 13 and 11
hours, respectively. The surface application rate achieved was about 0.52 cm/h. The SWS of the
test section on 3 May before the irrigation was 48.6 cm (82% of S) and it increased to 49.3 cm
(83% of S) on 4 May before the second irrigation application.

41

Sprinkler

Wetted Area

Figure 3.11: Irrigation setup and the wetted area.

Figure 3.12 presents cumulative percolation following the irrigation events.

While no

distinct breakthrough was observed within 45 days after the first irrigation test, after 45 days,
following a major precipitation event (~ 75 mm), significant percolation (~ 2 mm) was collected
during a relatively short period of time immediately after the storm even that occurred in late
July 2010. The S of storage layer was about 84% during the breakthrough. The slope of the
percolation plot during the 45 day period after the first irrigation event is relatively gradual and it
increases sharply following the storm event. The gradual increase is typical of capillary flow
through micropores. Whereas, the sharp change in the rate of flow is indicative of macropore
flow (Khire et al. 1997). Hydraulic conductivity of the clay was estimated by assuming unit
gradient through the cover when the breakthrough occurred in late July.
During the second irrigation test, the breakthrough occurred within 4.5 hours after the
irrigation started and continued for a period of about 4 days (Figure 3.12). Total 10.5 mm of
42

percolation was collected during this relatively short period. The average degree of saturation of
the cover was about the same as that during the first irrigation test (~ 84%). The rate of surface
application of water during the second irrigation test was less than that during the first test. Thus,
the relatively rapid breakthrough and relatively high magnitude of percolation are indicative of
formation of macropores in the clay cover during the service period. The estimated hydraulic
conductivities of the storage layer during breakthrough events are discussed in Section 3.6.

Cumulative Percolation (mm)

12

10

8

Test 2 (May 2013)
-6

Effective Field Hydraulic Conductivity (K ) ~ 9 × 10 cm/s
eff

6

Breakthrough within 4.5 hours
Test 1 (June 2010)

4

-6

Effective Field Hydraulic Conductivity ~ 1 × 10 cm/s
Breakthrough after 45 days

2

0

0

10

20

30

40

50

60

70

80

Days after irrigation events

Figure 3.12: Comparison of measured percolation during irrigation tests.

Figure 3.13 present the water contents before, during, and after the two irrigation tests.
Cumulative precipitation and irrigation and cumulative percolation are also presented on the
same plots to assess the effect of irrigation on the flow through the cover. In order for the plot to

43

be legible, water content measurements from all six depths within the cover are not shown in
Figure 6. During the third year of service, the WCR at the 15 cm depth went bad hence the water
contents at 15 cm depth was predicted using measured suction from co-located potential sensor
(HDMP) and in-situ soil water characteristic curve (SWCC) developed using sensor data. More
details on in-site SWCC are presented in Section 4.3.1 and Saravanathiiban and Khire (2014a).
The estimated water content is presented in Figure 3.13b.
During the first irrigation test, water content at 15 cm depth reached saturation while water
contents at 75 cm and 165 cm depths showed insignificant increase indicating only a little
portion of water applied infiltrate into storage layer. For about 45 days after the first irrigation
test, the water contents at 15 cm and 75 cm depth gradually decreased (test section was drying)
and the water content at 165 cm depth gradually increased as small amount of percolation was
collected by the lysimeter. Following a major precipitation event (~ 75 mm) in late July 2010,
the water content at 15 cm and 75 cm depths rapidly increased and significant percolation (~ 2
mm) was collected during a relatively short period of time. During the second irrigation test,
water contents at 75 cm depth showed slight increase, and water contents at 165 cm depth did not
change (Figure 3.13) while significant breakthrough occurred.

44

60

(a) Test 1 (Year 2010)

Depth:

0.35

45

75 cm

0.3

30
15 cm
Irrigation
15

0.25
165 cm

Aug/12

Jul/27

Jul/11

Jun/25

(b) Test 2 (Year 2013)

3

0

30
Depth:
15 cm

3

Irrigation

24

0.35

18

75 cm
0.3

12

Precipitaion
165 cm

0.25

6

May/14

Apr/26

Apr/18

Apr/9

0.2

Apr/1

Percolation

May/5

Volumetric Water Content (cm /cm )

0.4

Jun/9

0.2

May/24

Percolation

Cumulative Precipitation and Irrigation
or Percolation (cm)

3

Drying

Cumulative Precipitation and Irrigation
or Percolation (cm)

Precipitation

3

Volumetric Water Content (cm /cm )

0.4

0

Figure 3.13: Soil water content of cover, cumulative precipitation and irrigation during irrigation
Test-1 (a); and Test-2 (b).

45

Measured onset of percolation did not show any correlation with precipitation during the 1st
year except one event after about 5 weeks of drying (July 2010). However, water content at 15
cm depth increased after every major (> 3 mm/d) precipitation events. This observation suggests
that the storage layer had not significantly structurally changed since construction and the
hydraulic conductivity of the cover was close to as-built.

On the other hand, significant

percolation can be noticed immediately after every major storm event during the third year.
Moreover, percolation occurred within few hours after precipitation events. The rapid
percolation onsets after precipitation events suggest that the macropores were formed and were
relatively deep. After precipitation events, water content at 75 cm depth showed gradual increase
before dry period and it rose rapidly after the dry period (Figure 3.13). However, water content at
165 cm depth showed gradual increase during the 1st and the 3rd years. The gradual increase in
water contents at certain depths and rapid increase at other arbitrary depths suggest that the
macropore network was not densely connected. Because distanced macropores would allow
water flow bypassing the water content sensors at some depths and would flow through the
sensors at some other depths. Densely connected macropore network would increase the water
content at all depths (Albright et al. 2006a). Hence, rapid percolation events and rapid increase in
water contents at only few depths suggests relatively less dense but deep enough macropores in
the cover.

3.7 Effective Field Hydraulic Conductivity
Effective field hydraulic conductivity (keff) was estimated using measured percolation assuming
unit hydraulic gradient and steady-state flow. Similar approach was followed by Meerdink et al.
46

(1996) and Albright et al. (2006a, 2006b). Though the assumptions made in the calculation of keff
are simplistic, keff represents the in-service field hydraulic conductivity of the storage layer
(Albright et al. 2006b). In the absence of ponding, field ksat is expected to be greater than keff
(Albright et al. 2006b). If keff is estimated during similar field saturation level, then keff should be
relatively similar. Hence, keff can be used as an indicator to identify the change in field ksat.
keff of the test section was estimated when the breakthrough occurred after the irrigation
events. For the first irrigation test, the breakthrough occurred in late July 2010 was used to
estimate the keff. Estimated keff was 1 x 10-6 cm/s. Relatively steady percolation measured for 2
hours was considered for keff estimation. This keff corresponds to S of 84% based on the water
contents measured during breakthrough. Estimated keff during the second irrigation test was 9 x
10-6 cm/s which was about two times higher than the field ksat measured using SRI. Estimated keff
from the two irrigation tests correspond to the same S (~ 84%). After 3.5 years in service, k eff of
the clay has increased by a factor of nine. The increase in keff indicates the increase in field ksat
and confirms the presences of macropore flow in the cover.
Also, keff was estimated for the entire monitoring period using percolation measured after
major storm events. Minimum of 2 hours of percolation was identified and used for k eff
estimations. During the first year, only few percolation events can be noticed. During the second
year to the fourth year, frequent percolation breakthrough occurred after every major
precipitation event. Estimated keff and S of the storage layer are presented in Figure 3.14. The S
was mostly between 80-87% during the breakthroughs. However, keff showed an increasing trend
from the second year of service onwards. During the first year, average estimated keff was about
9 x 10-7 cm/s. The maximum keff of 4.3 x 10-5 cm/s was estimated during the second year. About
47

62 mm of rainfall was recorded within 8 days previous to the maximum keff event and the test
section was at its maximum S of 87.5%. Note that the second year was the wettest year in the 50
year history as mentioned previously. Except the two higher values of keff (above 2 x 10-5 cm/s)
during the second year, keff follows an increasing pattern from about 9 x 10-7 cm/s to about 1 x
10-5 cm/s. Hence, keff of the test section increased by at least an order of magnitude in about four
year of service due to the formation of macropores.
-4

90
85
80

10

S

-5

k

75
eff

70

10

65

-6

60
55
Year 1

Year 2

Year 3

Year 4

May/31/13

Nov/25/12

May/22/12

Nov/17/11

May/15/11

Nov/9/10

May/6/10

-7

Nov/1/09

10

50

Figure 3.14: Estimated keff and degree of saturation (S) of the storage layer.

48

Degree of Saturation, S (%)

S = 80% - 87%

eff

Effective Hydraulic Conductivity, k (cm/s.)

10

3.8 Gas Injection Tests and Identification of Macropores
Field gas injections test was carried out to evaluate the methane (CH4) oxidation capability of the
compacted clay cover and the gas test data was used to locate, confirm, and characterize
macropores. CH4 and carbon dioxide (CO2) in about 1:1 ratio by volume, were mixed and
injected into the lysimeter for about five months starting from 13 September 2012. Flow meters
with regulators were used to monitor and maintain the inflow rate of gas. The average flow rate
of CH4 and CO2 injected were about 660 g/day and 1100 g/day, respectively. Flux chambers and
gas probes were used to measure CH4 flux and the gas concentrations within the cover at various
depths. Figure 3.15 and Figure 3.16 present the cross sectional view and plan view, respectively,
of the test section with gas injection point, flux chambers, and gas probes. Flux chamber
measurements were conducted on top surface as well as side slopes. Nine flux chambers were
used to monitor CH4 flux and these chambers were moved to new locations from time to time to
increase the flux measurement locations. Flux measurements were carried out at total 17
locations. The flux chambers used were of two sizes: 0.65 m L × 0.65 m W × 0.20 m H and 1.0
m L × 1.0 m W x 0.3 m H. Collars were placed on the surface of the test section and sealed using
bentonite clay.

49

Flux
Chambers

Gas Probes

0.3 m thick topsoil

Subsurface
lateral
drainage
trench

0.3 m

1.5 m thick
Compacted
Clay

Subsurface

Gas lateral
Injection drainage
Point trench

Water content
and
water potential
sensors

Lysimeter
Filled with Pea Gravel

0.6 m
Liquid Level
Sensor

Geomembrane + GCL
8.5 m

Waste (MSW)
30 m

Figure 3.15: Schematic of the field test section with flux chambers and gas probes.

50

35 m

N
26 m

8.5 m

Gas Injection Area
(Lysimeter)

8.5 m

LEGEND
Surface Crack
Flux Chamber - Bigger

25 m
17 m

Active Macropore

Non-Active Macropore

Flux Chamber - Smaller

Gas Probes

Figure 3.16: Plan view of the field test section showing the flux chamber, gas probes, and observed macropore locations.
51

During sampling, three to six samples were collected from each chamber sequentially over a
20 – 30 minutes of period using 140 ml syringes and filled into Tedlar bags (Zefon International,
Inc.) or glass vials. Glass vials are commonly used to collect gas samples to test using gas
chromatography GC. Tedlar bags were used in this study to collect samples to test using portable
detectors. The bags allow supplying the sample in natural flow rate of the sampler pump in
detectors. 500 milliliters of Tedlar bags and 20 ml glass vials were used. The vials were tested
using GC (Gas Chromatography) to measure the concentrations of CH4 and CO2.
A dual detector portable total organic/inorganic vapor detection instrument, Toxic Vapor
Analyzer (TVA-1000B) was used for gas survey on the cover surface and to test the samples in
addition to GC. Gas survey on the ground surface of test section was carried out to locate the
macropores. The TVA uses an FID (Flame ionization detector) and a PID (photoionization
detector). FID is highly sensitive to hydrocarbons. However, FID does not differentiate different
gases rather measures bulk concentration of hydrocarbons in the sample. PID is sensitive to
aromatics, unsaturated hydrocarbons, and chlorinated hydrocarbons (Instruction Manual 2003).
Also, PID can measure some inorganic compounds (eg. ammonia, carbon disulfide, chloroform,
and ethylamine). However, PID does not detect methane. Hence, an elevated reading of FID and
no response in PID might indicate presence of methane (Instruction Manual 2003). FID and PID
were calibrated using standard gas cylinders before every field visit. Even though FID lacks the
specificity of gas chromatographs, it requires less operator proficiency, provides rapid results,
and is less expensive (Robbins et al., 1990). Duplicate samples were collected from flux
chambers to compare the concentration measured from FID and GC. About one hundred samples
were analyzed in GC and compared with TVA measurement which showed that the TVA
measurements are within ±10%. Gas surveying was carried out once every two weeks.
52

Those macropores identified on the ground surface were located and their surface
dimensions were measured. The typical shape of the macropores was elongated or cylindrical at
surface. The widths of the macropores were few mm to about 4 cm and the length of the
macropores ranged from 1 cm to 60 cm (Khire and Saravanathiiban 2013). Only one surface
crack was identified in the lysimeter area. Many elongated cracks (tension cracks) were located
towards the downslope side of the test section. While the most common reasons for the
formation of macropores in compacted clay covers is desiccation, freeze-thaw, root penetrations,
and animal or rodent burrows, the test section built at the site in this project may also have
developed tension cracks because of the slope and the location of the test section which hangs on
the side slope (1V:4H) of the landfill.
The gas surveying carried out using portable FID during CH4 gas injection tests were used to
locate unidentified macropores in the test section. Hot-spots of CH4 gas identified using portable
FID was marked as active macropores. Some of the cracks which were visible did not conduct
CH4 gas and those were identified as non-active macropores. This suggests that not all surface
cracks are deep enough or connected. Surface crack was not visible at few hot-spots indicated by
portable FID. Cracks could have formed in storage layer but those were not connected to surface
or at surface level the crack could have been filled over during rainfall or heavy wind.
Locations of flux chambers, gas probes, surface macropores, active macropores, and nonactive macropores are shown in Figure 3.16.

Five flux chambers were placed enclosing

macropores in order to measure the methane flux through macropores. Measured chamber fluxes
are shown in Figure 3.17. The average flux varied between 0.1 g/day/m 2 to 80 g/day/m2.
Micropore flux chambers (enclosed intact ground surface ) (F1, F4, F6, F8, F11, F12, F13, F15,
53

F18, and F19) had average chamber flux less than 1 g/day/m2 while surface crack chambers (F3,
F7, F9, F16, and F17) had average chamber flux between 8 to 80 g/day/m2. Flux chambers F2,
F5, F10, and F14 had flux between 3 to 10 g/day/m2. This higher flux compare to micropore flux
chambers suggests possibility of macropore presents in flux chambers F2, F5, F10, and F14.
Measured average methane fluxes through surface crack flux chambers were about two orders of
magnitude higher than average methane flux through micropore flux chambers. CH4 oxidization
capacity of cover was estimated using the procedure published in Chanton et al. (2008). The CH4
oxidization capacity of cover estimated for flux chambers enclosing macropore were zero for
three chambers and 1.8% and 3% for rest of the two chambers while CH4 oxidization capacity of
the intact portion of the cover was up to 30%. Two orders of magnitude higher methane flux and
almost zero CH4 oxidization capacity of flux chambers enclosing macropores indicates that those
macropores were well formed and reached deep enough to the lysimeter.

54

1000

4

2

CH Flux (g/m /day)

100

10

1
Observed
Surface
Crack

0.1

Macropore

0.01

F1
F4
F6
F8
F11
F12
F13
F15
F18
F19
F2
F5
F10
F14
F3
F7
F9
F16
F17

Micropore

Flux Chambers Locations

Figure 3.17: Measured chamber fluxes from October 12, 2012 to January 4, 2013

55

CHAPTER 4: WATER BALANCE MODELING

Water balance modeling was carried out using these two public domain numerical models: (1)
UNSAT-H; and (2) RZWQM.

4.1 UNSAT-H Model
UNSAT-H is a public domain, 1-D, finite difference, water flow and heat transfer model, and
widely used for earthen landfill final cover modeling (Fayer 2000, Khire et al. 1997, Khire et al.
2000, Benson et al. 2005, Bohnhoff et al. 2009, Mijares and Khire 2012). 1-D modeling is
considered adequate for modeling hydrology of earthen landfill covers (Khire et al. 1997).
Benson et al. (1994) indicated that error in the water balance by ignoring the lateral flow is
relatively small (less than 1.5%). Moreover, measured lateral flow in the field test section
confirmed that the flow is predominantly 1-D in the vertical direction (Mijares and Khire 2012).
Hence, 1-D UNSAT-H model was used to simulate percolation through the field test section.
In UNSAT-H model, infiltration is simulated in two stages. The first stage is flux-controlled,
that is the infiltration rate is controlled by water supply (Fayer 2000). During the second stage,
the infiltration rate is based on the ability of the soil profile to transmit water downward. The
infiltration rate is a function of the time from the onset of precipitation, the initial water content,
the hydraulic properties of the surface soil, and the hydraulic properties of the layers deeper
within the profile (Hillel, 1980, Fayer, 2000). Surface runoff is not explicitly simulated in
UNSAT-H, however if rainfall or irrigation rate is greater than the infiltration capacity, the
56

excess water results as surface runoff and removed from the system. UNSAT-H uses the
Richards’ equation for saturated and unsaturated water flow through soils (capillary pores). Soil
water retention can be described using commonly used mathematical functions such as the
Haverkamp function (Haverkamp et al. 1977), the van Genutchen function (van Genutchen
1980), and the Brooks-Corey function (Brooks-Corey 1964). Hysteresis in water retention
function can be included in the van Genutchen function. Options for defining the unsaturated
hydraulic conductivity functions include the Haverkamp model, the Mualem model, the Burdine
model, and the Brooks-Corey model. Initial condition is defined by soil suction. Transpiration is
simulated as a sink term in the Richard’s equation. Evaporation is computed using Fick’s law
when heat flow is modeled and evaporation is a function of potential evaporation when heat flow
is not modeled. Potential evapotranspiration is calculated using Penman type equation using
daily relative humidity, net solar radiation, wind speed, and daily minimum and maximum air
temperatures. UNSAT-H does not simulate snowmelt, freeze or thaw, and temperature
dependence of soil properties (Fayer 2000).
UNSAT-H was used to simulate percolation due to micropore flow through the test section
during the first three years of service. RZWQM predictions were compared with UNSAT-H
predictions.

4.2 Root Zone Water Quality Model (RZWQM)
RZWQM (version 2.42) is a 1-D, two-domain, finite difference, soil water flow and chemical
transport model that simulates saturated and unsaturated flow in vertical soil profiles (Ahuja et.
al. 2000). The model assumes that soil consist of two interacting domains: (1) micropore domain
57

representing micropores in the soil matrix; and (2) macropore domain representing cracks or
other larger pores where flow occurs due to gravity only. A schematic showing how RZWQM
computes the water balance is shown in Figure 4.1.
Flux Boundary Condition

Precipitation

If not raining
Evapotranspiration
(extended
Shuttleworth and Wallace model)
Infiltration
(Green-Ampt equation)

If overland flow >
If rainfall rate Macropore flow
Capacity
> Infiltration rate
Overland Flow
Runoff
Deadend
Macropores

Horizon 1: Top Soil

Continuous
Macropores
(Poiseuille’s
law)

Horizon 2: Storage Layer
Water Redistribution
(Richards’ equation)

Macropore to
Micropore
(Lateral Green-Ampt)

Unit Gradient Boundary
Percolation
(Micropore flow + Macropore flow)

Figure 4.1: Schematic showing conceptual model for RZWQM

Three grid systems are used to simulate water flow in RZWQM: (1) horizons/layers for
defining hydraulic properties, (2) layers for redistribution, and (3) 1 cm thick layers for
infiltration. The Green-Ampt approach (Green and Ampt 1911) is used to model infiltration into
the soil matrix (Eq. 4.1),
58

V  k sat

hc  H 0  Z wf
Z wf

Eq. 4.1

where V = the infiltration rate at any given time (cm/hr.); k sat = the effective saturated hydraulic
conductivity of the wetting zone (cm/hr., a correction factor of 0.5 is used to account for the
entrapped air and resulting viscous resistance);

hc = the suction head at the wetting front (cm)

estimated from unsaturated hydraulic conductivity function; Ho= the depth of surface ponding
(cm); and Zwf = the depth of the wetting front (cm).
Redistribution of water in the soil matrix following infiltration is simulated using the
Richards’ equation (Celia et al. 1990, Richards 1931) as follows,


 
h

 k (h, z )  k (h, z )  S r ( z, t )
z z 
z


Eq. 4.2

where  = volumetric water content (cm3/cm3); t = time (hr.); z = soil depth (cm); h = soil-water
pressure head (cm); k = unsaturated hydraulic conductivity (cm/hr.); and Sr = sink term for root
water uptake and tile drainage rates (hr-1). Soil hydraulic properties, SWCCs and unsaturated
hydraulic conductivity functions, for each layer is described using modified Brooks-Corey
function (Ahuja et. al. 2000). Initial condition can be defined by soil suction or water content for
each horizon.
RZWQM does not simulate surface runoff or ponding infiltration. However, the difference
between rainfall or irrigation rate and infiltration rate is resulted as overland flow when rainfall
or irrigation rate is more than infiltration rate. If macropores are present then the overland flow
re-routed into the macropores, and if the macropore flow capacity exceeds then the balance of
overland flow results as surface runoff and removed from the system. The surface boundary
59

condition is modeled as specified flux boundary similar to UNSAT-H. Evapotranspiration is
modeled as a function of potential evaporation until the soil dries to wilting point at which point
the upper boundary condition becomes constant head. Potential evapotranspiration is computed
using extended Shuttleworth and Wallace model (Farahani and Ahuja 1996). Water update by
plants is simulated as a sink term in the Richards’ equation using the approach of Nimah and
Hanks (1973). The lower boundary can be unit gradient or constant flux when water table is
present.

4.2.1 Macropore Flow Component
Overland flow generated at the soil surface flows into the macropores. The surface soil layer is
assumed to have cylindrical macropores while the layers beneath the surface layer can have
cylindrical or planar cracks. The Poiseuille law is used to calculate the macropore flow rate
capacity assuming gravity flow (unit gradient) as presented in Eq. 4.3 and Eq. 4.4,
For cylindrical holes:

N p gRp4
K

mac
8

Eq. 4.3

For planar cracks:

K

mac



Lc gw3
12

Eq. 4.4

where  = density of water (1.0017 g/cm3); g = the gravitational constant (1.27×1010 cm/hr2);
Rp = the radius of cylindrical holes (cm); w = the width of planar cracks (cm);  = the dynamic
viscosity of water (36.072 g/hr/cm); Np = the number of pores per unit area; and Lc = the total
60

length of planar cracks per unit area (cm). The number of macropores per unit area can be
calculated from total macroporosity and the radius of the single cylindrical macropore (Malone
et al. 2001),

Np 

Pmac
  R p2

Eq. 4.5

Np 

Pmac
w  Lc

Eq. 4.6

Similarly for planar cracks,

where Pmac is the macroporosity.
Continuous macropores are idealized to be vertical and extends to a point below the depth of
the soil profile modeled or to water table (if exists). However, horizontal macropores can be
defined as deadend macropores that branch off from vertical macropores in each soil horizon and
end within the soil horizons. However, the macropores are not physically modeled, only the
amount of water flow through continuous macropores, the amount of water stored in deadend
macropore, and the amount of water absorbed from continuous and deadend macropores to the
soil matrix is calculated. Macropore flow capacity limits the flow of ponded water into the
macropores. Macropore flow is sequentially routed downward through the continuous
macropores in 1 cm depth increments for each time step. During each depth increment,
macropore is allowed to absorb laterally by the soil matrix from deadend and continuous
macropores. Lateral/radial absorption rate is calculated by lateral Green-Ampt type equation. A
reduction factor is used to account for the influence of the compaction or an organic coating
surrounding the walls of macropores in lateral absorption (Malone et al. 2001). However, water
flow from the soil matrix to macropore is not considered in RZWQM.
61

4.3 Input Parameters
The input parameters used for water balance simulations using UNSAT-H and RZWQM are
summarized in Tables 4.1 and 4.2. Detail discussion on how the input parameters were measured
or estimated is presented in following subsections.

4.3.1 Field Measured Unsaturated Hydraulic Properties
In order to numerically model landfill earthen final covers, site-specific unsaturated hydraulic
properties of cover soils are required. The hydrology of earthen cover is strongly influenced by
hydraulic properties which include: saturated hydraulic conductivity, unsaturated hydraulic
conductivity function (the relationship between unsaturated hydraulic conductivities and soil
suction), and the soil water characteristics curve (SWCC) (the relationship between volumetric
water content and soil suction). Because cover soil is generally unsaturated, the unsaturated
hydraulic conductivity function and SWCC are required input for most water balance models.
Khire et al. (1995) presented the importance of using site-specific unsaturated hydraulic
conductivity function to simulate the hydrology of the earthen landfill final cover accurately. The
hydraulic properties are measured using laboratory methods or field methods. Field method is
recommended whenever possible because it is most representative of the field structure of the
soil which is often not replicated in the lab tests (Benson and Gribb, 1997).

62

Table 4.1: Input parameters for water balance modeling using UNSAT-H and RZWQM
Parameter

Input Value

Area (ha)
Elevation (m)
Aspect (deg.)
Latitude (deg.)
Longitude (deg.)
Slope (deg.)
Climate zone(rain>20 in/yr)

Soil type
Depth (cm)

Top Soil
Storage Layer
Top Soil
Storage Layer

Reference

Field Specifics
0.007 Measured
204 Google Earth
260 Google Earth
420 N Google Earth
-83.25 Google Earth
10 Measured
3 Measured
Soil Layer Description
Sandy Loam Mijares et al. 2012
Clay Loam Mijares et al. 2012
0-31 Design Values (Mijares et al. 2012)
32-182 Design Values (Mijares et al. 2012)

Model Application

RZWQM
RZWQM; UNSAT-H
RZWQM
RZWQM; UNSAT-H
RZWQM
RZWQM
RZWQM
RZWQM
RZWQM
RZWQM; UNSAT-H
RZWQM; UNSAT-H

Soil Hydraulic Properties
Brooks-Corey fitting
Top Soil
parameters
Brooks-Corey fitting
Storage Layer
parameters

Table 4.2 Lab Measured

RZWQM; UNSAT-H

Table 4.2 Saravanathiiban and Khire (2014b)

RZWQM; UNSAT-H

Meteorological and Hydrologic Data
Precipitation, air temperature, solar
radiation, wind speed
Dew point , cloud cover
Relative humidity

---

On-site measurement

-----

NOAA weather station
On-site measurement
63

RZWQM; UNSAT-H
UNSAT-H
RZWQM

Table 4.1 (cont’d)
Albedo
Albedo of dry soil
Albedo of wet soil

Initial condition
Upper boundary condition
Lower boundary condition

Total macroporosity (cc/cc)
Avg. radius of cylindrical macropores
(cm)
Fraction deadend macropores
Sorptivity factor for lateral infiltration

0.20 Chudnovskii (1966)
0.30 Ma et al. (2011)
0.20 Ma et al. (2011)
Initial and boundary conditions
Matric suction Field Measured
Water content Field Measured
Specified Flux Measured Precipitation
Unit gradient Fayer et al. (2000); Khire et al. (1997)

UNSAT-H
RZWQM
RZWQM
UNSAT-H
RZWQM
RZWQM; UNSAT-H
RZWQM; UNSAT-H

Macropore Parameters
9.50E-08 –
calibrated
1.37E-07

RZWQM

0.15 – 0.18

calibrated

RZWQM

0

Field observation

1

calibrated

RZWQM
RZWQM

64

Table 4.2: Saturated and unsaturated hydraulic input parameters used for modeling.

Model

Soil Layer

UNSAT-H

Topsoil
Storage Layer

RZWQM

Topsoil
Storage Layer

Saturated
Residual
VWC
VWC (r)
(s)
0.41
0.00
0.39
0.03

0.41
0.39

0.00
0.03

SWCC

Unsaturated conductivity function

he (cm)

B

63.1
29.7

3.33
8.3

he (cm)

A1

A2

63.1
29.7

0
0

0.30
0.12

ksat
(cm/s)
1×10-2
4.4×10-6
ksat
(cm/s)
1×10-2
4.4×10-6

hbk (cm)

B

l

63.1
12.0

3.33
-2.0

2
0.84

hbk (cm)

N1

N2

63.1
12.0

0.0
0.0

2.90
1.08

Table 4.3: In-situ saturated and unsaturated hydraulic properties.

Function

Curve

Saturated
VWC (s)

van
Genuchten

Lower
Average
Upper

0.39
0.39
0.39

Brooks-Corey

Average

0.39

Brooks-Corey

----

0.41

Note:

Storage Layer
SWCC
Residual

n
VWC (r)
cm
0.10
0.030
1.60
0.11
0.015
1.22
0.10
0.007
1.14
hb (cm)
λ
0.03
29.7
0.12
Topsoil
0.00
63.1
0.30

Original Brooks-Corey parameters

65

Unsaturated conductivity function
ksat

nk
l
(cm/s)
cm
-------0.02
1.9
-3.67
-------4.4×10-6
hbk (cm)
λk
12
-0.31
--1×10-2

63.1

0.30

----

In-situ water content and soil suction measured using co-located WCRs and HDMP sensors,
respectively, as shown in Figure 4.2 were used to develop field unsaturated hydraulic properties
(Saravanathiiban and Khire 2014b). A detail description on WCRs including the temperature
depended calibration is presented in Section 3.2.

Water potential
sensor
Length – 6 cm
Diameter – 1.5 cm
Water content
sensor
Rod length –
30 cm
Figure 4.2: Co-located water content and water potential sensors

HDMP sensor (Figure 4.2) is made up of a cylindrical-shaped porous ceramic tip with
heating element and a thermocouple inserted at its center. This sensor uses thermal conductivity
of the ceramic tip to measure soil water potential, indirectly, using a built in sensor-specific
calibration. When matric suction of the surrounding soil is in equilibrium with the ceramic tip,
soil suction is determined using the correlation between thermal conductivity and water content
of the ceramic tip that corresponds to matric suction. These sensors measure soil water potentials
66

ranging from -2,500 to -10 kPa. Because the air entry pressure of the ceramic tip is about -10
kPa, the sensors remain saturated until it reaches the air entry value. Hence, the sensor cannot
measure matric suctions between -10 kPa to 0 kPa. The sensors were calibrated individually
according to the procedure presented by Flint et al. (2002).
Field based SWCCs were developed using co-located water content and water potential
sensors. The data sets of volumetric water content (VWC) and water potential from all three
sensor nests measured at six depths in each nest were used. During the process of development
of SWCCs, a time lag was noticed between the changes in water content and corresponding
change in water potential. This time-lag is the result of unsaturated hydraulic conductivity of the
porous tip which causes a finite amount of time for the sensor potential to equilibrate with the
soil. Hence, the set of VWC and water potential data were manually chosen considering the time
lag in water potential measurements (if any) after drying events. The data set of VWC and water
potential from all three sensor nests is presented in Figure 4.3(a). The data shows some scatter.
Therefore, three SWCC curves were generated for each test section using the van Genuchten
fitting function (van Genuchten 1980) and the Brooks-Corey fitting function. The middle curve
represents the average of field properties and the other two curves show the upper and lower
bounds observed in the data. The van Genuchten fitting function is as follows,
( )

(

){

(

)

}

Eq.4.7

and
(for the Mualem model)
where

is the residual VWC;

is the saturated VWC; h is the matric suction; θ is the VWC at

suction h; and α, n, and m are the curve fitting parameters.
67

And, the SWCC curve according to Brooks-Corey function,

  s

when h < hb
Eq. 4.8

   r  hb 
 
s r  h 



when h≥ hb

where hb is the air entry pressure; and λ is the pore size distribution index. The fitting parameters
are listed in Table 4.3.
Benson and Gribb (1997) recommended the instantaneous profile method (Watson 1966,
Hillel et al. 1972, Khire et al. 1995, Meerdink et al. 1996) to measure unsaturated hydraulic
conductivities of fine-grained soils because large volume of soil can be tested. Moreover, they
recommend this method to measure the unsaturated hydraulic conductivities of landfill cover
soils because the boundary conditions are reasonably well defined when a monotonic drying is
observed in the field and percolation from the base of the cover being investigated is zero.
Hence, the instantaneous profile method is used to develop the field unsaturated hydraulic
conductivity functions using in-situ water content and soil suction. Data were selected for the
calculation based on nearly monotonic decreasing water contents and increasing matric suctions
noticed from in the data collected from the water content and water potential sensors,
respectively. Percolation was zero during the periods when the data was selected for unsaturated
hydraulic conductivity calculations. This ensured that the required key assumptions for the
instantaneous profile method were met: (1) flow is 1-D; and (2) flow is unidirectional. The
calculation was based on the assumptions that the flow was one-dimensionally orthogonal to the
test section and water was removed via evaporation only (i.e., the effect of transpiration is
negligible). Plant roots existed predominantly in the upper 30 cm thick vegetative layer and had
68

minimal penetration in the lower clay layer (more detail on plant roots is presented in Section
4.3.2). Hence, the assumption that evaporation was the key driver for the removal of water is
reasonable. A detailed step by step procedure followed in the instantaneous profile method to
calculate unsaturated hydraulic conductivity can be found in Khire et al. (1995) and Meerdink et
al. (1996). The suction values to calculate the hydraulic gradient and average suction were
inferred from the field based average SWCC (Figure 4.3). Preliminary numerical simulations of
water balance of the test section using UNSAT-H indicated that when the average SWCC is
input to the model, the water balance predictions are relatively accurate (Khire and
Saravanathiiban 2013).
Figure 4.3(b) presents the estimated unsaturated hydraulic conductivities. Also, predicted
and fitted unsaturated hydraulic conductivity curves using the van Genuchten-Mualem function
(van Genuchten 1980) and the Brooks-Corey function (Brooks and Corey 1964) are presented in
Figure 4.3(b) as a reference. The van Genuchten-Mualem conductivity function is as follows,

1  h 1  h   K
k ( h) 
2
n m

nm

1  h 

n ml

sat

Eq. 4.9

where l is the pore interaction term. Similarly, the Brooks-Corey conductivity function,

k (h)  k sat

when h < hbk
Eq. 4.10

h 
k (h)  k sat  bk 
 h 

k

when h≥ hbk

where hbk is the air entry pressure; and λk is the pore size distribution index for conductivity
function.
69

10

6

10

5

10

4

(a)

Matric Suction, h (cm)

Upper
Average
van Genutchen function
Lower
10

3

10

2

Brooks-Corey function

10

1
0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

Voumetric Water Content (VWC),  (cm /cm )

Hydraulic Conductivity (cm/sec.)

3

10

-4

10

-6

10

-8

3

(b)

10

-10

10

-12

10

-14

Measured using
Instantaneous Profile Method

Brooks-Corey (Fitted)

Brooks-Corey (Predicted)
van Genutchen-Mualem
(Predicted using l = 0.5)

1

10

100

1000

10000

10

5

10

6

Matric Suction, h (cm)

Figure 4.3: Estimated field-based SWCCs (a); and measured and predicted unsaturated hydraulic
conductivities (b) for storage layer.
70

The Brooks-Corey functions are defined in slightly different formats in UNSAT-H and
RZWQM. The SWCC curve according to UNSAT-H model:

  s

when h < hb
Eq. 4.11

   r  hb 
 
s r  h 

1

B

when h≥ hb

where B is the inverse of the original Brooks-Corey parameter λ. And, conductivity function
according to UNSAT-H model,

k (h)  k sat

when h < hbk
Eq. 4.12

h 
k (h)  k sat  bk 
 h 

2b

'

Bk

when h≥ hbk

where Bk is the inverse of λk; b’ represents 1+l for Burdine conductivity model (Burdine 1953)
where the exponent l is usually 2; and for Mualem model (Mualem 1976) b’ represents 2+l and l
is usually 0.5. The input parameters l, B, and Bk presented in Table 4.2 were calculated from the
original Brooks-Corey fitting parameters presented in Table 4.3.
Modified form of the Brooks-Corey function is used in RZWQM (Ahuja et al. 2000). Eq.
4.13 and Eq. 4.14 presents SWCC and conductivity function used in RZWQM.

   s  A1 (h)

when h < hb
Eq. 4.13

   r  B1 h A

2

when h≥ hb

where A1 and A2 are fitting parameters, and B1   s   r  A1  hb  hbA2 .

71

k (h)  k sat h

 N1

when h < hbk
Eq. 4.14

k (h)  C2 h

 N2

when h ≥ hbk

where N1 and N2 are fitting parameters (slopes of the two segments of the log (k) - log (h) curve
divided at air-entry pressure, hbk); and C2  k sat  hbkN2 . When A1 and N1 approach to zero, Eq.
4.13 and Eq. 4.14 reduce to original Brooks-Corey equations (Eq. 4.8 and Eq. 4.10) and the
parameters A2 and N2 are equal to λ and λk, respectively. The input parameters (B1, A2, and N2)
are presented in Table 4.2, which were calculated from the original Brooks-Corey fitting
parameters presented in Table 4.3.

4.3.2 Vegetative Data
Input data relevant to transpiration modeling for UNSAT-H includes root-length density
function, maximum root depth, and leaf area index (LAI), growing season, and percentage bare
area. The root-length density function was measured using samples collected from the field test
section. Two samples of 30 cm diameter and 50 cm deep were collected from two different
locations of the test section. Figure 4.4 shows a sample excavated from the test section.
The samples were cut into 5 cm slices (Figure 4.5). The plant roots were washed and the
weight of the plant roots from each slice was recorded. The average biomass at 2.5 cm and 50 cm
depths were about 2 g/m2 and 240 g/m2. The normalized root-length density was calculated by
dividing the biomass of each depth interval by total biomass for entire depth and is presented in

72

Figure 4.6. UNSAT-H uses the following equation for normalized root-length density (R)
function,
R  ae  zb  c

Eq. 4.15

where a, b, and c are the coefficients that fit the normalized root-length density data. The values
of a, b, and c for the measured data are 1.43, 0.3, 0.018, respectively. The maximum root depth
was used as 55 cm because negligible amount of roots were present at that depth.

Figure 4.4: Sample collected from the test section for root-length density measurement.

73

Healthy vegetation in the region can reach a maximum LAI of 4.0 to 4.5 according to
Schroeder et al. (1994). However, a conservative value of 3.0 was used for the maximum LAI
from visual observation. Average percentage bare area was set as 50% from visual observation
during site visits.

Figure 4.5: Measurement of root-length density function.

RZWQM has a plant growth module that simulate above and below ground biomass.
However, it requires calibration for specific plants. Hence, a simple known growth curve
approach is added to the model to simulate only water and nutrient uptake. Quick Turf module is
one of the growth curve approach modules in RZWQM that has Orchard grass built into it.
However, the growing season cannot be input to the model for Quick Turf module. The turf can
74

undergo dormancy at the beginning of the year only (Bartling et al. 2012). The turf is simulated
continuously according to planting and cutting dates. The initial plant height, plant height at cut,
initial LAI, LAI at cut, maximum root depth, and planting and cutting dates are the input
parameters relevant to water uptake for Quick Turf module. Planting and cutting dates were set
as the growth season similar to UNSAT-H model input. LAI at cut was set to 3.0 and initial LAI
was set to 0.1 to minimize the water uptake during early growing season. Preliminary simulation
results from both UNSAT-H and RZWQM showed that including plant makes relatively small
impact in simulated water balance, especially percolation. Hence, only results of water balance
modeling without plants are presented in Sections 4.4 and 4.5.

0

0.1

Normalized Root Biomass, 
rL
0.2
0.3
0.4
0.5
0.6

0.7

0

Depth, z (cm)

10

20
 = aexp(-bz)+c
rL

30

a =1.43, b=0.3, c=0.018

40

50

60

Figure 4.6: Measured normalized root-length density function.

75

0.8

4.3.3 Meteorological and Hydrologic Data
Meteorological data and soil albedo are used to calculate PET. Reflectance of total solar energy
is quantified by the soil albedo. A soil surface albedo of 0.2 was used for UNSAT-H
(Chudnovskii 1966; Khire et al. 1997). Meteorological data input to UNSAT-H includes daily
and hourly precipitation, daily maximum and minimum air temperature, average daily dew point,
daily total solar radiation, average daily wind speed, and average daily cloud cover. The data
recorded on site were used as input. Daily dew point temperatures were obtained from NOAA
weather station. Because net solar radiation was measured at the site, average cloud cover was
assumed to be zero. Hourly precipitation was applied as it was recorded by the tipping bucket
located at the site. However, the precipitation recorded at the site did not include snow. Hence,
snow received during winter recorded by the NOAA weather station located about 5 km from the
site was used as input. UNSAT-H does not simulate the snow melt. Hence, snow melting
approach presented in Kustas et al. (1994) was used to melt the precipitation in the form of snow
before supplied to the model. Estimated snow melt was applied as rainfall between 10:00 AM
and 6:00 PM when the solar radiation is typically higher and snow melt is observed. This
approach is consistent with the method used by Khire et al. (1997).
Hourly or daily meteorological data can be supplied to RZWQM. Hourly data of
precipitation, average air temperature, solar radiation, average wind speed, and average relative
humidity were supplied to RZWQM. RZWQM has snow module called PRMS (Precipitation
Runoff Modeling System) (Leavesley et al. 1983) to estimate the snow depth and snow melt
from the meteorological data. RZWQM requires albedo for wet and dry soil separately. The soil
albedo was set as 0.2 and 0.3 for wet and dry soil, respectively, for RZWQM (Ma et al. 2011).
76

4.3.4 Initial and Boundary Condition
Initial conditions were defined using soil suction values assigned to each node in UNSAT-H and
initial average water contents of each layer were specified in RZWQM. The simulations were
carried out individual years for the first three years of service. Field measured water contents and
soil suctions on October 31st 2009, October 31st 2010, and October 31st 2011 were used as the
initial conditions for the simulation of Year 1, Year 2, and Year 3, respectively. The initial soil
water storages for Year 1, Year 2, and Year 3 were 42.5 cm, 44 cm, and 58 cm, respectively.
The boundary conditions assigned in both the models were similar. The top boundary of the
model was assigned as a flux boundary which considers infiltration and evapotranspiration as
specified fluxes. The lower boundary was assigned as a unit gradient to simulate the vented
lower boundary of the lysimeter (Khire et. al. 1997; Benson et al. 2010).

4.3.5 Macropore Flow Data
Macropore flow module of RZWQM requires total macroporosity, average radius of cylindrical
macropore, fraction of dead end macropore, and sorptivity factor for lateral infiltration. One
macropore was identified in the field test section within the lysimeter area (Section 3.7) and the
flow rate was measured using irrigation tests (Section 3.5). Hence, the model was calibrated
using one cylindrical macropore for the lysimeter area. The average radius of the macropore was
calibrated to match the measured flow rate during the second irritation test. The fraction of dead
end macropores was set to zero because non-active macropores were not identified within the

77

lysimeter area. Also, sorptivity factor for lateral infiltration was set to 1. The sorptivity factor of
one allows lateral infiltration rate based on Green-Ampt equation without any reduction.

4.4 Evaluation of RZWQM for Water Balance Modeling of Earthen Landfill Cover
The simulations that considered only micropores were carried out using UNSAT-H and
RZWQM to quantify the micropore flow through the test section. Input parameters listed in
Table 4.1 and Table 4.2 were used for the simulations. Plants were not considered for these
simulations because the influence of plants was negligible on deep percolation, and all plant
related parameters were not measured except the root-length density function. Simulations were
carried out on yearly basis which allowed providing initial conditions at the beginning of each
year. From November 1st 2009 to October 31st 2010 was considered as Year 1. Similarly,
November 1st 2010 to October 31st 2011, November 1st 2011 to October 31st 2012, and
November 1st 2012 to October 31st 2013 were considered as Year 2, Year 3, and Year 4
respectively. However, numerical simulations were carried out only for the first three years
because the test section was dismantled in Year 4. Irrigation tests carried out during Year 1 and
Year 4 were simulated separately to obtain calibrated input parameters.

4.4.1 Simulations of Irrigation Tests
The first irrigation that was carried out during the first year after the test section was built. The
water balance of this test was simulated for two weeks starting four days before the irrigation
event. Measured and simulated cumulative percolation and cumulative precipitation and
78

irrigation are presented in Figure 4.7(a). The simulated percolation is relatively accurate. The
measured and simulated SWSs of the test section during irrigation test-1 are presented in Figure
4.7(b). The simulated SWSs are relatively similar to each other and within about 1cm compared
to the measured SWS.
Figure 4.8 presents simulated cumulative surface runoffs. Irrigation applied was about 21
cm and about 19 cm and 18 cm of the irrigation was resulted as surface runoff according to
RZWQM and UNSAT-H, respectively.
Hence, both the models simulated the water balance of the earthen cover relatively
accurately during the first irrigation test. The water balance predictions of both models are
comparable to each other.

79

300

(a)

0.8

240
Percolation
Micropore flow
simulations

Irrigation
0.6

180

UNSAT-H
0.4

120
RZWQM

Precipitation
0.2

60
Measured Percolaiton

0.0

Cumulative Precipitation and Irrigation (mm)

Cumulative Percolation (mm)

1.0

0
0

2

4

6

8

10

12

14

Day
30

(b)

Soil Water Storage, SWS (cm)

Irrigation

Simulated SWS
UNSAT-H

62

22.5

RZWQM

60

15

Measured SWS

58

7.5

Preciitaion

56

Cumulative Precipitaion and Irrigation (cm)

64

0
0

2

4

6

8

10

12

14

Day

Figure 4.7: Measured and simulated cumulative percolation (a); and soil water storage (b) during
the first irrigation test (2010).
80

30

(b)

Precipitation + Irrigation

RZWQM
UNSAT-H

Cumulative Runoff (cm)

20

24

Irrigation
15

18

10

12

5

6

0

0
0

2

4

6

8

10

12

Cumulative Precipitation (mm)

25

14

Day

Figure 4.8: Simulated cumulative surface runoff during the first irrigation test (2010).

Similar to the first irrigation test, water balance of the test section during the second
irrigation test was simulated for two weeks starting four days before the irrigation event.
Measured and simulated cumulative percolations, and cumulative precipitation and irrigation are
presented in Figure 4.9(a). The measured cumulative percolation was about 11 mm whereas the
simulated cumulative percolation was about 1 mm. The percolation measured during the second
irrigation test was influenced by macropore flow. UNSAT-H or RZWQM without macropore
flow module cannot simulate the percolation dominated by macropore flow in the fourth year of
service of the earthen cover. RZWQM simulations of water balance with macropore module are
presented in Section 4.5. Measured and simulated SWSs of the test section during the second
irrigation test are presented in Figure 4.9(b).

81

120
(a)
96

7.2

72
Measured Percolation

4.8

48
Percolation- Micropore flow
simulations
UNSAT-H RZWQM

2.4

24

0.0

0
0

2

4

6

Day

8

10

12

14

64

15

Soil Water Storage, SWS (cm)

(b)
Precipitation + Irrigation

62

10
Irrigation

UNSAT-H
RZWQM

60

5

Measured SWS
58

Cumulative Precipitation and Irrigation (cm)

Cumulative Percolation (mm)

Irrigation + Precipitation

9.6

Cumulative Precipitation and Irrigation (mm)

12.0

0
0

2

4

6

8
Day

10

12

14

Figure 4.9: Measured and simulated cumulative percolation (a); and soil water storage (b) during
the second irrigation test (2013).

82

Compare to UNSAT-H, RZWQM showed more drying, in other words more evaporation.
However, the SWSs simulated by both the models follow the trend observed in measured SWS.
Figure 4.10 presents simulated cumulative surface runoffs. UNSAT-H and RZWQM simulated
relatively similar surface runoff. The simulated cumulative surface runoff was about 5 cm while

6

15

Cumulative Surface Runoff (cm)

RZWQM
UNSAT-H
4

10
Irrigation + Precipitation

2

5

0

0
0

2

4

6

8
Day

10

12

Cumulative Precipitation and Irrigation (mm)

cumulative precipitation and irrigation was 11 cm.

14

Figure 4.10: Simulated surface runoff during the second irrigation test (2013).

4.4.2 Simulations of Year 1
Figure 4.11(a) presents the measured and simulated cumulative percolation, and recorded
cumulative precipitation and irrigation for the first year. During Year 1, cumulative water
applied to the model was about 993 mm. The simulated cumulative percolation exceeded the
83

measured percolation only by a few mm. Flow was predominantly through micropores during
the first year. Percolation simulated by UNSAT-H and RZWQM are comparable. Hence, both
models predicted percolation reasonable well when the flow was predominantly through
micropores. Figure 4.11(b) presents the measured and simulated SWSs of the test section. The
simulated SWSs followed similar trend as observed in the measured SWS except during Spring
2010. The simulated SWSs showed more drying than the measured SWS during March to May
2010. The soil was frozen to 15 cm depth during Winter 2010. The temperature of the soil
surface and the soil temperature at depth of 15 cm are presented in Figure 4.11(b). The models
do not simulate freeze-thaw process and do not account for the effect of frozen ground surface.
During Winter 2010, the model predictions were slightly different from each other because the
snow melt were simulated differently. The precipitation during Winter 2010 was applied during
early March 2010 for UNSAT-H model because the model does not have a snow melt algorithm.
RZWQM has an in-built snow melt module. After May 2010, the simulated SWSs were
identical.
Comparison of simulated cumulative surface runoff, and cumulative evaporation (E) and
cumulative potential evaporation (PE) during Year 1 are presented in Figure 4.12. Both models
produced runoff during the irrigation event only. Compared to UNSAT-H, RZWQM produced 2
cm of more runoff and 2 cm less infiltration during the irrigation event. Also, UNSAT-H
simulated 2 cm of cumulative E more compared to RZWQM.

84

10

1000
Precipitation + Irrigation

8

800
Irrigation
208 mm
Simulated Percolation
UNSAT-H

4

400

RZWQM

2

200

65

Nov/1/10

Sep/1/10

Jul/1/10

Mar/1/10

May/1/10

Measured Percolation
Jan/1/10

0

0

10

Soil Temperature
at 6"

60

0

Measured
SWS

o

Soil Surface
Temperature

55

-10

50

-20

45

-30
UNSAT -H

40

-40

Nov/1/10

Sep/1/10

Jul/1/10

May/1/10

Mar/1/10

Jan/1/10

35

RZWQM

Simulated SWS

Nov/1/09

Soil Water Storage,SWS (cm)

600

Soil Temperature ( C)

6

Cumulative Precipitation + Irrigation (mm)

1200

Nov/1/09

Cumulative Percolation (mm)

12

-50

Figure 4.11: Measured and simulated cumulative percolation (a); and soil water storage (b)
during Year 1.
85

20
(a)

Surface Runoff (cm)

RZWQM
15
UNSAT- H

10

5

Nov/1/10

Sep/1/10

Jul/1/10

May/1/10

Mar/1/10

(b) Thin line - UNSAT-H
Thick line - RZWQM

PE

100

80

Precipitation + Irrigation
E

60

40

Nov/1/10

Sep/1/10

Jul/1/10

May/1/10

Mar/1/10

0

Jan/1/10

20

Nov/1/09

Cumulative Precipitation or PE or E (cm)

120

Jan/1/10

0

Nov/1/09

Irrigation event

Figure 4.12: Simulated cumulative surface runoff (a); and cumulative evaporation (E) and
cumulative potential evaporation (PE) (b) during Year 1.
86

4.4.3 Simulations of Year 2
Figure 4.13(a) presents the measured and simulated cumulative percolation for Year 2
together with measured cumulative precipitation. During Year 2, the total water applied to the
model was about 1,038 mm and the measured percolation was 75 mm. Year 2 was the wettest
year in 50 year history (Figure 3.10a) which partly contributed to greater percolation. UNSAT-H
and RZWQM with only micropore flow simulated cumulative percolation of about 40 mm which
was about half of the measured cumulative percolation. In addition, the slope of the field
cumulative percolation was steeper than that simulated by UNSAT-H or RZWQM. Percolation
measured during Year 2 was influenced by macropore flow. Hence, those models that are based
on the Richard’s equation alone underestimated percolation. The macropore flow simulations
using RZWQM are presented in Section 4.5.
Figure 4.13(b) presents the measured and simulated SWSs of the test section. The simulated
SWSs followed similar trend as the measured SWS. However, the simulated SWSs showed less
storage than the measured SWS during Winter 2011. The soil temperature at 6 cm depth
indicated frozen ground during the Winter of 2011. The soil temperature at soil surface and the
soil temperature at depth of 15 cm are presented in Figure 4.13(b). Compared to Year 1,
simulations of Year 2 yield better match with the measured SWS. Effect of frozen ground was
less during Year 2 compared to Year 1. The frozen soil depth and duration were shorter during
Year 2 compared to Year 1. During Winter 2011, UNSAT-H and RZWQM predictions were
slightly different from each other because the snow melt were simulated differently. The snow
melt approach presented by Kustas et al. (1994) was used to melt the precipitation that was snow
before input to UNSAT-H whereas RZWQM used an in-built snow melt module PRSM
87

(Leavesley et al. 1983) to melt the snow according to the air temperature. After March 2011, the
models predictions were relatively accurate. From March 2011 to September 2011, the model
predictions of SWS were identical to each other. However, after September 2011, RZWQM
predicted slightly less SWS and more E compared to UNSAT-H. The simulated cumulative E
and PE are presented in Figure 4.14(a). Figure 4.14(b) presents the simulated cumulative surface
runoff during Year 2. The simulated surface runoff events between the models showed few
differences in Spring 2011 which could be because of the difference in snow melting approach
and RZWQM uses Green-Ampt approach to model infiltration. However, cumulative simulated
surface runoff at the end of Year 2 was similar to each other.

88

80

1200

60

900

UNSAT-H
40

600
RZWQM

70

Nov/1/11

Sep/1/11

Jul/1/11

Mar/1/11

May/1/11

300

Jan/1/11

20

0

10

Soil Temperature
at 6" depth

(b)
0

Soil Surface
Tempeature

Measured

60

-10

55

-20

UNSAT-H
50

-30

RZWQM

Nov/1/11

Sep/1/11

-60

Jul/1/11

35

May/1/11

-50

Mar/1/11

40

Jan/1/11

-40

Nov/1/10

45

o

65

Soil Temperature ( C)

0

Soil Water Storage, SWS (cm)

Precipitation

Cumulative Precipitation (mm)

Measured Percolation

Nov/1/10

Cumulative Percolation (mm)

(a)

Figure 4.13: Measured and simulated cumulative percolation (a); and soil water storage (b)
during Year 2.
89

120
Thick line - RZWQM
100

80

PE

60

E

Precipitation

40

Sep/1/11

Nov/1/11

Nov/1/11

Jul/1/11

May/1/11

Sep/1/11

8

Mar/1/11

0

Jan/1/11

20

Nov/1/10

Cumulative Precipitation or E or PE (cm)

(a) Thin line - UNSAT-H

(b)

6
5
4
UNSAT-H
3
2
1
May/1/11

Mar/1/11

Jan/1/11

Nov/1/10

0

RZWQM

Jul/1/11

Cumulative Surface Runoff (cm)

7

Figure 4.14: Simulated cumulative evaporation (E) and cumulative potential evaporation (ET)
(a); and cumulative surface runoff (b) during Year 2.
90

4.4.4 Simulations of Year 3
Figure 4.15(a) presents measured and simulated cumulative percolation for Year 3 together
with recorded cumulative precipitation. During Year 3, total water applied to models was about
784 mm and measured percolation was 137 mm. About 60% of total precipitation occurred
during Winter and Spring 2011/2012. Winter and spring precipitation can increase the
percolation through earthen cover because ET is relatively small during winter and spring
(Albright et al. 2010). Moreover, the test section was field saturated during Winter and Spring
2011/2012. Hence, initial saturated conditions, and higher precipitation during winter and spring
can cause higher percolation during Year 3. However, the slope of the field cumulative
percolation was steeper than that simulated by UNSAT-H or RZWQM indicating contribution
from macropore flow. The micropore only flow simulations using UNSAT-H and RZWQM
predicted cumulative percolation of about 12.5 cm and 9.5 cm, respectively, at the end of Year 3
which was about 90% and 70% of the measured cumulative percolation, respectively. However,
the simulated cumulative percolation using UNSAT-H at Line A and Line B shown in Figure
4.15(a) (late January 2012 and late March 2012) were about 37% and 60% of the measured
cumulative percolation, respectively. The simulated cumulative percolation using RZWQM was
even lower compared to UNSAT-H at the Line A and the Line B shown in Figure 4.15(a). The
time lag in simulated cumulative percolation compared to the measured cumulative percolation
indicates the influence of macropore flow during Year 3.
Figure 4.15(b) presents the measured and simulated SWSs of the test section. Also, the
temperature of the soil surface is presented in Figure 4.15(b). The simulated SWSs were slightly
more during Winter and Spring 2011/2012 and slightly less during Summer 2012 compared to
91

the measured SWS. Unlike during Year 1 and Year 2, ground did not freeze during Year 3 and
the simulated SWSs showed better match with the measured SWS. Also, the simulated SWSs
using UNSAT-H and RZWQM were comparable.
Comparisons of simulated cumulative surface runoffs, and simulated cumulative
evaporations and cumulative PEs during Year 3 are presented in Figure 4.16. Simulated
cumulative surface runoff using UNSAT-H and RZWQM were relatively similar. The UNSATH simulated about 5 cm less cumulative evaporation compare to RZWQM.

4.5 Macropore Flow Simulations using RZWQM
RZWQM was used to simulate macropore flow that was observed in the field test section. The
second irrigation test carried out during the fourth year of service indicated macropore flow
through the test section. Percolation measured during the second irrigation test was used to
calibrate the macropore flow component of RZWQM. The calibrated model parameters were
then used to simulate percolation during Year 2 and Year 3.

92

15

90
13.7cm

Measured

Cumulative Percolation (cm)

10.5cm

UNSAT-H
9.5cm

10
A

60

RZWQM
6.4cm

5.6cm

Precipitation

5

30

4.7cm

Cumulative Precipitation, P (cm)

12.5cm

B

2.1cm

Nov/1/12

Sep/1/12

Jul/1/12

May/1/12

Mar/1/12

Jan/1/12

1.5cm

Nov/1/11

0

0

10
Soil Surface
Temperature

70

0
o

UNSAT-H

Soil Temperature ( C)

Soil Water Storage, SWS (cm)

80

60

-10

RZWQM
50

-20

40

-30

Measured

Nov/1/12

Sep/1/12

Jul/1/12

May/1/12

Mar/1/12

Jan/1/12

-40

Nov/1/11

30

Figure 4.15: Measured and simulated cumulative percolation (a); and soil water storage (b)
during Year 3.
93

10

Cumulative Surface Runoff (cm)

UNSAT-H
RZWQM

8

6

4

Sep/1/12

Nov/1/12

Sep/1/12

Nov/1/12

Jul/1/12

May/1/12

Mar/1/12

Thin line - UNSAT-H
Thick line - RZWQM
PE

100

80

60

Precipitation
E

40

Jul/1/12

May/1/12

Mar/1/12

0

Jan/1/12

20

Nov/1/11

Cumulative Precipitation or E or PE (cm)

120

Jan/1/12

0

Nov/1/11

2

Figure 4.16: Simulated cumulative surface runoff (a); and cumulative evaporation (E) and
cumulative potential evaporation (ET) (b) during Year 3.
94

4.5.1 Calibration of RZWQM for Macropore Flow
Cylindrical macropores were conceptualized in RZWQM and calibrated for their input
parameters using field measured percolation during the second irrigation test carried out in the
fourth year. The number of macropores per unit area for a given diameter of macropore can be
calculated using Eq. 4.5. The spacing of macropores can be calculated using the number of
macropore per unit area. Several combinations of diameter and spacing of macropores were able
to simulate percolation during the irrigation test. Three sets of macropores used to simulate the
percolation are presented in Figure 4.17. Figure 4.17(a) presents the simulated cumulative
percolation with and without macropores and measured cumulative percolation. Diameter of
cylindrical macropores 1 mm, 3 mm, and 5 mm with spacing of 0.86 m, 8.6 m, and 24 m,
respectively were able to simulate the measured percolation relatively accurately. However, only
one active macropore was identified within the lysimeter through CH4 injection tests. Hence, the
spacing of macropores was assigned equal to the linear dimension of the lysimeter of 8.5 m. 3
mm diameter macropore with 8.6 m spacing was able to simulate the percolation measured
during the irrigation test.
Lateral infiltration from macropore to the soil matrix is fully allowed by setting the
sorptivity factor for lateral infiltration to 1.0. The sorptivity factor 1.0 allows the lateral
infiltration to occur according to the hydraulic gradient without any reduction. Also, a fraction of
the deadend macropores was set to zero because non-active macropores were not identified
within the lysimeter area during CH4 injection tests. The measured and simulated SWSs with and
without macropores are presented in Figure 4.17(b).

95

150

RZWQM: Cylindrical Macropore, D-Diameter, S-Spacing

D=1 mm & S=0.86 m
D=3 mm & S=8.6 m
D=5 mm & S=23.4m

12.0

120

Precipitation

9.0

90

Measured Percolation
6.0

60

Percolation- Micropore only
simulation

3.0

30

0.0

Cumulative Precipitation + Irrigation (mm)

Cumulative Percolation (mm)

15.0

0
0

2

4

6

8

10

12

14

Day

15
Precipitation + Irrigation

62

10

Micropore only
Macropore
D = 3 mm & 5 mm
Macropore
D = 1 mm

60

5

58

Cumulative Precipitation and Irrigation (cm)

Soil Water Storage, SWS (cm)

64

0
0

2

4

6

8
Day

10

12

14

Figure 4.17: Measured and simulated percolation (a); and soil water storage (b) during the
second irrigation test with micropores and macropores.

96

The simulated SWS using macropores having 3 mm and 5 mm diameter were similar.
However, when the diameter of macropore was 1 mm it increases the SWS by few mm. The
infiltration that enters macropores can be absorbed by the soil matrix depending on the degree of
saturation of the soil and the balance releases profile as percolation. When the number of
macropores is less and the diameter of the macropores large, percolation increases and relatively
small amount of water is absorbed into the soil matrix. When the diameter and spacing of
macropores decrease, the number of macropores increases. Hence, the total surface area of
macropores walls increases and the amount of lateral infiltration increases which leads to
increase in SWS. Figure 4.18 presents the simulated cumulative surface runoff with and without
macropores. The simulated surface runoff with macropore is less than the surface runoff
simulated without macropores. The cumulative surface runoff was reduced by about 1 cm
compared to micropore only simulation which was re-routed into macropores. For simulation
with macropores of diameter reduced to 1 mm, cumulative runoff was reduced further by few
mm because greater macropores wall area results in greater absorption of infiltration.

97

6

Micropore only

Cumulative Runoff (cm)

5

Macropore
D = 3 mm & 5 mm

4
Macropore
D = 1 mm

3

2

1

0
0

2

4

6

8
Day

10

12

14

Figure 4.18: Simulated surface runoff during the second irrigation test with and without
macropores.
4.5.2 Year 2 Simulation with Macropore Flow
The calibrated macropore parameters, 3 mm diameter cylindrical macropore at 8.6 m
spacing, were used to simulate percolation during Year 2. Figure 4.19(a) presents recorded
cumulative precipitation, and measured and simulated cumulative percolation for Year 2.
Combined micropore and macropore simulation using the calibrated macropore parameters
simulated cumulative percolation much better than micropore only simulation for Year 2.
However, macropore simulation predicted about 1cm less cumulative percolation (about 13%
less) compared to measured cumulative percolations. Year 2 was the wettest year in 50 years
history of the site. The estimated keff based on the percolation measurement showed higher keff
during Year 2 than during the second irrigation test (refer to Figure 3.14). Hence, macropore
parameter needed to be reassigned. Larger macropore diameter was input to match the measured
98

cumulative percolation during Year 2. Macropore diameter of 3.6 mm with spacing of 8.6 m
provided reasonably good prediction of cumulative percolation. The size of the macropore can
change in the field due to swelling or shrinkage of clay, plant root penetration and death, and
animal activity. Hence, calibration of macropore parameters needs to be done routinely to
capture the macropore flow accurately.
The measured and simulated SWSs with and without macropores are presented in Figure
4.19(b). The simulated SWS with macropores was comparable to the simulated SWS without
macropores.

99

80

Measured Percolation

D = 3 mm

40

70

Nov/1/11

Jul/1/11

May/1/11

Mar/1/11

Jan/1/11

Sep/1/11

Micropore only
Simulation

20

0

10

Soil Temperature
at 6" depth

0
Measured

60

-10

Macropore

55

-20

Micropore only

35

-60

Nov/1/11

-50

Sep/1/11

40

Jul/1/11

-40

May/1/11

45

Mar/1/11

-30

Jan/1/11

50

o

Soil Surface
Tempeature

Soil Temperature ( C)

65

Nov/1/10

Soil Water Storage, SWS (cm)

D = 3.6 mm

60

Nov/1/10

Cumulative Percolation (mm)

Macropore Spacing = 8.6 m

Figure 4.19: Measured and simulated percolation (a); and soil water storage (b) during Year 2
with and without macropores.
100

Simulated surface runoff during Year 2 with and without macropores is presented in Figure
4.20. As expected, macropores reduced simulated cumulative surface runoff compared to
micropore only simulation. The reduction in surface runoff due to macropores was very close to
the increase in simulated percolation. This shows that the model simulates all infiltration that
enters the macropores is either absorbed by the matrix or is routed as percolations.

Micropore only

7
6
5

Macropore D = 3 mm

4
Macropore D = 3.6 mm
3
2

Nov/1/11

Sep/1/11

Jul/1/11

May/1/11

Mar/1/11

0

Jan/1/11

1
Nov/1/10

Cummulative Surface Runoff (cm)

8

Figure 4.20: Simulated cumulative surface runoff during Year 2 with and without macropores.

101

4.5.3 Year 3 Simulation with Macropore Flow
Unlike the simulation of Year 2 with macropores where using the irrigation test data calibrated
macropore parameters underestimated cumulative percolation, simulation of year 3 with
calibrated macropore parameters provided excellent match with measured percolation. Figure
4.21(a) presents the measured and simulated cumulative percolation of Year 3. The estimated
maximum keff during Year 3 and Year 4 were comparable. During the winter and spring of the
Year 3 and Year 4, the test section was close to saturated compared to Year 2. When clay is
saturated, swelling can reduce the macropore cross section area. Hence, the calibrated size
simulated Year 3 better but not Year 2. Also, because the macropore size can be larger during
Year 2, Year 2 required larger macropore diameter than the calibrated value. These results show
that macropore geometries do change during service.
The measured and simulated SWSs with and without macropores are presented in Figure
4.21(b). The simulated SWSs with and without macropores were comparable. Figure 4.22
presents the simulated cumulative surface runoff during Year 3 with and without macropores.
Macropores reduced the cumulative surface runoff which is consistent with simulation for Year
2.

102

(a)
Macropore D = 3 mm
Spacing 8.6 m

Measured

10
Micropore only

70

(b)

Nov/1/12

Sep/1/12

Jul/1/12

Mar/1/12

Jan/1/12

0

May/1/12

5

Nov/1/11

Cumulative Percolation (cm)

15

Micropore only

Soil Water Storage, SWS (cm)

65
60

Macropore D = 3 mm
55
50
45

Measured
40
35

Nov/1/12

Sep/1/12

Jul/1/12

May/1/12

Mar/1/12

Jan/1/12

Nov/1/11

30

Figure 4.21: Measured and simulated cumulative percolation (a); and soil water storage (b)
during Year 3 with and without macropores.
103

10

8

6

Macropore D=3mm
4

Nov/1/12

Sep/1/12

Jul/1/12

May/1/12

Mar/1/12

0

Jan/1/12

2

Nov/1/11

Cumulative Surface Runoff (cm)

Micropore only

Figure 4.22: Simulated cumulative surface runoff during Year 3 with and without macropores.

104

CHAPTER 5: LABORATORY TESTING OF MACROPOROUS COMPACTED
CLAY SAMPLES

Testing macropore flow through cracks in fine-grained soil in laboratory is a challenging task
due to the self-healing nature of cracks. Increase effective stress, erosion of particles form
macropore wall clogging the cracks, and swelling of clay are attributed as major causes of this
self-healing nature of cracks (Eigenbrod 2003; Wang et al. 2013; Kim and Daniel 1992; Othman
and Benson 1993). Self-healing mitigates the macropore flow by partially closing the cracks.
However, evaluation of macropore flow in laboratory is needed to improve our understanding of
macropore flow phenomena in clay. Ahuja et al. (1995) used quick-setting cement to prevent the
closure of macropores during laboratory test in their experimental setup. A new technique was
developed in this study to create macropores in the lab. Crystal sugar supported by a thin highly
permeable insect mesh was inserted within the clay specimen during compaction. The sugar
dissolved during permeation and formed macropores where sugar once occupied the space.
Control tests were carried out to test the effect of sugar in the permeant on permeability of
compacted clay specimen. Double ring rigid wall permeameter was used to measure the
hydraulic conductivity of the compacted specimens. High resolution X-Ray CT scans were
performed on compacted clay specimens to develop 3D morphology of the macropores.
Macroporosity and morphology of macropores were analyzed using X-ray CT images. Also, the
X-ray CT images were used as direct input to LBM which is described in Chapter 6.

105

5.1 Self-healing of Cracks in Compacted Clay Soil
The self-healing behavior of clay used to construct the field test section (Section 3.1) was
evaluated. Particle size distribution of the specimen is presented in Figure 5.1(a). The uniformity
coefficient and coefficient of curvature of the clay are 70 and 0.7, respectively. Percentage of
fines in the clay soil is about 78%. The SP compaction test results are shown in Figure 5.1(b).
The optimum water content and maximum dry unit weight of the clay for SP effort was 13.4%
and 18.5 kN/m3, respectively. The LL and PI of the clay is 29 and 13, respectively. The SL of the
clay as predicted using the LL and the PI based on a procedure presented in Holtz and Kovacs
(2010) is 13.25. The clay is classified as CL according to USCS.
The first specimen, 30 cm diameter and 30 cm tall, to evaluate the self-healing behavior was
compacted at 4% dry of optimum water content in 6 cm thick lifts. Achieved dry unit weight was
about 98% of the maximum dry unit weight. The compaction criteria were similar to the field
compaction of the test section. The specimen was initially saturated and then dried for about
three months. Cracks were formed but the cracks were mostly shallow. Figure 5.2(a) presents the
picture of the cracked clay specimen. The specimen was wetted by gradual spraying. The cracks
partially self-healed and closed as shown in Figure 5.2(b). Similar observation was reported by
Greve et al. (2010).

106

(a)
Percentage Finer by Weight (%)

100.0

Sieve analysis data

80.0
Hydrometer data

60.0

40.0

20.0

0.0
100

10

1

0.1

0.01

0.001

0.0001

Particle Size (mm)
20.0
(b)

Standard Proctor test

3

Dry Density (kN/m )

19.0

Zero air void line
18.0

17.0

16.0
10

12

14
16
18
Molding Water Content (%)

20

Figure 5.1: Particle size distribution (a); and compaction curve (b) of the clay used for
macroporous specimens.
107

(a)

(b)

(c)

(d)

Figure 5.2: Cracked (a); and self-healed (b) sample-1; and cracked (c); and self-healed (d)
sample-2.

108

The second specimen was compacted in four lifts and two water content and water potential
sensors were placed in between each lifts. The specimen was compacted in a PVC specimen
holder. Compaction criteria used for the first specimen was followed. The specimen was initially
saturated with about 2 cm head of water ponded on the specimen and then the water level was
increased to about 90 cm to increase the gradient. Measured water contents indicated that the
specimen reached saturation in about a month. Measured water content is presented in Figure
5.3. The saturated specimen was dried for about a month. Cracks formed in the specimen as
shown in Figure 5.2(c). The specimen wetted by ponding of water. Initially ponding height was
about 2 cm and then the water level was increased to 90 cm. The measured water contents
increased suddenly after the water was ponded indicating the cracks were developed to at least
21 cm depth. The water contents dropped to a steady-state level in about a day indicating the
self-healing of cracks.

5.2 Macroporous Specimen Preparation
The specimens were compacted in 10 cm diameter Proctor compaction mold. Standard procedure
for SP compaction test was followed (ASTM 2007). However, food grade crystal sugar packed
in an insect mesh was placed within specimen at each lift as shown in Figure 5.4. The specimen
was compacted. The continuity of sugar column was maintained between lifts. Within a lift the
shape of the sugar packed space was straight or tortuous. Three macroporous specimens were
prepared namely S1, S2, and S3. S1 with two vertical columns of sugar, S2 with one vertical
column of sugar, and S3 with one tortuous column of sugar were prepared. After compaction the
specimens were retrieved using a hydraulic extruder.
109

H = 2 cm

3

Volumetric Water Content, VWC (cm /cm )

0.45
H = 2 cm

H = 90cm

H = 90cm

Drying

0.4

3

7 cm
depth

0.35

14 cm
depth

0.3
0.25

21 cm
depth

0.2
0.15
0.1
0.05

H - Ponded water level

Nov/18 Dec/2 Dec/16 Dec/30 Jan/13 Jan/27 Feb/10 Feb/24 Mar/10

3

3

Volumetric Water Content, VWC (cm /cm )

0.5

Depth:
14 cm

0.4

0.3

7 cm
0.2

0.1

21 cm
0
Jan/24

Jan/24

Jan/24

Jan/25

Jan/25

Jan/25

Jan/26

Jan/26

Figure 5.3: Measured water content of Sample 2 during drying and wetting.
110

(a)

(b)
Sugar
Insect
Nest

Proctor Mold
Figure 5.4: Macroporous specimen preparation: at the end of the first lift (a); and at the end of
compaction (b).

The specimens were tested in flexible wall permeameter (ASTM 2010a). Two control tests
were carried out on intact specimens using only DI water and DI water with dissolved sugar as
the permeant to test the effect of sugar on hydraulic conductivity of the clay specimens. The
concentration of sugar in the permeant (20 g/L) was decided based on the amount of sugar added
to and the pore volume of the macroporous specimens. Flexible wall permeameter results are
presented in Figures 5.5 and 5.6. The ksat of specimens without macropores (control) were 1×10-6
cm/s. and 5×10-7 cm/s. Ratio of outflow to inflow reached about one when one pore volume of
flow had passed. The cumulative volume change [Figure 5.6(b)] of the specimen shows
shrinkage at pore volume equal to 1.6 due to increase in confining pressure. The hydraulic
conductivity decreased when sugar solution was used as the permeant and the tests needed to be
stopped before they reached steady state duo to time limitations. Also, detailed analysis on how
111

sugar water influences the permeability of clay soil is not considered in this steady. However, the
tests confirmed that the sugar solution decreases the hydraulic conductivity. Hence, increase in
hydraulic conductivity of macroporous specimens would be due to the macropores and not due to
the sugar in the permeant.
Hydraulic conductivity of macroporous specimens were measured using flexible wall
permeameter. The influent, effluent, and confining pressures were maintained at 13.8 kPa, 3.4
kPa, and 20.7 kPa, respectively, for specimens S1 and S2. The confining pressure of specimen
S2 was increased to 27.6 kPa after the pore volume reached about 2.6 to study the effect of
confining pressure on macropores. For specimen S3, the influent, effluent, and confining
pressures were maintained at 3.4 kPa, 0.0 kPa, and 6.9 kPa, respectively. The sugar used to
create the macropores was flushed out once the influent was permeated. Figure 5.7(a) presents
the measured ksat of macroporous specimens using flexible wall permeameter. Measured ksat of
all three macroporous specimens was about 1×10-3 cm/s which is the permeability of the flexible
wall permeameter. Hence, the ksat of macroporous specimens can be ≥ 1×10-3 cm/s. The ksat of
macroporous specimens were measured using double ring rigid wall permeameter (discussed in
Section 5.3). The effect of increase in confining stress tested in S2 is presented in Figure 5.7(b).
Increase in confining pressure reduces ksat because increase in effective stress closes the crack
(Kim and Daniel 1992). The decrease in volume as the confining stress increased shown in
Figure 5.7(b) is an indication of healing of crack. Similar finding on natural macropores were
presented by Kim and Daniel (1992), Othman and Benson (1993), and Eigenbrod (2003). Hence,
the artificial macropores developed in this study behave similar to natural cracks when effective
stress increases.

112

Hydraulic Conductivity (cm/sec.)
Hydraulic Conductivity (cm/sec.)

10

-4

10

-5

Permeating Fluid
Open symbol: DI water
Closed symbol: sugar solution

Control Specimen-1
10

-6

10

-7

10

-8

Control Specimen-2

0

10

-4

10

-5

10

-6

10

-7

10

-8

50

100

150
Time (Hrs.)

200

250

Permeating Fluid
Open symbol: DI water
Closed symbol: sugar solution

Control Specimen-1

Control Specimen-2

0

0.4

0.8
1.2
Number of Pore Volume

1.6

Figure 5.5: Measured hydraulic conductivity of compacted clay intact samples with time (a); and
with number of pore volume (b).
113

2.5

(a) Permeating Fluid
Open symbol: DI water
Closed symbol: sugar solution

Ratio of Outflow to Inflow

2.0

Control Specimen-1
1.5

1.0

Control Specimen-2
0.50

0.0

0

0.4

Cumulative Volume Change of the Samples
(Exapansion +)

10

0.8
1.2
Number of Pore Volume

1.6

(b)

Control Specimen-2

5.0

0.0

-5.0

Control Specimen-1
-10

-15

Permeating Fluid
Open symbol: DI water
Closed symbol: sugar solution
0

0.4

0.8
1.2
Number of Pore Volume

1.6

Figure 5.6: Ratio of outflow to inflow (a); and cumulative volume change of the specimens (b)
during flexible wall permeability test.
114

10

-2

(cm/sec)

(a) S1, S2, and S3

Hydraulic Conductivity, k

sat

S3

10

-3

S2
S1

10

10

-4

0

0.5

1
1.5
Number of Pore Volume

2

2.5

-2

-2

3

(cm/sec)

Volume Change

-6

sat

Hydraulic Conductivity, k

Volume Change of the Sample (cm )
(Shrinkage)

(b) S2

10

-3

k

sat

-10

-14
10

-4

0

2

4

6

8

10

12

Number of Pore Volume

Figure 5.7: Measured hydraulic conductivity of macroporous specimen using flexible wall
permeameter.
115

5.3 Measurement of Macropore Flow
Because flexible wall permeameter with higher ksat measurement capability was not available, a
double ring rigid wall permeameter (DRRP) was used to measure ksat of the macroporous
specimens. Figure 5.8 shows the DRRP setup. Common issue with rigid wall permeameter on
measuring ksat of fine grain soil is side wall leakage (Boynton and Daniel 1985). The DRRP has
two rings on the bottom plate such that the areas enclosed by both the rings have the same area
as shown in Figure 5.8(a). The purpose of the outer ring is to collect water through the area of
specimen within the outer ring and evaluate side wall leakage (if any). Hence, water collected
from inner ring is not affected by the side wall leakage. The DRRP setup could hold standard
compaction specimens having 10.16 cm (4 inches) diameter and 11.64 cm (4.58 inches) height
(Figure 5.8(b)). Diameter of the inner ring is 7.18 cm. Geotextile was placed at the bottom within
the rings and at the top of the specimen. The specimens were pushed through the inner ring
divider as shown in Figure 5.8(c).
Constant head tests were conducted on macroporous specimens as shown in Figure 5.8(d).
The hydraulic gradient was maintained at about one. The specimens exhibited swelling during
permeation. Figure 5.9 shows the specimen before and after the test. Before the test, few mm of
spacing can be noticed between specimen perimeter and rigid wall cell. However, after the test,
the gap between the specimen perimeter and the rigid wall of the perimeter was closed or
reduced due to swelling of the clay during saturation. The permeability of the inner and outer
rings was measured separately. The macropores were within the inner ring on all three
specimens. The measured ksat for all three specimens are presented in Table 5.1. The measured
flow through the outer ring indicated some degree of side wall leakage except for specimen S1.
116

(a)

(c)

(b)

(d)

Figure 5.8: Double ring rigid wall permeameter (DRRP) setup: bottom plate (a); specimen inside
the rigid wall cell (b); bottom of specimen (c); and test setup (d).

The overall ksat of specimens were 0.06 cm/s, 0.05 cm/s, and 0.02 cm/s for S1, S2, and S3,
respectively. However, conductivity of the setup is 0.17 cm/s. Generally the conductivity of the
setup needs to be at least an order of magnitude higher than the conductivity of the specimen to
117

measure the conductivity accurately. Hence, except S1, the ksat of S2 and S3 could be greater
than listed in Table 5.1.

Figure 5.9: Macroporous specimen before test (a); and after test (b).

Table 5.1: Measured ksat of macroporous specimens
ksat (cm/s)

Specimen ID

Inner ring (macropore)

Outer ring

Overall

S1

0.12

0.0000

0.06

S2

0.10

0.0060

0.05

S3

0.03

0.0007

0.02

118

5.4 High-Resolution X-Ray CT Images
High resolution X-Ray CT scans of the compacted clay specimens were used to develop 3D
morphology of the macropore. The X-Ray CT facility at the Texas Transportation Institute,
Texas A&M University performed the X-Ray CT scans. Individual scan images are called slices.
An individual slice is shown in Figure 5.10(a). The vertical resolution (∆z) of the greyscale
images was 1.0 mm that means the vertical distance between two slices is 1.0 mm. The
horizontal resolutions (∆x and ∆y) can be calculated directly from the specimen diameter and the
pixel of the scan images. The horizontal resolution was 0.1688 mm/pixel in both directions. The
original grayscale image was converted to binary images using specific threshold where black
color was assigned for soil matrix and white color was assigned to macropores. The threshold
was selected based on the average value of variation of threshold along the macropore. Figure
5.10(b) presents the binary image developed using average threshold (50) and Figure 5.10(c)
presents variation of threshold along the cross section Line A-A’ shown in Figure 5.10(a).
However, the upper limit of threshold slightly varied depending on the orientation of the Line AA’ [Figure 5.10(a)]. Hence, sensitivity analysis was carried out to understand the effect of
threshold value on macropore parameter. The average threshold value used for images of
specimens S1, S2, and S3 were 55, 50, and 54, respectively.
XCAT image analysis software developed by Kutay et al. (2010) was used to develop 3-D
structure of the macropores within the specimens. The binary images developed were used as
input to XCAT software. The 3-D structures of macropores are presented in Figure 5.11. The
long continuous macropores that can be identified in the Figure 5.11 are the artificial
macropores. Isolated small macropores visible in Figure 5.11 are inter-clod voids.
119

(a)

(b)

A’

A

Scale:
10 mm
118 mm
(c)

Figure 5.10: X-ray CT scan image: original greyscale image (a); binary image (b); and variation
of threshold along the Line A-A’ (c).

120

(a) S1

(b) S2

Inter-clod voids

(c) S3

Macropore
formed by
crystal sugar

Figure 5.11: 3D structure of macropores: S1 (a); S2 (b); and S3 (c).

121

5.5 Macroporosity and Characterization of Macropores
The characterization of macropores consisted of measurement of the average diameter, minimum
width, and maximum width. High resolution X-ay CT scan images were used to characterize the
macropores and to measure the Macroporosity. The formed macropore was isolated from interclod voids [Figure 5.11(c)] in order to characterize the formed macropore. Areas showing
formed macropore and soil matrix on each slice were calculated by summing up the white color
pixels and black color pixels, respectively, on each slice. The micropores within the soil matrix
were not identified separately from soil matrix and hence not included in the macroporosity. The
total macroporosities of each specimen were estimated by adding the macropore areas from each
slice and diving by the total area of all slices. However, the resolution of the slices along the
vertical direction is larger than the resolution in the horizontal direction. Hence, the macropore
areas on each slice were multiplied by three in order to achieve similar resolution in all three
dimensions. Further details on construction of 3-D image structure with uniform resolution are
explained in Section 6.4. Estimated macroporosities are presented in Table 5.2. The specimen S1
had the maximum macroporosity of 6.84% followed by specimen S2 and S1 which had
macroporosity of 0.51% and 0.28%, respectively.
The equivalent radii were estimated assuming the macropore is cylindrical and having the
same cross sectional area as the estimated formed macropore area for each slice. The equivalent
radius for specimen S1 was calculated using the area of the only one continuous macropore. The
other macropores dead-ended and hence was not included in the analysis. Figure 5.12 presents
the frequency count of equivalent radii of macropores for each specimen. The minimum, mean,
mode, and maximum radius of macropores are tabulated in Table 5.2. The maximum equivalent
122

radius for specimen S1, S2, and S3 were 6.36 mm, 5.83 mm, and 5.01 mm, respectively.
Similarly, the minimum equivalent radius for specimen S1, S2, and S3 were 3.340 mm, 1.973
mm, and 1.966 mm, respectively.
Table 5.2: Characteristics of formed macropore and estimated macroporosity

Specimen Macroporosity
ID
(%)
Min.
S1
6.84
3.340

Radius (mm)
Mean Mode

Max.

4.610 4.096

6.365

S2

0.51

1.973

3.605 3.488

5.834

S3

0.28

1.966

2.604 2.574

5.010

S3

Frequency

41.0

31.0

S2

S1

21.0

11.0

1.0
1

2

3

4
5
Radius, R (mm)

6

7

Figure 5.12: Frequency distribution of equivalent radius of the formed macropore.

123

CHAPTER 6: LATTICE BOLTZMANN MODELS

Typically, the fluid flow simulations start from Navier-Stokes (NS) equations, which are
nonlinear partial differential equations for the conservation of mass and momentum. The NS
equations are discretized by finite differences, finite volumes, and finite elements, or spectral
methods and solved by common numerical methods. This approach is called ‘top-bottom’
approach (Wolf-Gladrow, 2000). The LBM follows a so-called ‘bottom-up’ approach, where it
starts from a discrete lattice model and reaches to the solution of NS equations as it evolves in
lattice space (Wolf-Gladrow, 2000). This chapter presents the LBM algorithms used for this
study. LBM algorithms developed in this study were first verified using analytical solutions, and
then validated using X-ray CT images of macroporous (cracked) clay specimens and ksat
measured in laboratory. A systematic analysis of morphology and tortuosity of macropores on
macropore flow rate was conducted. An equation was formulated to predict flow rate through
arbitrary shape and tortuous macropores using the flow rate of straight and vertical cylindrical
macropore which has equivalent sectional area. In addition, multiphase LBM was tried to model
unsaturated flow through macropores.

6.1 Saturated (single phase) LBM

In the LBM, the fluid phase is discretized according to the conservation of macroscopic
properties and the necessary symmetries required by NS equations (He and Luo, 1997). The
average movement of fluid particles are represented by distribution function, ( ), which is
124

defined for each lattice vector

(microscopic velocity in ith direction) at each site x. D2Q9 and

D3Q19 are two lattice models of DxQy series proposed by Qian et al. (1992), where x is the
space dimension and y is the number of discrete velocity vectors in the model. The microscopic
velocities for the lattice models D2Q9 and D3Q19 are illustrated in Figure 6.1.

e14=[-1,0,1]

e15=[0,1,1]
e5=[0,0,1]

e11=[1,0,1]

e18=[0,-1,1]

e6=[-1,1]

e2=[0,1]

e5=[1,1]
e8=[-1,1,0]

e3=[-1,0] e9=[0,0]

e2=[-1,0,0]

y

e7=[-1,-1]

e4=[0,-1]

e1=[1,0,0]

e4=[0,-1,0]

e8=[1,-1]
z

e7=[1,1,0]

e19=[0,0,0]

e1=[1,0]
e9=[-1,-1,0]

e3=[0,1,0]

e10=[1,-1,0]

e16=[0,1,-1]
e13=[-1,0,-1]
y

e6=[0,0,-1]
e17=[0,-1,-1]

e12=[1,0,-1]

x

x
(b)
(a)
Figure 6.1: Microscopic velocity directions of D2Q9 model (a); and D3Q19 model (b).

Similar to the implementation by Sukop and Thorne (2006), uniform unit particle mass (1
mu) is assumed, therefore the microscopic velocities and momenta are always equivalent. The
evolution of distribution function is done in two steps for each time step t: (i) one is advection
and (ii) the other one is collision. The advection of the distribution function occurs along the
microscopic velocity vectors and it can be represented as:

f i in ( x, t )  f i out ( x  ei , t  1)

Eq. 6.1

where fiin is the incoming distribution function from the neighboring nodes to the current node x
at time t; and fiout is the outgoing distribution function at the neighboring node x-ei at the previous
125

time step t-1. After the advection, a new distribution function, called equilibrium distribution
function ( f i eq ), is determined at each site to represent the effect of collision. The collision step
includes the ongoing collisions between particles by weighting the distribution function at a site
with its f i eq . The equilibrium distribution function depends on the fluid macroscopic velocity (u)
and density (ρ) at a node x, which can be computed as:

 ( x)  i f i in ( x)


u

i

f i in ei



Eq. 6.2

Eq. 6.3

f i eq for models D2Q9 and D3Q19 is:

9
3


2
f i eq  Wi 1  3ei .u   ei .u   u.u 
2
2



Eq. 6.4

where Wi is weight factor in ith direction. The weight factors differ based on the lattice models.
The weight factors of D2Q9 are: Wi =1/9 for i= 1 to 4 (face-connected particles); Wi =1/36 for i=
5 to 8 (edge-connected particles); and Wi =4/9 for i= 9 (rest particles). Similarly, the weight
factors of D3Q19 are: Wi =1/18 for i= 1 to 6 (face-connected particles); Wi =1/36 for i= 5 to 8
(edge-connected particles); and Wi =4/9 for i= 9 (rest particle). The numbering of the ith direction
is similar to the microscopic velocity numbering in Figure 6.1. More details on the derivation of
the weight factors can be found in Wolf-Gladrow (2000). Using the linear relaxation presented
by Bhatnagar et al. (1954), also known as BGK approximation, the collision step can be written
as:

126

f i ( x, t )  f i ( x, t ) 
out

in

f i in ( x, t )  f i eq ( x, t )



Eq. 6.5

where τ is the relaxation time which controls the rate at which particle distribution relaxes to
equilibrium. The relaxation time relates to the kinematic viscosity of the fluid (v) with
(

) in lattice units. The cs is the speed of sound in lattice units, which is equal to ⁄√

for D2Q9 and D3Q19 LB models.
The body force (BF) is applied to the system by adding the force term to the collision step
(Eq. 6.5) in the following form (Martys et al. 1998):

BFi  

Wi
 ei .g 
c s2

Eq. 6.6

where g is the gravitational acceleration; and  is density. In LBM used in this study, the
following form of the equation of state (EOS) defines the pressure (P) as follows:

P  c s2 

Eq. 6.7

Hence, pressure boundary is simulated by setting the boundary densities according to Eq. 6.7.

6.1.1 Real-numbered Solid Density Micropore and Macropore Flow Model
The permeable solid (i.e., the clay soil medium in this study) is implemented using the partialbounceback model proposed by Walsh et al. (2009). The partial-bounceback is implemented as a
secondary collision step in standard LBM as follows:

f i out ( x, t )  1  ns  f i out ( x, t )  ns f i in ( x, t )
127

Eq. 6.8

where ̅ indicates the opposite direction vector; and ns is the damping parameter related to solid
fraction. The value of ns ranges from 0 to 1. When ns is equal to 0, the model recovers the
standard LB model (i.e., equations Eq. 6.1 to Eq. 6.5). When ns is equal to 1, Eq. 6.8 reduces to
full-bounceback condition, i.e., the no-slip boundary condition.

6.2 Theoretical Verification of the D2Q9 and the D3Q19 Model
The D2Q9 permeable solid model was first verified using Poiseuille flow through a slit between
two parallel plates. The velocity at the plate boundaries are zero and follow a parabolic velocity
distribution between the plates. The velocity distribution between two parallel plates which is 2a
distance apart is given by,

uy 

Py
2

a

2

 x2



Eq. 6.9

where u y is the velocity in y direction; Py is the pressure gradient; µ is the dynamic viscosity of
the fluid; and

is the distance from centerline between parallel plates. Numerical simulations

were performed in a 50×50 lattice domain representing 5 mm each side. The relaxation time and
the density of the fluid were set to 0.8 and 1 in lattice units, respectively. No slip boundary
condition was applied at the wall boundaries using full-bounceback.
Pressure gradient Py of

⁄

was applied between top and bottom

boundary. Unknown distribution function components at the top and bottom boundaries were
calculated using the formulation given in Kutay et al. (2006). LBM simulation and analytical

128

solution are plotted in Figure 6.2 and the LBM simulation matches very well with the analytical
solutions.
0.50

y

u (mm/sec.)

0.40

0.30

0.20
LBM
Analytical

0.10

0.00
0

10

20

30

40

50

x (mm)

Figure 6.2: Verification of D2Q9 permeable solid model using Poiseuille flow simulation
between two parallel plates.

D3Q19 permeable solid model was verified by simulating flow through a cylindrical void
and comparing the velocity profile with the velocity profile given by the analytical solution (i.e.,
Poiseuille Law). A cylindrical cavity with a diameter (D) of 50 mm was utilized. Pixel resolution
(dx) and time step (ts) used were 1.0 mm/pixel and 0.1 sec/ts, respectively. No slip boundary
condition was applied at the boundary of the solid wall/water interface by assigning ns =1 for
solids located outside the cylinder. Pressure gradient (

) of 1.67

⁄

was

applied between the inlet and outlet boundaries. Unknown distribution function components at
129

the inlet and outlet boundaries were calculated using formulations given in Kutay et al. (2006).
The relaxation time and the density of the fluid were set to 0.8 and 1 in lattice units, respectively.
The Poiseuille’s analytical solution of velocity distribution of laminar flow through a cylinder is:
(
where

)

Eq. 6.10

is the velocity in the z-direction; and r is the distance from the centerline of the

cylinder in radial direction. Figure 6.3 shows the comparison of analytical solution and LB
simulation results where an excellent match is visible.

Figure 6.3: Verification of D3Q19 model using Poiseuille flow simulation through a cylindrical
pipe.

Dardis and McCloskey (1998) presented analytical solution for velocity distribution in y
direction (

) between parallel plates filled with solid scatters as follows:
130

( (

(

Where

(

)
)

)

Eq. 6.11

; ns = damping parameter related to solid fraction (see Eq. 6.8);

√ ⁄ ;

=kinematic viscosity; a= half of the gap between the plates; and x= distance from the
centerline. The derivation of this solution (Eq. 6.12) from the Navier-Stokes equations is
presented in Dardis and McCloskey (1998). Numerical simulations were performed in a 50×50
lattice domain representing 50 mm each side. Time step used was 0.1 sec/ts. Pressure gradient
(

⁄

) of 1.33

was applied between the inlet and outlet boundaries. Different

values of solid fraction were simulated. No-slip boundary condition was applied at the wall
boundaries using full-bounceback. Comparison of LBM simulation results and analytical
solution is shown in Figure 6.4 where an excellent match is visible.
y

-3

1.4 10

ns = 0.05

-3

1.2 10

x

LBM n = 0.05
s

Analytical n = 0.05

-3

y

u (mm/sec.)

1.0 10

s

LBM n = 0.20

y ns = 0.20

s

-4

Analytical n = 0.20

8.0 10

s

LBM n = 0.50
s

-4

6.0 10

x

Analytical n = 0.50
s

-4

4.0 10

y

ns = 0.50

-4

2.0 10

x
0.0
0

5

10
15
x (mm)

20

25

Figure 6.4: Velocity profiles between two parallel plates for different ns values
131

Since typically the hydraulic conductivity (k) is used to represent the permeability of the
soils, a relation between the ns used in the LB formulations and the k is needed. The analytical
expression that relates the intrinsic permeability (K) to ns was derived by Walsh et al. (2009) as
follows:
(

)

Eq. 6.12

where =kinematic viscosity. Using the relationship between the intrinsic permeability and the
hydraulic conductivity k,
Eq. 6.13
following equation can be derived,
(

)

Eq. 6.14

In order to verify this formulation and the accuracy of the LB algorithm, LB simulations
were conducted in a 3D domain of 50×50×50 (corresponding to a 0.5mm×0.5mm ×0.5mm
domain) that is assumed to have ns values ranging from 0.01 to 0.99. Time step used for
simulations was 1.0×10-5 sec/ts. Periodic boundary condition was applied to all the boundaries.
A gravitational body force was applied in the z-direction and resulting average velocity (

) in

the direction of gravity was used in Darcy’s formulation to compute the hydraulic conductivity
(i.e.,

, where i = hydraulic gradient = 1 because of gravity driven flow). Then this LB-

based k is compared with the k computed using the Eq. 6.14 for different ns values. As shown in
Figure 6.5, LB-based k and k computed using Eq. 6.14 match very well, verifying the algorithm
implementation as well as the formulation given by Walsh et al. 2009.
132

1

Hydraulic Conductivity, k (mm/sec.)

10

LBM
Analytical

0

10

10-1

-2

10

10-3

-4

10

0.00

0.20

0.40

0.60

0.80

1.00

ns

Figure 6.5: Hydraulic conductivity curve under gravity driven flow using analytical solution by
Walsh el al. (2009) and LB model prediction.

6.3 Laboratory Validation of the LBM
The laboratory measurement of ksat of macroporous specimens and X-ray CT scan images of
those specimens were used to validate the LBM. More detail on laboratory tests and X-ray CT
scan images can be found in Chapter 5. The binary images developed from X-ray CT scan
images were used as direct input to the LBM. The binary images were stacked to reconstruct the
3-D structure of the Specimens. However, the resolutions (scaling relation from LBM to physical
unit) of the images along the vertical and horizontal directions were not similar. Also, the
original horizontal resolution of 0.1688 mm/pixel would take longer computer time and memory
133

to simulate. Hence, the horizontal resolution of 0.3375 was achieved by resizing the image slices
using bilinear interpolation. Moreover, to achieve a uniform resolution in vertical and horizontal
directions, each slice were repeated three times while stacking to construct the 3-D structure of
the specimen as depicted in Figure 6.6. By repeating three times each slide, the original spacing
between slices was divided in to three times which results in a new spacing of about 0.3333 mm
(1/3 mm). The new spacing between slices is very close to the horizontal pixel spacing of 0.3375
mm. Similar approach to achieve a uniform resolution in all three directions of X-ray CT images
was followed by Sukop et al. (2013).

S1-035
0.3333 mm
S1-035
1.0 mm

0.3375 mm/pixel

S1-035

S1-036
0.3375 mm/pixel

Figure 6.6: Construction of 3-D structure of the specimen using image slices. Two image slices
of specimen S1 namely S1-035 and S1-035, and interpolation based on three repetitions of slice
S1-035.

134

The relaxation parameter () was set to 1.0 for all the simulations, which results in time step
of 0.019 sec/ts. The ns was set to zero which means soil matrix is impermeable. Assuming
impermeable soil matrix is reasonable since the flow through soil matrix is very small (ksat is
1×10-5 mm/sec.) compare to the flow through macropores. Pressure gradient was applied
between top and bottom boundary. The flux was calculated as sum of the velocities of all the
nodes on a slice. The difference in inflow and outflow of fluid mass was monitored. The
difference between inflow and outflow reaches zero when the system reaches steady-state. The
intrinsic permeability (K in pixel2) was calculated using Darcy’s law,

K  q

L
P

Eq. 6.15

where q is the Darcy flux (i.e., the flow rate per area in pixel/ts); and L is the length (in pixel)
between inlet and outlet of the domain. Since the pressure gradient was applied using density of
fluid at inflow and outflow, the average density (in mu.pixel3) was used for the intrinsic
permeability calculations. ∆P (in mu.lu-1.ts-2) can be calculated using the equation Eq. 6.7. The
intrinsic permeability calculated in LB units can be converted to physical unit (in mm 2) using the
2

scaling relation which is the resolution of the images. That is Kphysical

 L physical 
 where
= KLB× 
 LLB 

 L physical 

 is equal to the resolution of the image (0.3375 mm/pixel). The ksat can be calculated
 LLB 
using Eq. 6.13. Figure 6.7 presents the calculated ksat using LBM and measured ksat in laboratory.
The calculated ksat is at least one order of magnitude higher than the measured ksat. And the
difference increase with the increasing ksat. The measured ksat of specimen S1 and S2 can be
limited by the apparatus since the conductivity of the apparatus is very close to the measured k sat.
135

However, measured ksat of specimen S3 is about an order of magnitude lower than the
conductivity of the apparatus.
The occurrence of non-Darcy flow can reduce the ksat at higher hydraulic gradient (Sukop et
al. 2013). The laboratory measurement of ksat was carried out using about unit hydraulic gradient.
However, applying hydraulic gradient above 1×10-4 is not feasible in LBM on large domains.
Sukop et al. (2013) showed considerable reduction in ksat due to non-Darcy flow above hydraulic
gradient of 1×10-6. Effect of hydraulic gradient on calculated ksat of specimen S2 and S3 was
studied by varying hydraulic gradient from 1×10-9 to 1×10-4. Figure 6.8 presents the calculated
ksat for different hydraulic gradients. The calculated ksat did not show difference with varying
hydraulic gradient. The observation of calculated ksat insensitive to hydraulic gradient is
consistent for both specimens. Also, a simulation with lower resolution (0.1688 mm/pixel) was
carried out on specimen S3 to verify the effect of resolution on calculated k sat and presented in
Figure 6.8. The effect of resolution on calculated ksat was negligibly small.

136

100.0

using LBM

S1

Limited by
the apparatus

S2

Simulated k

sat

10.0

S3

1.0

Conductivity of the apparatus

0.1
0.1

1

10

Measured k

100

sat

Figure 6.7: Measured and calculated ksat of macroporous specimens.

S1, TH=56

Simulated k

sat

using LBM (mm/s.)

100.0

S2, TH=50

S3, TH=54,dx=0.1688

10.0

S3, TH=54
dx=0.3375
TH-Threshold
1.0 -9
10

-8

10

-7

10

-6

10

10

-5

10

-4

Hydraulic Gradient

Figure 6.8: Effect of hydraulic gradient and resolution on calculated ksat of macroporous
specimens using LBM.
137

Also, the effect of threshold was tested. The average number of variation of threshold along
the macropore was selected for each specimen (explained in section 5.4). However, the average
number can vary slightly. Hence, the ksat of specimen S3 was calculated using threshold of 49,
54, and 59 where 54 was the average threshold. The calculated ksat is presented in Figure 6.9.
The difference between the selected average threshold and 5 threshold upper or lower produced
only about ±15% of calculated ksat using selected average threshold.

Simulated k

sat

using LBM (mm/s.)

10.0

sample - S3
dx=0.3375
Hydraulic Gradient = 3.59E-07
1.0
48

50

52

54

56

58

60

Threshold

Figure 6.9: Effect of selected threshold on calculated ksat of macroporous specimens
using LBM.

Unlike macroporous rock (Sukop et al. 2013) or compacted aggregate (Kutay et al. 2006), the
size of macropores within compacted clay can change due to swelling of clay particles. The
138

specimens were relatively dry when the specimens were received from X-ray CT scan facility.
The swelling of the specimens during laboratory tests after the X-ray CT is presented in Figure
5.9 and discussed in section 5.3. The decrease in macropore size due to swelling of clay can
affect the comparison between measured vs. calculated ksat. Hence, the effect of swelling on
calculated ksat of macroporous specimens was studied through simulation flow through different
sizes of artificial vertical macropore. The shape of the artificial macropore was roughly vertical
cylinder. The difference between the sizes was in the increment of one pixel in radial direction of
artificial macropore. The resolution was assigned as 0.3375 mm/pixel and hydraulic gradient of
3.59×10-7 was applied which is similar to the typical simulations of the macroporous specimens.
Figure 6.10 presents the calculated ksat with varying diagonal distance of macropore. The
calculated ksat shows considerable decrease when a circle of pixels (shown in black color in
Figure 6.10) removed from the macropore (shown in gray color in Figure 6.10). This kind of
reduction in macropore area can occur due to swelling during laboratory measurement of k sat, the
due to the effect of X-ray CT scan resolution, and due to damage of specimen during
transportation form X-ray CT scan facility to laboratory to measure the ksat.

139

10.0

Width of black
area = 1 pixel

1.0

Simulated k

sat

using LBM (mm/s.)

100.0

0.1
Vertical Cylinder
dx=0.3375
Hydraulic Gradient = 3.59E-07
0.0
0

0.5

1

1.5

2

2.5

Diagonal Distance of Macropore (mm)

Figure 6.10: Effect of decrease or increase in macropore width in one pixel increment.

6.4 Effect of Morphology and Tortuosity of Macropore on Macropore Flow
In order to study the effect of macropore morphology on the overall flow regime in cracked
clays, several artificial macropore shapes were created in 3D domain. Synthetic 3D images were
developed using four different macropore cross sectional shapes: circle, diamond, ellipse, and
rectangle (see Figure 6.11 which shows the x-y plane), and three tortuosity levels in x-z plane
(see Figure 6.12). These shapes were created such that all the shapes have the same cross
sectional area. Two different set of shapes were used with different cross sectional areas: (i)
MP1= Macropore-1 and (ii) MP2= Macropore-2. The MP1 had an area of 0.784 mm2 and
macropore volume percentage of 23.9% (with respect to the overall domain size) and the MP2
140

had an area of 0.018 mm2 and macropore volume percentage of 5.2%. For each set of shapes, the
simulation domain area and volume were maintained constant for all the simulations.

Figure 6.11: Sectional shapes of macropores used to generate 3D image: Circle (a); diamond (b);
ellipse (c); and rectangle (d).

Figure 6.12: Tortuosity parameters and different tortuosity in same size of domain

Tortuosity is typically defined as the total length of the flow path divided by the shortest
distance between inlet and outlet as follows (see Figure 6.12):

141

∑

where L is the z-direction (vertical) distance between two macropore ends;

Eq. 6.16
is the length of the

ith segment of the macropore; n is the number of segment between the ends of the macropore;
and Le (effective distance) is the summation of li. Macropores with different tortuosity and the
tortuosity parameters are shown in Figure 6.12. The angle of macropore segment

was varied

from 18.4o to 45o. The increase in T has two important consequences: (1) the hydraulic gradient
on fluid decreases due to increase in effective distance, (2) change in the velocity direction from
the direction of gradient. The combined influence in flow due to these two effects is represented
by a tortuosity factor T2 in capillary tube models (Carman 1977; Bear 2013; Corey 1977). Since
gravity flow (unit gradient) is simulated in this study, increase in vertical distance alone will not
make change in flow rate. Hence, numerical simulations were carried out with inclined (angle α
from vertical) straight macropores and inclined bended macropores to study the effect of
inclination and multiple bends on macropore flow rates. Macropore shape and tortuosity were
varied systematically to investigate the influence of shape, perimeter of shape, aspect ratio of
shape, and tortuosity in macropore flow.
Table 6.1 presents macropore parameters such as dimensions of the shape (D1 and D2),
perimeter of the shape, aspect ratio (

⁄

), angle of macropore segment (α), tortuosity, and

macropore area used for numerical simulations. Pixel resolution (dx) and time step (ts) used for
MP1 simulations were 1.0×10-2 mm/pixel and 1.0×10-5 sec/ts, respectively. Similarly, for MP2
simulations the dx and ts were set to 5.0×10-3 mm/pixel and 2.5×10-6 sec/ts, respectively. For
MP1 and MP2 simulations the ns was set to 0.9 and 0.69, respectively. The ns was set such that
both the set of simulations give same micropore hydraulic conductivity of 5.45×10 -3 mm/sec in
142

physical units.

Gravity force of 9.81 m/s2 was applied. Periodic boundary condition was

implemented in all the boundaries. Relaxation parameter was set to 0.8 for all the simulations.
The MP2 simulation series from T2 to T4 were repeated with straight inclined macropores in
order to quantify the influence of inclined flow paths on macropore flow without the effect of
bends. During the steady gravity flow (fully developed laminar flow) the length of inclined
macropore does not influence the flow rate. Only the angle of inclination causes change in
hydraulic gradient on flow through inclined direction for gravity flow. Simulations were carried
out on high performance computing center (HPCC) of Michigan State University (MSU). A
typical simulation took about 10 hours of computer time using single processor.
Macropore flow was not affected significantly by incorporating micropore flow because the
simulations were carried out with a very small micropore hydraulic conductivity (5.45×10 -3
mm/sec). The smaller micropore hydraulic conductivity was selected in order to represent clay
soil. The increase in maximum velocity and macropore flow rate for cylindrical macropore due
to incorporating micropore flow was only 0.1% and 1%, respectively. However, incorporating
micropore flow allowed to simulate the actual physics at the boundary of micropore and
macropore region and to quantify the influence of the flow through micropore matrix compare to
macropore flow.

143

Table 6.1: Macropore parameters investigated
Simulation
ID

Macropore D1
D2
Aspect Perimeter α
Sectional
(mm) (mm) Ratio
(mm)
Shape

Tortuosity Macropore
(T)
Area
(mm2)

MP1_C1_T1 Circle

1.00

1.00

1.00

3.280

0

1.00

0.784

MP1_D1_T1 Diamond

1.50

1.04

1.44

3.825

0

1.00

0.784

MP1_E1_T1 Ellipse

1.50

0.66

2.27

3.691

0

1.00

0.784

MP1_R1_T1 Rectangle

1.00

0.78

1.28

3.517

0

1.00

0.784

MP2_C2_T1 Circle

0.15

0.15

1.00

0.497

0

1.00

0.018

MP2_D2_T1 Diamond

0.23

0.15

1.53

0.587

0

1.00

0.018

MP2_E2_T1 Ellipse

0.25

0.09

2.78

0.598

0

1.00

0.018

MP2_R2_T1 Rectangle

0.15

0.12

1.25

0.528

0

1.00

0.018

o

MP2_C2_T2 Circle

0.15

0.15

1.00

0.497

18.4

1.05

0.018

MP2_D2_T2 Diamond

0.23

0.15

1.53

0.587

18.4o 1.05

0.018

MP2_E2_T2 Ellipse

0.25

0.09

2.78

0.598

18.4o 1.05

0.018

MP2_R2_T2 Rectangle

0.15

0.12

1.25

0.528

18.4o 1.05

0.018

MP2_C2_T3 Circle

0.15

0.15

1.00

0.497

26.6o 1.12

0.018

MP2_D2_T3 Diamond

0.23

0.15

1.53

0.587

26.6o 1.12

0.018

MP2_E2_T3 Ellipse

0.25

0.09

2.78

0.598

26.6o 1.12

0.018

MP2_R2_T3 Rectangle

0.15

0.12

1.25

0.528

26.6o 1.12

0.018

MP2_C2_T4 Circle

0.15

0.15

1.00

0.497

45.0o 1.41

0.018

MP2_D2_T4 Diamond

0.23

0.15

1.53

0.587

45.0o 1.41

0.018

MP2_E2_T4 Ellipse

0.25

0.09

2.78

0.598

45.0o 1.41

0.018

MP2_R2_T4 Rectangle

0.15

0.12

1.25

144

0.528

45.0

o

1.41

0.018

6.4.1 Effect of Macropore Shape and Shape Parameters
In order to compare the magnitudes of flow rate through the micropores and macropores,
total flow rates through each zone was divided by their respective cross sectional area. In other
words, flow rate through the ‘white’ area in Figure 6.11 was divided by the cross sectional area
of the ‘white’ area to come up with a ‘normalized’ flow rate for macropore flow. Similarly, flow
rate through the ‘black’ area in Figure 6.11 was divided by the cross sectional areas of the
‘black’ portion of the overall cross section. Figure 6.13 shows the comparison of macropore flow
rate (per white area) to micropore flow rate (per black area) for two sets of macropores
simulated. As shown, the ratio of macropore flow rate per area to micropore flow rate per area
(RQ) varied from about 9800 to 13000 for MP1 macropores and it varied from 750 to 1100 for
MP2 macropores. In both cases the lowest ratio was relevant to the ellipse shape, which has the
highest aspect ratio.
Equivalent hydraulic conductivities (

) were calculated using Darcy’s law by dividing

the total flow rate at the bottom of the domain (flow rate through micropores + macropores) by
the total cross sectional area perpendicular to the gravity.

Then, the ratios of equivalent

hydraulic conductivities were defined as:
Eq. 6.17
where

is the equivalent hydraulic conductivity ratio.

145

RQ for MP1

(a) MP1: Percentage of Volume of Macropore - 23.9 %
1.4x10

4

1.2x10

4

1x10

4

8x10

3

6x10

3

4x10

3

2x10

3

MP1-T1

0
Circle

1.2x10

3

1x10

3

8x10

2

6x10

2

4x10

2

2x10

2

Diamond

Shape

Ellipse

Rectangular

(b) MP2: Percentage of Volume of Macropore - 5.2 %

RQ for MP2

MP2-T1

0
Circle

Diamond

Shape

Ellipse

Rectangular

Figure 6.13: Ratio of macropore flow rate per area to micropore flow rate per area of different
shape of macropores: MP1 (a); and MP2 (b).

Figure 6.14 presents
reduction in

of two set of macropore flow simulated. All the shapes show

as compared to circular shape and follow similar pattern for both the set of

macropore sectional areas (MP1 and MP2). Maximum reduction in
shape of MP2. Minimum reduction in

is about 35% for ellipse

is about 10% for rectangular shape of MP2. This
146

shows that the assumption of circular crack alone (not including other factors such as the
tortuosity) leads to errors between 10% and 35%. It should be noted that the shapes analyzed in
this paper (ellipse, diamond, rectangle) are simpler than the actual crack morphology in the field.
Therefore, more error should be expected because of complex angular cross section of the
cracks.
In order to analyze the influence of macropore parameters on the flow rate, ratios of
macropore flow rate were plotted against perimeter and aspect ratio of a given macropore cross
section. Ratio of macropore flow rate and ratio of perimeter were defined as:
Eq. 6.18

Eq. 6.19
where RQ is the flow rate ratio and RP is the perimeter ratio. Figure 6.15 presents the variation of
ratio of macropore flow rate with ratio of perimeter. Ellipse shape shows maximum reduction of
26% and 35% in flow rate for MP1 and MP2, respectively.

147

Normalized Equivalent
Hydraulic Conductivity (R

Keq

)

1.4
1.2

MP1-T1

MP2-T1

Area of Macropore:
2

MP1 = 0.784 mm

2

MP2 = 0.018 mm

1
0.8
0.6
0.4
0.2
0
Circle

Diamond
Ellipse
Shape of Macropore

Rectangular

Figure 6.14: Effect of macropore shape on the equivalent hydraulic conductivity

Generally there is a trend between RQ and RP, except for the diamond shape. Increase in
perimeter imply the increase in wall surface area, hence increase in wall friction. The reason for
higher flow rate of diamond shape compare to ellipse (even though perimeter length of diamond
is more than the perimeter length of ellipse) is that the diamond shape has larger eccentricity
(i.e., distance from its centroid to closest macropore wall is larger) than the ellipse shape. The
eccentricity (from the macropore wall) can be linked to the aspect ratio (see the D1 and D2
parameters in Table 6.1). The smaller aspect ratio represents larger eccentricity for the shapes
and shape parameters used in this study. The effect of viscous drag increases with decreasing
eccentricity (Panton 2006). In other words, as the centroid of the cross sectional area gets closer
to the wall, viscous drag increases, decreasing the overall flow rate (Panton 2006).

148

Area of Macropore:

MP1-T1
MP2-T1

Q

Ratio of Macropore Flow Rate (R )

1.1

2

MP1 = 0.784 mm

2

1

MP2 = 0.018 mm

0.9
Diamond
Rectangle
0.8

Ellipse

0.7

Ellipse
0.6
1

1.05

1.1

1.15

1.2

1.25

Ratio of Perimeter (R )
P

Figure 6.15: Effect of perimeter of macropore on macropore flow rate

The influence of aspect ratio on macropore flow rates (RQ) of macropores with same
⁄

sectional area but different aspect ratios (i.e.,

) is shown in Figure 6.16. The flow rate

decreases consistently with increasing aspect ratio. This is meaningful because as two solid (or
low permeable) sides get closer and closer, the overall velocity will decrease because of wall
friction. Both MP1 and MP2 follow the same line of reduction, indicating that this behavior is
not a function of the cross sectional area of the pore. The maximum reduction of 35% in flow
rate can be noticed for the maximum aspect ratio of 2.78 (MP2 ellipse).

149

Q

Ratio of Macropore Flow Rate (R )

1.1
MP1-T1
MP2-T1

1

Area of Macropore:
2

MP1 = 0.784 mm

2

MP2 = 0.018 mm

0.9

0.8

0.7

0.6
1

1.5

2

2.5

3

Aspect Ratio of Macropore Section (D1/D2)
Figure 6.16: Effect of aspect ratio of macropore on macropore flow rate

Figure 6.17 shows the 3D visualization of absolute value of the velocity distribution for
different cross sectional shapes. All of the subfigures are the simulations for tortuosity shape T4
illustrated in Figure 6.12. As shown in Figure 6.16, even though all the simulations were
performed in the same condition (i.e., gravity), the cross sectional areas are the same, the
velocity distribution exhibit significant difference for different cross sections. For example, the
ellipse shape has much lower velocity than the circle (and the others). The ellipse shape has the
highest aspect ratio and the flow rate through macropore is function of aspect ratio as explained
previously. It is noted that the scale of the colorbars for each subfigure is the same, where the
minimum (blue) shows 0 and the maximum (red) shows 7.28 mm/sec. The x-, y- and z- axis are
in lattice units and can be converted to physical units by multiplying each coordinate with the
150

resolution, which is 0.005 mm/lattice. The domain shown is 120x120x180 lattice units, which
corresponds to 0.60mm by 0.60mm by 0.90 mm domain size in physical units. The location of
maximum velocity is not at the center of the cross section in tortuous flow path, because the fluid
prefers to flow through the least resistance path (Clennell 1997).

Figure 6.17: 3D visualization of absolute value of velocity distribution for different cross
sections.

151

6.4.2 Effect of Tortuosity
Velocity vectors and streamlines are plotted in Figure 6.18 along the cross sectional plane
X-Z for circular cross section macropore. The solid nodes, velocity vectors (at every other node),
and streamlines (at interval of 4 nodes) are marked in gray color, red color arrows, and blue
lines, respectively. The zoomed rectangle shows the nature of the velocity profile along the
macropore of tortuosity of 1.41 (T4). The maximum velocity occurs away from the center line
(i.e., centroid of the cross section) of the macropore. The streamlines show that micropore flow
is minimal along convex side of the macropore bends and micropore flow prefers the concave
side of the macropore where less friction paths are found. The most efficient path for fluid flow
is not necessarily the shortest path (less tortuous) instead the path results in minimum energy
dissipation (Clennell 1997).
Figure 6.19 presents comparison of macropore flow rate through inclined straight
macropore and bended macropores. The flow rate through inclined straight macropores is
generally higher than the flow rate through bended macropores. The macropore flow rates of
inclined straight macropores and bended macropores were normalized with the vertical
macropore flow rate (i.e., when T=1). For a particular shape of macropore, the normalized
macropore flow rates are defined as:
Eq. 6.20

Eq. 6.21
where

is the normalized macropore flow rate for inclined straight macropores and

the normalized macropore flow rate for bended macropores.
152

is

T3

T2

T4

Figure 6.18: Velocity vectors and streamlines along the X-Z plane of circular cross section.
153

The

quantifies the influence of flow angle from the direction of gradient and the

influence of increase in effective length of flow path. The difference between

and

shows the effect of multiple turning or bends in flow path.

3

Macropore Flow Rate, Q (lu /ts)

3
Black color bars - Straight Inclined
White color bars - Bended

2.5
2
1.5
1
0.5
0

T2

T3

Circle

T4

T2

T3

T4

T2

T3

T4

Diamond
Ellipse
Shape of Macropore

T2

T3

T4

Rectangular

Figure 6.19: Comparison of flow rate through tortuous path and straight inclined path.

Figure 6.20 presents the normalized flow rates

and

with tortuosity. The

macropore flow rate decreases consistently with increasing tortuosity (increasing inclination
angle, effective length, and number of turnings/bends). The maximum effect of bends is about
25% on overall decrease of flow rate due to tortuosity for ellipse shape with three bents (T4).
Ellipse shape shows the maximum impact due to overall tortuosity and the macropore flow rate
reduces by about 70% for tortuosity of 1.41 (α = 45o). The lowest flow rate through ellipse shape
is in consistence with the straight vertical macropores which can be explained using aspect ratio.
154

Similarly, other shapes followed the same pattern in magnitude of flow rates through tortuous
macropores as straight vertical macropores. The effect of increase in effective length and
inclination angle shows more influence in decreasing flow rate than number of bends for the
tortuosity considered in this study. This shows that the tortuosity of a crack affects the flow rate
significantly and needs to be included in the analyses.
1.2

Area of Macropore:

Circle-R(QTS)
Diamond-R(QTS)
Ellipse-R(QTS)
Rectangle-R(QTS)
Circle-R(QTB)
Diamond-R(QTB)
Ellipse-R(QTB)
Rectangle-R(QTB)

QTS

QTB

Normalized Macropore Flow Rate
or R )
(R

2

MP2 = 0.018 mm
1

0.8

0.6

0.4

0.2
1

1.1

1.2

1.3

1.4

Tortuosity
Figure 6.20: Influence of tortuosity on macropore flow rate.

155

1.5

6.5 Development of a Simplified Predictive Equation for Macropore Flow
From the results presented on macropore flow rate through different shapes and tortuosity
(Figure 6.16 and Figure 6.20), we can come to a conclusion that the macropore flow rate is a
function of aspect ratio and tortuosity. Change in velocity direction, increase in effective length,
and bends/turnings are the influencing parameter in tortuosity. However, a particular pattern was
not noticed on the effect of bends for the tortuous paths considered in this study. Hence, an
equation can be formulated using the aspect ratio and tortuosity to predict the flow rate through
different shapes and tortuosity based on the flow rate of straight vertical cylinder. Such an
equation can be implemented in water balance models (e.g., RZWQM) in order to enhance the
macropore flow module to include different shapes and tortuous macropores without much
complication in programming. The flow rate of different shapes and tortuous macropore is
inversely proportion to the tortuosity and aspect ratio (Figure 6.16 and Figure 6.20). Carman
(1997) pointed out that the combined influence in flow due to change in velocity direction and
increase in effective length can be represented by the tortuosity factor T2. Figure 6.16 shows that
the effect of shape can be represented by

R A of the shape, where RA is aspect ratio. Hence,

following equation is formulated:
Qp 

where

Qs
T 2 RA

Eq. 6.22

is the predicted flow rate; and Qs is the straight vertical cylindrical macropore flow rate.

Figure 6.21 presents the predicted vs. simulated flow rates in lattice units. Predicted macropore
flow rates of straight inclined macropores match well with simulated macropore flow rate. For
the bended macropores, at lower tortuosity or lesser number of bends the prediction is
156

comparable with the prediction for straight inclined macropores. However, at higher tortuosity
(tortuosity more than about 1.3) the predicted flow rates are higher than the simulated flow rates
as expected. More detailed study is required to identify the influence of bends on macropore flow

3
Open Symbol - Straight Inclined Macropore
Solid Symbol - Bended Macropore

3

Predicted Macropore Flow Rate (lu /ts)

rate and to come up with a factor for prediction equation.

2.5
2
Circle
Diamond
Ellipse
Rectangle
Circle
Diamond
Ellipse
Rectangle

1.5
1
0.5

1:1 line

0
0

0.5

1

1.5

2

2.5

3

3

Simulated Macropore Flow Rate (lu /ts)
Figure 6.21: Comparison of simulated and predicted macropore flow rate.

6.6 Verification of the Simplified Predictive Equation for Macropore Flow
The prediction equation developed in section 6.6 was used to predict the ksat of macroporous
specimens using the X-ray CT scan images. The image analysis and macropore characterization
methodologies are presented in section 5.4 and section 5.5. The equivalent cylinders relevant to
157

minimum radii were used to predict the macropore flow rate through the specimens and hence
the ksat of specimens. The cross section of minimum radii slices are shown in Figure 6.22. The
parameters used to calculate the aspect ratio of the macropore shapes were marked in Figure
6.22. The tortuosity was calculated for each macropore according to the definition given in
section 6.4. The macropore shape parameters and tortuosity are listed in Table 6.2. Also, the
predicted ksat of macroporous specimens using equation Eq. 6.22 are presented in Table 6.2.
Calculated ksat using LBM are presented in Table 6.2 for comparison.
Figure 6.23 presents the comparison of calculated vs. predicted and measured ksat of
macroporous specimens. For specimens S1 and S2, the prediction was about 20% less than the
calculated ksat. However, for specimen S3 the prediction is slightly higher. Specimen S3 had the
tortuous macropore (tortuosity of 1.38) compare to specimen S1 and S2. The effect of number
bending on macropore flow is not included in the prediction equation. However, the predictions
of ksat using the equation for the macroporous specimens are comparable with the calculated ksat.
Hence, X-Ray CT images can be used to predict the ksat of the specimen using the prediction
equation.

158

S1

D1

D2

S2

D2

D1

S3
D1

D2

Figure 6.22: Cross section of macropore relevant to minimum equivalent radius (not to scale).
159

Table 6.2: Characteristics of macropore, and calculated and predicted ksat

Estimated
Macro- Min. k using
sat
Sample
porosity Radius Poiseuille Tortuosity
(-)
ID
law
(mm)
(%)
(mm/s.)

D1

Simulated
Predicted
Aspect ksat using
ksat using Eq.
LBM
D2 ratio
6.22 (mm/s.)
(mm/s.)

S1

6.84

3.340

118.4

1.10

8.78 5.06

1.49

98.9

80.3

S2

0.51

1.973

14.4

1.05

4.73 3.88

1.22

14.9

11.9

S3

0.28

1.966

14.2

1.38

2.70 1.20

2.25

3.9

5.0

100.0

Predicted or Measured k

sat

(mm/sec.)

S1

S2

Predicted

10.0

S3

Measured
1.0

0.1
0.1

1

10

Simulated k

sat

100

using LBM (mm/sec.)

Figure 6.23: Calculated vs. predicted and measured ksat of macroporous samples

160

6.7 Multiphase LBM
6.7.1 Shan-Chen (SC) Multiphase Model
The inclusion of an attractive force between neighbor particles, which leads to phase separation
is the principal characteristic that differentiate the Shan-Chen (SC) model from single phase
model (Sukop and Or 2004). The attractive force Fint is incorporated in LB algorithm as
proposed by Shan and Chen (1993). For the D2Q9 model the Fint is as follows:

Fint ( x, t )  G ( x, t )i Wi ( x  ei , t )ei

Eq. 6.23

where G is interaction strength constant that is negative for particle attraction; and the function

 is depend on density as follows:
 (  )   0 exp    0  



Eq. 6.24

where  0 and  0 are constants and the values of these constants can be varied arbitrarily (Sukop
and Thorne 2006).
The interaction forces (F) (fluid/fluid or fluid/solid) can be incorporated as velocity shifts as
follows (Sukop and Or 2004):

U u

F


Eq. 6.25

where U is the macroscopic velocity. The velocity U is used to compute fieq. The equations from
Eq. 6.1 to Eq. 6.6 and Eq. 6.24 to Eq. 6.25 form the single component multiphase (e.g., water
and vapor) model. This algorithm can be extended to multi component multiphase (e.g., oil and
161

water) model by adding individual distribution functions for each component ( f i where



indicates the component of the fluid). For multicomponent model, the f i (eq ) is calculated from
composite macroscopic velocity,
1

u

  

i



1

 

f i  ei

Eq. 6.26





Figure 6.24 presents the phase separation simulated using 2D single component multiphase
(liquid-vapor) SC model (blue and gray colors indicates water and vapor, respectively). A
domain of 200×200 pixel2 was used for the simulation with average initial density 200 in lattice
units and perturbed with random density of 1 in lattice units. Figure 6.24 shows that initially the
phase separation occurs rapidly and eventually it slows down due to the fact that large numbers
of small droplets start competing for the liquid particle in the vapor phase (Sukop and Thorne
2004). During the phase separation, interfaces formed to minimize their total length in 2D, a
consequence of free energy minimization and occurs in part by geometric rearrangement in to the
minimum surface area volume (a circle in 2D) (Sukop and Thorne 2004). Also, the smaller
droplets are more unstable compare to large droplets due to the surface tension effects on the
interface (Begum and Basit 2008).

162

ts =100

ts =500

ts =1000

ts =2000

ts =25000

ts =50000

Figure 6.24: Phase transition simulated using single component SC multiphase model

Equations Eq. 23 and Eq. 24 lead to a non-ideal equation of state (EOS) for SC model and
the EOS (SC model EOS) is given by,

p  cs2  

cs2
2
G  
2

Eq. 6.27

If there is no interaction force, that is function  reaches zero, the EOS (Eq. 6.27) reduces
to ideal gas EOS (Eq. 6.7). EOS defines the relationship between pressure, temperature and
volume (or density). In SC model EOS, the temperature is not explicitly defined but the
temperature is dictated by G.

163

The SC model EOS is qualitatively similar to the van der Waals (vdW) EOS. However, the
lack of a repulsive force in SC model EOS has led to a behavior where liquid phase is actually
more compressible than the vapor phase which is generally not physically correct (Sukop and
Throne 2006). Also, the spurious velocity (unphysical velocities occurs in two phase LBM at the
interface region) is a challenge in SC model since the large spurious velocity will make the
simulation unstable (Yuan and Schaefer 2006). The spurious velocity limits the density ratio
since increasing density ratio results in an increase in the spurious velocity (Yuan and Schaefer
2006). Yuan and Schaefer (2006) recommended that spurious velocity can be reduced and higher
density ratio can be achieved by changing the EOS. The EOS can be changed by changing the
form of function . The function  can be derived from Eq. 6.27 as follows:



2 p  cs2 
 ( ) 
cs2G



Eq. 6.28

Yuan and Schaefer (2006) studied different form of  namely the vdW EOS, CarnahanStarling (C-S) EOS, Redlich-Kwong (R-K), Redlich-Kwong Soave (RKS), and Peng-Robinson
(P-R). Huang and Lu (2009) used R-K EOS to study the relative permeability in steady- state
gas/liquid flow in porous media. Similar approach is followed in this study to simulate
percolation of water through unsaturated macropores. The R-K EOS is defined as,

RT
a 2
p

1  b
T 1  b 

Eq. 6.29

where p is the pressure; R is the gas constant; T is the temperature; a is the attraction parameter;
and b is the repulsion parameter. The function  for R-K EOS can be derived by replacing p in
Eq. 6.28 according to Eq. 6.29.
164

The adsorption force between fluid particles and solid surfaces ( Fads ) is introduced as
follows (Huang and Lu 2009):

Fads ( x, t )  G  ( x, t )i Wi  w s( x  ei )ei

Eq. 6.30

where  w is density of solid phase; and s is equal to 1 or 0 for a solid or pore respectively. The

 w is not true density of solid phase but it is a parameter used to tune different wall properties
(Huang and Lu 2009). When the  w changes between density of liquid (  l ) and density of air (

 a ), the contact angle changes between 0o to 180o.
Different contact angles were simulated using a domain of 200×200 pixel2. A rectangle fluid
was placed on the solid at the center of the domain. The parameters required for R-K EOS were
assigned according to Huang and Lu (2009). The parameters a and b were set to 2/49 and 2/21,
respectively. And, R is 1 and T is a function of Tc, T=0.85Tc where Tc is 0.1961. The liquid and
gas densities were set to 6.06 and 0.5, respectively. Figure 6.25 presents the three different
contact angles of 12o, 75o, and 180o simulated using 5.5, 3.0, and 0.8, respectively, for  w . Blue,
gray, and yellow colors indicate liquid, gas, and solid phase, respectively. The dx and dt were set
to 1.0 mm/pixel and 0.1667 sec./ts. The relaxation parameter was assigned as 1.0. The
simulations reached steady-state in about 3000 ts.

165

Θ=180o

v1
0.5
0.12

Θ=75o

50

0.4
max velocity, lu/ts

0.1
0.08

100
0.06
0.04

150

0.3

0.2

0.1
0.02
200
50

100

150

0

200

6

Θ=12o

5
50

0

3.5
3
2.5

4
2
100

3
1.5
2

150

200
50

100

150

200

1

0.5

0

0

Figure 6.25: Simulated different contact angles (color bar shows the density of fluid)
166

1

0

500

Figure 6.26 presents the simulation of infiltration of liquid through unsaturated macroporous
soil matrix using 2D single component multiphase SC model. The soil matrix was considered
impermeable. Hence, water can flow only through unsaturated macropores. Blue, gray, and
yellow colors indicate liquid, gas, and solid phase, respectively. The simulation parameters were
assigned similar to the previous simulations (contact angle simulations). The density ratio was
about 12. The contact angle between liquid and solid was set to 180o. Velocity boundary
condition was assigned at top and bottom boundary. Periodic boundary condition was assigned
for east and west boundary. The liquid infiltration during four time steps was shown in Figure
6.26. The liquid flows through bigger pores initially and then flow occurs through some smaller
pores. The model was able to simulate the flow of liquid qualitatively correctly. However, model
was unstable after 20000 ts. Also, model was unstable for gravity flow simulations.
ts =2000

ts =12000

ts =6000

ts =18000

Figure 6.26: Simulation of infiltration of liquid through macropores using SC model
167

6.7.2 Color Indices (CI) Model
Misici and Palpacelli (2005) presented a two-fluid LBM based on color indices. In their
algorithm, a color distribution function evolves together with the standard density distribution
function. Except the evolution of the color distribution function, the rest of the algorithm is same
as the single phase standard LBM algorithm (previously presented). Hence, only the evolution of
color distribution function is described in this section. The presence of two fluids (red and blue
fluid) is described by an indicator function  :


1 if   1

  0 if   0
 sin  0.5   1
i

otherwise
2


Eq. 6.31

where  is color distribution. The interface of the two fluids is those nodes for which
0   ( x, t )  1 and one of the two fluids is at nodes for which  ( x, t )  1and the other fluid is at

nodes for which  ( x, t )  0 . The color distribution (  ) is calculated as follows:

  i  i x, t 

Eq. 6.32

The color distribution function is defined as:

 i x, t    i x  ei , t  1    i x  ei , t  1   ieq x  ei , t  1

Eq. 6.33

where  is relaxation parameter of color distribution function and it is function of relaxation
parameters of individual fluids (  r for red fluid and  b for blue fluid):
168

  r  b 1   

Eq.6.34

The equilibrium color distribution function is given by:

 ieq 

  x, t  f i  x, t 


Eq. 6.35

Finally the surface tension term is added to the f i eq (of the standard LBM) in the following
form for a 2D model (the general term can be found in Misici and Palpacelli, 2005):
2
9  Gx  Gy 
Si  Wi  
 2 G
2 
G





 e  e 2  2 
iy
 ix
3 



Eq. 6.36

where σ is the surface tension strength parameter; and G is the gradient of the  .
Figure 6.27 presents the simulation of infiltration of water through macropores using 2D
color indices model. Actual values of density and kinematic viscosity of water and air (blue and
gray color in Figure 6.27, respectively) were used as input parameters for the simulation. Also,
the real gravitational acceleration was applied. Due to the instability of the model, very small
time step was used (dt = 1×10-6 sec/time step) which caused longer computer time. Regardless of
the very small time step used, the model showed instability near about 2.5×105 ts (Figure 6.27)
and the water area showed decreasing pattern. Refinement of the code is required for further
study.

169

ts =0

ts =150000

ts =250000

ts =300000

Figure 6.27: Simulation of infiltration of water through unsaturated macropores using CI model.

170

CHAPTER 7: SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS FOR
FURTHER STUDIES

7.1 Summary
The long-term change in hydraulic properties of earthen landfill final covers is one of the major
challenges in cover design as it impacts the long-term percolation into the waste and impacts the
long-term risk for ground water and surface water contamination. Formation of macropores and
resulting flow through macropores increases the percolation through earthen landfill covers
during service. Water flow in micropores can be downward or upward depending on the
hydraulic gradient due to the soil capillarity and gravity. However, water flow in macropores is
predominantly downward. Presently available models used for water balance modeling of
earthen covers are based on Richards’ equation that simulates only the micropore flow (capillary
flow). Hence, a model that can simulate micropore flow as well as macropore flow is required to
simulate long-term hydrology of earthen covers. The key objectives of this dissertation were to:
(1) quantify micropore flow vs. macropore flow through an engineered earthen landfill final
cover; (2) identify and validate numerical model capable of simulating micropore flow as well as
macropore flow through earthen landfill covers; and (3) develop a model based on lattice
Boltzmann method that can simulate saturated flow through micropores as well as macropores.

171

7.1.1 Field Data
A field-scale test section of an earthen cover was constructed at a landfill located near
Detroit, Michigan and monitored for about four years. The cover consisted of 1.5 m thick
compacted clay overlain by 0.3 m thick topsoil. The test section was instrumented to measure
water contents, water potentials and percolation. Measured soil water storage (SWS) and
percolation were analyzed for influence of macropore flow. The SWS calculated using point
measurement of water content could not show the influence of macropore flow. Hence, SWS
alone cannot be used to identify the presence of macropore flow in earthen covers.
Measured annual percolation for the first year was about 7 mm whereas the annual
percolations for the second and the third years were 75 mm and 136 mm, respectively. Annual
precipitations recorded for the 1st, the 2nd, and the 3rd years were 792 mm (35 percentile), 1,104
mm (98 percentile), and 828 mm (48 percentile), respectively. While the increase in annual
precipitation and higher winter and spring precipitation may have contributed to higher
percolation in Years 2 and 3, formation of macropores also contributed to the increase.
The test section was subjected to controlled irrigation during the 1st year and the 4th year of
service to evaluate macropore flow. The field test section reached similar degrees of saturation at
the end of each of the two irrigation tests. While the first irrigation test did not show distinct
breakthrough, the second test resulted in breakthrough within few hours confirming macropore
flow. Hence, controlled irrigation coupled with instrumentation can be used to assess evolution
of macropore flow. Percolation measured during the four years of service was used to estimate
effective hydraulic conductivity (keff) of the test section. The average degree of saturation of the
storage layer was about 80-87% during the observed breakthroughs. The keff increased by about
172

an order of magnitude during the four years of service. The increase in keff is an indication of the
increase in field ksat. Formation of macropores caused the increase in keff.
Methane (CH4) injection tests were carried out in the field to evaluate CH4 oxidization
capacity of earthen covers. The results aided in identification of macropores and their locations
and to distinguish the active and non-active macropores. The width of the macropores ranged
from few mm to about 4 cm and the length of the continuous segments of the macropores that
could be measured from the surface were from 1 cm to 60 cm. Only one visible surface crack
was identified in the lysimeter area. The average measured CH4 flux through macropores was
about two orders of magnitude higher than the average CH4 flux through intact ground surface.
The CH4 oxidization capacity of cover where a macropore present was almost zero while the
CH4 oxidization capacity of intact cover ranged from 0% to 30%. Two orders of magnitude
higher CH4 flux and almost zero CH4 oxidization capacity indicate that those macropores had
reached the lysimeter.

7.1.2 Water Balance Modeling
Field evaluation of RZWQM was carried out using the data collected from the field test
section. Numerical predictions using commonly used numerical model UNSAT-H were also
compared with RZWQM predictions. For the first year, both models simulated cumulative
percolation and SWS relatively accurately. However, both models are not capable of simulating
frozen ground and hence SWS predictions of both models were not as accurate during winter and
spring months. The numerically predicted annual percolation using UNSAT-H and RZWQM
173

were about 50% of measured annual percolation during the second year.

The measured

percolation during the second year was affected by macropore flow. UNSAT-H does not
simulate macropore flow and the macropore flow component was not activated in RZWQM in
the initial runs. UNSAT-H and RZWQM (without macropore flow) were able to simulate about
90% and 70% of cumulative percolation for the third year, respectively. However, the simulated
percolation had a delay of about four months compared to measured percolation. In the field, the
onset of percolation was four months early due to the relatively rapid flow through the
macropores.
The macropore parameters required for the macropore flow component of RZWQM were
calibrated using the flow rate measured during the second irrigation test and the number of
macropores identified using methane injection tests. Calibrated macropore parameters based on
the results of the second irrigation test and methane tests were used to simulate percolation
during the second and the third year of service. The micropore and macropore flow simulations
using RZWQM produced better percolation prediction compared to micropore only simulations
for the second year, but about 13% less than the measured cumulative percolation. However, the
micropore and macropore flow simulations using RZWQM for the third year simulated
cumulative percolation that matched very closely with measured cumulative percolation. The
size of the macropore can change in the field due to swelling or shrinkage of clay, plant roots
penetration and death, and animal activity. Hence, periodic field calibration of macropore flow
parameters may yield more accurate predictions of macropore flow.

174

7.1.3 Formation of macropore in Laboratory Specimen
Experimentally simulating macropore flow through cracks in clay in the laboratory is
challenging due to the self-healing nature of cracks in clay and the effect of setup walls on the
formation of cracks. A new technique was developed to create macropores in compacted clay in
the laboratory. Food grade sugar in crystal form wrapped into an insect mesh was placed within
clay specimen during compaction. During permeation, sugar dissolves while leaving a void
which is a macropore. Control tests were carried out to test the effect of sugar on hydraulic
conductivity of compacted clay specimen. Double ring rigid wall permeameter was used to
measure the hydraulic conductivity of compacted specimens. The macropore resulted in an
increase in the hydraulic conductivity by at least four orders of magnitude. High resolution Xray CT scans of compacted clay specimens were used to develop 3-D structure of macropore(s).
Macroporosity of macropores were estimated using the X-ray CT images.

7.1.4 Lattice Boltzmann Model (LBM)
A 3-D LBM simulating saturated flow through micropores and macropores was developed.
Verification of the LBM was carried out using analytical solutions. Simulated ksat of
macroporous specimens were about an order of magnitude higher than the laboratory
measurement of ksat. The X-ray CT scans were performed on partially dry specimens. Unlike
macropores in macroporous rock or compacted aggregate, the size of macropores within
compacted clay can change due to swelling of clay during wetting. Hence, saturated samples
need to be subjected to X-ray CT scans in order to develop structure of macropores.
175

Influence of macropore morphology defined by shape and shape related parameters, and
tortuosity were systematically investigated via numerical simulations (3-D LBM). Constant cross
sectional area and volume of macropores with various shapes were considered. The influence of
turnings/bends in flow paths, increase in effective length of the flow path, and inclination of flow
paths were analyzed. It was observed that the macropore flow rate decreases by about 10% to
70% as compared to a circular cross section when the tortuosity is increased by 40% or when the
shape of the cross section is altered. Increase in aspect ratio of sectional shape shows decrease in
macropore flow rate. Tortuosity plays a significant role in macropore flow where the macropore
flow rate reduces by about 70% for tortuosity of 1.41 as compared to a tortuosity of 1. About
25% of overall reduction in flow rate was due to bends in flow path for tortuosity of 1.41.
A prediction equation was formulated to estimate the flow rate of arbitrary shape and
tortuous macropore using the flow rate of straight vertical cylinders. The predicted flow rate of
different shapes and tortuous macropores based on straight cylindrical macropore gives good
match with simulated flow rate using LBM. However, more detailed study is required to predict
the flow rate of bended macropores with tortuosity higher than about 1.3. The prediction
equation was used to predict the ksat of macroporous specimens. Equivalent circular area of
macropores on each slices of X-ray CT images of macroporous compacted clay specimens were
calculated. The equivalent cylinders relevant to minimum radii were used to predict the
macropore flow rate through the macroporous specimens and hence the ksat of the specimens.
The macropore shape parameters required for prediction equation were estimated from the cross
section of macropore on minimum radius slice for each specimen. Also, the tortuosity was
calculated for each macropore from the 3-D structure of the macropores. Comparison of the
predicted ksat using the proposed formulation and simulated ksat using the LBM were showed
176

excellent match. Hence, X-Ray CT images can be used to predict the ksat of the macroporous
specimen with the help of the proposed prediction equation.

7.2 Conclusions
The key conclusions of this study are:
1. Formation of macropores alters the long-term hydraulic properties of earthen landfill
final covers. Field irrigation tests and estimation of keff indicated that the saturated
hydraulic conductivity of the field-scale earthen final cover test section increased by
about an order of magnitude due to the formation of macropores over the four years of
service life. Controlled irrigation tests coupled with instrumentation can be used to
quantify macropore flow through earthen landfill final cover.
2. RZWQM developed by Agricultural Research Services (ARS), USDA was evaluated to
simulate hydrology of earthen landfill final cover using the data collected from the fieldscale test section prior to the formation of macropores and RZWQM results were
compared with commonly used micropore flow model UNSAT-H. The simulation results
showed that RZWQM can be used to model long-term hydrology of earthen landfill final
covers.
3. Macropore flow through earthen landfill final cover can be modeled using RZWQM as
long as site specific macropore parameters are input. Macroporosity or number of
macropores and average diameter of macropores are the important calibration parameters
require to model the macropore flow.
177

4. Size of macropore can change in the field due to swelling and shrinkage, animal activity,
and plant root penetration. Hence, measurement or calibration of macropore parameters
needs to be done periodically.
5. Laboratory experiments on macropore flow through cracks in clay is challenging due to
the self-healing nature of the cracks. Hence, a new technique was developed to introduce
continuous macropores within lab-scale compacted clay samples and the technique
yielded consistent results.
6. 3-D LBM developed in this study can be used to simulate saturated flow through
micropores as well as macropores. Verification of the model confirmed the accuracy of
the model. Laboratory validation of the model was difficult due to the swelling nature of
the compacted clay specimen.
7. The LBM aided in quantification of the effect of macropore parameter and tortuosity in
macropore flow. Tortuosity and aspect ratio are the two key parameters that influence
flow through macropores of the same cross sectional area. The macropore flow rate
reduced by about 35% and 70% when the aspect ratio and tortuosity were increased to
about 2.8 and 1.4 from 1, respectively.
8. The prediction equation developed using LBM simulations on the effect of macropore
parameter and tortuosity can be used to predict the macropore flow rate as long as scan
images of macropores are available. The prediction equation requires only two
parameters, aspect ratio and tortuosity, other than the sectional area of macropore.
9. The prediction equation can be implemented in water balance models (e.g., RZWQM) in
order to enhance the macropore flow module to include different shapes and tortuous
macropores without much complication in programming.
178

7.3 Recommendations for Further Studies
The macropore flow component of RZWQM was calibrated using the irrigation test performed
on the test section during the 4th year of service. The simulated percolation using the calibrated
macropore parameters provided excellent match with measured percolation for the third year
field data but not for the second year field data. The size of the macropore can change in the field
due to swelling of clay, plant roots penetration and death, and animal activity. Hence, in the
further studies, it is recommended to perform controlled irrigation tests during different time
spans in service to create a data base for macropore parameters. Also, the calibrated macropore
parameters can be verified to simulate the macropore flow through earthen landfill final covers at
different climate zones in future studies.
The 3-D LBM was verified using analytical solutions. High resolution X-ray CT scans were
performed on compacted clay specimens to develop 3-D structure of macropore and to validate
the LBM using the measured ksat. Unlike macropores in macroporous rock or compacted
aggregate, the size of macropores within compacted clay can change due to swelling of clay
particles. Hence, saturated samples needs to be subjected to X-ray CT scans in order to develop
structure of macropores and to validate the LBM. In future studies, it is recommended to use
improved permeameter to measure the ksat of macroporous compacted clay samples and to
perform X-ray CT scan on saturated samples in order to validate the LBM to simulate macropore
flow through compacted clay samples.
The prediction equation developed to estimate the flow rate of arbitrary shape and tortuous
macropore using the flow rate of straight vertical cylinders was able to predict the flow rate
through macroporous compacted clay samples. The effects of sectional shape and tortuosity are
179

included in the equation. The effect of tortuosity can be due to change in velocity direction,
increase in distance, and bends/turnings. However, the effect of bends/turnings is not considered
in the equation. The LBM simulations showed that the effect of bends/turnings in macropore is
small compare to the effect of change in velocity direction and increase in distance. However, it
is recommended to study the effect of bends in detail in future studies. Also, the prediction
equation can be validated using samples collected from field in further studies.

180

REFERENCES

181

REFERENCES

Adu-Wusu, C. and Yanful, E. K. (2006), “Performance of engineered test covers on acidgenerating waste rock at Whistle mine, Ontario,” Can Geotech J 43:1–18
Ahuja, L. R., DeCoursey, D.G., Barnes, B.B. and Rojas, K.W. (1993), “Characteristics of
macropore transport studied with the ARS Root Zone Water Quality Model,”
Transactions of the ASAE, 36(2):369-380.
Ahuja, L.R., Johnson, K.E. and Heathman, G.C. (1995), “Macropore transport of a surfaceapplied bromide tracer: Model evaluation and refinement,” Soil Science Society of
America Journal, 59: 1234–1241.
Ahuja, L.R., Rojas, K.W., Hanson, J.D., Shaffer, J.J. and Ma, L. (Eds.) (2000), “The Root Zone
Water Quality Model,” LLC, Highlands Ranch, Colorado: Water Resources Publications.
Albrecht, B. and Benson, C. (2001). “Effect of desiccation on compacted natural clays.” J.
Geotech. & Geoenviron. Engrg., Vol. 127(1), 67–75.
Albright, W., Benson, C., Gee, G., Abichou, T., Roesler, A., and Rock, S. (2003), “Evaluating
the Alternatives”. Civil Engineering, 73(1), 70-75.
Albright, W., Benson, C., Gee, G., Abichou, T., McDonald, E.V., Tyler, S., and Rock. S.
(2006a). “Field performance of a compacted clay landfill final cover at a humid site.” J.
Geotech. & Geoenviron. Engrg., 132(11), 1393–1403
Albright, W., Benson, C., Gee, G., Abichou, T., Tyler, S., and Rock. S. (2006b). “Field
Performance of Three Compacted Clay Landfill Covers.” Vadose Zone Journal, 5, 11571171.
Albright, W., Benson, C., Gee, G., Roesler, A., Abichou, T., Apiwantragoon, P., Lyles, B. and
Rock, S. (2004), “Field water balance of landfill final covers,” Journal of Environmental
Quality, 33(6): 2317-2332.
Albright, W., Benson, C., Waugh, W.J. (2010). “Water balance covers for waste containment:
Principles and practice,” ASCE press, Reston, Virginia.
Apiwantragoon, P. (2007). “Field hydrologic evaluation of final covers for waste containment.”
Ph.D. dissertation, University of Wisconsin, Madison, Wisc.
ASTM D698. (2007), “Standard Test Methods for Laboratory Compaction Characteristics of
Soil Using Standard Effort,” West Conshohocken, Pennsylvania: ASTM International

182

ASTM D5084. (2010a), “Standard Test Methods for Measurement of Hydraulic Conductivity of
Saturated Porous Materials Using a Flexible Wall Permeameter,” West Conshohocken,
Pennsylvania: ASTM International
ASTM D5126. (2010b), Standard Guide for Comparison of Field Methods for Determining
Hydraulic Conductivity in Vadose Zone, West Conshohocken, Pennsylvania: ASTM International
Attar, E. and Korner, C. (2009), “Lattice Boltzmann method for dynamic wetting problems,”
Journal Colloid and Interface Science, 335, 84-93.
Bartling, P.N.S., Ma, L. and Ahuja, L. R. (2012), “RZWQM2 User Guide,” USDA-ARS
Agricultural Systems Research Unit, 2150 Centre Ave., Building D-200, Fort Collins,
CO, USA.
Bear, J. (1972), “Dynamics of Fluids in Porous Media,” Elsevier. New York.
Begum, R. and Basit, M.A. (2008), “Lattice Boltzmann method and its applications to fluid flow
problems,” European Journal of Scientific Research, 22(2), 216-231.
Benson, C.H., Abichou, T., Albright, W.H., Gee, G.W., and Albright, W.H. (2001), “Field
evaluation of alternative earthen final covers,” International Journal of
Phytoremediation, 3(1): 105–127.
Benson, C., Albright, W., Roesler, A., and Abichou, T. (2002), “Evaluation of Final Cover
Performance: Field Data from the Alternative Cover Assessment Program (ACAP),”
Proceedings of Waste Management 2002. Tucson, AZ.
Benson, C., Bohnhoff, G., Ogorzalek, A., Shackelford, C., Api-wantragoon, P. and Albright, W.
(2005), “Field data and model predictions for an alternative cover,” Waste Containment
and Remediation, GSP No. 142: 1-12. Reston, Virginia: ASCE
Benson, C., Bosscher, P., Lane, D., and Pliska, R. (1994), “Monitoring system for hydrologic
evaluation of landfill final covers,” Geotechnical Testing Journal, 17(2), 138 – 149.
Benson, C.H. and Daniel, D. (1990). “Influence of clods on the hydraulic conductivity of
compacted clay.” J. Geotech. Engrg., Vol. 116(8), 1231–1248.
Benson, C. H. and Gribb, M. M. (1997). “Measuring unsaturated hydraulic conductivity in the
laboratory and field.” Proceedings of sessions on unsaturated soils at Geo-Logan 1997,
July 15-19, Logan, Utah.
Benson, C. and Khire, M. (1995). “Earthen Covers for Semi-Arid and Arid Climates.” Landfill
Closures. ASCE GSP No. 53, 201-217.
Benson, C.H. and Othman, M. (1993). “Hydraulic conductivity of compacted clay frozen and
thawed in situ.” J. Geotech. Engrg., Vol. 119(2): 276–294.
183

Benson, C., Sawangsuriya, A., Trzebiatowski, B., and Albright, W. (2007), Post-construction
changes in the hydraulic properties of water balance cover soils, Journal of Geotechnical
and Geoenvironmental Engineering, 133(4): 349–359.
Benson, C. and Wang, X. (2006), “Temperature-compensating calibration procedure for water
content reflectometers,” Proc. TDR 2006, Purdue University, West Lafayette, USA, Sept.
2006, Paper ID 50.
Beven, K.J. and Germann, P.F. (1982), “Macropores and water flow in soils,” Water Resources
Research, 18: 1311–1325.
Bhatnagar, P.L., Gross, E.P., Krook, M. (1954), “A model for collision processes in gases. I.
Small amplitude processes in charged and neutral one-component systems,” Physical
review, 94(3):511.
Bohnhoff, G., Ogorzalek, A., Benson, C., Shackelford, C. and Apiwantragoon, P. (2009), “Field
data and water balance predictions for a monolithic cover in a semiarid climate,” Journal
of Geotechnical Engineering, 135(3): 333-348.
Boynton, S.S. and Daniel D.E. (1985), “Hydraulic conductivity tests on compacted clay,” J.
Geotech. Engrg., Vol. 111(4): 465–478
Brinkman, H.C. (1947), ‘‘A calculation of the viscous force exerted by a flowing fluid on a
dense swarm of particles,’’ Appl. Sci. Res., Sect. A1, 27.
Brooks, R. and Corey, A. (1964), “Hydraulic properties of porous media.” Hydro. Paper No. 3,
Civ. Engrg. Dept., Colorado State Univ., Fort Collins, Colorado.
Carman P.C. (1937), “Fluid flow through granular beds,” Transactions of the Institute of
Chemical Engineers, 50:150-166.
Cancelliere, A., Chang, C., Foti, E., Rothman, D. H. and Succi, S. (1990), “The permeability of a
random medium: comparison of simulation with theory,” Physics of Fluids A, 2,
2085−2088.
Celia, M. A., Bouloutas, E. T. and Zarba, R. L. (1990), “A general mass-conservative numerical
solution for the unsaturated flow equation,” Water Resour. Res. 26:1483-1496.
Chanton, J., Powelson, D., Abichou, T., and Hater, G. (2008). “Improved field methods to
quantify methane oxidation in landfill cover materials using stable carbon isotopes.”
Envir. Sci. & Tech., 42(3), 655-670.
Chin, J., Boek, E.S. and Coveney, P.V. (2002), “Lattice Boltzmann simulation of the flow of
binary immiscible fluids with different viscosities using the Shan-Chen microscopic
interaction model,” Philos. Trans. R. Soc. London, Ser. A, 360, 547–558.

184

Chudnovskii, A. (1966). “Plants and lght I. Radiant energy,” Fundamentals of agrophysics,
Israel Program of Scientific Translations, Jerusalem, Israel, 1-51.
Clennell, M.B. (1997), “Tortuosity: a guide through the maze,” In Lovetl, M.A., Harvey, P.K.
(eds). Developments in Petro-physics, Geological Society Special Publication, 122:299344.
Corey, A.T. (1977), “Mechanics of heterogeneous fluids in porous media,” Water Resources
Publication, Littleton, Colorado.
Daniel, D. E. (1984), “Predicting Hydraulic Conductivity of Clay Liners,” Journal of
Geotechnical Engineering ASCE, 110(2): 285-300.
Dardis, O. and McCloskey, J. (1998), “Lattice Boltzmann scheme with real numbered solid
density for the simulation of flow in porous media,” Physical Review E. 57(4):4834.
DeCoursey, D. G. and Rojas. K. W. (1990), “RZWQM- A model for simulating the movement of
water and solutes in the root zone,” In Proc. Int. Symp. on Water Quality Modeling of
Agricultural Nonpoint Sources, ed. D. G. DeCoursey, Part 2,813-821. USDA-Agri. Res.
Serv. ARS-81.
Dullien, F.A.L. (1992), “Porous Media: Fluid Transport and Pore Structure,” 2nd edition.
Academic Press, San Diego.
EPA. (1992). “Subtitle D Clarification.” 40 CFR 257 and 258. Federal Register.
EPA. (2007). “Municipal Solid Waste Generation, Recycling, and Disposal in the United States:
Facts and Figures for 2006.” EPA-530-F-07-030. November.
Eigenbrod, K.D. (2003), “Self-healing in fractured fine-grained soils,” Can. Geotech. J. 40: 435–
449.
Farahani, H. J., and Ahuja. L. R. (1996), “Evapotranspiration modeling of partial canopy/residue
covered field,” Trans. ASAE, 39(6): 2051-2064.
Fayer, M. (2000), UNSAT-H version 3.0: Unsaturated soil water and heat flow model – Theory,
user manual, and examples. PNNL-13249. Richland, Washington: Pacific Northwest
Laboratories
Fayer, M., Rockhold, M. and Campbell, M. (1992).” Hydraulic modeling of protective barriers:
Comparison of field data and simulation results,” Soil Science Society of America
Journal, 56: 690-700.
Ferreol, B. and Rothman, D.H (1995), Lattice-Boltzmann simulations of flow through
Fontainebleau sandstone, Transport Porous Media, 20: 3 –20.

185

Freed, D. M. (1998), Lattice-Boltzmann method for macroscopic porous media modeling.
International Journal of Modern Physics C, 9(08):1491-503.
Gerke, H.H. (2006), “Preferential flow descriptions for structured soils,” Journal of Plant
Nutrition and Soil Science, 169, 382–400
Glass, R.J., Steenhuis, T.S., and Parlange, J.Y., (1988). “Wetting front instability as a rapid and
far-reaching hydrologic process in the vadose zone,” Journal of Contaminant Hydrology,
3 (2–4), 207–226.
Green, W. H. and Ampt, G. A. (1911), “Studies on soil physics. 1: Flow of air and water through
soils,” J. Agr. Sci. 4:1-24.
Greve, A., Andersen, M.S. and Acworth, R.I.. (2010), “Investigations of soil cracking and
preferential flow in a weighing lysimeter filled with cracking clay soil,” Journal of
Hydrology, 393: 105–113.
Gunstensen, A. K., Rothman, D. H., Zaleski, S. and Zanetti, G. (1991), “Lattice Boltzmann
model of immiscible fluids,” Phys. Rev. A, 43(8), 4320– 4327.
Guyonnet, D., Amraoui, N. and Kara, R. 2000. Analysis of transient data from infiltrometer tests
in fine-grained soils. Ground Water 38 (3): 396-402.
Hardie, M. A., Cotching, W. E., Doyle, R. B., Holz G., Lisson, S., and Mattern, K. (2011),
“Effect of antecedent soil moisture on preferential flow in a texture-contrast soil,”
Journal of Hydrology, 398, 191–201
Haverkamp, R., Valcin, M., Touma, J., Wierenga, P., and Vauchaud, G. 1977. “A Comparison of
Numerical Simulation Models for One-dimensional Infiltration,” Soil Science Society of
America Journal, 41: 285-294.
Hauser, V. L., Weand B. L., and Gill, M. D. (2001), “Natural Covers for Landfills and Buried
Waste,” Journal of Environmental Engineering, 127(9), 768-775.
He X, and Luo L-S. (1997), “Theory of the lattice Boltzmann method: From the Boltzmann
equation to the lattice Boltzmann equation,” Physical Review E, 56(6):6811.
Hillel D. (1980), “Applications of soil physics,” Academic Press, Inc., New York
Hillel, D., Krentos, V. and Stylianou, Y. (1972). “Procedure and test of an internal drainage
method for measuring soil hydraulic characteristics in situ.” Soil Sci., Vol. 114, 5, 395400.
Holtz, R.D. and Kovacs, W.D. (2010), “An Introduction to Geotechnical Engineering,” Prentice
Hall, 2nd edition.

186

Huang, H. and Lu, X. (2009), “Relative permeabilities and coupling effects in steady-state gasliquid flow in porous media: A lattice Boltzmann study, “Physics of Fluids, 21, 092104.
Instruction Manual. (2003). “Toxic vapor analyzer (TVA-1000B).” Thermo electron corporation,
Environmental instruments, 27 Forge Parkway, Franklin, Massachusetts, 02038.
Jarvis, N.J. (2007), “A review of non-equilibrium water flow and solute transport in soil
macropores: principles, controlling factors and consequences for water quality,”
European Journal of Soil Science, 58 (3), 523–546.
Kang, Q, Zhang, D, and Chen, S (2002), “Unified lattice Boltzmann method for flow in
multiscale porous media,” Physical Review E. 66(5):056307.
Khire, M.V., Benson, C. and Bosscher, P. (1999), “Field data from a capillary barrier and model
predictions with UNSAT-H,” Journal of Geotechnical and Geoenvironmental
Engineering, 125: 518–528
Khire, M.V., Benson, C. and Bosscher, P. (1997). “Water balance modeling of earthen final
covers,” Journal of Geotechnical and Geoenvironmental Engineering, 123(8): 744-754.
Khire, M.V, Benson, C. and Bosscher, P. (2000), “Capillary barriers: Design variables and water
balance,” Journal of Geotechnical and Geoenvironmental Engineering, 126(8), 695-708.
Khire, M. V., Meerdink, J., Benson, C. H. and Bosscher, P. (1995). ‘‘Unsaturated hydraulic
conductivity and water balance predictions for earthen landfill final covers.’’ Soil Suction
Applications in Geotechnical Engineering Practice, W. Wray and S. Houston, eds.,
ASCE, New York, 38–57.
Khire, M.V. and Saravanathiiban, D.S. (2013), “Micropore vs. Macropore Flow: Implications for
Landfill Final Cover Design,” Coupled Phenomena in Environmental Geotechnics, 407.
Kim, W.H. and Daniel, D.E. (1992). “Effect of freezing on hydraulic conductivity of compacted
clay.” J. Geotech. & Geoenviron. Engrg., 118(7): 1083-1097.
Kohne, J.M., Kohne, S. and Simunek, J. (2009), “A review of model applications for structured
soils: a) Water flow and tracer transport,” Journal of Contaminant Hydrology, 4-35.
Körner, C., Thies, M. Hofmann, T.,Thürey, H., Rüde, U., (2005), “Lattice Boltzmann model for
free surface flow for modeling foaming,” J. Stat. Phys. 121(112): 179-196.
Kumar, A., Kanwar, R.S. and Ahuja, L.R. (1998), “Evaluation of preferential flow component of
RZWQM in simulating water and atrazine transport to subsurface drains,” Trans. ASAE,
41(3): 627-637.
Kung, K.-J.S. (1990), “Preferential flow in a sandy vadose zone. I. Field observation,”
Geoderma, 46: 51–58.
187

Kustas, W., Rango, A., and Uijlenhoet, R. (1994). “A simple energy budget algorithm for the
snow melt runoff model.” J. Geotech. Engr., 122(7): 565-576.
Kutay, M.E., Arambula, E., Gibson, N.H. and Youtcheff, J. (2010), “Three-Dimensional Image
Processing Methods to Identify and Characterize Aggregates in Compacted Asphalt
Mixtures,” International Journal of Pavement Engineering. 11(6): 511-528.
Kutay, M.E., Aydilek, A.H. and Masad, E. (2006), “Laboratory validation of lattice Boltzmann
method for modeling pore-scale flow in granular materials,” Computers and Geotechnics,
33(8): 381–395.
Leavesley, G. H., Lichty, R. W., Troutman, B. M. and Saindon, L. G. (1983), “Precipitationrunoff modeling system: User’s manual,” USGS Water Resources Investigations 834238, Denver, Colo., U.S. Geological Survey.
Ma, L., Ahuja, L. R., Nolan, B. T., Malone, R. W., Trout, T. J. and Qi, Z. ( 2012), “Root Zone
Water Quality Model (RZWQM2): Model Use, Calibration and Validation,” Trans.
ASABE, 55:1425-1446.
Ma, L., Ahuja, L. R., Saseendran, S. A., Malone, R. W., Green, T. R., Nolan, B. T., Bartling, P.
N. S., Flerchinger, G. N., Boote, K. J., and Hoogenboom, G. (2011), “A Protocol for
parameterization and calibration of RZWQM2 in field research,” In: Ahuja, L. R. and
Ma, L. (Eds). Methods of Introducing System Models into Agricultural Research, SSSA
book series, 1-64.
Malone, R. W., Shipitalo, M. J., Ma, L., Ahuja, L. R. and Rojas, K. W. (2001), “Macropore
component assessment of the Root Zone Water Quality Model (RZWQM) using no-till
soil blocks,” Trans. ASAE 44(4): 843-852.
Malone, R. W., Ma, L., Wauchope, R. D., Ahuja, L., Rojas, K., Ma, Q., Warner, R. and Byers,
M. (2004). “Modeling hydrology, metribuzin degradation, and metribuzin transport in
macroporous tilled and no-till silt loam soil using RZWQM,” Pest Mgmt. Sci. 60(3): 253266.
Martys, N. (2001), “Improved approximation of the Brinkman equation using a lattice
Boltzmann method,” Physics of fluids, 13(6), 1807-1810.
Martys, N. and Chen, H. (1996), “Simulation of multicomponent fluids in complex threedimensional geometries by the lattice Boltzmann method,” Phys. Rev. E, 53(1b), 743–
750.
Martys NS, Shan X, Chen H. (1998), “Evaluation of the external force term in the discrete
Boltzmann equation,” Physical Review E, 58(5):6855-6857.
MDE. (2003). “Estimated Costs of Landfill Closure Fact Sheet.” Maryland Department of the
Environment, Baltimore, MD.
188

Meerdink, J., Benson, C. and Khire, M.V. (1996), “Unsaturated hydraulic conductivity of two
compacted barrier soils,” Journal of Geotechnical Engineering, 122(7): 565–576.
Melchior, S. (1997). “In situ studies on the performance of landfill caps.” Proc., Int.
Containment Technology Conf., St. Petersburg, Fla., 365–373.
Mijares, R. G. (2011). “Hydraulic evaluation of lysimeters versus actual evapotranspirative
caps.” PhD. Dissertation, Michigan State University, East Lansing, MI.
Mijares, R. G. and Khire, M. (2010), “Soil water characteristic curves of compacted clay
subjected to multiple wetting and drying cycles,” ASCE Geotechnical Special Publication
No. 199: Advances in Analysis, Modeling and Design, Reston, VA., 400-409.
Mijares, R. and Khire, M.V. (2012), “Field Data and Numerical Modeling of Water Balance of
Lysimeter versus Actual Earthen Cap,” Journal of Geotechnical and Geoenvironmental
Engineering, 138(8): 889-897.
Mijares, R., Khire, M.V. and Johnson, T. (2012), “Field-scale evaluation of lysimeters versus
actual earthen covers,” Geotechnical Testing Journal, 35(1): 31–40.
Miller, C. J. and Mishra, M. (1989). “Modeling Leakage Through Cracked Clay Liners II: A new
Perspective,” Water Resources Bulletin, 25(3): 551-656.
Misici, L. and Palpacelli, S. (2005), “A lattice Boltzmann approach for immiscible fluids with
very different viscosities,” International Journal of Modern Physics C, 16(9), 1409-1435.
Montgomery, R. and Parsons, L. (1989). “The Omega Hills final cover test plot study: Three
year data summary.” Proc., 1989 Annual Meeting of the National Solid Waste
Management Association, Washington, D.C., 1–14.
Newmen, M.S., and Yin, X. (2013), “Lattice Boltzmann Simulation of Non-Darcy Flow in
Stochastically Generated 2D Porous Media Geometries,” SPE Journal, 18(1), 12–26.
Nimah, M., and Hanks, R. J. (1973), “Model for estimating soil water-plant-atmospheric
interrelation: I. Description and sensitivity,” SSSA Proc. 37(4): 522-527.
Nyhan, J., Schofield, T., and Starmer, R. (1997). “A Water Balance Study of Four Landfill
Cover Designs Varying in Slope for Semi-Arid Region.” Journal of Environmental
Quality, 26, 1385-1392.
Othman, M. and Benson, C. (1994). “Effect of freeze-thaw on the hydraulic conductivity and
morphology of compacted clay.” Can. Geotech. J., Vol. 30(2): 236–246.
Pan, C., Hilpert, M. and Miller, C. T. (2004a), “Lattice-Boltzmann simulation of two-phase flow
in porous media,” Water Resour. Res., 40, W01501, doi:10.1029/2003WR002120.
Panton, R.L. (2005), “Incompressible Flow,” Wiley. 3rd edition.
189

Phifer, M., Drumm, E. and Wilson, G. (1994). “Effects of post compaction water content
variation on saturated conductivity.” Hydraulic conductivity and waste contaminant
transport, STP 1142, D. Daniel, and S. Trautwein, eds., ASTM, Philadelphia, 318–334.
Qian Y, d'Humieres D, and Lallemand P. (1992), “Lattice BGK models for Navier-Stokes
equation,” EPL (Europhysics Letters), 17(6):479.
Richards LA. (1931), “Capillary conduction of liquids through porous mediums,” Physics 1:318333.
Robbins, G. A., Deyo, B. G., Temple, M. R., Stuart, J.D. and Lacy, M. J. (1990), “Soil-gas
surveying for subsurface gasoline contamination using total organic vapor detection
instruments: Part I. Theory and laboratory experimentation,” Ground Water Monit.
Remed., 10:122–131.
Saravanathiiban, D.S. and Khire, M.V. (2014a) “Controlled irrigation to estimate field-scale
hydraulic conductivity of a landfill final cover,” ASCE Geotechnical Special Publications
No. 234: Geo-Characterization and Modeling for Sustainability. Reston, VA
Saravanathiiban, D.S. and Khire, M.V. (2014b), “Field-scale unsaturated hydraulic properties of
compacted and uncompacted earthen covers,” ASCE Geotechnical Special Publications
No. 234: Geo-Characterization and Modeling for Sustainability. Reston, VA
Saravanathiiban, D.S., Kutay, M.E. and Khire, M.V. (2014c), “Effect of tortuosity and
macropore morphology on preferential flow through saturated soil: A lattice Boltzmann
study,” Computers and Geotechnics, 59: 44-53.
Scanlon, B., Reedy, R., Keese, K., and Dwyer, S. (2005). “Evaluation of Evapotranspiration
Covers for Waste Containment in Arid and Semiarid Regions in the Southwestern USA,”
Vadose Zone Journal, 4: 55-71.
Schroeder, P. R., Dozier, T.S., Zappi, P. A., McEnroe, B. M., Sjostrom, J.W., and Peyton, R. L.
(1994). "The Hydrologic Evaluation of Landfill Performance (HELP) Model:
Engineering Documentation for Version 3," EPA/600/R-94/168b, September 1994, U.S.
Environmental Protection Agency Office of Research and Development, Washington,
DC.
Shan, X. and Chen, H. (1993), “Lattice Boltzmann model for simulating flows with multiple
phases and components,” Phys. Rev. E, 47, 1815–1819.
Shan, X. and Chen, H. (1994), “Simulation of non-ideal gases and liquid-gas phase transitions by
the lattice Boltzmann equation,” Phys. Rev. E, 49(4), 2941– 2948.
Simunek, J., Sejna, M., Saito, H., Sakai, M. and van Genuchten, M.Th. (2008), “The HYDRUS1D software package for simulating the one-dimensional movement of water, heat, and
multiple solutes in variably saturated media,” Version 4.0, HYDRUS Software Series 4.
University of California Riverside, Riverside, CA, USA.
190

Simunek, J., Sejna, M. and van Genuchten, M.Th. (1998), “The HY-DRUS-1D software package
for simulating the one-dimensional movement of water, heat, and multiple solutes in
variably saturated media,” Version 2.0. IGWMC-TPS-70, School of Mines, Golden,
Colorado: International Ground Water Modeling Center
Simunek, J. and van Genuchten, M. Th. (2008), “Modeling nonequilibrium flow and transport
with HYDRUS,” Vadose Zone Journal, 7 (2): 782–797.
Song, Q. and Yanful, E. K. (2010), “Laboratory and numerical modeling of water balance in a
layered sloped soil cover with channel flow pathway over mine waste rock,” Environ.
Earth Sci., doi:10.1007/s12665-010-0488-4.
Spaid, M.A.A. and Phelan, F. R. (1997) ‘‘Lattice Boltzmann methods for modeling microscale
flow in fibrous porous media,’’ Phys. Fluids, 9, 2468.
Succi, S. (2001), “The Lattice Boltzmann Equation for Fluid Dynamics and Beyond,” Oxford:
Clarendon Press.
Succi S., Foti, E. and Higuera, F. (1989), “Three dimensional flows in complexgeometries with
the lattice Boltzmann method,” Europhysics Letter, 10: 433–438.
Sukop, M. C., Huang, H., Alvarez, P. F., Variano, E. A. and Cunningham, K. J. (2013),
“Evaluation of permeability and non-Darcy flow in vuggy macroporous limestone aquifer
sam-ples with lattice Boltzmann methods,” Water Resources Research, 49: 216–230.
Sukop, M.C., and Thorne, D.T. (2006), “Lattice Boltzmann Modeling: An Introduction for
Geoscientists and Engineers,” Springer, Heidelberg.
Sukop, M. C., and Or, D. (2004), “Lattice Boltzmann method for modelingliquid-vapor interface
configurations in porous media,” Water Resour. Res., 40, W01509, doi:
10.1029/2003WR002333.
Suter, G., Luxmoore, R. and Smith, E. (1993). “Compacted soil barriers at abandoned landfill
sites are likely to fail in the long term.” J. Environ. Qual., Vol. 22(2), 217–226.
Swift, M. R., Orlandini, E., Osborn, W. R. and Yeomans, J. M. (1996), “Lattice Boltzmann
simulations of liquid-gas and binary fluid systems, Phys,” Rev. E, 54(5), 5041–5052.
Taylor, G., Spain, A., Nefiodovas, A., Timms, G., Kuznetsov, V., and Benntt, J. (2003),
“Determination of the Reasons for Deterioration of the Rum Jungle Waste Rock Cover,”
Australian Center for Mining Environmental Research, Brisbane.
Trast, J., and Benson, C.H. (1995). “Estimating field hydraulic conductivity at various effective
stresses.” J. Geotech. Engrg., Vol. 121(10), 736–740.

191

Udawatta, R. P., Anderson, A. H., Motavalli, P. P. and Garett, H. E. (2011), “Calibration of a
water content reflectometer and soil water dynamics for an agroforestry practice,”
Agroforest Syst. 82: 61-75.
van Genuchten, M.Th. (1980). “A Closed-form Equation for Predicting the Hydraulic
Conductivity of Unsaturated Soils.” Soil Science Society of America Journal, 44: 892898.
Walsh, S. D, Burwinkle, H, and Saar, M. O (2009), “A new partial-bounceback latticeBoltzmann method for fluid flow through heterogeneous media,” Computers and
Geosciences, 35(6):1186-93.
Wang, J.J., Zhang, H.P. Zhang, L. and Liang, Y. (2013), “Experimental study on self-healing of
crack in clay seepage barrier,” Engineering Geology, 159: 31–35
Ward, A.L. and Gee, G.W. (1997). “Performance Evaluation of a Field-Scale Surface Barrier,”
Journal of Environmental Quality, 26, 694-705.
Watson, K.K. (1966). “An instantaneous profile method for determining the hydraulic
conductivity of unsaturated porous materials.” Water Resources Research, 2(4): 709–715
Wolf-Gladrow, D.A. (2000), “Lattice-Gas Automata and Lattice Boltzmann Models: An
Introduction.” Berlin: Springer.
Workman, S.R., and Skaggs, R.W. (1990), “PREFLO: a water man-agement model capable of
simulating preferential flow,” Transactions of the ASAE, 33(6): 1939–1948.
Yuan, P. and Schaefer, L. (2006), “ Equations of state in a lattice Boltzmann model,” Physics of
Fluids, 18.

192