Modeling the impact of wintertime rain events on the thermal regime of permafrost

Introduction Conclusions References


Introduction
Arctic permafrost areas represent a vast region which is expected to be strongly impacted by global warming in the coming decades to centuries.Simulations of future climate Correspondence to: S. Westermann (sebastian.westermann@geo.uio.no) using General Circulation Models (GCM's) suggest a warming of the near-surface air temperature of up to 10 K in the Arctic in the coming century, which is significantly more than the global average (Solomon et al., 2007).Based on the output of GCM's, a number of studies have modeled the thermal regime of permafrost soils, both for the entire permafrost domain (e.g.Delisle, 2007;Lawrence et al., 2008) and for selected regions (e.g.Etzelmüller et al., 2011).The results indicate a significant reduction of the total permafrost area and a pronounced deepening of the active layer in the remaining areas until 2100.The underlying process-based permafrost models used in all these studies assume heat conduction as the only process, through which energy is transferred in the ground and within the perennial snow pack.Advection of heat through vertical or horizontal water fluxes is generally neglected in thermal permafrost models (Riseborough et al., 2008).Significant advective heat transfer has been documented in spring, when meltwater infiltrates and subsequently refreezes in the frozen ground, causing a rapid increase of soil temperatures (e.g.Hinkel and Outcalt, 1995;Kane et al., 2001;Weismüller et al., 2011).However, this phenomenon only occurs during a very limited period of the year, so that it is generally not considered in permafrost models.
In this study, we demonstrate that advective heat transfer within the snow pack through infiltrating water can have a significant impact on the thermal regime of permafrost.During strong wintertime rain events, so-called "rain-on-snow" events, water can percolate to the bottom of the snow pack, where it refreezes, thus depositing considerable quantities of latent heat at the snow-soil interface and leading to the formation of basal ice layers (Woo and Heron, 1981).For a permafrost site in Svalbard, Putkonen and Roe (2003) show that few strong rain-on-snow events are sufficient to confine S. Westermann et al.: Modeling wintertime rain events the ground surface temperature (GST), i.e. the temperature at ground surface below the snow cover, to 0 • C for a large part of the winter season.Under such conditions, the "Bottom Temperature of Snow" (BTS) method (Haeberli, 1973), a commonly applied field technique to assess the occurrence and vitality of permafrost, can produce misleading results (Farbrot et al., 2007).
Winter rain events are common in permafrost areas in more maritime settings, like in Norway (e.g.visible in the data published in Isaksen et al., 2002) and Iceland (Etzelmüller et al., 2007;Farbrot et al., 2007).Furthermore, they are documented from Alaska, Arctic Canada (Grenfell and Putkonen, 2008) and the western-and easternmost parts of Siberia, e.g. on the Yamal and Chukotka peninsulas (Rennert et al., 2009;Bartsch et al., 2010).Using active microwave remote sensing, Bartsch (2010) estimates an average of one to three strong wintertime rain events per season in most areas on Svalbard, which is in an arctic maritime setting.Rain-on-snow events have received significant attention in biology and agricultural science, as the basal ice layers prevent ungulates from reaching plants and lichens under the snow, thus causing migration or starvation of the animals (Putkonen and Roe, 2003;Harding, 2003;Chan et al., 2005;Putkonen et al., 2009).Using field observations from W Svalbard and a simplified model of water infiltration in snow, we demonstrate that rain-on-snow events under certain conditions play a prominent role for the thermal regime of permafrost soils, so that they may deserve attention in model approaches targeting the permafrost evolution in a changing climate.

Study site
The study is performed on the Brøgger peninsula located at the west coast of Svalbard close to the village of Ny-Ålesund at 78 • 55 N, 11 • 50 E (Fig. 1).The area has a maritime climate influenced by a branch of the North Atlantic Current, so that the winters are relatively mild (average February air temperature −14 • C).The snow-free period typically lasts from July to September, but much longer (150 days) and much shorter (50 days) durations have been recorded (Winther et al., 2002).On average, about three quarters of the annual precipitation of 400 mm fall during the "winter" months from October to May, but the interannual variability of the winter precipitation is significant (Førland et al., 1997;eKlima, 2010).Strong wintertime rain events leading to the development of a bottom ice layer have been documented for the Brøgger peninsula for a number of years (Putkonen and Roe, 2003;Kohler and Aanes, 2004).
The Brøgger peninsula is located in the zone of continuous permafrost, with "lowland" permafrost being restricted to a 2 to 4 km wide strip between the Kongsfjorden and the locally glaciated mountain chain in the interior (Liestøl, 1977).The Bayelva climate and soil monitoring station about 2 km SW of the village of Ny-Ålesund (Fig. 1) has provided a record of climate and soil variables since 1998 (Roth and Boike, 2001;Boike et al., 2003Boike et al., , 2008;;Weismüller et al., 2011).It is located on top of Leirhaugen hill at an elevation of 25 m a.s.l., about half-way between the Kongsfjorden and the terminus of the glacier Brøggerbreen.The surface is covered by mud boils, a form of non-sorted circles, of diameters of about 0.5 to 1.0 m, so that exposed soil on the tops alternates with sparse vegetation in the depressions between the mudboils.The soil features a high mineral and low organic (volumetric fractions below 0.1) content, and the soil texture ranges from clay to silty loam (Boike et al., 2008).For this study, we use time series of outgoing long-wave radiation, from which the temperature of the ground or snow surface is calculated (Sect.3.3).Furthermore, snow depth, active layer temperatures and soil water contents through "Time-Domain-Reflectometry" (TDR) measurements (Roth et al., 1990) within the active layer (Table 1) are employed.As a record of GST, the soil temperature at 0.05 m depth below the surface in a depression between mudboils is used.The precipitation record from an unheated tipping bucket rain gauge, which only measures precipitation in the form of rain or slush, is employed to determine the occurrence of rain-onsnow events and to coarsely estimate the amount of rain.
In addition, we make use of measurements of the incoming longwave radiation at a station of the "Baseline Surface Radiation Network" (BSRN) located in the village of Ny-Ålesund (Ohmura et al., 1998).Furthermore, precipitation measurements from Ny-Ålesund (eKlima, 2010) are The Cryosphere, 5, 945-959, 2011 www.the-cryosphere.net/5/945/2011/employed, where regular maintenance of the instrumentation is guaranteed, unlike at the Bayelva station.Here, the total amount of precipitation is measured and the daily record is filed in the three classes "rain", "slush" and "snow" according to visual observations, which allows an independent determination of rain-on-snow events.However, the exact amount of liquid precipitation cannot be determined from either record, as the partitioning of slush in liquid and solid precipitation is unclear.

Model setup
The employed model is a thermal snow and soil model supplemented by a "cold-hydrology" scheme for percolation of rain water in snow.Unlike sophisticated snow schemes, such as SNOWPACK (Bartelt and Lehning, 2002;Lehning et al., 2002a,b), or fully coupled heat-and mass transfer models (Dall'Amico et al., 2011), such as COUP (Stähli et al., 1996;Gustafsson et al., 2004) or GEOtop (Zanotti et al., 2004;Rigon et al., 2006;Endrizzi et al., 2011), it does not include a comprehensive description of all natural processes in the snow and the soil.Instead, we only account for the processes that are most relevant for the formation of the thermal regime of the soil.As an example, water movement within the soil, which Weismüller et al. (2011) show to be of secondary importance for the thermal regime at the study site, is not included.

Soil thermal model
In the soil domain, we assume the temperature T to change over time t and depth z through heat conduction (Jury and Horton, 2004) as described by where k(z,T ) [Wm −1 K −1 ] denotes the thermal conductivity and c eff [Jm −3 K −1 ] the effective heat capacity (Jury and Horton, 2004), which accounts for the latent heat of freezing and melting of water/ice as L = 334 MJ m −3 denotes the specific volumetric latent heat of fusion of water.The first term is calculated from the volumetric fractions of the constituents as where θ α [-] and c α [Jm −3 K −1 ] represent the volumetric content and the specific volumetric heat capacity (following Hillel, 1982)  capacities as given in Table 2 are employed, which are in agreement with published values for the study site (Roth and Boike, 2001;Boike et al., 2008;Westermann et al., 2009).The resulting volumetric heat capacities for frozen and unfrozen soil are given in Table 2.The soil freezing characteristic θ w (T ) is determined by fitting the function to measurements of temperature and soil water content conducted at a depth of 0.89 m (Table 1), which is displayed in Fig. 2. The fit yields a value of δ = 0.17 • C, and a residual liquid water content in frozen soil of θ min w = 0.05 is assumed (see Roth and Boike, 2001 for a discussion of the accuracy of TDR measurements in this context).The ice content is then given by The thermal conductivity k is calculated following de Vries (1952) and Campbell et al. (1994) as where f α [-] denotes a weighting factor and k α the conductivities of the constituents water, ice, air and mineral (Hillel, 1982).The calculation of the f α is guided by the concept that one of the constituents occurs as interconnected "continuous phase" with thermal conductivity k c , while the other constituents are conceptualized as discontinuous phases, i.e. small domains intercepted by the continuous phase.Assuming spherical particles, the weighting factors can be calculated as (Campbell et al., 1994) www.the-cryosphere.net/5/945/2011/The Cryosphere, 5, 945-959, 2011 where k c denotes the conductivity of the continuous phase.
For unfrozen soil with fixed mineral, but varying water and air contents, Campbell et al. (1994) suggest a smooth transition from an air-dominated regime (k c = k a ) to a waterdominated regime (k c = k w ) by defining with The transition between air and water as continuous phase occurs at the "liquid recirculation cutoff" θ lrc [-], and the range, over which this transition takes place, is determined by the smoothing parameter s [-].Campbell et al. (1994) give values for both parameters for a range of soils that have been determined in laboratory experiments.We employ values of θ lrc = 0.1 and s = 3, which is in the range of the values found by Campbell et al. (1994) for silty and clayey soils (as found in the study area).We assume that the concept of Eq. ( 8) is applicable for air-ice systems (i.e.zero water content) and water-ice sys-tems (i.e.zero air content).In the absence of measurements or literature values, we employ the same parameters θ lrc and s for the air-ice system as for the air-water system.For the water-ice system, the choice of γ is uncritical, as the thermal conductivities of pure water and ice are not strongly different.We assume a linear interpolation between k w and k i according to the volumetric water and ice contents.Finally, the thermal conductivity of a soil with non-zero fractions of all constituents is obtained by interpolating between the three confining systems air-water, air-ice and water-ice, which span a three-dimensional space.
The resulting thermal conductivity for unfrozen soil is k thawed =1.45 Wm −1 K −1 , which is in good agreement with the published value of (1.3 ± 0.4) Wm −1 K −1 from measurements in the study area (Westermann et al., 2009).For freezing soil in the temperature range between −2 and −9 • C at the Bayelva station, Roth and Boike (2001) report a thermal diffusivity of k/c eff =8 × 10 −7 m 2 s −1 .From the slope of the freezing characteristic (Fig. 2), the volumetric fractions of all constituents and Eq. ( 4), c eff is estimated to be between 2.5 and 3.0 MJ m −3 K −1 for this temperature range.This results in thermal conductivities between 2.0 and 2.4 Wm −1 K −1 , so that the value of k frozen =2.5 Wm −1 K −1 obtained from the conductivity model for fully frozen soil is reasonable.
For depths of more than 10 m, we assume constant thermal properties (Table 2), as observations during a drilling campaign close to the Bayelva station in 2008 suggest that bedrock is found at this depth.We emphasize that the exact choice of the values does not strongly affect the simulations of GST on the considered timescales.

Snow thermal and hydrological model
We derive both the build-up and ablation of the snow cover from snow depth measurements at the Bayelva station, which in conjunction with measured precipitation rates determine the mass balance of the snow pack.If no liquid precipitation occurs, we assume heat transfer within the snow pack to be fully governed by conductive heat transfer.In case of a rain event, water freezes within the snow pack, which increases the volumetric ice content θ i and releases latent heat, thus increasing the snow temperatures towards the freezing point of water, 0 • C. Following Jury and Horton (2004), the governing equation of heat transfer within the snow pack becomes For the density ρ and the volumetric heat capacity of snow, we use published values from measurements of Westermann et al. (2009) at the study site (Table 3).To avoid inconsistencies due to mechanical effects during freezing of water, the density of ice is set equal to the density of water, so that a density of dry snow of 350 kg m −3 corresponds to a volumetric ice content of θ i =0.35.For simplicity, we choose thermal conductivities which linearly increase with time from the date of the snowfall for each grid cell.The confining conductivities k fresh and k old (Table 3) for fresh snow and for 280 days old snow (approximately corresponding to the length of the snow-covered period) are again derived from Westermann et al. (2009): while the snow density and thus the volumetric heat capacity are found to be relatively homogeneous in time and space, an increase of measured thermal diffusivities k snow /c snow at the bottom of the snow pack (which has not been strongly affected by rain events in the studied year) from 4.5 × 10 −7 m 2 s −1 to 7.0 × 10 −7 m 2 s −1 is recorded from December to March.This corresponds to thermal conductivities of 0.3 and 0.55 Wm −1 K −1 at the bottom of the snow pack.Considering that the timing of the first snowfall in the study period has been similar to the year studied by Westermann et al. (2009), the linear conductivity scheme with confining values k fresh and k old roughly reproduces the measured values for December and March.
Infiltration of rain water is described by an extended bucket scheme, which distributes the liquid precipitation rate P / t [m 3 m −2 s −1 ] over the snow pack.In addition to temperature, two further variables are used to characterize a grid cell in each time step, the volumetric water content θ w and the sum of volumetric ice and water contents, θ tot = θ w + θ i .For each snow cell, we calculate two maximum infiltration rates: ], leading to the infiltration of an amount of water, which refreezes and provides the energy to raise the temperature to or sustain it at 0 • C, plus infiltration  ], which allows infiltration until the sum of ice and water content is unity, i.e θ tot = 1.
The infiltration rates define the governing differential equations for θ tot and θ w and the term ∂θ i /∂t in Eq. ( 10), which couples the hydrological to the thermal scheme.Details on the employed equations are given in Appendix A. The infiltration scheme is schematically depicted in Fig. 3: the precipitation rate is distributed from top to bottom using I snow for each grid cell, until the entire precipitation rate is accounted for (Fig. 3a).When the sum over the entire snow pack of I snow is not sufficient to absorb the precipitation rate, water starts pooling at the bottom of the snow pack, which is treated by distributing I bottom from bottom to top (Fig. 3b) in the model.After the refreezing of this water, a bottom ice layer with thermal properties different from the snow forms (Table 3), for which literature values for ice are assigned (Hillel, 1982).Other than that, the thermal properties of the snow remain unchanged during and after an infiltration event in the model.Since both the exact amount of rain and the field capacity θ fc w can only be estimated, we use a fixed value of 0.01 for the field capacity and manually adjust the solid-to-liquid ratio of the precipitation to fit the modeled GST to the measurements.For 2005/2006, we assume that 60 % of the daily precipitation recorded at the Bayelva station (which measures rain and slush) occur in liquid form.In 2006/2007, where only precipitation data from Ny-Ålesund are available, only a single strong rain event occurs in late March.For this event, we consider the recorded precipitation as rain.

