Vol.2, No.4, 392-405 (2011)
opyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
Agricultural Scienc es
Spatial and temporal variability of soil freeze-thaw
cycling across Southern Alberta, Canada
Andrew J. Phillips1,2*, Nathaniel K. Newlands1
1Environmental Health, Agriculture and Agri-Food Canada, Lethbridge Research Centre, Lethbridge, Canada;
*Corresponding Author: ajphilli@ucalgary.ca
2Department of Geomatics Engineering, Schulich School of Engineering, University of Calgary, Calgary, Canada.
Received 9 September 2011; revised 21 October 2011; accepted 30 October 2011.
Soil freeze-thaw cycles play an important role in
all aspects of agro-ecosystems, such as crop
productivity, the evolution of the soil matrix,
including trace-gas emissions. In regions that
experience synoptic weather conditions through-
out the winter, freeze-thaw cycles generally oc-
cur in one of two categories; seasonal or winter
cycles. Current soil vegetation atmosphere mod-
els (SVAT’s) often include a heat-transport soil
freeze-thaw algorithm, but lack detail on complex
interactions between the main driving variables.
Boundary conditions for these models are often
based only on a few climate variables and typi-
cally lack regional context. A nested statistical
analysis was applied to identify the optimal set
of environmental variables (via a stepwise re-
gression selection procedure) to track soil
freeze-thaw dynamics. Historical data collected
between the years 2006-2009, for 17 long-term
climate stations distributed across southern
Alberta Canada was utilized. Cross-correlation
between wind speed and maximum air tem-
perature identified Chinook-driven freeze-thaw
events, with such interaction varying signifi-
cantly across the region and by soil depth. Cli-
mate-soil interactions were most significant
predictors of soil temperature during winter
months. The seasonal freeze-thaw cycle is es-
timated to vary between 112 - 131 days, consist-
ing of 12 - 20 winter cycles (1 cm depth), and 1-5
winter cycles (5 cm depth) with average lag time
of 26 - 112 days. Freeze-thaw prediction was
greatly improved when higher-order climate in-
teraction terms were considered. Our findings
highlight the importance of soil-water assume-
ptions within complex ecosystem SVAT models
in resolving regional-driven climatic trends. Alon g -
side improved representation of regional trends
aimed at reducing model-based uncertainty,
such efforts are expected to, in tandem, help
advance the geostatistical design, and imple-
mentation of agroenvironmental monitoring sys-
tems that combine in-situ and satellite/remote-
sensing derived estimates of near-surface soil
Keywords: Freeze-Thaw; Soil Temperature;
Agro-Ecosystem Modeling; Regional Climate;
Soil Science
Soil freeze-thaw cycles affect plant growth, greenhouse
gas emissions, soil structure, stability, microbe popula-
tions and nutrients. Freeze-thaw cycles have been shown
to increase grassland above-ground productivity, while
causing root damage [1]. Significant increases in nitrous
oxide (N2O) emissions, a major greenhouse gas that con-
tributes to global warming from agricultural systems, has
been documented during thawing events, with the great-
est emissions occurring during the spring thaw [2,3].
Freeze-thaw variability regulates the cumulative release
of N2O by controlling how water, in its different states,
physically blocks emission from the soil surface. At
depth, soil temperature, that regulates biological micro-
bial activity, is also affected by freeze-thaw variability.
Freeze-thaw cycles reduce soil aggregate size and make
soils more susceptible to erosion [4]. There is evidence
that soil microorganism populations adapt in order to
survive in regions that experience multiple winter freeze-
thaw cycles [5]. Recent findings indicate that climate
change may increase the frequency of freeze-thaw cycles
[6]. Alongside anticipated changes in the frequency of
freeze-thaw events, the amplitude or strength of influ-
ence is also expected to increase, given the longer ab-
sence of insulating snow cover (i.e., longer growing sea-
sons at polar latitudes) that accompanies a rise in global
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
mean air temperature due to global warming [7]. A better
understanding of freeze-thaw cycling behavior in rela-
tion to climate change, and its influence on soil condi-
tions can improve adaptive, multi-objective decision-
making that includes crop selection, use of herbicides,
irrigation and nutrient application timing, as well as
leaching and erosion management. There are a number of
soil-vegetation-atmosphere models (SVAT’s) that have
been developed to simulate soil water flow, soil tem-
perature and soil freezing and thawing. These models are
often driven by boundary conditions based on a few cli-
mate parameters and are not parameterized with a re-
gional-specific context using statistical relationships be-
tween observed data and relevant related variables.
Models that do not consider regional-specific aspects of
climate interactions may produce unrealistic model out-
put, underestimating the impact of freeze-thaw dynamics
on agricultural land and its productivity. Along with
modeling accuracy improvements, related advances in
rapid automated monitoring of soil and climate at fine
spatial scales and the integration of multiple sensor re-
sources may reduce measurement uncertainty.
Several existing SVAT models consider soil freeze-
thaw cycling, such as the Daily Century Model (Daycent)
[8], Ecosystem-Atmosphere Simulation Scheme (EASS)
[9], Simultaneous Heat And Water (SHAW) [10,11], and
HYDRUS [12] (Table 1). All these models aim to model
heat transport through a one-dimensional soil profile.
SHAW, HYDRUS and EASS are stochastic, while Day-
cent is deterministic. These models are all based on the
Table 1. A summary of leading parameters, assumptions and corresponding upper boundary conditions for current soil water freeze-
thaw algorithms.
Ecosystem Model Type Key Parameters Assumptions Upper Boundary Conditions
simulation scheme
Latent heat of fusion,
equilibrium total water
potential, ice potential,
absolute temperature
1-D soil layer model,
explicit thermal
separation of
vegetation & soil,
water & energy
balance coupled,
11, 0,1
,0 1
uns wspp
Iz l
 
