Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

6 Water Budget Record from Variable Infiltration Capacity(VIC) Model
Algorithm Theoretical Basis Document, Version 1.2

June 19, 2009
Huilin Gao1, Qiuhong Tang1, Xiaogang Shi1, Chunmei Zhu1, Ted Bohn 1, Fengge Su 1,
Justin Sheffield2, Ming Pan2, Dennis Lettenmaier1, Eric F. Wood2

1

Department of Civil and Environmental Engineering
University of Washington
Seattle, WA 98195
hgao@hydro.washington.edu
2

Department of Civil and Environmental Engineering
Princeton University
Princeton, NJ 08544

Abstract
This document describes the algorithms within the latest version of the variableinfiltration capacity (VIC) model. As a semi-distributed macroscale hydrological model,
VIC balances both the water and surface energy within the grid cell; and its sub-gridvariations are captured statistically. Distinguishing characteristics of the VIC model
include: subgrid variability in land surface vegetation classes; subgrid variability in thesoil moisture storage capacity; drainage from the lower soil moisture zone (base flow) as
a nonlinear recession; and the inclusion of topography that allows for orographic
precipitation and temperature lapse rates resulting in more realistic hydrology inmountainous regions. VIC uses a separate routing model based on a linear transfer
function to simulate the streamflow. Adaptations to the routing model are implemented in
VIC to allow representation of water management effects including reservoir operationand irrigation diversions and return flows. Since its existence, VIC has been wellcalibrated and validated in a number of large river basins over the continental US and theglobe. Applications using the VIC model cover a variety of research areas.

120

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Given the numerous improvements and updates of the VIC model through its nearlytwenty years of existence, this document serves as a general guideline for helping usersof the long term dataset to understand the fundamental VIC algorithms up to date.
Section 6.1 serves as an introduction, Section 6.2 gives a historical overview of the VIC
model development, and Section 6.3 explains the classic algorithms of the VIC model for
calculating the state variables, surface fluxes, and streamflow, as well as the newly
implemented algorithms for taking into account the water management. Section 6.4
describes the model forcings and model parameterizations, and Section 6.5 is about the
VIC calibration. Finally, Section 6.6 summarizes VIC validation and applications.

6.1

Introduction

The variable infiltration capacity (VIC) model (Liang et al., 1994, 1996), with a variety
of updates (Cherkauer et al, 2003; Bowling et al., 2004; Bowling and Lettenmaier, 2009),
has been extensively used in studies on topics ranging from water resources management
to land-atmosphere interactions and climate change. Throughout its existence, VIC has
played multiple roles, as both a hydrologic model and land surface scheme when coupled
to general circulation models. As a semi-distributed macroscale hydrological model, VIC
balances both the water and surface energy budgets within the grid cell; and its sub-gridvariations are captured statistically. Distinguishing characteristics of the VIC model
include: subgrid variability in land surface vegetation classes; subgrid variability in thesoil moisture storage capacity; drainage from the lower soil moisture zone (base flow) as
a nonlinear recession; inclusion of topography that allows for orographic precipitationand temperature lapse rates resulting in more realistic hydrology in mountainous regions.
To simulate streamflow, VIC results are typically post-processed with a separate routing
model (Lohmann, et al., 1996; 1998a; b) based on a linear transfer function to simulate
the streamflow. VIC has been adapted to allow representation of water management
effects (Haddeland et al, 2006a; b ; 2007) including reservoir operation and irrigation
diversions and return flows.
VIC has been well calibrated and applied in a number of large river basins over thecontinental US and the globe (Abdulla et al. 1996; Bowling et al. 2000; Lohmann et al.
1998b; Nijssen et al. 1997, 2001a; Shi et al., 2008; Su et al., 2005, 2006; Wood et al.
1997; Zhu and Lettenmaier, 2007). VIC has participated in the WCRP Intercomparison ofLand Surface Parameterization Schemes (PILPS) project and the North American Land
Data Assimilation System (NLDAS), where it has performed well relative to otherschemes and to available observations (Bowling et al, 2003a, b; Lohmann et al., 2004;
Nijssen et al. 2003; Wood et al., 1998). It has also been evaluated using soil moistureobservations in the U.S. (Maurer et al, 2002) and global snow cover extent data by(Nijssen et al, 2001b).

121

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Driven by high-quality meteorological forcings, VIC had been used to provide a longterm data record of land surface fluxes and states for the conterminous United States(1950-2000) (Maurer et al., 2002) and Mexico (1925-2004) (Zhu and Lettenmaier, 2007).
Applications using such a data record have covered many areas, such as: simulating
ensembles of streamflow and hydrologic variables for forecast purpose ( Hamlet andLettenmaier, 1999; Wood et al., 2002, 2005; Wood and Lettenmaier, 2006 );reconstructing and analyzing drought events ( Andreadis and Lettenmaier, 2006a;
Sheffield et al., 2004a; Sheffield and Wood, 2007; Wang et al., 2009); studying the North
American monsoon teleconnections (Zhu and Lettenmaier, 2005; Zhu et al., 2007, 2009);
drought prediction (Luo and Wood, 2007); conducting hydrologic studies over the Panarctic region (Bohn et al., 2007; Bowling et al., 2003c; Lettenmaier and Su, 2009; Slater
et al., 2007; Su et al., 2005, 2006); water management (Adam et al., 2007; Haddeland etal, 2006a, b, 2007); and many others (Section 6.6).
Given the numerous improvements and updates of the VIC model through its nearlytwenty years of existence, this document serves as a general guideline for helping usersof the long term dataset to understand the fundamental VIC algorithms. Section 6.2 gives
a historical overview of the VIC model development; Section 6.3 explains the classic
algorithms of VIC model for calculating the state variables, surface fluxes, andstreamflow, as well as the newly implemented algorithms for the water management.
Section 6.4 describes the model forcings and model parameterizations. Section 6.5 isabout the VIC calibration. And section 6.6 summarizes VIC validation and applications.

6.2

Historical Overview of the VIC Model

The VIC model was developed for incorporation in GCMs, aiming to improve therepresentation of horizontal resolution and subgrid heterogeneity in a simple way.
Employing the infiltration and surface runoff scheme in Xianjiang model (Zhao, 1980),
VIC was first described as a single soil layer model by Wood et al. (1992) and
implemented in the GFDL and Max-Planck-Institute (MPI) GCMs (Stamm et al. 1994).
The single soil layer model requires three parameters: an infiltration parameter, anevaporation parameter, and a base flow recession coefficient. In 1994, Liang et al. (1994)
generalized the two-layer VIC model (VIC-2L) to include the multiple soil layers andspatially varying vegetation and evaporation within a grid cell. In VIC-2L, infiltration,
drainage from the upper soil layer into the lower soil layer, surface and subsurface runoff
are calculated for each vegetation cover tile (in addition to the statistical parameterizationof heterogeneity of infiltration and runoff generation within a vegetation cover tile
present in the original VIC model). Therefore, the subgrid-scale heterogeneity isrepresented in soil moisture storage, evaporation, and runoff production. As a semidistributed land surface model, VIC calculates the sensible and latent heat fluxes
according to physical formulations, but it uses conceptual schemes to represent thesurface runoff and base flow. In 1996, Liang et al. (1996) found that the VIC-2L tends tounderestimate the evaporation due to the low soil moisture in its upper soil layer, and the122Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

main cause of this error is the lack of a mechanism for moving moisture from the lower
to the upper soil layer. VIC-2L was then modified to allow diffusion of moisture between
soil layers, and to have an additional 10cm thin soil layer on top of the previous upper
soil layer. In this way the three-layer VIC model (VIC-3L) was generated, and the VIC3L framework has been used ever since. The model currently allows for more than threesoil layers if desired.

A number of modifications to VIC have been made to improve the model such that it candeal with complicated hydrological processes. Since the VIC model does not represent
the geometry of the sub-grid variations, a separate routing model has been developed tosimulate the streamflow (Lohmann., et al., 1996, 1998a, 1998b). To represent the cold
land processes, the VIC model was upgraded to include a two-layer energy balance snow
model (Andreadis et al., 2009; Wigmosta et al., 1994; Storck et al., 1998), frozen soil andpermafrost algorithm (Cherkauer et al., 1999, 2003; Cherkauer and Lettenmaier, 2003),
and blowing snow algorithm (Bowling et al., 2004). To improve the simulations ofelevation-dependent components within a grid cell, elevation bands representing
topography were introduced (Nijssen et al., 2001b). With the evapotranspiration
algorithm, canopy responses to wind profile and surface radiation budget have been
incorporated (Wigmosta et al., 1994), and the leaf area index (LAI) and the vegetation
fraction were allowed to vary at each time step (Liang et al., 1996). The effects lakes oflake and wetlands on moisture storage and evaporation, which are particularly important
for runoff at high latitude, have been included ( Bowling et al, 2003c; Bowling andLettenmaier, 2009; Cherkauer et al., 2003 ). To simulate water management impacts, areservoir module has been implemented to the routing model and a sprinkle irrigation
scheme has been added to the soil moisture simulation (Haddeland et al., 2006a, 2006b,
2007).

Besides the above improvements to the water budget and energy balance processes in the
VIC model, efforts have been made to provide better meteorological forcings through thedata preprocessor. Using algorithms by Kimball et al. (1997), Thornton and Running(1999), and Bras (1990), a full suite of hydrologic variables is constructed from limited
observed driving data (precipitation, maximum and minimum air temperature, and windspeed) (Nijssen et al., 2001b).

6.3

VIC Model Description

6.3.1 Overview of VIC Model Processes

The overall VIC model framework has been described in detail in literature (Liang et al.
1994; Liang et al., 1996; Nijssen et al., 1997). The key characteristics of the grid-based123Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

VIC are the representation of vegetation heterogeneity, multiple soil layers with variableinfiltration, and non-linear base flow.
Figure 6.1 shows the schematic of the VIC model with a mosaic representation ofvegetation coverage and three soil layers. The surface of each grid cell is described by
N+1 land cover tiles, where n = 1, 2, … , N represents N different tiles of vegetation, andn = N+1 represents bare soil. For each vegetation tile, the vegetation characteristics, suchas LAI, albedo, minimum stomatal resistance, architectural resistance, roughness length,
relative fraction of roots in each soil layer, and displacement length (in the case of LAI)
are assigned. Evapotranspiration is calculated according to the Penman-Monteithequation, in which the evapotranspiration is a function of net radiation and vapor pressure
deficit. Total actual evapotranspiration is the sum of canopy evaporation and transpirationfrom each vegetation tile and bare soil evaporation from the bare soil tile, weighted by
the coverage fraction for each surface cover class. Associated with each land cover typeare a single canopy layer, and multiple soil layers (three layers are used for description inthis ATBD). The canopy layer intercepts rainfall according to a Biosphere-atmospheretransfer scheme (BATS) parameterization (Dickinson et al., 1986) as a function of LAI.
The top two soil layers are designed to represent the dynamic response of soil to theinfiltrated rainfall, with diffusion allowed from the middle layer to the upper layer whenthe middle layer is wetter. The bottom soil layer receives moisture from the middle layer
through gravity drainage, which is regulated by a Brooks-Corey relationship (Brooks andCorey, 1988) for the unsaturated hydraulic conductivity. The bottom soil layer
characterizes seasonal soil moisture behavior and it only responses to short-term rainfall
when the top soil layers are saturated. The runoff from the bottom soil layer is according
to the drainage described by the Arno model (Franchini and Pacciani, 1991). Moisture
can also be transported upward from the roots through evapotranspiration. Although
vegetation subgrid-scale variability is a critical feature for the VIC model, the soil
characteristics (such as soil texture, hydraulic conductivity, etc.) are held constant foreach grid cell. In the model, soil moisture distribution, infiltration, drainage between soil
layers, surface runoff, and subsurface runoff are all calculated for each land cover tile ateach time step. Then for each grid cell, the total heat fluxes (latent heat, sensible heat, andground heat), effective surface temperature, and the total surface and subsurface runoff
are obtained by summing over all the land cover tiles weighted by fractional coverage.
The VIC model can be run in either a water balance mode or a water-and-energy balance
mode. The water balance mode does not solve the surface energy balance. Instead, it
assumes that the soil surface temperature is equal to the air temperature for the current
time step. By eliminating the ground heat flux solution and the iterative processes
required to close the surface energy balance, the water balance mode requires
significantly less computational time than other model modes. These simplifications,combined with the daily time step that is typical of water balance mode simulations,yields a substantial savings in computational time. The exceptions to this are that thesnow algorithm and the frozen soil algorithm, both of which run at a sub-daily time step,
and which solve the surface energy balance to determine the fluxes needed to drive
accumulation and ablation processes, or to solve the frozen soil penetration, respectively(Andreadis et al., 2009; Bowling et al., 2004; Cherkauer and Lettenmaier 1999; Storck et

124Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Modelal., 1998). The full water-and-energy balance mode not only solves the complete waterbalance but also minimizes the surface energy balance error. The surface energy balanceis closed through an iterative process which tries to find the surface temperature thatyields surface energy fluxes (sensible heat, ground heat, ground heat storage, outgoing
longwave and indirectly latent heat) so that balance the incoming solar and longwave
radiation fluxes. This mode requires more computational time than the water balancemode as well as requiring a sub-daily simulation time step. However, it is critical forstudies in which the land-atmosphere interactions are of interest (e.g., coupling withclimate models).
In the VIC model, each grid cell is modeled independently without horizontal water flow.
The grid-based VIC model simulates the time series of runoff only for each grid cell,which is non-uniformly distributed within the cell. Therefore, a stand-alone routing
model (Lohmann., et al., 1996, 1998a) is employed to transport grid cell surface runoff
and base flow to the outlet of that grid cell then into the river system. In the routing
model, water is never allowed to flow from the channel back into the grid cell. Once itreaches the channel, it is no longer part of the water budget. Figure 6.2 shows theschematic of the routing model. A linear transfer function model characterized by its
internal impulse response function is used to calculated the within-cell routing. Then by
assuming all runoff exits a cell in a single flow direction, a channel routing based on thelinearized Saint-Venant equation is used to simulate the discharge at the basin outlet.

Figure 6.1 Schematic of the VIC-3L model with mosaic representation of vegetation coverage.

125Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Figure 6.2 Schematic of VIC network routing models.

6.3.2 Water balance