Initial and boundary conditions, numerical scheme
The method of lines (Schiesser, 1991) is employed to numerically solve the heat transfer equation (Eq. 1) for a soil domain with 100 m depth.The spatial derivatives are discretized using finite differences, leaving the time as continuous variable.The resulting system of ordinary differential equations (ODE) with time as variable is solved numerically in MATLAB with the ODE-solver "ode113" (Shampine and Gordon, 1975;Shampine and Reichelt, 1997), which continuously adjusts the integration time step to minimize computational costs.The grid spacing is increased with depth, with 0.02 m between the surface (defined as 0 m) and 1.6 m, 0.2 m between 1.6 m and 5.0 m, 0.5 m between 5.0 m and 20.0 m, 1.0 m between 20.0 m and 30.0 m, 5.0 m between 30.0 m and 50.0 m and 10.0 m between 50.0 m and 100.0 m.
If snow is present, additional grid cells with a grid spacing of 0.02 m are added on top of the soil.The position of the uppermost cell is determined from measurements of the snow depth, which are interpolated to the center positions of the snow grid cells.Within the snow pack, the numerical scheme is applied to the variables (T θ tot θ w ), which are defined by the system of coupled differential equations given by Eq. ( 10), Eqs.(A4) to (A8) and Eq.(A9).
As upper boundary condition, we use surface temperatures T surf based on measurements of outgoing and incoming long-wave radiation.For calculation, we use use Stefan-Boltzmann and Kirchhoff's Law, where σ S [ Wm −2 K −4 ] denotes the Stefan-Boltzmann constant and ε [-] the surface emissivity.As this study focuses on periods, where the ground is snow-covered, we use a constant surface emissivity of 0.985, which is in the range of typical values published for snow surfaces (e.g.Dozier and Warren, 1982;Hori et al., 2006).At the lower boundary, we prescribe a heat flux of 60 mW m −2 at a depth of 100 m, which is in the range of estimates for the geothermal heat flux on Svalbard (Liestøl, 1977;Isaksen et al., 2000;Van De Wal et al., 2002).
To a depth of 1.52 m, the initial condition is inferred from soil temperature measurements at the Bayelva station (Table 1), between which the temperatures are linearly interpolated.Below 1.52 m, no temperature measurements are available, so that the temperature distribution can only be estimated.For the season 2005/2006, we use the record of the lowermost temperature sensor at 1.52 m (which has been continuously in frozen ground) from July 2002 to June 2005 to generate the steady-state temperature distribution for this forcing, which is employed as initial condition below 1.52 m.This results in an initial temperature of −2.9 • C at 1.5 m, −3.8 • C at 3 m, −3.1 • C at 10 m and −3.1 • C at 20 m depth.Below, a stable gradient of 0.024 Km −1 (determined by the heat flux through the lower boundary and the conductivity of the bedrock, Table 2) forms, thus placing the base of the permafrost at 150 m depth, which is in agreement with estimates of permafrost thickness in coastal areas of Svalbard (Humlum, 2005) If the snow depth increases and a new grid cell must be added from one time step to the next, the temperature of the newly added grid cell, T + , must be specified (for decreasing snow depths, the uppermost grid cell is simply deleted).We use as defining differential equation, so that T + follows the upper boundary condition with a characteristic lag time τ + .We set corresponding to the timescale of heat diffusion through snow over a distance z (set to the grid spacing of 0.02 m).This choice ensures that the integration timesteps selected by the ode-solver are adequate both for the snow temperatures and for T + .The volumetric ice and water content θ tot of a snow cell is initialized to 0.35 according to snow density measurements (Table 3, Westermann et al., 2009), while the water content θ w is set to zero.