11, 0,1
= soil moisture flow fluxes between
layers l 1 and l; l and l + 1, re-
Pun = precipitation rate at ground surface
Es = evaporation from soil surface
ρw = density of water
Ps = percolate water supplied to underlying soil
Ip = potential infiltration rate before ponding
Ilim = limiting value to infiltration rate after
zp = depth of surface pond
Daily Century-Model
(DAYCENT) Deterministic
Diurnal temperature extremes,
depth, thermal diffusivity,
time step,
1-D, soil thermal
diffusion constant
with depth, diffusion
with no advection,
sinusoidal pattern of
soil surface
Rd = diurnal temperature range at depth
Rs = diurnal temperature range at soil surface
d = depth
α = thermal diffusivity
p = time step
HYDRUS 1D Stochastic
Apparent thermal conductivity of
soil, temperature, time, depth,
coefficient of apparent thermal
conductivity, apparent thermal
conductivity of water,
oxygen uptake rate
1-D, convection &
conduction heat
transport, coupled
energy, vapour
& water transport
sin 12
 
T0 = temperature at soil surface
T = average temperature at soil surface
A = amplitude of temperature sine wave
pt = time period of temperature cycle
Heat And Water
depth, soil thermal conductivity,
soil temperature, density of water
& ice, specific heat capacity of
water & ice, liquid vapour flux,
total solutes present per mass of
soil, volumetric heat capacity of
soil, time, density of ice, latent heat
of fusion, soil vapour density, latent
heat of vaporization, vapour flux
1-D layered system,
convective heat
transfer by liquid
and latent heat
transfer by vapour,
layered system
RHLEG 
H = sensible heat flux
Rn = net all-wave radiation
E = total evapotranspiration
Lv = latent heat flux
G = soil heat flux
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
physical properties of heat and water transport through
porous media. Cycling is represented through upper
boundary conditions that use climate information to
drive soil freeze-thaw cycling of the underlying soil lay-
ers. HYDRUS uses mean temperature and Daycent uses
diurnal temperature ranges in the boundary conditions as
the main climate drivers of soil freeze-thaw dynamics.
SHAW uses net-solar radiation balance and EASS uses
precipitation and infiltration rates as main climate bound-
ary condition driving variables. All of these models as-
sume that one or two climate variables are the main
drivers of freezing and thawing in the soil and other cli-
mate variable effects are negligible. The problems with
current soil freeze-thaw models is that they do not con-
sider interactions between climate variables and the
boundary conditions have little regional context to drive
the models.
Freeze-thaw cycles can be divided into two categories:
seasonal cycles and winter cycles. Seasonal cycles are
characterized by initial freezing of the soil water in the
late fall or early winter when maximum and minimum
air temperatures stay below 0˚C consistently, and like-
wise, by thawing when soil temperature stays above
freezing. The soil profile freezes from the top down, the
depth of freeze and rate are affected by the duration and
intensity of the cold weather as well as the insulation
capabilities of the soil and snow cover. The thawing por-
tion of a cycle begins when the air temperatures increase
to above 0˚C for long enough to thaw the soil water.
Thawing also occurs from the soil surface down; and if
the warm weather persists long enough, the entire frozen
portion will thaw, otherwise a return of cool weather
conditions will re-freeze the upper portion of the soil
profile. Such warm weather allows for partial thawing of
the frozen soil profile or even a complete thawing de-
pending of the extent of the initial freeze, duration and
intensity of the warm weather. While the physical state
of soil water is determined by soil temperature, other
variables such as air temperature, soil moisture content,
snow cover, precipitation and wind alter soil tempera-
Winter freeze-thaw cycles are events driven by synop-
tic weather conditions that bring unseasonably warm tem-
peratures to a region during a normally cold winter. Per-
sistent wind contributes greatly to the intensity and fre-
quency of winter freeze-thaw cycles. In Canada, such
wind-driven freeze-thaw cycles are most pronounced as
Chinooks that occur in Canada almost exclusively within
southern and central Alberta. Chinook winds belong to a
family of mountain winds that blow in regions where
long mountain chains are more or less perpendicular to
the prevailing wind. Chinooks are established when
moist Pacific air is blown westward by prevailing winds
[13]. As the air mass cools, it then rises over the Rocky
Mountains of British Columbia. The moist air mass con-
denses and precipitates over the mountain ranges. The
drier Pacific air mass reaches the leeward side of the
mountains and warms at the dry adiabatic lapse rate as it
descends onto the prairies with high wind speeds [13].
The warm moisture-deprived air has the ability to quickly
melt and evaporate any snow cover, while thawing the
upper layers of the soil. Temperatures during a Chinook
event can increase air temperatures by more than 20˚C in
a few hours [14]. A Chinook can last only a matter of
hours to more than 24 hours. After the Chinook dissi-
pates the region quickly returns to normal winter tem-
peratures. Climate-land interaction considerations pro-
vide an explanation of the occurrence probability of
dominant winter freeze-thaw events.
Freeze-thaw cycling is also influenced by long-term
climate teleconnections, as large-scale atmospheric and
oceanic circulations, such as the El Nino/Southern Os-
cillation (ENSO), Pacific Decadal Oscillation (PDO),
Pacific North America (PNA), and North Pacific (NP)
pattern. These large scale circulations have been shown
to account for a large portion of cold season temperature
trends in North America [15]. Teleconnections have an
influence on air temperature and fresh water freeze and
break-up dates in Canada [16]. Air temperature is a ma-
jor driving force of northern latitude freeze and thaw
cycles [17]. Teleconnections influencing seasonal air
temperature trends in North America will, in turn, also
have a driving effect on seasonal soil freeze-thaw cycles.
Nonetheless, the influence of teleconnections on re-
gional soil temperatures are difficult to quantify without
knowledge about details on timing of the events associ-
ated with regional impacts.
In this paper, we investigate the influence of climate-
soil interaction/s on the intensity and frequency of soil
freeze-thaw dynamics in time and space. For our study
region (southern Alberta, Western Canada) we apply a
nested, multivariate statistical modeling approach to
identify the best set of variables to predict freeze-thaw
cycling based on available data. In this way, our aim was
to devise statistical models that capture the influence of
interactions between climate, soil and other relevant en-
vironmental variables for improved tracking of freeze-
thaw variability in time and space. Our statistical predic-
tions of the timing, intensity, and inter-annual variability
of soil freeze-thaw cycling are important for improving
the reliability of SVAT models, the identification of most
significant driving variables and interaction terms. Also,
we obtain best-fit models that provide surface soil tem-
perature boundary conditions for soil water and solute
vertical transport models. Our findings illustrate the im-
portance of regional characterization of freeze-thaw
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
variability. We obtain correlation maps for identifying
where additional climate stations might be located to
improve the monitoring of freeze-thaw activity and its
potential impact on growing-season and winter cover-
crop productivity as well as agro-ecosystem resilience.
Many of the world’s most productive agricultural re-
gions are situated near mountainous areas. The magni-
tude and frequency of freeze-thaw cycling is anticipated
to increase as a result of future global projected climate
changes—yet, its potential regional variability is not
well understood [7]. For this reason, better tracking of
freeze-thaw activity in time and space, in relation to en-
vironmental variables, contributes to reducing obser-
vational and modeling uncertainty, enhancing the ability
of decision-makers to identify and gauge potential vari-
ability in freeze-thaw cycling. In turn, agricultural deci-
sion-makers may better assess both short- and long-term
risks to agricultural crops, in relation to freeze-thaw ac-
tivity and management practices for the use-efficiency of
irrigation, fertilizers and herbicides (i.e., amount and
timing of applications).
2.1. Data
The region of interest for our study was Southern Al-
berta, Canada, which experiences winter freeze-thaw
cycles due to warming effects of the Chinook winds, and
longer term seasonal freeze-thaw cycles that are strongly
influenced by teleconnections. Historical maximum,
minimum, mean air temperatures, wind speed and pre-
cipitation for 17 climate stations distributed across the
study region were obtained (Figure 1). These stations
recorded soil moisture and soil temperature at various
soil depths. A seasonal time series was obtained for each
station from 1 May 2006 to 31 May 2009 inclusive over
three winters (Figure 2). Station elevations ranged from
767 m to 1310 m above mean sea level. The climate sta-
tions are operated by Agriculture Drought Monitoring
(AGDM) and Agriculture and Agri-Food Canada (AAFC).
The datasets were inspected for errors, involving the
conversion of hourly to daily measurements for 4 sta-
tions (AAFC stations). Day lengths for each of the sta-
tions was computed using United States Naval Observa-
tory (USNO) day length calculator and the latitude and
longitude of each station.
Table 2 provides a brief summary of the data obtained
for each station. Due the coarse spatial resolution of the
climate stations, before considering more advanced geo-
spatial models, a statistical analysis was first under-
taken to identify leading variables, while awaiting more
fine scale data from the study region.
Figure 1. Long term climate monitoring stations distributed
within study region of Southern Alberta, Canada (n = 17).
Table 2. Summary of the data collected for each climate station in southern Alberta (n = 17).
Climate Stations Data Interval Variables Data Source
Barnwell, Bodo, Brocket, Champion,
Del Bonita, Foremost, Hussar, Morrin,
Olds, Oyen, Schuler, Stettler, Wrentham
Min, Max & Mean air temperature,
Soil Temperature & Moisture at
5, 20, 50, 100 cm, 10 m Wind speed
Agriculture Drought
Monitoring Station (AGDM)1
Lacombe, Lethbridge,
Onefour, Stavely hourly
Min, Max & Mean air temperature
Soil Temperature & Moisture at
5, 20, 50, 100 cm 10 m Wind speed
Agriculture Agri-Food
Canada (AAFC)2
All Climate stations daily Day Length United States Naval
Observatory (USNO)3
1Available at: www.agric.gov.ab.ca/app116/stationview.jsp; 2Available internally through AAFC; 3Available at: aa.usno.navy.mil/data/docs/RS_OneYear.php.
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
Figure 2. Average daily climate and soil conditions from May 1, 2006 to May 31, 2009, data was collected from
long term climate stations (n = 17) in southern Alberta. The 95% confidence interval is shown for maximum, mini-
mum air temperature, soil temperature at 5 cm, soil moisture at 5 cm and day length.
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. http://www.scirp.org/journal/AS/Openly accessible at
2.2. Methods
To identify the best or more representative statistical
model to capture freeze-thaw observed historical vari-
ability, we evaluated in a nested fashion, a wide array of
possible models comprising different explanatory, pre-
dictor variables and assumed interactions between them
(i.e., model and variable selection step). The daily mean
of each variable was calculated from the 17 station data-
sets. Generalized linear models (GLM) assuming a
Gaussian distribution for each random variable were
applied to the dataset to find the best explanatory (de-
pendent) variable and set of predictor variables (inde-
pendent). Soil temperature and moisture at 5 cm, accu-
mulated precipitation, mean air temperature, wind speed,
and day length were each tested as a dependent variable
in a GLM. The remaining variables were tested as inde-
pendent variables in the GLM, with no variable interac-
tions being considered. Based on the lowest Akaike In-
formation Criterion (AIC) and resultant residual plots,
soil temperature at 5 cm was selected as the best ex-
planatory variable (Table 3).
A linear regression model was constructed with soil
temperature at 5 cm as the explanatory variable. A nested
or step-wise model selection was employed that simpli-
fies a linear model by reducing the AIC to find the best
independent variables and coefficients. Mean air tem-
perature, soil moisture at 5 cm, accumulated precipita-
tion, wind speed, and day length (both with and without
interactions) were considered as predictors for soil sur-
face temperature (i.e., surface depth of 5 cm). Examina-
tion of the residual plots showed unsatisfactory results,
so a second stepwise model selection was performed
next using minimum and maximum air temperature in-
stead of the mean. The new selected independent vari-
ables were minimum, maximum air temperature, soil
moisture at 5 cm, accumulated precipitation, wind speed,
and day length, with and without interactions between
The values for each climate variable from November
1 to April 1 of each winter were extracted from the mean
seasonal dataset. Figure 3 shows the winter dataset for
three years (i.e., 2006-2007, 2007-2008, 2008-2009). A
stepwise linear model selection strategy was also employed
to select the best set of predictor and interaction terms
during the winter season for soil temperature at 5 cm. The
predictors tested were maximum, minimum air tempera-
ture, soil moisture at 5 cm, accumulated precipitation,
wind speed, and day length as well as their interaction.
To examine climate-soil interaction variability across
the study region, the statistical models were fitted to data
collected from all stations. Geostatistical spline-based
interpolation was then performed to interpolate between
stations and to generate a representative “correlation map”
of freeze-thaw variability across the study region. The
correlation between wind speed and maximum tempera-
ture was most representative of freeze-thaw activity and
was interpolated from best-fit predictions obtained at
each station location. High correlations between wind
speed and maximum air temperature were linked to
Chinook-driven freeze-thaw events, such that locations
with high correlation values between wind speed and
maximum air temperature would also experience more
freeze-thaw cycles relative to locations with lower cor-
relation values. Spline surfaces were created from the
potted correlation values using varied degrees of smooth-
ing. Using a spline surface smoothing tension of 0.7, a
leave-out-one validation was performed. Each station
was individually removed from the spline-surface (geo-
statistical) calculation, and the predicted correlation at
the location of each dropped or withheld climate station
was then recorded. Model cross-validation error of the
withheld correlation was calculated as the absolute dif-
ference between the actual correlation value and the
withheld correlation value. The leave-out-one validation
errors were then ranked in descending order to highlight
the stations with the greatest and least influence on the
correlation surface.
Table 3. Summary of the criterion used to judge the best dependent variable in a GLM, all models had 1121 residual deviance de-
grees of freedom and 1126 null deviance degrees of freedom.
Selected dependent variable (y),
n = 4 AIC BIC
y = Precipitation 4872 4907 4913 5860
y = Mean Air Temperature 6205 6240 16,038 137,995
y = Soil Moisture at 5 cm 9500 9535 298,494 440,386
y = Soil Temperature at 5 cm 5386 5421 7754 97,091
y = Wind Speed at 10 m 6709 6745 25,093 30,854
y = Day Length 3914 3949 2101.1 9537
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
Figure 3. The soil temperatures (black) at 5 cm depth for the 2006/2007, 2007/2008 and 2008/2009
winters (November 1 to April 1). Each winter is labeled and separated by a vertical gridline.
In order to understand the temporal variability be-
tween the three winters, the average soil freeze- thaw
characteristics across southern Alberta were com- puted
for each winter from November 1 to April 1. From this
data, a set of best-fit statistics was calculated for each
winter, this included the mean maximum, minimum air
temperature, mean soil temperature at 5 cm, total soil
freezing days, freeze-thaw cycles at 5 cm, estimated
freeze-thaw cycles at 1cm, the mean magnitude of soil
freeze, the maximum magnitude and day of soil freeze at
5 cm.
The number of soil freezing days, day of first freeze
and final thaw, length of seasonal freeze-thaw cycle,
mean lag between 5 cm freeze-thaw cycles and accumu-
lated precipitation were also computed to better charac-
terize inter-annual variation in freeze-thaw cycling be-
havior. The length of each winter’s seasonal freeze-thaw
cycle was computed by subtracting the time between the
initial freeze and final thaw of each winter. The mean
magnitude of soil freeze was computed as the average of
the negative 5 cm soil temperatures. The number of soil
freezing days each winter was the total of days each
winter where the average daily soil temperature at 5 cm
was below 0˚C. Soil freeze-thaw cycles at 5 cm depth
were calculated by counting the number of times the soil
temperature at 5 cm dropped below 0˚C then rose above
freezing again. The soil temperature at 5 cm was the
shallowest soil measurement available and therefore
shallow freeze-thaw cycles caused by daily temperature
extremes were not captured at this measurement depth.
Freeze-thaw cycles at the 1 cm soil depth were estimated
from daily maximum and minimum air temperatures.
One freeze-thaw cycle was defined by whenever the
daily maximum air temperature is greater than or equal
to 6.3˚C and the daily minimum is less than or equal to
–3.5˚C [18].
The seasonal time series for soil temperature at 5 cm
was filtered using Fourier domain smoothing and de-
noising with two different threshold frequencies (0.1 and
0.05) to remove the high-frequency (i.e., noise) variabil-
ity. The low frequency component of the filtered sea-
sonal time series represents the signal as the seasonal
trend in soil temperatures at 5 cm. Each filtered time
series of soil temperature was then fitted with a two
damped sinusoidal model. An average winter soil tem-
perature time series was created by averaging the 5 cm
soil temperature for each day between November 1 and
April 1 for the 2006/2007, 2007/2008 and 2008/2009
winters. The average winter 5 cm soil temperature time
series was fitted with a B-spline curve tracking winter
freeze-thaw variability.
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
Mean daily air temperature alone was unable to ade-
quately describe observed trends in near surface-soil tem-
perature, as Chinook events disrupt climate-soil interac-
tion via rapid increases in diurnal temperature change.
For this reason, minimum and maximum air temperature
better tracked temperature range and freeze-thaw cycling.
Figure 4 shows the distribution of residuals for the best-
fitting seasonal and winter models. The best fitting mod-
els for the seasonal and winter datasets both used maxi-
mum, minimum air temperature, soil moisture at 5 cm,
accumulated precipitation and wind speed at 10 m as
predictors for soil temperature at 5 cm. Best-fitting
models were obtained when interactions between inde-
pendent variables were considered (Table 4).
Figure 4. Statistical distribution quantiles (5%, 95%) and pattern of residuals for the seasonal model (left) and the winter model
(right). The standardized residuals lie within two standard deviations of the observations on the scale-location plot of both models.
Both the seasonal and winter models loose accuracy at extremes as shown in the quantile-quantile plots.
Table 4. Soil temperature predicted over a seasonal and winter time series, with and without interaction term considered, where p <
0.001 for all models. DF = Degrees of freedom.
Model-data partition Interactions AIC F-stat Terms DF Adjusted R²
seasonal no 5272 2404 6 1120 0.928
yes 4308 651 57 1169 0.971
winter no 1560 346 6 450 0.819
yes 1446 55 58 398 0.873
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
Nested/stepwise model selection for the winter dataset
did not improve the prediction skill of the winter model
with interactions as much as the seasonal model with
interactions. While more interactions were considered in
the best fit winter model, the greater number of predictor
terms in the best-fitting winter model is likely due to
stronger climate fluctuations during winter. The quan-
tile-quantile plot in Figure 4 shows that both the sea-
sonal and winter models need improved prediction ac-
curacy at soil temperature extremes. Eq.1 represents the
reduced form of the best fitting statistical model for the
seasonal time series. All predictor terms of soil tempera-
ture in Eq.1 had a p < 0.001. Eq.1 can be used to esti-
mate the soil temperature at 5 cm in the study region.
Eq.2 estimates winter soil temperatures between No-
vember 1 and April 1 in the study region and represents
the reduced form of the best fitting statistical model for
the winter time series. All the terms used for predicting
soil temperature in Eq.2 had a p < 0.01. Eqs.1 and 2 are
regional-specific and therefore capture complex interac-
tions between climate variables that drive soil tempera-
ture throughout the year and during the winter season.
max min
minmax min
5.4 0.23430.0474
 
