Non-stationary temporal characterization of the temperature proﬁle of a soil exposed to frost in south-eastern Canada

. The objective of this work was to compare time and frequency ﬂuctuations of air and soil temperatures (2-, 5-, 10-, 20- and 50-cm below the soil surface) using the continuous wavelet transform, with a particular emphasis on the daily cycle. The analysis of wavelet power spectra and cross power spectra provided detailed non-stationary accounts with respect to frequencies (or periods) and to time of the structure of the data and also of the relationships that exist between time series. For this particular application to the temperature proﬁle of a soil exposed to frost, both the air temperature and the 2-cm depth soil temperature time series exhibited a dominant power peak at 1-d periodicity, prominent from spring to autumn. This feature was gradually damped as it propagated deeper into the soil and was weak for the 20-cm depth. Inﬂuence of the incoming solar radiation was also revealed in the wavelet power spectra analysis by a weaker intensity of the 1-d peak. The principal divergence between air and soil temperatures, besides damping, occurred in winter from the latent heat release associated to the freezing of the soil water and the insulation effect of snowpack that cease the dependence of the soil temperature to the air temperature. Attenuation and phase-shifting of the 1-d periodicity


Introduction
The soil temperature has important biological, agricultural and climatic consequences (Grundstein, et al., 2005). Indeed, soil temperatures modulate the response of many soil biological and biochemical processes (Lloyd and Taylor, 1994;Quyang and Zheng, 2000;Luo et al. 2003;Beltrami and Kellman, 2003). For instance, soil organic matter decomposition and the resulting release of CO 2 to the atmosphere are dependent on soil temperatures (Kätterer et al. 1998;Kirschbaum, 2006). Coupled with air temperature, soil temperature in the upper few hundred meters of the earth's crust may also inform on past climatic conditions (Putman and Chapman, 1996;Schmidt et al., 2000;Bartlett et al., 2006;Pollack et al., 2005;Smerdon et al., 2006). Furthermore, soil temperatures below 0 • C influence snowmelt runoff and soil hydrology, especially by reducing soil permeability (Niu and Yang, 2006). Consequently, the modelling of such ecosystem processes often involves soil temperatures as a driving variable.
Soil temperature data are typically not available at a regional scale but has to be estimated, usually from empirical or semi-empirical relationships with standard meteorological data including soil properties and vegetation (e.g., Zheng et al. 1993;Kang et al. 2000). Although mainly influenced by anomalies in the atmospheric circulation (Hu and Feng, 2003), the soil temperature integrates the effects of air temperature, vegetation, snow cover, net radiation at the surface (Beltrami, 2001) and heat transfer through soils with different thermal properties in the vertical direction (Hu and Feng, 2003). Soil temperatures at various depths thus result from surface energy processes and from regional environmental forcing. This local response of soil temperature to atmospheric forcing varies both in time and depth.
Available observation networks of soil temperature have been examined by Zhang and al. (1996), Ouyang and Zheng (2000), Schmidt et al. (2000), Beltrami and Kell- man (2003), Hu and Fend (2003), Pollack et al. (2005), Niu and Yang (2006) and Smerdon et al. (2006). The soil thermal response in the field typically shows a diurnal pattern mostly driven by the temperature of the lower atmosphere. At the surface, Bartlett et al. (2006) reported air temperature explaining as much as 94% of the total variance of the soil temperature. Soil temperature in a vertical profile is governed by a heat diffusion process (Clauser, 1982). Putnam and Chapman (1996), Beltrami (2001), and Bartlett et al. (2006) showed that, for unfrozen ground, soil temperature tracks surface air temperature with amplitude attenuation and phase lag as a function of depth and soil properties. The incoming solar radiation also influences the soil temperature fluctuations. If surface soil temperature closely matches surface air temperature during the night, the soil temperature may be significantly warmer during the daylight hours, particularly in summer and autumn (Putnam and Chapman, 1996). Soil temperature time series are thus typically positively biased when compared to air temperature time series. The sensitivity of the offset caused by the incident solar radiation is largely attributable to the specific surface conditions at the site and may not be a feature of the soil-air tracking at all locations (Bartlett et al., 2006).
During winter, the soil is much less influenced by air temperature variations because freezing of the upper soil introduces a non-conductive mechanism (latent heat) that is effective until the upper layers of the soil are completely frozen (Beltrami, 2001); the latent heat released by freezing water warms up the surrounding soil (Luo et al., 2003). The presence of a snowpack may further reduce the conduction process at the surface (Beltrami and Mareschal, 1991). Both the annual number of snow cover days and the thermal resistance explain a considerable amount of the variance in the thermal offset between air temperatures and soil temperatures (Grundstein et al., 2005).
The objective of this work was to compare time and frequency fluctuations of air and soil temperatures at different depths using the continuous wavelet transform, with a par-ticular emphasis on the daily cycle. With wavelet analysis, a comprehensive description of the structure of soil temperature profile, its evolution over time, and how it relates to air temperature observations can be obtained without the stationary assumptions of standard Fourier frequency analysis. This may be particularly advantageous for the description of temperature time series of soils exposed to frost because of the non stationary behaviour of incoming solar radiation influences in high latitude, of the non-conductive mechanism that takes place only during freezing, and of the insulation properties of a fluctuating snowpack. Within the context of global warming, interest in soil temperature has been increasing because of its effect on the carbon release through soil organic matter decomposition -this release of carbon constitutes a dangerous positive feedback that leads to further warming. Several soil organic matter simulation models include a parameter that affects the decomposition rates as a function of soil temperature (e.g., Parton et al., 1987;Andrén and Katterer, 1997). These parameters are commonly estimated from standard weather station data such as mean daily or monthly air temperatures. Consequently, knowledge of the air and soil temperature relationships is valuable to soil carbon models.
In this study, soil temperature amplitude attenuation and time-lag with depth were evaluated as a function of the incoming solar radiation for a site exposed to frost in southeastern Canada. The energy exchange within the top 50 cm of the soil and the bottom 10 m of the atmosphere was examined using time series data for 1000-, 600-, and 150-cm high air temperature, soil temperature at 2-, 5-, 10-, 20-and 50-cm below the soil surface, incoming solar radiation and snowpack height.

