MULTIMEDIA EXPOSURE ASSESSMENT
MODEL (MULTIMED) FOR EVALUATING
THE LAND DISPOSAL OF WASTES--
MODEL THEORY
by
Atul M. Salhotra1
Phil Mineart1
Susan Sharp-Hansen2
Terry Allison3
Woodward-Clyde Consultants1
Oakland, CA 94607
AQUA TERRA Consultants2
Mountain View, CA 94043
Computer Sciences Corporation3
Athens, GA 30605-2720
EPA Contract No. 68-03-3513
and No. 68-03-6304
Project Monitor
Gerard Laniak
Environmental Research Laboratory
U.S. Environmental Protection Agency
Athens, GA 30605-2720
ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ATHENS, GEORGIA 30605-2720
-------
DISCLAIMER
The work presented in this document has been funded by the United States
Environmental Protection Agency. It has been subject to the Agency's peer and
administrative review, and has been approved as an EPA document. Mention of
trade names or commercial products does not constitute endorsement or
recommendation for use by the U.S. Environmental Protection Agency.
-------
FOREWORD
As environmental controls become more costly to implement and the penalties of
judgment errors become more severe, environmental quality management requires
more efficient management tools based on greater knowledge of the
environmental phenomena to be managed. As part of this Laboratory's research
on the occurrence, movement, transformation, impact and control of
environmental contaminants, the Assessment Branch develops management or
engineering tools to help pollution control officials assess the risk to human
health and the environment posed by land disposal of hazardous wastes.
This report describes a computer model for simulating the transport and
transformation of contaminants released from a hazardous waste disposal
facility into the multimedia environment. The MULTIMED model simulates
release to air and soil, including the unsaturated and saturated zones, and
possible interception of the subsurface contaminant plume by a surface stream.
It further simulates movement through the air, soil, groundwater and surface
water media to humans and other potentially affected species. MULTIMED is
intended for general exposure and risk assessments of waste facilities and for
analyses of the impacts of engineering and management controls.
Rosemarie C. Russo, Ph.D.
Director
Environmental Research Laboratory
Athens, GA
-------
ABSTRACT
The Environmental Protection Agency's Multimedia Exposure Assessment Model
(MULTIMED) simulates the movement of contaminants leaching from a waste
disposal facility. The model includes two options for simulating leachate
flux. Either the infiltration rate to the unsaturated or saturated zone can
be specified directly or a landfill module can be used to estimate the
infiltration rate. The landfill module is one-dimensional and steady-state,
and simulates the effect of precipitation, runoff, infiltration,
evapotranspiration, barrier layers (which can include flexible membrane
liners), and lateral drainage. A steady-state, one-dimensional, semi-
analytical module simulates flow in the unsaturated zone. The output from
this module, water saturation as a function of depth, is used as input to the
unsaturated zone transport module. The latter simulates transient, one-
dimensional (vertical) transport in the unsaturated zone and includes the
effects of longitudinal dispersion, linear adsorption, and first-order decay.
Output from the unsaturated zone modules--i.e., steady-state or time series
contaminant concentrations at the water table--is used to couple the
unsaturated zone transport module with the steady-state or transient, semi-
analytical saturated zone transport module. The latter includes one-
dimensional uniform flow, three-dimensional dispersion, linear adsorption,
first-order decay, and dilution due to direct infiltration into the
groundwater plume. Contamination of a surface stream due to the complete
interception of a steady-state saturated zone plume is simulated by the
surface water module. Finally, the air emissions and the atmosphere
dispersion modules simulate the movement of chemicals into the atmosphere.
The fate of contaminants in the various media depends on the chemical properties
of the contaminants as well as a number of media- and environment-specific
parameters. The uncertainty in these parameters is quantified using the Monte
Carlo simulation technique.
To enhance the user-friendly nature of the model, separate interactive pre- and
postprocessing software have been developed for use in creating and editing input
and in plotting model output.
This report provides the conceptual and theoretical details of the various
modules and the Monte Carlo simulation technique. A user's manual for
implementing MULTIMED is currently in progress. In addition, an application
manual describes the use of MULTIMED in modeling Subtitle D land disposal
facilities.
This report was submitted in partial fulfillment of Work Assignment Number 32,
Contract Number 68-03-3513 by Aqua Terra Consultants, under the sponsorship of
the U.S. Environmental Protection Agency. This report covers the period March
1990 to July 1990, and work was completed as of August 1990.
-------
TABLE OF CONTENTS
Page
Disclaimer ii
Foreword iii
Abstract iv
Figures viii
Tables x
Acknowledgements xi
1 . 0 INTRODUCTION 1
1.1 Overview of the Model 1
1.2 The Physical Scenario 1
1.3 Model Capabilities 2
1.4 Report Organization 3
2 . 0 THE LANDFILL MODULE 4
2.1 Introduction 4
2.2 Governing Equations and Solution Techniques 6
2.2.1 Surface Water Balance 6
2.2.1.1 Calculation of Runoff 7
2.2.1.2 Calculation of Evapotranspiration 10
2.2.2 Percolation through the Landfill 12
2.2.3 Landfill Control Features 12
2.2.3.1 Low-permeability Liners 13
2.2.3.2 The Drainage System 14
2.2.4 Solution Method for Identifying Steady-State
Leaching Rate 15
2 . 3 Assumptions and Limitations 17
2 . 4 Data Requirements 17
3 . 0 THE UNSATURATED ZONE FLOW MODULE 19
3.1 Introduction 19
3.2 Governing Equations and Solution Techniques 19
3.2.1 Spatial Discretization in the Unsaturated Flow
Module 23
3.3 Assumptions and Limitations 24
3 . 4 Data Requirements 24
4 . 0 UNSATURATED ZONE TRANSPORT MODULE 26
4.1 Introduction 26
4 . 2 Governing Equations 26
4.2.1 Unsteady-State Transport 26
4.2.2 Steady-State Transport 29
4.3 Assumptions and Limitations 31
4 . 4 Data Requirements 32
4.4.1 The Chemical Transformation Rate 32
4.4.2 The Distribution Coefficient 32
4.4.3 The Longitudinal Dispersion Coefficient 34
5 . 0 THE SATURATED ZONE TRANSPORT MODULE 36
5.1 Introduction 36
5 . 2 Governing Equations 36
5.2.1 Solution to the Gaussian Source Boundary
-------
Condition 41
5.2.2 Solution to the Patch Source Boundary Condition... 43
5.2.3 Receptor Well or Stream Location 46
5.2.4 Time Values for Saturated Zone Concentrations 48
5 . 3 Assumptions and Limitations 48
5.4 Coupling of the Unsaturated and the Saturated Zone
Modules 49
5.4.1 Steady-State Coupling 49
5.4.1.1 The Gaussian Source Boundary Condition.. 50
5.4.1.2 The Patch Source Boundary Condition 51
5.4.2 Unsteady-State Coupling 52
5 . 5 Data Requirements 52
5.5.1 Source-Specific Parameters 53
5.5.1.1 The Gaussian Source Boundary
Condition 53
5.5.1.2 The Patch Source Boundary Condition 57
5.5.2 Chemical-Specific Parameters 57
5.5.2.1 The Overall Decay Coefficient 57
5.5.2.2 The Distribution Coefficient 59
5.5.3 Aquifer-Specific Parameters 59
5.5.3.1 Retardation Coefficient 59
5.5.3.2 Porosity and Mean Particle Diameter 59
5.5.3.3 Bulk Density 59
5.5.3.4 Seepage Velocity 60
5.5.3.5 Hydraulic Conductivity 60
5.5.3.6 Dispersion Coefficients 61
5.5.3.7 Source Thickness (Mixing Zone Depth).... 62
6 . 0 THE SURFACE WATER MODULE 63
6.1 Introduction 63
6.2 Mathematical Description of the Surface Water Model 63
6.2.1 Computation of the In-Stream Decay Coefficient.... 66
6.2.2 Calculation of the Mean In-Stream Velocity 69
6.3 Exposure Due to Surface Water Contamination 69
6.3.1 Routes of Exposure 69
6.3.2 Human Exposure to Toxics through Drinking Water... 69
6.3.3 Human Exposure to Toxics Due to Fish Consumption.. 72
6.3.4 Exposure of Aquatic Organisms to Toxics 74
6 . 4 Assumptions and Limitations 74
6 . 5 Data Requirements 74
7 . 0 THE AIR EMISSIONS MODULE 77
7.1 Introduction 77
7 . 2 Governing Equations 77
7.2.1 The Air Emissions Diffusion Model 77
7.2.2 The Effect of Atmospheric Pressure Fluctuations... 78
7.2.3 Other Transport Processes 79
7.2.4 Computation of the Diffusion Coefficient 79
7.2.4.1 Effect of Engineering Controls 80
7.2.5 Computation of Vapor Concentration, Csl 80
7 . 3 Assumptions and Limitations 82
7 . 4 Data Requirements 82
8 . 0 AIR DISPERSION MODULE 84
8.1 Introduction 84
8 . 2 Governing Equations 84
8.2.1 The Gaussian Dispersion Equation 84
8.2.2 Area Source Approximations 89
8.2.3 Plume Rise 89
8.2.4 Estimation of Vertical Dispersion Coefficient 92
-------
8.2.5 Terrain Effects 93
. 3 Assumptions and Limitations 94
. 4 Data Requirements 94
-------
9 . 0 UNCERTAINTY ANALYSIS 96
9.1 Introduction 96
9.2 Statement of the Problem and Technical Approach 96
9.3 The Monte Carlo Analysis Technique 98
9 . 4 Uncertainty in the Input Variables 100
9.5 Description of the Parameter Distributions 101
9.5.1 Normal Distribution 102
9.5.2 Log-Normal Distribution 102
9.5.3 Uniform Distribution 103
9.5.4 Log-Uniform Distribution 103
9.5.5 Exponential Distribution 103
9.5.6 Empirical Distribution 104
9.5.7 The Johnson System of Distributions 104
9.6 The Random Number Generator 104
9 . 7 Analysis of the Model Output 105
9.8 Confidence Bounds for the Estimated Percentile 112
9. 9 The Number of Monte Carlo Simulation Runs 115
10.0 REFERENCES 116
11.0 APPENDIX A - Simplified Estimation for Mixing Zone Depth 124
-------
FIGURES
Page
2.1 Profile of a typical hazardous waste landfill 5
2.2 Depiction of a typical event duration 8
2.3 Soil moisture available for evapotranspiration 11
3.1 A schematic of the waste facility and leachate migration
through the unsaturated and saturated zones 20
4.1 A schematic of transport through the layered
unsaturated zone 30
5.1(a) A schematic diagram of the Gaussian source boundary
condition for the saturated zone transport module 38
5.1(b) A schematic diagram of the Patch source boundary
condition for the saturated zone transport module 40
5.2 A schematic of the source penetration and the well location... 47
6.1(a) Groundwater contaminant plume interception
by the surface stream 64
6.1(b) Downstream contaminant transport from the edge
of the initial mixing zone 64
6.2(a) Stages between failure of the waste containment facility
and human exposure due to drinking water 70
6.2(b) Stages between failure of the waste containment facility
and human exposure due to fish consumption 71
6.2(c) Stages between failure of the waste containment facility
and exposure to Aquatic Organisms 71
8.1 Wind direction sectors for Gaussian models 87
8.2 Virtual source approximation for long-term Gaussian models.... 90
8.3 Effective source area used by VALLEY 91
9.1 A schematic description of the Monte Carlo method
of uncertainty analysis 99
9.2 Selecting a Johnson distribution from skewness and kurtosis... 106
9.3 Comparison of the exact and the generated cumulative
frequency distribution for a normally distributed variable.... 109
9.4 Comparison of the exact and the generated cumulative
frequency distribution for a log-normally distributed
variable 109
9.5 Comparison of the exact and the generated cumulative
-------
frequency distribution for an exponentially distributed
variable 110
9.6 Comparison of the exact and the generated cumulative
frequency distribution for an empirically distributed
variable 110
9.7 Comparison of the exact and the generated cumulative
frequency distribution for a uniformly distributed variable... Ill
9.8 Typical results obtained using MULTIMED in the Monte
Carlo mode 113
-------
TABLES
2-1
3-1
4-1
4-2
5-1
5-2
Parameters Required for the Source Module
Parameters Required for the Unsaturated z
Parameters Required for Unsaturated Zone
Compilation of Field Dispersivity Values
Parameters Required for the Saturated Zon
one Flow Module
Transport Module
e Transport Module...
Additional Data Required to Compute Parameters for
the Saturated Zone Transport Module
Page
18
25
33
35
54
55
5-3 (a) Alternatives for Including Dispersivities in the
Saturated Zone Module ......................................... 62
5-3 (b) Probabilistic Representation of Longitudinal Dispersivity
6
7
7
7
8
8
8
9
9
-1
-1
-2
-3
-1
-2
-3
-1
-2
for Distance of 152.4 m
Parameters Required for the Surface Water Module
General Characteristics of Atmospheric Pressure Fluctuations..
State of Chemical Substances in Landfills and
Vapor Density Equilibrium Laws
Parameters Required for the Air Emissions Module
Atmospheric Stability Categories
Constants for Pasquill-Gif ford Curves for Each Stability Class
Parameters Required for the Air Dispersion Module
Qualitative Comparison of Uncertainty-Propagation Methods
(a) Results of Random Number Generator Test for 500 Values....
62
75
78
81
83
86
92
95
98
107
9-2 (b) Results of Random Number Generator Test for 1000 Values ..... 108
-------
ACKNOWLEDGEMENTS
A number of individuals were involved in the actual development and
implementation of the MULTIMED code. Key individuals include Jan Kool and
Peter Huyakorn of HydroGeoLogic Inc., Barry Lester of Geotrans Inc., Michael
Ungs of TetraTech, Inc., Bob Ambrose of U.S. EPA, Jack Kittle of AQUA TERRA
Consultants, Rob Schanz, Yvonne Meeks, and Peter Mangarella of Woodward-Clyde
Consultants, and Terry Allison of Computer Sciences Corporation.
The pre- and postprocessor for MULTIMED were developed by Jack Kittle and Paul
Hummel at AQUA TERRA Consultants. Paul Hummel and Constance Travers tested
and created tutorials for the pre- and postprocessors.
The draft documentation of the MULTIMED model theory was written by Atul M.
Salhotra and Phil Mineart at Woodward-Clyde Consultants in 1988 under EPA
Contract No. 68-03-6304. It was corrected and updated by Susan Sharp-Hansen
at AQUA TERRA Consultants in 1990 under EPA Contract No. 68-03-3513. At AQUA
TERRA, the revised document was reviewed by Constance Travers and Tony
Donigian. Word processing was performed by Dorothy Inahara and the figures
were prepared by Melissa Donigian. Section 9.5, the description of Monte
Carlo distributions, is taken, with slight modifications, from Volume 1 of the
RUSTIC documentation (Dean et al., 1989).
We thank the Project Monitors at the U.S. EPA Athens Laboratory, Lee Mulkey
and Gerard Laniak, for their continuous technical and management support
throughout the course of this project.
-------
SECTION 1
Introduction
1.1 Overview of the Model
This chapter provides an overview of the U.S. Environmental Protection
Agency's Multimedia Exposure Assessment Model (MULTIMED). The model simulates
the fate and transport of contaminants released from a waste disposal facility
into the multimedia environment. Release to either air or soil, including the
unsaturated and the saturated zone, and possible interception of the
subsurface contaminant plume by a surface stream are included in the model.
Thus, the model can be used as a technical and quantitative management tool to
address the problem of the land disposal of chemicals. At this time, the air
modules of the model are not linked to the other model modules. As a result,
the estimated release of contaminants to the air is independent of the
estimated contaminant release to the subsurface and surface water.
MULTIMED utilizes analytical and semi-analytical solution techniques to solve
the mathematical equations describing flow and transport. The simplifying
assumptions required to obtain the analytical solutions limit the complexity
of the systems which can be represented by MULTIMED. The model does not
account for site-specific spatial variability, the shape of the land disposal
facility, site-specific boundary conditions, multiple aquifers, or pumping
wells. Nor can MULTIMED simulate processes, such as flow in fractures and
chemical reactions between contaminants, which can have a significant affect
on the concentration of contaminants at a site. In more complex systems, it
may be beneficial to use MULTIMED as a "screening level" model which would
allow a user to obtain an understanding of the system. A numerical model
could then be used if there are sufficient data and necessity to justify the
use of a more complex model.
1.2 The Physical Scenario
The physical scenario simulated by the model is a land disposal facility that
releases pollutants into the air, soil, and/or groundwater. In response to a
number of complex physical, chemical, and biological fate and transport
processes, the pollutants move in the multimedia environment (air, water, and
soil), resulting in potential toxic exposure to humans and other receptors.
Note that the processes affecting air emissions are not linked in the model to
the processes affecting subsurface transport. In other words, the
concentration calculated in the one medium is not affected by the release of
the contaminant to the other medium.
The sources of pollutants considered in the current version are either the
leachate from a waste disposal facility or air emissions during the post-
closure period. Inadequate long-term functioning or failure of the facility's
engineering controls (i.e., the caps and liners) are assumed to occur after
closure and result in the release of leachate to soil or groundwater beneath
the facility and emission of vapor to the atmosphere. Note that the use of
the air-emissions module is most appropriate for high concentrations of waste
in the facility. Also, the model does not include fate processes, such as
complexation and solids precipitation, that affect the transport of metals.
-------
-------
1.3 MODEL CAPABILITIES
During the development of this model, emphasis has been placed on the creation
of a unified, user-friendly, software framework, with the capability to
perform uncertainty analysis, that can be easily enhanced by adding modules
and/or modifying existing modules.
The major functions currently performed by this model include:
Allocation of default values to some input parameters/variables.
Reading of the input data files.
Echo of input data to output files.
Derivation of some parameters, if specified by the user.
Depending on user-selected options:
> simulation of leachate flux emanating from the source
* simulation of unsaturated zone flow and transport
> simulation of saturated zone transport only
* computation of in-stream concentrations due to contaminant loading
assuming complete interception of a plume in the saturated zone
> computation of the rate of contaminant emission from the waste
disposal unit into the atmosphere
* simulation of dispersion of the contaminants into the atmosphere
Generation of random numbers for Monte Carlo simulations.
Performance of statistical analyses of Monte Carlo simulations.
Writing the concentrations at specified receptors to output for
deterministic runs. In the Monte Carlo mode, writing the cumulative
frequency distribution and selected percentiles of concentrations at
receptors to output.
Printing the values of randomly generated input parameters and the
computed concentration values for each Monte Carlo run.
The fate and transport of contaminants critically depends on a number of
media-specific parameters. Typically many of these parameters exhibit spatial
and temporal variability as well as variability due to measurement errors.
MULTIMED has the capability to analyze the impact of uncertainty and
variability in the model inputs on the model outputs (concentrations at
specified points in the multimedia environment), using the Monte Carlo
simulation technique.
To enhance the user-friendly nature of the model, separate interactive
preprocessing and postprocessing software has been developed, using the ANNIE
Interaction Development Environment (AIDE) (Kittle et al. , 1989), for use in
creating and editing input and in plotting model output. The pre- and
postprocessors have not been integrated with MULTIMED because of the size
limitations of desktop computers. Therefore, after using the preprocessor to
create or modify input, the model is run in batch mode. Afterwards, the
postprocessor can be used to produce plots of the Monte Carlo output or
concentration versus time.
Finally, model results can be used to manually 'back-calculate' the maximum
source concentration (for a steady-state, infinite contaminant source) of a
chemical which ensures protection of human health and/or the environment at a
-------
down-gradient point of exposure (see Section 9.6). Details of the back-
calculation procedure for an unsteady-state (finite contaminant source) case
have been discussed by Mulkey and Allison (1988) .
1.4 Report Generation
This report includes the theoretical background and other information
necessary to understand MULTIMED. Chapter 2 describes the landfill source
module, which can be used to estimate leachate rates from a waste disposal
facility. Chapter 3 contains a description of the steady-state, one-
dimensional, unsaturated zone flow module, which computes saturation as a
function of depth. This information is used by the landfill module and the
unsaturated zone transport module, which is described in Chapter 4. The
saturated zone transport module and its coupling with the unsaturated zone is
described in Chapter 5. Chapter 6 discusses the surface water module, which
computes in-stream concentrations of a contaminant based on the assumption of
complete interception of a steady-state saturated zone plume. The emission of
contaminants from the facility into the atmosphere and their transport and
dispersion are described in Chapters 7 and 8, respectively. Chapter 9
contains a discussion of the Monte Carlo simulation technique that is
incorporated into the model. Each of these chapters presents the necessary
mathematical equations and relationships, limitations, and major assumptions
of the modules and details of the required input parameters.
A user's manual for implementing MULTIMED is in progress at this time. In
addition, a MULTIMED Subtitle D application manual has been written (Sharp-
Hansen et al., 1990). The manual explains how to use the unsaturated and
saturated zone modules of MULTIMED to study and design Subtitle D land
disposal facilities. Installation of the code, the use of the pre- and
postprocessors, the format of the input and output, parameter estimation, and
example problems are addressed in the Subtitle D application manual.
-------
SECTION 2
The Landfill Module
2.1 Introduction
This section presents the details of the landfill module included in MULTIMED.
The landfill module provides a physically-based estimate of the amount of
leachate that may be generated at a RCRA Subtitle C or D landfill. The module
is one-dimensional and steady-state. It uses a water balance approach to
simulate the effects of hydrologic processes including precipitation,
evapotranspiration, runoff, infiltration, and lateral drainage. A schematic
diagram of a typical hazardous waste landfill is shown in Figure 2.1. The
design components shown in this figure have been used to develop the
conceptual basis of the landfill module. The user can adapt these components
to his specific problem by assigning varying layer thicknesses and hydraulic
properties, by rearranging layers, and by omitting or adding layers. Thus, a
variety of climatic conditions and typical landfill designs may be specified
by the user. Note that the landfill module was initially developed as a
"stand-alone" code and was later integrated into MULTIMED. Details of the
stand alone code are presented in Meeks et al. (1988).
In general, two options are available in the code for the simulation of
leachate flux from the source. Under option 1, the landfill module is not
used and the user specifies an infiltration rate and chemical concentration.
The unsaturated zone flow module (described in Section 3) then simulates the
movement of the leachate to the water table. Under option 2, the infiltration
rate is not specified. Instead, the leachate flux is calculated by the
landfill module described in this section. Because the leachate flux depends
strongly on underlying soil moisture conditions, the landfill module must be
coupled with the unsaturated zone flow module. These two modules are solved
iteratively in a sequential manner until the percolation rates estimated by
the two modules agree. The landfill module does not estimate the chemical
concentration in the leachate. Thus, a chemical concentration must be
specified by the user for either option.
This chapter is organized into four sections. Section 2.2 describes the
governing equations and the iterative solution technique used to estimate a
leachate flux that satisfies both the surface water balance and the
unsaturated zone flow module. The major assumptions and limitations inherent
in the module are identified in Section 2.3
The data required in the landfill module are described in Section 2.4.
2.2 GOVERNING EQUATIONS AND SOLUTION TECHNIQUES
The volume of leachate generated from a landfill is governed by:
The availability of water
Underlying soil and refuse conditions
The amount of water available for leachate generation from a landfill is
controlled by climatic conditions, surface conditions (e.g., soil and
vegetation type, soil moisture conditions), groundwater influx, and liquids
-------
inherently associated with the waste. In this module the contributions to
leachate from groundwater influx and liquids inherent in the waste are
considered negligible. The amount of water available from other sources is
estimated using an event-based surface water balance technique described in
Section 2.2.1.
The underlying soil and refuse conditions include 1) the hydraulic properties
of subsurface materials, and 2) the design and functioning of engineering
controls. These conditions affect the amount of infiltration that can
percolate through the landfill layers and the unsaturated zone. The
unsaturated zone flow module described in Section 3 calculates percolation
based on a semi-analytical solution to the Richard's Equation. The coupling
of the unsaturated zone flow module with the landfill module for the leachate
calculation is discussed in Section 2.2.2. Effects of engineering controls
(i.e., synthetic liners and lateral drains) are described in Section 2.2.3.
Under steady-state conditions, the water available for infiltration must equal
the lateral drainage flux plus the percolation rate through the unsaturated
zone below the drainage layer. That is, the water made available for leaching
must pass through each of the landfill layers and out of the lateral drains or
through the bottom of the landfill. This is necessitated by the fact that no
change in storage of water within the landfill can occur in the steady-state
formulation of the problem. Section 2.2.4 describes how the landfill module
uses the prescribed equality of these fluxes to calculate the steady-state
leachate rate.
The landfill module solves for the steady-state leachate rate based on
representative steady-state climatic data, which must be estimated from actual
data for precipitation, storm interval, etc. Because climatic variables, such
as evapotranspiration rates and precipitation characteristics, may have strong
seasonal variability, the code includes the option to consider seasonal
variability in precipitation events and evapotranspiration while retaining the
assumption of steady-state (see Section 2.2.1.2). This option should not be
used to estimate transient leaching rates because information is not passed
from one season to the following season.
2.2.1 Surface Water Balance
A surface water balance is used to estimate the water available for
infiltration into the landfill. The water balance is based on an event
approach. That is, the water balance is evaluated for the period of time from
the beginning of a precipitation event, and through the period of no
precipitation which follows. The calculation period ends at the onset of the
next precipitation event. Figure 2.2 depicts the event approach. The
algorithm used in the landfill module partitions precipitation into runoff,
evapotranspiration, and infiltration. Changes in surface water storage (e.g.,
snow melt) are not included. Because the module is steady-state, changes in
soil moisture are precluded and do not need to be considered in the water
balance. Infiltration (volume per unit area) into the soil during a
precipitation event can then be expressed by (Dass et al. , 1977):
INFIL = PRECIP - RO - ET (2.1)
where
INFIL = water available for infiltration during an event [cm]
PRECIP = storm precipitation [cm]
RO = runoff during an event [cm]
ET = evapotranspiration during the event interval [cm]
-------
-------
The procedures used to estimate the volume of runoff and of actual
evapotranspiration are described individually in the following sections.
Based on the calculated infiltration volume, the average steady-state
infiltration rate is given by:
I = INFIL/DEVENT (2.2)
where
I = infiltration rate [cm/d]
DEVENT = event duration [d]
2.2.1.1 Calculation of Runoff
Some fraction of the incident precipitation runs off the landfill and is lost
to overland flow before it has a chance to infiltrate. The landfill module
computes runoff by the SCS runoff curve number method in the manner described
in the HELP model documentation (Schroeder et al., 1984b). The SCS procedures
were developed from observed runoff-rainfall relationships for large storms on
small watersheds.
The relation between precipitation, runoff, and retention for a particular set
of environmental conditions was found to be:
_ (PRECIP- 0.2S)2
(PRECIP + 0.8S) '
where
RO = runoff during an event [cm]
PRECIP = storm precipitation [cm]
S = retention parameter [cm]
The retention parameter, S, for a given soil varies as a function of the soil
moisture in the underlying soil (Schroeder et al., 1984b):
., r. SM - WP,
s = s^--] (2-4)
where
Smx = maximum value of the retention parameter [cm]
SM = soil water content in the upper soil layer [dimensionless
fraction]
UL = upper limit moisture content in the upper soil layer
[dimensionless fraction]
WP = wilting point moisture content in the upper soil layer
[dimensionless fraction]
-------
The upper limit of soil moisture content is the moisture content at saturation
and is numerically equal to the effective porosity. The wilting point
moisture content is the lowest naturally occurring soil water content and is
equal to the effective porosity multiplied by the residual saturation. Since
soil water is not distributed uniformly throughout the soil profile and since
the soil moisture near the surface influences infiltration more strongly than
that located elsewhere, the retention parameter is depth-weighted in the
landfill module (Meeks et al., 1988). The soil profile in the uppermost soil
layer is divided into seven segments. The thickness of the top segment is set
equal to one thirty-sixth of the total upper layer thickness and the thickness
of the second segment is five thirty-sixths of the layer's thickness. The
thickness of each of the bottom five segments in the uppermost soil layer is
defined as one-sixth of the total layer thickness. The depth-weighted
retention parameter is computed with the following equation (Knisel, 1980):
SM. - WP
"1 (^W}] (2'5)
where
SMj = soil moisture content of segment j [fraction]
W.j = weighting factor for segment j [dimensionless]
The weighting factors decrease with the depth of the segment. In accordance
with the development of CREAMS (Knisel, 1980), the weighting factors for
segments 1, 2, 3, 4, 5, 6 and 7 are 0.111, 0.397, 0.254, 0.127, 0.063, 0.032,
and 0.016, respectively.
The maximum value of the retention parameter, Smx, is the value of S at the
lowest soil moisture content and is calculated from the SCS curve number for
average moisture conditions (CNIZ) .
2.2.1.2 Calculation of Evapotranspiration
Evapotranspiration from a landfill is a function of climatic conditions,
vegetation, soil moisture, and the ability of the soil to transmit water and
water vapor. A two-step approach is taken to calculate actual
evapotranspiration.
First, the potential evapotranspiration during the interval between storms is
calculated from the seasonal potential evapotranspiration input by the user.
This calculation takes the form:
where
EPET = event potential evapotranspiration [cm]
PET = seasonal potential evapotranspiration [cm]
DEVENT = event duration [d]
ISEA = number of seasons in a year [dimensionless]
-------
The second step involves estimating the availability of water, stored as soil
moisture in the shallow subsurface, to meet the evapotranspiration demand.
The available soil moisture is estimated as a function of the moisture content
above the wilting point in the upper soil layer. The user must take care in
defining the thickness of this layer so that its depth accurately reflects the
depth over which evapotranspiration takes place. In general, the zone of
evapotranspiration extends only a few feet or less from the surface. If the
top of the landfill is vegetated, it is suggested that the uppermost layer
thickness correspond to the root zone depth. The module assumes that at the
top of the uppermost layer the soil moisture can be depleted to the wilting
point and that at the bottom of the uppermost layer the soil moisture is not
affected by evapotranspiration. The soil moisture level above which soil
moisture is available is linearly interpolated between these two depths as
shown in Figure 2.3. That is, a triangular distribution is assumed from the
surface to the bottom of the uppermost layer with the maximum soil moisture
taken from near the surface. This approach is similar to that taken in the
PRZM model (Carsel et al., 1984). The available moisture for
evapotranspiration can be expressed as:
AW = 100 y. (SM. - [WP + ^ (SM-WP) ] ) Az. (2'7'
J THICK J
where
AW = water available for evapotranspiration [cm]
THICK = the thickness of the uppermost layer [m]
WP = wilting point moisture content in the upper soil layer
[dimensionless fraction]
SM = soil water content in the upper soil layer [dimensionless
fraction]
AZj = thickness of segment j [m]
SMj = soil water content in segment j [fraction]
z.j = depth to the center of segment j [m]
The potential evapotranspiration and the available soil moisture are compared
and the lesser of the two amounts is assigned as the actual
evapotranspiration. Thus, the actual evapotranspiration for an event during a
specific season is:
where
ET = min (EPET, AW) (2.8)
ET = evapotranspiration during an event in a specific season [cm]
10
-------
Finally, the evapotranspiration rates for all of the seasons are aggregated to
estimate a representative steady-state value.
2.2.2 Percolation through the Landfill
Percolation rates through the landfill are determined by treating the landfill
layers and the unsaturated soil layers as one continuous system and applying
the one-dimensional Richard's Equation for unsaturated flow. The governing
equations and solution procedure are described in Section 3.2. The results of
the unsaturated flow module calculations are the soil moisture and pressure
head profiles in the landfill and the unsaturated zone for a given percolation
rate. The soil moisture values are used as input into the runoff and
evapotranspiration calculations, described in Section 2.2.1. Thus, the soil
moisture values impact the surface water balance and the leaching rate
predicted by the landfill module. Section 2.2.4 provides a step-by-step
description of how the unsaturated zone module and the landfill module are
coupled.
2.2.3 Landfill Control Features
Engineering controls can be used in landfill design to reduce the amount of
leachate emanating from the base of the landfill. The two important types of
engineering controls that may be simulated using this module are:
Low-permeability synthetic liners
Lateral drainage systems
The user may indicate the presence of synthetic liners and a drainage layer at
any location within the landfill as shown in Figure 2.1. One synthetic liner
per landfill layer may be included; however, only one lateral drainage layer
can be simulated. In addition, the user can elect to specify failure of the
synthetic liners. The techniques used to implement the effects of liners and
drains are discussed in the following paragraphs.
2.2.3.1 Low-permeability Liners
Low-permeability liners are thin sheets of rubber or plastic material used as
barriers to vertical flow. Liners are generally installed immediately above
or within barrier soil layers and together with the barrier soil layer reduce
vertical flow. Because the liner itself is too thin to treat numerically as a
separate layer, its hydraulic properties are combined with those of the soil
layer below it and both units are treated as one layer by the module. That
is, the effective hydraulic conductivity of the soil layer below is reduced
when using a synthetic liner. According to Freeze and Cherry (1979), the
effective hydraulic conductivity of the liner and adjacent soil layer is:
(TL + Ts)
KLS ^~ (2.9)
1± + _1£
KL Ks
where
KLS = effective conductivity of the liner and soil layer [m/yr]
KL = hydraulic conductivity of the liner [m/yr]
Ks = hydraulic conductivity of the soil layer [m/yr]
TL = thickness of the liner [m]
Ts = thickness of the soil layer [m]
11
-------
The user can elect to specify the degree of failure of the synthetic liner
system. To implement a liner failure, the user selects a percent failure (0
to 100 percent). The hydraulic conductivity of the liner/soil system is
increased by the selected failure percentage. That is:
FPERC , , ,
*' = JC"'ToO <*«-**> (2-10>
where
FPERC = percent liner failure [dimensionless]
KF = hydraulic conductivity of failed liner and soil layer [m/yr]
Note that the maximum hydraulic conductivity (at 100 percent failure) is the
hydraulic conductivity of the soil layer alone.
2.2.3.2 The Drainage System
Lateral drains are commonly utilized in landfill design to remove excess water
which may accumulate above barrier soil layers. Therefore, they serve to
reduce percolation through the landfill. The landfill module is capable of
simulating the effect of one drainage system placed anywhere in the landfill
profile. Generally, such a drainage system would be placed immediately above
or below the refuse layer.
The presence of lateral drains affects the surface infiltration rate, the
percolation rate to the water table, and the rate of lateral drainage. These
three effects are simulated by an algorithm which is based on the assumption
that the steady-state pressure head at the drain is less than or equal to
zero. That is, the drains operate at full efficiency and remove all ponded
infiltration above the drainage layer.
The first step in the lateral drainage algorithm is the simulation of the
moisture distribution within the landfill and unsaturated zone without
considering the effect of the drainage system. The pressure head at the
intended drain location is checked. If the pressure head indicates
unsaturated conditions, the module assigns zero as a lateral drainage rate.
If the pressure head value indicates saturated conditions at the drain
location, the pressure head at the drain is reduced to zero, and the system is
resimulated. The zero pressure head at the drainage layer serves as a
boundary condition, that is:
l|/(LD) = 0.0 (2.11)
where
LD = location of the lateral drain [cm].
This boundary condition allows the landfill and the unsaturated zone to be
broken into two independent unsaturated flow systems--one above the drainage
layer and the other below the drainage layer.
The system above the drain is solved in the same manner as the unsaturated
flow system which was described in Sections 2.2.1 and 2.2.2. The only change
in the solution procedure is that the location of the bottom boundary
condition is changed from the water table to the drainage layer. The result
of this calculation is the infiltration rate into the top of the landfill.
12
-------
This rate is generally larger than it would be in the absence of the drainage
system.
The system below the drain is solved in a similar manner, but Equation 2.11 is
used as a top boundary condition. Therefore, the surface water balance
calculation described in Section 2.2.1 is not used in the calculation. That
is, the percolation rate under the drainage layer is not significantly
impacted by surface water hydraulics and the unsaturated flow equations in
Section 3 are solved directly. The calculated percolation rate underlying the
drainage system is less than or equal to what it would be in the absence of
the drainage system.
The rate of lateral drainage is calculated based on a mass balance approach.
The lateral drainage rate plus the rate at which percolation reaches the water
table (Q£) must be equal to the infiltration (I) rate into the top of the
landfill. Solving for the lateral drainage results in:
QLAT = I - Q, (2.12)
where
QLAT = lateral drainage rate [cm/d]
I = infiltration rate [cm/d]
Q£ = percolation rate [cm/d]
2.2.4 Solution Method for Identifying Steady-State Leaching Rate
Section 2.2.1 describes the estimation of the infiltration rate (I) into the
landfill, based on a water balance method. Section 2.2.2 describes the
estimation of the percolation rate (Q) through the landfill, using the
unsaturated zone module described in Section 3. Section 2.2.3 explains how
landfill control features (liners and drains) affect the surface infiltration
rate and percolation rates from the bottom of the landfill. At steady-state
the percolation rate from the bottom of the landfill plus the lateral drainage
rate must be exactly equal to the infiltration rate into the landfill as there
can be no net change in moisture storage within the landfill.
The infiltration and percolation rates are coupled through the soil moisture
content of the upper portion of the landfill. This section describes the
iterative procedure used to estimate the steady-state leachate rate (and the
soil moisture profile) starting with an initial estimate of percolation. If
the calculated infiltration rate is different than the estimated percolation,
the estimate of percolation is adjusted, and the moisture content profile and
infiltration rate are recalculated. This procedure is repeated until
sufficient agreement between the percolation and infiltration rates is
attained and a solution which satisfies the surface water budget and the
unsaturated flow equations is found.
The procedure to calculate leachate rates consists of 10 iterative steps:
1. Calculate the effective hydraulic conductivity of the synthetic
liner/soil layers. Find the layer with the smallest hydraulic
conductivity in the landfill profile. Estimate an initial
percolation rate (Q) which is equal to the saturated vertical
hydraulic conductivity of this layer. This initial estimate is not
critical to the final solution, but simply allows the module to
start computations with a reasonable value.
13
-------
14
-------
2. Calculate the pressure head in each of the landfill and unsaturated
layers between the water table and the soil surface. Convert the
calculated pressure heads to soil saturation. These calculations
are both performed in the unsaturated zone module described in
Section 3.
3. Calculate the soil moisture using the relationship that soil
moisture is equal to soil saturation multiplied by effective
porosity.
4. Calculate the depth weighted retention parameter from Equation 2.5.
Calculate the surface runoff (RO) using Equation 2.3.
5. Calculate the actual evapotranspiration (ET) using Equations 2.6 to
2.8.
6. Calculate the surface water infiltration rate (I) using Equations 2.1
and 2.2.
7. Compare the infiltration rate (I) and the percolation rate (Q) assumed
in Step 1.
a. If I - Q|/Q < the convergence criteria, proceed to Step 9. The
convergence criterion is 1CT4 cm/day.
b. If I > Q, increase Q. This corresponds to a case where
infiltration from precipitation is sufficient to create perched
water above the least conductive layer so that more percolation
is forced through the layer. Portions of the profile will
remain saturated under these conditions.
c. If I < Q, decrease Q. This corresponds to case where
infiltration from precipitation is so low that even the least
conductive layer becomes unsaturated and less percolation moves
through it.
8. With the new updated estimate of Q, apply Steps 2 through 8
iteratively until convergence is achieved or the maximum number of
iterations is exceeded. If the maximum number of iterations is
exceeded without meeting the convergence criteria in Step 7, check
the relative convergence between Q and I. If *I-Q*/Q < the relative
convergence criteria proceed to Step 9. The relative convergence
criteria is 1CT3. Otherwise, print a convergence error message and
stop.
9. Check the lateral drainage flag. If there is an active lateral
drain in the profile, recalculate the percolation rate using the
lateral drain option.
10. Repeat the entire procedure for each of the seasons. Output results
for each season and convert the calculated seasonal infiltration
rates to a steady-state average.
15
-------
2.3 Assumptions and Limitations
The important simplifying assumption made in developing the landfill module
include:
a) Lateral inflow, surface run-on, and transient climatic conditions
are considered negligible. Degradation or aging of design
components is not considered. However, the degree of failure of the
liners may be specified.
b) The seasonal potential evapotranspiration rates are assumed to
adequately characterize the dynamic variability of actual
evapotranspiration demand.
c) The event-averaged infiltration rate is an adequate representation
of the long-term average infiltration rate.
d) Groundwater is assumed to be below the bottom of the landfill.
e) No liquid is generated from the decomposition of waste or codisposal
of wastes and sludges.
f) Lateral drains are fully efficient and remove all ponded water above
the drainage layer.
2.4 Data Requirements
The landfill module has relatively modest data requirements. Table 2-1 lists
the parameters required.
16
-------
TABLE 2-1. PARAMETERS REQUIRED FOR THE LANDFILL MODULE
Parameters
Units
Climatic Data
For each season:
Typical storm precipitation amount
Typical interval between storms
Potential evapotranspiration
Design Specifications
SCS Curve Number II
Number of liners
Number of layers
Number of porous materials
Number of seasons
Location of the drainage layer
Thicknesses of the layers
For each liner:
Percent liner failure
Thickness of liner
Hydraulic conductivity of liner
For each porous material:
Saturated hydraulic conductivity
Porosity
Residual water content
Either:
van Genuchten alpha coefficient
van Genuchten beta coefficient
or
Brooks and Corey exponent
van Genuchten alpha coefficient
van Genuchten beta coefficient
[cm]
[d]
[cm]
[dimensionless]
[dimensionless]
[dimensionless]
[dimensionless]
[dimensionless]
[dimensionless]
[m]
[percent]
[m]
[m/yr]
[cm/hr]
[dimensionless]
[dimensionless]
[I/cm]
[dimensionless]
[dimensionless]
[I/cm]
[dimensionless]
17
-------
SECTION 3
The Unsaturated Zone Flow Module
3.1 Introduction
When the bottom of a waste disposal unit is located above the water table,
leachate can migrate through the unsaturated zone and into a saturated
aquifer. In such situations it is important to include the unsaturated zone
in the analysis of contaminant fate and transport. A schematic diagram of the
leachate migration is shown in Figure 3.1.
This chapter presents details of the semi-analytical unsaturated zone flow
module included in MULTIMED. The flow module computes the water saturation
values within the unsaturated zone which are used by the unsaturated zone
transport module (described in Section 4) to compute one-dimensional vertical
seepage velocities.
Two options, previously described in Section 2, are available for the
simulation of leachate percolation through the unsaturated zone. Under Option
1, the user specifies the leachate rate, and the unsaturated zone module
simulates the percolation of leachate to the water table. Under Option 2, the
leachate rate is not specified, rather it is calculated by the landfill module
(see Section 2). This calculation is dependent on the results of the
unsaturated zone module. Therefore, an iterative scheme is used to calculate
a leaching rate which satisfies both the surface water balance included in the
landfill module and the unsaturated zone water saturation profile estimated by
the unsaturated flow module.
This chapter is organized into four sections. Theoretical details of the flow
module and the underlying assumptions are presented in Sections 3.2 and 3.3.
Section 3.4 discusses the data requirements for this module. Note that the
unsaturated flow and transport modules can not be run independent of each
other. In addition, whenever the unsaturated flow and transport modules of
MULTIMED are used, the saturated transport module must also be used.
3.2 GOVERNING EQUATIONS AND SOLUTION TECHNIQUES
The unsaturated zone flow module simulates steady downward flow to the water
table. The governing equation is given by Darcy's law:
-Kvkrw (!L - i) = Q (3.D
18
-------
where
i|;
z
K^
kr
Q
the pressure head [m]
the depth coordinate which is taken positive downward [m]
the saturated hydraulic conductivity [m/yr]
the relative hydraulic conductivity [dimensionless]
the percolation rate [m/yr]
The boundary condition at the water table is:
= 0
(3.2)
where L, is the thickness of the unsaturated zone [m] .
To solve the above equations, it is necessary to specify the relationships of
relative hydraulic conductivity (krw) versus water saturation (SJ , and of
pressure head (i|/) versus water saturation. The relationship of pressure head
to water saturation is described in the model by the following equation (van
Genuchten, 1976; Mualum, 1976):
Se =
(3.3)
where
Swr = the residual water saturation [dimensionless fraction]
3,Y = soil-specific empirical parameters [dimensionless]
a = soil-specific empirical parameter [1/m]
i|/a = the air entry pressure head, subsequently assumed zero [m]
Se = the effective saturation [dimensionless fraction]
and where Se is related to the water saturation, Sw, as follows:
The parameters a, 3, and y are empirical coefficients defined by van
Genuchten. The parameters 3 and y are related through:
S -
(34)
(3.4)
Y = 1 - 1/3
(3.5)
19
-------
Hence, the parameter y need not be specified. Note that Equations 3.3 and 3.4
are not valid for 3 < 1.0 and y < 0.0 because the effective saturation, Se,
can not be greater than one.
Two alternative functional expressions can be specified by the user to
describe the relationship between relative hydraulic conductivity and water
saturation. The van Genuchten relationship is:
krw = Se1/2[l- (1 - Se1/Y)V (3.6)
Alternatively, the krw(Sw) relationship presented by Brooks and Corey (1966)
may be used. This relationship is given by:
krw=Sen (3.7)
where n = soil specific parameter [dimensionless] .
As a first step in the solution of Equations 3.1 and 3.2, the soil
constitutive relations in Equations 3.3 and 3.6 are combined. Using van
Genuchten Ts constitutive equations and assuming i|/a = 0 , leads to the following
expression for krw(i|/) :
C 1 i|; ;> 0
krw = \ (3.8)
,{l - (-ml/J^tl + (-ml/)"]1''1'-1}2
I))))))))))))))))))))))))))))))))) Hf < 0
Next, Equation 3.8 is substituted into Equation 3.1 and the derivative di\i/dz
is replaced by a backward finite difference approximation. This yields, after
some rearranging:
)) ()))))))))) + 1) - 1 = 0 Hf * 0
Q Az
(3.9)
K, {l - (-on|f ^[i + (-onlf)"]1"-1}2 ^z_iz - ^z
)) )))))))))))))))))!)))))))))))))) ()))))))))) + i) - i = o nf < o
Q [l + (-ml/)"] (1/2-1/2p) AZ
where i|; is the representative pressure head for the soil layer between z and
z-Az .
If the Brooks and Corey (1966) relationship is used, the expression for
relative hydraulic conductivity becomes:
r 1 i|; ;> 0
(3.10) .
krw = i
1(1 + (-on|/)p)-nv i|; < 0
20
-------
Substituting Equation 3.10 into Darcy's law (Equation 3.1), the resulting
expression equivalent to Equation 3.9 is:
)) ()))))))))) + 1) - 1 = 0, i|r * 0
(3.11)
Q Az
)) (l + (-c^)prnv ())))))))))) + l) - l = 0, Hf < 0
Q AZ
i|; can be written as a weighted average of i|/z and i|/z_iz:
if =
-------
distance necessary to obtain the proper depth. Thus, in the above
example, the distance between nodes 1 and 2 would be 0.4.
The minimum nodal spacing is 0.1 m. If the distance between nodes
1 and 2 is less than O.lm, then the total number of nodes is
decreased by one and the distance between nodes 1 and 2 is
increased by 1 m.
(c) If the depth is greater than 200 m, the number of nodes is set
equal to 200 and all layers are of equal thickness.
(d) If multiple layers are used, an additional node is inserted to
define the bottom of each layer.
Note that if the depth of the unsaturated zone is generated randomly for each
run, the unsaturated zone is considered homogenous and composed of only one
material. Thus, only one layer and one material can be specified in both the
unsaturated flow and unsaturated transport modules.
3.3 ASSUMPTIONS AND LIMITATIONS
The major assumptions on which the flow model is based are:
(a) Flow of the fluid phase is considered isothermal, one-dimensional,
and governed by Darcy's law.
(b) The flow field is considered to be steady.
(c) The simultaneous flow of the second phase (i.e., air) can be
disregarded.
(d) Hysteresis effects are neglected in the specification of the
characteristic curves.
3.4 DATA REQUIREMENTS
The data required by the unsaturated zone flow module are listed in Table 3-1.
Note that the van Genuchten parameters are required by the model to describe
the pressure head versus water saturation relationship. The user can specify
one of two options for describing the relationship between relative
permeability and water saturation: the van Genuchten or the Brooks and Corey
equation. If the Brooks and Corey option is used, a Brooks and Corey exponent
must be input in addition to the van Genuchten alpha and beta coefficients.
22
-------
TABLE 3-1. Parameters required for the Unsaturated Zone Flow Module
Parameter
Units
Source-Specific Data
Percolation rate from the facility
Unsaturated Zone Data
Number of layers
Number of porous materials
Thickness of each layer
Material associated with each layer
For each material:
Air entry pressure head
Porosity
Saturated hydraulic conductivity
Residual saturation (water content)
Either:
van Genuchten alpha coefficient
van Genuchten beta coefficient
or
Brooks and Corey exponent
van Genuchten alpha coefficient
van Genuchten beta coefficient
[m/yr]
[dimensionless]
[dimensionless]
[m]
[dimensionless]
[m]
[dimensionless]
[cm/hr]
[dimensionless]
[I/cm]
[dimensionless]
[dimensionless]
[I/cm]
[dimensionless]
Note: The model provides the option to use either van Genuchten's or Brooks
and Corey's constitutive relationship for relative permeability versus
water saturation. However, the relationship between pressure head and
water saturation is expressed in terms of van Genuchten parameters.
23
-------
SECTION 4
Unsaturated Zone Transport Module
4.1 Introduction
This section presents the details of the unsaturated zone transport module
included in MULTIMED. As mentioned in Section 3, transport within the
unsaturated zone is important only when the bottom of the waste disposal unit
is located above the water table. Also, when the unsaturated zone modules are
used, the saturated zone module must also be used.
The theoretical basis of the unsaturated zone transport module and the
underlying assumptions are discussed in Sections 4.2 and 4.3. Section 4.4
addresses the data requirements for this module.
4.2 Governing Equations
4.2.1 Unsteady-State Transport
The transport of contaminants in the unsaturated zone is treated as a one-
dimensional problem. Important fate and transport mechanisms considered by
the model include dispersion in the vertical direction, linear adsorption, and
first-order decay of the contaminant. With these constraints, the transport
equation can be expressed as:
where
C = the dissolved phase contaminant concentration in the
unsaturated zone [mg/f] ,
Dv = the dispersion coefficient in the unsaturated zone [m2/yr]
X,f = the first-order degradation rate within the unsaturated zone
[1/yr]
Rv = the unsaturated zone retardation factor [dimensionless]
vv = the steady-state unsaturated zone seepage velocity [m/yr]
t = time [yr]
z = the vertical coordinate which is positive downwards [m]
24
-------
The retardation factor in Equation 4.1 is computed using:
where
pbv = the bulk density of the unsaturated zone [g/cc]
Kdv = the contaminant distribution coefficient for the unsaturated
zone [cc/g]
6 = the porosity of the unsaturated zone [dimensionless fraction]
Sw = the saturation within the unsaturated zone [dimensionless
fraction]
The overall first-order degradation rate, X,f, which is calculated using
Equation 5.3 in Section 5, includes the effect of both biodegradation and
chemical hydrolysis reactions. The latter is discussed in Section 5.5.2.
Further, the unsaturated zone seepage velocity in Equation 4.1 is computed
using :
V =
QS
(4.3)
where Q is the steady-state percolation rate within the unsaturated zone.
Note that Q is assumed to be steady in MULTIMED. Also, the saturation, Sw, is
computed by the unsaturated zone flow module, as discussed in Section 3.
Solution of the above differential equation requires two boundary conditions.
The first boundary condition describes the source concentration and may be of
the following form:
C(0,t) = Co (4.4a)
or
C(0,t) = Co exp(-At) (4.4b)
or
C(0,t) = Co[l-s(t-T)] (4.4c)
25
-------
where
A = the source concentration decay rate [1/yr]
s(t-T) = the unit step function with a value of unity for t>T and zero
for t
-------
C± (l^t) = Ci+1 (0,t)
(4.10;
27
-------
where I = the thickness of a layer and the subscripts i and i+1 refer to
successive layers. Equation 4.10 implies that the source concentration at the
top of any layer i+1 is set equal to the concentration computed at the bottom
of the previous layer i. Note that the layers can be of different thickness.
The solution to the layered unsaturated zone is derived using Laplace
transform techniques to transform the governing partial differential equation
(Equation 4.1) and the boundary conditions to an ordinary differential
equation in the Laplace domain. The ordinary differential equation is solved
in the Laplace domain and then inverted using either the convolution theorem
or the Stehfest algorithm (Stehfest, 1970; Moench and Ogata, 1981). The
latter is a numerical inversion scheme. Both these solution schemes are
included in the model. In general, the Stehfest algorithm is computationally
faster. At very high Peclet numbers, however, there is a possibility that
this numerical solution may not converge. For such cases, the convolution
integration method may be used. Details of the solution scheme are presented
by Shamir and Harleman (1967) and Hederman (1980).
4.2.2 Steady-State Transport
For the case of a steady-state continuous contaminant source, the governing
Equation 4.1 can be simplified to yield:
For this case the boundary conditions are:
C(z=0) = Co (4.12a)
(z=°°) = 0 (4.12b)
G Z
The analytical solution to the above system of equations is:
V z
C(z) = Co exp {-- - z(XvRv/Dv + V2/4D2)1/2} (4.13a)
or
C(z) = Coexp{_ - -. (1 + ji/*} (4.13b)
-X LR
-^) (4.14;
28
-------
In the event that dispersion within the unsaturated zone is neglected, the
above equation reduces to:
where L = the depth of the unsaturated zone [m] .
For a layered unsaturated zone, Equation 4.14 can be expressed as:
(415,
where n = the number of homogenous layers within the unsaturated zone.
4.3 Assumptions and Limitations
The major assumptions on which the unsaturated zone transport model is based
are:
(a) The flow field within the unsaturated zone is at steady state.
(b) The seepage velocity and other model parameters (e.g., the
dispersion coefficient, partition coefficient) are uniform in each
layer. In other words, each layer is homogeneous and isotropic.
(c) Transport is assumed to be strictly one dimensional. Lateral and
transverse advection and dispersion are neglected.
(d) Adsorption and decay of the solute may be described by a linear
equilibrium isotherm and a first-order decay constant respectively.
The daughter products of chemical and biochemical decay are
neglected.
(e) Each layer is approximated as being infinite in thickness. This
assumption is valid and introduces negligible errors if the ratio
of dispersivity to the layer thickness is small (« 1.0).
4.4 Data Requirements
Table 4-1 lists the parameters required by the unsaturated zone transport
module. These include:
(a) Three source-specific parameters. Note that the model is linear
with respect to the source concentration. Thus, if the source
concentration is set to unity, the model computes normalized
downgradient well concentrations.
(b) Five chemical-specific parameters. Note that the overall decay
coefficient and the distribution coefficient for the unsaturated
zone can not be input directly. Rather they are computed from
other parameters as discussed is Sections 4.4.1 and 4.4.2.
(c) Seven unsaturated zone transport parameters. Remember that the
seepage velocity is computed using Equation 4.3, with the
saturation values supplied by the unsaturated zone flow module.
The computation of the dispersion coefficient by the code is
discussed in Section 4.4.3.
29
-------
(d) Two aquifer parameters. The pH and temperature of the unsaturated
zone is assumed to be the same as for the saturated zone.
4.4.1 The Chemical Transformation Rate
The liquid-phase and sorbed-phase chemical decay coefficient are computed
using the hydrolysis rate constants. The equations used are shown in Section
5.5.2.1. Note that the pH and temperature of the unsaturated zone can not be
input; rather, the corresponding values for the saturated zone are used. The
overall decay rate is then computed using Equation 5.3 in Section 5.
4.4.2 The Distribution Coefficient
The distribution coefficient is computed as the product of the normalized
distribution coefficient and the fractional organic carbon content. However,
the fractional organic carbon content is replaced by the percent organic
matter divided by a conversion factor, using the following relationship
(Enfield et al. , 1982) :
f~=T&74 (4'16)
where
foc = fractional organic carbon content [dimensionless]
fom = percent organic matter content [dimensionless]
172.4 = conversion factor from percent organic matter content to
fractional organic carbon content
30
-------
TABLE 4-1. Parameters Required for Unsaturated zone Transport Module
Parameters Units
Source-Specific Parameters
Source decay constant
Source concentration at top of unsaturated zone
Pulse duration (for unsteady state simulation)
Unsaturated Zone-Specific Parameters
Number of layers used to simulate transport
Porosity of the unsaturated zone (specified in
the unsaturated flow input)
For each layer:
Thickness
Longitudinal dispersivity
Bulk density of the soil
Biodegradation rate,
Percent organic matter
Saturated Zone-Specific Parameters
Temperature of the aquifer3
pH of the aquifer3
Chemical-Specific Parameters
[i.e., Koc)
Normalized distribution coefficient
Reference Temperature
Acid and base hydrolysis rates at reference
temperature
Neutral hydrolysis rate at reference temperature
[1/yr]
[mg/f]
[yr]
[dimensionless]
[cc/cc]
[m]
[m]
[g/cc]
[1/yr]
[dimensionless]
[pH units]
[cc/g]
[f/mole-yr]
[1/yr]
* Note that the temperature and pH used in calculating the unsaturated zone
overall decay rate are the temperature and pH specified for the aquifer.
4.4.3 The Longitudinal Dispersion Coefficient
Longitudinal dispersion is computed using the relationship:
DV = avvv (4.17)
where
Dv = the longitudinal dispersion coefficient [m2/yr]
Vv = the seepage velocity in the unsaturated zone [m/yr]
av = the longitudinal dispersivity [m]
The longitudinal dispersivity, av, can be either input directly by the user or
derived by the code. The equation used in the model to derive dispersivity
values is based on an analysis of data presented by Gelhar et al. (1985) and
shown in Table 4-2. Using regression analysis, the following relation was
developed:
31
-------
a = .02 + .022L, R2 =
where L is the thickness of the unsaturated zone. To avoid excessively high
values of dispersivity for deep unsaturated zones, a maximum dispersivity of
1.0 m is used. Thus, for all depths greater than 44.5 m, av is set equal to
1.0 m. Additional details of the above regression are presented in Salhotra
(1988) .
TABLE 4-2.
Compilation of Field Dispersivity Values (Gelhar et al . ,
1985)
Author
Yule and Gardner
(1978)
Hildebrand and
Himmelblau (1977)
Kirda et al .
(1973)
Gaudet et al .
(1977)
Brissaud et al .
(1983)
Warrick et al .
(1971)
Van de Pol et al .
(1977)
Biggar and Nielsen
(1976)
Kies (1981)
Jury et al . (1982)
Andersen et al .
(1968)
Oakes (1977)
Type of Vertical Scale
Experiment of Experiment (m)
Laboratory 0.23
Laboratory 0.79
Laboratory 0.60
Laboratory 0 . 94
Field 1.00
Field 1.20
Field 1.50
Field 1.83
Field 2.00
Field 2.00
Field 20.00
Field 20.00
Longi tudinal
Dispersivity
av (m)
0 . 0022
0.0018
0.004
0.01
0.0011,
0.002
0.027
0.0941
0.05
0.168
0.0945
0.70
0.20
32
-------
SECTION 5
The Saturated Zone Transport Module
5.1
Introduction
This chapter presents details of the module used to simulate contaminant fate
and transport within the saturated porous zone. Recall that the contaminant
can enter the saturated formation by direct leaching from a source (in the
absence of an unsaturated zone) or by percolation through the unsaturated
zone. MULTIMED allows either option to be specified. In both cases the
governing equations and the semi-analytical solution for transport in the
saturated zone are the same. Flow in the saturated zone is assumed to be
steady.
The following sections describe the governing equations, boundary and initial
conditions, model limitations, and parameters associated with the saturated
zone transport module.
5.2 Governing Equations
The three-dimensional solute transport equation on which the model is based
can be written as:
By
- v
x c
B B
where
x, y, z = spatial coordinates in the longitudinal, lateral and vertical
directions, respectively [m]
C = dissolved concentration of chemical [mg/f, g/m3]
Dx, Dy, Dz = dispersion coefficients in the x, y and z directions,
respectively [m2/yr]
Vs = one dimensional, uniform seepage velocity in the x direction
[m/yr]
Rs = retardation factor in the saturated zone [dimensionless]
t = elapsed time [yr]
Xs = effective first-order decay coefficient in the saturated zone
[1/yr]
q = net recharge outside the facility percolating directly into and
diluting the contaminant plume [m/yr]
B = the thickness of the saturated zone [m]
6 = effective porosity for the saturated zone [dimensionless]
In Equation 5.1, the retardation factor and the effective decay coefficient
are defined as:
R = 1 +
(5.2;
33
-------
and
where
pb = bulk density of the porous media [g/cc]
Kd = distribution coefficient [cc/g]
A! = first-order decay constant for dissolved phase [1/yr]
X2 = first-order decay constant for the sorbed phase [1/yr]
Xb = first-order lumped biodegradation rate in the saturated zone
[1/yr]
The flow domain is regarded as semi-infinite in the x direction (0 < x < <*>) ,
infinite in the y-direction (-<*> < y < <*>) and finite in the z-direction (0 < z
< B) .
Solution of Equation 5.1 requires two boundary conditions in each of the x, y,
and z directions and an initial condition. At the source (downgradient edge
of the waste disposal unit) the model allows a choice between two boundary
conditions with respect to the distribution of contaminant along the vertical
source plane. The first boundary condition specifies the contaminant
concentration as a gaussian distribution in the lateral direction and uniform
over the vertical mixing or source penetration depth, H. A schematic
description of the flow domain and this source boundary condition is shown in
Figure 5.1(a). Mathematically, the boundary condition can be expressed as:
rco exp [-y2/ (2a2,) ] 0
-------
Changes in the maximum dissolved concentration, C0, over time are modeled
using one of three options. The value of C0 can be constant for all time
(steady-state) :
C0 = C, t>0 (5.5)
or a finite pulse, Ts years in duration, of constant concentration:
C0 = C, 00 (5.7)
where
C, = user-specified leachate concentration [mg/f]
A = source decay rate constant [1/yr]
The second, alternative source boundary condition allowed by the model is a
rectangular patch source of thickness H and width 2y0. A schematic diagram of
this boundary condition is shown in Figure 5.1(b) . It can be mathematically
expressed as:
rco , -Yo^y^Yo and 0
-------
In addition to the source boundary condition, the initial and additional
boundary conditions used by the model can be expressed as:
C(x,y,z,t) = ^Cf(x,y,t) + ACp (x, y, z, t) (5.10]
Equation 5.9a implies that background concentrations of the contaminant in the
aquifer are zero. Equations 5.9d and 5.9e imply that there is zero flux of
contaminant at z = 0 and z = B.
5.2.1 Solution to the Gaussian Source Boundary Condition
Huyakorn et al. (1987) and USEPA (1985) have presented analytical solutions
for the system of Equations 5.1 to 5.4 and Equation 5.9. The general solution
can be expressed as:
where C£ and ACP are functions given by:
Cf(x,y,t) = ^tF(x, y, T) exp (~r\i] di (5.11)
. -, , . 2 £ v^ 1 , nnz. . , nnH.
AC (x, y, z, t) = - 2_, cos ( ) sin ( )
P n n.i n B B (5_12;
/ t F(x,y, T) exp (-3nT) di
J o
in which
(2nD*
F(x'y'T) =
C ax V'x
exP ( - 7) (5.141
'
36
-------
(5.15;
where
D*, D*, and D* = the retarded dispersion coefficients (DX/RS, Dy/Rs,
DZ/RS) in the x, y and z direction
V* = the retarded solute (seepage) velocity [vs* = VS/RS]
T = the variable of integration
Note that in the event that H = B (i.e., the source fully penetrates the
saturated formation), ACP = 0 in Equation 5.10. At any distance, x, from the
source, the maximum contaminant concentration occurs at the centerline of the
plume and can be represented as :
C(x, 0,0, t) = Cf(x, 0,t) + ACp(x, 0, 0, t) (5.17)
where C£(x,0,t) and ACP (x, 0, 0, t) are given by Equations 5.11 and 5.12 with
arguments y and z set equal to zero, and the function F(x,0,i) defined as:
exp
F(x,0,i) =
T3/2
As t approaches infinity, a steady-state condition is reached. The steady-
state concentration along the plume centerline can be expressed as:
C*(x, 0,0) = -| c; (x, 0) + ACp* (x, 0,0) (5.19)
where
C'(x, 0) = 5* f"exp [- °!ji! - x (?-%L + JL)1/2 ]du (5.20)
Jo 2 * *
X X
AC* (x, 0,0)= -^- £) sin
n £=i n B
X X
and
2 Co V*x
(5-22:
37
-------
The solution for the transient state (i.e., Equations 5.10 to 5.16) in
MULTIMED is achieved by incorporating a pre-existing code, EPATMOD (Huyakorn
et al., 1986). Similarly, the steady-state solution (i.e., Equations 5.19
through 5.22) in MULTIMED makes use of a pre-existing code named EPASMOD
(Huyakorn et al. , 1986).
Note that for large pulse durations, the transient solution is identical to
the steady-state solution. However, the steady-state solution is
significantly faster than the transient solution and should be used for
steady-state computations. Finally, note that the model uses the principle of
superposition to compute the plume concentration for a pulse source, i.e., a
contaminant source of finite duration, Ts.
5.2.2 Solution to the Patch Source Boundary Condition
Sudicky et al. (1988) have presented the analytical solution to the system of
Equations 5.1 through 5.3, 5.8, and 5.9. The general solution can be
expressed as:
C ( ~>c \r 7 i~}
Irrfr I
T3/2 Lcrrct^
- c x
0 1/2
X
Fprfr 1 '
r 21 / t i s , / . ,
y y y+ y
} rrfr { \ 1 r?T
rt (x-V*T)2 !
Jo 4D*T T '
x
> } -prfr i ^ y° ll V -1
(5.23a)
nnz nnz nllz n2n2D*T
[sin { } - sin { }] cos ( ) exp { - } di
' v
38
-------
where
z = 0
(5.23b)
and
=H
(5.23C)
If the solution is desired along the plane y
simplified by noting that the term
= 0, Equation 5.23a can be
erfc {
Y~Y
, 1/2
- erfc
1/2
(5.24a)
reduces to:
2-2 erfc {
2(Dy*T)1/2
(5.24b)
When z = 0, the term cos(nnz/B) becomes 1. A more substantial simplification
results in the case of a fully penetrating source (i.e., z{ = 0 and z2 = B) .
In this case, the entire second integral term in Equation 5.23 vanishes.
Similarly, for the steady-state case the solution is:
S~1 ( ^.r , - ~\ S~1 2 1 / \ i T7
x2 , (y-y"
2C x V'x ^ !
+ exp ( ) 2^ [si
*- ,-. _ # ^-~ ,-.
n
, 1/2
in
y
X D;
nnz nnz
) - sin ( ) ] cos
(5.25}
2* 22*21/2
X
M D*
N x
[
v2
i\ + y*//ip* + n2n2P*/R2) ( +
^/v T vs /t^x -" 11 vz/r> 1 \ D«
JC
, (y-y*)2
D*
y
(y-y*^ j
D' '
dy*
39
-------
where Kx is the modified Bessel function of the first kind and first order.
Note that the second integral term in Equation 5.25 vanishes again in the case
of a fully penetrating source.
An alternative formulation of the steady-state solution for the case where y =
0 is :
C(x, 0,z) = °Z\ZX (X + V2/4D;)1/2 exp (-) f" -i
2(nD')1/2B 2D' J° T3/2
exp { - T - ^-} erf (-^) di
(5.26;
V*x ^ ]_ nllZ2 nnz n z
exp ( r) / ., [sin (-) -sin(-)] cos (
n(nD')1/2 2D- ~ n ' B ' B ' B '
X X
(X + Vs*24D; + n'rfD.VBV" -1 exp { - T - } erf ( ) dt
with
(5.2?;
y
(5.28]
,X+V'2/4D + n2n2D7-B2,
- ^ - - - } (5.29)
As before, for the case of a fully penetrating source and/or when z = 0, this
solution can be further simplified. The above solutions, Equations 5.23 to
5.26, were earlier programmed in a code named PATCH3D. This code has been
incorporated into MULTIMED.
5.2.3 Receptor Well or Stream Location
A schematic of the receptor location relative to the waste facility is
presented in Figure 5.2. The location of a well is determined by specifying
in the input the radial distance to the well, the angle between the plume
centerline and the radial location of the well measured counterclockwise, and
the depth of penetration of the well. The code computes the cartesian
coordinates of the well location as:
xr = R Cos i}f (5.30)
40
-------
yr = R Sin ty (5.31)
where
R = the radial distance to the well [m]
i|; = the angle measured counterclockwise from the plume centerline
[degrees]
xr, yr = the cartesian coordinates of the well location [m]
The z coordinate of the well is specified in the input as a fraction of the
total depth of the aquifer, measured from the water table downward. Based on
the input, the code calculates the z-distance in meters from the water table
to the well. The well is assumed to have a single slot at that depth.
when the surface water module is used, a stream is the receptor instead of a
well. The stream location is determined from the radial distance to the
receptor alone (i.e., the angle off-center and z-distance are not used).
5.2.4 Time Values for Saturated Zone Concentrations
When operated in transient mode, MULTIMED requires that time values for
calculating the groundwater concentration at the well be entered by the user.
(The steady-state model does not require time values.) The source term for
the transient model can be either a square wave pulse, formed by a constant
leachate concentration of finite duration, or an exponentially decaying pulse,
formed by specifying a source decay rate. In either case, the concentration
at the well builds gradually, reaches some peak value, and declines as the
dispersed pulse passes the well. The user must enter the time values that
capture the various components of the pulse as it passes the well, such as the
break-through time or the time of maximum concentration. An analytical
technique to exactly predict when the pulse will reach the well or when the
maximum concentration will occur has not been developed. However, the
following equations may provide a reasonable estimate of the time, Tmc, when
maximum concentration will be reached:
Tmc = xRs/Vs + 0.5TS (5.32)
where
x = x-distance to the well
Rs = retardation factor in the saturated zone
vs = groundwater seepage velocity
Ts = duration of the leakage pulse
Using this as a starting time, the user can begin a search that should lead
quickly to the time values of interest. A good time step, Tstep, to use in the
search is given by:
Tstep = 0.4Tmc (5.33)
If the unsaturated zone is included in the simulation, Ts should be set to the
time when the unsaturated zone concentration has passed its maximum and is
about half way back to zero. This can be found by making an initial model run
and examining the resulting unsaturated zone time series found in the output
file VTRNSPT.OUT.
41
-------
5.3 Assumptions and Limitations
Following are the list of assumptions inherent in the saturated zone transport
module:
(a) A single aquifer with uniform thickness is modeled. The saturated,
porous medium properties are isotropic and homogeneous. The module
cannot be used to simulate transport in fractured media unless the
fractured medium is represented as an equivalent porous formation.
(b) The ground water flow velocity is steady and uniform. This implies
that the recharge through the facility and into the groundwater
plume is small compared to the natural (regional) flow.
(c) Contaminant degradation/transformation follows the first-order rate
law and is restricted to biodegradation and hydrolysis. The latter
is a second order process from which the first order rate is
obtained using the existing environmental pH. This assumption is
conservative since it neglects degradation due to other mechanisms
such as oxidation, reduction, etc. Further, the model does not
include by-products of degradation.
(d) Contaminant sorption follows a linear adsorption isotherm.
Adsorption takes place instantaneously and the adsorbed phase is in
local equilibrium.
(e) The initial contaminant concentration in the aquifer is zero.
Assumptions regarding the source boundary conditions and the extent
of the formation have been discussed in Section 5.2.
5.4 Coupling of the Unsaturated and Unsaturated Zone Modules
When modeling the transport of contaminants through the unsaturated and the
saturated zones, an important requirement is that the principle of
conservation of mass be satisfied (i.e., the mass flux of contaminant that
leaches from the bottom of the unsaturated zone (or out of the facility in the
absence of an unsaturated zone) must be equal to the mass flux that enters the
saturated zone). This mass flux consists of the sum of advective and
dispersive mass fluxes.
5.4.1 Steady-State Coupling
The mass that reaches the water table from the facility can be expressed as:
(5.34)
where
ML = the mass that leaches out of the facility [g/yr]
Aw = the area of the facility [m2]
Q£ = percolation rate [m/yr]
C, = concentration in the leachate from the facility [g/m3] if
attenuation within the unsaturated zone is neglected or the
unsaturated zone is absent. Alternatively, C, is the
concentration at the bottom of the unsaturated zone.
42
-------
The mass flux that is advected into the saturated zone is calculated by
integrating the source concentration in the y direction from -<*> to +<*> and over
the depth z = 0 to z = H. Thus the mass flux advected into the aquifer is:
Ma = IH f~ C(x=0,y, z) VBQ dydz (5.35)
J z=0 Jy=-
where
Ma = mass flux advected into the aquifer [g/yr]
C(x=0,y,z) = concentration as a function of y and z at the source [g/m3,
mg/{] as expressed by Equation 5.4 or 5.8
vs = the seepage velocity in the saturated zone [m/yr]
6 = effective porosity of the saturated zone (cc/cc)
Similarly, the mass flux that enters the saturated zone due to dispersion can
be expressed as:
where
Md = mass flux dispersed into the aquifer [g/m3]
Dx = dispersion coefficient in the x direction [m2/yr]
5.4.1.1 The Gaussian Source Boundary Condition--
Substituting Equation 5.4 into Equation 5.35 and integrating, with C0 assumed
uniform over the source depth H, yields:
(5.36)
Ungs (1987) has evaluated the integral in Equation 5.36 to yield:
Ma = (2n)1/2aVs6ffC (5.37)
Md = (2n)1/2aVs8ffCo[ -1/2 + 1/2 (1 + £_|_J£)i/2] (5.38)
where
Xs = the first-order decay coefficient [1/yr]
Rs = the linear retardation factor [dimensionless]
Dx = the longitudinal dispersion coefficient [m2/yr]
Note that in the event, that Dx = 0, the dispersive flux, Md, is zero.
The total flux into the saturated zone is given by the sum of advective
(Equation 5.37) and dispersive (Equation 5.38) fluxes:
MT = (2n)1/2 aV3QHCot,D (5.39)
43
-------
44
-------
where
Xs = the first-order decay coefficient [1/yr]
Rs = the linear retardation factor [dimensionless]
Dx = the longitudinal dispersion coefficient [m2/yr]
Note that if £D is set equal to unity, it implies that the dispersive flux is
neglected.
Equating Equations 5.34 and 5.40 yields the following expression of the mass
balance:
AwQfCt = (2n)1/2 aVs6HCoCD (5.41)
The above equation is used to couple the unsaturated and the saturated zone
models under steady-state conditions.
Equation 5.41 can be rearranged to yield:
C°= (2^QV3QHO^DC> (5'42)
or
Co = (NMF)Ct (5.43)
The factor NMF can be thought of as representing a near-field dilution effect
or the effect of mixing below the facility. Based on purely physical
considerations, this factor should be less than or at most equal to unity.
5.4.1.2 The Patch Source Boundary Condition
Substituting Equation 5.8 into Equation 5.35 and substituting Equation 5.23a
into Equation 5.36 and integrating, the total mass entering the saturated zone
can be expressed as:
MT = 2yoVs8ffCo + 2yoVsQHCo {-1/2 + 1/2 (1 + ^ *)1/2} (5.44;
or letting
45
-------
then,
MT = 2yoVs8ffCo
(5.46;
In Equation 5.44 the first term represents the advective mass flux and the
second term represents the dispersive mass flux entering the saturated zone.
Equating Equation 5.46 to Equation 5.34 yields the mass balance expression
used to couple the unsaturated and saturated zone (patch source) model under
steady-state conditions:
C =
(5.4?;
5.4.2 Unsteady-State Coupling
For the case of unsteady-state transport in the unsaturated zone, the mass
flux at the water table varies in time and the above approach for coupling the
unsaturated and the saturated zone is no longer valid. In the unsteady state,
concentrations in the saturated zone are determined using the convolution
integral approach that superimposes the effects of source changes over time as
follows:
ac*
C(x, y,z,t) = f fc 4r- T f U, y, z, t - T
JO OT
dT
(5.48)
where
C*(t) = the concentration at the water table at time t [mg/f]
f(x,y,z,t) = the normalized (with respect to source concentration)
solution of the saturated zone analytical solution [mg/f]
In Equation 5.48, the value of f(x,y,z,t) is obtained using Equations 5.10 to
5.16 (with C0 = 1) for the case of a gaussian source boundary condition or
Equation 5.23a (with C0 = 1) for a patch source boundary condition. In the
computer code, the integral is numerically evaluated using the trapezoidal
rule.
5.5 Data Requirements
Table 5-1 lists the primary input parameters required to compute the
contaminant concentrations in the saturated zone. These parameters are
classified into the following groups:
(1) Source-specific parameters
(2) Chemical-specific parameters
(3) Aquifer-specific parameters
A number of the parameters listed in Table 5-1 can be derived using other
variables (presented in Table 5-2) and a set of empirical, semi-empirical or
exact relationships. The derivation of parameters by the code is discussed
below.
5.5.1 Source-Specific Parameters
46
-------
The source-specific parameters which can be derived are the spread of the
source, the source width, the source length, and the source thickness (which
is actually classified as an aquifer parameter in the input). These are
discussed below in terms of the Gaussian and the Patch boundary conditions.
5.5.1.1 The Gaussian Source Boundary Condition
As described in Section 5.2, the gaussian source boundary condition is
uniquely specified when a (the standard deviation of the gaussian source), H
(the thickness of the source) and C0 (the maximum gaussian concentration) are
specified. Note that in the event of a finite duration source, an additional
parameter, Ts (the duration of the source), is necessary.
Source Thickness (Mixing Zone Depth) -- Percolation of water through the
facility (and unsaturated zone, if it exists) results in the development of a
mixing zone below the facility. This is shown in Figure 5.2. In the
saturated zone module, this mixing zone is the "source" of the contaminants,
hence the term source thickness is used interchangeably with the term mixing
zone depth. The thickness, H, of the mixing zone depends on the vertical
dispersivity of the media. If a value for H is not known, it can be derived
using the following relationship:
H = (20/,)172 + B(l - exp (- -L)) (5.49)
s
where
av = the vertical dispersivity [m]
L = the length scale of the facility--!.e., the dimension of the
facility parallel to the flow direction [m]
B = the thickness of the saturated zone [m]
In Equation 5.49, the first term represents the thickness of the mixing zone
due to vertical dispersion and the second term represents the thickness of the
mixing zone due to the vertical velocity below the facility resulting from
percolation. The derivation of the second term is presented in Appendix A.
While implementing this alternative, the code checks that the computed value
of the thickness of the source, H, is not greater than the thickness of the
aquifer, B.
The Spread of the Gaussian Source -- The standard deviation of the gaussian
source is a measure of the spread of the source. It can be derived by the
code using:
a = W/6
where
w = the width scale of the facility--!.e., the dimension of the
facility orthogonal to the groundwater flow direction
[m]Dividing by 6 implies that 99.86 percent of the area under
the gaussian source is flanked by the width of the facility.
The Length Scale and the width Scale--If the length, L, or the width, w, of
the source is not known, they can be derived by assuming that the waste
disposal facility has a square shape and taking the square root of the area:
47
-------
I 1/2
(5.51)
48
-------
5.5.1.2 The Patch Source Boundary Condition
The patch source boundary condition is uniquely specified when the width of
the source, 2y0, the source thickness, H, and the uniform concentration, C0,
are specified. The first two of these parameters can be derived.
TABLE 5-1. Primary Parameters Used in the Saturated Zone Transport Module
Parameters
Units
Source-Specific Parameters
Area of the land disposal facility
Leachate concentration at the waste facility
Either:
or
Standard deviation of the source (Gaussian)
width of source (Patch)
[m2]
[mg/1, g/m3]
[m]
[m]
Infiltration rate
Recharge rate into the plume
Duration of the pulse
Source decay constant
Aquifer-Specific Parameters
Porosity
Thickness of the aquifer
Thickness of source
Seepage velocity
Dispersivities (longitudinal, transverse, vertical)
Retardation coefficient
Radial distance from the site to the receptor
Angle between the plume center and the receptor
Well vertical distance
Time value at which concentration is required
Chemical-Specific Parameters
Effective first-order decay coefficient
Distribution coefficient
Overall Chemical Decay Coefficient
Solid phase decay coefficient
Dissolved phase decay coefficient
Bulk density
Distribution coefficient
Porosity
[m/yr]
[m/yr]
[yr]
[1/yr]
[cc/cc]
[m]
[m]
[m/yr]
[m]
[dimensionless]
[m]
[degrees]
[fraction]
[yr]
[1/yr]
[cc/g]
[1/yr]
[1/yr]
[g/cc]
[cc/g]
[cc/cc]
49
-------
TABLE 5-2. Parameters Used to Derive Other Saturated Zone Transport Module
Parameters
Parameters
Units
Solid and Dissolved Phase Decay Coefficients
Reference temperature
Aquifer temperature
Second-order acid-catalysis hydrolysis rate
constant at reference temperature
Second-order base-catalysis hydrolysis rate
constant at reference temperature
Neutral hydrolysis rate constant at reference
temperature
pH of the aquifer
Retardation Coefficient
Bulk density
Distribution (i.e.
Porosity
adsorption) coefficient
Bulk Density
Porosity
Porosity
Mean particle diameter of the porous medium
Particle Diameter
Porosity
Distribution Coefficient
Normalized distribution coefficient
for organic carbon, Koc
Fractional organic carbon content
Seepage Velocity
Hydraulic gradient
Hydraulic conductivity
Porosity
[f/mole-yr]
[f/mole-yr]
[1/yr]
[pH units]
[g/cc]
[cc/g]
[cc/cc]
[cc/cc]
[cm]
[cc/cc]
[cc/g]
[dimensionless]
[m/m]
[m/yr]
[cc/cc]
(continued...)
50
-------
TABLE 5-2. Parameters Used to Derive Other Saturated Zone Transport Module
Parameters (concluded)
Parameters Units
Hydraulic Conductivity
Porosity [cc/cc]
Mean particle diameter of the porous medium [cm]
Thickness of the Source (Mixing Zone Depth)
Length of the land disposal facility [m]
Thickness of the aquifer [m]
Seepage velocity [m/yr]
Porosity [cc/cc]
Infiltration rate through the facility [m/yr]
Vertical dispersivity [m]
Standard Deviation of the Source
width of the land disposal facility [m]
Length and Width of the Facility
Area of the land disposal facility [m2]
Dispersivities
Radial distance from the site to the receptor [m]
Source Thickness (Mixing Zone Depth)--If the source thickness is not specified
directly, it is derived in the same way as for the Gaussian source.
width of the Patch Source--The width of the patch source should be input as a
user-specified value. However, the width of the source can be derived as the
square root of the area using Equation 5.51.
5.5.2 Chemical-Specific Parameters
Two chemical-specific parameters shown in Table 5-1 can be derived from other
parameters: the overall chemical decay coefficient and the distribution
coefficient. Two additional chemical-specific parameters (shown in Table 5-2)
can be derived as well. They are the solid phase and the liquid phase decay
coefficients.
5.5.2.1 The Overall Chemical Decay Coefficient
The overall chemical decay coefficient for the saturated zone can be input
directly. When it is not known, it can be derived by the code using the first
term in Equation 5.3. If the solid-phase and the liquid-phase decay
coefficients are unknown, they too can be derived from temperature-corrected
values for hydrolysis rates, as described below. Note that chemical
degradation within the saturated zone is limited to hydrolysis and the by-
products of hydrolysis are assumed to be non-hazardous.
51
-------
The acid-catalyzed, neutral and base-catalyzed hydrolysis rates are all
influenced by groundwater temperature. This effect is quantified in the code
using the Arrhenius equation, that yields:
(5.52)
52
-------
where
K^b and Ka/b = the second-order acid- and base-catalysis hydrolysis rate
at temperature Tr and T respectively [f/mole-yr]
Kjf and K^ = the neutral hydrolysis rate at temperatures Tr and T
respectively [1/yr]
T = temperature of the groundwater [°C]
Tr = reference temperature [°C]
Rg = universal gas constant [1.987E-3 kcal/deg-mole]
Ea = Arrhenius activation energy [kcal/mole]
Note that, using the generic activation energy of 20 kcal/mole recommended by
Wolfe (1985), the factor Ea/Rg has a value of about 10,000.
The acid catalyzed, base catalyzed and neutral hydrolysis rate constants are
combined (Mill et al., 1981) to yield the composite, first order, dissolved
phase hydrolysis rate:
X± = fC/[ff+] + fC/ + fC/COff-] (5.53)
where
[H+] = the hydrogen ion concentration [mole/f]
[OH"1 = the hydroxyl ion concentration [mole/f]
Note that [H+] and [OH"1 are both computed from the pH of the aquifer, i.e.,
[ff+] = 10'PH (5.54)
[Off"] = io-(14-PH) (5.55)
For the case of sorbed phase hydrolysis, evidence suggests that base
neutralized hydrolysis can be neglected and that the acid neutralized
hydrolysis rate is enhanced by a factor of a. Thus, the effective sorbed
phase decay rate is expressed as:
A,, = afC/[ff+] + fC/ (5.56)
where a is the acid-catalysis hydrolysis rate enhancement factor for the
sorbed phase with a typical value of 10.0.
5.5.2.2 The Distribution Coefficient
The relationship most suited for relating the chemical distribution
coefficient, Kd, to soil or porous medium properties is discussed in detail by
Karickhoff (1984). In the absence of user-specified values, hydrophobic
binding is assumed to dominate the sorption process. For this case, the
distribution coefficient can be related directly to soil organic carbon
content using:
Ka = Kocfoc (5.57)
53
-------
where
Koc = normalized distribution coefficient for the chemical on organic
carbon [ml/g]
foc = organic carbon content in the saturated zone [dimensionless
fraction]
5.5.3 Aquifer-Specific Parameters
Of the 12 aquifer-specific input parameters shown in Table 5-1, six can be
either input directly or derived. In addition, some of the parameters used to
derive these primary parameters (see Table 5-2) can also be derived.
5.5.3.1 Retardation Coefficient
The retardation coefficient can be derived in the code using Equation 5.2.
5.5.3.2 Porosity and Mean Particle Diameter
In the absence of a user-supplied value for porosity, 6, it can be calculated
from the particle diameter using the following empirical relationship (Federal
Register, Vol. 51, No. 9, pp. 1649, 1986):
6 = 0.261 - 0.0385 ln(d) (5.58)
where
d = the mean particle diameter [cm].
Using the same relationship, the mean particle diameter can be derived from a
user-supplied value for porosity. Thus, only one or the other may be derived
in a given simulation.
5.5.3.3 Bulk Density
The soil bulk density directly influences the retardation of solutes and is
related to the soil structure. An exact relationship between the soil
porosity, particle density, and the bulk density can be derived (Freeze and
Cherry, 1979). This relationship is used in the code to derive bulk density.
Assuming the particle density to be 2.65 g/cc, the relationship is expressed
as:
pb = 2.65(1 - 8) (5.59)
where
pb = the bulk density of the soil [g/cc].
5.5.3.4 Seepage Velocity
The seepage velocity is related to aquifer properties through Darcy's Law.
Assuming a uniform, saturated porous medium, the seepage velocity can be
derived in the code using the following relationship:
54
-------
V =
KS
(5.60)
55
-------
where
K = the hydraulic conductivity of the formation [m/yr]
S = the hydraulic gradient [m/m]
6 = porosity [dimensionless]
Note that in general, the hydraulic gradient is a function of the local
topography, groundwater recharge volume and location, and the volume and
location of groundwater withdrawals. Further, it may also be related to the
porous media properties.
5.5.3.5 Hydraulic Conductivity
In the absence of site-specific measurements, a hydraulic conductivity value
can be derived using approximate functional relationships. One such
relationship, the Kozeny-Carman equation (Bear, 1979), is included in the
model:
K = £2. -? (5.61)
U (1-6)2 1-8
where
K = the hydraulic conductivity [cm/s]
p = the density of water [kg/m3]
g = acceleration due to gravity [m/s2]
u = the dynamic viscosity of water [N-s/m2]
d = mean particle diameter [cm]
In Equation 5.61 the constant 1.8 includes a unit conversion factor. Both the
density of water, p, and the dynamic viscosity of water, u, are functions of
temperature and are computed using regression equations presented in CRC
(1981). Note that at 15°C, the value of [pg/1.8u] is about 478. After using
Equation 5.61 to derive a value for hydraulic conductivity, the code converts
the value to units of meters per year.
The only use of the hydraulic conductivity parameter in MULTIMED is to compute
the seepage velocity. Therefore, if the user chooses to input the seepage
velocity directly, instead of specifying it as a derived parameter, the
hydraulic conductivity value is not needed.
5.5.3.6 Dispersion Coefficients
The model computes the longitudinal, lateral and vertical dispersion
coefficients as the product of the seepage velocity and longitudinal (aL) ,
transverse («T), and vertical (av) dispersivities. In the absence of user
specified values for dispersivities, the model allows two alternatives.
Alternative one is based on the values presented in Gelhar and Axness (1981).
Dispersivities are calculated as a fraction of the distance to the
downgradient receptor well, as follows:
0L = 0.1 xr (5.62)
56
-------
3.0
(5.63;
av = 0 . 056ar
(5.64)
where xr = the distance to the receptor well [m]. Option one is summarized in
Table 5-3 (a) .
Alternative two allows a probabilistic formulation for the longitudinal
dispersivity as shown in Tables 5-3 (a) and 5-3(b) [Gelhar (personal
communication), 1986]. The longitudinal dispersivity is assumed to be uniform
within each of the three intervals shown in Table 5-3(b). Note that the
values of longitudinal dispersivity shown are based on a receptor well
distance of 152.4 m. For other distances, the following equation is used:
QL(xr) =
= 152.4) (xr/152.4)'
(5.65)
TABLE 5-3(a). Alternatives for Including Dispersivities in the Saturated
Zone Flow Module
Dispersivity
«L [m]
«T [m]
av [m]
aL/aT
aL/av
TABLE 5-3 (b) .
Class
«L (m)
Probability
Cumulative
Alternative 1
Existing Values
0.1 xr
0.333 «L
0.056 aL
3
approx. 18
Probabilistic Representation
a Distance
1 2
0.1-1 1-10
0.1 0.6
0.1 0.7
Alternative 2
Gelhar 's Recommendation
Probabilistic Formulation
(See Table 5-3 (b) )
QijQ
aL/160
8
160
of Longitudinal Dispersivity for
of 152.4 m
3
10-100
0.3
1.0
Probability
The transverse and vertical dispersivity are assumed to have the following
values:
aT = aL/8
(5.66)
av = aL/160
(5.67)
57
-------
5.5.3.7 Source Thickness (Mixing Zone Depth)
The derivation of this aquifer-specific parameter is discussed in Section
5.5.1.1.
SECTION 6
The Surface Water Module
6.1 Introduction
This chapter presents details about the fate and transport of a toxic chemical
pollutant in a surface stream as included in MULTIMED. As illustrated in
Figure 6.1(a and b) , the module is based on the assumption that contamination
of the stream occurs due to the complete interception of a steady-state
groundwater plume. Interception of a finite pulse source can not be
simulated. Also, the case of partial penetration of the plume is not
considered. Further, the case of stream contamination due to direct overland
runoff from a 'failed' waste containment facility is not discussed. This
route of exposure has been analyzed and discussed in detail by Ambrose et al.
(1985) .
The technical approach presented in the following sections is based on the
surface water exposure assessment model developed by Ambrose et al. (1985),
and an enhanced version of that surface water model, termed SARAH2
(Vandergrift and Ambrose, 1988).
6.2 Mathematical Description of the Surface Water Module
The surface water model assumes complete interception of a contaminant plume
by a stream perpendicular to the stream as shown in Figures 6-1(a) and 6-1(b).
Further it is assumed that at x=0 (the downstream edge of the near-field
mixing region), the river is laterally as well as vertically mixed. With
these assumptions, the following mass balance equation can be written:
Mg = Cs (Qg + Qfl) (6.1)
where
Mg = contaminant mass flux entering the stream, [g/yr]
Cs = the near-field fully mixed contaminant concentration in the
stream [g/m3 or mg/f]
Qg = groundwater discharge [m3/yr]
Qs = stream discharge [m3/yr]
The stream discharge is supplied by the user. The groundwater discharge, Qg,
is calculated as:
Qg = V'Q (H + 2v/a^T) (3.920 + 2^0^T) (6.2)
58
-------
where
V* = the steady-state retarded seepage velocity [m/yr]
6 = the porosity [cc/cc]
H = the depth of penetration of the source [m]
av, aT = the vertical and transverse dispersivity [m]
xs = the distance of stream from the downgradient edge of the waste
disposal unit [m]
a = the standard deviation of the source [m]
Equation 6.2 calculates the volumetric flux from the groundwater, Qg, using a
steady-state retarded seepage velocity and an estimate of the cross-sectional
area of the plume that intercepts the stream. The cross-sectional area is
estimated based on the vertical and lateral spread of the plume under uniform
flow conditions. Note that in Equation 6.2, the factor 3.92 accounts for 95%
of the area of the initial gaussian source. In the event that the patch
source boundary for the saturated zone model is used, the term 3.92a in
Equation 6.2 is replaced by w, the width of the patch source.
For a continuous, steady source of contaminant at the waste disposal facility,
the mass loading to the stream, Mg, for steady-state conditions is:
where
v
Mg = ML exp (- -JLi) exp (- -S-J.) (6.3)
Xvd X x
JLi) exp (- -S-J.
ML = mass flux of contaminant leaching from the land disposal
facility [kg/yr] (also see Equation 5.34)
X, = the overall decay coefficient in the unsaturated zone [1/yr]
dv = the distance traveled by the contaminant within the unsaturated
zone [m]
V* = the retarded seepage velocity in the unsaturated zone [m/yr]
Xs, xs, V* = the corresponding quantities for the saturated zone
In the event that the unsaturated zone is simulated as a layered system
consisting of n layers, Equation 6.3 can be modified as:
^ X d X x
Mg = ML exp [- ^ (~^~}i ] 8XP [~ ~~^~] (6-4)
V S
Combining equations 6.1 and 6.3 or 6.4 gives the concentration of the
contaminant in the fully mixed portion of the stream at the downstream edge of
the near field mixing zone. At a downstream location within the stream, the
contaminant concentration [mg/f] is reduced due to degradation. Assuming a
first order degradation rate, the downstream concentration is:
CR(xR) = Csexp (i) (6.5)
59
-------
where
CR(xR) = contaminant concentration in the stream [mg/f]
AR = the in-stream decay rate [1/yr]
XR = the distance from the downstream edge of the mixing zone to the
point of interest [m]
VR = the mean in-stream velocity of flow [m/yr]
The model estimates the in-stream decay rate, XR, and the mean in-stream
velocity of flow, VR, by methods described below.
6.2.1 Computation of the In-Stream Decay Coefficient
The in-stream fate processes of volatilization and hydrolysis are included in
the model and assumed to be first order reactions. The overall decay
coefficient is represented as:
where
KHR = the hydrolysis rate constant at the in-stream temperature, T
[1/yr]
KVR = the volatilization rate constant [1/yr]
The factor 3 . 15E7 changes the units of XR from [1/s] to [1/yr] . The nominal
hydrolysis rate constant is calculated from the acid-catalyzed, neutral and
base-catalyzed hydrolysis rate constants (Burns et al . , 1982; Mill et al . ,
1981) :
where
K^r and Kjf = the second-order acid-catalysis and base-catalysis
hydrolysis rate constant at the reference temperature, Tr
[I/mole-yr]
Kjf = the neutral hydrolysis rate at temperature, Tr [1/yr]
[H+] = hydrogen ion concentration
a = sorption catalysis hydrolysis rate enhancement factor for the
sorbed compound [dimensionless]
Kn = neutral hydrolysis rate constant at the reference temperature
[1/yr]
fs = fraction of the chemical that is sorbed [dimensionless]
fD = fraction of the chemical that is dissolved [dimensionless]
[OH"] = the hydroxyl ion concentration [mole/f]
Note that [H+] and [OH"] can be computed from the pH of the stream.
The hydrolysis rate constants can be expressed as a function of temperature by
using the Arrhenius equation:
K^ = A exp [ -EjRgT] (6.8)
60
-------
where
Ea = the Arrhenius activation energy [kcal/mol]
Rg = the gas constant [1.987E-3 kcal/deg-mol]
T = the temperature [K]
The pre-exponential factor, A, has the same units as the hydrolysis rate
constant. Using the above equation, the hydrolysis rate constant is corrected
to ambient stream temperature by using the following expression:
(6.9)
where TE and T are the reference temperature and in-stream water temperature in
°C. Note for the case when the activation energy is 20 kcal/mole, the value
of Ea/Rg is 10,000. This value is used in the model.
The second transformation pathway considered is volatilization. The
volatilization rate constant is calculated from the Whitman, or two-resistance
model (Whitman, 1923; Burns et al . , 1982):
R + RG
s L G
where
KVR = the volatilization rate constant [1/s]
Ds = mean stream depth [m]
RL = liquid phase resistance [s/m]
RG = gas phase resistance [s/m]
The second term in Equation 6.10 represents the conductivity of the compound
through a liquid and a gas boundary layer at the water surface. The liquid
phase resistance to the compound is assumed to be proportional to the transfer
rate of oxygen, which is limited by the liquid phase only:
(6.11)
where
K02 = reaeration rate constant [1/s]
MW = molecular weight of the compound [g/mole]
32 = molecular weight of oxygen [g/mole]
The gas phase resistance to the compound is assumed to be proportional to the
transfer rate of water vapor, which is limited by the gas phase only:
G (H'/RT) (6.12)
61
-------
where
WAT = vapor exchange constant [m/sec]
HT = Henry's law constant [atm-m3/mole]
R = ideal gas constant [8.206 x 1(T5 m3-atm/mol°K]
T = stream water temperature [°K]
18 = the molecular weight of water [g/mole]
The reaeration and water vapor exchange constants will vary with stream reach
and time of year. They can be calculated using one of several empirical
formulations. In this model the reaeration rate constant is calculated by the
Covar method using stream velocity, VR, and depth, Ds, corrected for
temperature (Covar, 1976). The water vapor exchange constant is calculated
using wind speed and a regression relation proposed by Liss (1973):
WAT = 5.16 x 10"5 + 3 .156 x 10"3 W
where W = mean wind speed at 10 cm above surface [m/sec].
(6.13)
wind speed measured above 10 cm is adjusted to the 10 cm height assuming a
logarithmic velocity profile and a roughness height of 1 mm (Israelsen and
Hansen, 1962):
W= W log (0.1/0.001)/log(z/0.001)
(6.14)
where
wz = wind speed at height z [m/sec]
z = wind measurement height [m]
6.2.2 Calculation of the Mean In-Stream Velocity
The mean in-stream velocity, VE, is estimated in MULTIMED as:
VR = (Ds)2/3(SLOPE/n)1/2
(6.15)
where
Ds = mean stream depth [m]
SLOPE = channel slope [dimensionless]
n = Mannings roughness coefficient
[dimensionless]
6.3 Exposure Due to Surface Water Contamination
6.3.1 Routes of Exposure
Three routes of exposure to toxic substances in surface streams are considered
in the model. These are human exposure due to drinking contaminated water,
human exposure due to the consumption of fish exposed to the contaminated
water, and exposure to aquatic organisms. Figures 6.2(a), 6.2(b) and 6.2(c)
show flow charts of the various stages between failure of the waste
containment facility and these exposure routes. MULTIMED uses the
concentrations, CR and Cs, to predict the concentration of a contaminant in
62
-------
drinking water, fish, and/or aquatic organisms. Based on these results and
using procedures discussed in Section 9.6, the user can manually back-
calculate the maximum allowable concentration of the contaminant in the waste
disposal facility that is protective of human health and toxic aquatic
effects .
6.3.2 Human Exposure to Toxics through Drinking Water
Human exposure to dissolved chemicals in drinking water is assumed to occur
through water obtained from a treatment plant located downstream from the
initial mixing zone. The water is assumed to be treated by a primary settling
process allowing suspended solids and adsorbed chemicals to settle. The
influent concentration, CR is reduced to a concentration of CDB in the effluent
drinking water from the water treatment plant. This concentration is
represented by:
CDW = fD CR (6.16)
where fD is equal to the fraction of the contaminant that is dissolved. It can
be expressed as:
fn = (6.17)
where Cad = sorbed concentration [mg/f] . The dissolved aqueous concentration,
CR, and the sorbed concentration, Cad, can be related by:
(6.18)
where
S = sediment concentration [mg/f]
Kp = the equilibrium distribution coefficient [f/mg]
Further, it can be shown that for the sorption of hydrophobic organic
compounds :
(6.19)
where
Koc = the organic carbon partition coefficient [I/kg]
foc = organic carbon content of suspended solids [fraction]
Substituting Equations 6.18 and 6.19 into Equation 6.18, fD is expressed as:
f° = KocfolcS + 1 (6'2°)
63
-------
and the fraction of sorbed contaminant is (1 - fD) .
6.3.3 Human Exposure to Toxics Due to Fish Consumption
Dissolved neutral organic compounds in the water can be taken up by fish
through exchange across the gill and gut membranes, and through the skin.
Contaminated food can be ingested, resulting in further exchange of compounds
across the gut membrane. Concentration levels in the fish rise until the
activity of a compound in the blood equals the activity of that compound in
the water. This condition represents chemical equilibrium. Further uptake of
the compound, which results in higher blood concentrations, will lead to net
exchange out of the fish through the gill, gut, kidney and skin. The whole-
fish concentration of the pollutant can be expressed as:
CF = KFfDCB (6.21)
where KF equals the entire fish partition coefficient or the bioconcentration
factor. Note that the use of the near-field, in-stream concentration, Cs, in
Equation 6.21 rather than the downstream, attenuated concentration, CR, assumes
that the fish reside continuously within the most polluted reach of the stream
and, hence, is a conservative assumption.
If the fish is exposed to a steady aqueous concentration over a long period of
time, the distribution of the compound within the lipid and non-lipid tissues
of the fish will equilibrate so that:
Cl = KlCB (6.22)
and
Cnl = Km CB (6.23)
where
C, = the lipid phase biomass concentration [mg/kg]
CB = the concentration in the blood[mg/f]
K, = the lipid phase partition coefficient [I/kg]
Cn, = the non-lipid (blood, muscle) phase biomass concentration
[mg/kg]
Kn, = the non-lipid phase partition coefficient [I/kg]
Neglecting bioaccumulation, the equilibrium concentration in the blood will
not exceed the dissolved concentration in the river, i.e.,
CB = fDCs (6.24)
The average whole fish concentration, CF [mg/kg], can be expressed as the
weighted sum of the tissue concentrations:
CF = f, C, + (1-f,) Crf (6.25)
64
-------
where f, = lipid fraction of the fish biomass [dimensionless] .
Combining Equations 6.22 to 6.25 with Equation 6.21, CF can be expressed as:
Cp= [f,«i + (1-f,) Knl] fDCs (6.26)
Comparing Equation 6.21 and Equation 6.26, the entire fish partition
coefficient can be expressed as:
KF = [ffi + (1-f,) Knt] (6.27)
For less hydrophobic compounds, Kn, may contribute significantly to KF. Non-
lipid tissue is composed primarily of water along with protein and
carbohydrates. By replacing Kn, with the octanol-water partition coefficient
and assuming that partitioning to non-lipids is less than or equal to 1% of
the partitioning to lipids, the model calculates a conservative estimate of KF
as :
KF = KOW (ft + 0.01) (6.28)
The octanol water partition coefficient, Kow, is approximated using the
correlation with Koc developed by Karickhoff et al . (1979):
*oc=°-41^ (6.29)
6.3.4 Exposure of Aquatic Organisms to Toxics
Aquatic organisms are exposed to contaminants present in the stream. Only
dissolved species of a compound cross the membranes of aquatic organisms and
cause internal exposure. There is some evidence, however, that the presence
of suspended contaminant solids can enhance the rate of uptake and the
internal exposure to a compound. Also, the Criterion Continuous Concentration
(CCC) , which is set to protect against toxic effects to aquatic organisms, is
usually referenced to the total concentration of the contaminant in the
stream. Therefore, exposure to aquatic organisms is estimated by using the
total concentration of the compound in the stream.
6.4 Assumptions and Limitations
Following is a list of important assumptions on which the surface water module
is based:
(a) The surface stream model considers the case of complete plume
interception only.
(b) The module considers only the case of a steady-state continuous
source at the landfill.
(c) The stream is laterally and vertically well mixed. This implies
that the stream is relatively small.
65
-------
(d) The only in-stream transport process considered is lumped first
order decay due to volatilization and hydrolysis.
(e) The municipal water treatment plant assumes removal of contaminants
by adsorption onto particles that are removed by sedimentation or
filtration processes.
(f) The concentration of the contaminant in fish is assumed to be at
equilibrium with the near-field, in-stream concentration.
6.5 Data Requirements
Table 6-1 lists the surface water and chemical input parameters required to
compute the in-stream concentration, as well as concentrations in drinking
water, fish, or aquatic organisms. The surface water module must be run in
conjunction with the saturated transport module (and the unsaturated zone
modules, if needed). The saturated zone parameters are listed in Tables 5-1
and 5-2. Note that the location of the stream is specified using saturated
zone parameters (see Section 5.2.3).
TABLE 6-1. Parameters Required for the Surface Water Module
Parameters Units
Surface Stream Specific Data
Stream discharge
Distance to drinking water plant intake
Mean stream depth
Mannings roughness coefficient
Channel slope
Sediment concentration
Organic carbon fraction of suspended solids
pH of the stream
In-stream temperature
Fraction of fish that is lipid
Wind speed
Height at which wind speed is measured
Chemical Specific Data
Reference temperature
Second-order acid catalysis hydrolysis rate
constant at reference temperature
Second-order base catalysis hydrolysis rate
constant at reference temperature
Neutral hydrolysis rate constant at reference
temperature
Henry's law constant
Molecular weight of the contaminant
Normalized partition coefficient (i.e., Koc)
[m3/sec]
[m]
[m]
[dimensionless]
[dimensionless]
[mg/f]
[dimensionless]
[pH units]
[°C]
[dimensionless]
[m/s]
[m]
[f/mole-yr]
[f/mole/yr]
[1/yr]
[atm-m3/rnole]
[g/mole]
[ml/q]
For Monte Carlo simulations, the parameters shown in Table 6-1 may be input as
constant values or as distributions. The values of some of these parameters
can be computed indirectly using the parameter estimation methods discussed in
the MULTIMED user's manual.
66
-------
SECTION 7
The Air Emissions Module
7.1 Introduction
This chapter describes the algorithm used to estimate the air emission of
toxic substances from land-based waste disposal facilities. The current
version of the air emissions module includes an algorithm to simulate the
emissions from Subtitle C facilities. Algorithms to estimate emissions from
impoundments and Subtitle D facilities are currently not included in
MULTIMED's air emissions module.
The model accounts for diffusive transport through porous media based on
Pick's law. The effect of atmospheric pressure fluctuations and other
transport processes is accounted for by empirical factors. Approximate values
of these parameters, based on available numerical and analytical studies, are
included.
The current version of the air emissions module assumes that the wastes
contained in the facility are suitably segregated into cells so that there is
no chemical or biochemical activity within the facility. The case of co-
disposal with municipal or liquid wastes is not considered here. This
suggests that the possibility of gas generation within the landfill is highly
unlikely. Further, the waste is considered to be covered with soil. Thus,
the complex process of emissions from uncovered wastes is not included in the
model.
The model's air emissions estimates are independent of predictions of
unsaturated/saturated transport. Therefore, in the current version of the
code, the air emission and dispersion simulations are run separately from
subsurface transport simulations.
This chapter is restricted to the air-emissions module only and does not
include any discussion of near field or far-field dispersion or transport
processes. The latter is discussed in Section 8.
7.2 Governing Equations
7.2.1 The Air Emissions Diffusion Model
The model described below was first developed by Farmer et. al. (1978) for
computing the vapor flux of hexachlorobenzene through a dry soil cover from a
landfill. This method treats the waste volatilization or vapor loss of the
chemical from a landfill as a diffusion controlled process using Pick's Law
for steady-state diffusion. The diffusion into the atmosphere is assumed to
emanate from a plane surface with a constant concentration. The emission rate
is expressed as:
E. = De± (CB± ~ Cai)TSEC (7.1)
67
-------
where
E-L = the emission rate of chemical i [g/yr]
L = mean depth of the soil cover [cm]
Del = effective diffusion coefficient for chemical i in soil [cm2/s]
Csl = pore space concentration of chemical i [g/cc]
Cal = concentration of the chemical i in air [g/cc]
Aw = area of the land disposal facility [cm2]
TSEC = the number of seconds in a year [s]
7.2.2 The Effect of Atmospheric Pressure Fluctuations
The model, as presented above, accounts for only the diffusive transport of
vapors through the porous media. The effect of atmospheric pressure
fluctuations that 'pump' chemical vapors from the waste disposal facility is
not considered. Typical ranges of atmospheric pressure fluctuations are
presented in Table 7-1.
The occurrence of large barometric changes is associated with the passage of
weather fronts, typically 4 to 8 days in duration. Clements and Wilkening
(1974) investigated the effect of large scale atmospheric changes in 222Rn flux
across a soil-air interface. Field data show that pressure changes of 1 to 2
percent associated with the passage of frontal systems produce changes in the
222Rn flux from 20 to 60 percent, depending on the rate of change of pressure
and its duration. This finding was confirmed in a laboratory experiment using
a vertical column of 226Ra-bearing sand. Lu and Matuszek (1978) observed
similar atmospheric pressure-induced gas flow of tritiated compounds from
waste in a commercial radioactive-waste land burial site.
TABLE 7-1.
Parameter
Pressure
Amplitude
(mm Hg)
General Characteristics
Frontal Passage
10-20
of Atmospheric Pressure
Diurnal Variation
1-3
Fluctuations*
Local Gustiness
0.1-0.2
Duration 4-8 (d) 24 (h) 10-30 (s)
"Springer et al. (1986)
The pumping effect of diurnal barometric variations in extracting gases and
vapors from subsurface porous formations apparently has not been investigated
in detail. Fukuda (1955) investigated air and vapor movement in soil due to
wind gustiness. The soil depth to which air can penetrate as a result of wind
gustiness is very small. It was found that in sandy soil with mean particle
diameters of 0.25 to 0.5 mm, air penetrates only about 5 mm below the surface.
Models that simulate the effect of pressure fluctuations (Thibodeaux, 1982)
have not been sufficiently calibrated and verified. However, this effect has
been incorporated into MULTIMED by using an empirical factor, eBP. Thibodeaux
(1982) computed a value of 1.13 for this factor using a numerical simulation
model and two weeks of atmospheric pressure data.
7.2.3. Other Transport Processes
Thibodeaux et al. (1986) compared the flux estimated using Equation 7.1 with
experimental data (without atmospheric pressure fluctuations). Their results
indicate that Equation 7.1 underestimates fluxes by a factor of approximately
three . This anomaly between the measured and computed values has been
ascribed to the process of surface diffusion (Cussler, 1983) that occurs on
68
-------
the walls of soil grains. This is a relatively fast process and could
account, in part, for the enhanced fluxes.
In the absence of detailed physical understanding of this process, it is best
to incorporate its effect as an empirical enhancement factor, e£.
Incorporating the two empirical factors, e£ and eBP, into Equation 7.1, the
model equation can be expressed as:
E. = De± (Csi - Cai) eBp ef TSEC (7.2)
The concentration of the chemical in air, Cal, is assumed to be zero in the
code. This is valid for the case of high wind velocities at the surface that
rapidly transport vapor emissions away from the land disposal unit.
7.2.4 Computation of the Diffusion Coefficient
The effective diffusion coefficient for a chemical in soil, Del, is computed in
the code using a general, empirical relationship proposed by Currie (1961).
The relationship, which is based on a hydrogen diffusion experiment, accounts
for the effect of moisture content:
^ei = YDoiOr (9a)a (7.3)
where
D01 = diffusion coefficient for chemical i in air [cm2/s]
6a = air-filled porosity of soil [cc/cc]
6t = total porosity of soil [cc/cc]
a = exponent which is 4.0 for granular material subject to moisture
testing
Y = factor which varies with soil type (0.8 to 1.0)
u = exponent which varies with soil type (1.4 for spherical grains
(sand), 2.6 for kaolin (clay) and 11.0 for plate minerals)
The effect of temperature on the diffusion coefficient is considered in the
code using:
DOT = DOR (- J1'5 (7.4)
where DOR and DOT = the diffusion coefficients at the reference temperature, TE,
and temperature of the waste disposal unit, T, respectively.
7.2.4.1 Effect of Engineering Controls
The current version of MULTIMED can not account for the effect of engineering
controls. If the effect of composite engineering controls, such as vegetative
soil, a synthetic membrane, or a drain layer, is important to results,
Equation 7.2 can be solved by manual calculation using the resistances-in-
series concept. The total effective diffusion for the cover can be computed
as :
69
-------
n = ^2i
V^ , Ln , (7.5)
where
Ln = the thickness of each layer [cm]
Denl = the effective diffusion coefficient for chemical i in layer n
[cm2/sec]
N = total number of layers
Ltotai = sum of thickness of the N layers [cm]
7.2.5 Computation of Vapor Concentration, Csl
The accuracy of the emission estimates obtained using the above model
significantly depends on the accuracy with which the vapor concentration in
the pore spaces of the waste disposal unit is computed. This concentration
depends on the properties of the chemical species, the chemical quantity in
the cell, and the state of the chemical with respect to the other materials in
the cell. Table 7-2 (Thibodeaux, 1982) shows a categorization of the various
states of chemical substances in landfills and the equilibrium laws that can
be used to determine the pore space chemical concentration. Further details
about these methods is available in Groves et al. (1984).
Application of these eight different theories to compute the pore space
concentration requires detailed information that is typically not available
for a waste facility. MULTIMED air emission module assumes that the
information is not available and therefore adopts a conservative approach.
Maximum vapor phase chemical concentrations of solid wastes in pure form or as
mixtures of solid flakes and granules are considered. For such conditions, it
is safe to assume that a chemical will exert pure component vapor pressure.
Thus, the saturation vapor concentration for a chemical in the waste can be
determined using the ideal gas law:
X.P M.
where
Mi = mole weight of chemical i [gm/mole]
R = molar gas constant (mm Hg-cm3/°K-mole)
T = temperature of the landfill [°K]
X-L = mole fraction of chemical i in the mixture [dimensionless]
P01 = vapor pressure of chemical i [mm of mercury]
7.3 Assumptions and Limitations
Following is the list of assumptions on which the air emissions algorithm is
based:
(a) The air emissions algorithm is not applicable to the case where
landfill gas is generated within the facility.
70
-------
(b) Far field atmospheric dispersion of the contaminants is treated
separately as discussed in Section 8.
(c) Enhancement in the rate of emissions due to barometric pumping and
other, as yet undefined processes, is included using empirical
factors.
(d) The vapor phase contaminant concentration within the land disposal
unit is computed assuming that the unit contains chemicals in pure
solid state. This is a conservative approach in that it tends to
overpredict the emission of contaminants.
7.4 Data Requirements
The data required for the air emissions model are shown in Table 7-3.
TABLE 7-3. Parameters Required for the Air Emissions Module
Parameters
Units
Air Emission-Specific Parameters
Effective depth of the soil cover
Water content of the soil
Total porosity of the soil
Empirical coefficients to compute diffusion
coefficient in soil from a known diffusion
coefficient in air
Flux enhancement factor due to barometric
pressure fluctuations
Enhancement factor due to other transport
processes
Temperature of the landfill
Chemical-Specific Parameters
Mole weight of the chemical
Mole fraction of chemical in landfill
Vapor pressure of chemical
Diffusion coefficient for chemical in air
at a reference temperature
Reference temperature for air diffusion
Source-Specific Parameters
Area of the land disposal facility (converted
to cm2 in program)
[cm]
[cc/cc]
[cc/cc]
y, u, a
[dimensionless]
[dimensionless]
[dimensionless]
PC]
[g/mole]
[dimensionless]
[mm of mercury]
[cm2/s]
[m2]
71
-------
SECTION 8
Air Dispersion Module
8.1 Introduction
This section describes the air dispersion module used to calculate the
atmospheric transport of vapor emissions from a Subtitle C landfill. It must
be run in conjunction with the air emission model described in Section 7.
Pollutants emitted from landfills will be transported by advection due to the
mean wind field and by dispersion due to vertical and horizontal turbulent
wind fluctuations. During development of MULTIMED, over 60 air dispersion
models that simulate these transport processes were reviewed and classified
into the following five categories:
(i) Grid models
(ii) Trajectory models
(iii) Gaussian models
(iv) Screening level models
(v) Miscellaneous models
As will be discussed in subsequent sections, the model selected for inclusion
is a long-term gaussian plume model similar to the VALLEY model (U.S. EPA,
1977). This type of model has been verified through extensive use, accounts
for long-term variations in meteorological conditions, and involves a level of
computational effort appropriate for MULTIMED. In the current version of the
model, the air dispersion module has received only limited testing.
8.2 Governing Equations
8.2.1 The Gaussian Dispersion Equation
The general gaussian equation for the ground-level concentration of a
pollutant at a receptor located at a distance, x, from a point source (Wark
and Warner, 1981) is:
- <>''
where
C = ground level pollutant concentration at distance, x, from the
source [mg/t]
Q = source emission rate [g/s]
R = term that accounts for vertical plume dispersion and
reflection, computed using Equation 8.2
D = term that accounts for chemical decay, computed using
Equation 8.3
U = wind speed [m/s]
oy = horizontal (transverse) dispersion coefficient [m]
oz = vertical dispersion coefficient [m]
72
-------
y = transverse distance from the plume centerline [m]
This equation is derived by solving the one-dimensional, advection-dispersion
equation and assuming a gaussian distribution of pollutant in the transverse
and vertical directions. Note that the horizontal and vertical dispersion
coefficients are functions of atmospheric stability and the distance between
the source and the receptor.
The term, R, that accounts for vertical dispersion and the reflection of the
plume at the ground surface and at the top of the atmospheric mixing layer, is
computed using:
£! 2nL - H ! 2nL + H
{ exp [-( )2] + exp [-( )2]} (8.2)
where
Ha = height of the plume centerline above ground [m]
L = thickness of the atmospheric mixing layer [m]
Further, in Equation 8.1, D is an exponential decay term used to account for
the transformation and chemical degradation of pollutants:
D = exp (-k -|) (8.3)
where
k = a first-order decay rate constant (1/s)
Equation 8.1 assumes constant meteorological (wind speed, direction,
stability) and source emission conditions in level terrain, and is valid only
for calculating steady-state ground-level concentrations. To estimate the
long-term (i.e., seasonal, annual) average concentrations for periods during
which meteorological conditions change, a frequency-weighting approach, which
calculates concentrations for all possible combinations of wind direction,
wind speed, and stability class, is used.
The distribution of meteorological conditions over time is described by a
joint frequency distribution of wind direction, wind speed, and stability.
The range of possible wind directions in an area is represented by dividing
the areal x-y plane into 16 equal wind direction sectors each subtending an
angle of 2n/16 radians at the location of the point source (see Figure 8.1).
Wind speeds are similarly described by dividing the range of possible wind
speeds into six equal intervals, with each interval represented by a single
average wind speed. The six average wind speeds commonly used for reporting
meteorological data are supplied as defaults in the MULTIMED preprocessor.
They can be modified by the user, if desired.
Atmospheric stability is a measure of the turbulence in the atmosphere. The
two types of turbulence present in the environment are convective and
mechanical. Convective turbulence describes the turbulence caused by
differences between temperature at different heights (i.e., heat flux).
Mechanical turbulence describes the mixing caused by wind shear. Atmospheric
stability is represented by the Pasquill-Gifford (PG) classification scheme,
composed of the six discrete categories A through F defined in Table 8-1.
Stability Class A, the most unstable class, is dominated by convective
73
-------
activities typical of a sunny day with light winds. Stability Class D, or
neutral stability, describes a well-mixed atmosphere and is dominated by
mechanical turbulence. Stability Class F, the most stable class, describes a
stratified atmosphere where temperature increases with height and winds are
light. This is typical of a clear calm winter evening. Pollutant releases
disperse the most under unstable conditions and the least under stable
conditions.
To account for the variability of wind direction, pollutant concentrations are
averaged in the y-direction across each wind direction sector:
_ 16 QRD
c-' - = (8.4)
where
Ci = average ground level pollutant concentration in sector
xc = distance from the source to the receptor along the plume
centerline [m]
Note that the transverse dispersion terms in Equation 8.1 have been replaced
by the inverse of the sector width (2nxc /16) .
The long-term average concentration is then computed by weighting
concentrations for each combination of wind speed, wind direction, and
stability class by the corresponding joint frequency of occurrence. Thus the
long-term average concentration from a point source is expressed as (U.S. EPA,
1977) :
16 ^ f- f- r QRS±D
~^^ Z^ Z^ Z^ ^ijk L z__
Dj.azxcV/2n
fi^[7-^=] (8.5)
where
CL = long-term average ground-level pollutant concentration [mg/f]
fljk. = the joint frequency of occurrence of the ith wind direction,
jth wind speed, and kth stability class [dimensionless]
Si = a smoothing function computed using Equation 8.6
[dimensionless]
74
-------
TABLE 8-1. Atmospheric Stability Categories
Category Description
A Extremely unstable
B Moderately unstable
C Slightly unstable
D Neutral
E Slightly stable
F Moderately stable
In Equation 8.5, Si is a smoothing function used to eliminate discontinuities
in concentrations calculated near the boundaries of the wind direction
sectors. Frequencies for a wind direction sector may be very different from
those in the two adjacent sectors, resulting in unrealistic discontinuities at
the sector boundaries. To avoid this effect, the receptor concentrations are
calculated as weighted averages of concentrations in the wind direction sector
in which the receptor is located and the two adjacent sectors. The weighting
(or smoothing) factor is the normalized distance of the receptor from sector
centerlines as shown in Figure 8.1. Thus:
(8.6)
where
w = the distance between sector centerlines at the receptor [m]
y± = the distance of the receptor from the centerline of sector
i [m]
Equation 8.5 is the key equation used by the air dispersion module to
calculate ground-level concentrations at various distances from the waste
disposal facility. The following sections describe the application of the
gaussian equation to area sources, the estimation of plume rise, methods for
calculating vertical dispersion coefficients, and terrain corrections.
8.2.2 Area Source Approximations
Area sources are used to represent emissions from landfills and from areas
where point sources are too numerous to simulate individually. Emissions from
these sources are expressed as mass fluxes per unit area. Concentrations
resulting from area sources are calculated from the gaussian point source
dispersion equation using the virtual point source approximation. The
location of the virtual point source is determined such that the horizontal
plume dispersion at the centroid of the area source is equal to the width of
the area source. For long-term sector averaging models, the virtual point
source is located such that the distance between the wind direction sectors at
the centroid of the area source is equal to the width of the area source
(Figure 8.2).
75
-------
Thus, the distance, xc, used in Equation 8.5 is the distance from the virtual
point source to the receptor along the sector centerline. However, the
distance used in calculating the decay term (Equation 8.3) is the distance
from the centroid of the area source to the receptor (i.e., the actual
distance traveled by the pollutant).
The air dispersion module has an additional correction to the virtual source
approximation to account for receptors very close to the source. Receptors
near the source will not be affected by the entire area source, since
contributions from portions of the area will be carried past the receptor. In
these cases the effective source area is reduced as shown in Figure 8.3. The
emission rate used in Equation 8.5 is then the product of the effective source
area and the landfill emission rate per unit area.
8.2.3 Plume Rise
Gaussian plume models calculate dispersion after the plume has stabilized to
its initial height of rise. During emission the plume may rise vertically due
to initial buoyancy and momentum. However, emissions from landfills and other
area sources will often have negligible buoyancy or momentum. These emissions
occur by diffusion through the landfill cap material and will in general have
a temperature close to the ambient temperature. Area sources are therefore
usually modeled by assuming that the plume centerline does not rise above the
elevation of the source during emission. An option is available, however, to
estimate initial plume rise by assuming that the height of rise from area
sources is linearly proportional to wind speed (Irwin et al., 1985):
h = ( ) Ah5 (8.7)
where
h = height of plume rise [m]
Ah5 = height of rise at U = 5 m/s [m]
In the above relationship, Ah5 is a model input parameter which must be
determined by the user from field studies or other estimation methods.
8.2.4 Estimation of Vertical Dispersion Coefficient
The vertical dispersion coefficient determines the vertical distribution of
pollutant in the plume and is a critical parameter in estimating ground level
receptor concentrations. The vertical dispersion coefficient increases with
turbulence and with distance from the source. The model allows the user to
select between two options for estimating this parameter. The first is the
Pasquill-Gifford family of curves derived using the expression:
az = a xb + c (8.8)
where a, b, and c are experimentally derived constants associated with the six
meteorological stability classes. The values of a, b, and c for the six
stability classes, which are listed in Table 8-2, are included in the code as
data statements.
76
-------
TABLE 8-2. Constants for PASQUILL-GIFFORD Curves for each Stability Class
Stability =
Stability =
a
b
c
x > 1000 m
a
b
c
.001
1.89
9.6
.0476
1. 1
2.00
. 1192
1.91
-25
.615
5 .4
.5-
2.63
5 . 1
126.
3.6
5 . 14
-75.
100 m <. x 1000 m
BCD
a .001
b 1.89
C 9.6
Stability = A
.0476
1. 11
2.00
x <
B
.119
.915
-1.4
100 m
C
.187
.755
-1.1
D
.1345
.745
-2.7
E
.362
.55
F
.1742
.936
0
.1426
. 922
0
.1233
.905
0
.0804
.881
0
.06
.854
0
. 0434
. 814
0
Source: U.S. EPA (1977;
The Pasquill-Gifford curves were derived from a series of experiments on
dispersion in flat terrain, and may not be applicable to all locations. If
data on wind turbulence are available for the site, the vertical dispersion
coefficient can be calculated by the model from the standard deviation of the
wind elevation angle. This is the second option:
o, =
(8.9)
a (x - r
h
where
ae = the standard deviation of the wind elevation angle [radians]
r0 = radius of source [m]
Note that the standard deviation of the wind elevation angle, ae, is also
referred to as the vertical turbulence intensity and is defined as the
standard deviation of vertical wind velocity fluctuations divided by the mean
wind speed. The radius of the source is calculated as the square root of the
area of the facility divided by two.
8.2.5 Terrain Effects
No gaussian model explicitly simulates the effects of terrain on meteorology
or on the trajectory of the plume (i.e., trapping in valleys, obstruction by
hills, etc.). However, MULTIMED does account for the effect of elevation on
the height of the plume centerline above the ground surface. The plume
centerline is assumed to remain at a constant elevation; the effective height
of the plume is then modified for receptors higher or lower than the source:
77
-------
H. =
,rh - (zr - zs
10
(Zr - Zs)h
(8.10)
78
-------
where
Ha = effective height of plume centerline above receptor [m]
h = the height above ground of the plume centerline at the source
[m]
zr = elevation of the receptor [m]
zs = elevation of the ground at the source, assumed zero in the
model [m]
When (zr - Zs) > h, the model assumes that the plume centerline is at the
ground surface.
The model also has two options for using a terrain correction, which is a
function of the atmospheric stability. Under stable conditions, the terrain
correction is calculated as described above. For unstable and neutral
conditions the model assumes that the plume always remains at the initial
height above the ground surface. Under these conditions the plume centerline
essentially parallels the ground surface.
8.3 Assumptions and Limitations
The following assumptions and limitations apply to the air dispersion module
Of MULTIMED:
(a) The concentration profile of the pollutant is assumed to be
gaussian in the transverse and vertical directions.
(b) Calculated concentrations are long-term averages and will not
indicate the maximum concentrations which might result from extreme
meteorological conditions.
(c) The model does not account for the effects of terrain on the shape
of the plume or the trajectory of the plume.
(d) The pollutant is assumed to be a non-buoyant vapor which neither
rises nor settles out of the plume due to gravitational forces.
(e) Sources are assumed to be square areas with constant emission
rates.
(f) Chemical transformation of the pollutant is modeled by lumped
first-order decay; no mechanistic descriptions of processes such as
photolysis are included in the model.
8.4 Data Requirements
Input data for the air dispersion module include parameters describing the
receptor, the source, and the local meteorology. The data requirements are
listed in Table 8-3. Note that in MULTIMED the source emission rates are
estimated by the air emissions module, described previously in Section 7.
The code allows the user either to assign a constant wind and stability
condition or to use a wind-stability frequency-weighting approach. For the
first condition, Equation 8.5 is solved using a constant wind speed and
constant stability condition, supplied in the input. The direction of the
wind is assumed to be the same as the direction of the receptor.
79
-------
For the second approach, a separate input file of joint frequency
distributions is required to solve Equation 8.5. This joint frequency
distribution of stability, wind speed, and wind direction is used to compute
seasonal average concentrations. For the usual configuration of 16 wind
direction sectors, 6 wind speeds, and 6 stability classes, 576 (16 x 6 x 6)
joint frequency entries are required. Typically this joint frequency
distribution is available as STAR (Stability Array data) summaries for
airports but may be difficult to find for other locations.
TABLE 8-3. Parameters Required for the Air Dispersion Module
Parameters Units
Source-Specific Parameters
Area of the land disposal facility
(converted to cm2 in the code) [m2]
Air Dispersion-Specific Parameters
Height of plume rise at a wind speed of 5 m/s [m]
Mixing layer thickness [m]
Standard deviation of wind elevation angle [radians]
Receptor elevation [m]
Receptor distance from the source [m]
Receptor angle from the source [radians]
Decay coefficient for pollutant in air [s"1]
Either:
Constant wind speed [m/s]
Constant stability condition [dimensionless]
or
Wind speeds for 6 classes [m/s]
Joint frequency distribution for each wind
speed, direction, and stability class [dimensionless]
80
-------
SECTION 9
Uncertainty Analysis
9.1 Introduction
As described in the previous sections, MULTIMED simulates the movement of
contaminants emanating from a waste disposal facility. The model includes
algorithms that simulate the movement of the contaminant within the
unsaturated zone, the saturated zone, a surface stream, and the atmosphere
based on a number of user-specified parameters. These include chemical-
specific, media-specific, source-specific and receptor location-specific
parameters.
Typically the values of these parameters are not known exactly due to
measurement errors and/or inherent spatial and temporal variability.
Therefore, it is often more appropriate to express their value in terms of a
probability distribution rather than a single deterministic value and to use
an uncertainty propagation model to assess the effect of the variability on
the model output.
This section presents the uncertainty propagation method implemented in the
MULTIMED code. The method allows a quantitative estimate of the uncertainty
in the concentration at a downgradient receptor location due to uncertainty in
the model input parameters.
9.2 Statement of the Problem and Technical Approach
The objective of the uncertainty analysis/propagation approach is to estimate
the uncertainty in the receptor concentration given the uncertainty in the
input parameters. In other words, the objective is to estimate the cumulative
probability distribution of the downgradient well concentration given the
probability distribution of the input parameters. As an example, Cw represents
the downgradient well concentration:
CW=9(X) (9.1)
where X represents the vector of model inputs and g represents the
computational algorithms for transport in the unsaturated and the saturated
zones. Some or all of the components of X may vary in an uncertain way (i.e.,
FCf (C ') = Probability (Cw < C ') (9 . 2 )
they are random variables defined by cumulative probability distribution
functions). The goal is to calculate the cumulative distribution function,
Fc (Cj , given a probabilistic characterization of X. Note that Fc (cj is
destined as: w
where C,!, is a given downgradient well concentration.
81
-------
Five methods of evaluating Fc (C'w) were examined in order to select the most
appropriate method for MULTlrffeo (WCC, 1986). The methods are:
1. First-Order and First-Order-Second-Moment Analysis (FO, FOSM);
2. Monte Carlo Simulation (MC);
3. Discretization of Probability Distributions (DPD);
4. Response Surface Analysis (RS); and
5. Rackwitz-Fiessler Method and its variants (RF).
The selection criteria included:
1. Computational efficiency, measured by the number of response
calculations required to achieve a given level of precision in
estimation of an output statistic (in this case, the 90th percentile
of the output distribution).
2. Accuracy in evaluation of the output statistic (e.g., a specified
percentile value).
3. Generality of application, so that a number of modules and input
conditions, and all sources of uncertainty, can be accommodated by
the same uncertainty-propagation method.
4. Simplicity of usage, measured by the number of parameters that must
be specified by the user for each application.
5. Completeness of the information produced, which may include only the
mean and variance of the output distribution or may be the whole
distribution, and which may or may not contain information useful
for uncertainty decomposition.
6. Flexibility with respect to input distributions, so that the method
would be able to accommodate a number of different input
distributions.
Using the above criteria, a qualitative comparison of the various uncertainty-
propagation methods is included in Table 9-1. Based on the evaluation of the
methods and a knowledge of MULTIMED, the Monte Carlo simulation method was
selected. This approach is simple, unbiased and completely general. Further
the method is especially attractive when there are many input variables that
are randomly distributed, because the efficiency does not depend on the
dimensionality of the input vector. Further, because the model is analytical,
it is not very expensive to run a large number of independent executions of
the model to achieve satisfactory confidence limits on the downgradient well
concentration. Details of the Monte Carlo method are discussed below.
82
-------
TABLE 9-1. Qualitative Comparison of Uncertainty-Propogation Methods
Criterion FO, FOSM MC DPP RS RF
Computational
Efficiency *** ** -- ** *
Accuracy * * * ** **
Generality ** *** * * *
Simplicity *** *** *** ** *
Information
Produced ** * ** ** ***
Variation of Fx ** ^ ^ *** *
- criteria not satisfied
* - criteria partially satisfied
** - criteria satisfied in general
*** - criteria satisfied
9.3 The Monte Carlo Analysis Technique
Figure 9.1 illustrates the Monte Carlo method used in MULTIMED. Given a set
of deterministic values for each of the input parameters, X1( X2 . . . xn, the
composite model computes the downgradient receptor concentration, Cw, i.e.:
cw = 9 (X^.X^.X^ . . . Xn) (9.3)
Application of the Monte Carlo simulation procedure requires that at least one
of the input variables, X^ . . . Xn, be uncertain, with the uncertainty
represented by a cumulative probability distribution. The method involves the
repeated generation of pseudo-random values of the uncertain input
variable(s). The pseudo-random values are drawn from the specified
distribution and are within the range of any imposed bounds. Then the model
is applied, using these values, to generate a series of model responses (i.e.,
values of C₯). These responses are statistically analyzed to yield the
cumulative probability distribution of the model response. Thus, the various
steps for the application of the Monte Carlo simulation technique involve:
(a) Selection of representative cumulative probability distribution
functions for the relevant input variables.
(b) Generation of a pseudo-random number from the distributions
selected in (a). These values represent a possible set of values
for the input variables.
(c) Application of the model to compute the derived inputs and
output (s) .
(d) Repeated application of steps (b) and (c).
(e) Presentation of the series of output (random) values generated in
step (c) as a cumulative probability distribution function (CDF).
83
-------
(f) Further analysis and application of the cumulative probability
distribution as a tool for decision making.
84
-------
9.4 Uncertainty in the Input Variables
The variables required by MULTIMED can be broadly classified into two
different sets that exhibit different uncertainty characteristics. These are:
(a) Variables that describe the chemical, biochemical, and
toxicological properties of the hazardous constituent. Examples of
these variables include the Henry's law constant, the acid,
neutral, and base catalyzed hydrolysis rates, and the soil
adsorption coefficient.
(b) Variables that describe the environmental properties of the various
media and impact the fate and transport of the pollutant within
each medium. Examples of these variables include the hydraulic
conductivity, porosity, organic carbon content, and dispersivity
values.
Uncertainty in the first set of variables primarily arises due to laboratory
measurement errors or theoretical analysis used to estimate the numerical
values. In addition to experimental precision and accuracy, errors may arise
due to extrapolations from controlled (laboratory) measurement conditions to
uncontrolled environmental (field) conditions. Further, for some variables,
semi-empirical methods are used to estimate the values. In this case, errors
in using the empirical relationships also contribute to variability in the
model outputs.
Uncertainty in the second set of variables may include both measurement and
extrapolation errors. However, the dominant source of uncertainty in these is
their inherent natural (spatial and temporal) variability. This variability
can be interpreted as site-specific or within-site variation in the event that
the model is used to analyze exposure due to a specific land-disposal unit.
Alternatively it can represent a larger scale (regional/national) uncertainty
if the model is used to conduct exposure analysis for a specific chemical or
specific disposal technology on a generic, nation-wide or regional basis.
Note that the distributional properties of the variables may change
significantly depending upon the nature of the application.
Whatever the source of uncertainty, the uncertainty preprocessor developed for
the model requires that the uncertainty be quantified by the user. This
implies that for each input parameter deemed to be uncertain, the user selects
a distribution and specifies the parameters that describe the distribution.
Currently, the user can select one of the following distributions:
(a) Normal
(b) Lognormal
(c) Uniform
(d) Log uniform
(e) Exponential
(f) Empirical
(g) Johnson SB distribution.
The first two distributions require the user to specify the mean and the
variance. The third and the fourth require minimum and maximum values. The
fifth distribution requires only one parameter--the mean of the distribution.
For the empirical distribution, the user is required to input the coordinates
of the cumulative probability distribution function (minimum 2 pairs, maximum
20 pairs) which is subsequently treated as a piece-wise linear curve.
85
-------
Finally, the Johnson SB distribution requires four parameters--mean, variance,
the lower and upper bounds.
In all cases, MULTIMED requires the upper and lower bounds for all the
distribution types. If a value which is randomly generated during any Monte
Carlo run is outside the user-specified bounds, it is rejected by the code and
a new number is generated.
Of the seven distributions, the characteristics of the first six are readily
available in literature (Benjamin and Cornell, 1970). Details of the Johnson
system of distributions are presented in McGrath and Irving (1973) and Johnson
and Kotz (1970). The distribution types are briefly summarized in Section 9.5
below. Note that Section 9.5 is copied, with slight modifications, from
Volume 1 of the RUSTIC documentation (Dean et al., 1989).
9.5 Description of the Parameter Distributions
This description of the distributions which can be specified in the MULTIMED
Monte Carlo option includes: 1) the parameters of the distributions, 2) the
equations for the probability and cumulative density functions, and 3) a brief
discussion of the properties of each distribution.
9.5.1 Normal Distribution
The term "normal distribution" refers to the well known bell-shaped
probability distribution. Normal distributions are symmetrical about a mean
value and are unbounded, although values further from the mean occur less
frequently. The spread of the distribution is generally described by the
standard deviation. The normal distribution has only two parameters: the
mean and the standard deviation. Although the distribution is unbounded, the
user must enter minimum and maximum values for individual parameter
distributions. Generated values outside these bounds are not used by the
model.
The probability density function of x is given by:
I I x - mxl2!
exp irO.5 i)))))) i.i, (9.4)
L I Sx ) J
where Sx is the standard deviation, and mx is the mean of x. The cumulative
distribution is the integral of the probability density function:
x
Fn(x) = J fn(x)dx (9.5)
The above integration must be performed numerically, but tables of
numerically-integrated values of Fn(x) are widely available in the statistical
literature.
9.5.2 Log-Normal Distribution
The log-normal distribution is a skewed distribution in which the natural log
of variable x is normally distributed. Thus, if y is the natural log of x,
then the probability distribution of y is normal with mean, my, and standard
-------
deviation, Sy, and a probability density function similar to Equation 9.4. The
mean and standard deviation of x (mx and Sx) are related to the log-normal
parameters m,f and Sy as follows:
mx = exptmy +0.5 (Sy)2] (9.6)
Sx2 = mx2 [exp(Sy2) - 1] (9.7)
To preserve the observed mean and standard deviation of x, the parameters of
the log-normal distribution (m,f and Sy) are therefore selected such that the
above relationships are satisfied. Note that my and Sy do not equal the
natural logs of mx and Sx respectively. Log-normal distributions have an
absolute lower bound of 0.0 and no absolute upper bound, and are often used to
describe positive data with skewed observed probability distributions. As
stated above, within the absolute bounds, the user must input minimum and
maximum bounds for individual parameters.
Note that when a lognormal distribution is selected in MULTIMED, the user
should input the mean and standard deviation of the data (arithmetic space).
The transformation to the mean and standard deviation in lognormal space is
performed by the code.
9.5.3 Uniform Distribution
A uniform distribution is a symmetrical probability distribution in which all
values within a given range have an equal chance of occurrence. A uniform
distribution is completely described by two parameters: 1) the minimum value
(lower bound), A, and 2) the maximum value (upper bound), B. The equation for
the uniform probability density distribution of variable x is given by:
fu(x) = 1/(B - A) (9.8)
where fu(x) is the value of the probability density function for x. The
cumulative distribution F(x) is obtained by integrating Equation 9.8. This
yields the probability distribution:
Fu(x) = (x - A)/(B - A) (9.9)
where F(x) is the probability that a value less than or equal to x will occur.
9.5.4 Log-uniform Distribution
The log-uniform distribution type requires only a lower and an upper bound;
the mean and standard deviation are not used. This distribution results in
values between the lower and upper bounds having equal probability of
occurrence (as in a standard uniform distribution) but with relative
magnitudes that follow a logarithmic scale.
9.5.5 Exponential Distribution
The probability density function for an exponential distribution is described
by an exponential equation:
87
-------
fe(x) =
exp (-x/mj
(9.10)
-------
where mx is the mean of x. The cumulative distribution is given by:
Fe(x) = 1 - exp(-x/mx) (9.11)
The exponential distribution is bounded by zero; the probability density
function peaks at zero and decreases exponentially as x increases in
magnitude.
9.5.6 Empirical Distribution
At times it may be difficult to fit a standard statistical distribution to
observed data. In these cases it is more appropriate to use an empirical
piecewise-linear description of the observed cumulative distribution for the
variable of interest.
Cumulative probabilities can be estimated from observed data by ranking the
data from lowest (rank = 1) to highest (rank = number of samples) value. The
cumulative probability associated with a value of x is then calculated as a
function of the rank of x and the total number of samples. The cumulative
probabilities of values between observed data can be estimated by linear
interpolation.
9.5.7 The Johnson System of Distributions
The Johnson system involves three main distribution types, one of which, the
SB (Log-ratio or bounded), is included in MULTIMED. The Johnson SB
distribution basically represents a log-ratio or bounded transformation
applied to the random variable such that the transformed variable is normally
distributed.
SB: Y = to(-|gL) (9.12)
where
In = natural logarithm transformation
X = untransformed variable with limits of variation from A to B
Y = the transformed variable with a normal distribution
To determine if the Johnson SB distribution is appropriate for a sample data
set, the skewness and kurtosis of the sample data should be plotted as shown
in Figure 9.2. If the sample point is located in the region for the SB
distribution, it can be used for the sample data. For additional details
about the Johnson system of distributions, the reader is referred to McGrath
et al. (1973) and Johnson and Kotz (1970).
9.6 The Random Number Generator
Having selected a distribution for each of the input parameters to be Monte
Carloed, the model generates random values of these parameters. This requires
the use of pseudo-random number generating algorithms. Numerous non-
proprietary subroutines exist that can be used to generate random numbers. A
number of these are comparable in terms of their computational efficiency,
accuracy and precision. The specific routines included in this code are those
described by McGrath et al. (1973).
The performance of these algorithms has been checked to ensure that they
accurately reproduce the parameters of the distributions that are being
89
-------
sampled. In order to test the algorithms two sets of runs were made. For Run
1, 500 random numbers were generated; for Run 2, 1000 random numbers were
generated. For the five distributions tested, the input parameters and the
results are shown in Tables 9-2 (a) and (b) . In each case, the output
statistics for the randomly generated variables closely match the input
values .
For Run 2, the randomly generated variables were arranged in ascending order
and the cumulative probability distributions were plotted. The results are
shown in Figures 9.3 to 9.7. Visual inspection of these figures further
testifies to the accuracy of these algorithms.
Note that more rigorous statistical tests could be used to further test the
accuracy of the algorithms. However, the above simplified analysis and the
additional testing performed by Marin (1988) has provided sufficient proof of
the accuracy of the results and indicated that these algorithms satisfactorily
reproduce the input statistics and distributions of the variables.
9.7 Analysis of the Model Output
Using the randomly generated parameter values, the model is used to estimate
values of concentrations at a point located downgradient from the waste
facility. If Cw represents the normalized, steady-state concentration at the
receptor location calculated by the model when the leachate concentration at
the waste disposal facility is unity, and CT is the (health based maximum
allowable) threshold concentration for the chemical at the receptor, the
maximum allowable leachate concentration at the waste facility can be back-
calculated by hand using:
CT
Ct = - (9.13)
The maximum allowable leachate concentration defined by Equation 9.13 is the
leachate concentration for which the downgradient receptor concentration does
not exceed the threshold concentration.
Alternatively,
J_ = _C«_ (9. 14)
C C
^w T
Equation 9.14 states that the reciprocal of the computed normalized
concentration represents the maximum allowable ratio of leachate concentration
to the threshold concentration. Thus, for a simulated normalized
concentration of Cw = 0.05 mg/f, Equation 9.14 implies that the maximum
allowable leachate concentration from the landfill could be 20 times the
threshold value for the chemical. Note that both Cw and CT are chemical
specific .
Combining the above back calculation procedure with the Monte Carlo analysis
allows the maximum leachate concentration to be couched in a probabilistic
framework. For each chemical, application of the Monte Carlo method results
in an array of values for normalized concentration, each representing a
feasible result for the waste disposal facility-environmental scenario. These
values are statistically analyzed to derive the cumulative probability
90
-------
distribution function as shown in Figure 9.8. The cumulative probability
distribution, Fc (cj together with the allowable threshold value, CT, and the
back calculation^17 procedure (Equations 9.13 and 9.14) provide the information
necessary to calculate the maximum allowable leachate concentration. In
particular, the value of leachate concentration, C,, that leads to p% of the
realizations in compliance--i.e., the receptor well concentration is less than
or equal to the threshold concentration, is:
c = i
' CP
(9.15)
where Cp is the p-percentile concentration obtained from the cumulative
distribution function of the downgradient well concentration.
9.8 Confidence Bounds for the Estimated Percentile
As described above, the Monte Carlo simulation provides an estimate of Cp, the
p-percentile concentration obtained from a sample of n model simulations.
Since the sample size is finite, estimates of Cp will be uncertain, with the
degree of uncertainty decreasing with increasing sample size.
A quantitative estimate of the uncertainty in the estimate of Cp can be
obtained by computing a confidence interval around the estimate Cp. The upper
and the lower bounds of this confidence interval are defined such that:
Probability (CL < Cp < Cu) = 1 - a
(9.16)
where
CL = the lower bound of the confidence interval [mg/f]
Cn = the upper bound of the confidence interval [mg/f]
a = significance level [dimensionless]
The interval CL to Cn is usually referred to as the 100(1 - a) confidence
interval and implies that the true value of the estimate of the quantile Cp
lies within the interval CL to GU with a probability of 100(1 - a). The
confidence interval is estimated using the binomial distribution as described
below.
Define an indicator (Bernoulli) random variable, I, such that:
if CLCL
Also, define a random variable, K, equal to the number of trials for which I =
1 from the Monte Carlo simulations. The random variable K is then binomially
distributed with a mean of np and a variance of np(l - p), i.e.,
Probability^ I± = k} =
nl
k\ (n - k)
pk (1 - p)1
(9.18)
91
-------
where n is the number of independent realizations of Cw (by Monte Carlo
simulations) corresponding to n independent realizations of I.
The probability that K is less than a given positive integer is also the
probability that C(K) < Cp where C(K) is the Kth smallest simulated value of the
concentration. Thus a confidence interval on K, based on the binomial
distribution, can be used to establish the confidence interval, Cp.
Conover (1980) describes a procedure for obtaining a confidence interval for Cp
that essentially involves looking up values in a table of cumulative
probabilities of the binomial distribution (i.e., Prob {K <_ k}) to find the k
values corresponding to probabilities of a/2 and 1 - a/2. For sample sizes
greater than 20, he provides an alternate procedure based on the normal
approximation to binomial probabilities for large n. This simple
approximation requires the calculation of two values, r and s:
r = np + za/2[np(l - p)]1/2 (9.19)
s = np + z1.a/2 [np (1 - p) ] 1/2 (9.20)
where za/2 and zI_a/2 are quantiles of the standard normal (mean = 0, variance =
1) distribution (note that za/2 = -z1_a/2] . Conover recommends rounding up these
values of r and s to the next higher integers and estimate the corresponding
values of C(r) and C(s), as the rth and sth smallest values of Cw. The confidence
interval is then of the form:
Prob{C{r}
-------
An alternative criterion for specifying n might be the fraction of the range
of simulated ranks to be covered by the confidence interval. Thus,
(g
n n
or, solving for n:
n- (9.24)
93
-------
SECTION 10
References
Ambrose, R.B., L.A. Mulkey, and P.S. Huyakorn (1985), "A Methodology for
Assessing Surface Water Contamination due to Land Disposal of Hazardous
Waste." USEPA, Athens, GA.
Anatra, M., S. Bosch, P. Durkin, D.A. Gray, and D. Hohreiter (1986),
"Investigation of the Uncertainty in the Acceptable Daily Intake (ADI)."
Syracuse Research Corporation report submitted to Woodward-Clyde
Consultants, USEPA Project No. 68-03-6304, Cincinnati, Ohio.
Anderson, M.P. (1979), "Using Models to Simulate the Movement of Contaminants
through Groundwater Flow System." CRC Critical Reviews in Environmental
Control, Vol. 9, Issue 2, pp. 97-156.
Arthur D. Little, Inc. (1984), "Evaluation of Emission Controls for Hazardous
Waste Treatment, Storage and Disposal Facilities." EPA-650/3-84-017.
Office of Air Quality Planning and Standards, Research Triangle Park,
North Carolina.
Bear, J. (1979), Hydraulics of Groundwater, McGraw Hill, New York.
Benjamin, J.R., and C.A. Cornell (1970), Probability, Statistics, and Decision
for Civil Engineers. McGraw Hill Book Company.
Bicknell, Brian R. (1984), "Modeling Chemical Emissions from Lagoons and
Landfills." Anderson-Nichols and Co., Inc., Palo Alto, CA 94303.
Brooks, R.H., and A.T. Corey (1966), "Properties of Porous Media Affecting
Fluid Flow." ASCE J. Irrig. Drain. Div. 92 (2) -.61-68.
Burns, L.A., D.M. Cline, and R.R. Lassiter (1982), "Exposure Analysis Modeling
System (EXAMS): User Manual and System Documentation." USEPA, -600/3-
82-023, Athens, GA.
Carnahan, B., H.A. Luther, and J.O. wilkes (1969), "Applied Numerical
Methods", John Wiley.
Carsel, R.F., and R.S. Parrish (1988), "A Method for Developing Joint
Probability Distribution of Soil-Water Retention Characteristics."
Water Resources Research 24(5):755-769.
94
-------
Carsel, R.F., R.S. Parrish, R.L. Jones, J.L. Hansen, R.L. Lamb (1985),
"Characterizing the Uncertainty of Pesticide Leaching in Agricultural
Soils." Draft submitted to J. Env. Qual.
Carsel, R.F, et al. (1985), "Characterizing the Uncertainty of Pesticide
Leaching in Agricultural Soils." Draft Report U.S. EPA Environmental
Research Laboratory, Athens, Georgia 30613.
Carsel, R.F., C.N. Smith, L.A. Mulkey, J.D. Dean, and P. Jowise. (1984)
User's Manual for the Pesticide Root Zone Model (PRZM) Release 1.
EPA-600/3-84-109, 16-17.
Carsel, R.F., and R.S. Parrish. (1988) A Method for Developing Joint
Probability Distribution of Soil-Water Retention Characteristics. Water
Resource Research 24(5), 755-769.
Clements. W.F., and M.H. Wilkening (1974), "Atmospheric Pressure Effects on
Rn Transport across the Earth-Air Interface." J. Geophysical
Research, Vol. 79, No. 33.
Code of Federal Regulations. 40 CFR. Part 264.
Cohen, Y. (1986), "Organic Pollutant Transport. Improved Multimedia Modeling
Techniques Are the Key to Predicting the Environmental Fate of Organic
Pollutants." Environ. Sci. Technol., Vol. 20, No. 6.
Conover, W.J. (1980), Practical Nonparametric Statistics. 2nd ed., John Wiley
and Sons, New York, 493 pp.
Covar, A.P. (1976), "Selecting Proper Recreation Coefficient for Use in Water
Quality Models." Presented at the USEPA Conference on Environmental
Simulation and Modeling, April 19-22, 1976.
CRC (1981), Handbook of Chemistry and Physics, 62nd edition, CRC Press.
Currie, J.A. (1961), "Gaseous Diffusion in Porous Media. Part 3 - Wet
Granular Materials." British Journal of Applied Physics, June 1961.
Currie, J.A. (1960), "Gaseous Diffusion in Porous Media. Part 2 - Dry
Granular Materials." British Journal of Applied Physics, Vol. 11,
August 1960.
Cussler, E.L. (1983), "Diffusion: Mass Transfer in Fluid System," pp. 188-
189, 1st Ed., Cambridge University Press.
Dass, P., G.R. Tamke, and C.M. Stoffel. (1977) Leachate Production at
Sanitary Landfills. ASCE, Jour. Env. Eng. Div., 103 (EE6).
E.G. Jordan Co. (1987), Technical Memorandums dated June 2, 1987, and
September 1987, submitted to USEPA, OSW, Washington, D.C.
E.G. Jordan Co. (1985), "Analysis of Engineered Controls of Subtitle C
Facilities for Land Disposal Restrictions Determinations. Revised
Distribution of Leaching Rates." Draft Report ECJ Project No. 4756-01
prepared for Research Triangle Institute, North Carolina, and USEPA,
OSW, Washington, D.C.
95
-------
Electric Power Research Institute (1985), "A Review of Field Scale Physical
Solute Transport Processes in Saturated and Unsaturated Porous Media.'
EPRI EA-4190, Project 2485-5, Palo Alto, California.
-------
Enfield, C.G., et al. (1982), "Approximating Pollutant Transport to Ground
Water." Ground Water, Vol. 20, No. 6, pp. 711-722.
Farmer, W.J., et al. (1978), "Land Disposal of Hexachlorobenzene Wastes:
Controlling Vapor Movement in Soils." Proceedings of the 4th Annual
Research Symposium held at San Antonio, Texas, March 6-8.
Federal Register (1986), "Hazardous Waste Management System: Land Disposal
Restrictions." USEPA, Vol. 15, No. 9.
Freeze and Cherry (1979), "Groundwater." Prentice Hall, Englewoods Cliffs,
New Jersey.
Freeze, R.A., and J.A. Cherry. (1979) Groundwater. Prentice-Hall Inc.,
Englewood Cliffs, N.J., 30-36.
Fukuda, H. (1955), "Air and Vapor Movement in Soil Due to Wind Gustiness."
Soil Science 79:249-258.
Gelhar, L.W. and Axness, C.J. (1981), "Stochastic Analysis of Macro-Dispersion
in Three-Dimensionally Heterogeneous Aquifers." Report No. H-8,
Hydraulic Research Program. New Mexico Institute of Mining and
Technology, Soccorro, New Mexico. 140 p.
Groves, F.R., D.D. Reible, and L.J. Thibodeaux (1984), "Estimation of Physical
and Chemical Properties of Waste Organic Mixtures Associated with Land
Pollution." Section 3 in Exposure Assessment Involving Mixtures of
Environmental Pollutants, Ed. J.D. Dean et al. , Office of Health and
Environmental Assessment, DRD, USEPA, Washington, D.C.
Haderman, J. (1980), "Radionuclide Transport through Heterogenous Media."
Nuclear Technology 47:312-323, February 1980.
Hermann, D.J., and S.B. Tuttle. (1981) Performance Studies of Various
Landfill and Impoundment Liners. Fourth Annual Madison Conference of
Applied Research and Practice on Municipal and Industrial Waste.
University of Wisconsin, Madison, WI, 250-265.
Huyakorn, P.S., M.J. Ungs, L.A. Mulkey, and E.A. Sudicky (1987), "A Three-
Dimensional Analytical Method for Predicting Leachate Migration."
Groundwater, Vol. 25, No. 5. September-October 1982.
Huyakorn, P.S., M.J. Ungs, E.D. Sudicky, L.A. Mulkey, and T.D. Wadsworth.
1986. "RCRA Hazardous Waste Identification and Land Disposal Restriction
Groundwater Screening Procedure." USEPA, Washington, D.C.
Irwin, J.S., T. Chico, and J. Catalano. 1985. CDM 2.0 -- Climatological
Dispersion Model -- User's Guide. U.S. Environmental Protection Agency,
Research Triangle Park, North Carolina.
Israelsen, D.W., and V.E. Hansen (1962), "Irrigation Principles and
Practices." John Wiley and Sons, Inc., New York. 447 pages.
Johnson, N.L., and S. Kotz (1970), Distributions in Statistics: Continuous
Univariate Distributions, Houghton Mifflin Company, Boston.
97
-------
Karickhoff, S.W. (1984), "Organic Pollutant Sorption in Aquatic Systems."
ASCE J. Hyd. Div. 110(6):707-735.
Karickhoff, S.W., D.S. Brown, and T.A. Brown (1979), "Sorption of Hydrophobic
Pollutants on Natural Sediments." Water Res. 13:241-248.
98
-------
Knisel, W.J., Jr., Ed. (1980) CREAMS, A Field Scale Model for Chemical
Runoff and Erosion from Agricultural Management Systems. Vols. I, II,
and III, Draft Copy, USDA-SEA, AR, Cons. Res. Report 24.
Kulkarni R.B., G. Luster, G. Rao and A.M. Salhotra (1988) Enhanced Methods for
Characterizing Uncertainties in Numerical Models. Volume I:
Methodology Development and Validation and Volume II: User's Manual and
Programmers Guide. Report prepared for U.S. Environmental Protection
Agency, Office of Research and Development, Athens, GA Contract No. 68-
03-6309
Lester, B.H., P.S. Huyakorn, H.O. White, Jr., T.D. Washworth, and J.E. Buckley
(1986), "EPACML Composite Analytical-Numerical Model for Simulating
Leachate Migration in Unconfined Aquifers." Prepared by GeoTrans, Inc.
Liss, P.S. (1973), "Process of Gas Exchange across an Air-Water Interface."
Deep-Sea Research 20:221-238.
Lu, A.H., and J.M. Matuszek (1978), "Transport through a Trench Cover of
Gaseous Tritiated Compounds from Buried Radioactive Waste," presented at
Symposium on the Behavior of Tritium in the Environment, held at San
Francisco.
Lu, J.C.S., B. Eichenberger, and R.J. Stearns. (1985) Leachate from
Municipal Landfills Production and Management. Noyes Publications, Park
Ridge, N.J., 7-94.
Lyman et al. (1982), "Handbook of Chemical Property Estimation Methods.
Environmental Behavior of Organic Compounds." McGraw-Hill Book Company.
Marin, C. (1988), Personal Communication.
Marino, M.A. (1974), "Distribution of Contaminants in Porous Media Flow."
Water Resources Research 10 (5) : 1013-1018.
McGrath, E.J., and D.C. Irving (1973), Techniques for Efficient Monte Carlo
Simulation, Volume II. Random Number Generation for Selected
Probability Distributions. Report prepared for Office of Naval
Research. Project No. NR 364-074/1-5-72, Code 462.
McKoy and Associates (1986), "Predicting Volatile Organic Emissions through
Soil Covers." The Hazardous Waste Consultant, September 1986.
Meeks, Y, P. Mangarella, G. Palhegyi, and A. M. Salhotra (1988), "Landfill
Source Model." Report Prepared for U.S. Environmental Protection Agency
under EPA Contract No. 68-03-6304. Environmental Research Laboratory,
Athens, Georgia.
Mill, T., et al. (1981), "Laboratory Protocols for Evaluating the Fate of
Organic Chemicals in Air and Water." Final Draft, Prepared for U.S. EPA
Technology Development and Applications Branch under EPA Contract No.
68-03-2227. Environmental Research Laboratory, Athens, Georgia.
Millington, R.J., and J.P. Quirk (1961), "Permeability of Porous Solids."
Trans. Faraday Soc. 57:1200-1207.
99
-------
Moench, A.F., and A. Ogata (1981), "Numerical Inversion of the Laplace
Transform Solution to Radial Dispersion in a Porous Medium." Water
Resources Research 17 (1) -.250-252 .
Mood, A.M., F.A. Graybill, and D.C. Boes (1974), Introduction to the Theory of
Statistics, third edition. McGraw-Hill Book Co., New York. 564 pages.
100
-------
Mualem, Y. (1976), "A New Model for Predicting the Hydraulic Conductivity of
Unsaturated Porous Media." Water Resources Research 12(3):513-522.
Mulkey, L.A., and T. Allison (1988), "Transient versus Steady-State Land
Disposal Model Comparisons." Report prepared for EPA OSW Environmental
Research Laboratory.
Olkin, I., L.J. Gleser, and C. Derman (1978), Probability Models and
Applications. Macmillan, New York, 576 pp.
Perrier, E.R., and A.C. Gibson (1980), "Hydrologic Simulation on Solid Waste
Disposal Sites." SW-868, U.S. E.P.A, Cincinnati, OH.
Renyi, A. (1953), "On the Theory of Order Statistics." Acta Mathematica
Academiae Scientiarum Hungarical 4(3-4) :191-231.
Schanz, R., and A.M. Salhotra (1988), "A Water Treatment Plant Model for
Pollutant Exposure Assessment." Final Report submitted to Environmental
Research Laboratory, Office of Research and Development, Athens, GA.
Schroeder, P.R., J.M. Morgan, T.M. Waloki, and A.C. Gibson (1984), "The
Hydrologic Evaluation and Landfill Performance (HELP) Model." Interagency
Agreement Number AD-96-F-2-A140, U.S. Army Engineer Waterways Experiment
Station, Vicksburg, MS 39180.
Shamir, V.Y., and D.R.F. Harleman (1967), "Dispersion in Layered Porous
Media." Journal of Hydraulics Division, ASCC-HYS, pp. 237-260.
Shen, Thomas T. (1981), "Estimating Hazardous Air Emissions from Disposal
Sites." Pollution Engineering.
Shen, T.T. (1979), "Emission Estimation of Toxic Pollutants from Hazardous
Waste Disposal Sites." Unpublished report, Division of Air Resources,
NYS Department of Environmental Conservation.
Shen, Thomas T., and T.J. Tofflemire (1980), "Air Pollution Aspects of Land
Disposal of Toxic Wastes." ASCE, Vol. 106 EE1.
Salhotra, A. (1988), "Dispersivity Values for the Unsaturated Zone," Technical
Memo submitted to EPA/OSW dated February 10, 1988.
Schroeder, P.R., J.M. Morgan, T.M. Walski, and A.C. Gibson (1984a),
"Hydrologic Evaluation of Landfill Performance (HELP) Model: Volume I.
User's Guide for Version 1." Municipal Environmental Research
Laboratory, U.S. Environmental Protection Agency, Cincinnati, OH.
EPA/530-SW-84-009.
Schroeder, P.R., A.C. Gibson, and M.D. Smolen (1984b), "Hydrologic Evaluation
of Landfill Performance (HELP) Model: Volume II. Documentation for
Version 1." Municipal Environmental Research Laboratory, U.S.
Environmental Protection Agency, Cincinnati, OH. EPA/1530-SW-84-009.
Sharp-Hansen, S., C. Travers, P. Hummel, and T. Allison (1990), "A Subtitle D
Landfill Application Manual for the Multimedia Exposure Assessment Model
(MULTIMED)." Prepared for the Environmental Research Laboratory, U.S.
Environmental Protection Agency, Athens, GA.
101
-------
Spear, R.C. (1970a), "The Application of Kolmogorov-Renyi Statistics to
Problems of Parameter Uncertainty in Systems Design." International
Journal of Control 11 (5) : 771-778 .
102
-------
Spear, R.C. (1970b) , "Monte Carlo Method for Component Sizing." Journal of
Spacecraft and Rockets 7(9):1127-1129.
Springer, C., P.O. Lunrey, K.T. Valsaraj, L.J. Thibodeaux, and S.C. James
(1986), "Emissions of Hazardous Chemicals from Surface and Near Surface
Impoundments to Air," Part B Landfills, Project No. CR 808161-02, Office
of Research and Development, USEPA, Cincinnati, Ohio.
Stehfest, H. (1970), "Numerical Inversion of Laplace Transforms." Commun. ACM
13(1):47-49.
Sudicky, E.A., M.J. Ungs, and P.S. Huyakorn (1986), "Three Dimensional
Analytical Solution for Transport from Gaussian Source in Uniform Ground
Water Flow." Submitted to Water Resources Research.
Sudicky, E.A., T.D. Wadsworth, J.B. Kool, and P.S. Huyakorn (1988), "PATCH3D -
- Three Dimensional Analytical Solution for Transport in a Finite
Thickness Aquifer with First-Type Rectangular Patch Source." Report
prepared for Woodward-Clyde Consultants. January 1988.
Swann, R.L., Eschenroeder, A., eds. (1983), "Fate of Chemicals in the
Environment." ACS Symposium Series 225, American Chemical Society,
Washington, D.C.
Todd, O.K. (1970) The Water Encyclopedia: A Compendium of Useful
Information on Water Resources. Water Information Center, Port
Washington, N.Y.
Thibodeaux, L.J., C. Springer, and G. Hildebrand (1986), "Transport of
Chemical Vapours through Soil--a Landfill Cover Simulation Experiment,"
presented at 1986 Summer National American Institute of Chemical
Engineers, August 24-27, 1986, Boston, Massachusetts.
Thibodeaux, L.J., et. al. (1982), "Models of Mechanisms for the Vapor Phase
Emission of Hazardous Chemicals from Landfills." Jour, of Hazardous
Materials 7 (1982) : 63-74 .
Ungs, M. J. (1986), "Mathematical Analysis and Validation of the Monte Carlo
Uncertainty Analysis for the Surface Water Component of the Land
Disposal Restriction Determinations." Draft Report Prepared for U.S.
EPA Office of Solid Waste under EPA Contract No. 68-01-7266, Work
Assignment No. 13.
Ungs, M.J. (1987), Unpublished manuscript submitted to Woodward-Clyde
Consultants and U.S. EPA ERL Athens.
USDA, Soil Conservation Service. (1972) National Engineering Handbook,
Section 4, Hydrology. U.S. Government Printing Office, Washington, D.C.
USEPA (1985), "Development of Land Disposal Banning Decisions Under
Uncertainty." USEPA Environmental Research Laboratory, Athens, Georgia.
USEPA (1984), "Evaluation and Selection of Models for Estimating Air Emissions
from Hazardous Waste Treatment, Storage and Disposal Facilities." (EPA-
450/3-84-020) Office of Air Quality Planning and Standards, Research
Triangle Park, North Carolina.
103
-------
USEPA (1977), User's Guide to the VALLEY Model, Office of Air Quality Planning
and Standards, Research Triangle Park, North Carolina.
Vandergrift, S.B., and R.B. Ambrose, Jr. (1988), "SARAH2, A Near Field
Exposure Assessment Model for Surface Water." U.S. EPA, Athens, GA.
EPA/600/3-88/020.
Van Genuchten, M. (1976), "A Closed Form Equation for Predicting the Hydraulic
Conductivity of Unsaturated Soils." Soil 3d. Soc. J. 44 (5) : 892-898 .
Van Genuchten, M., and W.J. Alves (1982), "Analytical Solutions of the One-
dimensional Convective-Dispersive Solute Transport Equation." Technical
Bulletin No. 1611, United States Department of Agriculture.
Van Genuchten, M. Th., G.F. Pinder, and W.P. Saukin (1977), "Modeling of
Leachate and Soil Interactions in an Aquifer, Management of Gas and
Leachate in Landfills." Rep. EPA -600/9-77-026, pp. 95-102, USEPA,
Cincinnati, Ohio.
Veneziano, D., R. Kulkarni, G. Luster, G. Rao, and A. Salhotra (1987, in
press), "Improving the Efficiency of Monte Carlo Simulation for Ground-
Water Transport Models." Proceeding of DOE/AECL 87 Conference on
Geostatistical, Sensitivity, and Uncertainty Methods for Ground-Water
Flow and Radionuclide Transport Modeling, San Francisco, 1987.
Wark, K., and C.F. Warner (1981), Air Pollution -- Its Origin and Control.
Harper and Row. New York.
Whitman, R.G. (1923), "A Preliminary Experimental Confirmation of the Two Film
Theory of Gas Absorption." Chem. Metallurgy Eng. 29:146-148.
Wolfe, N.L. (1985), "Screening of Hydrolytic Reactivity of OSW Chemicals."
USEPA Athens Environmental Research Laboratory, Athens, Georgia.
Wolfe, N.L. (1987), Personal Communication.
Woodward-Clyde Consultants (1986), "Development of an Approach for Conducting
Uncertainty Analyses in Multimedia Modeling." Draft Report Prepared for
USEPA Technology Development and Applications Branch, Athens. Project
No. 68-03-6304. Work Assignment B-2.
Woodward-Clyde Consultants. (1988) SYNOP, Rainfall Analysis Algorithm.
Report under preparation.
104
-------
APPENDIX A
Simplified Estimation for the Mixing Zone Depth
The mixing zone depth of a solute plume under a waste facility can be
estimated by adding the contributions due to advection and dispersion, i.e.:
H = h*dv + hdiSP (A.I)
where
H = the total depth of penetration [m]
hadv = the vertically advected component of the penetration depth [m]
hdlsp = the vertically dispersed component of the penetration depth [m]
The advected depth, hadv, is the depth to which a particle would be transported
under the influence of vertical advection and is given by:
v^t (A. 2)
where
vz = the vertical seepage velocity [m/yr]
T = time of travel [yr]
If the vertical seepage velocity is constant with depth, then:
hadv = VZt (A. 3)
However, a better assumption is that the vertical seepage velocity varies
linearly with depth, with a maximum value at the water table and zero at the
bottom of the aquifer. This variation can be mathematically expressed as:
Vz = Vzo(l-z/B) (A. 4)
where
B = the saturated aquifer thickness [m]
z = the depth from the top of the water table [m]
Vzo = the maximum vertical seepage velocity [m/yr]
vzo can be estimated from the net vertical recharge rate.
Equation 2 cannot be integrated since Vz is not an explicit function of time.
Consider the following differential equation for the vertical seepage
velocity:
105
-------
dz
dt
= V (z)
(A.5)
106
-------
Rearrange terms in Equation 5 and integrate to depth hadv:
-, dz
Substitute Equation 4 into Equation 6 and integrate to get:
-^ ind-Auiv/B) = T (A. 7)
zo
Solve for hadv from Equation 7:
J»l (A.8)
h .. = B (1 - e B )
The time of travel, i, can be estimated as the time it takes for a particle to
be advected horizontally under the facility of length, L, i.e.,
where vx is the horizontal seepage velocity and is assumed to be a constant.
Prickett, Naymik, and Lonnquist (1981) estimate the magnitude of the effect of
dispersion on particle transport as:
Aio,, = ^/2^W (A. 10)
(A. 11)
107
-------
where
«L, av = the longitudinal and vertical dispersivities [m]
V = the magnitude of the seepage velocity [m/yr]
Along, Avert = the longitudinal and vertical dispersed distances that
correspond to one standard deviation of random transport
[m]
If the effect of the horizontal seepage velocity is assumed to be much larger
than that of the vertical, then the dispersed depth is estimated from Equation
11 as:
(A.12)
Hence, the total depth of penetration is the sum of the vertically advected
and dispersed components. Substitute Equations A.8 and A.12 into Equation A.I
to obtain the total estimated depth of penetration.
H = B(l-e
(A.13)
The solution to Equation A.13 needs to be checked when evaluating any
particular case so that a value of H greater than the aquifer thickness, B, is
not used. If the computed H is greater than B, set H equal to B.
Reference
Prickett, T., T. Naymik, and C. Lonnquist (1981), A Random-Walk Solute
Transport Model for Selected Groundwater Quality Evaluations. Bulletin
65, Illinois State Water Survey, Department of Energy and Natural
Resources, Champaign, Illinois. 103 pages.
108
-------
UJ
tl
cr
ft
Ll_
CO
w
c
u
D.
Q.
O
^
< k
UJ
u!
0
£C
C.
m
w
tr
5
0
i >
PRECIPITATION
I 1 1 EVAPOTRANSPIRATION
INFLATION ^|^ | ^^
lllill!ll!ll!ll!ll!llllllll!lll!!l!!!l!l!lll!l!!lilllll!l!!l!llll!lllllill!llllllli^ Vegetation
VEGETATIVE LAYER . |
cr
UJ
s
LATERAL DRAINAGE LAYER ^
0
o-
o
BARRIER SOIL LAYER
1
REFUSE LAYER
1
p
LATERAL DRAINAGE f
LATERAL DRAINAGE LAYER (FROM BASE OF LANDFILL) |
^ ^ , -^ ~~~ ^ -i i n:
O~' ^~~^ drain *~^~~ ^
:
BARRIER SOIL LAYER
i
__ . -]
1
i
PERCOLATION
(FROM BASE OF LANDFILL)
Figure 2.1 Profile of a typical hazardous waste landfill.
j
r
-------
Ul
tc
CL
z~
o
H
a
u
UJ
cc
STORM
DURATION
TIME
INTERVAL BETWEEN
STORMS
EVENT DURATION (DEVENT)
INFIL = PRECIP-RO-ET
, INF1L
DEVENT
where
NFIL
PRECIP
flO
DEVENT
I
ET
Water available for infiltration during the evert (cm)
Storm precipitation (cm)
Runolt^luring4he event (cm)
Even! duration (d)
Infiltration rate (cnVd)
Evapotranspiration during the event (cm)
Figure 2.2 Depiction of a typical event duration.
-------
SOIL MOISTURE CONTENT
Wilting Point
Available
Soil Moisture
Vegetation
Uppermost
Landfill Layer
Second
Landfill Layer
PLAN VIEW
SECTION VIEW
Monitoring
-.Well Ground Surface
Figure 3.1 A schematic of the waste facility and leachate migration
through the unsaturated and saturated zones.
-------
Flow
1
"a
Layer 1
Layer 2
Layer 3
z=0
z=l1+!2
Figure 4.1 A schematic of transport through the layered unsaturated
-------
Figure 5.1(b) A schematic diagram of the Patch source boundary
condition for the saturated zone transport module.
-------
X-Z PLANE
Figure 5.1 (a) A schematic diagram of the Gaussian source boundary
condition for the saturated zone transport module.
-------
:^>m^^^^^
Figure 6.1(a) Groundwater contaminant plume
interception by the surface stream.
«$8SSWS»
Plume Boundary,
. \ s *
' N«r Fkld I I
''M»tnaZenc> I
MtBurcmrfil
Point
Groundwttcr
Lewling
Figure 6.Kb) Downstream contaminant transport from
the edge of the initial mixing zone.
-------
STREAM aOW
WASTE \
FAILURE I
TRANSPORT IN
GROUNDWATER*
M1XNGAT
ENTRY POINT
TRANSPORT IN
STREAM
DRINKM3 WATER
PLANT
HUMAN EXPOSURE VIA DRINKING
WATER CONSUMPTION
' Unsalurated and saturated
Figure 6. 2 (a)
Stages between failure of the waste containment facility
and human exposure due to drinking water.
-------
Unciunictf *nd awnM
Figure 6.2(b)
Stages between failure of the waste containment
facility and human exposure' due to fish consumption.
EXPOSURE TO
AOUA1IC OflGANISMS
Figure 6.2(c)
Stages between failure of the waste containment
facility and exposure to Aquatic Organisms.
-------
RECEPTOR
Figure 8.1 Wind direction sectors for Gaussian models.
-------
VIRTUAL POINT
SOURCE
AREA
SOURCE
Figure 8.2 Virtual source approximation for long-term Gaussian models.
-------
EFFECTIVE
SOURCE
AREA
SECTOR
CENTERLINE
RECEPTOR
Figure 8.3 Effective source area used by VALLEY.
-------
Model Parameters/Data
EPAMM Model
w = g(x)
Model Output
INPUT VALUES
INPUT VALUES
2
INPUT VALUES
2
E
INPUT VALUES
OUTPUT VALUES
INPUT DISTRIBUTIONS
OUTPUT DISTRIBUTION
Figure 9.1 A schematic description of the Monte Carlo method of
uncertainty analysis.
-------
c: 4
Region for Johnson
Distribuiion
for Johnson SL Disiribution
(lognormal)
Line for
student t
distribution
Region for Johnson
Sn Distribuiion
Distribution
Impossible Region
. SKEWNESS
Source: McGrath at al. 1973
Figure 9.2 Selecting a Johnson distribution from skewness and kurtosis.
-------
ft
re
I
O.
o
5
i
Figure 9.3 Comparison of the exact and the generated cumulative frequency
distribution for a normally distributed variable.
I"
-8
D
O
.9.
.8.
.7
.6
.5
.4
.3
.2
.1.
10 11
Values
12
13
-15
Figure 9.4 Comparison of the exact and the generated cumulative frequency
distribution for a log-normally distributed variable.
-------
XI
re
X»
o
t»
>
30
50
ED
70 BO
Values
Figure 9.5 Comparison of the exact and the generated cumulative frequency
distribution for an exponentially distributed variable.
XI
n
o
£
120
Figure 9.6 Comparison of the exact and the generated cumulative frequency
distribution for an empirically distributed variable.
-------
10
Figure 9.7 Comparison of the exact and the generated cumulative frequency
distribution for a uniformly distributed variable.
-------
-------
o
u.
O
LL)
O
ID
tr
u.
ID
5
s
ID
O
i.o
0.9-
0.8-
0.7-
0.6-
0.5-
0.4-
0.3-
0.2-
0.1-
0.0
I 1 \ I I I I I 1
0.0 0.1 02 0.3 OX 05 0.6 0.7 0.8 0.9 1.0
NORMALIZED CONCENTRATION. Cw*
* Normalized with respect to source concentration
Figure 9.8 Typical results obtained using MULTIMED in the
Monte Carlo mode.
-------
TABLE 9-2(b). RESULTS OF RANDOM NUMBER GENERATOR TEST FOR 1000 VALUES
Input Statistics
Observed Output Statistics
Normal
LogNormal
Exponential
Empirical*
Uniform
mean std. dev.
10.00 1.00
10.00 1.00
10.00 10.00
18.855
10** 25***
mean std. dev. max min
10.10 0.95 12.90 7.40
10.03 1.00 13.40 7.50
10.00 9.66 81.20 0.00
19.73 27.33 99.90 0.10
17.38 - 25.00 10.10
Cumulative Probability 0.0
Values 0.1
Expected Mean
**Minimum Value
***Maximum Value
0.1 0.7 1.0
1.0 10.0 100.0
18.855
-------
TABLE 9-2(a). RESULTS OF RANDOM NUMBER GENERATOR TEST FOR 500 VALUES
Input Statistics
Observed Output Statistics
Normal
LogNonnal
Exponential
Empirical*
Uniform
mean std. dev.
10.00 1.00
10.00 1.00
10.00 10.00
18.855
10** 25***
mean std. dev. max min
10.00 1.05 13.40 6.90
9.97 0.98 13.20 7.60
9.80 9.67 53.70 0.00
18.54 25.54 99.20 0.10
17.4 - 24.9 10.1
*Cumulative Probability 0.0
Values 0.1
Expected Mean
**Minimum Value
***Maximum Value
0.1 6.7
1.0 10.0
18.855
1.0
100.0
-------
------- |