Estimating Snow Water Equivalent Using Snow Depth Data and Climate Classes
MATTHEW STURM,* BRIAN TARAS,
1
GLEN E. LISTON,
#
CHRIS DERKSEN,
@
TOBIAS JONAS,
&
AND JON LEA**
* U.S. Army Cold Regions Research and Engineering Laboratory, Ft. Wainwright, Alaska
1
Alaska Department of Fish and Game, Fairbanks, Alaska
#
Cooperative Institute for Research in the Atmosphere, Colorado State University, Fort Collins, Colorado
@
Climate Research Division, Environment Canada, Toronto, Ontario, Canada
&
WSL-Institute for Snow and Avalanche Research SLF, Davos, Switzerland
** National Resource Conservation Service, Oregon State Office, Portland, Oregon
(Manuscript received 13 August 2009, in final form 10 August 2010)
ABSTRACT
In many practical applications snow depth is known, but snow water equivalent (SWE) is needed as well.
Measuring SWE takes ;20 times as long as measuring depth, which in part is why depth measurements
outnumber SWE measurements worldwide. Here a method of estimating snow bulk density is presented and
then used to convert snow depth to SWE. The method is grounded in the fact that depth varies over a range
that is many times greater than that of bulk density. Consequently, estimates derived from measured depths
and modeled densities generally fall close to measured values of SWE. Knowledge of snow climate classes is
used to improve the accuracy of the estimation procedure. A statistical model based on a Bayesian analysis of
a set of 25 688 depth–density–SWE data collected in the United States, Canada, and Switzerland takes snow
depth, day of the year, and the climate class of snow at a selected location from which it produces a local bulk
density estimate. When converted to SWE and tested against two continental-scale datasets, 90% of the
computed SWE values fell within 68 cm of the measured values, with most estimates falling much closer.
1. Introduction
As global temperatures rise, the world’s snow resources
are predicted to change in significant ways (Hosaka et al.
2005; Christensen et al. 2007; Ra
¨
isa
¨
nen 2008; Deser et al.
2010). Long-term changes in global, regional, and local
snow depth (h
s
), snow water equivalent (SWE), and ex-
tent will ultimately have major ramifications for ecosys-
tem function, human utilization of snow resources, and
the climate itself through feedback mechanisms like snow
albedo (Barry 1996). Unfortunately, of the three snow
metrics listed above, only extent [i.e., snow cover area
(SCA)] is easily monitored using satellites. This moni-
toring, under way for several decades (Robinson 1993;
Frei and Gong 2005), has shown that global SCA has
been decreasing for the past 30 years (Lemke et al.
2007). SCA, however, is only an indirect measure of the
world’s snow water resources (Brown 2000; Brown et al.
2003). To fully understand global snow water trends, the
most fundamental metric to monitor is SWE, with depth
a close second . For monitor ing these metrics, two
potential methods are available: passive microwave
remote sensing and estimations based on direct mea-
surements. Many issues make estimating snow mass
from remote sensing problematic (Ko
¨
nig et al. 2001),
leaving field and station measurements currently the
primary means of inferring critical trends in snow re-
sources.
Of the two fundamental measurements, depth is
quicker and easier to measure in the field than SWE. No
detailed estimates of the total number of depth and
SWE measurements made worldwide is available, but
what is available suggests that considerably more depths
are collected than SWE measurements. For example,
Environment Canada operates 1556 snow depth sites,
but only 27 sites where SWE is measured. In the United
States, over 700 snow pillows (Beaumont 1965; Johnson
and Schaefer 2002), chiefly operated by the National
Resource Conservation Service (NRCS), are used to
make continuous SWE measurements, but conservative
estimates (based on sales of sonic sounders) suggest that
thousands of depth-monitoring stations are in operation.
Corresponding author address: Dr. Matthew Sturm, USA-
CRREL-Alaska, P.O. Box 35170, Ft. Wainwright, AK 99712.
E-mail: msturm@crrel.usace.army.mil
1380 JOURNAL OF HYDROMETEOROLOGY VOLUME 11
DOI: 10.1175/2010JHM1202.1
Ó 2010 American Meteorological Society
The U.S. Weather Service as well as the Swiss Service
measure more depths than water equivalents.
The two measurements are related but not inter-
changeable. The decision to make one type of mea-
surement versus the other (or a mix of both) is often
determined as much by tradition as an objective assess-
ment of accuracy, need, and available resources. To our
knowledge, no formal cost–benefit analyses of measuring
depth versus SWE have been conducted, but a rough es-
timate can be made of the relative time effort involved in
making each type of measurement. During field cam-
paigns we try to split our available time equally between
the two. A census of records for the past decade indicates
an approximate ratio of depth to SWE measurements of
about 30 to 1. In a more objective comparison, our stan-
dard field protocol calls for 201 depth measurements and
10 SWE measurements at each station. These are done in
the same amount of time giving a 20:1 ratio. Part of the
reason for this disparity is that we use an efficient self-
recording depth probe (U.S. Patent No. 5864059). In the
U.S. alone approximately 6000 man hours are expended
monthly during winter and spring in the collection of
;20 000 SWE data, but based on the 20:1 ratio, the same
work effort could produce nearly half a million depth
measurements.
What if there was a reliable method of converting
depth measurements into SWE? Such a method would
allow for a significant expansion in the number of loca-
tions where SWE could be estimated without incurring
much additional expense. This would improve our ability
to assess worldwide snow water resources. Here we de-
scribe a method that uses historical station data and
knowledge of snow climate classes (Sturm et al. 1995) to
estimate the bulk density of a snowpack. The bulk density
is then used to convert depth into SWE. Specifically, us-
ing a large (n 5 25 688) training set of depth, density, and
SWE measurements, we fit the data, including the climate
class of each datum (Sturm et al. 1995), with a nonlinear
analysis of covariance model (ANCOVA) using Bayes-
ian methods. From the analysis an equation is derived
that allows estimation of bulk density based on input of
location, snow depth, and the day of the year.Weexamine
the residuals to assess model fit and test the model against
a large independent dataset to provide the reader with an
understanding of model accuracy.
The method opens the possibility of converting thou-
sands of depth observations that have been, or currently
are being, collected in to SWE values. For example, it
could be used to convert Environmental Technical Ap-
plications Center monthly observed global snow depth
climatologies (Foster and Davy 1988) into SWE clima-
tologies. It is also likely to become even more useful in
the future if one promising new method of measuring
snow depth becomes operational: airborne LiDAR
(Hopkinson et al. 2004; Deems and Painter 2006; Deems
et al. 2006; Minoru and Hiroshi 2006; Schaffhauser et al.
2008). Airborne LiDAR snow surveys are conducted
much like gamma ray surveys (Peck et al. 1980) with a set
of measurements made prior to the start of the snow
season subtracted from a second set made after the snow
has accumulated. Unlike gamma ray surveys, however,
LiDAR surveys produce snow depth swath maps along
the flight path that can contain literally millions of depth
data. Accuracy depends on flight altitude, terrain, and
GPS navigation, but is typically about 615 cm in absolute
terms, with considerably higher accuracy possible when
frequent bench marks are included in the swath. Coupled
with the depth-to-SWE conversion method described
here, airborne LiDAR could greatly improve our ability
to estimate local-to-hemispheric snow resources.
2. Background
At a given point, snow depth (h
s
) is related to SWE by
the local bulk density (r
b
):
SWE 5 h
s
r
b
r
w
, (1)
where depth is measured in centimeters, density in
grams per centimeters cubed, r
w
is the density of water
(1 g cm
23
), and SWE is measured in centimeters of
water. It has long been observed (e.g., Dickinson and
Whitely 1972; Steppuhn 1976) that the natural range of
h
s
is many times greater than the range of r
b
. For a
northern Alaskan dataset (n 5 5323), 95% of all bulk
density values fell between 0.12 and 0.42 g cm
23
while
95% of all depths fell between 8 and 100 cm, a dynamic
range 4 times greater than that of density. For the
training dataset used in this paper the ratio is about 10:1.
Based on these ratios, estimating the more conservative
parameter (r
b
) while directly measuring the more dy-
namic (and easier to measure) parameter (h
s
) is the most
practical and potentially accurate method of estimating
SWE.
The form of Eq. (1), along with the large dynamic
range of h
s
, ensures that there will be a strong correla-
tion between SWE and depth, but it is mute about the
relationship between depth and bulk density. These two
parameters are functionally related, albeit in a complex
way. Steppuhn (1976) noted that Eq. (1) requires a co-
variance term between depth and density (C) when used
to compute mean areal values (indicted by overbars):
SWE 5 h
s
r
b
r
w
1 C. (2)
DECEMBER 2010 S T U R M E T A L . 1381
The covariance between depth and bulk density, how-
ever, tends to be negligible for a snow cover less than
80 cm deep (Pomeroy and Gray 1995), and at greater
depths is only about 2.5% of SWE (data from this study).
The error in neglecting C is therefore typically smaller
than other potential error sources in Eqs. (1) or (2), so
for simplicity, C is ignored.
Despite the fact that estimating r
b
using Eq. (1) makes it
possible to determine SWE from depth, there have been
relatively few formal attempts to do so. A handbook for
the NRCS field surveyors (Davis et al. 1970) provides
qualitative rules for adjusting the local mean snow den-
sity upward or downward depending on wind, exposure,
thaws, and the day of the year. Wilks and McKay (1996)
developed a power-law version of Eq. (1) with the bulk
density term (their ‘‘pseudodensity’’) varying with climate.
More recently, estimates of r
b
have been produced using
snow models that solve the surface energy balance and
evolve a snowpack layer by layer (Brun et al. 1989; Liston
and Elder 2006; Lehning et al. 2006), but these models
require extensive meteorological input as well as the cal-
culation of settlement of individual snow layers. They tend
to be difficult to apply to large basins and regions because
of their heavy data and computational demands.
3. Data and modeling
a. Data
The bulk density model was developed using a training
set of 25 668 records of snow depth, density, and SWE
from three countries and two continents (Table 1). These
data span the climate classes of seasonal snow defined by
Sturm et al. (1995): alpine, maritime, prairie, tundra, taiga,
and ephemeral. The classes are defined by the general
physical attributes of a snow cover (depth, density, type of
snow layers, etc.), attributes that tend to remain rela-
tively consistent within climate zones. The set combines
research surveys done in Canada, Switzerland, and Alaska
by the authors with data from monthly agency surveys
done by the NRCS on snow courses in the western United
States (Fig. 1a). Because few systematic measurements
have been made in ephemeral snow it had to be excluded
from the analysis. The data in the training set tend to be of
high quality because a high percentage were taken during
research campaigns.
At NRCS stations, data are collected on the first of the
month starting in November and continuing through
May or June each year. The Canadian and Alaskan data
were taken by the authors during over-snow traverses
conducted in March and April when conditions were
still well below freezing but the snowpack was at near-
maximum depth (Table 1). The Swiss data were taken at
TABLE 1. Sources of SWE, depth, and bulk density data for the training dataset.
Data source No. of data No. of sites Source
Alaska research studies 4977 10 traverses M. Sturm, USA-CRREL, unpublished data,
matthew.sturm@usace.army.mil
Canadian research studies 4934 4 sites, 2 traverses http://www.ccin.ca
C. Derksen, Environment Canada, unpublished data,
Chris.Derksen@ec.gc.ca
NRCS reporting stations 12 245 70 sites http://www.wcc.nrcs.usda.gov/snow/snowhist.html
Swiss research studies 3532 1 basin, various altitudes T. Jonas, SLF, Davos, unpublished data, jonas@slf.ch
Tot 25 688
FIG. 1. (a) The location of the sites in the training dataset
(Western Hemisphere only) used to develop the model. Swiss data
from a basin located at 478049N, 88439E are included in this set, and
(b) the location of the sites in the test dataset (Meteorological
Service of Canada 2000).
1382 JOURNAL OF HYDROMETEOROLOGY VOLUME 11
Alptal, Switzerland (Staehli and Gustafsson 2006) on
a weekly and biweekly schedule between November and
April. All data used in the analysis are available online
(see http://cdp.ucar.edu/cadis).
After developing the model, we located a large depth–
density–SWE dataset from Canada (i.e., the test dataset;
Fig. 1b; n 5 226 009) against which the model was tested.
This set, compiled in 2000 (Meteorological Service of
Canada 2000) and updated through 2004 by R. Brown,
comes from snow courses measured between 1935 and
2004 and is wholly different than the Canadian data used
in model development. The practices and equipment
used in collection of the large set varied widely and are
poorly documented. Considering this, and the age of the
data, they are probably not as accurate as the training set.
Additionally, the set is heavily biased toward sites clus-
tered along the U.S.–Canada border, many of which are
classified as maritime sites despite having markedly lower
depths than in the training set.
Both the training and test sets undoubtedly include
erroneous depth and SWE data (cf. Carroll 1995). SWE
measurements, in particular, are difficult to make accu-
rately. Using a corer that is either plunged or twisted into
the snow until it hits the ground, a snow core is removed
from the snowpack and weighed in the corer using
a spring balance calibrated to read out directly in SWE
units, or bagged and weighed later using a digital balance.
In either case, the measurement is gravimetric and re-
turns the bulk density of the snow; though in one case the
measurement is reported in SWE units. While the corer
is touching the ground, the depth is read off graduated
marks on the side of the corer with a nominal precision of
60.2 cm, but an accuracy that is closer to 61.0 cm.
At least four models of snow corer are in use in North
America [the Federal or Mt. Rose, the Adirondack, the
Eastern Snow Conference (ESC), and the Snow-Hydro]
and all suffer from similar problems, the most common
being a loss of snow. If the corer fails to retain a plug of
soil, vegetation, and/or ice, snow will fall out the bottom
of the corer when it is extracted from the snowpack
(Turcan and Loijens 1975). The corer can also hit an ice
layer, which will push snow out of the corer’s path, again
resulting in a light sample. We have tried to quantify the
magnitude of this undersampling by routinely comparing
core-based density values to those determined from
a density profile done in nearby snow pit. For measuring
the latter, a 3-cm high steel box cutter with a volume of
100 cm
3
is used with a digital balance. By integrating the
individual layer densities measured in the pit, a bulk
density can be computed and compared to the adjacent
values obtained by coring. Our experience suggests that
the pit measurements are more reliable and accurate,
albeit more labor-intensive, than the core measurements.
The results (spanning 8 yr, 5 field programs, and 3 types
of corers: Table 2) suggest core-based bulk densities av-
erage 7.1% lower than layer-integrated values.
Our results contradict prior studies (Work et al. 1965;
Peterson and Brown 1975; Goodison 1978; Farnes et al.
1982) that found that coring results in a high SWE bias.
This happens when excess snow is forced into a core tube,
a situation that is said to happen more frequently with
smaller diameter cutters (like the Federal Sampler with
across-sectionalareaof11vs30cm
2
for the ESC and
Snow-Hydro samplers), deeper, denser snow (Peterson
and Brown 1975), and poor cutter design (bad taper;
fewer and/or dull teeth). The type of snow (nature and
amount of depth hoar, number and thickness of ice layers)
also matters. SWE errors also arise when a spring balance
is used for weighing the cores (Bray 1973). Given our
findings, and those of the previous studies, we conclude
that within the training set (i) the SWE errors are small
and (ii) they are a function of the equipment used, the care
taken in collecting the measurement, and the nature of the
snowpack (none of which is routinely recorded). They are
effectively random in nature, and therefore cannot be
corrected easily, so no correction is applied.
b. Modeling
While it would be possible to develop a model relating
SWE directly to depth, it is more appropriate to model
bulk density (as detailed above) and use that relation-
ship to convert measured or assumed depths to SWE.
Bulk density is a complex function of snow depth (h
s
),
snow temperature (u), snow deposition history (t), and
the initial density of the individual snow layers (r
0
):
r
b
5 f (h
s
, t, u, r
0
), (3)
thus Eq. (1) can be written as
SWE 5 f (h
s
, t, u, r
0
)
h
s
r
w
. (4)
TABLE 2. Difference (%) between bulk density determined by
coring and obtained by integrating individual snow-layer densities.
The average value is weighted by the number of snowpits (n).
Yr Location n Error (%)
1996 Kuparuk basin, Alaska 37 10
2000 Ivotuk to Barrow, Alaska 60 23
2000 Barrow 14 25
2002 Ivotuk to Barrow, Alaska 26 29
2004 Nome to Barrow, Alaska 39 211
2005 Manitoba, Canada 101 210
2006 Manitoba, Canada 56 28
Tot 333 27.1
D
ECEMBER 2010 S T U R M E T A L . 1383
Without knowing explicitly the form of Eq. (3), i t is
clear that SWE is a complex and nonlinear function
of h
s
.
To model the bulk density we used Bayesian statisti-
cal methods to develop a nonlinear ANCOVA model.
Bayesian methods offer versatility when modeling com-
plex systems (Congdon 2003; Gelman et al. 2004; see also
http://www.bayesian.org/ and the journal Bayesian Anal-
ysis). When r
b
is plotted as a function of depth, it exhibits
nonlinearity in the mean response, a skewed error struc-
ture, and a nonco nstant error variance (heterosce dasticity) .
These complexities can be handled us ing Bayesian
methods, but are more problematic when using nonlinear
least squares. The Bayesian methods also produce useful
posterior probability distributions that define the uncer-
tainties associated with the model predictions, as well as
accounting for temporal and spatial dependencies. Three
easily obtained predictor variables: snow depth, DOY
(day of the year), and snow class, were used as input to the
model. The snow class, as discussed below, was intended
to capture environmental variables like temperature and
initial density, which are well known for impacting the
bulk density.
Many different functional forms of Eq. (3) (i.e., linear
models with and without quadratic and cubic functions of
depth and DOY; nonlinear functions) and all combina-
tions of predictor variables were tested. To choose which
form was best, as well as whether all three predictor
variables were needed, we used an objective method
of evaluation called the deviance information criterion
(DIC; Burnham and Anderson 2002; Spiegelhalter et al.
2002; Wheeler et al. 2010). DIC scores become lower (in
our case increasingly negative) as model accuracy goes
up, but if, and only if, the improvement in accuracy ex-
ceeds a penalty assessed for increasing model complex-
ity. A decrease in DIC score of 5–10 points implies a
significant improvement in model performance. The fi-
nal model that was selected had a DIC score ;1630
lower than its closest competitor, suggesting that it was
significantly better.
The variance in the bulk density data (when plotted
against snow depth, DOY, and subdivided into climate
classes) was modeled using both beta and normal distri-
butions, but in the end, a beta distribution was selected
because it could accommodate the skew in the data. In-
corporating a term to address the nonconstant variability
(heteroscedasticity) present in the training dataset im-
proved model performance (decreasing DIC by almost
2000 points) and was also incorporated into the final
model.
The software package WinBUGS (Lunn et al. 2000)
was used to run the various models and determine the
optimal result. Multiple Markov chains (using a range of
initial values) were run for more than 100 000 iterations
to simulate posterior distributions, with the first 50 000
iterations eliminated to allow for burn-in. We verified
convergence by examining run histories, autocorrelation
functions for simulated values, the Gelman–Rubin con-
vergence statistic (Gelman and Rubin 1992) as modified by
Brooks and Gelman (1998), and by performing additional
diagnostic tests using the CODA (Plummer et al. 2006)
libraries in Program R (R Development Core Team 2008).
Two potential types of errors affect the modeling. The
first relates to the use of snow climate classes as a proxy
for the physical processes that lead to the densification
of a snow cover. Snow bulk density is known to increase
with time and depth as the weight of the overlying snow
compacts underlying layers (Kojima 1966). It is also
moderated by temperature, solar radiation, the nature of
the initial snow deposit, infiltration of meltwater, and the
rate, sequence, and rapidity with which new snow layers
are added to the snow on the ground. Bilello (1969),
TABLE 3. Statistics for the training and test datasets with model errors also listed.
n
Depth
(cm)
Std dev
(cm)
Density
(g cm
23
)
Std dev
(g cm
23
)
SWE
(cm)
Std dev
(cm)
Density error
(g cm
23
) Std dev
SWE Error
(cm) Std dev
Training dataset
All data 25 688 108.4 111.3 0.312 0.093 38.7 48.7
Alpine 18% 130.0 82.5 0.335 0.086 46.2 35.5
Maritime 32% 176.6 149.9 0.343 0.101 68.4 67.1
Prairie 13% 88.5 72.8 0.312 0.085 30.0 31.4
Tundra 31% 43.8 24.7 0.284 0.075 12.5 8.2
Taiga 6% 59.6 22.7 0.217 0.056 12.8 6.0
Test dataset
All data 226 009 70.6 70.6 0.274 0.097 22.3 30.1 0.040 0.078 2.2 6.0
Alpine 24% 90.5 68.0 0.292 0.092 29.5 28.3 0.044 0.063 2.9 5.7
Maritime 47% 69.6 78.3 0.279 0.098 22.7 33.7 0.056 0.078 3.0 6.4
Prairie 7% 58.1 48.6 0.260 0.077 16.5 16.7 0.058 0.062 3.4 4.3
Tundra 6% 34.6 28.8 0.278 0.085 9.8 9.7 0.085 0.085 2.8 4.1
Taiga 16% 52.3 31.5 0.214 0.067 11.7 9.0 0.003 0.067 20.3 3.8
1384 JOURNAL OF HYDROMETEOROLOGY VOLUME 11
![FIG. 3. (a) SWE vs depth, (b) bulk density vs depth, and (c) SWE vs bulk density for the training dataset (n 5 25 688). In (a) and (c) a line has been fit to the data. In (b) an exponential curve [Eq. (4)] has been fit.](/figures/fig-3-a-swe-vs-depth-b-bulk-density-vs-depth-and-c-swe-vs-3pozuphp.webp)
![FIG. 8. Residual values [Eq. (7)] after applying the model to the training dataset.](/figures/fig-8-residual-values-eq-7-after-applying-the-model-to-the-32yg4t1y.webp)