Site and dataset description
The Mesonet-Montreal network is made up of of 43 meteorological stations in the Saint-Lawrence River valley. The Howick station that was used in this study is located at 45.0202 • N and 73.7566 • S, 43.9 m above sea level in southeastern Canada (Fig. 1). The station was situated in a cropped field. The soil at this site is classified as a Humic Gleysol, St-Urbain and Ste-Rosalie soil series (Canada Soil Survey Committee, 1978). The soil texture of the Ap horizon (≈0 to 30 cm) ranged from a clay loam to a silt clay loam with a soil organic matter content of 6.8%. The underlying Bg soil horizon (≈30 to 70 cm) was characterized by a texture ranging from clay to heavy clay.
Climate 1971-2000 normal values for the near-by Sainte-Martine automatic weather station are 6.3 • C, 813 mm of rainfall, and 171 cm of snowfall. The air-soil station was programmed to report observations every 5 min. All observations were averaged over the previous minute, except for rainfall that was 5-min total. For the purpose of the present study, data were averaged or summed in order to produce 3-h observations to reduce noise and to obtain Fourier periods of 6 h, or more.
Rainfall was measured with a tipping bucket gauge (TB3: Hydrological Services PTY), snowpack height with a sonic ranger sensor (SR50: Campbell Scientifics), incoming solar radiation with a pyranometer (CNR1: Kipp and Zonen), air temperature with a thermometer (HMP45C: Vaisala) 150 cm above ground level inside a Stevenson standard shelter and with thermometers (HMP45C: Vaisala) 600 cm and 1000 cm above ground level in naturally aspirated 10-plate radiation shields, and soil temperature with thermocouples 2, 5, 10, 20 and 50 cm below ground level (STP01: Hukseflux).
For a discrete sequence of observation, x n , the continuous wavelet transform W n is defined as the convolution product of x n with a scaled and translated wavelet ψ(η) that depends on a non-dimensional time parameter η as follows: Where n is the localized time index, s is the wavelet scale, δt is the sampling period, N is the number of points in the time series and the asterisk indicates the complex conjugate. Note that it is considerably faster to perform these convolution calculations in Fourier space (Torrence and Compo 1998).
Since a complex wavelet leads to complex continuous wavelet transforms, the wavelet power spectrum defined as W X n (s) 2 provides a convenient description of the fluctuation of the variance at different times and frequencies.
A number of functions may be used for wavelet transformation. An admissible wavelet function must have zero mean and be represented in both time and frequency domains (Farge, 1992). The selection of a particular wavelet function depends mainly on the experimenter's perspective and objectives. However, an appropriate wavelet should show patterns similar to the time series (Nakken, 1999). The Morlet wavelet shows peaks and troughs in a wavelike fashion suitable to air and soil temperature time series and provides a more acute definition in the spectral space than other common candidates such as the Mexican hat. The Morlet wavelet  ( Fig. 2) is a complex non-orthogonal wavelet consisting of a plane wave modulated by a Gaussian function as follows: Where ω 0 is the non-dimensional frequency equalled to 6 to satisfy the admissibility condition stated above. Spectral analysis often reveals that the variance of geophysical time series is organized along some preferential frequency bands that extend over a specific range of scales (e.g. Rajagopalan and Lall, 1998). The scale-averaged wavelet power allows to further examine fluctuations in variance over such bands of wavelet periods. The scale-averaged wavelet power is the weighted sum of the wavelet power spectrum over selected scales s 1 to s 2 (Torrence and Compo, 1998): Where δj is a factor fixing the scale resolution (chosen as 0.1) and C δ is a reconstruction factor specific to each wavelet form (C δ =0.776 for the Morlet wavelet). Similar to a cross-correlation analysis, cross wavelet phase angle evaluates the phase difference between the components of two time series. If W XY n (s) =W Y n (s) W Y * n (s) is the cross wavelet power spectrum from time series X and Y , then their phase angle is φ XY η (Labat et al., 2000) as follows: as computed from the real (ℜ) and imaginary (ℑ) parts of the cross wavelet power spectrum. In practice, W XY n (s) must be previously smoothed (Torrence and Webster, 1999;Grinsted et al., 2004;Maraun and Kurths, 2004;Maraun et al., 2007). The extent to which the two time series are out of step or the phase angle is expressed in radians between -π and π. If the phase angle is π/2 radians, the second series is lagging the first one by a quarter of a cycle. It may be convenient to convert the phase angle into a time-lag by means of T φ XY n /2π, where T is the wavelet period.