The water balance in the VIC model follows the continuous equation for each time-step:S PERtwhere dS/dt, P, E, and R are the change of water storage, precipitation, evapotranspiration,and runoff, respectively. Within the time step, all units of above variables are mm. Overvegetated areas, the precipitation is the througfall (Pt). The water balance equation in the
canopy layer (interception) is:Wi P  E c  Pt
twhere Wi is canopy intercepted water (mm), Ec is evaporation from canopy layer (mm),and Pt is througfall (mm).

6.3.2.1 Evapotranspiration

The VIC model considers three types of evaporation: evaporation from the canopy layer(Ec, mm) of each vegetation tile, transpiration (Et, mm) from each of the vegetation tiles,
and evaporation from the bare soil (E1, mm) (Liang et al. 1994). Total evapotranspiration
over a grid cell is computed as the sum of the above components, weighted by therespective surface cover area fractions. The formulation of the total evapotranspiration is:126Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) ModelN

E   C n  ( Ec ,n  Et ,n ) C N 1  E1
n 1

Where Cn is the vegetation fractional coverage for the nth vegetation tile, CN+1 is the bareN 1soil fraction, and  C n  1.
n 1

6.3.2.1.1 Canopy evaporation
When there is intercepted water on the canopy, the canopy evaporates at the maximum
value. The maximum canopy evaporation ( Ec* , mm) from each vegetation tile iscalculated using the following formulation:

Ec*  (

Wi 2 / 3
rw
) Ep
Wim
rw  roWhere Wim is the maximum amount of water the canopy can intercept (mm), which is 0.2
times LAI (Dickinson, 1984); the power of 2/3 is as described by Deardorff (1978). The
architectural resistance, r0, is caused by the variation of the humidity gradient between
the canopy and the overlying air (s m -1). In the model, r0 is assigned for each land covertype according to the vegetation library. The aerodynamic resistance, rw, represents thetransfer of heat and water vapor from the evaporating surface into the air above thecanopy (s m -1). Ep is the potential evapotranspiration (mm) that is calculated from thePenman-Monteith equation (Shuttleworth, 1993) with the canopy resistance set to zero,
which is:( Rn  G )   a c p (es  ea ) / raEp where λv is the latent heat of vaporation (J kg-1), Rn is the net radiation (W m-2), G is thesoil heat flux (W m-2), (es - ea) represents the vapor pressure deficit of the air (Pa), ρ a isthe density of air at constant pressure (kg m -3), cp is the specific heat of the air (J kg-1 K-1),
Δ represents the slope of the saturation vapor pressure temperature relationship (Pa K -1),and γ is the psychrometric constant (66 Pa K -1). The Penman-Monteith equation as
formulated above includes all parameters that govern the energy exchange andcorresponding latent heat flux (evapotranspiration) from uniform expanses of vegetation.
The aerodynamic resistance (rw, s m -1) is described as follows after Monteith andUnsworth (1990):
1
rw 
C wu z127

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

where u z is the wind speed (m s-1) at level z, and Cw is the transfer coefficient for water
which is estimated taking into account the atmospheric stability. The algorithm forcalculating Cw is based on Louis (1979).
When the continuous rainfall rate is lower than the canopy evaporation, the interceptedwater is not sufficient for meeting the atmospheric demand within one time step. In sucha case, the canopy evaporation (Ec, mm) is
Ec  f  E c*
where f is the fraction of the time step for canopy evaporation to exhaust the interceptedwater, and it is given by: W  P  t f  min1, i *E c  t 
6.3.2.1.2 Vegetation transpiration

The vegetation transpiration (Et, mm) is estimated using (Blondin, 1991; Ducoudre et al.,1993):
Wrw
Et  (1  ( i ) 2 / 3 ) E p
Wimrw  ro  rc
Where rc is the canopy resistance (s m -1) given by:
r0 c g T g vpd g PAR g sm
rc 
LAI
where r0c is the minimum canopy resistance (s m-1) according to the vegetation library,and g T, g vpd, g PAR , and g sm are the temperature factor, vapor pressure deficit factor,
photosynthetically active radiation flux (PAR) factor, and soil moisture factor,
respectively. Details about the four limiting factors are available through Wigmosta et al.
(1994).

When canopy evaporation happens only for a fraction of the time step (f) (see Section
6.3.2.1.1), the transpiration during that time step then has two parts as described byrw
Wrw
Et  (1  f ) E p
 f  (1  ( i ) 2 / 3 ) E p
rw  ro  rc
Wim
rw  ro  rc
where the first term represents the part of the time step when there is transpiration but nocanopy evaporation, and the second term represents the part of the time step when there isboth evaporation from the canopy and transpiration.
The vegetation transpiration from a certain vegetation tile is the total contribution fromall three soil layers, weighted by the fractions of roots in each layer.
6.3.2.1.3 Bare soil evaporation

128

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

The bare soil evaporation only occurs on the top thin layer. When the surface soil issaturated, it evaporates at the potential evaporation rate. When the top soil layer is not
saturated, its evaporation rate (E1) is calculated using the Arno formulation by Franchiniand Pacciani (1991). The infiltration capacity (i) uses the spatially heterogeneous
structure described by the Xianjiang Model (Zhao et al., 1980), which is expressed asi  im (1  (1  A)1 / bi )withi m  (1  bi )   S  zwhere im is the maximum infiltration capacity (mm), A is the fraction of area for which
the infiltration capacity is less than i, b i is the infiltration shape parameter, θs is the soil
porosity, and z is the soil depth (m). All these variables are for the top thin soil layer.
The bare soil evaporation is described as
1
 AS

i0


E1  E p  dA  dA
1 / bi
0

i
(
1

(
1

A
)
)
AS m


with As denoting the fraction of the bare soil that is saturated, and i0 representing thecorresponding point infiltration capacity.

6.3.2.2 Soil Moisture and runoff

The VIC model uses the variable infiltration curve (Zhao et al., 1980) to account for the
spatial heterogeneity of runoff generation. It assumes that surface runoff from the upper
two soil layers is generated by those areas for which precipitation, when added to soil
moisture storage at the end of the previous time step, exceeds the storage capacity of thesoil. The formulation of subsurface runoff follows the Arno model conceptualization(Franchini and Pacciani, 1991; Todini, 1996). The soil moisture and runoff algorithms for
the VIC-3L is explained with details in Liang et al. (1996).
Similar to the total evapotranspiration, the total runoff Q is expressed as:N 1

Q   C n  (Qd ,n  Qb ,n )n 1where Qd,n (mm) and Qb,n (mm) are the direct runoff (surface runoff) and base flow(subsurface runoff) for the nth land cover tile, respectively.

The VIC model assumes there is no lateral flow in the top two soil layers; therefore themovement of moisture can be characterized by the one-dimensional Richard’s equation:K ( )
 ( D( ) ) 
t z
z
z
where θ is the volumetric soil moisture content, D(θ) is the soil water diffusivity (mm2 d1), K(θ) is the hydraulic conductivity (mm d-1), and z is soil depth (m). By including the129

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Modelatmospheric forcing, the integrated soil moisture for the top two soil layers can be
described as (Mahrt and Pan, 1984):
 i
(i=1,2)
 z i  I  E  K ( ) |  D ( ) |ZItz  Z I
where I is the infiltration rate (mm d-1), z1 and z2 are soil depth for layer 1 and layer 2,
respectively. The infiltration rate I is the difference between the precipitation (orthroughfall if there is vegetation coverage) and the direct runoff Qd.
For the lower soil layer, an empirical formulation derived from large scale catchmenthydrology is used in which the drainage and subsurface drainage are lumped together asbase flow (Qb). The soil moisture for the soil layer is described by the water balance
equation including diffusion between soil layers as:
 3 ( z 3  z 2 )  K ( )|  D ( ) |  E  Qb
 Z2
tz  Z 2
If it is bare soil, the evapotranspiration term E is zero because there is no evaporationfrom the lower soil layer. Otherwise, if the vegetation roots go through into the lower soil
layer, the evapotranspiration term E needs to be considered.
Since the top thin soil layer has a very small water holding capacity within each time
step, the direct runoff (surface runoff, Qd) is calculated for the entire upper layer (layer 1
and layer 2) as (Liang et al., 1996):
i0  P 1 bi
P  i0  im
P  z 2  ( S   2 )  z 2   S  (1  i ) ,Qd  
m
P  z  (   ),P  i0  i m2
S2where the infiltration capacity associated terms (i0, im, θs, and bi) are explained in Section
6.3.2.1.3.

The formulation of base flow (sub surface runoff, Qb), which used the Arno model
formulation, (Franchini and Pacciani, 1991), is expressed as: DS Dm W    3 , 0   3  WS  S
 S S
Qb  
 DS Dm    ( D  DS Dm )( 3  WS  S ) 2 ,   W m
3
S S
 WS  S 3WS
 S  WS  S
where Dm is the maximum subsurface flow (mm d-1), DS is a fraction of Dm ,, and WS is
the fraction of maximum soil moisture (soil porosity) θs. The base flow recession curve islinear below a threshold (WS θs) and nonlinear above the threshold. The first derivative at
the transition from the linear to nonlinear drainage is continuous.

130

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

6.3.3 Energy Balance (without snow or frozen soil)

The energy balance and its components were explained with details in Liang et al. (1994).
Using the same scheme as most other land surface models do, the energy balance
equation at the land surface for each land cover type within the grid cell is described as
Rn  H   w Le E  G
where Rn is the net radiation (W m-2), H is the sensible heat flux (W m-2),  w v E is thelatent heat flux (W m-2) (ρ w is the density of liquid water, kg m-3; Le is the latent heat ofvaporization, J kg-1), and G is the ground heat flux (W m-2). When the land surface is flatand homogeneous, the energy balance for a layer of air adjacent to the ground surface canbe expressed as
Rn  H  Le E  G  H
with ΔH representing the change of the energy storage rate in that air layer (W m-2). The
net radiation, and sensible and latent heat fluxes are from the top surface of the air layer,
while the ground heat flux is from the bottom of the air layer. The rate of energy storagechange, when is considered to be significant, is a c p (TS  TS ) z aH 2twhere ρ a is density of air (kg m-3), cp is specific heat of air at constant pressure (J kg-1 K1), TS and TS are the surface temperature of the bottom of the air layer at the end andbeginning of a time step (K), respectively, and za is the height of the air layer (m).
The net radiation which is the total of shortwave and longwave radiation is given byRn  (1   ) RS   ( R L  Ts4 )where α is the surface albedo for the land cover type, RS is the downward shortwave
radiation (W m-2), ε is the surface emissivity of the land cover type, RL is the downward
longwave radiation (W m-2), and σ is the Stefan-Boltzmann constant (5.67×10 -8 W m-2
K-4).

The water and energy balances are linked through the latent heat flux, with E=Ec+Et forvegetated land type and E=E1 for bare soil (see Section 6.3.2.1 for more details).
The sensible heat is given byac p
H(Ts  Ta )rhwhere rh is the aerodynamic resistance to heat flow (s m -1), Ts and Ta are surface
temperature (K) and surface air temperature (K), respectively.
The ground heat flux for the top soil layer is described askG(TS  T1 )D1131

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

where k is the soil thermal conductivity (W m-1 K-1), T1 is the soil temperature between
the first and second soil layers (K), and D1 is the thickness of the first soil layer (m). The
calculation of T1 is described in Liang et al., 1999.

When snow is present, the surface energy balance is solved at the snow/air interface
instead of soil/air interface. Ground heat flux is still computed, but it is the flux from thesnowpack into the ground. The snowpack is treated as two layers, thermally.
6.3.4 Routing Model

The routing model is described in detail by Lohmann et al. (1996, 1998a). It essentiallycalculates the concentration time for runoff reaching the outlet of a grid cell as well as thechannel flow in the river network. It is assumed that most horizontal flow within the gridcell reaches the channel network within the grid cell before it crosses the border into aneighboring grid cell. Flow can exit each grid cell in eight possible directions but all flow
must exit in the same direction. The flow from each grid cell is weighted by the fractionof the grid cell that lies within the basin. Once water flows into the channel, it does notflow back out of the channel and therefore it is removed from the hydrological cycle of
the grid cells. The daily surface runoff and baseflow produced by the VIC model fromeach grid cell is first transported to the outlet of the cell using a triangular unithydrograph, and then routed to in the river network to the basin outlet.
Both parts of the routing scheme (within grid cell and river routing) are constructed as
simple linear transfer functions. The routing model extends the FDTF-ERUHDIT (First
Differenced Transfer Function-Excess Rainfall and Unit Hydrograph by a Deconvolution
Iterative Technique) approach (Duband et al., 1993) with a time scale separation and asimple linear river routing model. The model assumes that the runoff transport is linear,causal, stable, and time invariant. It also assumes the impulse response function is nevernegative. The following summarizes the within grid and river network routing
respectively according to the modeling algorithms cited from Lohmann et al. ( 1996;
1998a).

6.3.4.1 Routing within a Grid Cell

To simulate the in-grid-dynamic of the horizontal routing process, one first separates thefast and slow components of the measured discharge with the linear model described in
Duband et al., (1993):
dQ S (t ) k  Q S (t )  b  Q F (t )dtwhere QS(t) is the slow flow, QF(t) is the fast flow and Q(t) is the total flow with Q(t) =QS(t) + QF(t) .
For each river basin, the parameters b and k are assumed to be constant over the period of
calculation. The ratio of b over k represents the ratio of water in the slow flow over waterin the fast flow. The fast and slow components are analytically connected by:132

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) ModeltS

Q (t )  b exp( k (t   ))Q F ( )d  Q S (0) exp(kt )0

The equation shows that the initial condition QS(0) decays exponentially with the mean
residence time of water in the flow (1/k) and the half-life decay is T1/2=(ln2)/k. With
discrete data the discharge equation can be solved with:exp(  k  t ) Sb  tQ S (t ) Q (t  t ) _Q(t )1  b  t
1  b  t
Based on the assumption that there is a linear relationship between measured streamflow
and effective precipitation (Peff, the part of the precipitation that becomes streamflow ), it
is sufficient to find an impulse response function connecting the fast component, QF, andPeff, due to the analytical connection of the fast and slow components. This impulse
response function and Peff can be found by solving the following equation iteratively:t max

Q F (t )   UH F ( ) P eff (t   )d0

In the equation UHF (τ) is the impulse response function (also called unit hydrograph) for
the fast flow component and tmax is the time taken for all fast processes to decay. The
equation for QF can be expressed in its discrete format, in which there are n data points atthe time step of Δt, and tmax = (m-1) ·Δt. Starting with the measured precipitation, thefollowing discrete equation is solved iteratively for the calculation of UH iF .
 QmF   Pmeff  P1eff UH 0F 

 


      .   F   eff