Field data from the winter seasons 2005/2006 and 2006/2007
We investigate the winter seasons 2005/2006 and 2006/2007, which are characterized by near-surface air temperatures significantly above the long-term average (Isaksen et al., 2007;eKlima, 2010).In both years, the perennial snow cover forms in the course of September and lasts well into June of the following year.In 2005/2006, the snow depth is between 0.5 and 1.0 m for most of the winter period (Fig. 4).At the rain gauge at the Bayelva station, a total precipitation of 250 mm is recorded during the period when snow is present.At the same time, a total precipitation of 373 mm is measured in the village of Ny-Ålesund, of which 40 % are classified as snow, 47 % as slush and 13 % as rain.As strong differences in the precipitation between the two sites are not likely, these numbers suggest that the 250 mm measured at the Bayelva station contain a considerable fraction of solid precipitation.Therefore, the total amount of liquid precipitation can only be estimated within wide margins of error from the measurements.
The measurements at the Bayelva station yield an average surface temperature of −8.4 • C for the winter season of 2005/2006.The measured GST under the snow, however, is much warmer and does not fall significantly below −1 • C for most of the winter season.In December and January, rain events occur, after which the GST is confined to close to 0 • C for several weeks.About 175 mm of precipitation are recorded in Ny-Ålesund in this period, of which 40 mm are classified as rain and most of the remaining part as slush.At five days, the total measured precipitation exceeds 10 mm.In February, a stronger soil cooling is initiated: the minimum GST value of about −3 • C is reached in the end of March after a prolonged period with snow surface temperatures as low as −25 • C, before the temperature increases again in the course of spring (Fig. 4).The average GST during the period when the ground is covered by snow is −0.6 • C.
In the winter period of 2006/2007, the measured snow depth is slightly higher than in the previous season, with values between 0.7 and 1.1 m for most of the winter.Due to problems with the rain gauge at the Bayelva station in this winter, only the precipitation record from Ny-Ålesund is available.There, a total of 283 mm of precipitation is recorded, of which about 54 % are classified as snow and 42 % as slush, while rain occurs only at the beginning and at the end of the winter season.Only one strong precipitation event with a total of about 45 mm occurs in late March, which leads to a sharp increase of the measured GST towards 0 • C (Fig. 5), to where the temperature is confined for about ten days.The average measured GST during the period, when the ground is covered by snow, is −1.1 • C, while the average snow surface temperature is −9.8 • C.
From the yearly record of GST, the freezing and thawing indices, Fr and Th, are calculated as the accumulated sum of degree days for negative and positive values of GST.Using the thermal conductivity for frozen and thawed soil (Table 2), we can evaluate a simple criterion for permafrost occurrence (e.g.Carlson, 1952), www.the-cryosphere.net/5/945/2011/The Cryosphere, 5, 945-959, 2011 For the period from July 2005 to June 2006, this condition is clearly violated, with k thawed Th being about twice as large as k frozen Fr.For the following year, both terms are approximately equal.The weighting might be even more shifted towards the right side of Eq. ( 14), since the soil is not fully frozen during part of the winter, so that the average thermal conductivity is smaller than k frozen .In any case, the threshold, below which permafrost is sustainable on long timescales, is clearly exceeded for the environmental conditions encountered in the winter season 2005/2006.