Time series
The time series covered 862 d from 26 September 2003 to 1 February 2006 (Figs. 3 and 4). At the Howick air-soil station, annual air temperature fluctuations showed highs around 30 • C during summer (June to September), and lows around -30 • C during winter (December to March). The time series appeared to be noisy mainly due to diurnal patterns (Fig. 3). Differences between the three air temperature series for elevation ranging from 150 cm to 1000 cm were small altough the 150-cm air temperature was slightly colder in average (Table 1). The incoming solar radiation was observed from December 2004, with maximum values occurring at summer equinox. The incoming solar radiation annual peak preceded that of air temperature by about a month. Snowfall occurred between December and April. The sonic ranger was installed in December 2004 and rainfall (blue) was abundant the rest of the year.  As for air temperature, the soil temperature time series combined the influence of an annual component and of a smaller diurnal component (Fig. 4). However, soil temperature differed from air temperature in many aspects. In particular, the amplitude of the annual and diurnal components was attenuated with depth. For instance, the 2-cm series ranged from -10.7 to 31.2 • C, and the 50-cm series ranged from -2.7 to 21.0 • C ( Table 1). The diurnal component faded away gradually with depth and almost vanished at 50 cm. As shown in Fig. 4, the annual cycle was lagging with depth. However, it was not possible, from the same figure, to ascertain the lagging of the diurnal cycle. It is also noteworthy that mean soil temperatures are higher than the mean air temperatures (Table 1). In fact, on a monthly basis (results not shown) the 2-cm depth soil temperature is always higher than the air temperatures. This offset is caused, as stated before, by the combination of phenomena typical for soil exposed to frost in high latitude: incoming solar radiation heating during longer daylight periods, non-conductive mechanism associated to freezing, and snowpack insulation. Amplitude and phase relationships between air and soil temperatures may be visualised by a simple graph of soil temperature as a function of air temperature, what Beltrami (2001) named "phase-space plots". Departures from the interception line, produced by the sinusoidal oscillation of yearly cycles, yields qualitative information about the correlation of the temperature time series, in both amplitude and phase. In this instance, "phase-space plots" of 3-h temperature data shows the effect of time-lag across frequencies (Fig. 5). Near the surface (Fig. 5a), time-lags were small and, as long as soil temperature remained above freezing, the air-soil temperature couples followed the interception line. Below freezing, the air-soil couples changed slope under the combined influence of latent heat dissipation and snow cover insulation. Soil temperature remained constant at 0 • C for some time. At the 20-cm depth (Fig. 5b), more curly features developed as a result of increased daily time-lag, and the slope of the air-soil couples below freezing was influenced even more by latent heat dissipation and snow cover insulation. At the 50-cm depth (Fig. 5c) the interception ellipse was dominated by the annual component of the soil time series. The winter air temperature penetration into the soil was very limited.