eff
F 
 Qn   Pn  Pn  m 1 UH m 1 
After each of the iteration steps the following constraint is applied:m 11
with UH iF  0 i
UH iF 

b
i0
1
k
The constraint results from the fixed fraction of the water in the fast and slow component,

the fact that  UH (t ) dt  1 and the non-negative assumption of UH(t). The calculated0

UHF is then put into the following discrete equation to solve for Peff.
 QmF  UH mF1  UH 0F 0  0  eff  P1     0      





 







0


 
 eff 
 Q F   0  0 UH F  UH F  Pn  n  m 10
Again, after each iteration step the constraint that ( 0  Pi eff  Precipitation, i ) is applied
after solving above equation.

133

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

The newly calculated Peff is then put back into the first discrete equation and thedeconvolutions are repeated until convergence is reached. Grid cell impulse response
functions can be obtained via deconvolution of the catchment impulse response function
with the river network impulse response function belonging to that catchment (Lohmann
et al., 1996).

6.3.4.2 River Routing

The transport of water in channels is described using a simple linear river routing model,
which follows the linearized Saint-Venant equation. The model assumes that water istransported out of the grid box only in the form of river flow. The following is thelinearized Saint-Venant equation, where C and D are parameters denote wave velocity
and diffusivity respectively.
Q
 2Q
Q
 D 2 Ct
x
x
Either from measurements or by estimation from geographical data of the river bed, C
and D are regarded as effective parameters since there are often times more than one riverin one grid cell. This way each grid cell ultimately ends up with one C and one D value,
which characterize the water transport within the cell.
The Saint-Venant equation is solved with convolution integralst

Q( x, t )   U (t  s) h( x, s )ds0

where
h ( x, t ) x(Ct  x) 2
exp(
)
4 Dt
2t tD

is the impulse response function of the Saint-Venant equation with h(x,0)=0 when x>0
and h(0,t)=δ(t) for t≥0. Because this solution scheme is linear and numerically stable, theinfluence from human activities (e.g., dams, irrigation water use) can be easily
implemented in each node.

6.3.5 Snow Model and Frozen Soil Algorithm

The main processes represented in the VIC snow model are shown schematically in
Figure 6.3. The spatial resolution for macroscale models usually ranges from 10 to 100
km, which is larger than the characteristic scales of the modeled snow processes.
Therefore, subgrid variability in topography, land cover, and precipitation are modeled bya mosaic-type representation, wherein each grid cell is partitioned into elevation (snow)
bands each of which contains a number of land cover tiles. The snow model is thenapplied to each land cover/elevation tile separately, and the simulated energy and mass134

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

fluxes, as well as the state variables for each grid cell are calculated as the area-averagesof the tiles.
The following are the main sub-models for a number of processes described in the snow
model.

Figure 6.3 Schematic of snow accumulation and ablation processes in the VIC snow model.

6.3.5.1 Snowpack Accumulation and Ablation

The snow model in VIC represents the snowpack as a two-layer medium, and solves anenergy and mass balance for the ground surface snowpack in a manner similar to othercold land processes models ( Anderson, 1976; Wigmosta et al., 1994; Tarboton et al.,1995). Energy exchange between the atmosphere, forest canopy and snowpack occurs
only within the surface layer. The energy balance of the surface layer is (Andreadis et.
al., 2009):
dWT S
 wcs
 Q r  Q S  Ql  Q p  Q m
dt
where cs is the specific heat of ice (J kg-1K-1, ρ w is the density of water (kg m-3), W is thewater equivalent (mm), Ts is the temperature of the surface layer (°C), Qr is the net
radiation flux (W m-2), Qs is the sensible heat flux (W m-2), Ql is the latent heat flux (Wm-2), Qp is the energy flux advected to the snowpack by rain or snow (W m -2), and Qm isthe energy flux given to the pack due to liquid water refreezing or removed from the packduring melt (W m-2). The detailed processes were described in Andreadis et al. (2009).

135Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

As snow accumulates on the ground, it goes through a metamorphism process, whichcauses the snowpack to compact and increase its density over time (except for depthhoar). In addition to the change in density caused by metamorphism, gravitational settlingcaused by newly fallen snow also contributes to the densification process. Following asimilar approach to Anderson (1976), compaction is calculated as the sum of two
fractional compaction rates representing compaction due to metamorphism andoverburden, respectively. Snow albedo is assumed to decay with age, based onrelationships published by the US Army Corps of Engineers (1956):
0.58

a a  0.85tad
0.46

a m  0.85tmd
where αa, αm are the albedo during the accumulation and ablation seasons, td is the timesince the last snowfall (in days), λ a = 0.92, and λm = 0.70. Accumulation and ablationseasons are defined based on the absence and presence of liquid water in the snow
surface layer, respectively.
The latest official version of the model contains an option to use the algorithm of Sun etal. (1999) instead. Unlike the VIC implementation of the US Army Corps algorithm, thismethod determines whether the snowpack is in “accumulation” or “ablation” mode basedon snowpack cold content, rather than on prescribed dates, and is therefore appropriate
for simulations anywhere in the world.

6.3.5.2 Atmospheric Stability

The calculation of turbulent energy exchange is complicated by the stability of theatmospheric boundary layer. During snowmelt, the atmosphere immediately above thesnow surface is typically warmer. As parcels of cooler air near the snow surface aretransported upward by turbulent eddies, they tend to sink back toward the surface where
turbulent exchange is suppressed. In the presence of a snow cover, aerodynamic
resistance is typically corrected for atmospheric stability according to the bulk
Richardson’s number which is a dimensionless ratio relating the buoyant and mechanicalforces acting on a parcel of air (Anderson, 1976). While the bulk Richardson’s number
correction has the advantage of being straightforward to calculate based on observations
at only one level above the snow surface, previous investigators have noted that its usage
results in no turbulent exchange under common melt conditions and leads to anunderestimation of the latent and sensible heat fluxes to the snowpack (e.g. Jordan 1991;
Tarboton et al. 1995).

6.3.5.3 Snow Interception and Canopy Effects

The snow interception algorithm in the snow model represents canopy interception,
snowmelt, and mass release at the spatial scales of distributed hydrology models. During
each time step, snowfall is intercepted by the overstory up to the maximum interception
storage capacity according to:136

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) ModelI = fPs
where I is the water equivalent of snow intercepted during a time step (mm), Ps is thesnowfall over the time step (mm), and f is the efficiency of snow interception (taken as0.6) (Storck et al., 2002).
The maximum interception capacity (mm), B, is given by:
B = Lrm(LAI)
where LAI is the single-sided leaf area index of the canopy and m is determined based on
observations of maximum snow interception capacity (mm). Lr is the leaf area ratio
which is a function of temperature (Andreadis et al. 2009).
Snowmelt is calculated directly from a modified energy balance, similar to that applied
for the ground snowpack. Newly intercepted rainfall is calculated with respect to thewater holding capacity of the intercepted snow, which is given by the total capacity of thesnow. The bare excess rainfall then becomes throughfall. The intercepted snowpack cancontain both ice and liquid water. Snowmelt in excess of the liquid water holding
capacity of the snow results in meltwater drip. Mass release of snow from the canopy
occurs if sufficient snow is available and the ratio of 0.4 is derived from observations of
the ratio of mass release to meltwater drip (Storck et al., 2002).

6.3.5.4 Blowing Snow

The blowing snow algorithm was developed by Bowling et al. (2004) to estimate
topographically-induced sub-grid variability in wind speed, snow transport andsublimation. The blowing-snow algorithm is designed to work within the structure of theexisting VIC mass and energy-balance snow model. The algorithm accounts for theenergy advected by rainfall, throughfall, or drip (when overstory is present), as well as
net radiation, ground heat flux, and sensible and latent heat fluxes. Incoming shortwave
and longwave radiation and wind speed are attenuated through the canopy, if present. Ifsnow is present, it is assumed to cover the understory for purposes of radiation transfer.
For each vegetation fraction within the grid cell, the time rate of change of snow water(We) is:dWe P  M  p  Qv  Qe
dt
dW e
where
is the rate of snow water accumulation, P is precipitation, M is snowmelt
dt
and drainage, Qv is the sublimation from blowing snow, and Qe is evaporation andsublimation from the snowpack, for a time increment dt. All of the terms are in units ofmillimeters per time step. The spatial probability of occurrence of blowing snow, p, isunitless.
Along with the standard meteorological forcings, three additional parameters are needed
to run the blowing snow algorithm: standard deviation of terrain slope, standard deviationof terrain elevations, and the lag-one autocorrelations for each model grid cell. S

137

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

6.3.5.5 Snow Model Calibration

During the calibration for the VIC snow model, four parameters are adjusted for gridcells: 1) maximum air temperature at which snowfall occurs; 2) minimum air temperature
at which rainfall occurs; 3) the snow surface roughness; 4) the value of m, which controls
the maximum snow interception capacity as a function of LAI.
oo

Usually the first two parameters are set to 0.5 C and -0.5 C, respectively. In addition,
we suggest that the snow roughness parameter should be in the range from 0.001 m to0.03 m.
The VIC snow model is intended primarily for large-scale applications. It has been
incorporated as the standard snow scheme within the VIC model, which represents subgrid spatial variability by simulating state and fluxes in land cover/elevation tiles. Within
the VIC model, it is used in a real-time hydrologic forecast system for the western U.S.
(Wood and Lettenmaier, 2006), and has been used in numerous analyses, diagnoses, andpredictions of climate variability and change (e.g. Christensen and Lettenmaier, 2007).

6.3.6 Frozen Soil Algorithm

Over the cold regions, the soil ice content of the frozen soil directly affects infiltration,and indirectly affects the heat transfer to and from the overlying snowpacks. A frozen soil
algorithm (Cherkauer and Lettenmaier 1999, 2003; Cherkauer et al., 2003; Bowling et al.,2008) has been implemented into the VIC model to improve its modeling skills over thecold regions. The frozen soil algorithm uses the same soil moisture transport scheme as
described in Section 6.3.2.2, while the frozen soil penetration is calculated by solving thethermal fluxes through the soil column. For each time step, the thermal flux through thesoil column is solved first to determine the soil layer ice content. Moisture fluxes are thencomputed using the ice content.
The heat flux through the soil column is

T

T
Cs
 (k
)  i L f ( i )t z ztwhere k is the soil thermal conductivity (Wm-1K-1), Cs is the soil volumetric heat capacity(J m-3K-1), T is the soil temperature (°C), ρi is the ice density (kg m-3), Lf is the latent heatof fusion (J kg-1), θi is the ice content of the layer (m3 m-3), t is time (s), and z is the depth(m). The last term of the equation only applies when the soil is frozen.
Within the soil column, a number of nodes are specified by the user. There is a node at
the surface, a node at the bottom of layer 1, and a node in the middle of layer 1. There isa node at a user-specified maximum depth (specified in the soil parameter file, typically4m), but this node need not be at the bottom of the soil column. All remaining nodes arespaced evenly between the bottom of soil layer 1 and the user-specified maximum depth.
The model also has an option to space the nodes exponentially ( Adam, 2007), i.e. close
together near the surface and gradually further apart with depth, down to the userspecified maximum depth. This is good for simulating permafrost, for which it is often
necessary to specify a maximum depth of as much as 40m.

138

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

With the nodes specified, the soil temperature is then solved for numerically at an hourly
time step via an explicit finite difference approximation of the soil thermal flux equation(see Cherkauer and Lettenmaier 1999 for detailed numerical expressions).
Thermal conductivity (k) and volumetric heat capacity (Cs) of the soil layer are calculated
at each time step after the soil moisture and ice content are updated by the moisture flux
solution. The soil thermal conductivity (k) is computed after a modification of Farouki(1986) as:k  ( k sat  k dry )k e  k dry
where ksat and kdry are the thermal conductivity of saturated soil and dry soil (W m-1 K-1),
respectively; and ke is the kersten number which weighs the two soil conductivities.
The volumetric heat capacity Cs is computed by summing the volumetric heat capacitiesof the soil constituents (Flerchinger and Saxton, 1989):
C s    j c j jwhere ρj, cj, and θj are the density (kg m -3), specific heat capacity (J m -3K-1), andvolumetric fraction of the jth soil node, respectively.
By adding the ice content component in the heat flux equation, the impact of frozen soil
on moisture transport can be simulated by the moisture flux algorithm. The first way thatthe ice content in the frozen soil affects the moisture transport is through available
moisture storage. Each of the three soil layers is divided into thawed, frozen, andunfrozen sublayers. The thickness of these sublayers depends on the soil temperatures at
the nodes. When there is a frozen layer present, the ice content is based on the averagetemperature of the sublayer. The fraction of the unfrozen water as by Flerchinger andSaxton (1989) is;Bp 1 
Lf T



Wi  Wi  g e  T  273.16 
where Wi is the liquid water content of soil layer i (mm), Wic is the maximum watercontent of soil layer i (mm), g is acceleration due to gravity (m s-2), ψc is the air entry
potential (m), and Bp is the pore-size distribution.
The second way the ice content affects soil moisture transport is through its effect on
infiltration and drainage. When a soil layer has high ice content, on one hand, it will benearly saturate to the runoff calculations, but on the other hand there is little moisture(unfrozen) to be allowed to drain to the lower layer.
c

The frozen soil model also has an option to simulate excess ground ice (Adam, 2007), a
common feature of permafrost. The user specifies the amount of excess ground ice as an
optional soil parameter and then VIC computes new effective densities and porosities.
As this ice melts, the effective porosities and densities approach the non-excess-ice
values (and the melt water is added to the soil; any excess runs off as necessary).

139

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

