MODELING THE EFFECTS OF LAKES AND WETLANDS ON THE WATER
BALANCE OF ARCTIC ENVIRONMENTS
Laura C. Bowling, Department of Agronomy, Purdue University
Dennis P. Lettenmaier, Department of Civil and Environmental Engineering, University of Washington
Submitted on: July 12, 2008
Corresponding author address:
Lilly Hall of Life Sciences
915 West State Street
West Lafayette, IN 47907-2054
Lakes, ponds and wetlands are common features in many low-gradient arctic watersheds. Storage of
snowmelt runoff in lakes and wetlands exerts a strong influence on both the interannual and interseasonal
variability of northern rivers. This influence is often not well represented in hydrology models and the land
surface schemes used in climate models. In this paper, we describe an algorithm to represent the
evaporation and storage effects of lakes and wetlands within the Variable Infiltration Capacity (VIC)
macroscale hydrology model. The model is evaluated with respect to its ability to represent water
temperatures, net radiation, ice freeze/thaw and runoff production for a variety of high-latitude locations. It
is then used to investigate the influence of surface storage on the spatial and temporal distribution of water
and energy fluxes for the Putuligayuk River watershed, on the Alaskan Arctic coastal plain. SAR imagery,
which can successfully identify the surface water extent, is used to parameterize the VIC model.
Simulations of runoff from the Putuligayuk watershed using the parameterized relationship indicate that up
to 50% of snow melt water goes into storage each year and does not contribute to streamflow. Seasonally,
the area of open water increases on the order of 100 – 220% due to the influx of snowmelt each spring.
Watershed surface storage is depleted over the summer by evaporation from open water, which exceeds
summer precipitation by 9 to 46%.
Lakes and wetlands are particularly prevalent in much of the Arctic region, where the presence of
permafrost and modest relief impedes the subsurface drainage of water (UNESCO 1978). Thousands of
shallow thaw lakes cover up to 40 percent of the tundra on the Alaskan North Slope (Jeffries et al. 1999)
and 30 to 50 percent of the central Mackenzie River delta region (Marsh and Bigras 1988). Similarly, in
the western Siberia lowlands, tundra marshes compose up to 40 percent of the land area (Moskvin 1989).
The predominance of surface water storage controls many aspects of the land surface water balance in these
regions, through a tendency to suppress runoff and increase evaporation (UNESCO 1978). The presence of
many lakes is often cited as a reason for lower interannual and seasonal variability in runoff from the lakes
region of the Mackenzie River, relative to the large Eurasian arctic rivers (UNESCO 1978, Vuglinsky
1997, Cattle 1985).
Observations in arctic Canada and Alaska indicate that annual evaporation from lakes and wetlands
exceeds annual precipitation (Roulet and Woo 1986, Marsh and Bigras 1988, Rovansek et al. 1996, and
Mendez et al. 1998). Rovansek et al. (1996) found that evaporation is the major cause of the observed drop
in pond water levels during the snow-free period for a tundra wetland complex in northern Alaska. The
wetlands therefore depend on a supply of water from the surrounding uplands that comes primarily from
snowmelt (Bigras 1990, Rovansek et al. 1996) and in some cases by subsurface drainage and ground ice
melt (Moskvin 1989, Young and Woo 2000).
Runoff from northern wetlands occurs as intensive surface flow following snowmelt (Rovansek et al.
1996, Moskvin 1989). The high ice content of poorly drained wetland soils makes them slow to thaw, so
infiltration may be delayed relative to drier soils (Roulet and Woo 1986, Woo and Xia 1996). Surface
storage provides an initial control on runoff at the beginning of snowmelt (Glenn and Woo 1997, Rovansek
et al. 1996). Lake volumes in the Mackenzie Delta, N.W.T., Canada have been observed to increase
between 90 and 330 percent above mid-summer volumes with the influx of snow melt water (Bigras 1990).
Lake water levels typically drop rapidly following spring melt (Bigras 1990, Bowling et al. 2003a). In the
wetlands of northern Alaska and western Siberia, runoff often ceases in the summer and the active layer
dries. Lake and wetland storage is gradually decreased through evaporation and becomes disconnected
from the surrounding watershed (Moskvin 1989, Rovansek et al. 1996, Bowling et al. 2003a).
Because of the high heat capacity of lakes and the availability of water for evaporation, lakes and wetlands
act as surface heat reservoirs that can influence regional climate and the spatial distribution of energy
fluxes. During the summer in northern Alaska, Mendez et al. (1998) found that latent heat accounted for 46
percent of net radiation over wetlands, in contrast to 23 percent in a dryer upland tundra site. Jacobs and
Grondin (1988) found that strongly positive sensible heat fluxes from two large lakes on south-central
Baffin Island, N.W.T., Canada led to relatively warm (2o C above regional conditions) summer climate
conditions in the interior lowlands relative to dry regions of the interior. Jeffries et al. (1999) found that
conductive heat flow from water below the ice and snow layer near Barrow, Alaska provides a significant
source of winter heat flux to the atmosphere.
The importance of the disparate heat fluxes from land and open water areas has led to the introduction of
lake modules in regional climate models for the prediction of both past and present conditions (Goyette et
al. 2000, Hostetler et al. 2000). For example, Hostetler et al. (2000) used a lake model within a regional
climate model to simulate the effect of Lake Aggasiz on central North American climate 11,000 years ago.
Goyette et al. (2000) included a mixed layer lake model within the Canadian Regional Climate Model
(CRCM) to represent the evolution of water surface temperature and the initiation of lake ice cover, with
the motivation of improving prediction of lake-induced winter precipitation downwind of the Laurentian
Despite recent interest in representing the atmospheric effects of lakes in regional climate models, the
hydrologic effects of lakes were until recently largely neglected in the land surface schemes (LSS) in
GCMs used to predict global climate. This is due in part to the scale differential: at the scale of current-
generation GCM model grid cells, lakes rarely occupy an entire model grid cell, or even a major part of
one. Seven of the 21 LSS that participated in the Project for Intercomparison of Land-surface
Parameterization Schemes (PILPS) Phase 2(e) arctic model intercomparison (Bowling et al. 2003b)
calculated evaporation from open water areas within the model grid cell. Only two of the models (one was
the VIC model on which our work is based) represented the storage retardation effect of lakes on the runoff
In this paper, we describe in detail an algorithm to represent the evaporation and storage effects of lakes
and wetlands within the Variable Infiltration Capacity (VIC) macroscale hydrology model (see Liang et al
1994 and Cherkauer et al. 2003 for details). The model is evaluated with respect to its ability to represent
water temperatures, net radiation, ice freeze/thaw and runoff production for a variety of high-latitude
locations. It is then used to investigate the influence of ponds and wetlands on the spatial and temporal
distribution of water and energy fluxes for the Putuligayuk River watershed, on the Alaskan Arctic coastal
2. Model Description
The lake and wetland algorithm is intended for use within the framework of the VIC macroscale hydrology
model which is intended both for use in off-line simulations of large river basins (Nijssen et al 1997) and in
coupled land-atmosphere studies (Zhu et al 2008). VIC has been applied in many diverse environments,
including a global application at 2o latitude/longitude resolution (Nijssen et al. 2001) and across the
contiguous U.S. at 1/8o resolution (Maurer et al. 2002). Surface runoff is generated according to the
variable infiltration curve, which represents the dependence of infiltration capacity on the spatial
distribution of surface soil moisture (Liang et al. 1994). Baseflow is generated according to an empirically-
based nonlinear soil moisture relationship (Liang et al. 1994). The effects of seasonally and permanently
frozen soil on infiltration and runoff are represented using the method of Cherkauer and Lettenmaier
(1999). Once generated, surface runoff and baseflow are released from the model grid cell and are therefore
no longer available for use within the grid cell. An independent routing model is used to route the runoff
and baseflow to the basin outlet for prediction of streamflow hydrographs (Lohmann et al. 1998a and b).
The model is designed to be applied at scales sufficiently large such that subsurface transfer of moisture
between model grid cells can be neglected, and all transfer of moisture between grid cells occurs via the
runoff routing postprocessor.
The algorithm described here is intended to represent the effects of lakes and wetlands in the VIC model by
creating a surface wetland land class that can be added to the grid cell mosaic, in addition to the vegetation
and bare soil land classes. The wetland class represents seasonally flooded ground as well as permanent
water bodies. Water and energy components of the combined lake 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), and Patterson and Hamblin (1988), while that of the exposed wetland follows Cherkauer
and Lettenmaier (1999).
a. Lake Algorithm Description
Evaporation from the water surface is calculated in each time step by solving a surface energy balance in
the manner of Hostetler and Bartlien (1990) and Hostetler (1991) and summarized here. Energy exchange
with the atmosphere takes place within a surface water layer, which cannot exceed a user-specified depth
(zsurf), typically around 0.6 m. Absorption of solar radiation by the surface layer is assumed to follow
Beer’s law. A two-band system is assumed, so radiation intensity as a function of depth, h, is given by
(Patterson and Hamblin 1988):
I (h ) = I o [Av ⋅ exp (− λv h ) + ANIR ⋅ exp (− λ NIR h )] (1)
where Io is the net shortwave radiation at the water surface, Av and ANIR are the fractions 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 layers is resolved by solving a set of simultaneous
equations. Included in these equations are the effects of radiation absorption by each layer, the eddy
diffusion of heat from adjacent layers by molecular and wind-induced turbulent mixing and convective
mixing due to temperature instabilities (Hostetler and Bartlien 1990). The bottom boundary has a no flux
condition, meaning that energy is not exchanged with the sub-lake soils (this condition could be relaxed in
For the VIC model we have introduced dynamically varying solution layers, such that the layer thickness is
recalculated in each time step in response to variations in total lake liquid water depth, as follows:
d zsurf = z; d z = 0; l = 1 z < zsurf
d zsurf = z 2 ; d z = z 2 ; l = 2 zsurf <= z < 2zsurf (2)
d zsurf = zsurf ; d z =
(z − z )
; l = ( int ) z 2zsurf <= z <Nnodes zsurf
( l − 1) z surf
d zsurf = zsurf ; d z =
(z − z )
; l = N nodes z >= Nnodes zsurf
( N nodes − 1)
Where z is the current lake depth, l is the current number of solution layers, and dz is the thickness of all
solution layers excluding the surface layer. The current thickness of the surface solution layer is dzsurf and
zsurf is the maximum allowable thickness of the surface solution layer. All depths are in meters. 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 of Patterson and Hamblin (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 equation (1) is assumed for 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/water interface (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, the area 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.
b. Wetland Algorithm Description
Unique features of the VIC lakes and wetland algorithm include the interaction of the simulated lake with
the VIC model grid cell and the ability to represent wetlands of varying size. The algorithm can be
summarized as follows (see Figure 1):
• All open water areas within a VIC model grid cell are simulated together as an effective grid cell lake.
• A user-defined fraction of runoff from vegetated areas within the grid cell is diverted to the lake. This
represents the storage retardation effect of lakes on seasonal streamflow.
• Once the new lake level is calculated, runoff is released from the lake as a function of the lake level.
Baseflow is calculated from below the lake as a function of the liquid water content of saturated
• 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
For this work, we use a hydrologic definition of wetlands, specifically meaning areas that are saturated or
covered by water for some portion of the year (Zoltai 1979). The tendency of a region to flood 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 the subsequent 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, the lake
fraction does not distinguish between the pelagic open water zone and the benthic zone that may contain
emergent wetland vegetation.
As the stage of the simulated lake drops, the open water area is recalculated and additional wetland area is
exposed (Figure 2). 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 − Ew ⋅ f lake + Ev ⋅ (1 − f lake ) − Dlake
where ∆S is the change in soil moisture, lake water and ice and snow storage, P is precipitation, Ew and Ev
are evaporation from the open water and wetland vegetation, respectively. Dveg is the discharge (runoff and
baseflow) entering into the lake from the non-wetland portion of the gridcell and Dlake is the discharge out
of the lake. All of the runoff and baseflow generated by the exposed wetland area is assumed 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 of saturated soil (Figures 2b and c).
Likewise, as the lake area expands, some of the lake volume must go to saturating the newly inundated soil
(Figure 2d). The average soil moisture for the wetland fraction is therefore updated as follows:
W = (1 − f lake ) ⋅ Wv + flake ⋅ Wmax ′
f lake ≤ f lake (4)
W = (1 − f lake ) ⋅ Wv + f lake ⋅ Wmax
′ ′ ′
f lake > f lake
where f lake and f lake represent the new and old lake fraction, respectively. Wv is the wetland soil moisture
for the current timestep and Wmax is the total soil moisture storage capacity of the soil column. The volume
of water (as mm per unit area) from the expanding lake used to recharge the wetland soil moisture is
calculated as follows:
Recharge = ( f lake − f lake ) ⋅ (Wmax − Wv )
f lake > f lake (5)
Changes in lake stage are calculated via a water balance for the saturated lake area. Runoff into the lake is
composed of all runoff and baseflow from the exposed wetland and a fraction of the runoff and baseflow
from all other grid cell vegetation types, as follows:
Din = Dveg + Dwetland = α ⋅ ∑ ( Rv ( i ) + Bv ( i ) ) ⋅ Cv ( i ) + ( Rw + Bw ) ⋅ (1 − f lake ) ⋅ Cwet (6)
where Cv is the fraction of the grid cell occupied by vegetation type i, and Rv and Bv are the runoff and
baseflow from vegetation type i. α is a user-specified parameter controlling the portion of runoff from the
vegetated areas that flows into the lake. Rw and Bw are the runoff and baseflow from the exposed wetland
To avoid complications due to the variation of lake area with depth, lake volume is the state variable used
for the water balance. Lake depth is updated each time step by piecewise integration of the derived depth-
volume curve. Subsurface outflow from the lake 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 is reduced by the ice content of the exposed wetland soil profile in order to calculate the baseflow, as
described by Cherkauer and Lettenmaier (1999).
Surface outflow from the lake is calculated as a function of the new depth, based on the equation for flow
over a broad-crested weir, assuming that the velocity head is negligible:
Q = cd b g ( z − zmin ) 2 (7)
Where Q is the discharge in m3 s-1, b is the flow width (m), g is the acceleration due to gravity, z is the
current lake depth and zmin is the elevation above the lake bottom of the weir or lake outlet. The coefficient
of discharge, cd, is used to account for the velocity of approach, non-parallel streamlines over the crest, and
energy losses. It 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 (in m/time step) becomes:
1.6b ( z − zmin ) 2 ⋅ dt
Rlake ( z ) = (8)
Where dt is the time step length in seconds and A(z) is the lake surface area (m2) at depth, z. For natural
lakes and wetlands the width of the reservoir outlet is also likely to vary with water level. Assuming a
roughly circular lake, the flow width, b, can be expressed as a fraction of lake circumference:
b = 2 f π ⋅ A( z ) (9)
Where f is the fraction of the lake circumference. Surface runoff out of the lake (in m per time step)
0 z ≤ zmin
Rlake(z) = (10)
5.67 f ( z − zmin ) 2 ⋅ dt
A( z ) z > zmin
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 displacement of water by floating ice. When
the mass of ice exceeds the mass of liquid water, z is the depth of liquid water alone. The depth-area curve,
A(z); width fraction, f; and 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 calibration.
3. Lake Model Evaluation
The core lake thermal model was first evaluated with respect to water temperature, radiation balance and
ice cover using available datasets from several high latitude locations. Physical constants governing
radiation extinction and thermal conductivity were set during this evaluation, as summarized in Table 1.
These same values were used for all subsequent simulations.
Figures 3 and 4 compare model simulations and observations of net radiation and water temperature from
the largest of three thaw ponds monitored on the Arctic coastal plain at Prudhoe Bay, AK (148.8957o W,
70.2795o N) for the summers of 1999 and 2000. Other meteorological components needed to run the VIC
model, including downward solar and longwave radiation, wind speed and humidity are monitored at the
nearby Betty Pingo meteorological tower maintained by the Water and Environmental Research Center
(WERC) at the University of Alaska Fairbanks (UAF). Net radiation was monitored using a Q-7.1 net
radiometer from Campbell Scientific. Wind induced convective cooling as air moves over the sensors can
cause a reduction in the measured net radiation of up to 6 percent for positive fluxes and 1 percent for
negative fluxes. The observed net radiation was adjusted for this bias using observed wind speed from the
UAF meteorological tower. Since the ponds are shallow (1 to 1.5 m), they are not affected by stratification
and diffusive effects, and so provide a good evaluation of the surface energy balance alone. Figure 3 shows
that the energy balance calculation is able to represent both the daily average surface water temperature and
net radiation of the thaw pond. The mean diurnal cycle of net radiation is also well represented, although
the diurnal cycle of the simulated water temperatures is suppressed relative to the observations (Figure 4).
It is suspected that the large measured diurnal temperature variations are due to inadequate shielding of the
submerged temperature sensor in these clear, shallow waters.
Figure 5 illustrates the average (1993-1999) simulated and observed water temperature profile of Toolik
Lake, Alaska, for the ice-free period. Toolik Lake is a 1.5 km2 kettle lake located in the foothills of the
Brook’s Range (68o 38’ N, 149o 43’ W). Its maximum depth is 25 m (7 m average) and it is ice-free from
July through September. Periodic lake temperatures and ancillary meteorological data have been collected
by researchers at the Marine Biological Laboratory of the Woods Hole Oceanographic Institution since
1975 as part of the Arctic Long-Term Ecological Research project at the Toolik Field Station. Figure 5
illustrates that the VIC lake algorithm captures several important features of the observed water
temperature profile. One of these is the depth of the thermocline (approximately 5 to 7 m) which is
controlled in part by the depth of penetration of solar radiation. Another is the gradual downward migration
of heat driven by eddy diffusion. Finally, the influence of convective mixing is apparent in the spring and
fall, when the lake completely overturns and becomes isothermal in both the simulated and observed
profile. On average, the simulated profile becomes isothermal approximately 2 weeks later than observed.
The ability of the ice model to represent the depth and duration of ice cover for a range of arctic lakes was
evaluated using data from northern Sweden. Figure 6 compares simulated and observed ice thickness on
Torne Lake (Torneträsk; 68.1956o N, 19.9930o E) from January 1989 through December 1998. Torneträsk
is a large (330 km2) glacial lake in the high mountain region of Sweden with a mean depth of
approximately 150 meters. Ice thickness observations are collected weekly by the Swedish Meteorological
and Hydrological Institute (SMHI). The simulated data represents a fraction of the VIC model grid cell
centered at 19.875o E, 68.125o N. Meteorological data for this simulation came from the PILPS 2(e) model
intercomparison (Bowling et al. 2003b). As shown in Figure 6, the timing of ice melt is well-represented,
although maximum ice depth is over-predicted, apparently due to a bias in the initiation of ice-cover in the
fall. This bias is likely due to the inability of the simple one-dimensional model to represent the impact of
currents and mechanical action, which suppress continuous ice formation. Also shown in Figure 6 is an
alternative scenario in which ice cover was prevented from forming until January 1st of each year. In this
case, total ice thickness is better represented in the model.
The simulated dates of lake freeze-up and break-up for the period 1989 to 1998 on 14 lakes within the
Torne-Kalix River basin in northern Sweden are shown in Figure 7. These dates are compared to observed
dates of ice formation and break-up, which were collected by the SMHI. For the simulated ice cover, the
dates of freeze-up and break-up are defined, respectively, as the date on which the ice cover first exceeds or
falls below 75% (for the “effective lake” in the 1/4o VIC model grid cell corresponding to the observation
location). On average, the VIC algorithm appears to overpredict the duration of the ice-covered period for
these 14 lakes. In contrast to the Torneträsk simulations, in the small lakes the VIC model captures the date
of ice formation well, but there is a tendency for it to simulate a break-up date approximately two weeks
later than that observed. This delay is not particularly sensitive to the ice cover threshold used to define
break-up, but may be a result of a delay in the simulated snowmelt.
These four test scenarios have demonstrated the functionality of the water thermal algorithm, as
implemented in the VIC model. The remainder of this work will evaluate the dynamic wetland algorithm
with respect to a case study of the Putuligayuk River, AK.
4. Surface Water Storage in the Putuligayuk River Basin
The Putuligayuk River is an approximately 471 km2 watershed located entirely on the Arctic coastal plain
in northern Alaska (see Figure 8). The annual runoff hydrograph for the Putuligayuk River is dominated by
snowmelt and upwards of 90 percent of the annual discharge occurs in the two weeks following snowmelt.
The watershed contains numerous thaw ponds and other permafrost features that impede drainage. As a
result between 50 to 70 percent of the watershed is flooded seasonally (Bowling et al. 2003a). Bowling et
al. (2003a) investigated the role of surface storage in ponds and wetlands on the hydrologic response of the
Putuligayuk River and found that between 30 to 44% of the end of winter snow water equivalent (SWE) is
not immediately available for runoff due to surface impoundment. The maximum extent of surface water
was observed to increase by 175 to 338% over mid-summer values in 1999 and 2000 respectively, due to
the influence of snowmelt, as determined from classification of RADARSAT ScanSAR imagery (Bowling
et al. 2003a). Between 42 to 43% of the stored water is returned to the atmosphere by evaporation from
open water areas (Bowling et al. 2003a). Following the large reduction in surface water extent, the
drainage network becomes fragmented, and discharge from the watershed approaches zero.
a. Wetland Geometry
The physical description of the pond and wetland extent, A(z), for each model grid cell was based on
observations and consists of two primary components: the extent of permanent ponds and the extent of
seasonally flooded area, as illustrated in Figure 9. The bathymetry of three typical thaw ponds on the
coastal plain was determined using a laser-surveying station (total station) and inflatable raft. Three
transects were made of each pond. The surveyed cross-sections were used to generate elevation contours
within the Arc/Info geographical information system. The elevation contours were converted into a depth-
area relationship for all three ponds and the area was normalized by the maximum extent of each pond
(Amin in Figure 9). Although the three ponds differ in maximum area (6,000, 7,000 and 16,000 m2), the
maximum depths are all between 1 and 1.5 meters, due to the influence of the permafrost layer and the
scaled depth-area relationships are all similar (open circles in Figure 10). It was therefore assumed that the
fitted polynomial profile of pond fractional area versus depth is representative of the collective lakes and
ponds of the coastal plain, with zmin=1.25.
The depth-area relationship for the wetland portion of the watershed was derived using a time-series of
RADARSAT ScanSAR images from the summers of 1999 and 2000 classified into water and non-water
areas, as described by Bowling et al (2003a). The classified time series (shown in Figure 11), shows the
rapid decline in inundated area following snowmelt in both years, before the water level approaches zmin.
At this point, streamflow decreases rapidly and further reduction in pond depth is due to evaporation and
limited subsurface drainage. For each year (i) and gridcell (j), the maximum open water extent ( Amax ) was
identified and the open water extent at the end of rapid runoff ( Amin ) was inferred as the inflection point of
the saturated area time series (Figure 11). The volume of streamflow between these two dates (expressed
in meters) is assumed to come entirely from this change in stored water (Vdischarge in Figure 9). Assuming a
linear relationship between area and depth for simplicity, the slope dA ( z ) dz can be calculated from:
dA ( z )
max − Amin )( Amax + Amin )
i, j i, j i, j
dz 2 ⋅Vdisch arg e ⋅ Acell
Where Acell is the area of each grid cell in square kilometers. Combining this with the derived curve for
the pond bathymetry yields a theoretical depth-area relationship for the lake and wetland portions for each
grid cell of the Putuligayuk catchment, as follows:
Aj z A < Amin
min z min
Aj (z) = (12)
A j + dA ( z ) z − zA > Amin
min ( min )
where A(z) is the lake area, and dA ( z ) dz is the average of the slopes calculated from 1999 and 2000.
The fractional area curve that is input to the VIC model is obtained by dividing equation (7) by the grid cell
area, Acell . Maximum and minimum inundated fractions (Amax and Amin) were determined for each grid cell
based on the classified RADARSAT images. The maximum saturated area possible in the VIC model was
set equal to 5% more than that observed by RADARSAT in 2000, while Amin was set to the minimum value
for 1999. This approach assumes that the wetland storage available in each grid cell is constant, and
therefore the surface slope must be steeper in grid cells that experience less change in inundated area. In
reality, it is likely that both volume and slope vary by grid cell, but there are insufficient data to derive such
As described above, the depth corresponding to Amin, zmin was assumed constant for each grid cell. Using
Amin, Amax and zmin, zmax was calculated for each grid cell from equation (7). The observed wetland areas at
maximum extent varied from 20 to 84 percent of the model grid cells, while at minimum extent wetland
area varied from 6 to 29 percent.
b. Model Implementation
The Putuligayuk River was simulated using the VIC model with the lakes and wetlands algorithm at a grid
cell resolution of 1/8o latitude by longitude (approximately 13.9 km by 4.7 km) and an hourly time step.
The simulation period for the Putuligayuk was August 1987 through July 2001, with the first simulated
year not used for analysis. Gridded meteorological data, including hourly vapor pressure, air temperature
and wind speed and daily precipitation was prepared as described in Bowling et al. (2004) for the adjacent
Kuparuk River basin, and summarized here. Daily total precipitation was obtained from two National
Climatic Data Center (NCDC) Co-operative observer stations located at the Prudhoe Bay and Kuparuk air
fields and one Natural Resource Conservation Service snow telemetry site, located at Sagwon (see Figure
8). NCDC precipitation was corrected for gauge undercatch using the method of Yang et al. (1998) for an
unshielded National Weather Service 8 inch standard gauge. The corrections for this wind swept
environment resulted in a 70 percent increase in cumulative precipitation over the simulation period and are
a major source of uncertainty in winter snow accumulation. Hourly relative humidity, air temperature and
wind speed were obtained from the four nearest WERC meteorological stations. Downward shortwave
and longwave radiation are not monitored in this remote location during the winter months, so they were
estimated internally to the VIC model, using the method of Thornton and Running (1999).
Soil freezing and thawing is represented with a constant temperature lower boundary condition set to -10o
C at a depth of 4 m using the algorithm of Cherkauer and Lettenmaier (1999). The boundary condition was
set based on reasonable reproduction of the observed soil temperature profile at the WERC Betty Pingo site
(Figure 12). Sublimation from blowing snow is also represented, using the parameters from Bowling et al.
(2004). Decay in snow albedo over time is represented using the algorithm of Sun et al. (1999). New snow
albedo and the snow water holding capacity were adjusted in comparison with snow surveys (Figure 12), as
summarized in Table 2.
Stream discharge measurements were re-initiated by UAF in 1999 after being discontinued by the USGS in
the 1980’s. Observed discharge in the Putuligayuk River is currently available for 1999 through 2001. The
period from September 1998 through September 1999 was used for calibration of the streamflow
hydrograph, while observed saturated extent for 1999 and 2000 from the RADARSAT-based estimates of
Bowling et al. (2003a) were used for calibration of seasonally-flooded area. Calibration of streamflow with
the VIC model typically involves adjustment of the variable infiltration parameter, bi, the three baseflow
parameters (Ws, Ds and Dsmax) and sometimes the thickness of the second and third soil layers. The lake
model adds three potential calibration parameters, the width fraction, f, the fraction of upland vegetation
that drains through the lake, α, and zmin, the minimum depth for surface flow from the lake. As summarized
in Table 3, a limited single factor sensitivity analysis was conducted to assess the influence of these
parameters (excluding α which was always set equal to 1.0) on the annual peak discharge and maximum
flooded area. Due to the presence of continuous permafrost which limits infiltration and subsurface flow
contribution, the simulated spring snowmelt hydrograph is relatively insensitive to the traditional VIC
calibration parameters. Both the peak discharge and peak flooded area are somewhat inversely correlated
with parameters controlling baseflow, Ws, Ds, Dsmax and Depth 3, so that as summer baseflow is increased,
the subsequent spring flood is reduced. As expected, the simulated streamflow and flooded area are most
sensitive to the new lake parameters, as illustrated in Figure 12. In this application, zmin was fixed in the
derived wetland geometry described above and f was adjusted during manual calibration of the streamflow
hydrograph and saturated area curves.
c. Model Evaluation
As shown in Figure 12, the timing of snowmelt and subsequent active layer development are well-
represented for this coastal plain location. Observations of SWE in Figure 12 are from snow surveys at the
Betty Pingo monitored wetland complex adjacent to the Putuligayuk catchment. Simulated SWE is from
the model grid cell centered at 148.9375o W, 70.3125o N from August 1996 through August 2000. The
model shows both under- and over-prediction of maximum SWE for these years, with observed maximums
of 133, 99 and 84 mm, respectively for 1997, 1999 and 2000. Simulated maximum SWE for these years
was 149, 81 and 86 mm. The predicted SWE is based on the catch-corrected and interpolated winter
precipitation less sublimation/condensation from blowing snow and the surface snow pack. Errors in any
of these quantities may account for the discrepancy in simulated SWE.
Snowmelt typically occurs in late May and early June with the active layer beginning to thaw once the
snow disappears. The observed 0o isotherm calculated from temperature profiles measured by WERC at
both the wetland site and an adjacent upland site are also shown in Figure 12. Both sites thaw at
approximately the same rate for the first 30 to 50 cm. Below this depth, the wetland site thaws more slowly
due to the greater phase change energy required to thaw the saturated soils. The maximum depth of the
active layer in the upland site is unknown for each of the three years shown, since the temperature exceeds
zero at the lowest measurement point (60 cm). The active layer at the wetland site is typically shallower,
and varies between 50 and 70 cm (maximum measurement depth 75 cm). The VIC-simulated 0o isotherm
is an average of both the wetland and the upland portions of the grid cell and is consistent with the
observations in both timing and the depth of thaw. The model appears to be more similar to the observed
wetland profile, which is reasonable since 75% of this grid cell is considered to be wetland.
Detailed observations of the duration of ice cover on the lakes and ponds are not available. The duration of
ice cover varies with lake size and depth, with ice cover typically disappearing in the latter half of June (D.
Kane, personal communication). The simulated last day of ice cover averaged over all cells varied between
14 June and 12 July between 1989-2001. The mean ice-free date was 25 June, and in most years the ice
was gone by the end of June. For the simulated years, the difference in ice-free date between the first and
the last grid cell to melt was as few as 5 days to as many as 25 days. Simulated mean freeze-up dates (first
date with 100 percent ice cover) varied between 8 September and 8 October during the 1988 -2000 test
period, with a mean date of 22 September.
The rapid snowmelt on top of frozen soils illustrated by Figure 13 results in extensive flooding of the low-
gradient coastal plain. The simulated and observed saturated extents for 1999 and 2000 from the
RADARSAT-based estimates of Bowling et al. (2003a) are shown in Figure 11. The magnitude of
maximum saturated extent for both years matches well with the RADARSAT estimates, with an
underprediction in 2000 (52% of simulated area versus 60% observed). The rate of recession in saturated
area in 1999 is not as steep as the observations, possibly due to the presence of ice. In addition, in 1999
both the simulated saturated area and the simulated discharge (see Figure 14) are too responsive to late
summer rain events.
In general, the model reflects the seasonal cycle inferred from previous field studies, with the saturated
extent reaching an annual maximum a few days following snowmelt. The maximum surface water extent
for 1989-2001 varied between approximately 37 and 58 percent of the total simulated area. This represents
a 103% to 218% increase over mid-summer values. This is less than the 175 to 338% increase estimated
from RADARSAT for 1999 and 2000, and reflects in part an underprediction of draw-down during the
summer months, resulting in an overprediction of minimum flooded area.
Observed discharge for the Putuligayuk River is shown in Figure 14, along with the discharge predicted by
the VIC model both with and without the wetland algorithm. For both scenarios and in both years, the
predicted snowmelt hydrographs peak substantially lower than observed. The total streamflow volume for
the period of observations in each year, is overpredicted in 1999 (80 mm versus 64 mm observed) and
underpredicted in 2000 (113 mm versus 121 mm observed). These differences seem to be related to
discrepancies in the gridded precipitation inputs, as discussed below. The primary influence of surface
water storage in the Putuligayuk basin is illustrated by the difference between the two simulated
hydrographs. The VIC wetland algorithm is better able to capture the delay in streamflow and confines the
spring freshet to a single peak with greater magnitude, rather than multiple small peaks. The wetland
algorithm also buffers the streamflow response in the late summer, relative to the simulation without
a. Precipitation Inputs
In order to produce a continuous precipitation record for the duration of the simulations presented here,
precipitation was taken from NCDC Cooperative observer stations and the NRCS SNOTEL site at Sagwon.
Bowling et al. (2003a) utilized summer precipitation measured by WERC at the wetland study site. The
summer precipitation used in the simulations was 11 mm (11%) higher in 1999 than the WERC values and
13 mm (22%) lower in 2000. Such discrepancies further highlight the difficulties in precipitation
measurement in remote areas such as the Alaska North Slope.
The comparisons with end-of-winter SWE at the Betty Pingo site shown in Figure 12 suggest that on
average, the catch-corrected, gridded precipitation dataset reasonably represents winter precipitation at this
location, but perhaps not for locations further from the coast. The observed watershed average end-of-
winter SWE measurements were 87 mm in 1999 and 124 mm in 2000 (Bowling et al. 2003a), while at
Betty Pingo, the end-of-winter SWE was 99 and 84 mm, respectively. The simulated watershed average
end-of-winter SWE was 107 and 84 mm, respectively. Due to the position of the available NCDC stations,
the watershed simulations are strongly influenced by the weather conditions recorded at the coastal plain
wetland site, and may not adequately reflect differences in precipitation patterns further from the coast,
where precipitation tends to increase (see e.g., Bowling et al. 2004). This suggests that the observed
watershed average SWE values better reflect the watershed average precipitation inputs, than do the point
measurements reflected in Figure 12, and the gridded precipitation data used in this study.
b. Surface Water and Energy Fluxes
The role of wetland storage in the water balance of this low-gradient arctic river was further explored by
investigating components of the simulated water balance for water years 1988 through 2000. The percent
of maximum snow water equivalent entering storage (rather than contributing to direct runoff) averages
20% on an annual basis and varied between 0 and 51% over the test period. For this purpose, direct runoff
was computed as the sum of all surface runoff before the hydrograph inflection point, identified as the point
when the rate of change of surface runoff drops below 1 mm day-1. This simulated rate of recharge
compares well with the 30 to 44% inferred from the basin water balance for 1999-2001.
Over the summer, some of this stored water is returned to the atmosphere via evaporation from lake and
wetland surfaces, which on average exceeds precipitation by 24 mm (between 8 and 37 mm yr-1 for
individual grid cells). Evaporation from the lakes and wetlands increases the prediction of seasonal (July,
August, September) evaporation from individual grid cells by between 4 and 17 mm depending on the total
lake area of the grid cell, or by 7 to 39 percent of the seasonal evaporation without lakes and wetlands.
Including the representation of open water evaporation in the VIC model only increases the watershed
average annual evaporation rate by 8%, since evaporation is generally not water limited in this
environment. What may have more influence is the change in spatial heterogeneity in evaporation and
latent heat, which can influence the generation of summer convective storms and the position of the Arctic
frontal zone in regional weather and climate modeling.
This paper describes an algorithm to represent the hydrologic effects of ponds and wetlands within a
macroscale hydrology model. The algorithm allows for the growth and reduction of saturated area over
time. Although our focus here is on the representation of arctic environments, the algorithm is a potentially
useful tool for studying the interaction of water fluxes between wetlands and adjacent surface water bodies
in other climates and regions as well. Our evaluation with respect to observed conditions for lakes in
Alaska and Sweden showed that:
• The one-dimensional model was able to reasonably approximate the temperature and ice conditions
of high-latitude lakes.
• The model does not take into account the mechanical action of wind, waves and currents, and as a
result over-predicts the ice-covered period somewhat.
To evaluate the wetland algorithm, the Putuligayuk River watershed in northern Alaska was simulated from
1987 to 2000. This watershed is a dynamic system, in which the total basin saturated area varies by over
50 percent during the snow-free period. The combined influences of permafrost, low topographic
gradients, snow cover and numerous ponds and wetlands make it a difficult system to represent with a
traditional hydrology model. The depth-area relationship for the wetland portion of the watershed was
derived from a time-series of classified RADARSAT ScanSAR images. The bathymetry of the lakes
themselves was assumed based on surveys of three typical thaw ponds on the coastal plain. The wetland
algorithm was able to reproduce the observed change in basin saturated area on a seasonal cycle and
improve the timing and shape of the snow-melt induced hydrograph. The simulated results suggest the
following about the role of surface storage in this system:
• Up to 50% of snow melt water goes into storage each year and does not contribute to streamflow.
• Seasonally, the area of open water increases on the order of 100 – 220% due to the influx of snowmelt
• Wetland storage is reduced by evaporation from open water, which exceeds summer precipitation by 9
The watershed-average SWE and total streamflow are underpredicted, particularly in 2000, suggesting that
the model forcing data does not correctly capture an increasing precipitation gradient away from the coast.
This makes it difficult to explore the absolute values of simulated water balance components in this
environment. On going research will focus on evaluating the wetland algorithm in temperate
environments, as well as relaxing the thermal boundary condition at the water / soil interface, in order to
better explore the interaction of the soil moisture and thermal regimes in permafrost environments.
We would like to acknowledge the assistance of Keith Cherkauer and Ted Bohn for their role in merging
the lake algorithm and its updates with the VIC model. We are also indebted to Doug Kane and the Water
and Environmental Research Center at the University of Alaska, Fairbanks for access to field sites and
much of the data used for model simulation and evaluation. The early work was funded in part by a NASA
Earth System Science fellowship to the first author and NASA Grant NAG5-10092 to the University of
Washington. More recently, this work was supported by NASA Grant NNG05GD17G to the University of
Washington, and Purdue University via sub-contract.
Bigras, S.C., 1990: Hydrological regime of lakes in the Mackenzie delta, Northwest Territories, Canada.
Arctic and Alpine Research, 22, 163-174.
Bowling, L.C., D.L. Kane, R.E. Gieck, L.D. Hinzman, and D.P. Lettenmaier, 2003a: The role of surface
storage in a low-gradient Arctic watershed. Water Resour. Res., 39 (4), 1087-1099,
Bowling, L.C., D.P. Lettenmaier, B. Nijssen, L.P. Graham, D.B. Clark, M. El Maayar, R. Essery, S. Goers,
Y.M. Gusev, F. Habets, B. van den Hurk, J. Jin, D. Kahan, D. Lohmann, X. Ma, S. Mahanama, D. Mocko,
O. Nasonova, G. Niu, P. Samuelsson, A.B. Shmakin, K. Takata, D. Verseghy, P. Viterbo, Y. Xia, Y. Xue
and Z. Yang, 2003b: Simulation of high-latitude processes in the Torne-Kalix basin: PILPS Phase 2(e) 1:
Experiment description and summary intercomparisons. J. Global and Planetary Change, 38 (1-2), 1-30.
Bowling, L.C., J.W. Pomeroy and D.P. Lettenmaier, 2004: Parameterization of blowing snow sublimation
in a macroscale hydrology model, J. Hydrometeorology. 5 (5), 745-762. doi: 10.1175/1525-7541.
Brutsaert, W., 1982: Evaporation into the atmosphere. D. Reidel, Dordrecht.
Cattle, H., 1985: Diverting Soviet rivers: some possible repercussions for the Arctic Ocean. Polar Record,
Cherkauer, K.A. and D.P. Lettenmaier, 1999: Hydrologic effects of frozen soils in the Upper Mississippi
River basin. J. Geophys. Res., 104d, 19599-19610.
Cherkauer, K. A., L. C. Bowling and D. P. Lettenmaier, 2003: Variable Infiltration Capacity (VIC) Cold
Land Process Model Updates. J. Global and Planetary Change, 38(1-2), 151-159.
Glenn, M.S. and M.K. Woo, 1997: Spring and summer hydrology of a valley-bottom wetland, Ellesmere
Island, Northwest Territories, Canada. Wetlands, 17 (2), 321-329.
Goyette, S., N.A. McFarlane, and G.M. Flato, 2000: Application of the Canadian Regional Climate Model
to the Laurentian Great Lakes region: Implementation of a lake model. Atmos.-Ocean, 38 (3), 481-503.
Hamill, L., 2001: Understanding Hydraulics. 2nd ed, Palgrave, 632p.
Hostetler, S.W., 1991: Simulation of lake ice and its effect on the late-Pleistocene evaporation rate of Lake
Lahontan,.Climate Dyn., 6, 43-48.
Hostetler, S.W. and P.J. Bartlein, 1990: Simulation of lake evaporation with application to modeling lake
level variations of Harney-Malheur Lake, Oregon. Water Resour. Res., 26 (10), 2603-2612.
Hostetler, S.W., P.J. Barlein, P.U. Clarke, E.E. Small and A.M. Soloman, 2000: Simulated influences of
Lake Agassiz on the climate of central North America 11,000 years ago. Nature, 405, 334-337.
Jacobs, J.D. and L.D. Grondin, 1988: The influence of an Arctic large-lakes system on mesoclimate in
south-central Baffin Island, N.W.T., Canada. Arctic and Alpine Research, 20 (2), 212-219.
Jeffries, M.O., T. Zhang, K. Frey and N. Kozlenko, 1999: Estimating late-winter heat flow to the
atmosphere from the lake-dominated Alaskan North Slope. J. Glaciol., 45 (150): 315-324.
Liang, X., D.P. Lettenmaier, and E.F. Wood, 1996: One-dimensional statistical dynamic representation of
subgrid spatial variability of precipitation in the two-layer Variable Infiltration Capacity model. J.
Geophys. Res., 101d, 21403-21422.
Liang, X., Lettenmaier, D.P., Wood, E.F. and Burges, S.J., 1994: A simple hydrologically based model of
land surface water and energy fluxes for general circulation models. J. Geophys. Res., 99 d, 415-428.
Lohmann, D., Raschke, E., Nijssen, B. and Lettenmaier, D.P., 1998a: Regional scale hydrology: I.
Formulation of the VIC-2L model coupled to a routing model. Hydrological Sciences Journal, 43, 131-141.
Lohmann, D., Raschke, E., Nijssen, B. and Lettenmaier, D.P., 1998b: Regional scale hydrology: II.
Application of the VIC-2L model to the Weser River, Germany. Hydrological Sciences Journal, 43, 143-
Maurer, E.P., A.W. Wood, J.C. Adam, D.P. Lettenmaier and B. Nijssen, 2002: A long-term hydroligally-
based data set of land surface fluxes and states for the continental United States. J. Climate, 15, 3237-
Marsh, P. and Bigras, S.C., 1988: Evaporation from Mackenzie delta lakes, N.W.T., Canada. Arctic and
Alpine Research, 20 (2), 220-229.
Mendez, J., L.D. Hinzman and D.L. Kane, 1998: Evapotranspiration from a wetland complex on the Arctic
coastal plain of Alaska. Nordic Hydrology, 29 (4/5), 303-330.
Moskvin, Y.P., 1989: Runoff from hummocky marshes of western Siberia. Soviet Meteorology and
Hydrology, 3, 88-94.
Nijssen, B., D. P. Lettenmaier, X. Liang, S. W. Wetzel, and E. F. Wood, 1997: Streamflow simulation for
continental-scale river basins. Water Resour. Res., 33, 711-724.
Nijssen, B., R. Schnur, and D.P. Lettenmaier, 2001: Global retrospective estimation of soil moisture using
the Variable Infiltration Capacity land surface model, 1980-1993. J. Climate, 14, 1790-1808.
Patterson, J.C. and P.F. Hamblin, 1988: Thermal simulation of a lake with winter ice cover. Limnol.
Oceanogr., 33(3), 323-338.
Roulet, N.T. and M.K. Woo, 1986: Wetland and lake evaporation in the low Arctic. Arctic and
Alpine Research, 18: 195-200.
Stele, M., 1996: A simple model of the Arctic Ocean freshwater balance, 1979-1985. J. Geophys. Res.,
Rovansek, R.J., Hinzman, L.D. and Kane, D.L., 1996: Hydrology of a tundra wetland complex in the
Alaskan Arctic coastal plain, U.S.A.. Arctic and Alpine Research, 28, 311-317.
Thornton, P.E., and S.W. Running, 1999: An improved algorithm for estimating incident daily solar
radiation from measurements of temperature, humidity and precipitation. Agricultural and Forest
Meteorology, 93, 211-228.
UNESCO, 1978: World water balance and water resources of the earth, United Nations Educational,
Scientific and Cultural Organization.
Vuglinsky, V.S., 1997: River water inflow to the Arctic Ocean – conditions of formation, time variability
and forecasts, Proceedings, Conference on Polar Processes and Global Climate, Part II of II, 275-276.
Woo, M.K. and Xia, Z., 1996: Effects of hydrology on the thermal conditions of the active layer. Nordic
Hydrology, 27, 129-142.
Young, K.L. and M.K. Woo, 2000: Hydrological response of a patchy high Arctic wetland. Nordic
Hydrology, 31 (4/5), 317-338.
Zhu C.M., L. R. Leung, D. Gochis, Y. Qian and D. P. Lettenmaier, 2008: 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-Discussion (in review).
Zoltai, S.C., 1979: An outline of the wetland regions of Canada, in Proc. workshop on Canadian wetlands
12, Environ. Can., Lands Dir., Ecol. Land Classif. Ser., pp. 1-8.
Figure 1. Schematic of the VIC lake algorithm. I: Evaporation from the lake is calculated 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.
Figure 2. Schematic for the wetland algorithm: a) when the lake is at its maximum extent the soil
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
Figure 3. Simulated and observed daily average a) water temperature and b) net radiation for a
thaw pond on the Alaskan Arctic coastal plain, June – August, 2000
Figure 4. Simulated and observed mean diurnal cycle for a) net radiation and b) water
temperature for a thaw pond on the Alaskan Arctic coastal plain, averaged for the period June –
Figure 5. Simulated (left) and observed (right) mean ice-free season water temperature profile
(1993-1997) for Toolik Lake in the foothills of the Brooks Range in northern Alaska.
Figure 6 Simulated and observed ice thickness for Torne Lake (68.1956o N, 19.9930o E),
northern Sweden, for the period 1989-1998.
Figure 7. Simulated and observed day of a) lake ice break-up and b) lake freeze-up for the period
1989-1998 for 14 lakes in the Torne-Kalix basin, Sweden.
Figure 8. Location of the Putuligayuk River watershed on the Alaskan Arctic coastal plain (inset)
and routing network for the VIC model at 1/8o resolution.
Figure 9. Schematic illustrating the parameterization of stored surface water volume and
inundated area for the Putuligayuk River application
Figure 10. Average depth-area relationship (bold line) using RADARSAT observations of water
extent in individual grid cells (black diamonds), and surveyed depths for ponds (open circles).
Figure 11. Simulated total saturated extent versus time for the region including the Putuligayuk
River watershed.for 1999 through 2000 (solid line); total inundated area on specific dates derived
from RADARSAT ScanSAR imagery (Bowling et al. 2002a) is shown by black circles.
Figure 12. Sensitivity of simulated maximum annual flooded area (left) and discharge (right) to
the new model parameters of width fraction (top) and zmin (bottom). Annual values (1995-2000)
are shown by dashed lines, with the multiyear mean shown as solid black.
Figure 13. Simulated and observed SWE (top) and active layer depth (bottom) for the wetland
study site adjacent to the Putuligayuk catchment (1997-2000).
Figure 14. Simulated and observed discharge for the Putuligayuk catchment, 1999-2001.
Table 1. Lake physical parameters used in evaluation
Coefficients for radiation Visible NIR
attenuation (m-1) (m-1)
Ice, λice 1.5 20
Water, λwater 0.3 1.4
Snow, λsnow 6 20
Thermal conductivity of ice 2.3 W/m·K
Thermal conductivity of snow 0.7 W/m·K
Table 2. Calibration parameters for the Putuligayuk River
New snow albedo 0.95
Max. liquid storage in snow 0.1 * SWE
Min. SWE for new albedo 7 mm/day
Lake hydraulic parameters
Width scaling ratio, f 0.015
Min. depth for runoff, Zmin 1.25 m
Max. allowable depth, Zmax 1.9 m to 2.1 m
Table 3. Relative sensitivity of simulated snowmelt-induced flooding extent and maximum
discharge rates to VIC model parameters
Parameter Range of Values Change in average Change in average
max. flooded area 1 peak discharge1
bi 0.00001 – 0.4 2.0 % 4.5 %
Ws 0.5 – 0.99 9.1 % 6.4 %
Depth 2 0.15 – 1.0 m 7.2 % 11.9 %
Ds 0.001 – 0.1 20.9 % 12.8 %
Depth 3 0.4 – 2.0 m 11.0 % 28.8 %
Dsmax 10.6 – 100 mm/day 28.1 % 26.3 %
zmin 1.12 – 1.38 m 46.8 % 11.8 %
f 0.001 – 0.5 92.8 % 46.3 %
Represents the absolute value of percent change and does not imply a direction of change. The table is
meant to illustrate relative sensitivity, rather than direction since in some cases the change is not monotonic
across the full range of parameters.
Figure 1. Schematic of the VIC lake algorithm. I: Evaporation from the lake is calculated 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.
Figure 2. Schematic for the wetland algorithm: a) when the lake is at its maximum extent the soil 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
Figure 3. Simulated and observed daily average a) water temperature and b) net radiation for a thaw
pond on the Alaskan Arctic coastal plain, June – August, 2000
Figure 4. Simulated and observed mean diurnal cycle for a) net radiation and b) water temperature for a
thaw pond on the Alaskan Arctic coastal plain, averaged for the period June – August, 2000
Figure 5. Simulated (left) and observed (right) mean ice-free season water temperature profile (1993-
1997) for Toolik Lake in the foothills of the Brooks Range in northern Alaska.
Figure 6 Simulated and observed ice thickness for Torne Lake (68.1956o N, 19.9930o E), northern
Sweden, for the period 1989-1998.
Figure 7. Simulated and observed day of a) lake ice break-up and b) lake freeze-up for the period 1989-
1998 for 14 lakes in the Torne-Kalix basin, Sweden.
Figure 8. Location of the Putuligayuk River watershed on the Alaskan Arctic coastal plain (inset), and
routing network for the VIC model at 1/8o resolution.
Figure 9. Schematic illustrating the parameterization of stored surface water volume and inundated area
for the Putuligayuk River application
Figure 10. Average depth-area relationship (bold line) using RADARSAT observations of water extent in
individual grid cells (black diamonds), and surveyed depths for ponds (open circles).
Figure 11. Simulated total saturated extent versus time for the region including the Putuligayuk River
watershed from 1999 through 2000 (solid line); and total inundated area on specific dates as derived
from RADARSAT ScanSAR imagery (Bowling et al. 2002a) (black circles).
Figure 12. Sensitivity of simulated maximum annual flooded area (left) and discharge (right) to the new
model parameters of width fraction (top) and zmin (bottom). Annual values (1995-2000) are shown by
dashed lines, with the multiyear mean shown as solid black.
Figure 13. Simulated and observed SWE (top) and active layer depth (bottom) for the wetland study site
adjacent to the Putuligayuk catchment (1997-2000). .
Figure 14. Simulated and observed discharge for the Putuligayuk catchment, 1999-2001.