Wavelet transform analysis
The continuous wavelet transform, proposed for the first time to our knowledge for the non-stationary temporal characterization of these series, assumes changes in spectral power properties over time and can simultaneously account for both time and frequency domains, in contrast with the Fourier transform (Daubechies, 1992).
In the Morlet power spectra of air temperatures (150 cm) and soil temperatures (2 and 20 cm) time series (Fig. 6), the abscissa is time and the ordinate is the equivalent Fourier period related to wavelet scale. Interpretations may be extended to other air temperature elevations and soil temperature depths. Colour contours are the variance in excess of 0.1, 1.6, 6.4, and 25.6 • C 2 . This selection graphically depicts the 1-d variance. The colour associated with each variance  bands is displayed in the legend along with their occurrence. The red curve gives the cone of influence beyond which the edge effects may be of concern (Torrence and Compo, 1998). Any feature beyond 64 d was not considered in this study in line with the duration of the series, i.e. 862 d. The annual cycle was thus not formally identified. The Morlet power spectrum of air temperature exhibits a dominant power peak at 1-d periodicity (Fig. 6a). This nonstationary feature is particularly prominent from spring to autumn. For periods less than one day, much lower variance occasionally displays features at 0.5-d periodicity. Variance is also low for periods up to 2 to 4 d, indicating no recurrent meteorological phenomenon. The sharp white edges at high frequencies (variance less than 0.1 • C 2 ) are all caused by missing observations of air temperature that were linearly filled before analysis. They shall be neglected hereafter.
Soil temperature follows air temperature very closely (Beltrami, 2001). Wavelet analysis diagnoses similitude or dissimilitude as a function of periodicity. In the Morlet power spectrum of soil temperatures at the 2-cm depth (Fig. 6b), a dominant power peak occurs again at 1-d periodicity. However, time span is more limited than for air temperature series. The large white patches (Fig. 6b) extending essentially over the 4-month episode from early December to the end of March, and over periodicity up to 4 to 8 d, refer to variance smaller than 0.1 • C 2 (There were no missing data in the soil temperature series). Hence, the soil does not follow the winter air temperature time series, most likely because of the introduction of a non-conductive mechanism (latent heat) and the insulation effect of snowpack. Furthermore, the soil damping effect at the 2-cm depth temperatures (Fig. 6b) showed a less pronounced peak at 1-d periodicity, while the low variance band for periods longer than one day extended up to 4 to 8 days. Finally, the influence of incoming solar radiation produced weaker intensity at the 1-d peak during the fall, because of shorter daylight periods and of a lower sun at solar noon.
All features described by the Morlet power spectrum of the soil temperature at 2-cm depth were more pronounced at the 20-cm depth (Fig. 6c). The soil damping effect on temperature considerably weakened the 1-d peak and further extended the low variance band to periods up to 8 to 16 d. For instance, the variance for periodicity less than about 0.75 d was always less than 0.1 • C 2 and the winter episode increased in both time and periodicity.
The scale-averaged power spectrum allows further examination of the fluctuations in variance of the 1-d periodicity, more specifically over the 0.8-to 1.25-d wavelet period band. The point here is to better assess the non stationary behaviour of the time series of the dominating frequency band. The air temperature variance was indeed non-stationary and lower in winter than in the rest of the year (Fig. 7a). These modulations are only partly transferred to the 2-cm soil temperatures, where there is even less stationarity than for air temperatures (Fig. 7b). The 2-cm soil temperature variance is basically absent in the winter, active in spring and summer, and gradually reduced during autumn. Except in winter, this phenomenon may to some extent be linked to the incoming solar radiation (Fig. 3d) that is nearly stable in spring and summer, and drops gradually during autumn. The attenuation of fluctuations at the 20-cm depth was as expected stronger than at the 2-cm depth (Fig. 7c).
In addition to the amplitude of the attenuation, there is a time-lag to the air temperature variations with soil depth. As mentioned before, the time-lag associated to the heat diffusion process is easily computed from the wavelet phase angle of two time series. The time-lag for 1-d oscillations (more specifically for the band from 0.8 to 1.25 d) of the data at hand was computed from the phase angle that originates from the cross wavelet spectra of air temperature at 150-cm and soil temperature at 2-, 5-, 10-, and 20-cm depths (Fig. 8) note that the variance associated to other frequency bands is too low to lead to reliable time-lag estimations. As expected, the time-lags increased with depth. The median value was 1.4 h at the 2-cm depth, 2.4 h at the 5-cm depth, 3.8 h at the 10-cm depth, and 6.9 h at the 20-cm depth (10% and 90% quantiles are given in Table 2). No time-lag results are provided for the winter season, again because the variance is very low at that time (see Fig. 7). These results cannot be generalized to other sites since the heat diffusion process depends on the local soil properties. However, the results confirm the potential of the methodology for non stationary soil temperature time series, which may be useful to anyone exploring soil temperature models..