6.3.7 Lake and Wetland ModelIn the VIC model, the effects of lakes and wetlands are simulated by creating alake/wetland tile that can be added to the grid cell mosaic, in addition to the vegetation
and bare soil tiles (Bowling and Lettenmaier, 2009). The lake/wetland tile represents
seasonally flooded ground as well as permanent water bodies. The tile contains a body ofopen water (lake) whose areal extent is allowed to change in response to the lake waterbalance. The wetland portion of the tile is the (time-varying) remaining portion of thelake/wetland tile not covered by the lake. Water and energy components of the combinedlake and wetland are resolved at each model time step. The energy balance of the lake
component builds on the work of Hostetler and Bartlien (1990), Hostetler (1991), andPatterson and Hamblin (1988), while that of the exposed wetland follows Cherkauer andLettenmaier (1999).
There are a few limitations with the current version of the lake and wetland model. First,
the vegetation in the wetland portion is prescribed to be shrubs that are
typical of tundra vegetation. Second, the model only simulates lakes that receive all of
their inflows from within the same grid cell, i.e. no channel inflows from other grid cells.
This restricts the lakes to just those small lakes whose drainage basins are contained
within the current grid cell. Third, the wetlands in the current version are essentiallyuplands that could potentially be flooded by lake expansion, with no special wetland
processes considered. Only wetlands formed by seasonal flooding due to local
precipitation and snowmelt and poor drainage are modeled.
The description of this section is primarily cited from Bowling and Lettenmaier (2009).

6.3.7.1 Lake Algorithm

Evaporation from the water surface is calculated in each time step by solving a surface
energy balance using the formulations by Hostetler and Bartlien (1990), and Hostetler(1991). The energy exchange with the atmosphere occurs within the surface water layer,
which is limited to a user-specified depth (zsurf), typically around 0.6 m. The absorptionof solar radiation by the surface water layer is assumed to follow Beer’s law. The
radiation intensity at the depth h is assumed to be a two-band system and expressed as(Patterson and Hamblin 1988):
I (h)  I o Av  exp(v h)  ANIR  exp( NIR h)
where Io is the net shortwave radiation at the water surface (Wm-2), A v and ANIR are thefractions of total radiation in the visible and near-infrared bands, respectively, and λv andλNIR are the attenuation coefficients of the two bands. Av and ANIR are set to 0.7 and 0.3,
respectively.
In deeper lakes, the average temperature for additional water layers is resolved by solving
a set of simultaneous equations. Included in these equations are the effects of radiation
absorption by each water layer, the eddy diffusion of heat from adjacent layers bymolecular diffusion, wind-induced turbulent mixing, and convective mixing due to140Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Modeltemperature instabilities (Hostetler and Bartlien 1990). The bottom boundary has a noflux condition, meaning that energy is not exchanged with the sub-lake soils (thiscondition could be relaxed in the future).
The lake layer thickness is recalculated in each time step in response to variations in totallake liquid water depth, as follows:d zsurf  z; d z  0; l  1
d zsurf  z / 2; d z  z / 2; l  2
d zsurf  z surf ; d z  z  z surf  l  1; l  int z / z surfz  z surfz surf  z  2 z surf
2 z surf  z  N nodes z surfd zsurf  zsurf ; d z  z  zsurf  N nodes  1; l  N nodesz  N nodes z surf

Where z is the current lake depth (m), l is the current number of solution layers, and dz isthe thickness of all solution layers (m) excluding the surface layer. The current thicknessof the surface solution layer is dzsurf (m), and zsurf is the maximum allowable thickness ofthe surface solution layer (m). The maximum number of computational layers in the lake,
Nnodes, is a user-specified parameter.
Freezing and thawing of the lake ice are represented using the method by Patterson andHamblin (1988). Snow accumulation and melt over the frozen surface is solved using the
VIC two-layer energy balance snow model (Cherkauer et al. 2003). In this case, heat flux
out of the lake takes the place of the ground heat flux. A two-band solar radiation
absorption model similar to the equation for the unfrozen lake case is applied to the snow
and ice layers. Heat flux through the ice is driven by the temperature gradient, and must
balance the heat flux from the lake and the energy of ice formation at the ice/waterinterface (Patterson and Hamblin 1988). In the VIC model, the water equivalent of lake
ice is a state variable, and the available liquid water volume is checked before new ice
can form (which is important in shallow wetland systems). For stability, lake ice must
exceed a user-specified minimum thickness (usually taken to be 10 cm). As ice melts, thearea of lake ice is adjusted to maintain this minimum thickness, resulting in fractional ice
coverage if ice area is less than the surface area of liquid water.

6.3.7.2 Lake/Wetland dynamics

Unique features of the VIC lakes and wetland algorithm include the interaction of thesimulated lake within the VIC model grid cell and the ability to represent wetlands ofvarying size. The algorithm can be summarized as follows (see Figure 6.4):
All open water areas within a VIC model grid cell are simulated together as aneffective grid cell lake.

A user-defined fraction of runoff from vegetated areas within the grid cell isdiverted to the lake. This represents the storage retardation effect of lakes on
seasonal streamflow.

141Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Once the new lake level is calculated, runoff is released from the lake as afunction of the lake level. Base flow is calculated from below the lake as afunction of the liquid water content of saturated wetland soils.

Specification of a variable depth-area relationship allows for the representation of
the reduction in surface water extent and the emergence of wetland vegetation
following drainage of seasonally flooded wetlands.

Figure 6.4 Schematic of the VIC lake and wetland algorithm. I: Evaporation from the lake iscalculated via energy balance, II. Runoff enters the lake from the land surface, III: Runoff out of
the lake is calculated based on the new stage, and IV: The stage is re-calculated.

In surface hydrology, the term wetland specifically means areas that are saturated or
covered by water for some portion of the year (Zoltai 1979). The tendency of a region toflood periodically can be represented within the VIC model by a user-input depth-area
relationship, A(z), for the maximum inundated fraction of the grid cell. For clarity of thesubsequent discussion, wetland fraction, Cwet, will be used to refer to the maximum
fraction of the VIC model grid cell that can be flooded, while lake fraction, flake, refers to
the fraction of Cwet that is inundated for a given time step. In this context, therefore, the142Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Modellake fraction does not distinguish between the pelagic open water zone and the benthic
zone that may contain emergent wetland vegetation.

Figure 6.5 Schematic for the wetland algorithm: a) when the lake is at its maximum extent thesoil column is saturated, b) as the lake shrinks runoff from the land surface enters the lake and c)evaporation from the land surface depletes soil moisture, d) as the lake grows, water from the lake
recharges the wetland soil moisture

As the stage of the simulated lake drops, the open water area is recalculated andadditional wetland area is exposed (Figure 6.5). The energy and water balance of the
newly exposed wetland is solved as an additional vegetation tile. The water balance of
the wetland fraction (Cwet) can be represented as follows:S  P  Dveg  E w  f lake  Ev  (1  f lake ) Dlake
where ΔS is the change in soil moisture, lake water and ice and snow storage (mm); P isprecipitation (mm); E w and Ev are evaporation from the open water and wetland
vegetation (mm) respectively. Dveg is the discharge (runoff and baseflow) entering into
the lake from the non-wetland portion of the grid cell (mm) and Dlake is the discharge outof the lake (mm). All of the runoff and baseflow generated by the exposed wetland area isassumed to enter the grid cell lake, so this internal transfer is not needed in the wetland
water balance.
The soil column under the lake is assumed to be saturated. As the lake area is reduced,
the soil moisture of the non-lake area is updated to include the newly exposed fraction ofsaturated soil (Figure 6.5b and Figure 6.5c). Likewise, as the lake area expands, some of
the lake volume must go to saturating the newly inundated soil (Figure 6.5d). The
average soil moisture for the wetland fraction is therefore updated as follows:  f lakeW  (1  f lake )  Wv  f lake  Wmax f lake )  Wv  f lake  WmaxW  (1  f lake  f lakef lake143

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

where f’lake and f lake represent the new and old lake fractions, respectively. Wv is thewetland soil moisture for the current time step (mm) and Wmax is the total soil moisturestorage capacity of the soil column (mm). The volume of water (as mm per unit area)from the expanding lake used to recharge the wetland soil moisture is calculated asfollows:
  f lake )(Wmax  Wv ) f lake  f lake