Model results
The strong differences between the snow surface and the ground surface temperature encountered in both winter seasons form as a combination of two factors.First, the thermal insulation of the snow, which has a much lower thermal conductivity than the underlying soil, delays the refreezing of the active layer, so that GST is sustained at higher values.Second, wintertime rain events cause a considerable warming of the snow and soil and can sustain GST at values close to 0 • C for prolonged periods.To separate both effects, we use the permafrost model described in Sect.3, which can simulate the thermal effect of rainwater infiltration on the soil.In addition to the model run, a control run is performed, where the infiltration routine is deactivated and the temperature distribution is only determined by heat conduction.
Figure 4 depicts the modeled temperature distribution within the snow pack and the uppermost meter of the soil.In addition, the modeled ground surface temperatures (taken as the temperature at 0.05 m depth, see Sect. 2) from both the model and the control run are displayed.At the beginning of the season, the model slightly underestimates GST.A possible explanation is that the snow depth at the measurement location in the depression between mudboils is higher than inferred from measurements with the sonic ranging sensor (Table 1), which averages over a larger footprint area.The impact of infiltrating rain is reflected in the position of the −0.01 • C isotherm (chosen since the modeled temperatures never reach exactly 0 • C for computational reasons): after a rain event, the snow temperatures close to the surface rapidly cool, while the temperatures in deeper snow layers remain around 0 • C, until all water is frozen.During that time, the heat transfer to the snow surface is impeded by the overlying snow layers, so that the freezing of the infiltrated water in the snow occurs slowly.If the infiltration is sufficiently strong that water pools at the bottom of the snow pack and large quantities of water are stored in the bottom snow layers (due to the maximum infiltration rate I bottom , see Fig. 3b), the effect is particularly pronounced and confines the modeled GST to close to 0 • C for a prolonged period.This situation occurs in January and February 2006, where the model run reproduces the measured GST with reasonable accuracy, while the control run shows considerably cold-biased temperatures (Fig. 4), which persist for the rest of the winter season.In contrast, the rain events between September and December 2005 have only minor impact on the modeled GST, so that the results of the model and the control run do not deviate strongly.As a result of the strong rain events, a 0.1 m thick bottom ice layer forms.Ice layers of similar thickness are documented for the area around Ny-Ålesund (Kohler and Aanes, 2004).This relatively conductive ice layer is still overlain by more than 0.5 m of snow (Fig. 4), so that the heat conduction through the snow pack is not significantly enhanced after the rain-on-snow event.
The model run yields an average GST of −0.8 • C for the period, when the ground is covered by snow, which is close to the measured value of −0.6 • C. The modeled average GST in the control run is −1.5 • C, which still appears to be a good approximation of the measured GST.However, the freezing index Fr (which is proportional to the average GST in • C) obtained from the control run is almost twice as large compared to the the model run.Therefore, the permafrost is not sustainable according to Eq. ( 14) in the model run, while the control run suggests that the left and right side of Eq. ( 14) are approximately equal.
In the winter season 2006/2007, the GST obtained from model and control run do not deviate until March, despite of a few small rain events (Fig. 5).The impact of the strong rain event in late March on GST is well reproduced in the model run, while the control run again yields cold-biased values of GST.Compared to the previous winter, the deviation between model and control run is smaller, with an average GST of −1.6 • C in the control and −1.3 • C in the model run, compared to a measured average GST of −1.1 • C. The rain event results in the formation of a 0.04 m thick bottom ice layer.
In both seasons, the modeled GST is considerably colder than the measured GST at the end of the winter period, from end of April in 2006 and from mid of May in 2007 (Figs. 4, 5).The step-like increase of measured GST is most likely caused by infiltration of water from melting snow, which is not accounted for in the current scheme.During this period, the agreement could most likely be improved, if nearsurface melt rates are calculated from a surface energy balance scheme (Boike et al., 2003).The infiltration of the melt water could then be treated similar to rain water.

