Intro To Hydrology
Intro To Hydrology
Intro To Hydrology
MARGULIS
Introduction to
Hydrology
Cover photo credit: NASA Goddard Space Flight Center Image by Reto Stckli (land surface, shallow water, clouds). Enhancements by Robert Simmon (ocean
color, compositing, 3D globes, animation). Data and technical support: MODIS Land Group; MODIS Science Data Support Team; MODIS Atmosphere Group;
MODIS Ocean Group Additional data: USGS EROS Data Center (topography); USGS Terrestrial Remote Sensing Flagstaff Field Center (Antarctica); Defense Meteorological Satellite Program (city lights).
Throughout the text I have attempted to provide credit to the author/creator/reference from which any graphics
were obtained. Most of the graphics that are not original creations for this text were taken from the publiclyavailable COMET (MetEd) program. Others were taken from the public-domain and/or from credited sources. I
apologize for any omissions and would greatly appreciate being informed of them for correction in future versions.
Preface
This open-access e-textbook is the first iteration in an
experiment designed to provide a novel tool for learning the
basic concepts in hydrologic science. I hope it is a tool that
will be useful to the hydrologic community at large. It is
expected that the book and accompanying codes will be a
dynamic entity that will continue to evolve going forward. The
book is optimized for use as an electronic book on an iPad (or
other iBooks platform) to take full advantage of multimedia,
search capability, web links, etc. However, for those without
access to those platforms we also oer a PDF version (which
will not have the full functionality of the iBooks version).
Note: To be able to play the embedded multimedia
links in the PDF you will need to view the document
using the freely available Adobe Reader program
(www.adobe.com).
The material in the book is primarily derived from notes and
other material used in the undergraduate Civil &
Environmental Engineering (CEE) 150: Introduction to
Hydrology course oered in the Department of Civil and
Environmental Engineering at UCLA. CEE 150 is a 10-week
course designed to provide a survey of the hydrologic cycle.
The format of the course and topical coverage owes much of
its genesis to a similar course taught at MIT by Dara
Entekhabi. I have taught the course at UCLA for over 10
years and it has undergone a significant evolution over those
years.
Table of Contents
Preface ..........................................................i
Chapter 1: The Hydrologic Cycle .................7
Section 1: Learning Objectives ..........................................8
Section 2: Motivation ........................................................9
Section 3: The Hydrologic Cycle ......................................12
Section 4: Unique Properties of Water .............................14
Section 5: Mass balance, fluxes, and units .......................18
Section 6: Global Hydrologic Cycle and Average
Mass Balance.........................................................22
Section 7: Watershed Mass Balance .................................26
Section 8: MOD-WET Codes ...........................................32
Section 9: Conceptual Questions ......................................33
Section 10: Sample Problems ...........................................34
1:
2:
3:
4:
5:
6:
7:
8:
Section
Section
Section
Section
Section
Section
Section
Section
1:
2:
3:
4:
5:
6:
7:
8:
Model ................................................368
vi
Chapter 1
The
Hydrologic
Cycle
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
1. Provide a basic definition of hydrology
2. List some of the key motivations for studying hydrology
3. Describe the primary reservoirs in the global hydrologic
cycle and their relative sizes
4. Describe the primary water fluxes connecting the
reservoirs in the global hydrologic cycle
5. Describe and dierentiate the unique properties of water
relevant to hydrology
6. Convert between mass, volume, and energy fluxes and flux
densities
7. Write down and apply the mass balance equation for a
particular control volume
8. Identify the average rate of the global hydrologic cycle and
the corresponding average residence times of water in
dierent global hydrological reservoirs
9. Define a watershed and explain its relevance to hydrology
S ECTION 2
Motivation
Before beginning our discussion of hydrology we should
define what it is and why its study is important. Simply put,
hydrology is the study of water. In particular, hydrologic
science is a branch of geoscience that aims to diagnose and
predict: 1) the spatial and temporal characteristics of water in
its various storage reservoirs (terrestrial, atmospheric,
oceanic) and 2) the corresponding fluxes of water between
these reservoirs.
There are many motivations for studying hydrology. The
fundamental motivation is that water is required for life,
making it of paramount societal importance. While the overall
amount of water is fixed, population is increasing and
therefore having an understanding of water and its storage
and movement is key to our survival and will likely dictate
how and where we will live in the future.
There is also a significant scientific motivation for
studying water because it plays a key role in the Earth system
as a whole, including weather and climate processes,
landscape evolution, and biogeochemical processes. For
example, in terms of climate, we would like to know how and
why water varies from one location to another. A relevant
example in terms of annual average precipitation distribution
over California is shown in Figure 1.1. The map shows clear
gradients in precipitation from North to South as well as a
F IGURE 1.1 Map of annual average precipitation over California (from www.ocs.orst.edu/prism).
10
11
S ECTION 3
VOLUME
(km3)
% TOTAL WATER
1,322,000,000
97.2
29,199,700
2.1
Groundwater (near-surface)
4,171,400
0.31
130,700
0.017
Soil Moisture
66,700
0.005
Atmosphere
12,900
0.0009
feeds surface water bodies via lateral flow and runo. The
atmospheric water is replenished via evaporation from the soil
and open water surfaces and transpiration from vegetation.
Fluxes and storage are directly linked via mass balance as
described in more detail below.
A key aspect of the hydrologic cycle is the fact that it is
driven by energy inputs (primarily from the sun; Figure 1.3).
At the global scale, the system is essentially closed with
respect to water; negligible water is entering or leaving the
system. In other words, there is no external forcing in terms of
a water flux. Systems with no external forcing will generally
eventually come to an equilibrium state. So what makes the
hydrologic cycle so dynamic? The solar radiative energy input,
which is external to the system, drives the hydrologic cycle.
Averaged over the globe, 342 W m-2 of solar radiative energy
is being continuously input to the system at the top of the
atmosphere. This energy input must be dissipated, and this is
done, to a large extent, via the hydrologic cycle. Due to this
fact, the study of hydrology is not isolated to the study of
water storage and movement, but also must often include
study of energy storage and movements.
13
S ECTION 4
msduncanchem.com/Unit_11/phase_diagrams_ws_files/image001.gif).
SPECIFIC HEAT
TEMPERATURE
DENSITY
(C)
(kg m-3)
999.87
4216
15
999.13
4184
30
995.67
4177
CAPACITY
(J kg-1 K-1)
15
E XAMPLE 1.4.1
How much energy would be released if all of the
atmospheric water vapor were condensed?
According to Table 1.1, the volume of water in the
atmosphere is approximately equal to 12,900 km3. The
equivalent amount of mass (using a density of water of
1000 kg/m3) is then given by:
(1000)3 m 3
Mass = (12,900 km )(1000 kg m )
1 km 3
= 1.29 1016 kg
3
-3
E XAMPLE 1.4.2
If the energy released in the previous example
were absorbed by the global oceans, how much of
a temperature change would there be?
The specific heat capacity of a substance dictates how its
temperature will change (per unit mass) given an energy
input. According to Table 1.1, the volume of water in the
global oceans is approximately equal to 1,322,000,000
km3. The equivalent amount of mass (using a nominal
density of water of 1000 kg/m3) is then given by:
(1000)3 m 3
Mass = (1,322, 000, 000 km )(1000 kg m )
1 km 3
= 1.322 10 21 kg
3
-3
Temperature change =
(3.225 10 22 J)/(1.322 10 23 kg)/(4216 J kg -1K-1 )
= 5.8 10 5 K
Due to the large heat capacity of the ocean (and its large
mass) there would be an imperceptible change in
temperature due to this large energy input. This is
illustrative of the large buering capacity of water.
17
S ECTION 5
Properties of the fluid are denoted as either extensive or
intensive. Extensive properties are those that are in some way
related to the total mass of the system, including mass (m),
momentum (mV), and energy (E). Intensive properties are
generally normalized by mass (e.g., momentum per unit mass,
which is simply velocity, energy per unit mass, etc.). In the
control volume approach, we can represent an extensive
property as B and an intensive property as B (dB/dm). The
general control volume equation can be expressed as (Mays,
2005):
dB d
=
B dV + CS BV dA
dt
dt CV
(1.5.1)
dB
total rate of change of extensive property
dt
d
of change of extensive property in control
B dV rate
dt CV
volume (where is density of fluid and dV is a
differential volume)
rate of outflow of extensive property across
BV dA netcontrol
surface (where V is velocity and dA
CS
The mass balance (or continuity) equation is derived
from the Reynolds transport theorem by setting B to mass
(making B = 1) and noting that for mass conservation the
left-hand-side term equals zero, which yields:
d
dM
dV
=
= V dA
dt CV
dt
CS
(1.5.2)
18
d
dV = CS V dA
dt CV
(1.5.3)
kg 1
m3
=
;
s
w
s
dS
= I i Oi
dt
i
i
(1.5.4)
I = O
i
In the context of mass-balance (or other) applications,
several forms of fluxes are used and it is important to be able
to dierentiate between them and the units associated with
them. Mass fluxes (i.e. used in Equation (1.5.2)) have units of
mass/time (e.g., kg s-1). Volume fluxes (i.e. used in Equation
(1.5.4)) have units of volume/time (e.g., m3 s-1). To convert
between the two, the density of the fluid can be used, i.e. for
water:
(1.5.5)
w = 1000 kg m 3
kg 1
kg
2
s A
ms
m3 1
m
mm m
Volume:
(or
,
, etc.)
s
A
s
d
y
Mass:
19
As will be discussed in much more detail later, and has
been alluded to above, water and energy fluxes are often
directly connected. The most relevant example is that the
evaporation flux is tied directly to a phase change energy flux.
For example, the latent heat flux associated with phase
change is simply the mass flux density multiplied by the latent
heat of vaporization:
J kg
= Wm 2
2
kg m s
where the units of W m-2 are the most commonly used for
energy fluxes (actually flux densities) in hydrology.
E XAMPLE 1.5.1
A golf course has requested a permit to install a
40,000 ft2 pond to enhance the beauty of its
facilities. It is hypothesized that due to high
evaporation rates, the water in the pond will
have to be supplemented with pumped
groundwater. There is a small creek that
discharges an average of 0.0005 m3/s into the
pond. The outlet valve from the pond releases an
average rate of 0.0003 m3/s to keep the pond
from getting stagnant. Precipitation on the pond
is 260 mm/year, and the annual evaporation is
estimated to be 105 inches/year. The pond will
where the mass balance for the pond can be written as:
dS pond
dt
20
Apond
! 1m $
2
2
= 40000 ft #
& = 3718 m
" 3.28 ft %
mm 1 m
3718 m 2 = 967 m 3 /y
y 1000 mm
in 2.54 cm 1 m
E = 105
3718 m 2 = 9916 m 3 /y
y 1 in 100 cm
P = 260
Qout
m3
Qin = (0.0003 0.0005)
s
(3600 24 365) s
= 6307 m 3 /y
1y
21
S ECTION 6
dS
= Pland Eland + Pocean Eocean
dt
(1.6.1)
X =
1
T
t +T
X dt
(1.6.2)
equation:
dS
= Pland Eland + Pocean Eocean
dt
(1.6.3)
dS
0 = Pland Eland + Pocean Eocean = Pglobal E global
dt
(1.6.4)
Pglobal = E global
(1.6.5)
Pglobal 1 m yr -1 = E global
Residence Time =
Volume of reservoir
Average volumetric flux rate
(1.6.6)
VOLUME
(km3)
RESIDENCE
TIME
1,322,000,000
2500 years
29,199,700
Groundwater (near-surface)
4,171,400
8 years
130,700
88 days
Soil Moisture
66,700
47 days
Atmosphere
12,900
9 days
F IGURE 1.9 Illustrative photo of the arid terrain of the Atacama desert (from en.wikipedia.org/wiki/File:Paranal_360-degree_
Panoramic.jpg).
25
S ECTION 7
26
27
dS
= I i Oi = P E Q
dt
i
i
(1.7.1)
dS
0 = P E Q
dt
(1.7.2)
so that over the long-term the average basin yield (i.e. the
28
Q = P E
(1.7.3)
en.wikipedia.org/wiki/File:California_water_system.jpg).
dS
= P E +Q
dt
where the storage and fluxes (precipitation, evaporation,
and inflow respectively) can be expressed in terms of
dh
= P E +Q
dt
where the fluxes on the right-hand-side can be written in
terms of volume flux densities, which are just the volume
fluxes per unit area:
P = 1.0 m y -1
! 1 m $ ! 365 d $
-1
E = 1.1 mm d #
&#
& = 0.4 m yr
" 1000 mm % " 1 y %
!
$ ! 365 d $ ! 1 $
3
-1 86400 s
-1
Q = 1.75 m s #
&#
& # 8 & = 0.55 m yr
" 1 d % " 1 y % " 10 m %
-1
dh
= (1.0 0.4 + 0.55) m y -1 = 1.05 m y -1
dt
Based on this rate, with no outflows, the reservoir will
take 2.9 years to fill up to a depth of 3 meters.
If there is a reservoir release (Qout) the (steady-state)
30
dS
= P E +Q Qout = 0
dt
Qout = P E +Q = 1.05 m y -1
which, not surprisingly, is the same as the filling rate if
there is no outflow. The above flux can be converted
back to volumetric flux units by multiplying by the
reservoir area, which yields: Qout = 3.3 m3 s-1.
In practice, most reservoirs are primarily responsible for
smoothing out natural variability so that in wet years
the reservoir may add storage and in dry years storage
may be depleted. The size of the reservoir is a key design
variable that needs to be large enough to meet demands
or attenuate floods with a specified reliability criteria.
Note that in this simple example the reservoir area is
constant, whereas in most cases, the area will increase as
the reservoir fills.
31
S ECTION 8
MOD-WET Codes
Throughout the course, several of the equations and/or
concepts that can be applied numerically are implemented for
student use in MATLAB code. Together these functions are
referred to as the Modular Distributed Watershed
Educational Toolbox (MOD-WET). The MOD-WET codes
are available as a companion to the book and are provided to
help cement basic applied concepts for students. They are
generally implemented as functions (i.e. with inputs/outputs)
rather than scripts. This makes them modular in structure so
they can be used in conjunction with other MOD-WET codes
or additional codes you develop on your own. They are listed
at the end of each chapter where the relevant concepts are
introduced. The actual functions provide documentation
about the inputs/outputs/units and other details. The user
should carefully examine each code before applying it so that
you are aware of the details behind the calculations.
Relevant functions based on concepts introduced in this
chapter include:
Watershed delineation code:
watershed_area_and_stream_delineation.m
32
S ECTION 9
Conceptual Questions
33
S ECTION 10
Sample Problems
Problem 1.1. Based on this chapter and background
reading references (e.g the chapter Tapping into a Planetary
Cycle from the book Introduction to Water in California
[Carle, 2004]) answer the following questions:
a) Name Earths hydrologic reservoirs where the great water
wheel takes place? What powers the water cycle (i.e. what is
the energy input that drives the hydrologic cycle)? What
fraction of the Earths freshwater is considered active and
which hydrologic reservoirs are associated with this water?
What fluxes transport water through the hydrologic cycle and
connect one reservoir to another?
b) What is the average annual precipitation rate in Los
Angeles, CA? In Eureka, CA? What period of the year does
California receive most of its precipitation? Which direction
do winds typically blow in California? How are these winds
related to the dry region of the Mojave Desert?
c) Why is the Sierra Nevada snowpack considered a water
reservoir? What technique has been used by water supply
planners to forecast water availability during spring-summer?
d) What technique has been used by scientists to reconstruct
periods of drought or wet years when a historic rainfall record
is not available? Using this method, experts have been able to
e) Suppose you are making yourself 0.50 liters of tea but its
temperature is 75C, which is too hot for you to drink. In
order to cool its temperature, you immerse three ice cubes of
the same dimension as the one in part c) in your tea. If the
energy used to melt the ice came from the internal energy (i.e.
temperature) of the tea, what is the lowered temperature of
your tea?
f ) Suppose you do not drink all of the tea, and you leave 50
ml in the cup. How much energy would be required to fully
evaporate the remaining tea?If this energy came from the air
in a room surrounding the cup of tea, how would the air
temperature change, i.e. by what amount and via an increase
ordecrease in air temperature? Assume the room has
dimensions of 5 m x 5 m x 3 m and the air density is 1.2 kg/
m3.
Problem 1.4. The reservoir shown in cross-section below
receives inflow from a river with annual mean discharge Q =
2.55 m3s-1. Observations from a nearby precipitation station
YANGTZE
RIVER
(CHINA)
OB RIVER
(RUSSIA)
WATERSHED
AREA (km2)
7,180,000
1,970,000
2,950,000
PRECIPITATION
(mm/yr)
1,777
1,120
563
STREAMFLOW
(m3/s)
190,000
35,000
12,500
36
37
Chapter 2
Atmospheric
Thermodynamics
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
1. List the key constituents of the atmosphere, specifically
including the radiatively active components
2. Write down the ideal gas law (equation of state) for moist
air in terms of air temperature, density, pressure, and
water vapor content
3. Use the ideal gas law to compute the density of air from
the other states
4. Define the concepts of saturated vapor pressure and
saturated specific humidity and how they are dependent
on temperature (i.e. the Clausius-Clapeyron equation)
5. List and define the typical descriptors used to identify the
actual water vapor content of air and be able to convert
between the dierent metrics
6. Describe the dierence between specific humidity or vapor
pressure and saturated specific humidity or saturated vapor
pressure.
7. Use the hydrostatic equation to relate atmospheric
pressure and density as a function of height
39
S ECTION 2
Atmospheric Composition
The atmosphere is a gaseous mixture of several
constituents, some of which are relatively constant, while
others are quite variable. The composition of the atmosphere
is shown in Table 2.1. By volume, the air in the atmosphere is
composed almost entirely of nitrogen and oxygen (together
equaling 99%). These two constituents and many others
(argon, neon, helium, krypton, hydrogen, and nitrous oxide)
represent the vast majority of the atmosphere and are
relatively constant in their concentration. Several components
vary as a result of natural or anthropogenic processes. The
most variable constituent is water vapor, which varies
naturally as a result of the hydrologic cycle. Carbon dioxide
varies both naturally (due to photosynthesis and respiration)
and anthropogenically (due to the burning of fossil fuels). Of
particular importance to climate and the hydrologic cycle are
those gases that are radiatively active (i.e. water vapor,
ozone, methane, carbon dioxide, and to a lesser extent
oxygen).
A radiatively active gas is one whose molecular structure
is such that it absorbs (and emits) radiative energy (this will
be discussed in much more detail in the next chapter).
Radiatively active gases are often referred to as greenhouse
gases. The primary hypothesis for anthropogenic climate
change has to do with the combined facts that 1) humans
% BY VOLUME
Nitrogen (N2)
78.084
Oxygen (O2)*
20.946
Argon (A)
0.934
0.033
Neon (Ne)
0.00182
Helium (He)
0.000524
Methane (CH4)*
0.00016
Krypton (Kr)
0.0014
Hydrogen (H2)
0.00005
0.000035
Ozone (O3)*
Water (H2O)*
0 - 0.000007
0-4
S ECTION 3
Atmospheric States
The state of the atmosphere is the set of variables that
can be used to quantitatively describe it. The thermodynamic
states include temperature, pressure, density, and some metric
of water (vapor) concentration. It is important to keep in
mind that the atmosphere is simply a moving fluid (with the
above characteristics varying in space and time) so another
key state variable is the wind speed (which has components in
three orthogonal directions). Knowledge of the state of the
atmosphere at a given location/time will ultimately provide a
mechanism for quantifying many of the water and energy
movements in the atmosphere and between the atmosphere
and land surface. Some key questions that will be addressed in
this and other sections of this chapter include:
1. How are the variables that comprise the thermodynamic
state of the atmosphere related (at a given location)?
2. What metrics for water vapor concentration are used, and
related to some of these metrics, what is the upper limit to
the concentration of water vapor?
3. How do key variables vary with height and what
relationships govern this variation?
To start with the first question, it is worth pointing out
that the atmosphere behaves like an ideal gas (something
pi = i RiT
(2.3.1)
= d + v
(2.3.2)
42
pd
p
p
e
+ v = d +
RdT RvT RdT RvT
(2.3.3)
where pd and pv are the partial pressures of dry air and vapor
respectively, and Rd and Rv are the dry air gas constant (287 J
kg-1 K-1) and the vapor gas constant (461 J kg-1 K-1)
respectively. In many texts related to hydrology, the partial
pressure for water vapor (i.e. vapor pressure) pv is replaced
by e, which is what we will do here. Equation (2.3.3) can be
rearranged noting that the total pressure is equal to the sum
of the two partial pressures (i.e. p = pd + e) to form a more
compact equation of state (i.e. ideal gas law) for moist air:
Rd
p "
e%
1
(1
)
;
=
= 0.622
$
'
RdT #
p&
Rv
(2.3.4)
p = RdTv ; where Tv =
e
1
(1
(2.3.5)
Tv =
"
e%
$1 (1 ) '
p&
#
= 300.99 K
(25 + 273.15) K
"
2500 Pa %
$1 (1 0.622)
'
100000 Pa &
#
43
p
100000 Pa
=
= 1.16 kg m -3
-1
-1
RdTv (287 J kg K )(300.99 K)
44
S ECTION 4
w=
v
e
=
d
p e
(2.4.1)
q=
v
e
=
v + d
p (1 )e
(2.4.2)
q w
e
p
(2.4.3)
Typical values for specific humidity (or mixing ratio) near the
surface are 5-1510-3 kg kg-1, where the small values simply
represent the small fraction of vapor in the atmosphere (Table
2.1). Because of the small values, units of g/kg are often used,
in which case the typical values near the surface are 5-15 g
H2O/kg Air.
Several other metrics exist that are relative in nature.
Before defining these, the upper limit (water vapor holding
capacity) of the atmosphere must be discussed. From
thermodynamic principles one can define the saturated vapor
45
bottle with some air space above the liquid water. Given the
temperature of the system in the bottle, the equilibrium
amount of vapor in the air space would have a vapor pressure
equal to the saturated vapor pressure. In truth, some water
molecules will escape the liquid state and others will leave the
vapor state and condense on the liquid. The saturated vapor
pressure is the equilibrium condition where the rate of
vaporization of liquid water molecules is equal to the rate of
condensation of vapor molecules. An analogous quantity can
be defined over bulk ice and yields:
des Lv es
=
dT Rv T 2
L
ei (T ) = es 0 exp s
Rv
(2.4.4)
1 1
T0 T
(2.4.6)
L
es (T ) = es 0 exp v
Rv
In practice the saturated vapor pressure is to a very
good approximation the upper limit to the actual vapor
pressure in air, i.e.:
1 1
T0 T
(2.4.5)
0 e es (T )
(2.4.7)
qs
e ; with 0 q q s (T, p)
p s
(2.4.8)
46
L
e = es (T = Td ) = es 0 exp v
Rv
1
1
T T
0
d
(2.4.10)
These newly defined saturation quantities allow for the
definition of other metrics of vapor concentration in air:
1 Rv e
Td =
ln
T
L
es 0
0
v
RH =
e
q
=
es q s
(2.4.9)
Td T
(2.4.11)
(2.4.12)
VPD = e = es e
(2.4.13)
T Tw
L
= v
q s (p,Tw ) q c p
(2.4.14)
E XAMPLE 2.4.1
For a parcel of air with a vapor pressure of 1500
Pa, a temperature of 20C, and a pressure of
98000 Pa, determine the following alternate
humidity metrics for the parcel:
a) vapor density, b) specific humidity, c) relative
humidity, d) dew point temperature, and e)
vapor pressure deficit
a) The vapor density is given by the ideal gas law for
water vapor (Equation (2.3.1)), i.e.:
v =
e
1500 Pa
=
RvT (461 J/kg/K)(293.15 K)
q =
=
e
p (1 )e
1500 Pa
98000 Pa (0.378)(1500 Pa)
48
) 2.5 10 6 J/kg #
&,
1
1
es = (611 Pa)exp +
%
(.
461
J/kg/K
273.16
K
293.15
K
$
'*
= 2368 Pa
where the larger the vapor pressure deficit the drier the
air is (i.e. relative to saturation). Air with a relative
humidity of 100% will have a vapor pressure deficit equal
to zero.
)
# 1500 Pa &,
1
461 J/kg/K
Td = +
ln
%
(.
6
* 273.16 K 2.5 10 J/kg $ 611 Pa '= 286.1 K = 13. C
49
S ECTION 5
T(z) = T0 z,
T0 surface temperature
(2.5.1)
dp
= g
dz
(2.5.2)
F IGURE 2.3 Typical profiles of specific humidity in the troposphere at different latitudes.
51
To characterize the bulk amount of vapor in the
atmospheric column, we often use the so-called total
precipitable water, which is simply the integration of water
vapor mass (per unit volume, i.e. vapor density) with height:
Wp =
dz
(2.5.3)
Wp =
1
v dz = q dz =
g
0
1
q dp = g
ps
ps
q dp
(2.5.4)
Wp = 1.12 exp(0.0614Td )
(2.5.5)
E XAMPLE 2.5.1
Estimate the precipitable water from the tropical
specific humidity profile in Figure 2.3. Assume a
standard atmospheric pressure profile. How does
that estimate compare to the one from the
empirical expression shown in Equation (2.5.5)?
Without data for the specific profile, it can be
approximated via a piecewise linear profile where at
altitudes of 0, 3, 5, 7, and 9 km, the corresponding
specific humidities are approximately given by 16, 5.5,
2.5, 1.5, and 0.5 g/kg respectively. The atmospheric
pressures corresponding to each level can be roughly
estimated from the standard atmosphere shown in Figure
2.2 as 101300, 70120, 54050, 41110, and 30800 Pa
respectively for the altitudes listed above.
The integral in Equation (2.5.4) can then be evaluated
numerically (i.e. using the trapezoidal rule):
1
Wp =
g
ps
q dp g 2 (q
p
i+1
+ qi )(pi+1 pi )
An illustrative example of the variability in total
precipitable water on the timescale of days is shown in Movie
2.1. The key features are that the highest precipitable water
occurs over equatorial regions (where oceans are prevalent and
0.622
The dewpoint temperature is then given by:
1
)
# 2524 Pa &,
1
461 J/kg/K
Td = +
ln %
(.
6
273.16
K
611
Pa
2.5
10
J/kg
$
'*
= 294.2 K = 21. C
evolution (over the ocean) as measured from the AMSU satellite (from sos.noaa.gov/videos/PrecipitableWater.mov).
53
M OVIE 2.2 Illustration of monthly averaged precipitable water over the globe from July 2002-July 2012 (from
earthobservatory.nasa.gov/GlobalMaps/).
54
S ECTION 6
MOD-WET Codes
Relevant functions based on concepts introduced in
this chapter include:
Air density:
air_density.m
Bisection (root finding) function:
bisect.m
Vapor pressure (from dew point temperature):
dew_point_temperature_to_vp.m
Script with definitions of key physical constants used in
other functions:
physical_constants.m
55
S ECTION 7
Conceptual Questions
1. Name at least two radiatively active gases in the
atmosphere.
2. What gas in the atmosphere is most variable?
3. What is the general trend of carbon dioxide in the
atmosphere. What are its two key modes of variability?
Explain.
4. Name the atmospheric state variables that are related to
each other via the equation of state for moist air (i.e. ideal
gas law).
5. With all else equal, is an air parcel with more vapor
concentration less dense or more dense than a parcel with
less vapor?
S ECTION 8
Sample Problems
Problem 2.1. Several atmospheric constituents, most notably
water vapor, greatly aect many of the processes we will
study. Therefore, it is important to be able to estimate how
much water is in the atmosphere over an area of interest.
Additionally, it is important to become familiar with various
temperature and humidity relationships. Consider the
following:
a) A parcel of moist air is at 870 mb, with a temperature of
12C and a relative humidity of 95%, what is the moist
density of the air in kg m-3?
b) What is the saturated vapor pressure of air (in Pa) at a
temperature of 37C? Explain the physical meaning of
saturated vapor pressure. How does it change as the
temperature increases? How does it change as the humidity
level of air changes?
c) What is the specific humidity (in g/kg) of air in part a)?
d) What temperature would a dry parcel of air need to have
in order to have the same density of the air in part a)?
e) What is the dew point temperature (in C) of the air in
part a)? What is the meaning of the dew point temperature?
that takes in air from the closed room, removes water vapor
by condensing the vapor into liquid, and then discharges the
dry air back into the closed room. You run the dehumidifier
until the specific humidity is lowered to the required value.
a) What mass of water (in kilograms) will need to be
condensed out to reach the required humidity level?
b) What is the amount of energy released (in Joules) as a
result of condensing the humidity?
c) Assuming the energy computed in part b) goes into
warming the air, what is the expected rise in air temperature
due to the dehumidification. Note: The specific heat capacity
of air is 1004 J/kg/K. For your calculation, you may assume
that the change in air density (or change in air mass) is
negligible as a result of the heating/dehumidification.
d) Will the air temperature still be in the necessary range as
a result of the dehumidification or will it need to be
additionally cooled?
Problem 2.4. In winter-time, heated indoor air reduces the
relative humidity inside buildings. This is especially a problem
inside hospitals in cold climates where the air temperature
and relative humidity must be kept at 69F (20.6C) and 75%
respectively for health and hygiene reasons.
Consider a cold winter day in Chicago when the outside air
temperature is 20F (-6.7C) and the relative humidity is 90%.
(Assume air density is 1.2 kg m-3 and air pressure is 1000 mb).
# g
p(z) = p0 exp %
$ Rd
dz &
T (z)(
'
0
v
b) Derive the expression for the pressure profile for the case
with a linear virtual temperature profile:
Tv (z) = Tv 0 z
which is a reasonable approximation to the profile in the bulk
troposphere (i.e. with a constant temperature lapse rate).
Problem 2.6. Several atmospheric constituents, most notably
water vapor not only impact water fluxes, but greatly aect
the radiative properties of the atmosphere. Therefore, it is
important to be able to estimate how much water is in the
atmospheric column over an area of interest. The
58
T = 5 K km -1
" z%
v (z) = v 0 exp $ '
# H&
Wp =
(z)dz
0
59
Chapter 3
Radiation
Processes
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
S ECTION 2
Basics of Radiation
Radiation is a form of energy that plays a crucial role in
hydrology. Radiative energy provided by the sun is the
external forcing of the system that drives the global
hydrologic cycle and atmospheric and oceanic motions. The
Earths surface and atmosphere reflect and absorb energy
from the sun as well as emit radiative energy of their own.
The absorbed energy and energy fluxes tends to drive water
fluxes and transformations (i.e. phase changes). Due to the
importance of radiative fluxes we would like to understand the
mechanisms responsible for the variability in radiative forcing
and develop simple models for these fluxes. We will start with
a quick refresher on basic radiation physics. For a much more
detailed treatment of radiation processes, the reader is
referred to Liou (1992, 2002).
Radiation is a form of energy carried by electromagnetic
(EM) waves (or photons). These waves travel at the speed of
light and can be characterized by either wavelength (the
length of a single period of the wave) or frequency (the
number of periods of the wave passing a given point per unit
time). The relationship between the two is given by:
= c
wavelength [m]
c speed of light [3.0 10 8 m/s]
(3.2.1)
frequency [s 1 = Hz]
62
R = B =
2 hc 2
(3.2.2)
5 ehc/(K T ) 1
R = B d
(3.2.3)
R = T 4 [Wm 2 ]
(3.2.4)
R = B ;
spectral emissivity []
(3.2.5)
63
R = T 4
broadband emissivity []
(3.2.6)
64
S ECTION 3
Shortwave: Rs =
=4 m
=0.1 m
R d
when you see something, what your eye is really sensing is the
solar (visible) radiation being reflected from that object. The
objects color corresponds to the wavelengths of the visible
light that are reflected (rather than absorbed).
The terrestrial spectrum covers two dierent parts of the
EM spectrum: infrared (thermal infrared and far infrared) and
microwave (in the tail of the distribution beyond 50 microns).
The peak occurs in the thermal infrared (around 12 microns).
The basic principle of night-vision goggles takes advantage of
this fact where the sensors are tuned to be sensitive to the
thermal infrared region. So even though little visible light may
be available, the thermal infrared signature is still apparent.
65
Due to these dierences, we often treat the two fluxes
separately, in terms of both measurement and modeling. The
solar fluxes are referred to as shortwave radiation (spanning
approximately 0.1-4 microns) and terrestrial fluxes are
referred to as longwave radiation (spanning approximately
4-100 microns). Formally, the integrated fluxes can be written
as:
Shortwave: Rs =
Longwave: Rl =
=4 m
=0.1 m
=100 m
=4 m
R d
(3.3.1)
R d
(3.3.2)
66
S ECTION 4
Radiative Properties of
Media
As shown above, an emitted radiative flux from a
particular body depends on its physical temperature and its
emissivity. The broadband emissivity strongly depends on the
media characteristics. For land surfaces, the surface type
largely dictates the broadband emissivity (see Table 3.1) and
typically ranges between 0.9-1.0. The broadband emissivity of
the atmosphere strongly depends on the concentration of
radiatively active gases, i.e. conceptually:
EMISSIVITY
Water
0.92-0.97
Snow (old)
0.82-0.89
Snow (fresh)
0.90-0.99
Ice
0.92-0.97
Bare soil
0.84-0.97
Grass (long, 1 m)
0.90
0.95
Agricultural crops
0.90-0.99
Forests
0.97-0.99
t + a + r = 1
(3.4.1)
ts + as + rs = 1
(3.4.2)
tl + al + rl = 1
the sun and earth and b) absorption as a function of wavelength and molecular absorption (from homepages.ius.edu/
kforinas/ClassRefs/ElectroMagnetic.html).
tl + al 1 tl 1 al
(3.4.4)
as + rs 1 as 1
(3.4.7)
al 1
(3.4.8)
surface albedo
The same conservation principle can be applied to the
land surface for shortwave and longwave fluxes:
ts + as + rs = 1
(3.4.5)
tl + al + rl = 1
(3.4.6)
69
S ECTION 5
on Earths TOA relative to the incoming solar beam and the resulting solar zenith angle.
4
Tsun
(5.67 10 8 W m 2 K4 )(5780 K)4 = 6.33 10 7 W m 2
I 0 = 1367 W m 2
(i.e. the solar constant -- here scaled by the ratio of the intercepted area to the global surface area) is distributed over larger
areas as one moves further from the location normal to the incoming beam.
M OVIE 3.1 Animation showing the minute-by-minute evolution of the day-night terminator for the winter solstice (from
sos.noaa.gov/videos/oneday.mov).
72
Figure 3.9 provides an illustration of the diurnal and
seasonal evolution of the solar position for a given location (in
this case 40 North latitude). The details of how this position
is calculated is given below; here only the qualitative aspects
are discussed. The sun is generally lowest in the sky at the
winter solstice: only 26.5 above the south horizon (i.e. a solar
zenith angle of 63.5) at solar noon. Solar noon is simply the
time corresponding to the highest point in the sky on a given
day (and is generally dierent that local noon which is
impacted by time zones, daylight saving time, etc.). At this
location the sun is 73.5 above the south horizon (i.e. a solar
zenith angle of 16.5) at solar noon on the summer solstice.
For the Spring and Fall equinoxes, the sun is 50 above the
south horizon (i.e. a solar zenith angle of 40) at solar noon
since there is no relative tilt between the Eartha axis and the
axis of the orbital plane on those days.
Based on the above factors, the amount of incoming
shortwave radiation at any location/time at the TOA is a
function of Earth-Sun geometry which is completely defined
by: i) Latitude (i.e. location), ii) Hour of day (due to rotation
of the Earth), and iii) Day of year (due to tilted axis and
elliptical orbit around Sun). Several models for the TOA flux
based on these inputs are available at varying levels of
precision. A relatively simple model for the flux at the
location perpendicular to the TOA is:
cos 0
I0
, daytime: 0 90
2
Rs 0 =
d
0,
nightime
(3.5.1)
F IGURE 3.9 Schematic showing the seasonal evolution of solar position in the sky at a latitude of 40 degrees North.
(3.5.2)
# 2
&
23.45
declination angle =
cos %
(172 DOY )(
180
$ 365
'
latitude
T 12
hour angle = 2 h
24
# 2
&
d = 1 + 0.017 cos %
(186 DOY )(
$ 365
'
73
F IGURE 3.10 Seasonal evolution of the solar declination angle (from physicalgeography.net/fundamentals/images/declination.jpg).
" 2
%
23.45
cos $
(172 172)' = 0.409 rad
180
# 365
&
75
= 2
12 12
= 0 rad
24
" 2
%
d = 1 + 0.017 cos $
(186 172)' = 1.0165
# 365
&
The cosine of the solar zenith angle is then given by
(using radians for all angles, where the latitude is equal
to 0.593 radians):
cos 0 = sin(0.409)sin(0.593) +
cos(0.409)cos(0.593)cos(0)
Rs 0 = (1367 W m -2 )
0.983
-2
=
1300.5
W
m
(1.0165)2
" 2
%
23.45
cos $
(172 355)' = 0.409 rad
180
# 365
&
12 12
= 0 rad
24
" 2
%
d = 1 + 0.017 cos $
(186 355)' = 0.9835 rad
# 365
&
= 2
cos 0 = sin(0.409)sin(0.593) +
cos(0.409)cos(0.593)cos(0)
Rs 0 = (1367 W m -2 )
0.539
-2
=
761.7
W
m
(0.9835)2
77
S ECTION 6
F IGURE 3.13 Schematic showing how solar direct beam radiation is absorbed and scattered in the atmosphere (from
powerfromthesun.net/Book/chapter02/chapter02.html).
78
Rs = Rs,dir
+ Rs,dif
+ Rs,bs
Rs,dir
downwelling direct beam flux
Rs,dif
downwelling diffuse flux
Rs,bs
backscattered flux
(3.6.1)
F IGURE 3.14 Photograph showing diffuse radiation being reflected off of clouds near sunset.
In terms of estimating these fluxes, models of varying
complexity can be used. Of high complexity are those models
generally referred to as radiative transfer codes. These require
a full accounting of the profiles of atmospheric constituents
and surface conditions and attempt to model the absorption,
scattering and transmission as radiation interacts with each
discretized layer of the atmosphere. These types of models are
commonly used as modules within global climate models,
where the atmospheric state profiles are explicitly modeled at
every time step. Measurements of such profiles are not
generally available with sucient sampling density in space/
79
time and so for oine hydrologic analysis, simpler semiempirical bulk models are often used instead. A general
example of such a bulk model is given by:
Rs,dir
= tsRs 0
s,dif
= Rs 0
(3.6.3)
Rs,bs
= (Rs,dir
+ Rs,dif
)
(3.6.4)
Rs = ts + + ts + 2 Rs 0
(3.6.5)
SURFACE TYPE
ALBEDO
0.03-0.10
0.10-1.0
Snow (old)
0.40-0.70
Snow (fresh)
0.45-0.95
Ice
0.20-0.45
Bare soil
0.05-0.40
0.16
0.26
Agricultural crops
0.10-0.25
Forests
0.05-0.20
ts = sa dust
(3.6.6)
80
The model shown above is applicable to clear sky
conditions. To take into account clouds, an additional
attenuation must be included. For the purposes of generality,
Equation (3.6.5) could be rewritten as:
(3.6.7)
(3.6.9)
as = 0.0363 0.0084Wp
bs = 0.0572 0.0173Wp
M opt = 1 / cos
Rs = fsc (cloud)"#ts + + ts + 2 $% Rs 0
(3.6.8)
fsc (C ) = 1 0.65C 2
(3.6.10)
fsc (C ) = 1 (1 ksc )C
(3.6.11)
where,
(3.6.12)
E XAMPLE 3.6.1
The slope (S) is the angle in the vertical plane between
the tangent plane in the steepest descent direction and a
horizontal plane. Based on this definition, a horizontal surface
has a slope of 0 degrees while a vertical surface has a slope of
90 degrees. The aspect (A) angle is the azimuth angle in the
horizontal plane between the steepest descent direction and an
arbitrarily defined reference direction (usually North). So a
point that is facing due North, East, South, or West is
generally labeled with an aspect of 0, 90, 180, or 270 degrees
respectively. Based on the slope and aspect angles, one can
define the local illumination (zenith) angle, which is simply
the zenith angle of the solar direct beam relative to the sloped
surface:
(3.6.13)
F IGURE 3.15 Schematic showing how slope/aspect can impact incident shortwave radiation.
In most applications, the slope and aspect are generally
computed from gridded DEM data. Several methods can be
used to compute the slope and aspect angles. The most
commonly applied computes the slope at a point by fitting a
plane to the neighboring nine cells (inclusive of the point of
83
dz dz
tan(S) = +
dx dy
dz dx
tan(A) =
dz dy
(3.6.14)
(3.6.15)
84
F IGURE 3.18 Distribution of aspect (in degrees) of the terrain represented by the DEM in Figure 3.16.
85
To quantify the impact of topography and shading, the
modified geometry can be accounted for to correct the
horizontal plane fluxes. Muller and Scherer (2005) formulated
the slope/aspect and shade eects into a single correction
term for direct beam flux received on a given unit area as:
s,dir
"
cos s
1 %
= tsRs 0 $maskshade
'
cos
cos(S)
#
&
0
(3.6.16)
86
Rs,dif
= Rs 0SVF + Rs (1 SVF)
(3.6.17)
s,dir
" cos(20.6 )
1 %
= (0.6)(1300.5 W m )$(1)
'
87
s,dir
" cos(0.6 )
1 %
= (0.6)(1300.5 W m )$(1)
'
cos(10.6
)
cos(10
)
#
&
-2
88
(3.7.2)
Rl = aTa4
S ECTION 7
(3.7.1)
While the atmospheric emissivity is dependent on all
greenhouse gas concentrations, the most variable is water
vapor, so several empirical models have been developed to
estimate the emissivity based on humidity measurements.
Examples include that by Idso (1981):
a = 0.74 + 0.0049ea
(3.7.3)
by Satterlund (1979):
(T 2016)
a = 1.08(1 exp(ea a
));
[ea ] = mb;
[Ta ] = K
(3.7.4)
by Brutsaert (1975):
0.14
!e $
a = 1.24 ## a &&
"Ta %
[ea ] = mb;
[Ta ] = K
(3.7.5)
a = 0.605 + 0.048ea0.5 ;
[ea ] = mb
(3.7.6)
The above model is for clear-sky downwelling radiation.
To account for clouds, a similar approach to that used in the
shortwave case can be used, i.e.:
Rl = flcaTa4
(3.7.7)
flc = (1 + 0.22C 2 )
(3.7.8)
Rl = sTs4
(3.7.10)
(3.7.9)
Ts surface temperature
90
E XAMPLE 3.7.1
! 15 mb $
a = 1.24 #
&
" 293.15 K %
= 0.81
91
92
(3.8.2)
S ECTION 8
Rl Rl
The previous two sections provided physical descriptions
and relatively simple models for the key shortwave and
longwave fluxes occurring at the land surface. While it is
convenient to model each flux individually based on the
dierent physical relationships, the key radiative input is the
net (or integrated) radiation absorbed by the surface. It is this
net flux that will drive hydrologic processes at the surface.
The net radiation can be written simply as:
Rs Rs = Rs Rs = Rs(1 )
(3.8.1)
Putting the two together we can write the net radiation
at the surface as:
Rn (x,y,t) = Rs(1 ) + Rl Rl
(3.8.3)
S ECTION 9
MOD-WET Codes
Relevant functions based on concepts introduced in
this chapter include:
Calculation of diuse shortwave scattering coecient:
diffuse_sw_scattering_coefficient.m
Calculation of direct beam shortwave transmissivity:
direct_sw_transmissivity.m
Computation of incoming solar radiation over complex
terrain:
disaggregate_SW.m
Conversion between easting/northing UTM vectors and
lat/lon vectors:
easting_northing_to_lat_lon.m
Slope and aspect determination for a given DEM:
generate_slope_and_aspect_from_DEM.m
Idso model for atmospheric emissivity:
idso_emiss.m
Satterlund model for atmospheric emissivity:
satterlund_emiss.m
compute_shade_lookup_table_and_SVF.m
96
S ECTION 10
Conceptual Questions
1. Name the respective sources of shortwave and longwave
radiation fluxes. What part of the spectrum is covered by
each? At what wavelength does the peak for each occur?
2. According to the Planck function, name the physical
variable controls the magnitude (and distribution) of
emitted (blackbody) radiation. Name the additional
physical parameter that describes how close an actual body
is to a blackbody.
3. Name at least two atmospheric radiatively active gases that
absorb solar radiation. Name at least two atmospheric
radiatively active gases that absorb longwave radiation.
4. Define in words what surface albedo represents. What land
surface types generally have the smallest and largest
broadband albedo and what are their typical values?
5. Name the variables are needed to determine the TOA
incoming shortwave flux.
6. Define in words what the declination angle represents.
7. What is the maximum possible value of incoming shortwave
radiation at the TOA? Where on the globe will this
generally occur?
S ECTION 11
Sample Problems
Problem 3.1. In this problem, you are asked to examine the
influence of the atmosphere on the radiative flux incident at
the earths surface. Incoming solar radiation in the shortwavelength spectra (shortwave) at the top of the atmosphere
(TOA) is partially absorbed and scattered (reflected) by
atmospheric molecules and particulates. As a result, only a
fraction of the TOA flux actually reaches the surface. Answer
the following questions using the equations provided in
Sections 5 and 6.
a) Compute declination angle, zenith angle, and incoming
top-of-atmosphere (TOA) solar radiative flux for solar noon
on December 21st, 2011 (winter solstice), and June 21st, 2012
(summer solstice) in Hilo, Hawaii (19.73 deg. N). For this
latitude and days of the year (and time of day), is the sun in
the southern sky, northern sky, or directly overhead? [Hint:
You can answer this by comparing the latitude to the
declination angle]. Is your answer generally true throughout
the year? Justify your answer.
b) What is the predicted shortwave direct beam
transmissivity of the atmosphere on June 21st based on the
precipitable water value you estimated in problem #1 part c)?
Assume attenuation due to dust Based on your answer,
approximately what fraction of the incoming direct beam solar
For each of the three times listed above, answer the following
questions. Note: In reality the atmospheric parameters will
vary throughout the day due to the solar flux passing through
more or less atmosphere. Below we will assume they are
constant for simplicity.
b) Assuming that a transmissivity of 0.52 is representative for
this day over the Sierras, what is the magnitude of the direct
beam flux incident on a horizontal plane at the surface?
c) Assuming a scattering coecient of 0.185, what is the
magnitude of the diuse flux incident on a horizontal plane at
the surface?
d) Assuming the average surface albedo is 0.25 in this region,
what is the magnitude of the backscattered shortwave flux
incident on a horizontal plane at the surface?
e) What is the total shortwave radiation incident on a
horizontal plane at the surface? How do the three components
compare in magnitude?
f ) What is the net shortwave radiation absorbed by the
surface?
g) Describe conceptually what locations (latitudes) on the
globe would have increased solar radiation on south facing
slopes vs. those that would have increased solar radiation on
north facing slopes. [Hint: This is related to the declination
angle.]
= 0.3 0.1s
Simultaneous satellite measurements indicate that the areaaveraged outgoing longwave radiation from the entire region is
536 W m-2. A ground based measurement indicates that the
dry area surface temperature is 315 K. Incoming solar
radiation at this time is equal to 900 W m-2 and average air
temperature and relative humidity over the entire region is
312 K and 35% respectively. Assume a representative value of
0.97 can be used for the land surface emissivity in this area.
a) What is the surface temperature of the moist area?
b) What is the net radiation over the dry and moist areas?
Explain why they are dierent.
c) What is the area-averaged net radiation?
Problem 3.7. To a large degree, the climate of Earth is
controlled by radiative fluxes. Radiative equilibrium consists
of the balance of net solar (shortwave) radiation, which is an
external input to the system, and net terrestrial (longwave)
fluxes, which generally depend on the Earths surface and its
atmosphere. A schematic picture of the global radiative energy
balance is shown in the figure below. Radiative equilibrium
refers to the condition where the net radiative energy is equal
100
102
Chapter 4
Atmospheric
Circulation
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
1. Sketch a conceptual picture of the annual average net
radiation absorbed at the TOA by the Earth system
2. Describe how and why the net radiation distribution
causes atmospheric motions
3. Sketch a conceptual picture of the general circulation of
the atmosphere
4. List the key characteristics of the general circulation and
discuss its impact on climate patterns (including
precipitation patterns)
5. Identify how high/low pressure cells are connected to
circulation and vertical motion
6. Identify which circulation feature is primarily responsible
for equator to mid-latitude energy transport
7. Identify which circulation feature is primarily responsible
for mid-latitude to pole energy transport
104
S ECTION 2
Rn 0 (x,y) = Rs 0 (1 0 ) Rl0
(4.2.1)
0 planetary albedo
M OVIE 4.1 Illustration of monthly evolution of TOA net radiation (from earthobservatory.nasa.gov/GlobalMaps/
view.php?d1=CERES_NETFLUX_M).
" R2 %
2
E
'
1367 Wm $$
=
342
W
m
2 '
# 4 RE &
2
106
from global net radiation distribution along with estimated contribution by the atmosphere and ocean circulation. PW = 1015
W (used with permission from Marshall and Plumb, 2007).
The heat flux necessary to yield equilibrium can be
inferred from the imbalance in net radiation. An estimate of
107
S ECTION 3
Atmospheric Motions
Driven by Latitudinal Energy Imbalance
Given the implied heat transport discussed in the
previous section, the primary question is how does the
atmospheric transport work? Here we will focus primarily on a
qualitative description of the key characteristics of the
atmospheric motion, i.e. the so-called general circulation. It
is important to note that the general circulation we will be
discussing is representative of the long-term average. Hence,
the general circulation is responsible for much of the
climatology on Earth, while individual weather events
(including extremes) can be quite dierent than the mean.
The instantaneous state and motion of the atmosphere can
dier markedly from the general circulation.
The primary impact of the net radiation distribution
shown in Figure 4.2 is that, within the troposphere, warmer
air exists in the tropics relative to the poles (Figure 4.4).
Thermal expansion of the tropical air columns (relative to
those at the poles) leads to a sloping pressure surface aloft
(also seen in the temperature isotherms) which would be
expected to drive air toward the poles aloft. Additionally, air
near the surface in the tropics would be expected to be
warmest, which due to buoyancy would generate an uplift
temperature as a function of height. Red arrows indicate expected vertical motions at equator and poles.
(warm air rises) with cold air aloft at the top of the
troposphere generating sinking motion at the poles. If these
were the only factors, they would imply a single circulation
cell as shown in Figure 4.5, with surface flows from pole to
equator, upward flow at the equator, return flow aloft, and
downward flow near the poles. This simple convection cell
model would have the desired trait of transporting warm air
toward the pole and cool air to the equator, eectively
redistributing the net radiation imbalance. In fact, early
conceptual models of the Earths atmospheric flow proposed
this model (Marshall and Plumb, 2007). However, an
108
F IGURE 4.5 Schematic of expected single-cell thermallydriven atmospheric general circulation on a non-rotating
planet. (from http://www.ux1.eiu.edu/~cfjps/1400/circulation.html).
109
F IGURE 4.8 Mean annual surface pressure distribution showing areas of typically high and low pressure (used with permission
from Marshall and Plumb, 2007).
110
F IGURE 4.10 Schematic showing frontal systems at midlatitudes that result in mixing of cold and warm air via horizontal
eddies and frontal uplift (used with permission from Marshall and
Plumb, 2007).
113
S ECTION 4
F IGURE 4.11 Schematic showing key general circulation features aloft (top panel) and at the surface (bottom panel) (used
with permission from Marshall and Plumb, 2007).
In the Northern hemisphere, the high surface pressure
zones are associated with anti-cyclonic (clockwise) flow, while
the low surface pressure zones are associated with cyclonic
(counter clockwise) flow (Figure 4.12). All of these features
are relatively persistent throughout the year, but can
114
These so-called jet streams result from the geostrophic
balance (discussed in more detail in Section 5) between
latitudinal pressure/temperature gradients and the Coriolis
force and generally occur in regions of strong temperature
contrast (i.e. where circulation cells meet; Figure 4.13). These
jet streams circle the globe at fast speeds and are responsible
for much of the advection of storms. The annual average zonal
wind is shown in Figure 4.14. Note that since this is an annual
average and not a snapshot in time, the feature seen near the
top of the troposphere is really indicative of both the
subtropical and midlatitude jets which migrate latitudinally
throughout the year. In particular, the polar (or midlatitude)
line showing the key circulation features including areas of uplift, downdraft, overturning and jets (from theairlinepilots.com/
forumarchive/met/atmospherecirculation.jpg).
115
showing large jet streams (used with permission from Marshall and
Plumb, 2007).
A point worth reiterating is that the atmosphere is a
highly dynamic fluid and the above discussion has primarily
focused on its mean behavior. Other atmospheric features like
Monsoons, El Nino/La Nina, hurricanes/typhoons, etc. can be
as important as the general circulation, especially in localized
regions. This should be kept in mind. Movie 4.4 shows a
simulation of atmospheric dynamics in an attempt to more
F IGURE 4.16 Conceptual picture showing the key mechanism responsible for global heat and momentum transport
(used with permission from Marshall and Plumb, 2007).
117
with key features and the implication of these features on precipitation and evaporation distribution with latitude (from
sonoma.edu/users/f/freidel/global/207lec2images.htm).
118
M OVIE 4.4 Animation of model output showing seasonal evolution of precipitable water (white colors) and precipitation (orange) (from vets.ucar.edu/vg/CCM3T170/index.shtml).
119
S ECTION 5
Fundamental Equations of
Atmospheric Motion
The general circulation of the atmosphere described
qualitatively in the earlier sections is the product of the
governing equations of motion for the atmosphere. The full set
of governing equations are the coupled set of equations
consisting of: 1) the equation of state, 2) conservation of
momentum, 3) conservation of air mass, 4) conservation of
water mass, 5) conservation of energy. The equation of state
was already defined in Equation (2.3.4). In addition to that
form, it can be expressed in terms of specific humidity rather
than vapor pressure as:
(6.5.1)
V = [u,v,w]T
(6.5.2)
u=
dx
;
dt
v=
dy
;
dt
w=
dz
dt
(6.5.3)
= (x(t),y(t),z(t),t)
(6.5.4)
which states that the variable depends on space and time, but
that the spatial coordinates are also changing with time due
to position changes. The so-called total or material
derivative is given by:
d
dx dy dz
( ) =
+
+
+
dt
t x dt y dt z dt
=
+u
+v
+w
t
x
y
z
(6.5.5)
d
()
() =
+ (V i )()
dt
t
(6.5.6)
follow.
The conservation of momentum is simply an expression
of Newtons 2nd law, i.e.:
" du dt %
F = a = dV = $$ dv dt ''
m
dt $
'
# dw dt &
(6.5.7)
V
t visc
"
2
2
2
2
2
2
$ u x + u y + u z
$
= $ 2 v x 2 + 2 v y 2 + 2 v z 2
$ 2
2
2
2
2
2
$# w x + w y + w z
(
(
(
= 2 V;
where F is the set of forces acting on the fluid, m is the air
parcel mass, and a is the acceleration. Note that the
derivative in Equation (6.5.7) is the total derivative so that
using Equation (6.5.6):
dV V
=
+ (V i )V
dt
t
(6.5.8)
dV V
V
V
V
=
+ (V i )V =
+
+
dt
t
t press t visc t grav
(6.5.11)
kinematic viscosity
(6.5.12)
which acts only in the vertical direction. Hence, in a nonrotating reference frame the governing momentum equation
can be written concisely as:
dV V
1
=
+ (V i )V = p + 2 V g
dt
t
(6.5.13)
(6.5.9)
In a rotating reference frame an additional apparent force,
i.e. the so-called Coriolis force, becomes relevant. The rotation
is characterized by the rotation vector:
V
1
1 $ p p p '
= p = & , , )
t press
% x y z (
)
)
)
%
'
'
'
'
'&
(6.5.10)
= = 2 rad/day
(6.5.14)
V
t coriolis
$ v sin
&
= 2 V 2 & u sin
&%
0
$ fv '
'
&
)
)
=
2
& fu )
)
&
)
)(
% 0 (
(6.5.14)
Ug =
V
1
+ (V i )V = p + 2 V g 2 V
t
(6.5.15)
2 U g =
1
p
(6.5.16)
1
k h p
f
(6.5.17)
Ug =
g
k hZ g
f
(6,5.18)
An example is shown in Figure 4.19 in terms of
geopotential height (at 400 mb in the atmosphere) over the
U.S. At any given location, the direction of the geostrophic
flow can be diagnosed by the cross product. At the sample
point shown in Figure 4.19, the vertical vector is out of the
page, the geopotential gradient is southward and therefore the
122
1 p
p
g = 0
= g
z
z
(6.5.19)
The conservation of air mass equation can be written as
the balance between local mass change and convergence of
mass flux:
= i ( V) = (V i ) i V
t
(6.5.20)
123
d
= i V
dt
(6.5.21)
In the atmosphere it is often valid to assume that the lefthand-side is much smaller than the right-hand-side so that the
air mass conservation equation could be replaced by:
i V = [u x + v y + w z ] = 0
(6.5.22)
dq
= body sources/sinks
dt
(6.5.23)
" 2q 2q 2q % Sevap S
dq
= q $ 2 + 2 + 2 ' +
cond
dt
y
z &
# x
(6.5.24)
Sevap Scond
dql
=
+
dt
(6.5.25)
T d
= 2 Rn LvSevap + LvScond ;
dt
R /c
#p& d p
= T %% ((
; p0 = 100000 Pa
p
$ 0'
thermal conductivity of air
c p
(6.5.26)
125
S ECTION 6
Conceptual Questions
1. Sketch the basic picture of net radiation at the top of
atmosphere as a function of latitude. (Make sure your axes
are labeled.)
2. What part (i.e. latitudes) of the globe absorbs the most net
radiation?
3. If averaged over the whole globe, what is the net radiation
at the TOA?
4. Which is ultimately more responsible for heat
redistribution: atmospheric or oceanic motions?
126
S ECTION 7
Sample Problems
Problem 4.1. Model output from a ten-year climate model
simulation yields the average values for outgoing TOA
longwave radiative fluxes and planetary albedo shown below.
a) Plot the planetary albedo as a function of latitude. Explain
the physical reason for the latitudinal dependence of planetary
albedo. Note: Planetary albedo is the bulk (combined eects
of surface and atmosphere) reflectivity representing what
fraction of incoming TOA shortwave radiation will be reflected
back into space.
b) Plot the outgoing longwave as a function of latitude.
Explain the physical reason for the latitudinal dependence of
outgoing TOA longwave flux. Note: The outgoing TOA
longwave flux is the combined flux emitted by the atmosphere
and the flux emitted by the surface that is transmitted to the
TOA.
c) Compute the annual mean TOA incident solar radiation for
all the latitudes tabulated below and plot the annual solar
incident radiation as a function of latitude. The MOD-WET
code TOA_incoming_solar.m may be useful.
d) Compute and plot the net TOA shortwave as a function of
latitude. On the same figure plot the net TOA longwave. Also
compute, and plot on the same figure, the total net TOA
LATITUDE
PLANETARY
ALBEDO (%)
OUTGOING
LONGWAVE
(W M-2)
-90
58.61
148.1
-84
63.39
150.9
-76
64.85
160.0
-68
57.92
186.6
-60
49.30
201.5
-52
43.54
209.9
-44
38.25
220.6
-36
31.46
235.3
-28
27.62
248.3
-20
24.79
254.6
-12
23.02
254.3
-4
20.50
250.6
19.39
248.9
12
21.26
250.5
20
25.50
254.6
28
29.47
249.5
36
33.26
235.7
44
37.75
225.0
52
40.63
213.8
60
44.05
206.1
68
48.04
197.3
76
58.19
191.6
84
60.92
182.9
90
59.47
192.6
127
128
Chapter 5
Precipitation
Processes
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
2. Define the dry adiabatic lapse rate and state its value
3. Define the moist adiabatic lapse rate and describe how/
why its typical value corresponds to the dry adiabatic
lapse rate
4. List and describe the two primary growth mechanisms for
water droplets in clouds
5. Use the Marshall-Palmer (or other) size distribution
function to compute the number of droplets in a cloud and
the liquid water content in clouds
6. Draw and label a schematic of the various processes
related to cloud/precipitation physics
7. Compute the precipitation reaching the surface using the
integral relating drop size distribution, uplift and terminal
velocities, and critical drop diameter
8. List and describe the four primary lifting mechanisms
involved in the meteorology of precipitation
130
S ECTION 2
Thermodynamics of Cloud
Formation
Precipitation generally refers to the flux of water
reaching the surface and is a result of a set of processes that
ultimately convert atmospheric water vapor to hydrometeors
(either liquid droplets or ice crystals) that can fall through the
atmosphere. Clouds are just a collection of hydrometeors that
are not large enough to fall to the surface. The processes
required to form clouds (and ultimately precipitation) include
both large scale meteorological phenomena as well as
microphysical processes on the molecular/cloud scale. For a
more detailed discussion of cloud and precipitation processes
the reader is referred to Yau and Rogers (1984) and
Pruppacher and Klett (2010).
The first necessary requirement for cloud formation is a
condensation mechanism whereby water vapor molecules are
converted to liquid or ice. To understand water vapor
transformations it is convenient to use parcel theory, which
is simply a conceptual representation used to follow parcels
of air as they move through the atmosphere. To a good
approximation, moving air parcels carry vapor around with
them. In other words, the specific humidity of an air parcel (q)
is conserved (i.e. q = constant) as a parcel moves, provided
there is no condensation. Given this fact, vapor will be carried
around by a parcel right up until the parcel becomes
dp
= pg
dz
(5.2.1)
dTp
Tp
Rd dp
cp p
1
dp = 0 (adiabatic)
p
(5.2.2)
(5.2.3)
131
3) Combining Equations (5.2.1) and (5.2.2) yields the socalled dry adiabatic lapse rate, which is defined as:
dTp
g
9.8 m s 2
1
d =
=
=
9.8
K
km
dz
c p 1004 J kg -1K-1
(5.2.4)
Tp (z) = Ts d z
(5.2.5)
T(z) = Tsurf z
p(z) = psurf
!T(z) $Rd
#
&
#T &
" surf %
133
q = q s (z LCL ) =
es (Tp (z LCL ))
p(z LCL )
#L
es 0 exp % v
%R
$ v
#1
&&
1
%
((
%T T z ((
$ 0
surf
d LCL ''
g
psurf
#T z &Rd
LCL
% surf
(
%
(
T
$
'
surf
the air were drier (i.e. it would have a higher LCL). For
colder air, the parcel will also be closer to saturation and
therefore have a lower LCL. For example, if the surface
air had a temperature of 20 degrees Celsius, the LCL
would be reduced to approximately 950 m.
(9.8 10 3 K/m)
s =
(2.5 10 6 J/kg)
1+
(6.84 10 4 kg/kg/K)
(1004 J/kg/K)
9.81 m s -2
%
((
$ 461 J/kg/K $ 273.16 K 281.1 K ''
(281.1 K)2
= 6.84 10 4 kg kg -1 K-1
135
S ECTION 3
Cloud Microphysics
While cooling via adiabatic uplift is generally necessary
for cloud formation, it is not sucient. Several microphysical
processes are responsible for the formation, growth, and
ultimate conversion of vapor to precipitation.
It turns out that while condensation is generally
observed when an air parcel reaches saturation, water vapor
alone is often insucient to form hydrometeors. Rather,
particulates that are pervasive in the atmosphere (aerosols,
ash, salt, etc.), which are often referred to as cloud
condensation nuclei (CCN), are at the heart of the initial
conversion of water vapor to hydrometeors. These CCNs are
typically on the order of 0.2 micrometers in diameter
compared to a water molecule which is on the order of 0.3
nanometers (i.e. three orders of magnitude smaller). The
process of water vapor molecules spontaneously joining
together via condensation to form a liquid droplet (or ice
crystal) is called homogeneous nucleation. While this is
theoretically possible, the small size of a water molecule,
combined with the surface tension eects as a function of
radius (i.e. higher surface tension for a smaller radius), makes
such droplet formation relatively rare and/or inecient.
Instead, CCNs provide a nucleation site for the water vapor
molecules to condense onto via what is referred to as
heterogeneous nucleation (Figure 5.2). The larger nuclei makes
weathergamut.com/wp-content/uploads/2001/10/nuclei.jpg).
F IGURE 5.3 Illustration of typical relative sizes of various hydrometeors in a warm cloud (from geog.ucbs.edu/~joel/g110_w08/
lecture_notes/precip_processes/agburt07_02.jpg).
Once droplets
of sucient size exist
in the cloud, some
start to fall relative
to the updraft.
Bigger droplets will
fall faster than
smaller ones,
overtaking them
(Figure 5.4). There
will generally be
some collisions
between droplets and
some coalescence.
This collision/
coalescence growth
mechanism yields a
much faster growth
rate and is generally
responsible for the
137
unit volume of air per unit diameter bin. Figure 5.3 provides
some indication of this distribution, with small cloud droplets
being in large numbers (e.g., 106 droplets per liter of air in
this example) compared to much fewer large cloud droplets
(e.g., 103 per liter) and even fewer raindrops (e.g., 1 droplet
per liter). These values are just an example and the actual
distribution varies considerably from one cloud to another.
The most common DSD used to describe precipitating (warm)
clouds is the so-called Marshall-Palmer (1948) distribution:
N(D) = N 0 exp(cD)
(5.3.1)
c=
1
;
D
(5.3.2)
(5.3.3)
D3
LWC = w
N(D)dD
6
0
c=
(5.3.4)
1
1
=
= 2 mm -1
D 0.5 mm
Total # of drops =
N 0e cD dD =
0
N 0 cD
e
0
c
N0
N
[0 1] = 0
c
c
4
4 (100 cm)
(0.08 cm )
4
1
m
=
= 8000 drops m -3
1000 mm
2 mm 1
1m
=
139
LWC =
N0
D3
w
N 0e cD dD = w
6
6
w N 0 3! w N 0
=
6 c4
c4
De
3 cD
dD
Vt = D
(100 cm)4
(1000 kg m ) (0.08 cm )
1 m 4 1000 g
=
4
1 kg
-1 4 (1000 mm)
(2 mm )
1 m4
= 1.57 g m -3
-3
Vup = K Z ;
The DSD can also be used in the calculation of the
precipitation rate leaving the cloud base by considering the
fact that the relative fall velocity of a drop is given by the
dierence between the terminal velocity and the updraft
velocity: (Vt -Vup) and that only drops that have a positive fall
velocity (i.e. that can overcome the updraft) will actually
contribute to the flux leaving the cloud base. Given this, the
precipitation leaving the cloud base can be expressed as:
D3
Pb =
N(D)[Vt (D) Vup ]dD
6
Dup
Z = cloud thickness
(5.3.7)
(5.3.6)
(5.3.5)
To compute the precipitation that reaches the surface
one needs to account for the sub-cloud evaporation from
falling drops. This has two impacts: i) a critical diameter (Dc)
exists below which all drops will evaporate and ii) a fraction
of mass from the droplets with D > Dc will also be lost (due
to partial evaporation). To account for these, a modified
integral can be used to get the surface precipitation (P):
P=
Dmin
D3
X(D)
N(D)[Vt (D) Vup ]dD
6
(5.3.8)
(5.3.9)
E XAMPLE 5.3.2
Estimate the cloud-base precipitation rate (in
mm/hr) for the cloud described in Example 5.3.1.
First assume the prevailing updraft velocity is
negligible and for simplicity that the terminal
velocity has a linear form as in Equation (5.3.6)
with a multiplicative parameter: 1 x 103 s-1.
Qualitatively, how would the answer change if the
prevailing updraft velocity were 0.5 m/s?
Qualitatively, how would the precipitation rate at
the surface dier from that at the cloud base?
The precipitation rate for negligible updraft velocity is
given by:
Pb =
D3
cD
N 0e [D]dD = N 0 e cDD 4 dD
6
6
0
4!
1
N 0 5 = 4 N 0 5
6
c
c
4
4 (100 cm)
= 4 (0.08 cm )
(1000 s 1 )
4
1m
1
5
1 5 (1000 mm)
(2 mm )
1 m5
1000 mm 3600 s
= 3.14 10 6 m s 1
= 11.3 mm h -1
1m
1h
=
141
0.5 m s -1
Dup =
=
= 0.5 mm
1000 s -1
The cloud base precipitation would then be obtained by:
D3
Pb =
N 0e cD [D Vup ]dD
6
Dup
It should be noted that the above description strictly
applies to warm cloud processes. In cold clouds (where
temperatures reach freezing) the basic concepts are the same,
with a few modifications. First, a cold cloud generally consists
of (at least initially) a mixture of ice crystals, super-cooled
liquid droplets, and vapor. The ice crystals can take on a
variety of shapes/sizes as discussed in more detail in the next
chapter. The presence of super-cooled liquid occurs as a result
of the ineciency of homogeneous freezing of existing liquid
droplets that are lifted into the cold region of the cloud. An
additional process occurs in cold clouds as a result of
dierences in the saturated vapor pressure over bulk water vs.
bulk ice as given by their respective Clausius-Clapeyron
equations (Equations (2.4.5) and (2.4.6) respectively). Droplet
or crystal growth/decay is driven by gradients in vapor
pressure between the surface of the hydrometeor and the
ambient cloud air. For sub-freezing temperatures the
saturated vapor pressure over liquid water is higher than that
over ice. So for the same ambient vapor pressure (and same
size liquid and ice hydrometeors), the vapor pressure gradient
will often be away from the liquid droplet and toward the ice
crystal. The end-result is that liquid droplets will decay via
diusion and ice crystals will grow. Depending on the time
scale, this will ultimately lead to growth of ice crystals at the
expense of supercooled liquid droplets. For cold clouds, falling
crystals will collide and form conglomerated crystals (e.g.
snowflakes), and will melt and/or sublimate in the subsaturated cloud layers.
142
S ECTION 4
Meteorology of Precipitation
The previous two sections highlighted the need for an
uplift mechanism to initiate cloud/precipitation processes.
These mechanisms are generally the result of large-scale
meteorological processes, some of which are directly connected
to the general circulation features discussed in Chapter 4.
Here we will focus on the qualitative description of the four
main mechanisms that provide vertical uplift in the
atmosphere.
The first mechanism involves large-scale horizontal
convergence (Figure 5.6). This consists of large-scale air
masses of similar thermal characteristics coming together at
the surface and being forced upward as a result of mass
conservation. The best example of this is the ITCZ which is
the coming together of the surface branches of the Hadley
cells. As mentioned in Chapter 5, the Hadley cells are a strong
and persistent feature of the atmospheric general circulation.
This uplift is the chief mechanism of tropical clouds near the
equator, is global in extent, and occurs throughout the year
(with seasonal shifts north/south of the equator). Of
particular note is the fact that water vapor is plentiful in the
tropics due to the large amount of net radiation near the
equator that drives evaporation. This combined with the
F IGURE 5.6 Schematic showing large-scale horizontal convergence leading to uplift. This mechanism is most likely to occur
at the ITCZ and/or in tropical low-pressure systems (copyright
Pearson Prentice Hall, 2005).
F IGURE 5.10 Schematic showing typical horizontal and vertical scales of a warm front (from sci.uidaho.edu/scripter/geog100/lect/
05-atmos-water-wx/05-part-7-atmos-lifting-fronts/05-30-front-warm-diag.j
pg).
Frontal systems are often shown schematically on surface
weather maps (Figures 5.12 and 5.13), where dierent colors/
symbols are used to denote the type of front. Note that the
systems tend to form via cyclogenesis where cyclonic/anticyclonic motions centered around highs or lows initiate
movement of air masses. In addition to warm and cold fronts,
stationary fronts (where there is no relative motion of one air
mass to the other) or occluded fronts (where one front is
overtaken by another) often occur before the frontal storm
146
148
150
E XAMPLE 5.4.1
Three islands have peak altitudes of 510 m, 1020
m, and 1530 m corresponding respectively to
pressures of 950, 900, and 850 mb. Each island
experiences a prevailing wind with surface air
having a temperature of 28 degrees Celsius and
specific humidity of 12 g/kg. The air is lifted up
and over the peak of the island before returning
to sea level on the leeward side. Which, if any, of
the three islands would be expected to have
orographic clouds?
151
es (Tpeak1 )
p peak1
# 2.5 10 6 J/kg #
&&
1
1
(611 Pa)exp %
%
((
461
J/kg/K
273.16
K
296.2
K
$
''
$
= 0.622
95000 Pa
= 18.7 g/kg
q s,peak 2 =
es (Tpeak 2 )
p peak 2
# 2.5 10 6 J/kg #
&&
1
1
(611 Pa)exp %
%
((
$ 461 J/kg/K $ 273.16 K 291.2 K ''
= 0.622
90000 Pa
= 14.5 g/kg
q s,peak 3 =
es (Tpeak 3 )
p peak 3
= 0.622
# 2.5 10 6 J/kg #
&&
1
1
(611 Pa)exp %
%
((
$ 461 J/kg/K $ 273.16 K 286.2 K ''
= 11.0 g/kg
85000 Pa
152
S ECTION 5
Precipitation Climatology
and Extremes
Of particular interest to hydrologists is the expected
(climatological) precipitation, which relates directly to water
supply issues, and extreme precipitation (both positive and
negative), which relates directly to flooding or drought. Here
we focus on a brief summary of climatology and extremes.
The climatology of precipitation is simply the long-term
average. As such it tends to simply be an expression of the
typical mechanisms described above that occur over a given
region. Movie 5.1 shows an animation of global monthly
average water vapor (left panel) and precipitation (right
panel) from 2002-2012 from various data sources. First and
foremost, precipitation requires the presence of water vapor as
a source of condensate. This is clearly shown as there is a high
correlation between the two maps. However, precipitation also
requires the presence of an uplift mechanism (described in the
previous section) and other microphysical processes. Some of
these large-scale mechanisms can clearly be seen, most
notably in the form of the ITCZ. It should be noted that the
animation and other climatologies described below, especially
at global scales, rely on models and remote sensing data for
many areas since the in-situ network is sparse.
Figure 5.20 shows a 30-year climatology for annual
average precipitation based on satellite-borne measurements
F IGURE 5.21 Zonal average of satellite-derived global precipiF IGURE 5.20 Global precipitation climatology as diagnosed
from satellite data from 1979-2010 (from NASA).
large precipitation events), the physics are the same, but often
involve some rare set of circumstances leading to above
average water vapor, uplift velocities, and hydrometeor
growth. In the case of droughts, the mechanisms are less well
known. In some cases they can be simply the result of the
natural variability in the system, but are also thought to be
related to feedbacks in the climate system, where a given state
(i.e. dry) is self-sustaining. A simple example of a proposed
feedback is that if there is initially a negative soil moisture
anomaly (i.e. below average conditions), this will lead to less
available water for evaporation. The reduced evaporation will
lead to reduced water vapor, which will lead to less
precipitation. Finally, the reduced precipitation will lead to
less soil moisture. Such a positive feedback mechanism is one
hypothesis for such sustained drought events. In general,
extreme events (both flood and drought) are much less
predictable than climatology.
In describing extreme precipitation events, three
parameters are generally of interest: precipitation intensity,
storm duration, and frequency of recurrence. It should also be
noted that extreme implies some relative comparison to
climatology (which varies in space/time), i.e. what might be
an extreme event at one location may not be at another.
Analysis of extreme events is therefore generally tied to a
particular region.
To put extreme events into some context, Figure 5.23
shows an example of extreme precipitation events from across
the globe in terms of varying duration and intensity. In
particular, the data points are those events with the
155
R = 425D 0.47
(5.5.1)
F IGURE 5.23 Maximum amounts of recorded rainfall vs. duration (adapted from Dingman, 2008). Blue dots indicate actual data
points from locations across the globe and the red line is based
on Equation (5.5.1).
For design purposes, most local codes specify a particular
hypothetical design storm with a characteristic duration
and frequency (return period). For a given region, these
specifications define a particular intensity for which to use in
the design. Example maps can be found in Dingman (2005)
which provides so-called intensity-duration-frequency (IDF)
maps of intensities and durations for dierent return period
(frequency) storms. The return period (i.e. T = 100-year
storm) indicates the average recurrence interval (i.e. once
every T years). Such design storms can then be used in
rainfall-runo models (Chapter 10) for predictions of design
runo hydrographs or peak flows resulting from these design
events.
156
S ECTION 6
Precipitation Measurement
Climatology maps and IDF curves used for design
mentioned in the previous section are ultimately determined
via measurements of precipitation. The traditional method of
measuring precipitation is via a rain gauge, which is simply a
vessel open to the atmosphere that periodically or
continuously records the quantity of water it collects. Figure
5.24 shows an example of a rain gauge. Gauges come in a
variety of forms including: i) non-recording gauges (i.e.
collects cumulative amount that must be measured manually),
ii) weighing-recording gauges (vessel sits on a scale that
measures the cumulative amount of water via its weight), iii)
tipping-bucket gauges (water is routed to a small vessel of
known volume that tips when filled and records each tipping),
iv) optical gauges (measures disturbance of an infrared beam
as a result of passing drops), as well as other varieties. These
gauges are of varying complexity (and therefore cost) and
provide varying information. The first two gauge types listed
above provide cumulative amounts (one manual and one
automated), but do not provide rainfall intensities (i.e. rain
per unit time), while the tipping bucket, optical, and other
types provide the actual rainfall intensities. In using an in-situ
gauge one does not want the gauge to interfere with the
variable it is trying to measure. For this reason, many gauges
have wind shields in order to minimize the prevention of
raindrop collection due to turbulent flow around the gage
(Figure 5.24).
A key aspect to remember is that a rain gauge is a pointscale measurement. What is generally of most interest is the
spatial field, i.e. the spatial distribution of precipitation, or
the spatial average (mean areal precipitation). If a network of
gauges are deployed with sucient spatial density, then it
would be expected to provide an accurate representation of
the spatial field, or at least its mean. Sucient density in this
context is not an independent variable, but is tightly coupled
to the variability of the underlying field. For example, if a
field has no spatial variability (i.e. is constant in space) then
157
Figure 5.25 shows the global rain gauge network. The
key point is that the vast majority of gauges are in developed
countries (U.S., Europe, Japan, Australia) due to costs. So
when large-scale hydrology problems are of interest, this lack
of in-situ gauges can be a factor in characterizing
precipitation. While the U.S. has far and away the largest
network, what is not obvious from the figure is that even in
the U.S. the gage network is insucient to fully characterize
precipitation. Note for example that the number of gauges in
A classical problem in hydrology is the calculation of
mean areal precipitation (MAP), i.e. over a watershed, from a
set of discrete gauge measurements (Figure 5.26). This can be
defined mathematically as:
F IGURE 5.26 Illustrative example of Thiessen polygon construction -- basin with gauge locations.
158
= 1 p(x,y,t)dx dy
P(t)
A
A
(5.6.1)
= w p (t)
P(t)
g g
g =1
(5.6.2)
w
g=1
= 1 and wg 1
(5.6.3)
= 1
P(t)
G
p (t)
g =1
(5.6.4)
159
F IGURE 5.29 Illustrative example of Thiessen polygon construction -- perpendicular bisector construction.
160
F IGURE 5.30 Illustrative example of Thiessen polygon construction -- final Thiessen polygons.
ag
G
1
=
P(t)
A pg (t) = A ag pg (t)
g =1
g =1
G
(5.6.5)
F IGURE 5.31 Illustrative example of Thiessen polygon construction -- polygon area computation.
161
1 !
P =
(18)(15) + (63)(10) + (29)(14) + (61)(11) + (52)(10)#$
"
223
= 11.2 mm
The arithmetic average is instead:
1
P = !"(15) + (10) + (14) + (11) + (10)#$
5
= 12. mm
In this case there is about a 7% dierence between the
two estimates. One would expect the Thiessen polygon
method to be a better one as it accounts for the relative
contribution of rainfall measured at each gage. In the
case where the Thiessen polygons have the same areas,
the two estimates would be the same.
So far the discussion of the measurement of precipitation
has been limited to in-situ measurements via gauges and
estimation of MAP using the Thiessen polygon method. As
alluded to above, the high degree of spatial variability in
precipitation fields makes it dicult for a sparse gauge
network to accurately sample the fields (even using the
Thiessen polygon method). A lack of representativeness in
gauge measurements will end up propagating to the MAP
estimates via Thiessen polygon or other methods. This has led
to a significant amount of research into how to better estimate
precipitation and its spatial structure.
Chief among these newer methods are those under the
category of remote sensing. Remote sensing is becoming
pervasive in hydrology as an enabling technology beyond just
precipitation estimation. It refers to methods used to measure
163
164
Z = D 6N(D)dD
(5.6.6)
Z(x,y) = aP(x,y)b
(5.6.7)
To overcome the issues associated with ground-based
sensors (either gauge or RADAR) another method becoming
more common to hydrology is the use of satellite-based remote
sensing. Satellite-based sensors measure electromagnetic
radiation (either emitted by or reflected by the atmosphere
and surface) in various wavelengths. In principle, the same
concept used in NEXRAD can be used by mounting a
RADAR on a satellite, but to date only one such research
satellite (the Tropical Rainfall Measurement Mission
[TRMM]) exists. Instead, more common sensors are those
measuring passive microwave or infrared radiation emitted by
166
167
Z=
D N e
6
cD
dD = N 0
6!
c7
E XAMPLE 5.6.2
If a cloud is characterized by the MarshallPalmer (MP) drop size distribution, derive a
reflectivity vs. precipitation relationship in the
form of Equation (5.6.7) in terms of the
Marshall-Palmer parameters. In the derivation
use the cloud-base precipitation (which is what is
typically measured by the RADAR) and assume
updraft velocities are relatively small.
The cloud base precipitation under the assumptions
mentioned above is given by (from Example 5.3.2):
1
P = Pb = 4 N 0 5
c
Z = N0
6!
4 N 0
7 5
P 7 5
720N 0
(4 N )
7 5
P
7 5
i.e. where:
Z = N0
6!
(4 N )
7 5
P 7 5
;
7 5
b=
= aP b
with:
a=
720N 0
(4 N )
0
7
5
!
1$
c = # 4 N 0 &
P%
"
= 4 N 0
15
P 1 5
S ECTION 7
MOD-WET Codes
Relevant functions based on concepts introduced in
this chapter include:
Computation of mean areal precipitation using the Thiessen polygons:
compute_mean_precip_from_thiessen.m
Computation of saturated adiabatic lapse rate:
sat_adiabatic_lapse_rate.m
Automated Thiessen polygon construction:
thiessen.m
169
S ECTION 8
Conceptual Questions
1. If an air parcel is lifted upward by 1 km (without
condensation), by how much will its temperature change?
Will it be cooler or warmer? What is the rate at which it
changes temperature with height called?
2. Is the moist adiabatic lapse rate generally larger or smaller
than the dry adiabatic lapse rate?
3. What type of surface air parcel characteristics (in terms of
humidity and temperature) will be expected to have the
highest cloud base when lifted?
170
S ECTION 9
Sample Problems
Problem 5.1. Clouds typically contain a spectrum of drops
of varying diameters as a result of condensation and droplet
growth processes. The distribution of the number of droplets
within clouds as a function of drop size (diameter) is often
described by the Marshall-Palmer drop size distribution
(DSD).
a) Sketch the Marshall-Palmer DSD. According to this DSD,
what are the implied most probable and least probable drop
sizes? Is the location of the actual peak physically realistic?
b) Derive an expression for the total (integrated) number of
drops in a cloud as a function of the Marshall-Palmer DSD
parameter c. Based on your answer, what are the implications
on the total number of drops in a cloud as the average drop
diameter in the cloud increases?
c) Derive an expression for the total liquid water content in a
cloud (i.e. mass of liquid water per volume of air) assuming
spherical cloud droplets. This requires writing down the total
mass of drops of a given size, which is the product of the drop
size distribution and the mass of a single drop and then
integrating across all possible drop sizes. Hint: A useful
integral is:
n x
x
e dx =
0
n!
n+1
P=
Dmin
D3
X(D)N(D)
[VT (D) Vup ]dD
6
a) Write the full expression for the precipitation flux for the
simplified case of no sub-cloud evaporation, negligible updraft
velocity and a Marshall-Palmer DSD.
STATION
1.1
1.7
1.3
1.8
1.9
2.3
hand and counting squares for each polygon area. Make sure
to be very precise in the construction of the polygons.
Consider only the area inside the basin for the weights.
c) Discuss the dierence in the two annual averages you find
for the watershed. Why are they dierent?
d) Which method do you think provides the most
representative mean areal precipitation estimate over the
basin? Why? So which method should generally be employed
for estimating mean areal precipitation?
Problem 5.4. Ground-based RADAR is now becoming more
commonly used in the U.S. to estimate precipitation. The
basis for this remote sensing technology is that microwave
pulses can be sent out into the atmosphere and are then
backscattered based on how many cloud/rain drops there are
in the field of view. From this backscattered power the rain
D3
P = N(D)
[VT (D)]dD
6
0
which corresponds to the case where updraft velocity is small
relative to terminal velocity and sub-cloud evaporation is
negligible.
Using the Marshall-Palmer distribution and assuming the
terminal velocity is directly proportional to diameter, i.e.:
VT (D) = D
along with the simplifying assumption that terminal velocity
is much larger than updraft velocity (VT>> Vup) integrate the
above precipitation equation to express precipitation (P) in
terms of all other parameters. Rewrite the equation to express
c as a function of all other variables.
173
Z = aP b
d) What is the theoretical value of the coecient b? What
terms make up the coecient a? For a cloud droplet terminal
velocity parameter of 3000 s-1, determine the constant a in
units of mm8/5 hr7/5. You should use the standard value for
the Marshall-Palmer distribution parameter N0.
e) Demonstrate that a=325 and b=7/5, if the units for Z are
mm6 m-3 and the input units of P are mm hr-1. In other words,
show that the following is true:
VT (D) = D;
= 3000 s 1
Chapter 6
Snow
Processes
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
1. Define the surface energy balance and name and define the
various fluxes involved
2. Define and be able to compute snow porosity, liquid water
content, snow density, and snow water equivalent
S ECTION 2
Rn = LE + H + G
(6.2.1)
It is often said that the net radiation gets partitioned
into the other three fluxes on the right-hand-side. How it gets
partitioned is very much a function of the surface conditions
178
LE = LvE or LE = LsE
(6.2.2)
180
S ECTION 3
Snowpack Characteristics
In areas with sub-freezing temperatures, the processes
described in Chapter 5 can lead to solid-phase precipitation.
Such snowfall events generally require cold cloud processes
aloft along with sub-freezing air that extends all the way
down to the surface. Exceptions to this include hail storms
where there are cold cloud processes aloft (often due to
towering cumulus clouds), but above-freezing air in the subcloud layer. In such cases, the hail stones are generally large
enough (and fall fast enough) that the above-freezing subcloud temperatures are not sucient to fully melt them before
they reach the surface. In the context of this chapter we will
mostly be discussing the former case of snowfall leading to a
variety of snowpack processes. For a more detailed discussion
of snow processes, the reader is referred to Male and Gray
(1981) and Armstrong and Brun (2010).
Snowfall generally consists of snow grains (or
snowflakes), rather than droplets, that are often made up of a
large collection of tiny ice crystals (Figure 6.2). A snowpack is
simply an accumulation of snow grains on the ground. These
grains come in various sizes and morphologies depending on
the atmospheric conditions in which they were generated.
Figure 6.3 shows some typical types (habits) of snow grain
crystals as a function of atmospheric temperature. Unlike rain
droplets which are essentially the density of liquid water, the
(6.3.1)
182
s =
Va +Vw
Vs
(6.3.2)
Vw
Vs
(6.3.3)
where the liquid water content is zero for a dry snowpack and
positive for a wet snowpack.
The snow density is given by the ratio of mass of solid
ice and liquid water (where air mass is neglected) to the total
volume of snow, i.e.:
s =
M i + M w iVi + wVw
=
= i (1 s ) + w
Vs
Vs
ret
s2
s
4
= 0.0735 + 2.67 10 ;
w
w
(6.3.5)
[ ] = kg m -3
Finally, and perhaps most importantly, we can define the
variable that is typically of most interest to hydrologists,
namely the snow water equivalent (SWE). This is simply
defined as the equivalent amount of liquid water mass that
would result if the entire snowpack were melted and is given
by:
SWE =
( )Vs + (1 s )Vs i w
A
(6.3.6)
(6.3.4)
where the numerator represents the liquid water in the
snowpack (first term) plus the equivalent liquid water bound
183
SWE = + (1 s ) i w ds = ds
w
(6.3.7)
SWE = ( s
n
w )dz n
(6.3.8)
F IGURE 6.5 Examples of variation in snow profiles as a function of location and time of season. Each panel corresponds to
a different site. Bars correspond to density on the top axis and
dotted lines correspond to temperature on the bottom axis.
Black and grey correspond to March and April conditions respectively (from Filippa et al, 2010).
! 200 kg m -3 $
SWE = #
&(2 m) = 0.4 m = 40 cm
-3
" 1000 kg m %
The liquid water content is known to be equal to 0 since
the snowpack is sub-freezing. Liquid water only exists in
a snowpack when it is at freezing. The snow porosity is
given by (rearranging Equation (6.3.4)):
s
200 kg m -3
s = 1 = 1
= 0.78
i
917 kg m -3
Hence the snowpack at the measurement time is 78%
pore space and 22% solid ice.
185
S ECTION 4
d(SWE)
= P E M
dt
(6.4.1)
d(SWE)
= P SWE =
dt
P dt
2
ei (r) = ei ()exp
rRv wT
(6.4.2)
Density however is not as straightforward as SWE. New
fallen snow generally has a relatively low density (100-200 kg
m-3). As soon as snow begins to accumulate, gravity and other
metamorphism factors begin to take hold with the end result
being a densification (i.e. increase in density) over time
(Figure 6.7). The primary mechanisms at work include:
1. Gravitational settling: This mechanism causes a
densification of the snowpack as a result of the weight of
overlying snow, which compacts the snowpack, reducing the
porosity.
F IGURE 6.7 Typical snow depth (upper panel) and snow density (lower panel) time series during accumulation season for a
high-elevation location with seasonal snowpack (in this case in
the Sierra Nevada).
187
L
e(z) = ei (T(z)) = es 0 exp s
Rv
1
1
T0 T(z)
(6.4.3)
188
E XAMPLE 6.4.2
Suppose that during the accumulation season a 2
meter snowpack has a linear temperature profile
of: T(z) = mz; z 0
where m is the temperature lapse rate in the
snowpack, with a value of -2 K/m. Derive an
expression for the vapor pressure profile in the
snowpack and describe the implied vapor flux.
The vapor pressure is, to a good approximation, equal to
the saturated vapor pressure (with respect to ice) since
the air is in equilibrium with the solid ice matrix. The
vapor pressure gradient can thus be written using the
chain rule in terms of the temperature gradient:
de dei dei dT
=
=
dz dz dT dz
The first term on the right-hand-side is simply equal to
the Clausius-Clapeyron equation (Equation (2.4.4), but
for bulk ice) and the second term is equal to the
temperature gradient from the temperature profile, so
that the vapor pressure gradient is equal to:
#L
Ls
de Ls ei T
=
=
e exp %% s
2
2 s0
dz Rv T z RvT(z)
$ Rv
#1
1 && T
%%
((((
$T0 T(z) '' z
190
"L
mLs
de Ls ei
=
m=
e exp $$ s
2
2 s0
dz Rv T
Rv (mz)
# Rv
"1
1 %%
$$
''''
#T0 mz &&
Another key snow property related to metamorphism is
snow albedo. In general, snow has a very high reflectivity (in
the visible part of the spectrum) compared to almost all other
media, with albedos as high as 90% in some cases. Since
downwelling shortwave is the largest energy flux in many
cases, the albedo can strongly regulate the overall heat
absorption by the snowpack. For snow however, the albedo is
not a constant, but evolves, primarily due to grain size
evolution in the snowpack and addition of lower albedo
contaminants (i.e. dirt). Figure 6.11 shows an illustrative
191
where ki is the thermal conductivity (in this case for ice). The
energy balance can then be written as:
(iciT )
= (Heat flux) + q
t
(6.4.4)
As described above, the metamorphism in the snowpack
is to a large extent driven by temperature profiles (and
specifically gradients). Since most of the transport of energy
(and vapor) occur in the vertical, the energy budget for the
snowpack can be written schematically (i.e. for a given layer)
as:
(iciT ) 2(kiT )
=
+q
2
t
z
Heat flux = ki (T z)
(6.4.5)
G = Rn LE H = Rs(1 ) + Rl Rl LE H
(6.4.6)
From Equation (6.4.6), it becomes clear that the primary
heating terms of the snowpack are the net shortwave radiation
(where albedo plays a significant role), the downwelling
longwave radiation, and in cases where the air temperature is
warmer than the snowpack surface, the sensible heat flux. The
primary cooling mechanisms are the outgoing longwave
radiation, the latent heat flux, and in cases where the
snowpack temperature is warmer than the overlying air, the
sensible heat flux. These fluxes drive the energy balance of the
snowpack, which in turn drives the various metamorphism
processes described above (Figure 6.12). During the
accumulation season the net radiation is often negative (i.e.
193
S ECTION 5
Snowmelt
In an idealized snowpack, the snowmelt season is that
period generally following the accumulation season. In truth,
the two seasons are generally not completely distinct and
depend on the climatology of the region. At certain elevations,
snowpack is quite transient, where many accumulation and
melt events can occur throughout the winter and spring. In
any case, the melt period generally begins when the surface
net radiation transitions from negative to positive (Rn > 0).
Melt itself will not occur immediately when this transition in
net radiation occurs due to the thermal inertia of the
snowpack. At the time net radiation becomes positive, the
snowpack will generally be at sub-freezing temperatures. The
positive net radiation implies an energy input to the
snowpack. The specific impact this energy input has defines
three commonly described phases of the snowmelt season.
The first phase is the so-called warming phase. This
phase consists of the period of time when the positive net
radiation goes into increasing the temperature of the
snowpack toward the melting point of water. At any given
time, the amount of energy required to raise the average snow
temperature to the melting point is generally referred to as
the cold content of the snowpack, which can be defined as:
Qcc = ci wSWE(Tsnow Tm )
(6.5.1)
194
The second and third phases of melt have to do with the
actual conversion of solid ice to liquid water. Keep in mind
that phase transition of a bulk media occurs at an isothermal
temperature (0C for water) so that during the last two
phases the temperature of the snowpack is constant
(isothermal) and at the melting point. By definition, the
snowpack is wet during these phases. The second phase is
referred to as the ripening phase and begins when the
snowpack becomes isothermal. Any cumulative energy input
beyond Qm1 will start to melt the snow. However early on in
the process, water is retained in the pores (i.e. none leaves the
snowpack as runo). The amount of water retained is that
given by Equation (6.3.5). The net energy input required to
complete the ripening phase (Phase 2) is simply:
Qm 2 = retds wLf
(6.5.2)
(6.5.3)
Qm 2 + Qm 3 = SWE wLf
(6.5.4)
It is important to keep in mind that the above described
accumulation and melt of snow is highly seasonal over much of
the globe. Movie 6.1 shows a monthly time series of images of
snow cover over the northern hemisphere over North America
(from January 2004-December 2004). Accumulation and melt
E XAMPLE 6.5.1
= 13.4 10 6 J/m 2
The required energy input for the melting phase is given
196
197
S ECTION 6
Impact of Vegetation
The impacts of vegetation on snow processes (and land
surface processes in general) are, among others, those related
to energy balance at the surface and those related to mass
balance. Each is discussed briefly below.
In terms of radiative energy balance at the snow surface,
the primary impact of vegetation is the attenuation and/or
augmentation of radiative energy fluxes. The clear- and
cloudy-sky shortwave fluxes discussed in Chapter 3 dealt with
those fluxes reaching the surface. In the case of bare soil/snow
these fluxes would directly contribute to the SEB. When
vegetation is present, those fluxes would correspond to those
reaching the surface just above the canopy. The downwelling
above-canopy shortwave radiation will interact with the
vegetation canopy, with the vegetation having media
properties as defined in Chapter 3 including transmissivity,
absorptivity, and reflectivity. The vegetation reflectivity
(albedo) is generally much lower than for snow. In the case of
a dense vegetation canopy (i.e. forest), the transmissivity will
also be relatively low. Together this means that much of the
above-canopy radiation will be attenuated (absorbed) by the
vegetation.
The interaction of solar radiation with canopy is
relatively complicated. To avoid complex radiative transfer
calculations, simple empirical expressions have been developed
(6.6.1)
fsv = exp(3.91F)
(6.6.2)
The impact of vegetation on mass balance of the
snowpack on the ground (or the surface in general) primarily
has to do with the so-called interception of precipitation by
the canopy. The leaves of a dense canopy can hold a nonnegligible amount of snow (or rain) that is at least initially
removed from the snow/soil surface mass balance. An example
of this in the case of snow in an evergreen forest is shown in
Figure 6.13. The interception is most often treated as a small
reservoir that captures any precipitation coming in contact
with the canopy and can hold up to a specified volume (or
(6.6.3)
d = (1 F)flc (cloud)a + F
(6.6.4)
Rl = dTa4
(6.6.5)
199
depth, i.e. volume per unit area) of water before the water
begins to drain or slough o (in the case of snow) the canopy
onto the surface. The amount of precipitation that falls
through the spaces in the canopy combined with the drainage
from the canopy is usually referred to as throughfall. The
primary impact of the intercepted water is that it can then
evaporate/sublimate back into the atmosphere before ever
reaching the surface. Hence, some of the water that would
have contributed to the snowpack and/or soil mass balance in
a bare surface case, will be removed from the surface mass
balance when vegetation is present. The evaporation (or
interception) loss from the canopy is discussed in more detail
in Chapter 8.
Both of the energy and mass balance implications
associated with vegetation cover discussed above play a
significant role in the spatial and temporal variability seen in
snowpacks. Moreover, they extend to non-snow covered
regimes, where the impact on energy/mass budgets are
equally important. This includes a significant impact on
evapotranspiration at the surface, which is discussed in more
detail in Chapter 8.
200
S ECTION 7
Snow Climatology
The climatology of snow is largely just a climatology of
precipitation, which drives snow accumulation, and a
climatology of radiation, which drives snow melt. Together
these will determine various climatological factors related to
snow including: average SWE, duration of accumulation and
melt seasons, types of snow, etc. Snowfall occurrence is largely
a product of latitudinal and elevational dependence. Areas at
high latitudes and/or at high elevation will experience air
temperatures cold enough to turn any appreciable
precipitation into snowfall.
Figure 6.14 shows the major mountain ranges across the
globe which experience significant seasonal snowfall and
accumulation as a result of elevational lapse rate eects
(despite being at midlatitudes). Much of the remaining
portions of snow covered areas of the globe are a result of less
radiation at higher latitudes. The combination of varying
elevation, latitude, and other factors leads to varying
climatological classification of snow types. A common
classification involves six dierent types of snow (Sturm et al.,
1995):
Tundra: Thin, cold, wind-blown snow usually found above or
poleward of tree line
F IGURE 6.16 Map of mean annual total snowfall over the continental U.S. (in inches).
202
F IGURE 6.17 Contour map showing the average annual number of days with snow cover.
S ECTION 8
Snow Measurement
In describing snow measurement one must distinguish
between snowfall measurement and snowpack measurement.
The former attempts to measure precipitation in the form of
snow. Hence, the same techniques discussed in Chapter 5 can
theoretically be used to measure snowfall, including snowfall
gauges, RADAR, and satellite-based techniques. The primary
complication with snowfall gauges have to do with the fact
that wind eects can be more amplified for low density
snowflakes (compared to raindrops) in the sense that the
turbulence caused by the presence of the gauge itself can
much more easily deflect snow such that it is not captured by
the gauge. Hence gauge under-catch, which is prevalent in the
case of rain, is generally exaggerated in the case of snowfall
measurement. It is not uncommon for a gauge to catch only
50% of the snowfall that might actually be falling, which can
lead to large (negative) biases in snowfall estimates.
Additionally, snow can freeze or block the gauge opening, so
often gauges designed to measure snowfall are equipped with a
heating unit to melt the snow which can then flow into the
tipping bucket or water reservoir.
The RADAR-based techniques used for snowfall
estimation also have diculty, mainly caused by the much
more complicated relationship between snow hydrometeors
and reflectivity. While liquid raindrops are usually assumed to
F IGURE 6.18 Example of snow depth measurement (from National Weather Service, Wilmington Ohio).
204
F IGURE 6.21 Photo of a snow pillow site prior during the dry
season.
Another in situ point-scale measurement that is
impractical for operational implementation, but often used for
more detailed measurement of snowpack stratigraphy is the
digging of snow pits. This technique involves digging a pit to
ground level to expose a face of the snowpack (Figure 6.23).
Once exposed, various profile measurements including
temperature, density, SWE, hardness, grain size, etc. at a
206
snow pit where measurements of temperature, density, hardness, grain size profiles are made.
207
regular interval across the entire profile are taken from the
face of the snow pit.
While the remote sensing of snowfall is dicult as
described above, other aspects of snow are routinely measured
via satellite-based remote sensing. Satellite-based remote
sensing generally measures electromagnetic radiation in a
particular and relevant part of the spectrum. Snow generally
has two relevant wavelength bands that provide useful
information: visible/near-infrared (Vis/NIR) and microwave.
In the Vis/NIR (i.e. shortwave) part of the spectrum
snow is highly reflective compared to other media (Table 3.2).
Satellite sensors measuring in this part of the spectrum are
208
209
S ECTION 9
MOD-WET Codes
Relevant functions based on concepts introduced in
this chapter include:
Disaggregation of meteorological station forcings:
distribute_met_forcings.m
Computation of surface pressure over complex terrain:
disaggregate_press.m
Computation of near-surface specific humidity over complex terrain:
disaggregate_qair.m
Computation of air temperature over complex terrain:
disaggregate_Tair.m
A simple 1-layer snow model with mass/energy balance:
snow_model.m
A snow albedo model based on the USACE formulation:
albedo_usace.m
210
S ECTION 10
Conceptual Questions
1. Name the four fluxes in the surface energy balance
equation. What are the units on the fluxes?
2. Describe what the porosity of the snowpack represent.
3. What is the liquid water content of a dry snowpack. What
(order of magnitude) is the maximum liquid water content
of snowpack?
4. Describe the meaning of the liquid water retention capacity
of a snowpack.
S ECTION 11
Sample Problems
Problem 6.1. Measurements of snowpack in three locations
in a small basin yield the following measurements:
LOCATION
DEPTH (CM)
95
102
105
SWE (CM)
27
27
29
TEMPERATURE
(C)
-3
-5
TEMPERATURE
(C)
-1
50
-3
100
-2
Chapter 7
Unsaturated
Flow and
Infiltration
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
1. Define the meaning of a porous media
2. Use the soil triangle to define a particular soil type/
texture
3. Define the volumetric water content and relative soil
saturation and how they are related to each other
4. Describe what factors (i.e. pressure head, elevation head,
velocity head) control flow in porous media
5. Define matric head and hydraulic conductivity and how
they depend on soil moisture and soil type
6. Define and explain Darcys Law for flow in unsaturated
conditions
7. Understand and describe the two terms that in particular
control vertical flow in unsaturated soils
215
S ECTION 2
or relative saturation:
In terms of characterizing porous soil media, many of the
same or analogous definitions that were used for snow are
used for soil. The mineral particle density is given by:
m =
Mass of minerals
Volume of minerals
(7.2.1)
b =
Mass of minerals
Volume of soil
(7.2.2)
and varies depending on the soil type and porosity. The soil
porosity has the analogous definition used in snow (in this
case volume of pore space per unit volume of soil) and can be
related to the above densities:
s = 1 b m
(7.2.3)
Volume of water
Volume of soil
(7.2.4)
s=
(7.2.5)
Most of the above defined properties depend on soil
type. However this has not yet been rigorously defined. The
common way of defining soil type (or texture) is via a grain
size distribution. This involves taking a soil sample and
sorting it with progressively smaller sieve openings to
characterize the distribution of soil grains in various size bins.
For the types of soil commonly found in the vadose zone,
three end-member types are defined: sand (size range between
0.06 and 2 mm), silt (size range between 0.002 and 0.06 mm)
and clay (sizes less than 0.002 mm). For a given soil sample
one can determine the fractional weight in each of these end217
ture triangle for determining soil type from fraction of clay, silt,
and sand.
POROSITY
Sand
0.395 (0.056)
Loamy sand
0.410 (0.068)
Sandy Loam
0.435 (0.086)
Silt Loam
0.485 (0.059)
Loam
0.451 (0.078)
0.420 (0.059)
0.477 (0.057)
Clay loam
0.476 (0.053)
Sandy clay
0.426 (0.057)
Silty clay
0.492 (0.064)
Clay
0.482 (0.050)
218
pores are still filled is called the field capacity. The typical
situation is shown where a mix of air and water is stored in
the pores, where water may be held up against gravity due to
capillary forces (discussed more below). Finally, if additional
drying is experienced (i.e. by evaporation or root uptake), the
soil moisture condition will eventually reach the wilting point.
In this case, only immobile water that is essentially adhered to
grains is left behind. In this condition, a plants system will
not be able to uptake moisture and may wilt. The amount of
water corresponding to each of these states is dependent upon
the soil type as shown in Figure 7.5. The reason for this has
219
E XAMPLE 7.2.1
A soil sample taken from the field with a 10 cm
long by 5 cm diameter cylindrical tube has a field
weight of 306 g and an oven-dried weight of 245
g. For the soil sample, calculate the: a) bulk
density and porosity and b) volumetric soil
moisture and relative saturation.
a) The volume of the (cylindrical) soil sample is given
by:
Vsoil =
b =
245 g
= 1.25 g cm 3
3
196.4 cm
220
61 cm 3
=
= 0.31
196.4 cm 3
The relative saturation is then given by:
s=
0.31
= 0.59
0.53
221
S ECTION 3
1 WV 2 pW
Total Energy = Wz +
+
2 g
w g
(7.3.1)
V2
p
H =z+
+
2g w g
(7.3.2)
H z+
p
= h(x,y,z,t)
w g
(7.3.3)
F IGURE 7.6 Illustration of capillary rise as a function of capillary tube diameter. Capillary rise is proportional to suction in
water. Soil suction in soil pores act analogously with smaller
pores holding water at higher suction pressures. The capillary
rise is marked by a curved meniscus due to surface tension.
h = +z
(7.3.5)
The flow (flux of water) in porous media is often
described by the volumetric flow (Q; [Q]=L3/T) per unit
cross-sectional area (A; [A]=L2) and is given by Darcys
223
(1856) Law:
q=
Q
= Kh
A
(7.3.6)
() () ()
+
+
x y z
h
=
+1
z
z
q z = K z
h
z
h(z,t) = (z,t) + z
(7.3.8)
() =
(7.3.7)
(7.3.9)
Based on Equations (7.3.7) and (7.3.9), flow in the
unsaturated zone depends on the hydraulic conductivity of the
soil and matric head of the soil. In particular these two
variables depend on two factors: i) Soil type (which is fixed in
time, but varies in space) and ii) soil moisture (which varies in
both space and time). To quantify these relationships,
experiments have been performed to develop constitutive
relationships between the soil type, moisture and conductivity
and matric head.
Several models exist, but a commonly used one is the
Brooks-Corey model (1966). The water retention curve (i.e.
soil matric head vs. soil moisture) is given by:
( ) = s
s
(7.3.10)
b Brooks-Corey parameter
K( ) = K s
s
(7.3.11)
c = 2b + 3
HYDRAULIC
CONDUCTIVITY
[cm/hr]
SAT.
MATRIC
HEAD
[cm]
BROOKSCOREY
PARAMETER
Sand
63.36
-12.1 (14.3)
4.05 (1.78)
Loamy sand
56.16
-9.0 (12.4)
4.38 (1.47)
Sandy Loam
12.49
-21.8 (31.0)
4.90 (1.75)
Silt Loam
2.59
-78.6 (51.2)
5.30 (1.96)
Loam
2.50
-47.8 (51.2)
5.39 (1.87)
Sandy clay
loam
2.27
-29.9 (37.8)
7.12 (2.43)
0.61
-35.6 (37.8)
7.75 (2.77)
Clay loam
0.88
-63.0 (51.0)
8.52 (3.44)
Sandy clay
0.78
-15.3 (17.3)
10.4 (1.64)
Silty clay
0.37
-49.0 (62.1)
10.4 (4.45)
Clay
0.46
-40.5 (39.7)
11.4 (3.7)
Silty clay
loam
Given these relationships, they can be used in Darcys
Law to express flux in terms of the soil moisture, i.e.:
225
(b+1)
" %
d
b
= s $$ ''
d
s # s &
(7.3.13)
where while the magnitude can vary with soil moisture, the
sign of Equation (7.3.13) is always positive (Figure 7.7).
F IGURE 7.7 Plot of matric head (blue curve/left axis) and hydraulic conductivity (green curve/right axis) as a function of
volumetric soil moisture content.
d
q z = K( )
+ 1 = K( )
K( )
z
d
(7.3.12)
The key point illustrated by Equation (7.3.12) is that the
flux at any depth in the soil is driven by two factors: i)
capillary forces driven by soil moisture gradients (first term)
and ii) gravity drainage (second term). The gravity term is
always acting downward (negative). The matric head term,
however, can act either upward or downward, depending on
the sign of the soil moisture gradient. If the soil moisture
gradient is positive (i.e. increasing soil moisture as you go up
in the soil column), then the matric head (gradient) term is
negative and reinforces the gravity term. If the soil moisture
gradient is negative (i.e. decreasing soil moisture as z increases
upward in the soil profile), then the matric head (gradient)
term is negative and acts against gravity. In both cases, the
matric head gradient term acts to drive flow from wetter soil
locations to drier soil locations. Depending on the conditions,
if the gradient is negative and large in magnitude, the first
term may be larger than the second so that flow will actually
be upward in the soil. The relative magnitude of the terms
and the sign of the matric head term will vary (significantly)
in time. During infiltration events, the soil moisture in the
surface layer will become moister than the soil below, driving
flow downward in conjunction with gravity. During a drydown period, where evaporation removes moisture from the
soil, the gradient in soil moisture may reverse and become
226
340 cm
fc = s
1
b
(7.3.14)
wp
15000 cm
= s
1
b
(7.3.15)
E XAMPLE 7.3.1
Soil moisture sensors in the field show that the
volumetric soil moisture profile follows the
following linear relationship: (z) = surf mz
where z is the depth below the surface (i.e. is
zero at the surface and negative below the
surface) and m is a positive constant. Determine
an expression for the Darcy flux in the soil for
this soil moisture profile. Estimate the magnitude
and direction of the flux at a depth of 0.5 m if
the soil is a sandy loam, the slope m is equal to
0.05 m-1 and the surface volumetric soil moisture
is equal to 0.4.
If the soil moisture varies with depth, then so will the
matric head and hydraulic conductivity. The matric
(pressure) head in the soil is given by:
b
b
! mz $
! (z) $
&
(z) = s ##
&& = s ## surf
&
s
" s %
"
%
c
! mz $
! (z) $
&
K(z) = K s ##
&& = K s ## surf
&
s
" s %
"
%
# &
#d &
q z = K( )% + 1( = K( )%
+ 1(
$ z
'
$ d z
'
(b+1)
= 1.87 cm/hr
bm # surf mz &
(
=
%
(
z
s s %$
s
'
where
" (z) %
d
b
= s $$
'
d
s # s '&
(2(4.9)+3)
(b+1)
" mz %
b
'
= s $$ surf
'
s #
s
&
and
(4.9)(.0005 cm 1 )
=
(21.8 cm)
(0.435)
(4.9+1)
= m
z
" mz %
'
q z = K s $$ surf
'
#
&
s
(b+1)
(
+
"
%
mz
*bm $ surf
'
+
1
s
$
'
*
s
#
&
s
)
,
" mz %
'
K(z = 0.5 m) = K s $$ surf
'
#
&
s
228
= q
t
S ECTION 4
Modeling Unsaturated
Zone Flow Dynamics
= (q z ) = K( )
+ 1
t
z
z
z
The modeling of unsaturated flow (or the unsaturated
flow equation) does not consist solely of Darcys Law, but
requires consideration of water going into or out of storage,
which changes the flow properties (matric head, hydraulic
conductivity). This can be done via application of the
dierential form of mass balance:
M
= ( wq dx dy dz)
t
(7.4.1)
M = w dx dy dz
(7.4.3)
(7.4.2)
(7.4.4)
=
K(
)
+
1
t z
d
(7.4.5)
D( ) = K( )
d
d
(7.4.6)
In the case where the soil has vegetation roots within the
unsaturated zone, then a sink term S can be added:
"
(7.4.7)
The root uptake sink term (S) can vary with depth and time
depending on the root density structure and available soil
moisture. In the case of no vegetation or root uptake, S = 0
and Equation 7.4.7 reduces to 7.4.5. This becomes the
229
2
=D 2
t
z
(7.4.8)
" z
%
'
(z,t) = 0 + (i 0 )erf $
1/2
$# 2(Dt) '&
(7.4.10)
erf(y) =
e x dx
(7.4.11)
#D &
q z (z,t) = D
= (i 0 ) % (
z
$ t '
# z
&
(
exp %
% 2(Dt)1/2 (
$
'
(7.4.12)
" ,
$
=# i
$% 0 ,
z 0,t = 0
(7.4.9)
z = 0,t > 0
230
E XAMPLE 7.4.1
E XAMPLE 7.4.1
" mz %
'
q z = K s $$ surf
'
#
&
s
(b+1)
(
+
"
%
mz
bm
surf
* $
'
+ 1'
* s$
s
#
&
s
)
,
(b+2)+c
# mz &
bcm
(
= K s 2 s %% surf
(
t
s
$
'
s
2
K sb(b + 1)m 2
s2
c1
cm # surf mz &
%
(
+ Ks
(
s %$
s
'
(b+2)+c
# mz &
(
s %% surf
(
$
'
s
(b+2)+c
# mz &
bm
(
= (c + b + 1)K s 2 s %% surf
(
s
$
'
s
2
c1
cm # surf mz &
%
(
Ks
(
s %$
s
'
#
&,
)
= (q z ) = +K(z)% (z) + 1(.
t
z
z *
$ z
'#
&
2
=
K(z) % (z) + 1( + K(z) 2 (z)
z
z
$ z
'
c1 #
(b+1)
&
)
,
)
,
mz
mz
cm
bm
. % s + surf
.
=
K s ++ surf
+ 1( +
.
+
.
%
(
s
s
s
*
- $ s
*
'
) mz ,
.
K s ++ surf
.
*
s
(b+2) &
#
2
)
,
mz
%b(b + 1)m + surf
(
.
s+
2
.
%
(
s
s
*
$
'
231
d
f (t) = K( [z = 0,t ])
d z
S ECTION 5
Infiltration
A complete modeling of unsaturated zone dynamics
requires the solution of Richards equation. If the soil moisture
profile evolution is known, then everything about storage and
flux in the entire unsaturated zone is known (or can be
determined). The one aspect of unsaturated flow that is of
most interest to hydrologists is infiltration, which is the flux
at the surface resulting from a precipitation event. The prime
motivation for modeling infiltration is that, if infiltration and
precipitation are known, then the surface runo is known. The
infiltration rate can be defined as:
z =0
(7.5.1)
+ K( [z = 0,t ])
(7.5.3)
z =0,t
F(t) =
f (t )dt f (t) =
0
dF
dt
(7.5.4)
f (t) = K
+ 1
z
z =0
which can be written more explicitly as:
(7.5.2)
f = P;
t tp
At t = t p :
(z = 0,t p ) = s ;
K( [z = 0,t p ]) = K s ;
0<
=
z z
z =0,t =t p
(z = 0,t > t p ) = s ;
K( [z = 0,t > t p ]) = K s ;
0<
=
z z
<
t >t p
t =t p
Recall from Equation (7.5.3) that the first term in the
infiltration is proportional to the soil moisture gradient (at the
surface), while the second term is equal to the hydraulic
conductivity (at the surface). Since the gravity term at the
surface is constant after ponding, the implication is that as
the soil moisture gradient decreases over time the infiltration
rate will also decrease over time. Figure 7.9 shows a schematic
of the actual infiltration rate as a function of time, which can
be written mathematically as:
P, t t t
0
p
f (t) = *
f (t), t p t tr
(7.5.5)
f (t) = f *(t),
t 0 = t p t tr
(7.5.6)
f (t) = P,
t 0 t tr
(7.5.7)
235
S ECTION 6
Infiltration Capacity
Models
In the discussion of infiltration in the previous section it
was mentioned that after ponding the infiltration rate should
be described by a monotonically decaying function f*(t). To
gain some insight into that function, and ultimately derive it,
we can start by examining models for infiltration
capacity (or potential infiltration): fc(t), which is the
infiltration rate that would occur under immediately ponded
conditions. Analytical models for the infiltration capacity can
be developed under the following simplifying conditions:
1. Ponded upper boundary condition, i.e.:
(z = 0,t) = s
2. The unsaturated zone is infinitely deep, i.e. the
groundwater table is not near the surface
3. The initial condition is uniform with depth, i.e.:
(z,t 0 ) = 0
Under these conditions, two commonly used models for
infiltration capacity can be developed: the Philip solution and
the Green-Ampt model. Note that the analytical model for
flow shown in Equation (7.4.12) was derived under similar
fc (t) = q z
z =0
(7.6.1)
The Philip (1960) solution keeps the first two terms, which is
expected to be valid under short times:
fc (t) =
Sp
2
(7.6.2)
t 1/2 + K s
2b + 3
S p = (s 0 )K s s
b
+
3
1/2
(7.6.3)
This solution predicts that the decay occurs to the order t -1/2,
and asymptotes to Ks as expected. Interestingly, the constant
diusivity case shown in Equation (7.4.12) (when evaluated at
the surface: z = 0) shows the same time dependence. In some
versions of the model, Ks is replaced by a fitting parameter.
Also note the sorptivity is only a function of soil properties
and the (uniform) initial condition (before infiltration begins).
Based on the above model, the cumulative infiltration
capacity can be derived directly:
Fc (t) =
f (t )dt = S t
t0
1/2
+ K st
(7.6.4)
236
The second model often used for infiltration capacity is
the so-called Green-Ampt (1911) model. This model uses the
same assumptions as those above, with one additional
simplification: the wetting front is a discontinuous front at
depth z = -Lf (Figure 7.10). The Green-Ampt model applies
Darcys Law across the wetting front:
qz = K s
h hsurf
dh
K s f
dz
z f z surf
(7.6.5)
h f = f + z f = f Lf
(7.6.6)
2b + 3
f
b + 3 s
(7.6.7)
(7.6.8)
fc (t) = q z = K s
Lf + f
Lf
= K s 1 +
Lf
(7.6.9)
This is not the final solution however since the length of the
wetting front (Lf) itself depends on the amount of cumulative
infiltration, i.e.:
Fc (t) = Lf (s 0 )
F IGURE 7.10 Conceptualization used to represent the wetting front in the Green-Ampt model for infiltration capacity
(adapted from Mays, 2005).
(7.6.10)
Lf =
Fc (t)
(s 0 )
(7.6.11)
237
f (s 0 ) dF
= c
fc (t) = K s 1 +
Fc
dt
(7.6.12)
f (s 0 )
1
t=
Fc (t) + f (s 0 )ln
Ks
Fc (t) + f (s 0 )
(7.6.13)
238
S ECTION 7
P, t 0 t t p
f (t) =
fc (t tc ), t p t tr
(7.7.1)
using the time-compression approximation. The dotted line represents an infiltration capacity model (i.e. under immediately
ponded conditions) which is shifted by the compression time.
239
fc (t tc )
t =t p
= fc (t p tc ) = P(t p )
(7.7.2)
tp
t p tc
P(t)dt =
fc (t)dt
(7.7.3)
Pt p =
t p tc
fc (t)dt
(7.7.4)
If applying the Philip solution (and assuming the storm
intensity P is constant), using Equations (7.7.2) and (7.7.4), it
can be shown that:
tp =
S p2
Ks
1
+
2P(P K s )
2(P K s )
Sp
tc = t p
2(P
K
)
(7.7.5)
(7.7.6)
where the actual infiltration rate over the entire storm is then
given by:
P, t 0 t t p
f (t) = S
p
(t tc )1/2 + K s , t p t tr
(7.7.7)
Similarly, if using the Green-Ampt model, it can be
shown that:
tp =
Ks
( 0 )
P(P K s ) f s
( )
1
f
s
0
tc = t p
Pt p + f (s 0 )ln
Ks
Pt
+
p
f
s
0
(7.7.8)
(7.7.9)
241
where the actual infiltration rate over the entire storm is then
given by:
) P,
+
+
# ( )&
f (t) = *
f
s
0
(,
+ K s %1 +
%
Fc (t tc ) ('
+,
$
t0 t t p
t p t tr
(7.7.10)
P
it can be seen that for either model, the time to ponding and
compression time both go to zero. This simply reduces the
actual infiltration rate to the infiltration capacity rate (as
would be expected since it is consistent with the original
assumption of instantaneous ponding).
2. In the case where the precipitation rate is less than the
saturated hydraulic conductivity of the soil, then no ponding
will ever occur (since both infiltration capacity curves
asymptote to the hydraulic conductivity). In this case the
actual infiltration rate is simply given by: f(t) = P and the
above computation of ponding and compression times can be
skipped altogether. So it is important to first check whether P
> Ks when computing actual infiltration. Otherwise, using the
242
1/2
(
" 2(1) + 3 %+
S p = *(0.35 0.5(0.35))(11 mm/h) 500 mm $
'(1)
+
3
#
&,
)
34.7 mm h 1 2
! 2(1) + 3 $
f = #
&(500 mm) = 625 mm
1
+
3
"
%
"
%
(34.7 mm h 1/2 )2
11 mm/h
tp =
$1 +
'
2(30 mm/h)(30 11) mm/h # 2(30 11) mm/h &
(s 0 )=0.175;
= 1.36 h
2
tr
0
tc = 2.1 hr
1
"(30 mm/h)(2.1 h) + (109.375 mm)
11 mm/h #
11 mm/h
(109.375 mm)
(30 mm/h)(30 11) mm/h
= 2.1 hr
tp =
11)
mm/h
#
&
F(t = tr ) =
f (s 0 ) = 109.375 mm
Fc (t = tr tc = 7.1 h) = 187 mm
which is about 6% dierent than predicted by the Philip
solution.
244
S ECTION 8
MOD-WET Codes
Relevant functions based on concepts introduced in
this chapter include:
Brook-Corey retention curve model:
brooks_corey_PSI_and_K.m
Field capacity (using Brooks-Corey model):
field_capacity.m
Green-Ampt infiltration capacity model:
green_ampt.m
Philip infiltration capacity model:
philip.m
Soil sorptivity:
sorptivity.m
Time-compression approximation for actual infiltration:
TCA_infiltration.m
Permanent wilting point (using Brooks-Corey model):
wilting_point.m
245
S ECTION 9
Conceptual Questions
1. Define the meaning of porosity. How is it used to transform
between volumetric water content and relative saturation in
a soil?
2. What are the three end-member soil types in the SCS soil
triangle?
3. Which has a higher porosity, sand or clay? Which has a
higher field capacity? Which has a higher wilting point?
4. What are the two terms in the piezometric head? Describe
what each term physically represents.
5. What is the sign of matric head. Explain why.
6. How does the hydraulic conductivity and matric head
change as soil moisture decreases.
7. The Darcy flux dictates that the flow (per unit area) is
proportional to the piezometric head gradient. What soil
property is the proportionality constant in the Darcy
equation? What two factors generally cause that parameter
to vary in unsaturated soil (i.e. what does it depend on)?
8. When fully expanded, what are the two terms in the Darcy
flux in an unsaturated soil? Explain the physical meaning of
each term.
246
S ECTION 10
Sample Problems
Problem 7.1. After collecting multiple soil samples in the
field, you return to the laboratory and conduct soils tests to
classify the soil types present in your study basin. Sieve
analysis indicates that soil #1 has a composition of 34% sand
and 46% silt, while soil #2 has a composition of 81% sand
and 15% silt.
a) Classify both soils (i.e. determine the soil texture) using
the soil texture triangle. Based on the soil parameters listed in
Tables 7.1 and 7.2, plot the matric head and conductivity vs.
volumetric soil moisture curves for these two soils. Discuss the
relationship between soil moisture and matric head and
hydraulic conductivity. How do these relationships vary with
soil type?
Use the soil parameters from part a) for soil #1 to answer the
following:
b) Soil moisture sensors installed at location #1 (i.e. the
sampling site of soil #1) indicate that the volumetric soil
moisture profile at a given snapshot in time is:
(z) = 0 + az;
2 m z 0
247
q(z)
=
t
z
s = 50 cm; K s = 25 mm hr -1; b = 1
Chapter 8
Evaporation
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
14. Compute the latent and sensible heat flux using the
mass-transfer approach
S ECTION 2
Basics of Evapotranspiration
Evapotranspiration (ET) refers to all the processes by
which water in liquid phase at or near the Earths surface
becomes atmospheric water vapor. In particular it generally
refers to the flux of water from the surface into the
atmosphere. Evaporation usually describes the direct
vaporization of water from either open water surfaces or bare
soil. Transpiration refers to water loss from within the leaves
of plants. Together the two fluxes are referred to as
evapotranspiration. Sublimation is a special case whereby
vaporization occurs directly from ice. A more comprehensive
treatment of evaporation is given in Brutsaert (1982).
Evapotranspiration is important for several reasons
including:
1. The long-term water balance, which depends on
evapotranspiration (as well as precipitation and runo)
determines the available water for human use.
2. Much of the food supply is grown via irrigated agriculture
and ecient irrigation practices requires knowledge of
evapotranspiration.
F IGURE 8.1 Photo of a standard evaporation pan which provides an estimate of potential evaporation.
252
dS
= P E Q 0
dt
(8.2.1)
E = P Q
(8.2.2)
G = Rn LE H 0
(8.2.3)
LE = Rn H
E Ep
(8.2.4)
253
S ECTION 3
Evaporation can be thought of as a diusive process,
whereby a flux is driven by concentration gradients. In this
case, the concentration refers to the concentration of water
vapor molecules, which can be expressed in terms of the vapor
pressure. Over a liquid water surface, the air in immediate
contact with the liquid has a saturated vapor pressure. So
when the air above the surface layer is sub-saturated, there is
a gradient in vapor concentration, which would be expected to
correspond to a vapor flux away from the surface. It is
important to remember that the saturated vapor pressure is a
function solely of temperature, hence if there is a temperature
gradient between the surface (Ts) and the overlying air (Ta),
this may lead to a vapor pressure gradient even if both are
saturated. In Figure 8.3 the air in contact with the liquid has
a vapor pressure equal to es(Ts), while that above is at
ea<es(Ts). Under these conditions, there would be expected to
be an evaporative flux away from the surface. In general,
neither the surface or overlying air may be saturated, however
the gradients will still drive the net flux.
Speaking more generally of a diusive flux, the rate of
transfer (flux) of constituent X in direction z can be written
(using Ficks 1st Law):
Fz (X ) = DX
dC(X )
dz
(8.3.1)
E = Fz (vapor) = Dv
d v
dq
de
= Dv
= Dv
dz
dz
p dz
(8.3.2)
H = Fz (heat) = DH
d( c pTa )
dz
= c pDH
dTa
dz
(8.3.3)
V (z) =
u* " z d %
ln $
';
# z0 &
z > z0 + d
(8.3.4)
z 0 0.1h
d 0.7h
(8.3.5)
V
>0
z
(8.3.6)
So the diusivities in the flux equations for E and H are
related to turbulence. Here we will focus on the case over bare
soil/open water and will extend to vegetated surfaces in the
next section. If it is assumed that DH = DV, it can be shown
that the impact of turbulence can be embedded in E and H
models (using a finite dierence approximation):
2V22
(q 2 q1 )
2
( " z d %+ (V2 V1 )
** ln $ 2
'-z
) # 0 &,
(8.3.7)
256
(8.3.8)
2V22
(T2 T1 )
H c p
2
( " z d %+ (V2 V1 )
** ln $ 2
'-z
) # 0 &,
T1 = Tsurf ,
q1 = q surf , V2 = V ,
T2 = Ta ,
raN
h = m2 = (1 15RiB )1/2 ;
h = m = (1 5RiB )1;
RiB 0 (unstable)
0 RiB < 0.2 (stable)
(8.3.11)
where RiB is the bulk Richardson number, which is a nondimensional ratio of the thermal stability to wind shear:
q 2 = qa
( " z d %+
** ln $
'-) # z 0 &,
=
2V
(8.3.10)
ra = raN mh
RiB =
(g /Tsurf )T z
(V z)
(8.3.12)
(8.3.9)
RiB = 0;
m = h = 1
(8.3.13)
stability correction factors are less than 1.0 (which reduces the
resistance; enhanced turbulence), while for stable conditions
they are greater than 1.0 (which increases the resistance;
suppressed turbulence).
Using this notation the evaporation flux model can be
written as:
E=
(q surf qa )
ra
(8.3.14)
LE = Lv
(q surf qa )
ra
(8.3.15)
H = c p
(Tsurf Ta )
ra
(8.3.16)
Over well-watered surfaces (i.e. open water or moist soil),
the air in contact with the surface is expected to be saturated.
In this special case one can generally assume: qsurf = qs(Tsurf)
which is equivalent to E = Ep (i.e. potential evaporation).
Based on this, one model used for potential evaporation is:
Ep =
(q s (Tsurf ) qa )
ra
(8.3.17)
In the more general soil case, soil moisture limits qsurf (i.e. qsurf
< qs(Tsurf)), and E < Ep (i.e. actual evaporation is less than
the potential evaporation). In such cases, the soil moisture
impact is often modeled as:
q surf = ( )q s (Tsurf )
(8.3.18)
E = ( )E p
(8.3.19)
#
1,
fc
%
%%
wp
, wp fc
=$
% fc
wp
%
0,
wp
%&
(8.3.20)
259
q surf = q s (Tsurf ) =
es (Tsurf )
p
q surf = 0.622
36.2 mb
= 23.1 10 3 kg/kg
973 mb
E=
(q surf qa )
ra
es (Ta )
p
36.9 mb
= (0.69)(0.622)
= 16.3 g/kg
973 mb
qa = RH q s (Ta ) = RH
raN
p
RdT[1 + 0.608qa ]
97300 Pa
(287J/kg/K)(300.35K)[1 + 0.608(16.3 10 3 )]
= 1.12 kg/m 3
so that the evaporation rate is:
= 1.49 10 4 kg/m 2 /s
The latent heat flux is determined by multiplying the
above mass flux by the latent heat of vaporization:
260
E[mm/day] = (1.49 10
1 m3
kg/m /s)
1000 kg
2
1000 mm 86400 s
= 12.9 mm/day
1m
1 day
261
S ECTION 4
Transpiration
The models presented in the previous section are relevant
for open water surfaces or bare soil surfaces. However in many
applications, vegetation is present and is expected to play a
key role in ET. Transpiration by plants occurs via the
vascular system of the plant structure (Figure 8.4). It is
helpful to first start with the motivation underlying the
plants role in water loss from the surface. In this context,
photosynthesis is the key driver. Plants use photosynthesis to
build plant structures (leaves, stems, roots, etc.) in order to
grow. The photosynthetic reaction can be written as:
H 2O+CO2 CH 2O+O2
(8.4.1)
F IGURE 8.4 Schematic of key structures of a plant as relevant to transpiration flux (by Laurel Jules from
http://en.wikipedia.org/wiki/File:Transpiration_Overview.svg).
dioxide into the leaves or the flow of water vapor out of the
leaves. The reason for this regulation has to do with the fact
that photosynthesis requires all of the ingredients listed above.
When the ingredients are unavailable (i.e. nighttime when
PAR = 0, or when water is limited due to soil moisture
shortage, etc.), there is little benefit to having the stoma open
and needlessly lose the water (which will diuse out into the
air). So in this context the transpiration is essentially a water
loss term that happens at the leaves of the plant.
263
The evaporative loss occurs because the air inside the
stomatal cavity is essentially saturated (i.e. q = qs), while that
outside the leaf is generally sub-saturated (i.e. q < qs). This
gradient in vapor will drive a diusive flux as described in the
previous section. The guard cells attempt to regulate the
carbon dioxide influx in a way that minimizes environmental
stresses on the plant (chiefly water loss). This ability to
regulate transpiration will add an additional resistance to
the flux equations developed in the previous section. This can
be modeled in terms of an extra stomatal or vegetation
canopy resistance (Figure 8.7). The details of this
physiological behavior are complicated, but the way it is
handled in the modeling of ET is generally simplified. The
resistance for the entire vegetation canopy (i.e. all of the
leaves combined) is often scaled from the single-leaf scale
using the so-called leaf area index (LAI), which is defined as
the leaf area per unit area of ground surface. So dense
vegetation with a large leaf area to ground surface area may
have values of LAI of approximately 6.0 or greater. This
parameter has both a strong spatial and seasonal variation
(Movie 8.1).
F IGURE 8.7 Conceptualization of controls on evapotranspiration process using a resistance analogy approach (from
(8.4.2)
media.wiley.com/mrw_images/els/articles/a0003206/image_n/
nfgz002.gif).
earthobservatory.nasa.gov/Features/LAI/LAI3a.php).
E=
(q surf qa )
ra + rc
(8.4.3)
As mentioned above, the stomatal openings open/close
depending on environmental stresses being experienced by
the plant. So far these eects have not been included
explicitly in the model. The four stresses that are usually
accounted for are the incoming shortwave radiation (PAR), air
humidity, air temperature, and rootzone soil moisture. These
eects are usually modeled in the stomatal resistance
mentioned above as:
rsmin
rs =
(8.4.4)
fR fT fe f
s
rz
fR (R ) =
s
1.105Rs
1.007Rs + 104.4
0 R 1100 W m
-2
(8.4.5)
fe (e) = 1 0.000238e,
0 e 4200 Pa
(8.4.6)
Ta (40 Ta )1.18
fT (Ta ) =
,
a
690
0 Ta 40 C
(8.4.7)
#
1,
%
%%
wp
,
f (rz ) = $
rz
% fc
wp
%
0,
%&
fc
wp fc
(8.4.8)
wp
F IGURE 8.8 Illustrative example of environmental stress factors used in canopy resistance model.
low, the stress factor for soil moisture will approach zero,
which will increase the stomatal resistance, and ultimately
make the transpiration equal to zero. The specific forms of
these functions are generally empirically-based. One example
of such equations is the set given by (Dingman, 2008):
For the incoming shortwave function, the stress function
goes to zero (i.e. stoma close) when the radiation is zero (at
night) and it goes to 1.0 when the stoma are fully open during
a sunny day with a significant amount of PAR. The air
temperature stress function is generally a function with a
peak or plateau region, which represents the optimal
temperature conditions for the plant. If the temperature gets
too cold or too hot the stoma will close completely. The vapor
pressure deficit (VPD) function is 1.0 when the VPD is zero
and decays as the VPD increases. This is because evaporation
losses will generally be high when the VPD is high and low
when the VPD is low. For soil moisture, when the rootzone
266
E XAMPLE 8.4.1
A vegetated surface is adjacent to the bare soil
surface described in Example 8.3.1. The
vegetation consists of a grass surface with a
characteristic height of 5 cm, a leaf area index of
2.5 and a minimum stomatal resistance of 70 s/
m. The incident shortwave radiation at the
surface is measured to be 900 W m-2. Assume the
meteorological conditions and surface conditions
for the vegetated surface are the same as that in
Example 8.3.1. Compute the evaporation rate
from the vegetated surface.
The primary dierence for a vegetated surface is the
additional resistance due to stomatal control on vapor
loss. This requires computation of each of the
environmental stress factors in the stomatal resistance
function. Based on the saturated soil and meteorological
conditions and using functions illustrated in Figure 8.8
(Equations (8.4.5)-(8.4.8)), the stress factors can be
estimated as:
1.105(900 W m 2 )
fR (R ) =
= 0.98
2
s
1.007(900 W m ) + 104.4
267
= 219 W m 2
which is less than the bare soil evaporation.
rc =
rsmin
(70 s/m)
=
LAI(fR fT fe f ) (2.5)(0.98 0.80 0.73 1)
s
rz
= 49 s/m
The aerodynamic resistance also changes due to the
increased roughness of the vegetated surface, i.e.:
2
raN
= 0.88 10 4 kg/m 2 /s
268
S ECTION 5
Additional ET Models
c p (T2 T1 )
B=
Lv (q 2 q1 )
The mass-transfer models described in the last two
sections provide instantaneous (i.e. over 10-30 min.) models
for LE and H (over bare soil/open water or vegetated surface
respectively), but require measurements/estimates of Tsurf,
qsurf, Ta, qa, and V (as well as characterization of other surface
properties). Other methods/models have been developed for
cases where dierent measurements may be available.
The first approach discussed in this section is the socalled Energy Balance Bowen Ratio (EBBR) method. The
Bowen ratio is defined as the ratio of the sensible heat to the
latent heat flux: B = H/LE. Note that the Bowen ratio is a
parameter indicative of the state of the surface and its role in
the partitioning of energy. When the surface is very dry, most
of the available energy will go into sensible heat (high Bowen
ratio). When the surface is wet, most of the available energy
will go into latent heat flux (low Bowen ratio). From the
original diusion analogy models in Equations (8.3.2) and
(8.3.3), this can be written as:
c p T / z
B=
Lv q / z
(8.5.1)
(8.5.2)
Rn G = LE + H = LE(1 + B)
(8.5.3)
LE =
Rn G
1+B
(8.5.4)
H = B LE
(8.5.5)
E XAMPLE 8.5.1
A meteorological tower takes temperature and
humidity measurements at two levels (where level
1 is at 5 m and level 2 is 10 m above the surface)
and available energy measurements:
Air temperature (level 1) =27.2C
Air temperature (level 2) =24.2C
Air specific humidity (level 1) =12 g/kg
Air specific humidity (level 2) =11 g/kg
Available energy (Rn - G) = 500 W m-2
Estimate the latent and sensible heat fluxes from
the surface.
From the two-level measurements, the Bowen ratio can
be estimated as:
B=
LE =
500 Wm
= 227 Wm 2
1 + (1.2)
2
The second model described here is the so-called Penman
model. Penman (1948) combined the mass-transfer models
with surface energy balance to derive a model for potential
evaporation. The primary development of the model is as
follows:
1. Assume that the surface humidity in the mass-transfer
model for evaporation is equal to the saturated specific
humidity (this implies potential evaporation):
q surf = q s (Tsurf )
(8.5.6)
Keeping only the linear terms (i.e. the first two terms on the
right-hand-side) yields:
q s (Tsurf ) q s (Ta ) +
dq s
(T Ta )
dT T surf
(8.5.7)
where by definition:
dq s des Lv es
=
=
dT p dT p Rv T 2
(8.5.8)
H = (1.2)(227 Wm 2 ) = 272 Wm 2
both of which are away from the surface.
dq s des
=
=
dT p dT p
(8.5.9)
(Tsurf Ta ) =
Hra
c p
(8.5.10)
(Rn G) + v q
ra
LE p =
1+
psychrometric constant =
q = qs (Ta ) qa
(8.5.11)
pc p
Lv
The Penman model provides potential evaporation rates
given estimates of available energy, reference-level air
LE = ( )LE p
(8.5.12)
E XAMPLE 8.5.2
A reservoir has a nearby meteorological station
that takes the following measurements:
Air temperature =29.2C
Air relative humidity =76%
Surface air pressure = 980 mb
Windspeed (at 2 m height) = 10.3 m/s
Available energy (Rn - G) = 383.2 W m-2
Estimate the evaporation rate (in mm/day).
Assume the characteristic wave height on the
reservoir surface is 3 cm.
Note that in this case there is no data characterizing the
surface temperature. (If the surface temperature was
known, we could reasonably assume that the surface
specific humidity was equal to the saturated specific
humidity at that temperature.) Hence the mass-transfer
271
(Rn G) + v q
ra
LE p =
1+
%
((
$ 461 J/kg/K $ 273.16 K 302.35 K ''
= 41.57 mb
ea = (0.76) ( 41.57 mb) = 31.59 mb
q =
(0.622)
e =
(4157 3159 Pa) = 6.3103 kg/kg
p
(98000 Pa)
= 246.6 Pa/K
(98000 Pa)(1004 J/kg/K)
=
= 63.3 Pa/K
(0.622)(2.5 10 6 J/kg)
246.6 Pa/K
=
= 3.9
63.3 Pa/K
Putting all terms together yields the latent heat flux:
(3.9)(383.2 W/m 2 )
LE p =
+
1 + 3.9
(1.11 kg/m 3 )(2.5 10 6 J/kg)
(6.3 10 3 )
(25.6 s/m)
1 + 3.9
= 444 W/m 2
The equivalent evaporation rate is then given by:
444 W/m 2
1 m3
E[mm/day] =
= 15.3 mm/day
1m
1 day
272
Other models have followed directly from Penman. For
example over a large moist surface with minimal advection the
air will become saturated under continued evaporation, i.e.:
q 0
Under these conditions one can define the equilibrium
evaporation, i.e. that driven solely by energy considerations:
(R G)
n
LEe =
1+
(Rn G) + v q
ra
LE =
r
1+ + c
ra
(8.5.15)
Finally, it must be noted that the Penman model does
not consider vegetation cover. To account for the impact of
vegetation, the Penman-Monteith model is an extension which
follows the same recipe as above, but considers the canopy
The set of models described in this and the previous
sections, take various approaches to predict ET based on
various inputs. The model chosen depends on the surface type
as well as the measurements that are available. The specific
(8.5.13)
(Rn G)
LE p = LEe =
1+
(8.5.14)
273
(3.9)(383.2 W/m 2 )
LE =
+
50 s/m
1 + 3.9 +
20 s/m
(1.11 kg/m 3 )(2.5 10 6 J/kg)
(6.3 10 3 )
(20 s/m)
50 s/m
1 + 3.9 +
20 s/m
= 326 W/m 2
which is a lesser value than the potential rate from the
open-water surface in Example 8.5.2 due to the
additional stomatal control, despite having a lower ra.
274
S ECTION 6
MOD-WET Codes
Relevant functions based on concepts introduced in
this chapter include:
Aerodynamic resistance to turbulent transport:
aero_resistance.m
psychrometric_constant.m
S ECTION 7
Conceptual Questions
1. Describe the physical meaning of potential evaporation.
How does it generally compare to actual evaporation?
2. Is pan evaporation a measure of actual evaporation or
potential evaporation?
3. For a watershed mass balance, write the long-term average
evaporation as a function of other long-term average mass
fluxes.
4. Using an energy balance, write the long-term average
evaporation as a function of other long-term average energy
fluxes.
5. How do you convert between latent heat flux (in units of W
m-2) and evaporation mass flux (in units of kg m-2 s-1)? How
do you convert between an evaporation mass flux to units
of depth per unit time (i.e. mm/day)?
6. In the mass-transfer model for evaporation from bare soil,
what time-varying gradient is evaporation proportional to
(i.e. the driver of the flux). What gradient is sensible heat
flux proportional to?
7. What type of function describes the form of the mean
horizontal wind profile (i.e. as a function of z) in the neutral
surface layer?
276
18. Identify whether the EBBR, Penman, and PenmanMonteith models give actual or potential
evapotranspiration. Explain.
277
S ECTION 8
Sample Problems
Problem 8.1. Meteorological and surface measurements
collected at noon from a bare soil field site are:
Air temperature = 26.8C
Air relative humidity = 69%
Surface air pressure = 973 mb
Wind speed at 2 m height = 3.0 m s-1
Soil surface temperature = 26.2 C
The measurements are taken on the first sunny day after a
rainy period (you can assume the soil is saturated). The soil
roughness elements have a characteristic height (h) of 4 cm.
a) Using only the measurements/information provided, which
model could be used to estimate the instantaneous surface
evaporation? Which model/s cannot be used if you only have
the data provided? Explain why.
b) What additional information/input data would you need if
the soil were not saturated?
c) What additional information would you need if the surface
was vegetated?
d) In general, qsurf over bare soil can be modeled as:
q surf = q s (Tsurf );
0 ( ) 1
278
Problem 8.3.
a) Describe the meaning of potential evaporation. What
conditions are necessary for potential evaporation from a
surface? How does potential evaporation dier from actual
evaporation?
b) Evaporation from a reservoir can be obtained using the
Penman equation. If the available energy at the surface of the
reservoir is 220 W m-2, aerodynamic resistance is 30 s m-1, and
air specific humidity is 8 g/kg, what is the change in
evaporation (in mm/day) from the surface if the air
temperature is reduced from 20C to 15C as a cold air mass
passes over the reservoir? You may neglect any changes in net
radiation at the surface as a result of the dierent air
characteristics.
c) Qualitatively, explain how an increase in the surface
temperature would aect the evaporation rate from the
reservoir.
d) Qualitatively, explain how a lower humidity would impact
the evaporation rate from the reservoir.
Problem 8.4. Suppose the following measurements are made
over a bare soil surface (with a characteristic height of 2 cm
for surface roughness elements):
Reference-level (i.e. 2 meter height) wind-speed: 3 m/s
Reference-level air temperature: 25.9C
Reference-level air vapor pressure: 2217 Pa
Surface soil temperature: 37.9 C
CROP
STOMATAL
RESISTANCE (S/M)
VEGETATION
HEIGHT (CM)
30
30
60
100
280
Chapter 9
Groundwater
Flow
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
1. Write down Darcys Law for specific discharge in
saturated flow; show how it is a special case of the
expression used in unsaturated flow
2. Compute the average fluid velocity through an aquifer
based on the Darcy flux and soil porosity
3. List the key dierences between unconfined and confined
aquifers
4. Describe the dierent mechanisms by which water is
stored in unconfined and confined aquifers
5. Define the respective storage parameter relevant to
unconfined and confined aquifers
6. Write down the general 3D groundwater flow equation for
confined aquifers
7. State the key dependent variable that is being solved for
in any groundwater flow problem (i.e. once this variable is
known anything else can be solved for)
S ECTION 2
F IGURE 9.1 Illustration of various layers in subsurface including vadose zone, an unconfined and confined aquifer and an
aquiclude.
h=
p
+z
g
(9.2.1)
q = Kh
(9.2.2)
F IGURE 9.2 Typical aquifer configuration showing piezometric surface for unconfined and confined aquifers, which can be
measured via monitoring wells (from
artinaid.com/2013/04/types-of-aquifers) .
qx K x
h
,
x
qy = K y
h
,
y
q z = K z
h
z
(9.2.3)
v=
Q
q
=
sA s
(9.2.4)
284
F IGURE 9.3 Map showing the large-scale groundwater aquifers across the U.S. (from USGS;
pubs.usgs.gov/ha/ha730/ch_a/gif/A004_us.gif).
E XAMPLE 9.2.1
The
q = Kh
where the gradient must be estimated from the head
distribution. The head gradient is mathematically in the
direction of steepest ascent, which graphically is
perpendicular to head contours. The negative sign in the
Darcy flux indicates that flow is in the down-gradient
direction. The down gradient direction at x = 3 km, y =
2 km is shown in the annotated figure below. The
gradient can be estimated using a finite dierence:
x =3km,y=2km
h
5m
= 0.01
L x =3km,y=2km 500 m
a) On the map,
sketch
the flowof
direction
water starting
from point
The
pattern
head of
shows
a general
peak(3km,
near3km)
theto the edge of the island. in
(2 pts.)
b) Where is the Darcy flux the maximum and what is its approximate value? Explain your reasoning for
choosing where the flux is largest. (4 pts.)
The steady-state head field for a confined aquifer under a 5 km x 5 km square island (with a
homogeneous and isotropic hydraulicqconductivity
of 1 m/day) is plotted in the contour map below. The
0.01 m/d
v
=
=
= the
0.029
datum is such that the head on all sides of the boundary of
islandm/d
is 0 m.
0.35
+"
a) On the map, sketch the flow direction of water starting from point (3km, 3km) to the edge of the island.
Note that the Darcy flux and velocity will vary in space.
(2 pts.)
b) Where is the Darcy flux the maximum and what is its approximate value? Explain your reasoning for
contours
closest
(i.e. largest gradient), which occurs
choosing where
the flux isare
largest.
(4 pts.)
drainage ditches which allow water to be removed from the soil through lateral groundwater flow. The
water in the drainage ditches is maintained at a depth H0. The irrigation rate I is applied at the top of the
soil profile. The water which is not evaporated via transpiration by the crops recharges (at a rate R) the
unconfined aquifer below (which is bounded by a horizontal impermeable layer at depth D). The crops
being grown by the farmer have a rooting depth of d and have an average evapotranspiration rate E. The
farmer must decide the irrigation rate such that the rootzone remains unsaturated (i.e. the water table
287
S ECTION 3
Ss
(9.3.3)
h
= (q)
t
The groundwater flow equation is analogous to the flow
equation for unsaturated flow, but under saturated conditions.
Here we will start with the general 3D flow equation using
mass balance and then simplify it to 1D and 2D versions as
special cases.
The statement of mass balance can be written as (same
starting point as in unsaturated flow):
M
= ( wq dx dy dz)
t
(9.3.1)
M
h
= wSs dx dy dz
t
t
(9.3.2)
Ss
Ss
q
q
q
h
= x y z
t
x
y
z
h
h
h
=
K
+
K
+
t x x x y y y z
(9.3.4)
h
K
z z
(9.3.5)
K x K y K z
K i K i (x,y,z)
=
=
=0
x
y
z
Second, an aquifer is said to be isotropic if there is no
directionality to the conductivity field, i.e.:
K x = Ky = K z = K
Finally, if the flow is in steady-state, then there are no head
variations in time, i.e.:
h
=0
t
Given these special conditions, one can derive simplified forms
of the 3D flow equation. For example, for a homogeneous/
isotropic aquifer, we can write the governing equation as:
2 h 2 h 2 h
h
Ss
=K 2 + 2 + 2
t
y
z
x
(9.3.6)
2 h 2 h 2 h
2
2 + 2 + 2=h =0
y
z
x
(9.3.7)
In many natural groundwater systems, the flow is largely
horizontal (i.e. in the x- and y-directions) with little to no
vertical flow. As such, in many cases solving the 2D GW flow
equation is sucient. Here we develop the 2D flow equation
from the general 3D equation shown above. The usual
mechanism for doing so if via the so-called Dupuit
approximation, which essentially involves integrating the 3D
equation with respect to z, assuming no vertical flow inside
the aquifer (qz = 0) and that the aquifer has a horizontal lower
boundary. Together these mean that h=h(x, y, t). Note that
the lack of dependence on z implies that h is constant in any
vertical plane. This is consistent with the no vertical flow
assumption (i.e. if h varied with z then there would be a
gradient and therefore a Darcy flux). In such cases, as depth
increases into the aquifer, there are one-to-one tradeos
between the pressure and elevation head terms so that the
head is constant in the vertical. In other words, at the surface
the pressure is zero so the piezometric head is equal to the
elevation head. As one moves down vertically into the aquifer
the pressure increases such that the pressure head increases by
the exact amount that the elevation head decreases.
To develop the 2D GW flow equation we integrate
Equation (9.3.5) from z = 0 to z = H, where H = h (i.e. the
water table) in the case of an unconfined aquifer and H = b
(i.e. the aquifer thickness) in the case of an confined aquifer:
H
H
h
h
h
S
dz
=
K
+
K
0 s t
0 x x x y y y + z
h
K z z dz
289
Ss
0
h
h
dz =
S dz
t
t 0 s
h
[unconfined aquifers]
t
h
h
=Ss b
= S ; [confined aquifers]
t
t
= Sy
K
dz
=
0 x x x x
h
h
0 K x x dz = x K x x z
0 z
h
K
h
[unconfined aquifer]
x x x
h
h
K
b
=
T
[confined aquifer]
x x x x x x
h
h
h
0 y Ky y dz = y 0 Ky y dz = y Ky y z
h
=
K
h
[unconfined aquifer]
y y y
h
h
K
b
=
T
[confined aquifer]
y y y y y y
h
h
K
z z dz = K z z
= Kz
h
h
Kz
=R0
z H
z 0
h
h h
=
T
+
T
+R
t x x x y y y
(9.3.8)
Sy
h
h
h
=
K
h
+
K
h
+R
t x x x y y y
(9.3.9)
h 1 2
=
(h );
x 2 x
h 1 2
=
(h )
y 2 y
290
Sy
S 2
h
h
Sy
= y
(h )
t
h0 2h0 t
Sy h 2
1 h 2
1 h 2
=
K
+
K
+R
2h0 t
x x 2 x y y 2 y
d 2 R
0= 2 +
; subject to BCs
C2
dx
(9.3.12)
(9.3.10)
2 2 R
C1
= 2+ 2+
t x
y C 2
(9.3.11)
= h 2 2 [unconfined]; = h [confined]
C 1 = Sy (Kh0 ) [unconfined]; = S T [confined]
C 2 = K [unconfined]; = T [confined]
In many applications, especially those types of problems that
can be solved analytically, we will be dealing with steady-state
and/or 1D problems. In the case of steady-1D problems, we
can use the above unified equation as a starting point. A 1D
flow problem involves flow only in one direction (i.e. the x-
(x) =
R 2
x + a 0x + a 1
2C 2
(9.3.13)
291
E XAMPLE 9.3.1
b) The steady-state irrigated system yields a onedimensional unconfined aquifer which is governed by
(from Equation (9.3.11)):
d 2h 2
2R
=
K
dx 2
Integrating twice yields:
h 2(x) =
R 2
x + a 0x + a 1
K
h 2(x = 0) =
R 2
(0) + a 0 (0) + a1 = H 02
K
a1 = H 02
R
(L)2 + a 0 (L) + H 02 = H 02
K
RL
a0 =
K
h 2(x = L) =
R " L % RL " L %
2
2
2
h (x = L 2) = $ ' +
$ ' + H 0 = (D d)
K #2&
K #2&
RL2
+ H 02 = (D d)2
4K
(I E)L2
+ H 02 = (D d)2
4K
where solving for the irrigation rate yields:
I =E+
4K
(D d)2 H 02
2
L
h 2(x) =
R 2 RL
x +
x + H 02
K
K
293
S ECTION 4
Groundwater Flow to
Pumping Wells
As mentioned above, groundwater is a primary source of
freshwater. To extract this supply from an aquifer generally
requires a pumping well (Figure 9.4) whereby the well taps
into the existing (natural) flow conditions. A pumping well
generally consists of a lined bore hole (i.e. with a pipe) that
has a screened end within the aquifer and a pump on the
other end of the pipe. Turning on the pump takes water out of
the aquifer and by mass balance induces flow toward the well
as will be described in more detail below. Note that a
monitoring well is quite dierent, as it does not have a pump
and is generally used solely for monitoring the piezometric
surface and/or constituents in the water (i.e. is passive rather
than an active part in the flow system). A pumping well may
be drilled into a near-surface unconfined aquifer or more
deeply into a confined aquifer. In this section we focus on the
hydraulics associated with flow toward pumping wells.
For mathematical tractability we will focus on a single
pumping well in an aquifer of infinite horizontal extent and for
conditions where the original piezometric surface (i.e. prior to
pumping) was horizontal. Extensions to multiple pumps and
more realistic cases can be built up from these single-well
solutions via superposition as discussed in the next section.
The basic schematic for such single well configurations is
F IGURE 9.4 Illustration of impact of pumping on the regional natural groundwater system.
around a single pumping well in a horizontally infinite confined aquifer (from Mays, 2005).
s(r,t) = h0 h(r,t)
(9.4.1)
around a single pumping well in a horizontally infinite unconfined aquifer (from Mays, 2005).
1 R
C1
=
r
+
t r r r C 2
(9.4.2)
C1
1
=
r
t r r r
(9.4.3)
S h 1 h
=
r
T t r r r
(9.4.4)
s(r,t) = h0 h(r,t) =
Q
W (u)
4T
(9.4.5)
e u
W (u) = du
u
u
(9.4.6)
Sr 2
u=
4Tt
(9.4.7)
0.219
0.049
0.013
3.8
10-3
1.1
10-3
3.6
10-1
1.2
10-4
3.8
10-5
1.2
10-5
10-1
1.82
1.22
0.91
0.7
0.56
0.45
0.37
0.31
0.26
10-2
4.04
3.35
2.96
2.68
2.47
2.3
2.15
2.03
1.92
10-3
6.33
5.64
5.23
4.95
4.73
4.54
4.39
4.26
4.14
10-4
8.63
7.94
7.53
7.25
7.02
6.84
6.69
6.55
6.44
10-5
10.94
10.24
9.84
9.55
9.33
9.14
8.99
8.86
8.74
10-6
13.24
12.55
12.14
11.85
11.63
11.45
11.29
11.16
11.04
10-7
15.54
14.85
14.44
14.15
13.93
13.75
13.6
13.46
13.34
10-8
17.84
17.15
16.74
16.46
16.23
16.05
15.9
15.76
15.65
10-9
20.15
19.45
19.05
18.76
18.54
18.35
18.2
18.07
17.95
10-10
22.45
21.76
21.35
21.06
20.84
20.66
20.5
20.37
20.25
10-11
24.75
24.06
23.65
23.36
23.14
22.96
22.81
22.67
22.55
10-12
27.05
26.36
25.96
25.67
25.44
25.26
25.11
24.97
24.86
10-13
29.36
28.66
28.26
27.97
27.75
27.56
27.41
27.28
27.16
10-14
31.66
30.97
30.56
30.27
30.05
29.87
29.71
29.58
29.46
10-15
33.96
33.27
32.86
32.58
32.35
32.17
32.02
31.88
31.76
296
S h 1 h
=
r
=0
T t r r r
(9.4.9)
d dh
r
=0
dr dr
h(r) = c 0 ln(r) + c1
(9.4.10)
dh
Q = 2 rb K
dr
(9.4.11)
c
Q
Q
Q = 2 rb K 0 c 0 =
=
2 Kb 2T
r
(9.4.12)
h(r1 ) = h1 =
Q
Q
ln(r1 ) + c1 c1 = h1
ln(r1 )
2T
2T
(9.4.13)
!r $
Q
h(r) =
ln # & + h
2T " r1 % 1
(9.4.14)
h(r) =
r
Q
ln + h0 ;
2T R
rw r R
(9.4.15)
r
Q
s(r) = h0 h(r) =
ln ;
2T R
rw r R
(9.4.16)
T =
Q ln[r2 r 1 ]
2 h2 h1
(9.4.17)
S h 2 1 1 h 2
=
r
=0
Kh0 t
r r 2 r
(9.4.18)
d 1 dh 2
r
=0
dr 2 dr
which can be integrated (twice) and using the same procedure
for identifying integration constants as shown above (for the
confined aquifer case) yields:
298
r
Q
h (r) =
ln + h12 ;
K r1
2
rw r R
(9.4.19)
s(r) = h0 h 2(r);
rw r R
(9.4.20)
K =
Q ln[r2 r1 ]
(h22 h12 )
(9.4.21)
So to summarize, we have defined the solutions to the
groundwater flow problem for three single-well infinite aquifer
cases: i) transient flow in a confined aquifer, ii) steady-state
flow in a confined aquifer, and iii) steady-state flow in an
unconfined aquifer. Additionally, while we do not have an
analytical solution for the fourth case of transient flow in an
unconfined aquifer, the Theis solution can be used if the
unconfined aquifer is relatively thick compared to the
drawdown. Finally, it should be noted that the above
development defined a positive pumping rate as one that
extracts flow from the aquifer. These cases yield a cone of
depression based on the given pumping rate. All of the above
(0.0001)(25 m)2
u=
= 0.00003
4(50 m 2 /day)(10 day)
From Table 9.1 the Well function at this value is: W(u)
= 9.84. The drawdown is:
(10 m 3 /day)
s(r,t) =
(9.84) = 0.16 m
2
4 (50 m /day)
299
E XAMPLE 9.4.2
A single well pumps from a horizontally infinite
confined aquifer at a rate of 5 m3/day until
steady-state is reached. The well is of radius 5 cm
and shows a drawdown of 1 m. At a monitoring
well located 100 meters from the pumping well
the drawdown is 0.1 m. Estimate the
transmissivity of the aquifer. Based on the
estimate, what is the expected steady-state
drawdown at a radius of 50 meters?
The two dierent head measurements can be used to
estimate the transmissivity:
" 50 m %
(5 m 3 /day)
s(r) =
ln
$
' = 0.19 m
2 (6.1 m 2 /day) # 215 m &
300
S ECTION 5
To start, the easiest example of superposition is the case
of a confined aquifer of infinite extent with multiple wells.
Conceptually, we know that a single well generates a
symmetric drawdown around the well. Multiple wells will have
multiple cones of depression that can form a complicated
head/drawdown field. The linearity of the single well solution
(with respect to h) allows for the determination of the total
drawdown as simply the summation of drawdown from each
well in the well field. (Note for an unconfined aquifer,
superposition does not strictly apply since the governing
equation is not linear in h). So the actual drawdown (s) is
given by:
(9.5.1)
where this assumes there are N wells. The x- and ycoordinates are used rather than r simply because each well
will have its own radial coordinate system. In practice, the
drawdown at a given location can be determined by finding
the distance from that point to each pumping well and using
that distance as the argument in the single-well solution (i.e.
in Equations (9.4.5) or (9.4.16)) to get each of the individual
drawdown predictions in Equation (9.5.1).
Figure 9.7 shows an example of a case with 2 wells with
pumping rates Q1 and Q2. In general, the pumping rates may
dier. Even for this relatively simple case, the composite
drawdown curve is complicated. For example, there is a local
maxima in between the two wells, with local minima at the
two wells. The local maxima corresponds to the location
where the head gradient is zero (implying no flow). It
301
E XAMPLE 9.5.1
It is required to de-water (i.e. lower the water
table) below a construction site that 80 m by 80
m in area as shown in plan below.
F IGURE 9.7 Drawdown solution for multiple wells via superposition (solid line) from single well solutions (dashed line).
302
Q =
2Ts(r)
2 (15 m/day)(100 m))(0.375 m)
=
ln(r / R")
ln(56.7 m / 600 m)
= 1498 m 3 /day
Building upon the example of multiple wells, we can
apply superposition to solve problems that involve non-infinite
aquifers. In this context, non-infinite means there is some
feature which impacts the flow to the well that is within its
radius of influence. There are generally two examples of this:
no-flux and fixed head BCs.
The first case we will discuss involves a no-flux boundary
within the radius of influence of the well. An example of such
a system is shown in the top panel of Figure 9.8.
303
(9.5.2)
F IGURE 9.9 Contour map of a single well solution in the presence of an impervious boundary. Note this is obviously incorrect as the gradient is non-zero at the boundary.
the two solutions over the entire domain provides the actual
cone of depression for this case (shown with the solid line).
Another example is shown in Figures 9.9 - 9.11 in terms
of the drawdown fields. Figure 9.9 illustrates that a single well
solution for the real system violates the no-flux BC. Figure
9.10 shows the contours for the drawdown fields predicted by
each of the (real and image) well drawdown solutions. Figure
9.11 shows the summation of the two drawdown fields, which
305
The other type of non-infinite aquifer condition that
often must be dealt with is that of a stream boundary. For
simplicity, the stream is often modeled as a fixed-head BC
(Figure 9.12). The primary consequences of the fixed-head
boundary is that the boundary acts as a water source that
would not exist in an infinite aquifer, so that less water is
required to be supplied by the rest of the aquifer leading to
non-symmetric and generally less drawdown. Superposition
can again be used to solve this problem. We again choose to
add an image well at the same distance as the real well, but
on the other side of the fixed-head boundary. However, rather
than making the image well a pumping well (to represent the
water loss and no-flux boundary condition provided by the
impervious boundary) we make it a recharge well (to represent
the extra source of water provided by the stream boundary).
By symmetry, since the real and image wells are the same
distance and the same magnitude of pumping rate, the
amount of drawdown predicted by the single-well solution for
the real well is exactly balanced by the same amount of
build-up as a result of recharge. Therefore, everywhere along
the fixed-head boundary the drawdown is zero. This is shown
in the bottom panel of Figure 9.12. Moreover we see that
generally there is less overall drawdown everywhere in the
aquifer due to the extra source provided by the stream.
306
Another example analogous to that shown in Figures
9.9-9.11 for the no-flux BC is shown in Figure 9.13 for the
fixed-head BC. The setup is exactly the same as before except
that the image well involves recharge so that there is a sign
change in the image well drawdown contours. In this case the
drawdown from the image well is negative (build-up). By
summing up the two single-well solutions, we get the result
shown in Figure 9.13. As can be seen, the drawdown is
head boundary using the image well approach. The top panel is
the real system and the bottom panel is the conceptual model
using superposition with a recharge image well to mimic the
no-flow boundary The dashed lines in the bottom image are the
single well solutions (one discharge and one recharge) which
when superposed yield the real cone of depression (adapted from
Mays, 2005).
307
E XAMPLE 9.4.2
The drawdown due to the pumping well if there were no
boundary would be (using the Theis solution):
Q
W (u1 );
4 T
Sr12
(10 4 )(50 m)2
u1 =
=
= 3.3 10 5
2
4Tt 4(950 m /day)(2 day)
W (u1 ) = 9.76
s1 = s(r = 50 m, t = 2 days) =
s1 =
86400 s
)
1d
9.76 = 4.24 m
4 (950 m 2 /day)
(0.06 m 3 s 1
Q
W (u2 );
4 T
Sr12
(10 4 )(150 m)2
u1 =
=
= 2.96 10 4
2
4Tt 4(950 m /day)(2 day)
W (u1 ) = 7.55
s2 =
86400 s
)
1d
7.55 = 3.28 m
2
4 (950 m /day)
(0.06 m 3 s 1
308
These examples show some of the power of the
superposition principle to solve relatively complex problems in
a relatively easy way. Even these examples/concepts of
superposition could be used to build up solutions for more
complicated problems. For example, if a pumping well is in a
river valley with a horizontal head distribution prior to
pumping, the aquifer may have a stream BC condition on one
side and a impervious boundary on the other. Depending on
the radius of influence of the well (i.e. given its pumping rate)
both boundaries may impact the flow patterns. In such a case
one may use two image wells and have an additive solution
with three terms instead of two. Or even more complicated
yet, the original head distribution (prior to pumping) may not
be horizontal, but itself the solution of a 1D flow problem
between the two BCs with recharge. In this case the solution
could be the summation of the original piezometric head
distribution combined with the drawdown from the real and
imaginary wells. Another example would be the case with
multiple real wells with a no-flux or fixed-head BC (or both).
In this case each real well would have a corresponding image
well. In these types of problems the hardest part is breaking
down the real problem into its constituent parts (each of
which has an analytical solution), with the easy part being
determining the solution via superposition of the individual
solutions.
As groundwater flow problems become more complex,
one may instead need to solve the problem numerically. In
such a case, the aquifer is discretized into small dierential
units and the GW flow equations are solved on the finite
309
310
S ECTION 6
Conceptual Questions
1. Name the two types of groundwater aquifers. Describe how
they dier.
2. What type of aquifer generally occurs closest to the land
surface (i.e. is shallowest in the subsurface)?
3. In an unconfined aquifer, what physical boundary does the
piezometric surface correspond to?
4. In an unconfined aquifer, is the piezometric surface above
or below the upper boundary of the aquifer? Explain.
311
S ECTION 7
Sample Problems
t(r) =
dr
v
!
314
Chapter 10
Runoff and
Streamflow
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
S ECTION 2
F IGURE 10.2 Flow accumulation patterns across terrain derived from DEM in Figure 10.1. Areas of high flow accumulation (red) correspond to stream channels while areas of low
flow accumulation area (dark blue) represent watershed
boundaries.
319
S ECTION 3
Runo Generation
Mechanisms
Runo is generated as a result of a water flux being
applied to the surface. This flux can either be in the form of a
rainstorm or snowmelt (Figure 10.4). Storm events can vary in
terms of intensity and duration, both of which should impact
the runo. For conceptual purposes, one can think of
F IGURE 10.4 Conceptual picture of runoff resulting from either rainfall or snowmelt.
F IGURE 10.5 Components of various runoff and related processes contributing to streamflow.
320
P, t 0 t t p
f (t) =
fc (t tc ), t p t tr
M OVIE 10.2 Animation of infiltration excess runoff mechanism at a given location on a hillslope.
(10.3.1)
t0
ie
(t) =
tr
tr
[P f (t)]dt = [P f (t t )]dt
t0
tp
(10.3.2)
322
E XAMPLE 10.3.1
Compute the cumulative infiltration excess runo
for the case described in Example 7.7.1.
The cumulative infiltration excess runo can be
computed via the integral in Equation (10.3.2).
Alternatively if the cumulative infiltration is already
computed (as in Example 7.7.1), then the runo is
simply the dierence between the cumulative rainfall and
the cumulative infiltration. So for the Philip solution the
estimated runo would be:
tr
t0
ie
(t) =
tr
t0
= 63 mm
t0
ie
(t) =
tr
t0
= 53 mm
The second overland flow generation mechanism is
referred to as saturation excess runo. It is sometimes referred
to as Dunne runo in reference to the work that first
identified the mechanism (Dunne, 1975). This mechanism
occurs when the groundwater table (due to recharge from
above) saturates the soil from below. Movie 10.5 shows a
conceptual representation of this phenomenon. In this case all
323
F IGURE 10.7 Illustration of variable contributing areas growing and contracting during and after a storm. These areas are
the locations expected to generate saturation excess runoff
(adapted from Dingman, 2008).
324
The final subsurface runo mechanism is generally
referred to as baseflow. It is simply lateral groundwater flow
into the stream channel network (Figures 10.8 and 10.9). The
flow can be from either unconfined or confined aquifers. In
Chapter 9, streams were treated as boundary conditions in the
groundwater flow problem. As seen in those examples the flux
is generally into the stream channels, which is exactly what
baseflow is. In the steady-state groundwater flow problems
325
326
S ECTION 4
Streamflow Hydrographs
The mechanisms described in the previous section
ultimately lead to runo that contributes to streamflow. One
could envision measuring streamflow at a given point in the
stream network. If measured, the amount of flow crossing that
point (i.e. cubic meters per second) would simply be the result
of all upstream flow processes (both runo to the closest
downstream channel location and flow down the stream to the
measurement point) as a function of time. The streamflow
crossing a given point over time is generally referred to as a
hydrograph, where one can conceptualize the total flow (Q)
as:
Baseflow corresponds to groundwater flow, while other components together are often referred to as stormflow. The fractional components depend on the storm and basin characteristics (adapted from Mays, 2005).
(10.4.1)
F IGURE 10.13 Impact of rainfall timing on stormflow hydrograph (adapted from Mays, 2005).
highest intensity rainfall occurs late in the storm and the right
panel shows a case where the highest intensity rainfall occurs
earlier in the storm. In the former case the rising limb is less
steep with a later peak.
Figure 10.14 illustrates how particular basin
characteristics can impact the hypothetical storm hydrograph.
The top panels show the response for basins of diering slope;
steeper basins will have hydrographs that respond and recede
more quickly. The bottom panels show the response for basins
of diering roughness; basins with higher roughness will
respond and recede less quickly. Another factor that is not
shown is the level of imperviousness, which is especially
relevant in urbanized basins. Basins with high levels of
urbanization will have more infiltration excess runo and
therefore the hydrograph will respond and generally recede
more quickly. Additionally, the available storage in the basin
can have an impact; basins with little storage will respond and
recede more quickly.
328
F IGURE 10.14 Impact of basin slope and roughness on stormflow hydrograph (adapted from Mays, 2005).
329
330
S ECTION 5
Rainfall-Runo Modeling:
Unit Hydrograph
The modeling of streamflow for a basin is not only of
implicit hydrologic interest, but has deep roots in the
engineering practice of flood forecasting and the need for
design flood estimation. In either of these applications, the
peak flow and the timing of the flow are key input parameters
to the design/forecast. Models that can be used for predicting
hydrographs span a large range of complexity. A detailed
treatment of modeling runo is given in Beven (2012).
In this section we will focus on a classical method used
in hydrology called the unit hydrograph (UH) method, which
is specifically designed to predict the stormflow hydrograph
for a given storm. Some of the underlying choices in this
method include: treating the basin as a lumped response unit,
using historical (archived) data to develop the model response
function, and making some simplifying assumptions including
using a systems (black-box) response approach. The primary
benefit of these choices is that the model is very
computationally ecient. In the next section we will focus on
a more physically-based approach that is becoming more
commonly applied in research and/or practice.
The UH approach is an empirically-based method to
predict stormflow response for a given storm. Conceptually,
F IGURE 10.16 Illustrative example of a streamflow hydrograph that can be used to derive a unit hydrograph.
333
N days = 1.21A0.2 ;
(10.5.1)
F IGURE 10.19 Plot of the stormflow hydrograph with baseflow removed from the original streamflow hydrograph.
Pe =
1
(Q B)dt
A
F IGURE 10.20 Conceptual picture of the process of converting total volume of excess runoff to the equivalent amount of
uniform runoff depth over the basin.
(10.5.2)
ht =
Qt Bt
Pe
(10.5.3)
where the units of the UH are flow units per unit depth of
stormflow, i.e. m3/s/cm. The UH simply represents the
response function that would occur if 1 cm of eective rainfall
occurred over D hours.
335
is converted to the unit hydrograph by dividing by the equivalent depth of runoff. Note the units of the unit hydrograph are
volumetric flow per unit stormflow depth (e.g. m3/s/cm).
336
E XAMPLE 10.5.1
The hydrograph shown below was measured at
the outlet of a 1 km2 watershed in response to a
storm event that corresponded to 2 hours of
excess precipitation. Construct the 2-hour UH
from this hydrograph. For simplicity assume the
flow at the start of the hydrograph is
representative of the baseflow.
TIME
(HOURS)
FLOW
(m3/s)
0.05
0.10
0.15
0.13
12
0.10
15
0.08
18
0.05
TIME
(HOURS)
STORMFLOW
(m3/s)
0.0
0.05
0.10
0.08
12
0.05
15
0.03
18
0.0
UH (m3/s/cm)
0.0
0.074
0.148
0.118
12
0.074
15
0.044
18
0.0
The end result of the UH construction process is the
determination of the unit response of the basin (i.e. to 1 cm of
excess rainfall) over a D-hour duration. Once constructed, the
UH can be used in applications, including the prediction of
response to other D-hour events of varying intensity as well as
response to storms of dierent excess precipitation intensity
and duration. This is done via superposition based on the
assumed linear response of the system. Three examples of
applications of the UH method are presented here including:
prediction of stormflow from a dierent D-hour storm,
construction of an (nD)-hour UH from a D-hour UH (where n
is an integer multiple), and prediction of response from a more
complicated excess precipitation event.
The simplest application of the UH method involves
using a D-hour UH to predict the response of a dierent Dhour event. For example, suppose for design or flood
prediction purposes you want to predict the resulting
hydrograph for a 6-hour eective rainfall event with Pe of
excess precipitation. The entire basis of the UH approach is
that the system behaves linearly so that the shape should be
invariant, but the magnitude can be scaled. So in this case, if
you have a 6-hour UH, you would simply scale it by Pe:
Qt = Pe ht
(10.5.4)
339
from two 1-hour UHs). This could then be used to predict the
response to any 2-hour event.
The final application discussed here is the case where
there is an excess precipitation event of longer duration than
the available UH and of varying excess precipitation
throughout the event. An example is that shown in Figure
10.25, which shows an excess precipitation time series over a
24-hour period. For the purposes of illustration, we will
assume we have already constructed a 6-hour UH for this
basin. One approach would be to construct the 24-hour UH
from four lagged 6-hour UHs as described above and then
341
E XAMPLE 10.5.2
TIME
(HR)
UH FLOW
(m3/s/cm)
10
342
While in the UH construction the excess precipitation is
determined directly from the historical hydrograph,
applications to other events require some estimate of Pe. In
reality it is determined by the runo generation mechanisms
occurring in the basin, which depend on the precipitation
intensity, soil type, and antecedent conditions (i.e. initial soil
moisture). In practice it is often estimated empirically. One
commonly applied method is the so-called SCS method
developed originally by the U.S. Soil Conservation Service
(now known as the Natural Resources Conservation Service
(NRCS)). A schematic of the assumed processes is shown in
Figure 10.28, where Ia is the initial abstraction (i.e. where all
water infiltrates), Fa is the continuing abstraction (potentially
decaying infiltration rate), P is the total precipitation and Pe
is the excess precipitation. Based on empirical evidence from
many small experimental watersheds, the excess precipitation
can be estimated by:
(P 0.2S)2
Pe =
P + 0.8S
(10.5.5)
S=
1000
10
CN
(10.5.6)
344
S ECTION 6
Rainfall-Runo Modeling:
Physically-based Models
The unit hydrograph approach to rainfall-runo
modeling described in the previous section has been and
continues to be widely used. If one needs only basin outlet
runo predictions and the underlying assumptions of the UH
method are reasonably valid for the basin of interest, then
UH-based predictions may be sucient (and often are, e.g. for
design purposes). The empirical nature of the approach
however includes limitations.
The primary alternative to empirical modeling is
modeling using physically-based approaches. This simply
means that the processes within the basin are modeled using
the physics that have been the primary basis of this book. As
has been made clear in earlier chapters, the primary drivers of
watershed processes are precipitation and net radiation. Given
these inputs, a set of processes ensue that include infiltration,
evaporation, unsaturated zone moisture redistribution,
recharge, groundwater flow, and runo. Each of these is
governed by physical processes that can be expressed in terms
of models. If tied together into an integrated unit, the model
becomes a physically-based hydrologic watershed model. One
should be plainly aware of the tradeos between models.
Physically-based models have the potential for a more robust
modeling framework (that may include modeling of interior
F IGURE 10.33 TIN-based distributed schematic representaF IGURE 10.32 Distributed schematic representation of basin
zone storage and net radiation into: surface sensible and latent
heat fluxes, ground heat flux and soil energy storage.
Subsurface processes that are modeled may include: one
dimensional flow in the unsaturated zone (including
redistribution and groundwater recharge), and groundwater
mass balance. Lateral fluxes modeled at the pixel-scale include
overland flow (infiltration excess and saturation excess) and
subsurface runo (interflow and baseflow). These runo
generation mechanisms yield an outflow hydrograph for each
F IGURE 10.35 Illustration of saturation excess runoff generation predicted by the tRIBS distributed model (from Enrique Vivoni personal communication).
S ECTION 7
0=
d
dV + CS V dA
dt CV
(10.7.1)
inlet
V dA = (Q + qdx)
(10.7.2)
351
V
dA
=
(Q
+
dx)
(10.7.3)
x
outlet
where the second term on the right-hand-side represents the
increment in flow due to lateral inflow and/or any storage
changes. Finally, the rate of change of mass stored in the
control volume is given by:
dV
=
( Adx)
dt CV
t
(10.7.4)
( Adx)
Q
(Q + qdx) + (Q +
dx) = 0
t
x
(10.7.5)
A Q
+
q = 0
t
x
(10.7.6)
F =
d
V dV + CS V V dA
dt CV
(10.7.7)
F = F
+ Ff + Fe + Fp
(10.7.8)
y
dx
x
(10.7.9)
352
V V dA = (VQ + v qdx);
inlet
(10.7.10)
momentum coefficient
V
dA
=
VQ
+
dx)
(10.7.11)
x
outlet
where the last term represents a change in momentum flux
(either as a result of the lateral influx of momentum or due to
a storage change). Finally, the rate of change of momentum
storage is given by:
d
Q
V
dV
=
dx
dt CV
t
(10.7.12)
# y
&
Q (Q 2 A)
+
+ gA % S 0 + S f + Se ( qvx = 0
t
x
$ x
'
A
Q
+
q =0 = 0
t =0 x
which simply states that Q = constant from one cross-section
to the next (i.e. Q1 = Q2 or V1A1 = V2A2), which is the
commonly applied steady-state mass balance. Under the same
assumptions, the momentum equation simplifies as well:
Q
t
# y
&
(Q 2 A)
+
+ gA % S 0 + S f + Se ( q =0vx = 0
x
$ x
'
=0
(10.7.13)
353
dy
dh
S0 =
dx
dx
(10.7.14)
Fr =
gD
D=
A
;
B
B top width =
dA
dy
(10.7.17)
dh d V 2
+
= S f Se
dx dx 2g
d V
2 dy
=
F
r
dx 2g
dx
2
(10.7.18)
(10.7.15)
dy d V 2
+
= S0 S f
dx dx 2g
dy S 0 S f
=
dx
1 Fr2
(10.7.16)
Q
t
=0
y
(Q 2 A)
+
+ gA
x
x
=0
=0
S 0 S f Se q =0vx = 0
which yields:
S 0 = S f + Se
(10.7.19)
S0 = S f
(10.7.20)
The special case of uniform flow is one of particular
interest and worthy of additional discussion. From dimensional
analysis one can show that the bed shear stress associated
with friction is given by:
V =
V2
0 = Cf
= gRS 0
2g
(10.7.21)
V =
2g
RS 0 = C RS 0
Cf
(10.7.22)
1 23 12
R S0
n
(10.7.24)
Q=
1
AR 2 3S 01 2
n
(10.7.25)
C =
1 16
R
n
(10.7.23)
355
A
(1.5 m 2 )
R= =
= 0.375 m
P (2(0.5 m) + 3 m)
1
Q=
(1.5 m 2 )(0.375 m)2 3(0.001)1 2 = 1.2 m 3 s 1
(0.02)
2 3
! (3 m)y $
((3 m)y) #
&
" 2y + (3 m) %
= (0.1 m 3 s 1 )(0.02)(0.001)1 2
= 0.0632
E XAMPLE 10.7.2
P = B + 2y B, if B << y
! By $
1
12
Q = (By) #
& S0
n
" 2y + B %
which shows that even for a simple cross-sectional
geometry, the nonlinearity of the Manning equation
requires an iterative solution for the depth y. Specifically,
2 3
" By %
1
1
Q (By) $ ' S 01 2 = By 5/3S 01 2
n
n
#B &
y = (QnS 01 2B 1 )3/5
For the conditions shown above this would yield a depth
estimate of 0.099 m, which is a close approximation to
the real solution. This approximation is only valid under
the wide-channel assumption mentioned above.
356
The Saint-Venant equations are sometimes referred to as
the dynamic wave equation because they fully describe the
dynamics of 1D unsteady flow in a channel. In this regard
they can be used to model any number of phenomena
including floods, tides, nonuniform flow, uniform flow, etc.
The price for this generality is one of computational demand
since the Saint-Venant equations are expensive to solve.
The kinematic wave equation is another special case of
the Saint-Venant equations that is generally more easily
solved. By definition a kinematic wave (as opposed to a
dynamic wave) is one where the acceleration terms, pressure
term, and lateral influx of momentum are all negligible, i.e.
The result of the above simplifications is essentially a
reduction of the two equation governing system to a single
governing equation. For uniform flow the momentum equation
can also be expressed in the general form relating A and Q as:
A = aQ b Q = cAd
A Q
+
= q(x,t)
t
x
(10.7.26)
S0 = S f
(10.7.27)
3 5
3
;b =
5
or
12
1 S0
c=
;
n P2 3
d=
5
3
(10.7.29)
A = aQ b
(10.7.28)
nP 2 3
a=
S0
Q (Q A)
y
+
gA
qvx 0
t
x
x
A dA Q
Q
=
= abQ b1
t dQ t
t
(10.7.30)
abQ b1
Q Q
+
= q(x,t)
t
x
(10.7.31)
Q = cAd
Q dQ A
A
=
= cdAd 1
x dA x
x
(10.7.32)
A
A
+ cdAd 1
= q(x,t)
t
x
(10.7.33)
358
S ECTION 8
F IGURE 10.39 Stream channel reach storage conceptualization used in the Muskingum hydrologic streamflow routing
method (adapted from Mays, 2005).
Hydrologic routing starts with a lumped mass balance
equation for a specified reach of channel (Figure 10.39):
dS
= I(t) Q(t)
dt
(10.8.1)
(10.8.2)
359
S j +1 S j
t
I j +1 + I j
2
Q j +1 + Q j
2
(10.8.3)
(10.8.4)
Q j +1 = C 1I j +1 + C 2I j + C 3Q j
(10.8.5)
C1 =
t 2KX
2K(1 X ) + t
(10.8.6)
C2 =
t + 2KX
2K(1 X ) + t
(10.8.7)
C3 =
2K(1 X ) t
2K(1 X ) + t
(10.8.8)
where the C1, C2, and C3 coecients sum to 1.0, meaning that
Equation (10.8.5) is simply a weighted average of the righthand-side terms. Provided the inflow hydrograph (as a
function of time), the initial outflow, and the parameters K
and X, Equation (10.8.5) can be applied recursively to yield
the outflow hydrograph. Methods for determining the
parameters generally require historical inflow and outflow
hydrographs or estimation via other means (Mays, 2005).
Cunge (1969) provided a connection between the
Muskingum method and the kinematic wave method. Recall
that the kinematic wave model is a hydraulic routing method
that provides a distributed estimate of discharge (i.e. at
discretized locations along the channel). One can express the
discretized solution of the kinematic wave model as:
j+1
j
Qi+1
= C 1Qij+1 +C 2Qij +C 3Qi+1
(10.8.9)
S ECTION 9
Measurement of Streamflow
In the preceding sections it was assumed in various
locations that hydrograph data was available, i.e. at a basin
outlet. Here we briefly outline the primary ways those
measurements are made.
As with many hydrologic variables, streamflow is
generally not measured directly, but estimated indirectly from
measurements that are more straightforward to make. In a
stream, the measurements which are easiest to make include
water depth and velocity. In general, streamflow is estimated
via a rating curve, which is a derived or known from a
relationship between streamflow depth and discharge (flow).
Examples where a rating curve is known theoretically (or
empirically) include flow over or through certain constructed
structures (often called weirs or flumes; Mays, 2005; WMO,
2010). Figure 10.40 shows an example of a weir which is
constructed into a channel reach. The primary point of such
structures is to force the flow to go through a known flow
state transition, i.e. critical flow, which is a transition
between supercritical and subcritical flow conditions (see
Mays, 2005 for a thorough discussion of critical flow and
depth). This is usually accomplished via a constriction in the
flow geometry. When the flow is critical and the structure
geometry is known, there are empirical relationships between
flow and the critical depth or measured head. In such cases,
stalled (bottom panel) v-notch weir in a stream channel for estimating streamflow (from WMO, 2010).
Another mechanism for deriving rating curves and
streamflow measurements is via the manual measurement of
mean velocities in a channel cross-section. Velocity in a stream
can be measured using vertical axis mechanical current meters
(WMO, 2010; Figure 10.41), which are devices placed
perpendicular to the flow. The flow velocity turns a propellor
or anemometer (attached to a wading rod) which can then be
used to back out the velocity at the measurement point. It
should be noted that flow in a stream is generally turbulent,
which means that the velocity field can vary significantly (and
somewhat randomly) in space and time. To get an accurate
measurement of the mean velocity, measurements should be
taken over a long enough period to average over the turbulent
eddies. While flow is often written as: Q = VA, the V is the
average velocity across the entire stream cross-section, so
Q = V dA
A
(10.9.1)
Q = Viyi wi
(10.9.2)
i =1
streamflow from stream cross-section area and multiple velocity measurements (adapted from Mays, 2005).
362
363
S ECTION 10
Conceptual Questions
S ECTION 11
Sample Problems
Problem 10.1. You are hired to assess the flash flood
potential of a given basin, where flash floods are generally
associated with short duration (high intensity) storms that
generate infiltration excess runo. The basin is composed of a
homogeneous silt loam soil. To comply with regulations for
the region, dierent downstream infrastructure must be
designed for various return-period storm events, where the
return-period is associated with the probability (and therefore
magnitude) of the storm event. The 1-hour duration design
storm events for the 10-year and 25-year return periods in this
region are:
STORM RETURN PERIOD
10-year
25-year
PRECIP. INTENSITY
2 cm/hour
7 cm/hour
TIME
(HOUR)
FLOW (m3/s/cm)
0.5
1.0
10
1.5
15
2.0
2.5
3.0
TIME
(HOUR)
DISCHARGE
(m3/s)
3.11
3.45
6.51
16.36
18.25
12.28
8.29
5.72
4.53
10
3.11
1-HR UH
TIME
(HOUR)
UH FLOW
(m3/s/cm)
b) Find and plot the 2-hour unit hydrograph for the basin.
10
367
Chapter 11
A Simple
Watershed
Model
S ECTION 1
Learning Objectives
By the time you finish this chapter you should be able to:
1. Understand the key components (i.e. inputs, governing
equations, and outputs) of a distributed watershed model
2. Describe the basic idea behind the TOPMODEL approach
3. Understand how analytical mass/energy balance equations
can be discretized and solved numerically
4. Setup and run watershed simulations using the provided
MOD-WET numerical watershed model code
5. Reconcile model simulation output (i.e. watershed
response) with your hydrologic understanding developed
in previous chapters.
369
S ECTION 2
! a $
i = ln # i &
" tan Si %
(11.2.1)
K s = K 0 exp(z i / m)
(11.2.2)
T(z i ) =
K(z)dz K m exp(z
zi
/ m) = T0 exp(z i / m)
(11.2.3)
! a
$
i
i = ln #
&
T
tan
S
" 0
i%
(11.2.4)
372
dTs
= CT [Rn LE H ] C d (Ts Td )
dt
(11.2.5)
dSWE
= P E M
dt
(11.2.6)
C snow
dTsnow
= Rn LE H wLf M
dt
(11.2.7)
where Csnow is the snow-layer heat capacity (J m-2 K-1) and the
right-hand-side terms are the net radiation, latent and
sensible heat fluxes, and latent heating due to melt. Advected
energy by rain could also be included.
For simplicity, the soil surface energy balance is only
solved for prognostically when snow disappears. In doing so,
the soil surface temperature is set equal to the snow
temperature (i.e. 0C) just before the snow disappears. This is
obviously a simplification as it ignores soil dynamics under
373
dSrz
= f E qdrain ;
dt
(11.2.8)
dSuz
= qdrain qv ;
dt
0 Suz SD
(11.2.9)
qv = K 0 exp(SD / m)
(11.2.10)
qb = T0 exp(SD / m)tan S
(11.2.11)
d SD
dt
= Qb Qv
(11.2.12)
SD = SD + m "# i $%
(11.2.13)
B = Acc
(11.2.14)
!
tan S
n = n0 #
# tan S
"
1/3
$
& ; S mean slope
&
%
K=
dx
ck
(11.2.17)
(11.2.15)
Aside from the spatial discretization of the river network
(dx) and temporal discretization (dt), the C coecients
depend on the dynamic parameters K and X. The parameter
K is given by:
(11.2.16)
ck =
dQ
dA
(11.2.18)
Q=
1
1
AR 2 3S 01 2 = A5/3P 2 3S 01 2
n
n
(11.2.19)
ck =
dQ
5 2/3 2 3 1 2
=
A P S0
dA 3n
3/5
1/2
5 " S0 %
= $ 2/3 ' Q 2/5
3 $#nP '&
(11.2.20)
(11.2.21)
376
P = B + 2y B (for B >> y)
(11.2.22)
"
%
Q
X = 0.5 $$1
''
Bc
tan(S)dx
#
&
k
(11.2.23)
For a given time step in the model, Equations (11.2.14)(11.2.23) can be applied to each reach of the stream network
from upstream nodes all the way downstream to the basin
outlet. In some cases the flow velocity may be larger enough
that it can cover the distance between two pixels in less than
the specified time step. A dynamic time step can be used to
reduce the routing time step for such cases. The streamflow at
the basin outlet is the outlet hydrograph.
It should be reiterated that the version of the MODWET model described herein does not contain explicit
representation of vegetation. As such, it is best suited to nonvegetated basins. We would expect vegetation to modify the
processes in a variety of ways including: rainfall/snowfall
canopy interception, additional control on evapotranspiration
via stomatal resistance, attenuation of solar radiation reaching
377
S ECTION 3
Numerical Implementation
of the MOD-WET Model
The MOD-WET model defined in Section 2 needs to be
implemented numerically since analytical solutions are not
generally available. Numerical solutions are obtained by
discretizing the equations in space and time. Here the spatial
discretization is chosen to coincide with the underlying DEM
or other spatially-distributed fields that are necessary inputs
for the model. So, for example, a typical spatial discretization
may be on the order of 30 m x 30 m using readily available
DEMs. Alternatively, the DEMs may be coarsened to a lower
resolution to reduce computational expense. The time
discretization is done via a time step dt. The choice of the
time step is generally governed by the dynamics of the system.
To resolve the surface energy balance, a time step of less than
an hour may be necessary, while for more slowly varying
dynamics a time step as large as a day may be acceptable.
For numerical implementation, the above equations are
posed by representing all storage (or deficit) states in terms of
equivalent depth of water and the energy state at the surface
is represented by surface temperature (either snow or soil). In
the discretization used below, the indices i and t are used to
represent spatial and temporal dependence respectively. In the
solution of all dierential equations, an explicit forward finite
dierence scheme is used for computational eciency. This is
Ts (i,t + 1) = Ts (i,t) +
(11.3.1)
For snow-covered pixels the discretized single-layer snow
mass balance and snow surface energy balances are given by
(based on Equation 11.2.6):
(11.3.2)
dt
[R (i,t) LE(i,t) H(i,t) wLf M(i,t)]
C snow n
(11.3.3)
378
(11.3.4)
(11.3.5)
(11.3.6)
(11.3.7)
(11.3.8)
(11.3.9)
(11.3.10)
(11.3.11)
as given by:
(11.3.12)
1
Qv (t) =
N
q (i,t)
i=1
(11.3.13)
dt 1
Qb (t) =
Lx N
q (i,t)
i=1
(11.3.14)
SD !(i,t + 1) = SD (t + 1) + m #$ (i)%&
(11.3.15)
where
!
$
a(i)
(i) = ln #
&
T
(i)tan
S(i)
" 0
%
(11.3.16)
!
$
$
a(i)
1 N !
a(i)
= ln #
& = ln #
&dA(i) (11.3.17)
T
(i)tan
S(i)
A
T
(i)tan
S(i)
" 0
%
" 0
%
i=1
"
$ q ! (i,t) + SD !(i,t + 1) , SD !(i,t + 1) < 0
q se (i,t) = # se
q se! (i,t),
otherwise
$%
(11.3.18)
(11.3.20)
Initialization of the model can be done by specifying an
initial condition of the basin-averaged saturation deficit and
then applying Equation (11.3.15). Other initial states that
need to be specified include the rootzone and unsaturated
zone storage values. A summary of the states and parameters
(in addition to the DEM and meteorological inputs) that need
to be specified for this numerical implementation of the model
are shown below in Table 11.1. Many other parmeters must
also be specified. In the default implementation of the model,
the soil parameters (roughness, hydraulic, thermal, and
VARIABLE
SD(t 0 )
Suz (i,t 0 )
Srz (i,t 0 )
Srz max
Srz min
T0
the fluxes that occurred over the time step. Finally, the raw
data is averaged or saved at the specified temporal resolution.
Outputs are stored in structured arrays based on variable
type, i.e.: state_maps, flux_maps, and time_series, where
each contains the key variables and the units of each.
The MOD-WET model is constructed in an attempt to
run reasonably eciently, but also in a way that makes for
easy student learning. The forward finite dierence
formulations for the energy mass/budgets are simple to
understand and allows for vectorization of the calculations.
This simply means that the change in any model state (e.g.
temperature) over the course of one time step can be
computed simultaneously for all pixels across the domain.
Alternatively, one could loop over pixels in the domain, but
loops tend to be very slow relative to vectorized calculations
and hence this saves a significant amount of time. The
primary drawback of these simple finite dierence schemes is
the need for a small time step (e.g. 15 min.), which tends to
be most necessary for the surface energy budget calculations.
More sophisticated finite dierence schemes (implicit/
iterative) may be more stable and therefore allow a larger
time step, but are generally less straightforward to vectorize.
The primary computational (CPU) expense is proportional to
the number of pixels being simulated (i.e. size of domain and
resolution of DEM) and the length of the simulation (number
of time steps required). Hence, short-duration simulations
with a small domain at coarse resolution will take the least
amount of CPU time. Long-duration simulations over a large
domain at high resolution will take the most CPU time. The
382
383
S ECTION 4
Example MOD-WET
Model Applications
The model developed in the last two sections is done
generally so that it can be applied to basins of varying sizes
and characteristics. The key inputs to the model include a
DEM, soil parameters, meteorological data, etc. Here a simple
toy basin is used for qualitative demonstration of the model.
Figures 11.3 and 11.4 shows the DEM for the toy basin.
Geographically, its location (i.e. easting/northing coordinates)
F IGURE 11.4 Surface plot of the toy basin DEM showing the
varying slopes along the valley hillslope cross-section.
The basin consists of a simple geometry with a
symmetric valley that has a stream that runs from East to
West with a specified bed slope. The hillslopes of the valley
have dierent slopes: a shallower slope in the lower part of the
valley (i.e. a flood-plain) with a steeper slope up to the ridge
lines. The basin outlet has an elevation around 3000 m with
the highest ridge line ranging up to 3050 m. The toy basin
configuration is chosen in order to highlight some of the key
points made earlier in the book including variability in
radiation, snowmelt, runo contributing areas, etc. Figure
11.5 shows the computed slope and aspect maps for the basin.
384
F IGURE 11.6 Map of the topographic index for the toy basin.
385
F IGURE 11.8 Saturation deficit map (in meters) over the toy
basin domain after 300 days of drainage.
A second experiment which is illustrated in this section
is a full-year simulation over the basin using realistic
meteorological forcing data. Figure 11.9 shows the dailyaveraged precipitation, incoming shortwave radiation, and air
temperature over a water year (WY; i.e. October 1st September 30th) corresponding to a meteorological station at
an elevation of 3300 m. It should be reiterated that the model
is run at a 15-minute time step to provide a stable solution for
the surface energy balance. The data is representative of the
Southern Sierra Nevada climate, showing the majority of
variables for full year simulation: precipitation (top panel), incoming shortwave radiation (middle panel), air temperature
(bottom panel). These forcing variables from a gage at a specified elevation and distributed spatially using topography.
387
Figure 11.12 shows the daily-averaged surface energy
fluxes. Net radiation shows a decreasing trend early in the
simulation in response to decreasing incident shortwave
radiation and air temperature (which reduces incoming
longwave radiation) as shown in Figure 11.9. It is relatively
low during the snow season and then increases during the
summer. The daily-averaged sensible heat flux is generally
positive during the summer when the surface is warmer than
the overlying air and is generally negative when the surface is
389
Figure 11.13 shows the daily-averaged runo terms
from the simulation. Saturation excess runo is shown to be
the largest contributor to runo for this basin and model
inputs. In fact, the infiltration excess runo is zero throughout
the simulation. This is because the specified saturated
hydraulic conductivity for the basin is larger than the
precipitation (or melt) rate experienced throughout the
simulation. Most of the saturation excess runo is in response
to intermittent melt events during the snow accumulation
season as well as to the sustained melt during Days 200-275 of
the Water Year.
saturation and infiltration excess runoff (top panel) and baseflow (bottom panel).
F IGURE 11.14 SWE map (in meters) over the toy basin domain on Day 260 of the Water Year.
390
over the toy basin domain on Day 260 of the Water Year.
the toy basin for each day of the Water Year (DOWY).
391
M OVIE 11.2 Animation of daily-average rootzone soil moisture (in meters) over the toy basin for each day of the Water
Year (DOWY).
F IGURE 11.16 Saturation deficit map (in meters) over the toy
basin domain on Day 260 of the Water Year.
392
meters) over the toy basin for each day of the Water Year
(DOWY).
M OVIE 11.4 Animation of daily-average saturation excess runoff (in m3/h) over the toy basin for each day of water year
(DOWY).
The runo from Tokopah is primarily snowmelt driven.
For the precipitation and/or snowmelt rates and soil
parameters used, infiltration excess runo is zero throughout
the simulation. The saturation excess runo is shown in Movie
11.9. The results show that the saturation excess runo is zero
over most of the domain and is intermittently non-zero early
in the WY around the main channel during intermittent
rainfall/snowmelt events. The saturation excess runo is most
persistent in this same region during snowmelt. Note that
much of this predicted runo is exfiltration resulting from
upstream snowmelt infiltration that raises the water table.
The fact that most of the runo is the result of saturation
S ECTION 5
Sensitivity of Hydrograph
Response to Watershed
Characteristics
This section uses the MOD-WET model to illustrate how
hydrologic response (in terms of the runo hydrograph at the
outlet) is impacted by various basin characteristics. As a
starting point, the toy basin used in Section 4 is chosen as a
baseline (nominal) watershed (Figure 11.20; upper left
panel). The basin is a simple rectangular shape with the main
stream channel moving from east to west with the outlet at
the far left (west) center of the domain. The nominal domain
is 15 x 10 pixels on a DEM of 30 m resolution so that the
nominal basin area is 0.135 km2. The flow patterns are
therefore relatively straightforward with hillslope flow
primarily in the north-south direction and main stream
channel flow in the east-to-west direction. To assess the
impact of shape and size, three other basins are synthesized.
The longer basin (Figure 11.20; upper right panel) is the
same width as the nominal, but is four times longer (basin
area of 0.54 km2). The wider basin (Figure 11.20; lower left
panel) is the same length as the nominal, but is four times
wider (basin area of 0.54 km2). The larger basin (Figure
11.20; lower right panel) is both four times longer and wider
than the nominal (basin area of 2.16 km2).
let) of simple synthetic rectangular watersheds: nominal (upper left), longer (upper right), wider (lower left), and
larger (lower right). The outlet of each basin is at the far left/
center of each rectangular domain. White areas are outside of
each basin and all basins are plotted using the same horizontal
scale.
For simplicity, results from a Hortonian or infiltration
excess runo (Chapter 10, Section 3) event (occurring
uniformly over the basin) are used for illustration. Such events
provide a much more straightforward way to gain insight into
the sensitivity of the runo hydrograph to basin
398
401
S ECTION 6
MOD-WET Codes
Relevant new functions based on concepts introduced in this chapter include:
Distributed model driver function:
MOD_WET_model_driver.m
Specification of model input parameters:
MOD_WET_model_static_and_control_parameters.m
Initialization of model states and distributed parameters:
initialize_model.m
Topographic index:
topo_index.m
Topographic-soil index:
topo_soil_index.m
TOPMODEL wrapper for evolving subsurface states and
fluxes:
TOPMODEL.m
Construction of flow network:
flow_network.m
402
S ECTION 7
Conceptual Questions
1. Provide a basic description of the topographic index and
what it represents in the TOPMODEL framework.
2. In the TOPMODEL framework, what is typically assumed
about the saturated soil hydraulic conductivity as a
function of depth?
3. In the TOPMODEL framework, what must the value of
saturation deficit (SD) be at a given pixel for it to generate
saturation excess runo?
4. In the TOPMODEL framework, what variables control the
baseflow from a given pixel?
5. Describe the meaning of the wide channel approximation
and why it is typically invoked?
6. What static channel parameters and dynamic model
variables are needed in the Muskingum-Cunge routing
equations?
403
S ECTION 8
Sample Problems
Chapter 12
MATLAB
Basics
It is generally assumed that the reader has had a course
in basic MATLAB or another programming language. The idea
of this section is not to provide a full primer or reference, but
highlight some key points to help in using MOD-WET
functions and provide pointers to more general references.
First, it should be noted that MATLAB is capable of much
more than covered here, has an extremely useful built-in help
system, and Mathworks provides their own rather extensive
primer for beginners located at: http://www.mathworks.com/
help/pdf_doc/matlab/getstart.pdf. You should consult these
resources and keep them handy for reference. There are also
many textbooks designed to provide more thorough
introductions to MATLAB (e.g. Chapman, 2007; Pratap, 2009).
The best way to find a solution to a problem you have is to
find an example of how such things are done via the help
command, using the primer, or other online resources. You
should also be aware that none of the applications used in this
book need to be solved using MATLAB, it is simply chosen as a
numerical tool since it is readily available and is relatively
easy to learn.
In its simplest form, MATLAB can be used as a calculator
to perform basic arithmetic operations on values stored in
variables, vectors, or multi-dimensional arrays. As a first step,
you need to be comfortable with how variables are stored in
arrays, array operations, indexing in arrays, etc. (see the
Language Fundamentals chapter in the MATLAB primer).
The idea of specifying a search path is particularly important when using the MOD-WET functions. Paths can be
manually added to MATLAB, but the approach described below automates this process with a simple line that you can
add to your own code. MOD-WET has several subdirectories.
Assume that you have the following directory structure:
/Users/username/CEE150/MOD_WET/chapter1
/Users/username/CEE150/MOD_WET/chapter2
/Users/username/CEE150/MOD_WET/chapter3
406
/Users/username/CEE150/MOD_WET/chapter5
/Users/username/CEE150/MOD_WET/chapter6
/Users/username/CEE150/MOD_WET/chapter7
/Users/username/CEE150/MOD_WET/chapter8
/Users/username/CEE150/MOD_WET/chapter11
MATLAB provides several colorbar options, with jet
being the default. Other colorbars, or ones you design yourself
may be more useful in other circumstances. For example, the
hsv colormap is useful for plotting maps of aspect or wind
direction (circular/angular directions) because angles of both
0 and 360 degrees should be represented by the same color
since they are actually the same angle. High and low values
are represented with the same color in this colormap which
can be invoked via: colormap(hsv).
Finally, invariably you will need to include graphical
figures into homework assignments or lab reports being
written in a word processing program. There are multiple
ways to save figures. Saving can be done from within a script
or function using the saveas command, e.g.
408
saveas(gcf,'figure_name.png','png')
409
Chapter 13
References
412
414