Recharge  ( f lake
Changes in lake stage are calculated via a water balance for the saturated lake area.
Runoff into the lake is composed of all of the runoff and baseflow from the exposed
wetland and a fraction of the runoff and baseflow from all other grid cell vegetation tiles.
To avoid complications due to the variation of lake area with depth, lake volume is thestate variable used for the water balance. Lake depth is updated each time step bypiecewise integration of the derived depth volume curve. Subsurface outflow from thelake is calculated using the VIC model Arno baseflow curve. Since the sub-lake soil
thermal regime is not resolved, the maximum moisture storage in the bottom soil layer isreduced by the ice content of the exposed wetland soil profile in order to calculate thebaseflow, as described by Cherkauer and Lettenmaier (1999).
Surface outflow from the lake is calculated as a function of the new depth, based on theequation for flow over a broad-crested weir, assuming that the velocity head is negligible:
32
Q  c d b g ( ( z  z min ) 2 )3
Where Q is the discharge (m3 s-1), b is the flow width (m), g is the acceleration due togravity, z is the current lake depth (m) and zmin is the elevation above the lake bottom of
the weir or lake outlet (m). The coefficient of discharge, cd, is used to account for thevelocity of approach, non-parallel streamlines over the crest, and energy losses. cd varies
between about 0.8 and 1.2 and frequently has a value of about 0.94, which was adopted
here (Hamill, 2001). Surface runoff out of the lake (m) becomes:
32

1.6b( z  z min )  dt
A( z )
Where dt is the time step length in seconds and A(z) is the lake surface area (m2) at depthz (m). For natural lakes and wetlands the width of the reservoir outlet is also likely tovary with water level. Assuming a roughly circular lake, the flow width, b, can be
expressed as a fraction of the lake circumference:b  2 f   A( z )where f is the fraction of the lake circumference. Surface runoff out of the lake becomes:
Rlake ( z ) 

0

3/ 2

Rlake  5.67 f z  z min  dt

Az 


z  z minz  z min144

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Lake depth, z, is calculated as the depth of both liquid water and lake ice water equivalent
when liquid water exceeds the lake ice water equivalent, to account for the displacementof water by floating ice. When the mass of ice exceeds the mass of liquid water, z is thedepth of liquid water alone. The depth-area curve (A(z)), width fraction (f); and minimum
allowable depth of the lake (zmin)are input parameters to the lake model. The width
fraction, f, can take on values from 0 - 0.5, and is typically adjusted during thecalibration.

6.3.8 Irrigation Scheme and Reservoir Module

The VIC irrigation modeling framework was developed to represent the effect ofirrigation on the water balance of large continental rivers. It is coupled with the VIC
model by including a sprinkle irrigation scheme (Haddeland et al., 2006a; 2006b). The
sprinkle irrigation is based on a standardized method of irrigation scheduling andinformation about growing season and irrigation intensity given by the United Nations
Food and Agriculture Organization’s (FAO) database AQUASTAT (FAO, 2003). The
reservoir model uses the VIC model as its centerpiece by adding a reservoir module(Haddeland et al., 2006a; 2006b) in the Lohmann et al. (1996, 1998a) routing model. The
main references for this section are Haddeland et al. (2006a; 2006b).
The energy balance mode of the VIC model is utilized, which means that the model
iterates for the surface temperature to reach closure of the surface energy and water
budgets at each time step. Required minimum input data for the model are daily
precipitation, and maximum and minimum daily temperatures. When radiation and vapor
pressure data are not supplied to the model, VIC calculates these variables based on daily
precipitation and daily minimum and maximum temperatures, using algorithms
developed by Thornton and Running (1999), and Kimball et al. (1997) as described in
Nijssen et al. (2001b). If wind speed or atmospheric air pressure are not provided, themodel uses default values (1.5 m s -1 and 95.5 kPa).
The main purpose of irrigation is to avoid vegetation stress caused by limited soil
moisture availability. The VIC model was therefore modified to allow for irrigation wateruse, based on the model’s predicted soil moisture deficit. Irrigation starts when soil
moisture drops below the level where transpiration becomes limited, and continues until
soil moisture reaches field capacity. Grid cells in which irrigation occurs are partitioned
into an irrigated part and a non-irrigated part, based on the fractional area irrigated within
the cell (Siebert et al., 2002).
Crop characteristics are determined according to FAO’s guidelines for computing crop
evapotranspiration (FAO, 1998a). Reference crop evapotranspiration is first calculated
within each model grid cell based on the Penman–Monteith equation ( Shuttleworth,
1993, see Section 6.3.2 for details). Crop coefficients and heights specified by FAO are
thereafter used to calculate LAI values throughout the growing season. The crop
coefficients have already taken into account soil evaporation as part of the waterrequirements. Crops with crop coefficients calculated in this way are assigned to the145Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

irrigated part of the grid cell, and the remaining vegetation is assigned to the nonirrigated part.
Storage in reservoirs can affect the streamflow significantly, and for this project areservoir module was developed and included in the Lohmann et al. (1996, 1998a)routing model. Reservoir characteristics and operating purposes were taken from theinternational commission on large dams (ICOLD) (ICOLD, 2003). An optimization
scheme based on the SCEM-UA algorithm ( Vrugt et al., 2003) was used to calculate
optimal releases given reservoir inflow, storage capacity, and downstream water or power
demands (see also Table 6.1 and Figure 6.6). Irrigation water can be extracted from river
runoff locally, or, in periods of water scarcity, from reservoirs or any other prescribed
point in the river basin. In this case, irrigation is restricted by water availability.
Alternatively, irrigation water is assumed to be freely available, and the model simulates
irrigation water requirements. In this case, irrigation is not restricted by water
availability, and it is hence possible that more water is used for irrigation than is available
in the river basin. The VIC model, like most land surface schemes, does not represent
groundwater in a way suitable for modeling groundwater withdrawals. A single-reservoir
algorithm is used — that is, it does not consider the simultaneous operation of multiple
reservoirs in a river basin. The reservoir model was run at a daily time step. However,water demands were calculated on a monthly basis, and within each month releases were
kept constant if possible. The economic value of reservoir releases for hydropower andwater supply was assumed to be constant throughout the year.
Table 6.1 Objective Functions Used in the Reservoir Model a

PurposeIrrigation

Objective Function365

min  (Qd i  Qri ), Qd  Qr
i 1

365

Flood control
min  (Qri  Q flood ) 2 , Qr  Q flood
i 1
365

Hydropower1
i 1 Qri  hg

min 
365

Water supply, navigationmin  (Qri  Qmean )i 1a

Qd: water demands (km3); Qr: reservoir releases (km3); Qflood: mean annual flood (km3),
calculated based on simulated naturalized discharge; Q mean: mean annual flow (km3); r:
density of water (kg m -3); η: efficiency of the power generating system; h: hydrostaticpressure head (KPa); g: acceleration due to gravity (9.8 m s -²).
Irrigation demands are calculated based on simulated irrigation water requirements
downstream of the reservoir — that is, the grid cell elevation must be lower than that of
the reservoir grid cell and lower than the five maximum grid cells from the reservoir’s
downstream river course. If there are multiple reservoirs upstream of an irrigated area,
but the reservoirs themselves are located in separate tributaries, demands are divided
based on reservoir capacity. For reservoirs located on the same river course, irrigation

146

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Modeldemands for shared downstream areas are used to represent water demands for all
reservoirs. Flood damages are expected when river discharge exceeds bankfull discharge,
which has a recurrence interval on the order of once in 1.5 to 10 years ( Mosley and
McKerchar, 1993). The mean annual flood (the mean of the annual maximum daily
discharges) may be used as a rough approximation of bankfull discharge. For hydropower
reservoirs, the optimization scheme is used to maximize hydropower production. When areservoir has multiple purposes, irrigation demands are given priority, followed by flood
control. Any excess water is used to maximize hydropower production, if applicable.
The reservoir module is retrospective ─ that is, it assumes perfect knowledge of future
reservoir inflows. At the beginning of each operational year, the next 12 month’s inflows
are used to determine reservoir releases. The start of the operational year is defined as thetime when mean monthly simulated naturalized streamflow shifts from being higher thanthe mean annual flow to being lower than the mean annual flow, following theconvention of Hanasaki et al. (2006). Minimum release (Q min, km3) can be set as theseven-day consecutive low flow with a ten year recurrence period, and is calculated basedon simulated naturalized flow at the reservoir location. The maximum volume of water(Qmax, km3) released for the current day (i) can be written as:
365
365
365


Qmax i  min ( S i 1  Qini ), ( S i 1  S end   Qinday   Qmin   E resday )
day  iday i 1day iwhere Si-1 is reservoir storage at the end of previous day (km 3), Send is storage at the endof the operational year (km3), Qin is simulated inflow to the reservoir (km3), and Eres isreservoir evaporation (km3), which is calculated using the Penman equation. S end varies
between 60 and 80 percent of maximum reservoir capacity, depending on water demands
during the current 12-month simulation period.
The model can be run assuming water availability is not a limiting factor, in which case it
calculates irrigation water requirements (which is defined as the water required inaddition to water from precipitation (soil moisture) for optimal plant growth during thegrowing season, (i.e., the difference between potential evapotranspiration and actual
evapotranspiration). When the modeling scheme takes into account how much water isactually available (locally or in upstream reservoir(s) built for irrigation purposes),consumptive irrigation water use is calculated. For the dam datasets, although thereservoir capacity information is complete, the reservoir's surface area is lacking over
22% of the dams. For these reservoirs, the relationship between reservoir volume andsurface area according to Takeuchi (1997) is used:V  9.208 A1.114
where A is the surface are (km2), and V is the capacity of the reservoir in 10×6 m3. Allreservoirs are assumed to have rectangular cross-sections, which is used to calculate
height-storage relationships for the reservoirs built for hydropower purposes.
In order to run the VIC model with the reservoir module implemented, the model first has
to be run to simulate natural conditions, i.e. without the reservoir module and irrigation
scheme implemented. Based on the naturalized simulations, the seven-day consecutive

147

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

low flow with a ten year recurrence period (Q7day), the mean of the annual maximum
daily discharges (Qflood), and the start of the hydrological year are calculated. The
streamflow information is used as input to the reservoir model. In addition, downstream
irrigation water requirements are required as an input to the reservoir model. The VIC
model is run for all grid cells upstream of the reservoir. The streamflow (Q in) is thenrouted to the reservoir, before optimal releases from the reservoir are calculated. When areservoir is built for multiple purposes, a combination approach is taken. Water
withdrawals from reservoirs are based on simple rules intended for implementation in any
river basin. As described above, the elevation of the grid cell in need of water has to belower than the elevation of the reservoir. Water is only extracted from the reservoir when
there is not enough water available locally. Another generalization is that upstream
locations are given priority at the cost of the possible needs of downstream locations. The
irrigation scheme is illustrated in Figure 6.6.

Figure 6.6 Schematic representation of the VIC irrigation scheme. The model grids and routing
network are shown on the left, and an example grid cell is shown on the right. Water is extracted
from the river and reservoir, and applied to the irrigated part of the cell.The excess water returns
to the river system.

Validation runs have shown that the model simulates irrigation water requirements thatare close to the reported ones, and the model is able to capture the main hydrologic
effects of reservoir operations and irrigation water withdrawals (Haddeland et al., 2006a,b). The reservoir model and irrigation scheme are described in more detail in Haddeland
et al. (2006a, b, 2007).

6.4

VIC Model Parameters and Forcings

Land surface characteristics required by the VIC model include soil data, topography, andvegetation characteristics. These parameters as well as the model forcing data aredescribed in this section. The locations of these parameters and forcings, as well as themodes of operation (e.g., energy balance mode or not), are described by a global
parameter file.

148

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

6.4.1 Soil Parameters

The soil parameter file used by VIC describes the unique soil properties for each grid cellin the model domain. It is also the main file that identifies which grid cells will besimulated, and what their latitudes and longitudes are (which is used to find the forcing
files for the grid cells).
Soil texture information and soil bulk densities were derived from the 5-min Food andAgriculture Organization data set (FAO, 1998b). The soil parameters generally fall into
two categories. The first category of soil parameters is not adjusted once it is determinedfrom the FAO data. These parameters include porosity θs (m3m-3), saturated soil potentialψs (m), saturated hydraulic conductivity ksat (ms-1), and the exponent B for unsaturated
flow (which was based on Cosby et al., 1984).
Another category of soil parameters is subject to calibration based on the agreement
between simulated and observed hydrographs. Parameters in this category include thethickness of each soil layer, d i; the exponent of the infiltration capacity curve, bi; and thethree parameters in the baseflow scheme: Dm, Ds, and Ws. The soil texture is based on a
5-min Food and Agriculture Organization dataset (FAO 1998). The specific soil
characteristics (e.g., field capacity, wilting point, and saturated hydraulic conductivity)
were obtained from algorithms by Cosby et al. (1984), Rawls et al. (1998), and Reynolds
et al. (2000) for each soil texture type.
Soil hydrologic and thermal parameters needed for the different VIC model set-ups arelisted below in Table 6.2. All columns of the input file must be filled, but certain
parameters need only be defined if the full energy or frozen soil models are activated.
Table 6.2 Soil parameters

VariableNamegridcellatloninfilt
Ds
Dsmax
Ws

c

Units

Number DescriptionofValuesN/A
1
1 = Run Grid Cell, 0 = Do Not
Run
N/A
1
Grid cell number
degrees 1
Latitude of grid cell
degrees 1
Longitude of grid cell
N/A
1
Variable infiltration curveparameter (bi)fraction 1
Fraction of Dsmax where nonlinear baseflow begins
mm/day 1
Maximum velocity of baseflow
fraction 1
Fraction of maximum soil
moisture where non-linear
baseflow occurs
N/A
1
Exponent used in baseflow curve,
normally set to 2

149

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

expt

N/A

NlayerKsat
phi_s
init_moist
elev
depth

mm/day
mm/mm
mm
m
m

NlayerNlayerNlayer
1
Nlayer

avg_T

C

1

dp

m

1

bubblequartzbulk_density
soil_densitycm
fractionkg/m 3kg/m 3

NlayerNlayerNlayerNlayer

off_gmt
Wcr_FRACT

hours
fraction

1
NlayerWpwp_FRACT fraction

Nlayerroughsnow_roughannual_prec
resid_moistmm
mm
fraction

1
1
1
Nlayerfs_active

1 or 0

1

Parameter describing the variationof Ksat with soil moisture
Saturated hydrologic conductivity
Soil moisture diffusion parameter
Initial layer moisture content
Average elevation of grid cell
Thickness of each soil moisturelayer
Average soil temperature, used as
the bottom boundary for soil heat
flux solutions
Soil thermal damping depth (depth
at which soil temperature remains
constant through the year, ~4 m)
Bubbling pressure of soil
Quartz content of soil
Bulk density of soil layer
Soil particle density, normally
2685 kg/m3
Time zone offset from GMT
Fractional soil moisture content at
the critical point (~70% of field
capacity) (fraction of maximum
moisture)
Fractional soil moisture content at
the wilting point (fraction of
maximum moisture)
Surface roughness of bare soil
Surface roughness of snowpack
Average annual precipitation.
Soil moisture layer residual
moisture.
If set to 1, then frozen soil
algorithm is activated for the gridcell. A 0 indicates that frozen soils
are not computed even if soil
temperatures fall below 0°C.

OPTIONAL (include if JULY_TAVG_SUPPLIED = TRUE in global parameter file):
VariableName
July_Tavg

Units
C

Number DescriptionofValues1
Average July soil temperature,

150

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

used for treeline computations.

6.4.2 Vegetation Parameters

Land cover characterization was based on the University of Maryland global vegetation
classifications described by Hansen et al. (2000), which has a spatial resolution of 1 km,
and a total of 14 different land cover classes. From these global data, the land cover types
present in each grid cell in the model domain and the proportion of the grid cell occupied
by each are identified, as described by Maurer et al. (2001). The primary characteristic of
the land cover that affects the hydrologic fluxes simulated by the VIC model is LAI. The
LAI is derived from the gridded monthly global LAI database of Myneni et al. (1997),
which is inverted using the Hansen et al (2000) land cover classification to derive
monthly mean LAIs for each vegetation tile for each grid cell. The LAI values do notchange from year to year in this implementation of VIC; hence, interannual variations invegetation characteristics are ignored. Furthermore, the Myneni et al. (1997) LAI values
to which the method is tied are based on averages over the period 1981–1994, which may
not be representative of the entire simulation period. Rooting depth is specified for eachland use type so that shorter crops and grasses draw moisture from the upper soil layers,
and tree roots from the deeper soil layer (e.g., Jackson et al. 1996). Additional parameters
for each vegetation tile were assembled based on several sources, including roughness
length and displacement height (Calder, 1993), architectural resistance (Ducoudré et al.,1993), and minimum stomatal resistance (DeFries and Townshend, 1994).
The vegetation parameter file describes the vegetative composition of each grid cell, anduses the same grid cell numbering as the soil file (latitudes and longitudes are notincluded in the file). This file cross-indexes each vegetation tile (from any land-coverclassification scheme) to the classes listed in the vegetation library. Vegetation
parameters and vegetation library for the VIC model are listed below in Table 6.3 andTable 6.4, respectively.
Table 6.3 Vegetation Parameters

Variable Name

Units

Description

vegetat_type_num
N/A
Number of vegetation tiles in a grid cell
Repeats for each vegetation tile in the grid cell:
Variable Name
Units
Description
veg_class

N/A

Vegetation class identification number (reference
index to vegetation library)Cv
fraction
Fraction of grid cell covered by vegetation type
Repeats for each defined root zone, within the vegetation tile:
Variable Name
Units
Description

151

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Variable NameUnitsroot_depthmDescription
Root zone thickness (sum of depths is total
depth of root penetration)

root_fract
fraction
Fraction of root in the current root zone.
OPTIONAL (include if BLOWING_SNOW = TRUE in global parameter file) - Include
for each vegetation tile:
Variable Name
Units
Description
sigma_slope

N/A

Standard deviation of terrain slope for eachvegetation class

lag_oneN/A

Lag one gradient autocorrelation of terrain
slope

fetch

m

Average fetch length for each vegetation class

OPTIONAL (include if GLOBAL_LAI = TRUE in global parameter file) - Include foreach vegetation tile:
Variable Name
Units
Description
GLOBAL_LAI

N/A

Leaf Area Index, one per month

Table 6.4 Vegetation Library

Variable
Name

Units

Numberof Values

veg_class

N/A

1

Vegetation class identification number(reference index for library table)

overstory

N/A

1

Flag to indicate whether or not the current
vegetation type has an overstory (TRUEfor overstory present [e.g. trees], FALSE
for overstory not present [e.g. grass])

rarc

s/m

1

Architectural resistance of vegetation type(~2 s/m)

rmin

s/m

1

Minimum stomatal resistance ofvegetation type (~100 s/m)

12

Leaf-area index of vegetation type

LAI

Description

albedofraction 12

Shortwave albedo for vegetation type

rough

M

12

Vegetation roughness length (typically
0.123 * vegetation height)

displaceme Mnt

12Vegetation displacement height (typically
0.67 * vegetation height)

wind_h1

Height at which wind speed is measured.

M

152

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Variable
Name

Units

RGL

W/m^2 1Minimum incoming shortwave radiation atwhich there will be transpiration. For trees
this is about 30 W/m^2, for crops about
100 W/m^2.

rad_atten

fract

1

Radiation attenuation factor. Normally set
to 0.5, though may need to be adjusted forhigh latitudes.

wind_atten fract

1

Wind speed attenuation through theoverstory. The default value has been 0.5.

trunk_ratio fract

1

Ratio of total tree height that is trunk (nobranches). The default value has been 0.2.

comment

1

Comment block for vegetation type.
Model skips end of line so spaces are valid
entrys.

N/A

Numberof Values

Description

6.4.3 Elevation Band

This file contains information needed to define the properties of each elevation band usedby the snow model. It is only needed when the snow mode is set “TRUE” in the global
parameter file. Snow elevation bands are used to improve the model's performance inchanging topography, especially mountainous regions where the effects of elevation on
snowpack accumulation and ablation might be lost in a large grid cell. The number ofsnow elevation bands to be used within the model is defined in the model control file. This
file is only read if the number of snow elevation bands is greater than 1. Parameters for theelevation band are listed in Table 6.5.
Table 6.5 Elevation Band parameters

Variable Name

Units

Number of
ValuesDescription

Cellnum

N/A

1

Grid cell number (shouldmatch numbers assigned insoil parameter file)

AreaFract

fraction

SNOW_BAND Fraction of grid cell coveredby each elevation band. Sumof the fractions must equal 1.

153

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Variable Name

Units

Number of
ValuesDescription

elevationm

SNOW_BAND Mean (or median) elevation ofelevation band. This is used tocompute the change in air
temperature from the grid cellmean elevation.

Pfactorfraction

SNOW_BAND Fraction of cell precipitationthat falls on each elevation
band. Total must equal 1. To
ignore effects of elevation on
precipitation, set thesefractions equal to the area
fractions.

6.4.4 Lake parameters

The lake parameter file is also optional ─ it is only needed when “LAKES = TRUE” inthe global parameter file. In addition to the information about inflows and outflows, theuser must specify information about lake depth and area. Because the lake area isallowed to vary with the lake volume, the shape of the lake basin must be specified. The
area of the lake basin at each node can be either calculated empirically as a function of
maximum depth and number of nodes or defined (along with the depth) at each node by
the user. Note: when specifying the lake basin shape empirically, the number of points
required is equal to the number of lake thermal nodes, even if the lake does not initiallyoccupy the entire basin.
Table 6.6 Lake parameters

Variable Name

Units

Number of
ValuesDescription

numnodN/A

1

Number of lake thermal nodes(also = number of lake basin
profile points)

mindepthm1Minimum allowable lake
depth

wfrac

N/A

1

Outflow channel width, as
fraction of lake perimeter(outflow is modeled as flow
over a broad-crested weir, ofwidth=wfrac*perimeter)

depth_inm1

Initial lake depth

maxdepthm1

Maximum lake depth154

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Variable Name

rpercentUnitsfraction

Number of
Values1

Description

Fraction of the grid cell
runoff routed through the lake

The next set of parameters is specified in pairs. For LAKE_PROFILE=TRUE in the
global parameter file, only one pair is needed (VIC computes the depth profile as aparabola between area 0 and the specified surface area). Otherwise, numnod pairs mustbe specified, giving the depth and area of the lake basin, one for each node:

Variable Name

Units

Number of
ValuesDescription

depth

m

1 (per pair)

Depth of lake basin at eachnode

surfacem2

1 (per pair)

Area of lake basin at eachnode

6.4.5 Meteorological and Radiative Forcings

The VIC model is forced with observed surface meteorological data which include
precipitation, temperature, wind, vapor pressure, incoming longwave and shortwave
radiation, and air pressure. The forcings data are over the land areas of the globe, at 3hourly, 1 degree resolution for 1948-2006. The forcing algorithms are explained indetails in the forcing ATBD.

6.5

Calibration

Like most physically based hydrologic models, the VIC model has many parameters thatmust be specified (about 20, depending on how the term “parameter” is defined).
However, most of the parameters can be derived from in situ measurement and remotesensing observation. The usual implementation approach (see e.g. Nijssen et al. 1997)
involves calibration of six parameters: a) the infiltration parameter (b i), which controls
the partitioning of rainfall (or snowmelt) into infiltration and direct runoff (a higher valueof b i gives lower infiltration and yields higher surface runoff); b) D2 and D3, which are
the second and third soil layer thicknesses ( D1, the top soil layer depth, is usuallyspecified a priori) and affect the water available for transpiration and baseflow
respectively (thicker soil depths have slower runoff response ─ baseflow dominated─with higher evapotranspiration, but result in longer retention of soil moisture and higher
baseflow in wet seasons); c) Dsmax, Ds, and Ws, which are baseflow parameters and also
are estimated via calibration. Dsmax is the maximum baseflow velocity, Ds is the fractionof maximum baseflow velocity, and Ws is the fraction of maximum soil moisture content