Long-term impact of repeated winter rain events on the ground thermal regime
It is evident from the measured freezing and thawing indices, that the permafrost at the study site is not sustainable for the environmental conditions encountered in the winter season 2005/2006 (Sect.4.1).To investigate the speed of a potential degradation and assess the contribution of the wintertime rain events, we force the model with the data set from July 2005 to June 2006 for ten consecutive years.As initial condition, we again prescribe the soil temperature distribution of 1 July 2005, as detailed in Sect.3.3.We emphasize, that The Cryosphere, 5, 945-959, 2011 the permafrost temperatures are comparatively low at this time, with temperatures of less than of −3 • C between 2 and 20 m depth.As the shape of the soil freezing characteristic has a pronounced impact on subsurface temperatures (Romanovsky and Osterkamp, 2000), we perform the simulation for two different parameters δ determining the amount of unfrozen water at subzero temperatures.On the one hand, we use δ = 0.17 • C, as estimated for the silty to clayey soil at the Bayelva station (Fig. 2) and employed for the simulations in the previous section.On the other hand, the simulation is performed for a steeper soil freezing characteristic (δ = 0.04 • C), which corresponds to a soil of a higher sand and gravel content.
The results of the model and the control runs are displayed in Fig. 6, which shows the fraction of unfrozen soil water relative to the potentially freezable amount, i.e. (θ w − θ min w )/(θ max w − θ min w ), for the uppermost four meters of the soil.While a zone with a considerable fraction of unfrozen soil water develops in all simulations, the degradation occurs considerably faster for the model run, which includes the effect of rain water infiltration, compared to the control run.For the flat soil freezing characteristic, the final temperature and soil water distribution of the control run after ten years is reached after roughly four years in the model run.Simultaneously, the maximum thaw depth increases by more than half a meter within three to four years in the model run, while the increase is slower in the control run.For the steep soil freezing characteristic, the speed of degradation is even higher, with a zone with constantly unfrozen soil water forming after only three to four years in the model run, which after five years already extends over more than half a meter of soil.Again, the process is slower in the control run, although a zone with constantly unfrozen soil water develops at the end of the ten-year period.Sandy and gravelly areas, which most likely feature a steeper freezing characteristic, exist in the vicinity of the Bayelva station.A study by Westermann et al. (2010) demonstrates thaw depths of up to two meters in 2008 for a gravel plain located approximately 200 m from the Bayelva station, which might at least partly be a consequence of the warm winters investigated in this study.
While the initialization below 1.52 m depth (Sect.3.3) and the assumptions on the thermal conductivity of the bedrock (Sect.3.1) introduce some degree of uncertainty for the ten ten-year simulation, all conclusions remain fully valid for slightly perturbed input values.

