Transfer function models to quantify the delay between air and ground temperatures in thawed active layers

Abstract. Air temperatures influence ground temperatures with a certain delay, which increases with depth. Borehole temperatures measured at 0.5 m depth in Alpine permafrost and air temperatures measured at or near the boreholes have been used to model this dependency. Statistical transfer function models have been fitted to the daily difference series of air and ground temperatures measured at seven different permafrost sites in the Swiss Alps. The relation between air and ground temperature is influenced by various factors such as ground surface cover, snow depth, water or ground ice content. To avoid complications induced by the insulating properties of the snow cover and by phase changes in the ground, only the mostly snow-free summer period when the ground at 0.5 m depth is thawed has been considered here. All summers from 2006 to 2009 have been analysed, with the main focus on summer 2006. The results reveal that in summer 2006 daily air temperature changes influence ground temperatures at 0.5 m depth with a delay ranging from one to six days, depending on the site. The fastest response times are found for a very coarse grained, blocky rock glacier site whereas slower response times are found for blocky scree slopes with smaller grain sizes.


Transfer function models to quantify the delay between air and ground temperatures in thawed active layers E. Zenklusen Mutter
1 , J. Blanchet 2 , and M. Phillips

Introduction
Ground temperature depends on air temperature, ground surface characteristics (e.g. the presence of vegetation or scree), snow cover and ground properties such as grain size and water/ice content (Goodrich, 1982;Harris and Pedersen, 1998;Gorbunov et al., 2004;Luetschg et al., 2004;Zhang, 2005;Gruber and Hoelzle, 2008).
Air temperature changes will sooner or later induce changes in the thermal state of the ground.The active layer is the uppermost part of permafrost, and is subject to seasonal thawing and refreezing (Burn, 1998) and therefore reacts first to changing weather and climate.Increasing air temperatures can lead to active layer deepening and ground instabilities, which can cause problems for infrastructure (Phillips and Margreth, 2008;Bommer et al., 2010) and in mountainous terrain lead to mass movements like rock falls or debris flows (Noetzli et al., 2003;Gruber et al., 2004b;Rabatel et al., 2008;Harris et al., 2009).It is therefore of interest to know how permafrost reacts to changing environmental conditions and in particular to changes in air temperature.
Various approaches have been applied to investigate the relation between air and ground temperature (e.g.Geiger, 1965;Harlan and Nixon, 1978;Beltrami, 1996;Frauenfeld et al., 2004;Smerdon et al., 2004Smerdon et al., , 2009)).For example Harlan and Nixon (1978) described the influence of the accumulated thawing temperatures, soil properties and moisture content on the active layer thickness.This approach was later applied in different studies, in both high latitude and high altitude permafrost regions (e.g. Brown et al., 2000;Hinkel and Nelson, 2003;Christiansen, 2004;Mazhitova et al., 2004;Åkerman and Johansson, 2008;Smith et al., 2009;Zenklusen Mutter and Phillips, 2011).
In the Alps the heterogeneity of the surface and subsurface, as well as the complex topography complicate the modelling of ground temperature (Gruber et al., 2004a;Noetzli et al., 2007).It is for example still a big challenge to model precipitation in mountainous terrain and hence estimate its influence on permafrost (Frei et al., 2003;Salzmann et al., 2007).Furthermore, there are still too few borehole sites in the Alps for a complete overview of the complex thermal conditions reigning in the ground.In order to analyze the future evolution of the thermal state and ground ice content in permafrost, empirical, analytical and numerical physically-based model approaches of different complexity are usually taken (Riseborough et al., 2008;Harris et al., 2009;Dall'Amico et al., 2011).In a recent study Engelhardt et al. (2010) performed sensitivity analyses on the influence of atmospheric forcing parameters (temperature and precipitation) on mountain permafrost evolution with a one-dimensional coupled heat and mass transfer model.They showed that at the permafrost site Schilthorn in the Introduction

Conclusions References
Tables Figures

Back Close
Full Swiss Alps, summer temperature changes during the snow-free period and the timing of the first snow cover in autumn have the largest impact on ground temperatures.Unlike these physically-based approaches, this study presents a statistical model to identify the temporal relation between day-to-day changes in the air temperature and those in the ground temperature of Alpine permafrost soils.Using a transfer function model, ground temperature changes measured in permafrost boreholes are linked to corresponding changes in air temperature.The goal is to quantify the temporal dependence of ground temperature changes on air temperature changes and find out whether this dependence is similar for all the different sites or if site-specific properties such as lithology, hydrology or ground-ice content alter the relation, and if so, to what extent.As the snow cover and thawing/freezing processes decouple ground temperatures from air temperatures (Luetschg et al., 2004;Zhang, 2005), only a three-month period in summer 2006 during which the active layer was established is here considered in detail.Ground temperatures measured in boreholes at 0.5 m depth in the active layer and air temperatures measured at or near the borehole sites are modelled statistically.Furthermore, all summers from 2006 to 2009 have been analysed in order to place the 2006 results in a larger context.

Data
Daily ground temperatures used for this work were measured in seven boreholes located in mountain permafrost at elevations between 2400 m a.s.l. and 3295 m a.s.l. in the Swiss Alps (Fig. 1 and Table 1).
Yellow Springs Instruments thermistors (YSI 44008) are used for the borehole temperature measurements, in combination with Campbell CR10X or CR1000 data loggers, with a calibrated precision of 0.01-0.02• C.

Conclusions References
Tables Figures

Back Close
Full are of different lengths and the depth of the boreholes at the study sites varies from about 20 m to more than 50 m.For the analyses carried out in this study mean daily ground temperatures in the active layer (at 0.5 m depth) and mean daily air temperatures either from a meteorological station at the borehole site or nearby (IMIS network, WSL Institute for Snow and Avalanche Research SLF) have been used (see Fig. 1 and Table 2).An overview of the mean daily air and ground temperature time series for summer 2006 (1 July to 30 September) measured at the different sites is given in Figs. 2 and 3.

Statistical methods and results
Statistical transfer function models have been applied to model the influence of air temperature on ground temperature.In the following, the procedure will be presented with example data from site A2 and summer 2006.The general framework follows the procedure of chapter 11.4 of Cryer and Chan (2010); see also Chapt. 11 in Box et al. (2008), Chapt.5.7 in Shumway and Stoffer (2011).
As there is no weather station close to borehole A2, air temperatures have been taken from the station Met.Ab which is located 9.2 km away at a comparable elevation (Tables 1 and 2).The courses of air and ground temperature show clear similarities (Fig. 3).However, there is a period of 23 days in the middle of the summer where discrepancies are visible.These are due to the presence of snow which attenuates the relation between both parameters.However, for simplicity these discrepancies have not been taken into consideration.
To determine the influence of air temperature on ground temperature, transfer functions were used to model the first differences of both temperature series (Figs. 4 and  5, top graphs).The idea behind this is to explain to which extent and with what delay, for example, a 10 • C change in air temperature from one day to another is represented in the ground temperature in summer.The processes resulting from differencing have zero mean.The application of the Phillips-Perron test (Phillips and Perron, 1988) and Introduction

Conclusions References
Tables Figures

Back Close
Full the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test (Kwiatkowski et al., 1992) reveals that the air and ground temperature difference series from site A2 can be assumed to be stationary.They are however autocorrelated (Fig. 4 and Fig. 5, bottom graphs) and mutually dependent, as are the initial ground and air temperature series themselves.Our goal is to quantify how changes in ground temperature are influenced by present and past changes in air temperature.According to Cryer and Chan (2010), the following general regression model relating the two time series can be considered: (1) In Eq. ( 1) parameter X t describes the input data (i.e. in our case the air temperature difference series), Y t the output data (i.e. in our case ground temperature difference series) and Z t corresponds to an error term, assumed to be independent of X t .Index t labels the time (i.e.days here) ranging from 1 to n. Parameter β k labels an infinite number of regression coefficients describing the dependency of the output on the input time series.In real applications however, a finite sum of coefficients is used and the model simplifies to Equation ( 2) is just a linear model for the ground temperature differences with the lagged air temperature differences as a covariate.It is known variously as the transfer function model, the finite distributed lag model, or the dynamic regression model.In our case, the input series X t (air temperature changes) influences the output Y t (ground temperature changes) but not the other way round, so m 1 and m 2 are expected to be non-negative.The values of the regression coefficients β k determine how quickly and efficiently ground temperature reacts to a change in air temperature.Roughly speaking, we expect to have large coefficients β k for small k (e.g.k = 0 or 1) in ground that reacts Introduction

Conclusions References
Tables Figures

Back Close
Full efficiently and quickly to air temperature, and moderately large β k up to a quite large lag (e.g.k = 5 or 6) in ground that has a diffuse response to air temperature.These coefficients β k need to be estimated from the data.There are two possibilities for this.
The first is to assume some structure on the relation of the coefficients β k .The most important structure is the Almon lag model (Almon, 1965), see chapter 15 of Judge et al. (1985) for a review.The second possibility is to impose no structure on the β k and let the data decide the shape of the lag structure.This is the approach we follow here.In Eq. ( 2), the coefficients β k are thus (m 2 − m 1 +1) unstructured regression coefficients to be estimated.Equation ( 2) can be written in vector form as where As the input X t and the output Y t are stationary autocorrelated (Figs. 4 and 5 for the example site A2), the same is expected for the error term Z t .Hence, we assumed in the regression model Eq. ( 2) that Z t follows some ARMA(p,q).Given arbitrary values for m 1 , m 2 , p and q, standard methods for dealing with correlated errors in linear models can be used to fit the transfer model Eq. ( 2).However, in practice m 1 , m 2 , p and q are unknown and reasonable values have to be chosen before being able to estimate the model.
We chose m 1 and m 2 according to the procedure described in Cryer and Chan (2010).The selection is made on a transformed version of model Eq. ( 2), obtained as follows.The stationary autocorrelated X t are assumed to follow some ARMA(p ,q ) model.The series X t being observed, we can use standard techniques of time series to choose reasonable values of p ,q and estimate the corresponding ARMA model for is then a white noise series.Applying the filter π(B) to both sides of Eq. ( 2) now gives where Because X is a white noise and X is independent of Z, the theoretical cross-correlation between the transformed processes X and Ỹ at lag k, Corr( Xt−k , Ỹt ), is then proportional to the regression coefficient β k .Thus in the cross-correlation plot (CCF plot) of Xt and Ỹt CCF values should be almost 0 for k < m 1 and k > m 2 , or at least nonsignificant.Therefore, the CCF plot between Xt and Ỹt can be used as a graphical tool for the selection of m 1 and m 2 .In Fig. 6 the CCF-Plot for site A2 is depicted, revealing m 1 = 0 and m 2 = 4 as possible candidates.However, when applied to the data of all seven sites, lags ranging from m 2 = 1 (site A3) to m 2 = 6 (site M2) are revealed (not shown).Considering different values of m 1 and m 2 for the different boreholes would impede intra-site comparisons of the estimated regression coefficients, and thus would make the intra-site comparison of the reaction speed of ground temperature changes to air temperature changes difficult.Another possibility therefore is to assume a fixed number of lags m 1 and m 2 in the regression model Eq. ( 2) for all sites and to distinguish the significant coefficients β k from the non-significant ones in the end.Following the cross-correlation plots of the prewhitened data from all sites (not shown), we decided to use m 1 = 0 and m 2 = 8.

Conclusions References
Tables Figures

Back Close
Full Once m 1 and m 2 have been selected, the choice of the orders p and q for the ARMA model for the residuals Z t in Eq. ( 2) remains.The difficulty arises because Z t is not observed.A first possibility would be to choose appropriate orders p and q based on the residuals obtained by fitting model Eq. ( 2) under the improper assumption of uncorrelated residuals, and then to use these orders to properly fit model Eq. ( 2).
However, we preferred at this point to have a more automatic criterion in order to be able to apply the procedure automatically to all our boreholes.The choice we made was to fit different models Eq. ( 2) with m 1 = 0 and m 2 = 8 and for different orders p and q, chosen to be relatively small to restrict the number of free parameters.We allowed p ranging from 0 to 8 and q ranging from 0 to 2. This leads to 27 different models Eq. ( 2) with 9 to 19 free parameters.All these models were fitted using the R function arima, with argument "xreg" set to the appropriate matrix X of Eq. ( 3).It uses the well-known Cochrane-Orcutt procedure (Cochrane and Orcutt, 1949) for fitting linear models with correlated errors.Finally the model with the lowest AIC value was selected as the best one (Akaike, 1973).
In summary, the applied procedure for fitting the transfer function model to the difference series X t and Y t can be described as follows: 1. Find an appropriate ARMA(p ,q ) for the X t .Prewhiten the input data X t with the corresponding filter, leading to Xt (Eq.4).
2. Transform the output data Y t by applying the same filter as in step 1 (Eq.5).Call the transformed output data Ỹt .
3. Create a cross-correlation plot between Xt and Ỹt .Deduce from it candidates for m 1 and m 2 .
4. Fit model Eq. ( 2) for this choice of m 1 and m 2 , with different orders p and q for the ARMA(p,q) model of the residuals Z t .Choose the one which best fits the data (lowest AIC).

TCD Introduction Conclusions References
Tables Figures

Back Close
Full For the example site A2 an ARMA(0,1) (i.e. an MA(1)) model has been selected for the residuals Z t (Table 3).In total this selected model Eq. ( 2) has 10 parameters.The residual analysis plots show that the residuals in the MA(1) model for Z t are uncorrelated and normally distributed (Fig. 7).The application of the Ljung-Box test (Ljung and Box, 1978) for site A2 and different lags ranging from 1 to 15 reveals P values from 0.5 (lag 9) to almost 1 (lag 3) and confirms that the residuals are "white noise".This is consistent with the theory of the transfer model.The comparison of the fitted with the measured difference data also shows a relatively good fit of the transfer model (Fig. 8, upper graph).However, positive differences seem to be slightly underestimated whereas negative differences are rather overestimated, leading to a little less fluctuation than observed.
Although model Eq. ( 2) applies to the difference data and not to the air and ground temperatures directly, one can now use these fitted values to reconstruct the time series of ground temperatures at day t + 1, given the ground temperature at day t.Let Y t be the observed ground temperature at day t.Hence, the predicted value Ŷ t+1 can now be obtained for t ranging from m 2 + 1 to n − 1 from the value Y t and the predicted first difference Ŷt of Eq. ( 2) by This so-called "one-day ahead prediction" can be used for example if there is a missing value at day t + 1 in the ground temperature series, but the air temperature is available.The comparison of these one-day ahead predicted and measured ground temperatures confirms the good fit of the transfer function model (Fig. 8, lower graph).
The estimated regression coefficients β k of model Eq. ( 2) for summer 2006 in A2 are shown in Fig. 9, dotted in black.They reach their maximum at lag k = 2, meaning that air temperature changes occur at 0.5 m depth with a delay of around two days at site A2, confirming the cross-correlation plot shown in Fig. 6.For site A2 all nine coefficients are significant (i.e. the absolute value of the estimated coefficient is larger than its standard error).The procedure described above was then applied to the data Introduction

Conclusions References
Tables Figures

Back Close
Full from the other six borehole sites.The estimated regression coefficients β k are plotted in Fig. 10.Comparing these coefficients for all sites reveals that depending on the site, the delay between air temperature changes and ground temperature changes at 0.5 m depth varies from about one to six days in summer 2006 (Figs. 9, 10 and Table 4).These differences can be explained by the different ground properties at the borehole sites.
The fastest response occurs at site A3, with a delay of one day (Table 4).The borehole also shows much higher values of the regression coefficients β k than the other sites (Figs. 9 and 10).This fast and efficient response can be explained by the very coarse, blocky surface cover at A3 (Fig. 11c) which allows an efficient relation between air and ground temperature changes.For site A2, which is located on a weatherexposed ridge (Fig. 11b), the delay is around two days (Fig. 9).On the contrary, scree slope sites A1, F1, M1 and M2 show more delayed and prolonged response times (Fig. 10, Table 4).The surface at these sites is not as coarse-grained as at site A3 (Fig. 11a,c,e,f) and the scree still contains air-filled voids, providing good thermal insulation.Although sites M1 and M2 are located only 50 m away from each other, differences in snow cover distribution (M1 lies between avalanche defence structures, which cause delayed snow melt, see also Fig. 11f) and hydrology (Rist and Phillips, 2005) lead to a slightly longer response time at M1 than at M2 (Table 4).Finally, in the artificially modified finer-grained terrain site A4, located at a chairlift midway station (Fig. 11d), similar responses as for the scree slopes are found, although less diffuse (Fig. 10c).
In addition to summer 2006, we applied the same procedure to data acquired in summers 2007, 2008 and 2009.The grey area in Figs. 9 and 10 depicts the enveloping curves of the estimated regression coefficients for these four summers.The overall picture is quite similar for all summers but some variability exists.A wider envelope is found for site A2 (Fig. 9), due to large variability in summer conditions here, and in particular due to the presence of snow in summer 2008.As already mentioned, snow cover can attenuate or even interrupt the relation between air and ground temperatures Introduction

Conclusions References
Tables Figures

Back Close
Full (see also Fig. 3) and a weak relation between air and ground temperature results in small regression coefficients β k .The presence of snow in summer 2008 in A2 was due to the presence of a huge snow cornice from the previous winter (K.Lauber, personal communication, 2008), causing the low model coefficients at the lower border of the grey area in Fig. 9. On the other hand, sites A3, A4 and F1 show very low variability over the four summers.At site A4 the borehole is located at the chairlift station where the regular artificial maintenance of the snow cover impedes the establishment of seasonally varying snow depths at the end of the winter.Seasonal differences at A4 are therefore restrained to summer snow events which have no long temporal implications.At A3 and F1 the narrowness of the envelope is probably simply due to similar conditions for the summers of the period 2006 to 2009.

Discussion and outlook
Transfer function models have been fitted to air and ground temperature data measured at 0.5 m depth at seven different permafrost sites in the Swiss Alps.Due to the influence of snow cover and ground ice on the relation between ground temperature and air temperature, the model has only been fitted for the mostly snow-free summer period when the ground at 0.5 m depth is thawed.Estimated model coefficients show that the ground temperature changes at 0.5 m depth at different permafrost sites depend on the air temperature changes which occurred about one to six days earlier, depending on the site characteristics.Very coarsegrained, blocky ground surfaces lead to faster response times.Smaller-scale blocks in scree slopes insulate the ground below and therefore induce longer response times.For some sites the relation between air and ground temperature -and therefore the coefficient estimates -were influenced by short periods with summer snow cover, which was not taken into account for this study.The study confirms the well-known influence of ground characteristics and snow cover on the thermal regime of the active layer and allows the quantification of the lag between air and ground temperature changes in different types of terrain.

Conclusions References
Tables Figures

Back Close
Full A subdivision of the year into snow-free and snow-covered, frozen and thawed time periods could be helpful to determine the role of various disturbing factors and is planned in future.Furthermore, the derivation of physical parameters such as heat capacity, thermal conductivity and hence thermal diffusivity from the model coefficients might be an interesting and important step to follow.Other questions of interest in this context are, for example, to estimate to what depth ground temperatures are influenced by daily or seasonal air temperature changes, or to determine the optimal vertical thermistor distribution and temporal resolution for borehole temperature measurements.

Conclusions References
Tables Figures

Fig. 1 .
Fig. 1.Locations of the borehole sites (open circles) and meteorological stations (stars) in Switzerland.

Fig. 2 .Fig. 3 .
Fig. 2. Air and ground temperature time series at 0.5 m depth for summer 2006 at the different sites (for site A2 see Fig. 3).

Fig. 4 .
Fig. 4. Air temperature day-to-day difference series measured at Met.Ab (top) and corresponding autocorrelation (lower left) and partial autocorrelation (lower right) plots.

Fig. 5 .
Fig. 5.As in Fig. 4 but for ground temperature day-to-day difference series measured at 0.5 m depth at A2.

Fig. 10 .
Fig. 10.As in Fig. 9 but for the remaining six sites (note the different scales on the y-axes).Significant (non-significant) coefficients are plotted as filled (empty) black dots.

Table 1 .
Investigated sites and their particular characteristics.

Table 2 .
Meteorological stations used in the study.

Table 3 .
ARMA(p,q)selected with AIC to fit the transfer model for the different sites.