155

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Modelof the third soil layer at which non-linear baseflow occurs. These three baseflow
parameters determine how quickly the water stored in the third soil layer is evacuated asbaseflow (Liang et al. 1994). The three baseflow parameters and the third soil layer depth(d3) (Nijssen et al., 2001a, Su et al., 2005) are used with only minor adjustment during
the calibration, while the infiltration parameter (bi) and the second soil depth ( d2) aretargeted for intensive calibration. Parameters bi and d 2 are calibrated independently.
Here are the general guidelines to VIC model calibration:
1) Ds - [>0 to 1] This is the fraction of Dsmax where non-linear (rapidly increasing)baseflow begins. With a higher value of Ds, the baseflow will be higher at lower
water content in the lowest soil layer.
2) Dsmax - [>0 to ~30, depends on hydraulic conductivity] This is the maximum
baseflow that can occur from the lowest soil layer (in mm/day).
3) Ws - [>0 to 1] This is the fraction of the maximum soil moisture (of the lowest
soil layer) where non-linear baseflow occurs. This is analogous to Ds. A higher
value of Ws will raise the water content required for rapidly increasing, non-linear
baseflow, which will tend to delay runoff peaks.
4) bi - [>0 to ~0.4] This parameter defines the shape of the Variable Infiltration
Capacity curve. It describes the amount of available infiltration capacity as afunction of relative saturated gridcell area. A higher value of b i gives lower
infiltration and yields higher surface runoff.
5) Soil Depth (of each layer) - [typically 0.1 to 1.5 meters] Soil depth effects manymodel variables. In general, for runoff considerations, thicker soil depths slow
down (baseflow dominated) seasonal peak flows and increase the loss due toevapotranspiration. The maximum soil moisture storage capacity is dynamicallydetermined by the change of soil thickness. The thicker the soil depths are(resulting in more soil moisture stored in the soil layers), the less runoff isgenerated.

The calibration of these parameters is conducted via a trial and error procedure that leads
to an acceptable match of model-predicted discharge with observations. Besides visualcomparison of monthly simulated and observed hydrographs, two objective functions areoften used. One is the Nash-Sutcliffe efficiency (Ef) which describes the prediction skillof the modeled streamflow as compared to the observed value. The other is the relative
error (Er) between simulated and observed mean annual runoff. Ef and Er are calculated
as:
N

 (Q
mod, i Qobs ,i ) 2

E f  1  i N1

 (Qobs ,i Qobs ) 2

i 1

156

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) ModelE r  (Qmod  Qobs ) / Qobs
where Qmod,i is the monthly modeled streamflow for month i, Qobs is the monthly
observed streamflow for month i, N is the number of months, and Qmod and Qobs are themean of the monthly modeled and observed streamflows, respectively. When Ef equals
1.0 it means that the model perfectly predicts the observations. In addition, Shi et al.
(2008) has successfully coupled the VIC model with the Multi Objective COMplex
evolution (MOCOM-UA) algorithm of Yapo et al. (1998) and implemented thisautomated parameter estimation algorithm for eight river basins across the western
United States.
Since its existence, VIC has been well calibrated in a number of large river basins overthe continental US and the globe (Abdulla et al., 1996; Bowling et al., 2000; Crow et al.,
2003; Lohmann et al., 1998b; Maurer et al. 2001; Nijssen et al. 1997, 2001a; Su et al.
2005; Troy et al., 2008; Wood et al. 1997; Zhu and Lettenmaier, 2007). The procedure ofmatching the simulated and observed streamflows through calibration ensures that
evapotranspiration is realistically estimated over a sufficiently long enough time. This isbecause the change in surface storage is relatively small compared to other accumulated
variables in the water balance system. On this basis, and given the physically basedmodel parameterizations of the soil moisture and energy fluxes calculation, the othersurface fluxes and state variables such as soil moisture should represent observations
reasonably well, at least in the aggregate. Nijssen et al. (2001a) selected 9 basins acrossthe globe to calibrate, with each basin representing a unique climate zone (and there is atleast one basin on each continent). The calibration was successful in the arctic andtemperate climate zones, although it did not work as well in the tropical climate zone. Su
et al. (2005) partitioned the pan-Arctic drainage basin system into 12 regions andcalibrated 9 of them. The VIC model did a good job of reproducing observed streamflow
in the coldest areas of the domain that were mostly underlain by permafrost, while
problems remained in areas of discontinuous permafrost where simulated streamflow was
mostly overestimated. Zhu and Lettenmaier (2007) selected 14 comparatively small
basins (less than 10 000 km 2) over the whole Mexico to calibrate (with the basinsrepresenting different climate zones). Generally, VIC did a good job of capturing thepeak time and temporal pattern of streamflow for both arid and wet regions. The obvious
problem is the great overestimation or underestimation of peak flows in some years
especially for arid basins.

6.6

Validation and Applications

The VIC model has been validated at large scales by participating in large projects suchas PIPLS and NLDAS. For PILPS, VIC water and energy fluxes were evaluated together
with 15 other models over the Arkansas-Red river basin (Liang et al., 1998; Lohmann etal., 1998c; Wood et al., 1998), in tropical forests (Pitman et al., 1999), and in the Tomeand Kalix river systems in the cold region (Bowling et al, 2003a; 2003b; Nijssen et al.,2003). For the NLDAS, the validated variables include streamflow and water budget(Lohmann et al., 2004), soil moisture and surface temperature (Robock et al., 2003;

157

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Schaake et al., 2004), energy budget (Robock et al., 2003), and snow cover and snowpackcontent (Sheffield et al., 2003; Pan et al., 2003).The NLDAS validation results by
Mitchell et al. (2004) suggest that VIC performed well as compared to other land surface
schemes. VIC has also been evaluated using soil moisture observations in the U.S.
(Maurer et al, 2002) and global snow cover extent data (Nijssen et al, 2001b).
Throughout its existence, the VIC model has been used in many research areas, such as
meteorology and atmospheric sciences, water resources, geosciences, environmental
sciences, remote sensing, etc. So far VIC has been cited about 1700 times from all over
the world. Figure 6.7 shows the number of citations in each year, suggesting a continuous
positive trend of VIC popularity. Table 6.7 lists a selected number of papers directly
using the VIC model or VIC model results.

Figure 6.7 VIC model citations in each year

Table 6.7 List of research papers using VIC model directly

Authors

Studied region

Specific topic
Reservoir influence on
streamflow

Adam et. (2007)

Eurasian arctic

Andreadis andLettenmaier (2006a)

United States

Arora and Boer (2006)

GlobalArora (2001)
Berbery et al (2003)
Cherkauer andLettenmaier (2003)
Cosgrove et al. (2003)
Crow et al. (2003)
Demaria et al. (2007)
Feng et al. (2008)

Russian lowlands
Mississippi basin

Temporal variability of soil
moisture
Water balance
Hydrological cycle

Minnesota River

Snow and frozen soil

North America
Southern Great Plains
Selected US basins
Colorado

Model spin-up behavior
Model calibration
Model parameter sensitivity
Snow simulation

Drought analysis

158

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Gao et al. (2006)

North America

Gao et al. (2007)

Southern Great Plains

Gao et al. (2004)

Southern Great Plains

Huang et al. (2003)
Li et al. (2007)

Selected US basins
Rio Grande

Liang et al. (2003)

Pennsylvania

Liang et al. (1994)
Lucas-Picher et al.
(2003)
Luo et al. (2005)
Luo and Wood (2007)
Maurer et al. (2001)
Meng and Quiring (2008)
Miguez-Macho et al.
(2008)Mo (2008)

Kansas

Surface soil moistureremote sensing
Data assimilation
Surface soil moistureremote sensing
Model parameter transfer
Water cycle
Surface and ground waterinteraction
VIC model

North America

Routing scheme validation

Columbia and Colorado basin
United States
Mississippi river basin
Texas and Maryland

Nijssen et al. (2001a)

Global basins

Nijssen et al. (2001b)
O’Donnell et al. (2000)
Pan and Wood (2006)
Pan et al. (2003)
Pan et al. (2008)
Parada and Liang (2008)

GlobalOhio river basin
Southern Great Plains
North America
Arkansas-Red river basin
Southern Great Plains

Rhoads et al. (2001)

Arkansas-Red river basin

Water cycle
Drought analysis
Water budget evaluation
Soil moisture modeling
Water table and Soil
moisture modeling
Drought indicesHydrologic sensitivity ofriversSoil moisture modeling
Water and energy balance
Data assimilation
Snow modeling
Data assimilation
Data assimilation
Surface temperaturemodeling

Sheffield and Wood
(2007)Sheffield et al. (2004a)
Sheffield et al. (2004b)
Sheffield et al. (2003)
Stamm et al. (1994)
Su et al. (2006)
Su et al. (2005)
Wang et al. (2008)

North America
United States

GlobalDrought analysis

United States
GlobalNorth America
Globalpan-Arctic
ArcticDrought analysis
Precipitation correction
Snow modeling
Global climate sensitivity
Surface water flux
Streamflow
Integration to climate
model

Global159

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Warrach-Sagi et al.
(2008)
Wojcik et al. (2008)
Wood et al. (2002)
Wood et al. (1992)
Wu et al. (2007)
Xie et al. (2007)
Zhou et al. (2006)
Zhu and Lettenmair(2007)
Abdulla and Lettenmaier(1997a)Abdulla and Lettenmaier(1997b)Andreadis andLettenmaier (2006b)

Germany

Streamflow

Colorado
United States
North Carolina
China
China
Baohe river, China

Snow modelingHydrologic forcasting
VIC model
Soil moisture modeling
Parameter estimation
Hydrological cycle
Surface hydrology andenergy flux

MexicoArkansas-Red river basin

Parameter estimation

Arkansas-Red river basin

Water balance

Snake river basin

Data assimilation

Guo et al. (2004)

Illinois river watershed

Huang and Liang (2006)
Hurkmans et al. (2008)
Lakshmi and Wood
(1998)
Liang et al. (2004)
Liang and Xie (2003)
Lobmeyr et al. (1999)
Lohmann et al. (1998b)
Matheussen et al. (2000)
Mengelkamp et al.
(1999)
Nijssen et al. (1997)

Selected watersheds in US
Rhine river basin

Precipitation impact onwater budget
Parameter estimation
River discharge

King’s creak catchment

evaporation

Blue river watershed
Pennsylvania
Elbe river, Germany
West river, Germany
Columbia river basin

Water flux
RunoffWater balanceStreamflow
Streamflow
Surface energy and waterbalance
Streamflow

Shaman et al. (2002)

Netherland, Arkansas
Major rivers in US
Sleepers river and Black Rockcatchment

Storm flow

Silberstein et al. (2002)

Western Australia

Sivapalan et al. (1997)
Sivapalan and Woods(1995)Su and Xie (2003)
te Linde et al. (2008)

Australia

Water logging and ground
water
RunoffAustralia

Water balance fluxes

China
Rhine river basin

RunoffDischarge

160

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Troy et al. (2008)
VanShaar et al. (2002)
Wooldridge SA, Kalma(2001)
Yang and Xie (2003)
Yuan et al. (2004)
Liang et al. (1996)
Slater et al. (2007)
Han and Li (2008)
Hillard et al. (2003)

United states
Columbia river basin

Parameter estimation
RunoffEastern Australia

Streamflow

China
Hanjiang river basin (China)
PILPS sites
pan-Arctic
California
Upper Mississippi river basin,
central Canada