max minmaxmin
0.1237 0.1877
Ts = soil temperature at 5 cm
Tmax = daily max temperature
Tmin =daily min temperature
M = soil moisture at 5 cm
W =wind speed at 10 m
P =precipitation
D =day length
Maximum and minimum air temperature, as individ-
ual predictors of soil temperature in both the winter and
seasonal models, were significant (p < 0.01), with their
interaction also significant (winter p < 0.001, seasonal p
< 0.05). The higher significance level of climate-soil
interactions during winter indicates that diurnal tem-
perature range is critical for predicting soil temperature
during the winter season. In the seasonal model the only
significant predictors of soil temperature were: day
length (p < 0.001), soil moisture (p < 0.01) and precipi-
tation (p < 0.05). The high significance level of day
length in the seasonal model was expected because in
Northern mid-latitudes the length of day has been pre-
viously shown to drive seasonal climate variation. Southern
Alberta receives very little winter precipitation; the av-
erage accumulated precipitation for the region was 138
mm in total over the three winters. Whereas the average
accumulated precipitation between May 1, 2006 and
May 31, 2009 for the study region was 1078 mm. In the
winter model, precipitation did not significantly (p > 0.1)
contribute to the prediction of soil temperature, but pre-
cipitation was significant (p < 0.05) in predicting soil
temperature in the seasonal model. Soil moisture was
only a significant (p < 0.01) predictor of soil temperature
in the seasonal model. Likely precipitation and soil mois-
ture individually play a more significant role in driving
soil temperature during the spring and summer when
precipitation events are more frequent and soil water is
unfrozen. Despite the low amount of precipitation during
the winter, precipitation and soil moisture still had sig-
nificant predictive ability for soil temperature when in-
teracting with other climate variables in the winter model.
Interacting terms were significant predictors of soil
temperature in the winter model, indicating that winter
soil temperatures are driven by complex interactions
between climate variables. A notable significant interac-
tion term in the winter model was maximum air tem-
perature, minimum air temperature and wind speed (p <
0.05), this interaction had a positive influence on soil
temperature. The maximum, minimum air temperature
and wind speed interaction may be capturing the Chi-
nook events; these events can increase daily maximum,
minimum air temperatures associated with high wind
speed [13]. The specific interaction terms that may be
related to long term synoptic climate events, such as
teleconnections nonetheless require more investigation
as the complex nature of how teleconnections influence
climate and freeze-thaw events over large areas.
Between November 1 and April 1, wind speed and soil
moisture were strongly correlated across all stations,
with correlation ranging between 0.81 at Barnwell and
Stavely to 0.96 at Oyen, with a mean of 0.88. The strong
winter co-occurrence of soil moisture and wind speed in
southern Alberta is best explained as driven by Chinook
events. Chinooks bring warm-dry-high-speed winds that
have the ability to melt snow and thaw the upper layers
of soil [14], snow melt that is not evaporated can infil-
trate into the thawed upper soil layers. Maximum air
temperature and wind speed were not as strongly corre-
lated as wind speed and soil moisture during the winter,
0.66 at Onefour to 0.89 at Oyen with an average correla-
tion of 0.71 (Table 5). Figure 5 shows significant spatial
variability in the correlation between maximum air tem-
perature and wind speed. The various fitted spline sur-
faces of this correlation pattern of climate-soil variability
reveal a similar spatial pattern to the observed regional
pattern of Alberta’s Chinook signal from 1951-1990 [13].
The similar spatial pattern between Alberta’s average
Chinook signal and the maximum air temperature and
wind speed correlation surface supports the idea that a
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
Table 5. Best-fit model cross-validation (leave-one-out) statistical results: Maximum air temperature and wind speed correlations at
each station before and after each station is withheld from the spatial interpolation, with the absolute error of predicted correlation
value error and error rank of each station. The maximum air temperature and wind speed correlation values were used as an indicator
potential freeze-thaw cycling due the relationship between Chinook events, unseasonably warm temperatures and wind speed.
Correlation of wind speed and
max air temperature
Predicted correlation at
withheld station
Absolute error of predicted
Error rank of predicted
Barnwell 0.70 0.700 0.000 17
Bodo 0.72 0.924 0.204 1
Brocket 0.76 0.778 0.018 14
Champion 0.70 0.744 0.044 6
Del Bonita 0.76 0.718 0.042 9
Foremost 0.66 0.602 0.058 5
Hussar 0.66 0.704 0.044 7
Lacombe 0.68 0.681 0.001 16
Lethbridge 0.74 0.734 0.007 15
Morrin 0.71 0.667 0.043 8
Olds 0.69 0.708 0.018 13
Onefour 0.60 0.636 0.036 10
Oyen 0.89 0.707 0.183 2
Schuler 0.71 0.771 0.067 3
Stavely 0.76 0.727 0.033 12
Stettler 0.68 0.713 0.033 11
Wrentham 0.64 0.699 0.057 4
Mean error: 0.052
Figure 5. Predicted strength of spatial occurrence of freeze-
thaw, based on station monitoring data (n = 17) for maximum
air temperature and wind speed (10 m) using thin-plate spline
interpolation under varying degrees of smoothing.
Chinook event is a complex interaction between two or
more climate variables. Locations with high correlation
values between wind speed and maximum air tempera-
ture were expected to have an increased occurrence of
Chinook events and freeze-thaw cycles. Tab le 5 shows
the results of the leave-one-out cross-validation, the pre-
dicted wind speed and maximum temperature correlation
value at the location of each withheld station was re-
corded, the absolute difference between observed corre-
lation value and the predicted values were calculated and
ranked. Barnwell, Lethbridge and Lacombe had the low-
est absolute difference between observed and predicted
correlation values, this can be explained by these sta-
tions’ relative close proximity to surrounding stations
with similar wind speed and maximum temperature cor-
relation values (Tabl e 5 ). The results of the leave-one-
out validation show that Oyen, Bodo and Schuler had the
largest errors (Ta b l e 5) and therefore were the most in-
fluential stations in calculation of the spline surfaces
(Figure 6). Oyen, Bodo and Schuler are all in the most
eastern part of the study region with a large spatial gap
between these monitoring stations and the nearest sta-
tions to the west. For example, there is a larger variation
in correlation pattern at Hussar (0.66) and Oyen (0.89).
In order to better track climatic influences on agro-eco-
systems, additional monitoring stations likely need to be
placed near such locations where higher observed corre-
lation between climate and soil interaction is detected.
The area between Hussar and to the east Oyen, Bodo and
Schuler is a prime candidate for an additional climate
monitoring station that would improve the tracking of
Chinook related climate interactions in eastern Alberta.
The winter of 2006/2007 was the mildest of the three
winters; it had the warmest mean soil temperature, mean
maximum and minimum air temperatures (Table 6). The
average and the maximum magnitude of soil freeze were
warmer than the following winters. Consistently warmer
air temperatures would have led to repetitive thawing
and refreezing of the upper 1 cm of the soil leading to an
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
Figure 6. The results of the leave-one-validation for the three
stations with the largest absolute error of the predicted strength
of freeze-thaw spatial occurrence (spline tension = 0.7).
estimated 18 freeze-thaw cycles at 1 cm soil depth. One
explanation for there being only 4 freeze-thaw cycles at
5 cm depth is soil temperature at 1 cm depth will be sen-
sitive to daily fluctuations in air temperature, whereas
soil at 5 cm depth will need prolonged exposure to warm
temperatures to thaw [19]. The difference between the
length of the seasonal freeze-thaw cycle and soil freez-
ing days within the cycle show that when the soil thawed
at 5 cm depth it remained unfrozen for multiple days.
The 2006/2007 winter received the second most winter
precipitation that may have provided some insulation to
the soil [20] preventing the soil from freezing to the 5
cm depth. The 2006/2007 winter had the earliest date of
initial soil freeze and final thaw at 5 cm, cold tempera-
tures may have come early with warmer temperatures
coming earlier as well. The date of maximum magnitude
of soil freeze being in mid-January instead of later in the
month provides evidence of an earlier winter. The winter
of 2007/2008 had the lowest amount of accumulated
winter precipitation; this may have contributed to the
low average soil temperature, the coldest minimum soil
temperature and the low number of freeze-thaw cycles at
5 cm depth (Tab le 6). The lack of winter precipitation
may have led to minimal snow cover to insulate the soil.
Without much snow cover the soil surface would be
more susceptible to air temperature extremes [20]. The
2007/2008 winter experienced the largest range between
the mean minimum air temperature and the mean maxi-
mum air temperature; this would explain the larger num-
ber of 1 cm depth freeze-thaw cycles when compared to
the 2008/2009 winter. A contributing factor to the lack of
freeze-thaw cycles at 5 cm depth during the 2007/2008
winter may be that the positive temperatures lacked per-
sistence and intensity to thaw the soil much below 1 cm
depth, leading to an estimated 20 cycles near the surface
and yet only 1 at the 5 cm depth. An interesting aspect of
the 2007/2008 winter is that the soil freezing days and
length of the seasonal cycle are the same, supporting the
idea that the warm weather was never able to thaw the
soil to 5 cm. The date of coldest soil temperature, initial
soil freeze and final thaw indicate that the 2007/2008
winter was shorter, with a late start and earlier end than
the preceding and following winters. The coldest of the
winters was 2008/2009; the mean maximum and mini-
mum air temperatures were the lowest of the three win-
ters. The 2008/2009 winter had a warmer mean soil tem-
perature, mean and maximum magnitude of soil freeze at
Table 6. The average climate variability of the southern Albertan winters (November 1 to April 1) for each winter. Each reported
variable is an average value computed using all the southern Alberta climates stations (n = 17).
Variables related to freeze-thaw 2006/2007 2007/2008 2008/2009
mean max air temperature (˚C) 1.302 0.533 –0.716
mean min air temperature (˚C) –9.923 –11.557 –12.172
mean soil temperature at 5 cm depth (˚C) –2.013 –2.894 –2.867
soil freezing days 114 112 127
Freeze-thaw cycles at 5 cm depth 4 1 5
estimated freeze-thaw cycles at 1 cm depth 18 20 12
average magnitude of soil freeze at 5 cm depth (˚C) –3.253 –4.469 –3.916
coldest soil temperature at 5 cm depth (˚C) –7.848 –10.325 –9.546
day of coldest soil temperature at 5 cm depth 12-Jan-07 29-Jan-08 26-Jan-09
first freeze 02-Nov-06 19-Nov-07 20-Nov-08
final thaw 10-Mar-07 10-Mar-08 31-Mar-09
length of winter freeze-thaw cycle 128 112 131
average lag time between freeze-thaw cycles 32 112 26
accumulated precipitation (mm) 52.415 36.943 69.118
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
Figure 7. B-spline curve fitted to average soil temperature (5 cm depth, R² = 0.99) in southern Alberta
between November 1 and April 1. Fitted soil temperature B-spline curve represents the boundary condi-
tion calibrated with a regional-specific context. All fitted soil temperature fit within the observed soil
temperature at 5 cm depth 95th percentiles.
5 cm than the previous winter (Tab le 6). More accumu-
lated winter precipitation is likely the main reason for
warmer soils, snow cover may have insulated the soil
from colder air temperatures [20]. The 2008/2009 winter
had the least number of soil freeze-thaw cycles at 1 cm
(12) and the most at 5 cm (5); this is likely caused by a
few intense warming events with a long enough duration
thaw deeper soil layers instead of many short-duration
low-intensity warming events. The 2008/2009 winter
had the latest day of initial soil freeze and final thaw at 5
cm but was still the longest winter as indicated by the
soil freezing days (Table 6).
Figure 7 shows the fitted B-spline curve fitted to the
average winter soil temperature time series for southern
Alberta. The B-spline fitted curve had a R² = 0.99 when
compared to the actual average winter soil temperature
at 5 cm depth and was used as the winter season regional-
specific portion of Eq.3. The B-spline cure fitted tends
to overestimate the peaks and under estimate troughs
observed in the average winter soil temperature time
series. The regional-specific soil surface boundary con-
dition is represented by Eq.3 and can be used to drive
freeze-thaw algorithms using regional-specific climate
data. Fitting sinusoidal curves to a temporally filtered
seasonal signal provides the regional-specific seasonal
boundary condition for seasonal soil temperature trends
in Eq.3. Fitting the B-spline curve the average winter
time series provides a regional-specific winter boundary
condition that can estimate winter freeze-thaw cycles.
The proposed model was developed using regional-spe-
cific climate data instead of generic boundary conditions
that require only one or two climate variable inputs.
13.247sin 0.0030.244
,Apr. 1Nov. 1
Apr. 1Nov. 1,
Nov. 1Apr. 1
Tt t
 