Modeling wintertime rain events -challenges and improvements
The presented scheme makes use of a simplified bucket model of water infiltration in the snow, instead of solving the governing physical equations (e.g.Colbeck, 1972Colbeck, , 1979;;Illangasekare et al., 1990).Nevertheless, it incorporates two important physical constraints: -Before water infiltration beyond a snow layer is possible, the temperature of the layer must be raised to 0 • C through freezing of the corresponding amount of water.
-A certain water content must be established and maintained to increase the hydraulic permeability of the snow to a level that facilitates flow of water.
In the employed model, the two constraints are satisfied by choosing sufficiently small values for the parameters τ 1 and τ 2 (Eqs.A1, A2), which control the timescales of the temperature increase to 0 • C and the saturation to the field capacity within a grid cell.We employ values on the order of seconds, so that (at least for realistic precipitation rates) significant infiltration in a grid cell can only occur when both conditions are satisfied for the grid cell above.Therefore, the model produces a progressing wetting front in the snow, which is fundamental for correctly describing water infiltration in dry snow (Colbeck, 1976).In addition to the precipitation rate, the progression of the wetting front depends on the temperature and one free parameter, the field capacity θ fc w of the snow, which can be adjusted to control the infiltration dynamics in dry snow.The physical interpretation of θ fc w is the water content, at which the hydraulic permeability switches from a very small value, which essentially prevents a flow, to a very large value, where flow is unobstructed and water movement occurs instantly.Accordingly, a precipitation signal instantly contributes to the water pool at the bottom of the snow pack, when the two conditions are fulfilled for the entire snow pack (i.e. the snow is "wet").In this point, the employed model differs from more sophisticated approaches describing the delay and dispersion of a precipitation pulse in wet snow by a hydraulic equation (e.g.Conway and Benedict, 1994).While the infiltration dynamics is not reproduced, the employed model still allows a good approximation of the amount of water reaching the bottom of the snow pack, which is the relevant quantity for the thermal regime of the soil.Furthermore, the timescale of infiltration in wet snow is on the order of minutes to hours (e.g.Singh et al., 1997), while water can persist at the bottom of the snow pack for several weeks (Sect.4.1), so that the infiltration dynamics is of minor relevance for the thermal regime of the permafrost.
In the model, we use a constant of value of θ fc w = 0.01, independent of snow properties, which is at the lower end of observed field capacities (denoted water holding capacities in some studies).While (Conway and Benedict, 1994) and Singh et al. (1997) report field capacities of more than 0.05, Kattelmann (1986) outlines that field capacities of less than 0.02 are common.Furthermore, water infiltration in snow features a high degree of spatial variability, as it preferentially occurs along localized flow finger or through fissures in ice layers in the snow pack (Marsh and Woo, 1984a,b;Conway and Raymond, 1993), so that a spatially averaged value like θ fc w (as employed in a 1-D-model) could be lower than suggested by measurements within a preferential flow path.
As exact rates of liquid precipitation are not available in this study, we cannot independently determine θ fc w by fitting the model output to observed GST values, so that it remains unclear whether the use of a constant value is sufficient for the purpose of modeling the thermal regime of permafrost.However, more sophisticated modeling schemes that take into account a possible dependence of infiltration on snow properties (e.g.Colbeck, 1979) generally require specification of additional parameters, for which a comprehensive set of field measurements is not available.Therefore, it seems questionable that the performance of a permafrost model would benefit from such improvements, particularly considering the strong uncertainties associated with both measurements (Førland and Hanssen-Bauer, 2000) and model predictions of precipitation (e.g.Serreze and Hurst, 2000) under arctic conditions.