Groundwater table
Streamflow
VIC model
Hydrologic processes
Data assimilation

161

Snow modelingChapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

6.7

References

Abdulla, F. A., et al. (1996), Application of a macroscale hydrologic model to estimate
the water balance of the Arkansas Red River basin, J. Geophys. Res.-Atmos., 101(D3),7449-7459.
Abdulla, F. A., and D. P. Lettenmaier (1997a), Development of regional parameter
estimation equations for a macroscale hydrologic model, Journal of Hydrology, 197(1-4),
230-257.
Abdulla, F. A., and D. P. Lettenmaier (1997b), Application of regional parameter
estimation schemes to simulate the water balance of a large continental river, Journal of
Hydrology, 197(1-4), 258-285.
Adam, J. C. (2007), Understanding the Causes of Streamflow Changes in the Eurasian
Arctic, 156 pp, University of Washington
Adam, J. C., et al. (2007), Simulation of reservoir influences on annual and seasonal
streamflow changes for the Lena, Yenisei, and Ob' rivers, J. Geophys. Res.-Atmos.,
112(D24), -.
Anderson, E. A. (1976), A Point Energy and Mass Balance Model of a Snow Cover,
NOAA Technical Report.
Andreadis, K., et al. (2009), Modeling snow accumulation and ablation processes inforested environments, Water Resour Res, in review.
Andreadis, K. M., and D. P. Lettenmaier (2006a), Trends in 20th century drought overthe continental United States, Geophys Res Lett, 33(10).
Andreadis, K. M., and D. P. Lettenmaier (2006b), Assimilating remotely sensed snow
observations into a macroscale hydrology model, Advances in Water Resources , 29(6),872-886.
Arora, V. K. (2001), Assessment of simulated water balance for continental-scale river
basins in an AMIP 2 simulation, J. Geophys. Res.-Atmos., 106(D14), 14827-14842.
Arora, V. K., and G. J. Boer (2006), The temporal variability of soil moisture and surfacehydrological quantities in a climate model, J. Clim., 19(22), 5875-5888.
Berbery, E. H., et al. (2003), Eta model estimated land surface processes and thehydrologic cycle of the Mississippi basin, J. Geophys. Res.-Atmos., 108(D22).
Blondin, C. (1991), Parameterization of land-surface processes in numerical weather
prediction, in Land Surface Evaporation: Measurements and Parameterization , edited byT. J. Schmugge and J. C. Andre, pp. 31-54, Springer-Verlag, New York.

162

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Bohn, T. J., et al. (2007), Methane emissions from western Siberian wetlands:
heterogeneity and sensitivity to climate change, Environ. Res. Lett., 2(4), 9.
Bowling, L., et al. (2008), Current capabilities in soil thermal representations within alarge scale hydrology model, in 9th International Conference on Permafrost , edited,
Fairbanks, AK.
Bowling, L. C., et al. (2000), Hydrologic effects of logging in western Washington,
United States, Water Resour Res, 36(11), 3223-3240.
Bowling, L. C., et al. (2003a), Simulation of high-latitude hydrological processes in theTorne-Kalix basin: PILPS phase 2(e) - 1: Experiment description and summaryintercomparisons, Global Planet Change, 38(1-2), 1-30.
Bowling, L. C., et al. (2003b), Simulation of high-latitude hydrological processes in theTorne-Kalix basin: PILPS phase 2(e) - 3: Equivalent model representation and sensitivity
experiments, Global Planet Change, 38(1-2), 55-71.
Bowling, L. C., et al. (2003c), The role of surface storage in a low-gradient Arctic
watershed, Water Resour Res, 39(4), -.
Bowling, L. C., et al. (2004), Parameterization of blowing-snow sublimation in amacroscale hydrology model, Amer Meteorological Soc.
Bowling, L. C., and D. P. Lettenmaier (2009), Modeling the effects of lakes and wetlandson the water balance of Arctic environments, Journal of Hydrometeorology , submitted.
Bras, R. A. (1990), Hydrology, an Introduction to Hydrologic Science , 643 pp., AddisonWesley.
Brooks, R. H., A. H. Corey (1988), Hydraulic properties of porous media, Hydrol. Pap.,
Colorado State University, 3.
Calder, I. R. (1993), Hydrological effects of land-use change. Handbook of Hydrology,
D.R. Maidment, Ed., McGraw-Hill, 3.1-3.50.
Cherkauer, K. A., and D. P. Lettenmaier (1999), Hydrologic effects of frozen soils in theupper Mississippi River basin, J. Geophys. Res.-Atmos., 104(D16), 19599-19610.
Cherkauer, K. A., et al. (2003), Variable infiltration capacity cold land process model
updates, Global Planet Change, 38(1-2), 151-159.
Cherkauer, K. A., and D. P. Lettenmaier (2003), Simulation of spatial variability in snow
and frozen soil, J. Geophys. Res.-Atmos., 108(D22), 14.
Christensen, N. S., and D. P. Lettenmaier (2007), A multimodel ensemble approach toassessment of climate change impacts on the hydrology and water resources of theColorado River Basin, Hydrology and Earth System Sciences , 11(4), 1417-1434.

163

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Cosby, B. J., et al. (1984), A statistical exploration of the relationships of soil moisturecharacteristics to the physical properties of soils, Water Resour Res, 20(6), 682-690.
Cosgrove, B. A., et al. (2003), Land surface model spin-up behavior in the North
American Land Data Assimilation System (NLDAS), J. Geophys. Res.-Atmos., 108(D22),-.
Crow, W. T., et al. (2003), Multiobjective calibration of land surface model
evapotranspiration predictions using streamflow observations and spaceborne surface
radiometric temperature retrievals, J. Geophys. Res.-Atmos., 108(D23), -.
Deardorff, J. W. (1978), Efficient prediction of ground surface-temperature and moisture,with inclusion of a layer of vegetation, Journal of Geophysical Research-Oceans andAtmospheres, 83(NC4), 1889-1903.
Defries, R. S., and J. R. G. Townshend (1994), NDVI-derived land cover classification at
global scales, International Journal of Remote Sensing , 15(17), 3567-3586.
Demaria, E. M., et al. (2007), Monte Carlo sensitivity analysis of land surface parameters
using the Variable Infiltration Capacity model, J. Geophys. Res.-Atmos., 112(D11).
Dickinson, R. E. (1984), Modeling evapotranspiration for three-dimentional global
climate models, in Climate Processes and Climate Sensitivity, Monogr. Ser., edited by J.
E. Hansen and T. Takahashi, pp. 58-72, Washington, D.C.
Dickinson, R. E., A. Henderson-Sellers, P. J. Kennedy, M. F. Wilson (1986), Biosphereatmosphertera nsfer scheme( BATS) for the NCAR community climate model, NCAR
Tech. Note TN-275 +STR.
Duband, D., et al. (1993), Unit-hydrograph revised - an alternate iterative approach to UH
and effective precipitation identification, Journal of Hydrology, 150(1), 115-149.
Ducoudre, N. I., et al. (1993a), SECHIBA, a new set of parameterizations of thehydrologic exchanges at the land-atmosphere interface within the LMD atmosphericgeneral circulation model, J. Clim., 6(2), 248-273.
Ducoudre, N. I., et al. (1993b), Sechiba, a New Set of Parameterizations of theHydrologic Exchanges at the Land Atmosphere Interface within the Lmd Atmospheric
General-Circulation Model, J. Clim., 6(2), 248-273.
FAO (2003), AQUASTAT: FAO’s Information System on Water and Agriculture. Foodand Agriculture Organization of the United Nations, Rome, Italy. .
Farouki, O. T. (1986), Thermal Properties of Soils, in Rock and Soil Mechanics, Vol 11,
pp. 102-119, Trans-Tech Publications, Claushal-Zellerfeld, Germany.
Feng, X., et al. (2008), The Impact of Snow Model Complexity at Three CLPX Sites,
Journal of Hydrometeorology , 9(6), 1464-1481.

164

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Flerchinger, G. N., and K. E. Saxton (1989), SIMULTANEOUS HEAT AND WATERMODEL OF A FREEZING SNOW-RESIDUE-SOIL SYSTEM .1. THEORY AND
DEVELOPMENT, Transactions of the Asae, 32(2), 565-571.
Franchini, M., and M. Pacciani (1991), Comparative-analysis of several conceptual
rainfall runoff models, Journal of Hydrology, 122(1-4), 161-219.
Gao, H., et al. (2006), Using TRMM/TMI to retrieve surface soil moisture over thesouthern United States from 1998 to 2002, Journal of Hydrometeorology , 7(1), 23-38.
Gao, H. L., et al. (2004), Using a microwave emission model to estimate soil moisturefrom ESTAR observations during SGP99, Journal of Hydrometeorology , 5(1), 49-63.
Gao, H. L., et al. (2007), Copula-derived observation operators for assimilating TMI andAMSR-E retrieved soil moisture into land surface models, Journal of Hydrometeorology ,8(3), 413-429.
Guo, J. Z., et al. (2004), Impacts of different precipitation data sources on water budgets,
Journal of Hydrology, 298(1-4), 311-334.
Haddeland, I., et al. (2006a), Effects of irrigation on the water and energy balances of theColorado and Mekong river basins, Journal of Hydrology, 324(1-4), 210-223.
Haddeland, I., et al. (2006b), Anthropogenic impacts on continental surface water fluxes,
Geophys Res Lett, 33(8), -.
Haddeland, I., et al. (2007), Hydrologic effects of land and water management in North
America and Asia: 1700-1992, Hydrology and Earth System Sciences , 11(2), 1035-1045.
Hamill, L. (2001), Understanding Hydraulics , 2 ed., Palgrave.
Hamlet, A. F., and D. P. Lettenmaier (1999), Columbia River streamflow forecasting
based on ENSO and PDO climate signals, J. Water Resour. Plan. Manage.-ASCE, 125(6),333-341.
Han, X., and X. Li (2008), An evaluation of the nonlinear/non-Gaussian filters for thesequential data assimilation, Remote Sensing of Environment, 112(4), 1434-1449.
Hanasaki, N., et al. (2006), A reservoir operation scheme for global river routing models,
Journal of Hydrology, 327(1-2), 22-41.
Hansen, M. C., et al. (2000), Global land cover classification at 1km spatial resolution
using a classification tree approach, International Journal of Remote Sensing , 21(6-7),
1331-1364.
Hillard, U., et al. (2003), Assessing snowmelt dynamics with NASA scatterometer(NSCAT) data and a hydrologic process model, Remote Sensing of Environment, 86(1),52-69.

165

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Hostetler, S. W., and P. J. Bartlein (1990), SIMULATION OF LAKE EVAPORATIONWITH APPLICATION TO MODELING LAKE LEVEL VARIATIONS OF HARNEYMALHEUR LAKE, OREGON, Water Resour Res, 26(10), 2603-2612.
Hostetler, S. W. (1991), Simulation of lake ice and its effect on the late-Pleistocene
evaporation rate of Lake Lahontan Climate Dynamics, 6, 43-48.
Huang, M., et al. (2003), A transferability study of model parameters for the variableinfiltration capacity land surface scheme, J. Geophys. Res.-Atmos., 108(D22).
Huang, M. Y., and X. Liang (2006), On the assessment of the impact of reducing
parameters and identification of parameter uncertainties for a hydrologic model withapplications to ungauged basins.
Hurkmans, R., et al. (2008), Water balance versus land surface model in the simulation ofRhine river discharges, Water Resour Res, 44(1).
ICOLD (2003), World Register of Dams 2003. International Commission on Large Dams
(ICOLD), Paris, France.
Jackson, R. B., et al. (1996), A global analysis of root distributions for terrestrial biomes,
Oecologia, 108(3), 389-411.
Jordan, R. (1991), A one-dimensional temperature model for a snow cover: Technicaldocumentation for sntherm89, US Army Corps of Engineers Cold Regions Research andEngineering Laboratory, Hanover, NH.
.
Kimball, J. S., et al. (1997), An improved method for estimating surface humidity fromdaily minimum temperature, Agric. For. Meteorol., 85(1-2), 87-98.
Lakshmi, V., and E. F. Wood (1998), Diurnal cycles of evaporation using a two-layerhydrological model, Journal of Hydrology, 204(1-4), 37-51.
Lettenmaier, D. P., and F. Su (2009), in ARCTIC Climate Change-The ACSYS Decadeand Beyond, edited, in press.
Li, J., et al. (2007), Modeling and analysis of the variability of the water cycle in theupper Rio Grande basin at high resolution, Journal of Hydrometeorology , 8(4), 805-824.
Liang, X., et al. (1994), A SIMPLE HYDROLOGICALLY BASED MODEL OF LANDSURFACE WATER AND ENERGY FLUXES FOR GENERAL-CIRCULATIONMODELS, J. Geophys. Res.-Atmos., 99(D7), 14415-14428.
Liang, X., et al. (1996), Surface soil moisture parameterization of the VIC-2L model:
Evaluation and modification, Global Planet Change, 13(1-4), 195-206.

166

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Liang, X., and Z. H. Xie (2003), Important factors in land-atmosphere interactions:surface runoff generations and interactions between surface and groundwater, Global
Planet Change, 38(1-2), 101-114.
Liang, X., et al. (2003), A new parameterization for surface and groundwater interactions
and its impact on water budgets with the variable infiltration capacity (VIC) land surface
model, J. Geophys. Res.-Atmos., 108(D16).
Liang, X., et al. (2004), Assessment of the effects of spatial resolutions on daily water
flux simulations, Journal of Hydrology, 298(1-4), 287-310.
Lobmeyr, M., et al. (1999), An application of a large scale conceptual hydrological modelover the Elbe region, Hydrology and Earth System Sciences , 3(3), 363-374.
Lohmann, D., et al. (1996), A large scale horizontal routing model to be coupled to land
surface parameterization schemes, Tellus(48A), 708-721.
Lohmann, D., et al. (1998a), Regional scale hydrology: I. Formulation of the VIC-2Lmodel coupled to a routing model, Hydrol. Sci. J.-J. Sci. Hydrol., 43(1), 131-141.
Lohmann, D., et al. (1998b), Regional scale hydrology: II. Application of the VIC-2Lmodel to the Weser River, Germany, Hydrol. Sci. J.-J. Sci. Hydrol., 43(1), 143-158.
Lohmann, D., et al. (2004), Streamflow and water balance intercomparisons of four land
surface models in the North American Land Data Assimilation System project, J.
Geophys. Res.-Atmos., 109(D7), 22.
Louis, J. F. (1979), Parametric Model of Vertical Eddy Fluxes in the Atmosphere, BoundLay Meteorol, 17(2), 187-202.
Lucas-Picher, P., et al. (2003), Implementation of a large-scale variable velocity river
flow routing algorithm in the Canadian Regional Climate Model (CRCM), AtmosphereOcean, 41(2), 139-153.
Luo, L. F., and E. F. Wood (2007), Monitoring and predicting the 2007 U.S. drought,
Geophys Res Lett, 34(22), -.
Luo, Y., et al. (2005), The operational eta model precipitation and surface hydrologic
cycle of the Columbia and Colorado basins, Journal of Hydrometeorology , 6(4), 341-370.
Mahrt, L., and H. Pan (1984), A 2-Layer Model of Soil Hydrology, Bound-Lay Meteorol,
29(1), 1-20.
Matheussen, B., et al. (2000), Effects of land cover change on streamflow in the interior
Columbia River Basin (USA and Canada), Hydrol. Process., 14(5), 867-885.