T = temperature [˚C]
t = time [s]
λ = decay rate
n = degree of B-spline
Pi = knot (control point)
bi = B-spline coefficient
Our findings illustrate the importance of regional
context in modeling freeze-thaw cycles. Our results in-
dicate that the interaction between air temperature and
wind drives variability in soil temperatures and if in-
cluded when specifying boundary conditions or cali-
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. Openly accessible at http://www.scirp.org/journal/AS/
brating to historical climate, may significantly improve
the prediction skills of models framed at many different
spatial and temporal scales, e.g. soil water and solute
vertical transport or ecosystem scale models. By consid-
ering a nested statistical modeling approach to capture
regional-specific freeze-thaw variability, we utilize a set
of best-fit measures (leave-one-out cross validation), and
geostatistical spine interpolation, to better identify the
influence of measurement and model-based uncertainty
on our model predictions of freeze-thaw variability over
our study region and both annual-scale and winter time-
The spatial analysis of the correlation of wind speed
and maximum air temperature revealed areas of southern
Alberta where the correlation between these two variables
varies most strongly. This kind of spatial analysis that
reveals the most influential climate monitoring stations
on model predictions and can be used to infer where
additional climate monitoring stations could be located
to better track climate signals and freeze-thaw dynamics
within the study region (southern Alberta). The interwinter
variability revealed that freeze-thaw cycling behaviour
varies differently in the winter months compared to the
full annual seasonal cycle. From this result, we infer that,
to more reliably predict and track freeze-thaw activity,
its spatial and temporal variability needs to be better
represented. Often historical data is limited, but none-
theless statistical modeling can be performed and used to
guide how new monitoring data can, in the future, be
collected to improve SVAT model prediction skill.
Specifically, our findings suggest that freeze-thaw sig-
nalling within southern Alberta is strongest during winter
months at several key monitoring stations. Where
climate-soil correlation is most statistically significant
are guiding regions where monitoring activity could be
sampled at higher frequency and at greater resolution
using networks of wireless soil moisture sensors.
Our findings also provide a statistical soil surface
boundary condition that captures both the seasonal and
winter-scale variability of soil temperature. This boundary
condition can be incorporated in SVAT and soil water/
solute vertical transport models. Our work illustrates how
models can be given regional context, by calibrating them
to a representative statistical model that captures main
climatic trends, regional variation, and climate-soil
interaction terms. In such a statistical representative
model, we show here that freeze-thaw cycling is one ex-
ample of a climate-soil interaction term. Using this
statistical approach, regional calibration of models that
require tracking of soil water movement can be improved,
especially in regions where freeze-thaw activity is strong.
Use of a single, simplified statistical model also potentially
offers a way to simplify climate-related inputs to complex
SVAT and agro-ecosystem models.
Future modeling work will explore how the spatial
pattern and inter-annual variation of soil freeze-thaw
cycling dynamics is regulated by climate teleconnections
and synoptic climate events. Related work will look at
the integration of fine scale climate and soil data with
multiple sensor resources to help advance the geostatis-
tical design, and implementation of agro-environmental
monitoring systems that combine in-situ and satellite/
remote-sensing derived estimates of near-surface soil
The authors would to acknowledge R. Scott Erickson (Agriculture
and Agri-Food Canada) and Prof. Steven Liang (University of Calgary)
for helpful editorial comments that improved an earlier version of this
manuscript. This work was made possible by funding from the Sus-
tainable Agriculture Environmental Systems (SAGES) Growing For-
ward Program, Agriculture and Agri-Food Canada (Federal Govern-
ment of Canada).
[1] Kreyling, J., Beierkuhnlein, C., Pritsch, K., Schloter, M.
and Jentsch, A. (2008) Recurrent soil freeze-thaw cycles
enhance grassland productivity. New Phytologist, 177,
938-945. doi:10.1111/j.1469-8137.2007.02309.x
[2] Chang, C. and Hao, X. (2001) Source of N2O emission
from a soil during freezing and thawing. Phyton, 41,
[3] Matzner, E. and Borken, W. (2008) Do freeze-thaw
events enhance C and N losses from soils of different
ecosystems? A review. European Journal of Soil Science,
59, 274-284. doi:10.1111/j.1365-2389.2007.00992.x
[4] Bullock, M.S., Larney, F.J., Izaurralde, R.C.S. and Feng,
Y. (2001) Overwinter changes in wind erodibility of clay
loam soils in Southern Alberta. Soil Science Society of
America, 65, 423-430. doi:10.2136/sssaj2001.652423x
[5] Walker, V.K., Palmer, G.R. and Voordouw, G. (2006)
Freeze-thaw tolerance and clues to the winter survival of
a soil community. Applied and Environmental Microbiology,
73, 1784-1792. doi:10.1128/AEM.72.3.1784-1792.2006
[6] Henry, H.A.L. (2008) Climate change and soil freezing
dynamics: Historical trends and projected changes. Cli-
matic Change, 87, 421-434.
[7] IPCC (2007) Climate change 2007: Synthesis report. In:
Pachauri, R.K and Reisinger, A., Eds., 4th Assessment
Report of the Intergovernmental Panel on Climate
Change IPCC, Geneva, 104.
[8] Eitzinger, J., Parton, W.J. and Hartman, M. (2000) Im-
provment and validation of a daily soil temperature sub-
model for freezing/thawing. Soil Science, 165, 525-534.
doi: 10.1007/s10584-007-9322-8
[9] Chen, B., Chen, J.M. and Ju, W. (2007) Remote sensing-
based ecosystem-atmosphere simulation scheme (EASS)—
Model formulation and test with multiple-year data.
A. J. Phillips et al. / Agricultural Sciences 2 (2011) 392-405
Copyright © 2011 SciRes. http://www.scirp.org/journal/AS/Openly accessible at
Ecological Modelling, 209, 277-300.
doi: 10.1016/j.ecolmodel.2007.06.032
[10] Flerchinger, G.N. and Saxton, K.E. (1989) Simultaneous
heat and water model of a freezing snow-residue-soil
system. I. Theory and development. Transactions of the
American Society of Agricultural Engineers, 32, 565-
[11] Flerchinger, G.N. (2000) The simultaneous heat and wa-
ter (SHAW) Model: Technical documentation. Northwest
Watershed Research Center, USDA Agricultural Research
Service, Boise.
[12] Šimůnek, J., Sejna, M., Saito, H., Sakai, M. and Genuch-
ten, M.T.V. (2009) The HYDRUS-1D software package
for simulating the one-dimensional movement of water,
heat, and multiple solutes in variably-saturated media.
PC Progress, Prague.
[13] Nkemdirim, L.C. (1996) Canada’s chinook belt. Inter-
national Journal of Climatology, 16, 441-462.
[14] Nkemdirim, L.C. (1991) Chinooks and winter evapora-
tion. Theoretical and Applied Climatology, 43, 129-136.
doi: 10.1007/BF00867470
[15] Hurrell, J.W. (1996) Influence of variations in extra-
tropical wintertime teleconnections on northern hemis-
phere temperature. Geophysical Research Letters, 23,
665-668. doi: 10.1029/96GL00459
[16] Bonsal, B.R., Prowse, T.D., Duguay, C.R. and Lacroix,
M. P. (2006) Impacts of large-scale teleconnections on
freshwater-ice break/freeze-up dates over Canada. Jour-
nal of Hydrology, 330, 340-353.
[17] Smith, N.V., Saatchi, S.S. and Randerson, J.T. (2004)
Trends in high northern latitude soil freeze thaw cycles
from 1988 to 2002. Journal of Geophysical Research,
109, 1-14. doi:10.1029/2003/2003JD004472
[18] Baker, D.G. and Ruschy, D.L. (1995). Calculated and
measured air and soil freeze thaw frequencies. Journal of
Applied Meteorology, 34, 2197-2205.
[19] Henry, H.A.L. (2007) Soil freeze-thaw cycle experiments:
Trends, methodological weaknesses and suggested im-
provements. Soil Biology and Biochemistry, 39, 977-986.
doi: 10.1016/j.soilbio.2006.11.017
[20] Groffman, P.M., Driscoll, C.T., Fahey, T.J., Hardy, J.P.,
Fitzhugh, R.D. and Tierney, G.L. (2001) Colder soils in a
warmer world: A snow manipulation study in a northern
hardwood forest ecosystem. Biogeochemistry, 56, 135-
150. doi:10.1023/A:1013039830323