Impact assessment under different environmental conditions
The model runs presented in Sect. 3 illustrate that rain events only have a strong impact on the soil temperatures, if water percolates to the bottom of the snow pack.This occurs if the amount of rain exceeds a threshold related to both temperature and water holding capacity of the snow.The energy required to cause an increase of the temperature of a snow pack of 1 m depth with the thermal properties given in Table 3 by 1 K corresponds to the latent heat released by the freezing of about 2.3 mm of water.Furthermore, for a field capacity of 0.01 as assumed in the model runs, 10 mm of water can be stored in the snow pack.For an average snow temperature of −5 • C, the threshold of liquid precipitation can thus be estimated to about 20 mm, which is only exceeded by few rain events in the two investigated winter seasons.After rain water has percolated to the bottom of the snow pack, the freezing rate is controlled by the heat flux through the snow pack.For a constant snow surface temperature of −10 • C, a snow depth of 1 m and a thermal conductivity of 0.5 Wm −1 K −1 , a heat flux of 5 Wm −2 is sustained, so that 10 mm of water at the bottom of the snow pack would freeze within eight days.For a snow depth of 0.5 m, the freezing would occur within four days.In reality, the freezing would occur faster than this simple estimate suggests since the downward heat flux in the ground, which leads to a warming of the underlying soil, dissipates some energy in addition to the upward heat flux through the snow pack.
Although large snow depths increase the threshold that must be exceeded for water to reach the bottom of the snow pack, they also delay the freezing of the infiltrated water and thus increase the impact on soil temperatures.Therefore, rain events have the strongest impact on the soil temperatures if the snow depth is high and large amounts of rain fall within short periods, as it has been the case for the two events in January 2006 and March 2007.Furthermore, if repeated rain events occur within short periods of time, the threshold can be considerably lowered, if the snow temperatures are close to 0 • C and water is still present in the snow pack.

Accelerated permafrost degradation through wintertime rain events?
The multi-annual simulations presented in Sect.4.3 suggest that rain-on-snow events can significantly accelerate the warming of soil temperatures in permafrost areas.As the long-term climate change is superimposed by considerable climate variability on shorter timescales, such non-linear processes can modify the response to the long-term warming signal: the occurrence of wintertime rain events may amplify the warming of permafrost temperatures in case of positive temperature anomalies, while only the much slower process of heat conduction leads to subsequent cooling in case of negative anomalies, thus resulting in a stronger net warming.Furthermore the model results highlight the possibility, that an initially stable permafrost system can quickly enter a state of degradation if environmental conditions as encountered in the winter 2005/2006 persist for several consecutive years.Thus, climate and weather extremes leading to an increased frequency of wintertime rain events may increase the vulnerability of permafrost to climate change.
Wintertime rain events play a larger role in permafrost regions with a maritime climate compared to continental areas, where they are only reported in spring (Groisman et al., 2003).Thus, the described processes are especially relevant for Svalbard and for mountain permafrost areas in Northern Europe.In Iceland, winter rain events are common at almost all elevation levels, producing iso-thermal snow conditions (Etzelmüller et al., 2007) and even contributing to the partial or total melting of the snow pack (Farbrot et al., 2007).In Norway, winter rain events seems common over the entire mountain area and have recently been recorded at elevations of more than 1600 m (Farbrot et al., 2011).Regional climate projection for Scandinavia indicate an increase of winter precipitation, but also an increase of precipitation intensity and air temperature for the mountain areas (e.g.Benestad, 2005;Beldring et al., 2008).This suggests that winter rain events will increase in frequency, so that they may accelerate the projected degradation of permafrost in these regions (Hipp et al., 2011).
Furthermore, the pronounced changes in the arctic and sub-arctic climate system projected for the next 100 yr (Solomon et al., 2007) may lead to the occurrence of wintertime rain events in permafrost area, where they have not been recorded previously or where such events have been very rare.Rennert et al. (2009) have scanned the output of GCMs for synoptic conditions correlated to the occurrence of winter rain events.Their analysis suggests an increase of both the frequency and the areal extent over the next 50 yr, with northwestern Canada, Alaska, and the pacific regions of Siberia being among the most affected areas.
Wintertime rain events are not accounted for in the models used to predict the future extent and thermal regime of permafrost.This study suggests that incorporating the effects of such events could improve predictions at least for regions where they occur regularly.While the presented model could in principle be driven with precipitation data sets obtained from climate models, it is questionable whether such coarsescale models can sufficiently capture the threshold nature of wintertime rain events.Furthermore, climate models deliver average precipitation rates for large grid cells, which cannot reproduce the subgrid variability typical for strong precipitation events.Due to the non-linearity of the rainwater infiltration in the snow, statistical (Wilby et al., 1998) or stochastic (Bates et al., 1998) downscaling algorithms may be required in order to correctly reproduce the net effect of wintertime rain events on long timescales.

Conclusions
We present measurements of the ground surface temperature conducted during two winter seasons on W Svalbard, where a number of wintertime rain events have occured.A simplified model of rainwater infiltration coupled to a transient permafrost model is employed to separate the impact of the rain events on soil temperatures from the thermal effects of the insulating snow cover.From this study the following conclusions can be drawn: -Small amounts of rain freeze close to the snow surface and thus have negligible impact on the soil temperatures.
-Strong rain events can modify the ground surface temperature over prolonged periods.During such periods, a permafrost model based solely on conductive heat transfer cannot reproduce the measured ground surface temperatures and rain water infiltration must be taken into account.
-In one of the investigated winter seasons, an average ground surface temperature of −0.6 -In a model simulation, where the environmental conditions encountered in this winter are repeatedly applied for several consecutive winters, clear evidence of permafrost degradation is visible after three to five years, as a zone develops, where only a minor part of the soil water refreezes.In a control run, where rain water infiltration is not accounted for, twice the time is required to reach a similar state.
-Strong wintertime rain events have a distinct impact on soil temperatures.As a process of highly non-linear nature, they have the potential to generate an accelerated warming in case of climate extremes, during which they occur with increased frequency.In areas with a maritime climate, they deserve attention in modeling schemes targeting the future thermal conditions of permafrost soils.