167

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Maurer, E. P., et al. (2001), Evaluation of the land surface water budget in NCEP/NCAR
and NCEP/DOE reanalyses using an off-line hydrologic model, J. Geophys. Res.-Atmos.,
106(D16), 17841-17862.
Maurer, E. P., et al. (2002), A long-term hydrologically based dataset of land surface
fluxes and states for the conterminous United States, J. Clim., 15(22), 3237-3251.
Meng, L., and S. M. Quiring (2008), A comparison of soil moisture models using Soil
Climate Analysis Network observations, Journal of Hydrometeorology , 9(4), 641-659.
Mengelkamp, H. T., et al. (1999), SEWAB - a parameterization of the Surface Energy
and Water Balance for atmospheric and hydrologic models, Advances in Water Resources ,23(2), 165-175.
Miguez-Macho, G., et al. (2008), Simulated water table and soil moisture climatology
over North America, Bulletin of the American Meteorological Society, 89(5), 663-+.
Mo, K. C. (2008), Model-Based Drought Indices over the United States, Journal of
Hydrometeorology, 9(6), 1212-1230.
Monteith, J. L., and M. H. Unsworth (1990), Principles of Environmental Physics, 2nded., Routledge, Chapman and Hall, New York.
Mosley, M. P., and A. I. McKerchar (1993), Streamflow, in Handbook ofHydrology, edited by D. R. Maidment, p. 8.3, McGraw-Hill, New York.
Myneni, R. B., et al. (1997), Estimation of global leaf area index and absorbed par using
radiative transfer models, Ieee Transactions on Geoscience and Remote Sensing , 35(6),1380-1393.
Nijssen, B., et al. (1997), Streamflow simulation for continental-scale river basins, Water
Resour Res, 33(4), 711-724.
Nijssen, B., et al. (2001a), Predicting the discharge of global rivers, J. Clim., 14(15),3307-3323.
Nijssen, B., et al. (2001b), Global retrospective estimation of soil moisture using thevariable infiltration capacity land surface model, 1980-93, J. Clim., 14(8), 1790-1808.
Nijssen, B., et al. (2003), Simulation of high latitude hydrological processes in the TorneKalix basin: PILPS phase 2(e) - 2: Comparison of model results with observations,Global Planet Change, 38(1-2), 31-53.
O'Donnell, G. M., et al. (2000), Macroscale hydrological modeling using remotely sensed
inputs: Application to the Ohio River basin, J. Geophys. Res.-Atmos., 105(D10), 1249912516.

168

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Pan, M., et al. (2003), Snow process modeling in the North American Land Data
Assimilation System (NLDAS): 2. Evaluation of model simulated snow water equivalent,
J. Geophys. Res.-Atmos., 108(D22), 14.
Pan, M., and E. F. Wood (2006), Data assimilation for estimating the terrestrial water
budget using a constrained ensemble Kalman filter, Journal of Hydrometeorology , 7(3),534-547.
Pan, M., et al. (2008), Estimation of regional terrestrial water cycle using multi-sensor
remote sensing observations and data assimilation, Remote Sensing of Environment ,112(4), 1282-1294.
Parada, L. M., and X. Liang (2008), Impacts of spatial resolutions and data quality on soil
moisture data assimilation, J. Geophys. Res.-Atmos., 113(D10).
Patterson, J. C., and P. F. Hamblin (1988), THERMAL SIMULATION OF A LAKE
WITH WINTER ICE COVER, Limnology and Oceanography , 33(3), 323-338.
Rawls, W. J., et al. (1998), Use of soil texture, bulk density, and slope of the water
retention curve to predict saturated hydraulic conductivity, Transactions of the Asae ,41(4), 983-988.
Reynolds, C. A., et al. (2000), Estimating soil water-holding capacities by linking theFood and Agriculture Organization soil map of the world with global pedon databases
and continuous pedotransfer functions, Water Resour Res, 36(12), 3653-3662.
Rhoads, J., et al. (2001), Validation of land surface models using satellite-derived surface
temperature, J. Geophys. Res.-Atmos., 106(D17), 20085-20099.
Shaman, J., et al. (2002), Representation of subsurface storm flow and a more responsive
water table in a TOPMODEL-based hydrology model, Water Resour Res, 38(8).
Sheffield, J., et al. (2003), Snow process modeling in the North American Land Data
Assimilation System (NLDAS): 1. Evaluation of model-simulated snow cover extent, J.
Geophys. Res.-Atmos., 108(D22), 13.
Sheffield, J., et al. (2004a), A simulated soil moisture based drought analysis for theUnited States, J. Geophys. Res.-Atmos., 109(D24), 19.
Sheffield, J., et al. (2004b), Correction of the high-latitude rain day anomaly in the
NCEP-NCAR reanalysis for land surface hydrological modeling, J. Clim., 17(19), 38143828.
Sheffield, J., and E. F. Wood (2007), Characteristics of global and regional drought,
1950-2000: Analysis of soil moisture data from off-line simulation of the terrestrialhydrologic cycle, J. Geophys. Res.-Atmos., 112(D17), -.

169

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Shi, X. G., et al. (2008), How Essential is Hydrologic Model Calibration to Seasonal
Streamflow Forecasting?, Journal of Hydrometeorology , 9(6), 1350-1363.
Shuttleworth, W. J. (1993), Evaporation, in Handbook of Hydrology , edited by D. R.
Maidment, pp. 4.1-4.53, McGraw-Hill, Inc., New York.
Siebert, S., et al. (2002), Global Map of Irrigated Areas Version 2.1. Center for
Environmental Systems Research, University of Kassel, Germany/Food and Agriculture
Organization of the United Nations, Rome, Italy.
Silberstein, R. P., et al. (2002), Mechanisms and control of water logging andgroundwater flow in the 'Ucarro' sub-catchment.
Sivapalan, M., and R. A. Woods (1995), EVALUATION OF THE EFFECTS OF
GENERAL-CIRCULATIONMODELSSUBGRID
VARIABILITYAND
PATCHINESS OF RAINFALL AND SOIL-MOISTURE ON LAND-SURFACEWATER-BALANCE FLUXES.
Sivapalan, M., et al. (1997), Variable bucket representation of TOPMODEL andinvestigation of the effects of rainfall heterogeneity, Hydrol. Process., 11(9), 1307-1330.
Slater, A. G., et al. (2007), A multimodel simulation of pan-Arctic hydrology, J Geophys
Res-Biogeo, 112(G4), -.
Stamm, J. F., et al. (1994), SENSITIVITY OF A GCM SIMULATION OF GLOBAL
CLIMATE TO THE REPRESENTATION OF LAND-SURFACE HYDROLOGY, J.
Clim., 7(8), 1218-1239.
Storck, P., et al. (1998), Application of a GIS-based distributed hydrology model for
prediction of forest harvest effects on peak stream flow in the Pacific Northwest, Hydrol.
Process., 12(6), 889-904.
Storck, P., et al. (2002), Measurement of snow interception and canopy effects on snow
accumulation and melt in a mountainous maritime climate, Oregon, United States, Water
Resour Res, 38(11), 16.
Su, F. G., and Z. H. Xie (2003), A model for assessing effects of climate change on
runoff in China, Progress in Natural Science , 13(9), 701-707.
Su, F. G., et al. (2005), Streamflow simulations of the terrestrial Arctic domain, J.
Geophys. Res.-Atmos., 110(D8), -.
Su, F. G., et al. (2006), Evaluation of surface water fluxes of the pan-Arctic land region
with a land surface model and ERA-40 reanalysis, J. Geophys. Res.-Atmos., 111(D5), -.
Takeuchi, K. (1997), Least marginal environmental impact rule for reservoir
development, Hydrol. Sci. J.-J. Sci. Hydrol., 42(4), 583-597.

170

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Tarboton, D. G., et al. (1995), A spatially distributed energy balance snowmelt model, in
Biogeochemistry of Seasonally Snow Covered Catchments , edited by K. A. Tonneson, etal., pp. 141–155, International Association of Hydrological Sciences.
te Linde, A. H., et al. (2008), Comparing model performance of two rainfall-runoff
models in the Rhine basin using different atmospheric forcing data sets, Hydrology and
Earth System Sciences, 12(3), 943-957.
Thornton, P. E., and S. W. Running (1999), An improved algorithm for estimating
incident daily solar radiation from measurements of temperature, humidity, andprecipitation, Agric. For. Meteorol., 93(4), 211-228.
Troy, T. J., et al. (2008), An efficient calibration method for continental-scale land
surface modeling, Water Resour Res, 44(9), 13.
VanShaar, J. R., et al. (2002), Effects of land-cover changes on the hydrological responseof interior Columbia River basin forested catchments, Hydrol. Process., 16(13), 24992520.
Vrugt, J. A., et al. (2003), A Shuffled Complex Evolution Metropolis algorithm for
optimization and uncertainty assessment of hydrologic model parameters, Water Resour
Res, 39(8).
Wang, A., et al. (2008), Integration of the variable infiltration capacity model soilhydrology scheme into the community land model, J. Geophys. Res.-Atmos., 113(D9), -.
Wang, A., et al. (2009), Multimodel ensemble reconstruction of drought over thecontinental United States, J. Clim., accepted.
Warrach-Sagi, K., et al. (2008), Streamflow simulations reveal the impact of the soil
parameterization, Meteorologische Zeitschrift, 17(6), 751-762.
Wigmosta, M. S., et al. (1994), A DISTRIBUTED HYDROLOGY-VEGETATIONMODEL FOR COMPLEX TERRAIN, Water Resour Res, 30(6), 1665-1679.
Wojcik, R., et al. (2008), Multimodel Estimation of Snow Microwave Emission during
CLPX 2003 Using Operational Parameterization of Microphysical Snow Characteristics,
Journal of Hydrometeorology , 9(6), 1491-1505.
Wood, A. W., et al. (2002), Long-range experimental hydrologic forecasting for theeastern United States, J. Geophys. Res.-Atmos., 107(D20), 15.
Wood, A. W., et al. (2005), A retrospective assessment of National Centers forEnvironmental Prediction climate model-based ensemble hydrologic forecasting in thewestern United States, J. Geophys. Res.-Atmos., 110(D4), 16.

171

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Wood, A. W., and D. P. Lettenmaier (2006), A test bed for new seasonal hydrologic
forecasting approaches in the western United States, Bulletin of the American
Meteorological Society, 87(12), 1699-+.
Wood, E. F., et al. (1992), A Land-Surface Hydrology Parameterization with Subgrid
Variability for General-Circulation Models, J. Geophys. Res.-Atmos., 97(D3), 2717-2728.
Wood, E. F., et al. (1997), Hydrological modeling of continental-scale basins, Annu. Rev.
Earth Planet. Sci., 25, 279-&.
Wood, E. F., et al. (1998), The Project for Intercomparison of Land-surface
Parameterization Schemes (PILPS) phase 2(c) Red-Arkansas River basin experiment: 1.
Experiment description and summary intercomparisons, Elsevier Science Bv.
Wooldridge, S. A., and J. D. Kalma (2001), Regional-scale hydrological modelling using
multiple-parameter landscape zones and a quasi-distributed water balance model,
Hydrology and Earth System Sciences , 5(1), 59-74.
Wu, Z. Y., et al. (2007), Thirty-five year (1971-2005) simulation of daily soil moistureusing the variable infiltration capacity model over China, Atmosphere-Ocean, 45(1), 3745.
Xie, Z. H., et al. (2007), Regional parameter estimation of the VIC land surface model:
methodology and application to river basins in China, Journal of Hydrometeorology , 8(3),447-468.
Yang, H. W., and Z. H. Xie (2003), A new method to dynamically simulate groundwatertable in land surface model VIC, Progress in Natural Science , 13(11), 819-825.
Yapo, P. O., et al. (1998), Multi-objective global optimization for hydrologic models,
Journal of Hydrology, 204(1-4), 83-97.
Yuan, F., et al. (2004), An application of the VIC-3L land surface model and remotesensing data in simulating streamflow for the Hanjiang River basin, Canadian Journal of
Remote Sensing, 30(5), 680-690.
Zhao, R.-J., et al. (1980), The Xinanjiang model, Hydrological Forecasting Proceedings
Oxford Symposium, IASH 129, 351-356.
Zhou, S. Q., et al. (2006), Effects of heterogeneous vegetation on the surfacehydrological cycle, Advances in Atmospheric Sciences , 23(3), 391-404.
Zhu, C. M., et al. (2005), Role of antecedent land surface conditions on North American
monsoon rainfall variability, J. Clim., 18(16), 3104-3121.
Zhu, C. M., et al. (2007), Role of antecedent land surface conditions in warm season
precipitation over northwestern Mexico, J. Clim., 20(9), 1774-1791.

172

Chapter 6 Water Budget Record from Variable Infiltration Capacity (VIC) Model

Zhu, C. M., and D. P. Lettenmaier (2007), Long-term climate and derived surfacehydrology and energy flux data for Mexico: 1925-2004, J. Clim., 20(9), 1936-1946.
Zhu, C. M., et al. (2009), Evaluating the influence of antecedent soil moisture on
variability of the North American Monsoon precipitation in the coupled MM5/VIC
modeling system, The Journal of Advances in Modeling Earth Systems , submitted.
Zoltai, S. C. (1979), An outline of the wetland regions of Canada, paper presented at Proc.
workshop on Canadian wetlands.

173