Conclusions
For this particular application to the temperature profile of a soil exposed to frost, when the characterization was calculated for a periodicity between 6 h and 64 d, both the air temperature and the near-surface soil temperature (2-cm depth) time series exhibited a dominant power peak at 1-d periodicity. A peak that was particularly prominent from spring to autumn. This feature gradually damped as it propagated deeper into the soil and was weak for the 20-cm depth. Influence of the incoming solar radiation was also revealed in the Morlet power spectra analysis by a weaker intensity of the 1-d peak during the fall because of shorter daylight periods and of a lower sun at solar noon. The principal divergence between air and soil temperatures, besides damping, occurred in winter from the latent heat release associated to the freezing of the soil water and the insulation effect of snowpack that cease the dependence of the soil temperature to the air temperature. Attenuation and phase-shifting of the 1-d periodicity could be quantified through scale-averaged power spectra and timelag estimations over the 0.8-to 1.25-d wavelet period band. Air temperature variance was only partly transferred to the 2-cm soil temperature time series and much less so to the 20cm soil depth. Unfrozen soil time-lags were about 1.4 h at 2 cm and 6.9 h at 20 cm.
Wavelet analyses are descriptive techniques that nicely complement the arsenal of tools employed for the diagnostic analysis of time series. The analysis of wavelet power spectra and cross power spectra provided detailed non-stationary accounts with respect to frequencies (or periods) and to time of the structure of the data and also of the relationships that exist between time series. These results are encouraging and suggest that this type of analysis, if expanded on a large number of sites, could enable us to provide more general conclusions with respect to soil temperature profiles at a regional scale. At a time when adaptation to and attenuation of climate change are prominent environmental and social issues, the carbon source/sink relationships between soil and atmosphere cannot be overlooked anymore. Climate is the major determinant of soil biological activity, including the decomposition rates of soil organic carbon. Better assessment of the air and soil temperature relationships, at a regional scale, is one of the components that can improve soil carbon balance models.