Governing equations of snow thermal and hydrological regime
As a standard ODE-solver is employed for the numerical time integration, the model system is fully defined if the time derivatives of the three governing variables, (T θ tot θ w ), are given for each grid cell (see Sect. 3.3).In addition to Eq. (10), which defines the time derivative of T , the ODEs for θ tot and θ w are required, plus an additional equation for ∂θ i /∂t, through which the hydrological is coupled to thermal scheme (see Eq. 10).The maximum infiltration rates (see Sect. 3.2) are defined by where T f is the freezing temperature of water (set to 0 • C).In the model runs, we use τ 1 = τ 2 = 10 s, so that the timescale of infiltration in a grid cell is short compared to the timescale of conductive heat transfer.If θ tot = θ w + θ i reaches one for a grid cell, both I snow and I bottom are set to zero, thus preventing further infiltration.Note, that this condition is inherent in the definition of I bottom , but not in I snow .The derivative of θ tot with respect to time is calculated with a recursive scheme.The snow grid cells are indexed from top to bottom as j = 1,...,N , with N being the grid cell above the snowsoil interface (Fig. 3).Using finite differences (see Sect. 3.3), Eqs. (A1) and (A2) are discretized in space with a grid spacing z j for the j -th cell, so that we obtain the infiltration rates I snow This yields the derivative of the ice content θ i with respect to time as which is employed in Eq. ( 10), where this term describes the effect of refreezing water on the temperatures in the snow pack.Thus, the defining system of ODEs, from which the model can be fully reproduced, are Eq.(10) for T , Eqs. (A4) to (A7) for θ tot , Eq. (A8) for θ w and Eq.(A9) for θ i .
We briefly illustrate the application of the scheme for a few relevant cases: in case the infiltration rate for a grid cell is equal to the maximum infiltration rate, i.e. ∂θ tot /∂t = I snow , the derivative of the ice content becomes Similarly, the liquid water content exponentially saturates towards the field capacity θ fc w with the governing differential equation In the following, the temperature is sustained at T f , as the amount of freezing water balances the conductive heat flux to the neighboring grid cells, which is determined by the gradient of the heat flux, ∂/∂z(k snow ∂T /∂z).When the infiltration process ceases, ∂θ tot /∂t becomes zero.For a grid cell with θ w > 0 and T = T f , the water content decreases according to as heat is flowing out into neighboring colder grid cells.With ∂θ i /∂t = −∂θ w /∂t, the grid cell is sustained at T = T f until all water is refrozen (θ w = 0).Then, ∂θ i /∂t becomes zero and normal heat conduction resumes (compare Eq. 10).

Fig. 1 .
Fig. 1.Map of Svalbard with the location of Ny-Ålesund.The inset shows the area around Ny-Ålesund with the location of the Bayelva station (bold lines: roads; thin lines: 10 m contour lines).

Fig. 3 .
Fig.3.Schematic description of the infiltration scheme, with snow grid cells numbered from j = 1,...,N .In case of light precipitation, the precipitation rate is distributed from top to bottom (a), while water starts pooling at the bottom for strong precipitation (b).Note that the maximum infiltration rates I snow and I bottom of each grid cell vary depending on temperature, water content and total water and ice content.See text.
. For the season 2006/2007, the initial condition below 1.52 m is obtained by forcing the 2002-2005 steady-state conditions with measured 1.52 m-temperatures from July 2005 to June 2006.

Fig. 4 .
Fig. 4. Top: modeled temperature distribution (in • C) within the snow pack and the first meter of soil for the winter season 2005/2006.The −0.01 • C isotherm is depicted in blue.Bottom: measured and modeled (including control run) ground surface temperature for the same time interval.In addition, the snow depth in m (right axis) and daily rainfall in mm d −1 (as employed in the model runs) are displayed.

Fig. 5 .
Fig. 5. Top: modeled temperature distribution (in • C) within the snow pack and the first meter of soil for the winter season 2006/2007.The −0.01 • C isotherm is depicted in blue.Bottom: measured and modeled (including control run) ground surface temperature for the same time interval.In addition, the snow depth in m (right axis) and daily rainfall in mm d −1 (as employed in the model runs) are displayed.

Fig. 6 .
Fig. 6. Results of scenario runs, where the forcing of the period July 2005 to June 2006 is applied for ten consecutive years.Fraction of unfrozen soil water relative to the potentially freezable amount (0 %: θ w = θ min w ; 100 %: θ w = θ max w , see Eq. 4).The model is initialized with measured (0 to 1.5 m depth) and estimated (>1.5 m depth) soil temperatures of 1 July 2005.Left: model run including infiltration of rain water in the snow.Right: control run.Top: flat soil freezing characteristic leading to a high content of unfrozen soil water at subzero temperatures (δ = 0.17 • C, as assumed for the study site).Bottom: steep soil freezing characteristic leading to a low content of unfrozen soil water at subzero temperatures (δ = 0.04 • C).
case depicted in Fig.3a, ∂θ tot /∂t of the first and the n-th grid cell (n = 2,...,N ) are calculated recursively from top to bottom as Eq.A3) is not fulfilled, i.e. for the case depicted in Fig.3b, ∂θ tot /∂t of the N-th and n-th (n = N −1,...,1) grid cells are calculated recursively from bottom to top as spatial derivatives instead of finite differences are given in the following.The derivative of the volumetric water content θ w with respect to time is calculated as temperature exponentially approaches T f with time constant τ 1 asT (t) = T f + T (t 0 )e −t/τ 1 .(A12)

Table 1 .
Data sets employed in this study.If not mentioned otherwise, the installations are located at the Bayelva station.

Table 3 .
Snow and ice properties employed in the model.