The Cryosphere Modeling the thermal dynamics of the active layer at two contrasting permafrost sites on Svalbard and on the Tibetan Plateau

Employing a one-dimensional, coupled thermal and hydraulic numerical model, we quantitatively analyze high-resolution, multi-year data from the active layers at two contrasting permafrost sites. The model implements heat conduction with the de Vries parameterization, heat convection with water and vapor flow, freeze-thaw transition parameterized with a heuristic soil-freezing characteristic, and liquid water flow with the Mualem-van Genuchten parameterization. The model is driven by measured temperatures and water contents at the upper and lower boundary with all required material properties deduced from the measured data. The aims are (i) to ascertain the applicability of such a rather simple model, (ii) to quantify the dominating processes, and (iii) to discuss possible causes of remaining deviations. Soil temperatures and water contents as well as characteristic quantities like thaw depth and duration of the isothermal plateau could be reproduced. Heat conduction is found to be the dominant process by far at both sites, with nonconductive transport contributing a maximum of some 3 % to the mean heat flux at the Spitsbergen site, most of the time very much less, and practically negligible at the Tibetan site. Hypotheses discussed to explain the remaining deviations between measured and simulated state variables include, besides some technical issues, infiltration of snow Correspondence to: K. Roth (kurt.roth@iup.uni-heidelberg.de) melt, dry cracking with associated vapor condensation, and mechanical soil expansion in detail.


Introduction
Permafrost is present on about one quarter of the Northern Hemisphere's land surface.With respect to hydrology, biology, and mass movement, the active layer with its annual freeze-thaw cycles is the most important part of permafrost landscapes.Aspects that depend crucially on the thermal and hydraulic dynamics of the active layer range from erosion and local hydrology all the way to the global climate, which crucially depends on the regional hydrology in northern regions, on their emission of greenhouse gases, and on the balance between different energy fluxes.Depending on the perspective, a quantitative representation is required at different scales that range from local, with extents of a few meters, for the detailed understanding of the pertinent processes, to global for the exploration of scenarios for future developments of our physical environment.
Here, we focus on the local scale, where the underlying processes are understood in sufficient detail to allow trustworthy predictions, e.g. for the development of the thermal and hydraulic dynamics at a particular site with changing forcing.The significance of such studies reaches well beyond their primary scale in that they provide structurally Published by Copernicus Publications on behalf of the European Geosciences Union.
correct templates for the representation of permafrost at larger scales.
Already a superficial consideration of the dynamics of the active layer reveals a disturbing spectrum of processes that may have to be taken into account, several of them highly nonlinear, several of them strongly coupled with each other.This remains true even if notoriously difficult aspects like vegetation and soil-atmosphere coupling or mechanical deformations and the associated pattern formation are neglected.This was summarized in the recent review by Riseborough et al. (2008).Since direct measurements for most processes and in the natural environment are not feasible with current technology, the approach is to model them and to compare the corresponding representations with available data, typically time-series of some state variables at a few locations.To this end, a wide range of analytical, semianalytical and numerical models have been used already.
The most widely applied analytical approach to soil freezing and thawing is based on Stefan's model (Carslaw and Jaeger, 1959).It presumes that (i) heat is transported from the soil surface to the thawing front through conduction only, (ii) the temperature gradient is constant in space, hence temperature relaxes instantaneously from any disturbance, and (iii) the transported heat is consumed completely at the thawing front for melting more water.This model yields a sharp thawing front whose penetration is in rough qualitative agreement with observations.Hence the two processes, heat conduction and solid-liquid phase change, were identified as the dominating ones.Nakano and Brown (1971) suggested the incorporation of a finite freezing zone instead of a sharp freezing front into this model to account for the unfrozen water that may be present at sub-zero temperatures.They included latent heat as an effective heat capacity within the freezing zone, and presented a numerical scheme for the solution.Harlan (1973) coupled Richards' equation with the heat equation based on Fourier's law and solved it using a finite-difference scheme in order to reproduce the movement of water and heat from warm to cold regions.Guymon and Luthin (1974) presented a very similar model and solved it using a finite-element scheme.Engelmark and Svensson (1993) proposed a solution scheme based on finite-volumes that directly yields the ice and water content and successfully checked the results against laboratory measurements.Hinzman et al. (1998) used a onedimensional subsurface model for simulating heat conduction and coupled it with a surface heat balance model to predict the thaw depth of the Kuparuk River watershed on the North Slope of Alaska from meteorologic data at seven stations.Harris et al. (2001) developed an analytical approximation for the heat transport problem in porous media, explicitly accounting for non-equilibrium processes occurring during rapid phase change.
While many studies assume purely conductive heat transport, Kane et al. (2001) examined different non-conductive heat transport processes that can be expected in freezing soils, and finds that infiltration of liquid water can contribute significantly to soil thawing in spring.Daanen et al. (2007) worked with a three-dimensional finite-difference model based on the heat equation coupled with Richards' equation and used this to examine liquid water movement during freezing.Romanovsky and Osterkamp (1997) compared a model similar to the one by Guymon and Luthin (1974) with two other commonly used models, one based on Stefan's problem and the other based on a finite-difference scheme with front tracking.
This paper aims at quantifying the extent to which the thermal and hydraulic dynamics in the active layer of permafrost soils can be represented by a rather simple model that combines (i) effective heat conduction, (ii) liquid-solid phase transition described by a static and non-hysteretic soilfreezing characteristic, and (iii) Richards' equation for soil water flow.Conversely, this also yields an estimate for the contribution of non-conductive processes to the movement of heat.Our primary focus is on the overall representation of all the phases during the annual cycle.This focus is relevant for a large-scale perspective, where simplified representations of the active layer are required that do not sacrifice qualitatively important phenomena.The secondary focus is on deviations between model and data, and on the discussion of possible causes.We emphasize right at the outset that obviously our findings are strictly applicable only to the two sites studied here.However, they are representative for a large range of permafrost landscapes.Important classes of sites that are not represented by them are hill slopes and shallow active layers in wet regions with strong radiative forcing.

Study sites
We consider two contrasting environments, one in the Arctic, the Bayelva site on Spitsbergen, and the other one on the Western Tibetan Plateau, the Tianshuihai site in the desert region of Aksai Qin, Xinjiang Uigur Autonomous Region, China (Fig. 1).These two sites differ with respect to radiative forcing, precipitation, and soil material.

The Bayelva site
This site and the corresponding data were described earlier, e.g.Roth and Boike (2001), Boike et al. (2003) and Boike et al. (2008), and we just repeat the pertinent features here.
The Bayelva site is located about 3 km west of Ny Ålesund at 78 • 55 N, 11 • 50 E.The active layer is underlain by continuous permafrost to several hundred meters depth.Mean annual precipitation is about 400 mm, mostly falling as snow from September throughout May.
The land surface at the study site consists of unsorted, circular mud boils with diameters of about 1 m and a height of about 0.15 m.Their centers are bare and low vegetation grows between them (Fig. 2).In this study, we use daily averaged data from temperature sensors and time-domain reflectometry (TDR) probes that have been installed in the center of one of these mud boils at depths as shown in Fig. 5.The accuracy of the temperature sensors was estimated as 0.02 • C over the whole temperature range observed, while the TDR measurements return soil volumetric water contents with an estimated error of 0.02 to 0.05.There is an additional sensor pair in 0.25 m depth, but the measured temperatures start drifting, so it was omitted here.The temperature sensor in 0.77 m also starts drifting by 0.1 to 0.2 • C in 2001, but has been included nevertheless due to the small extent of the deviation.We work with data starting on 14 September 1998 and ending on 2 December 2001.This data set includes the range presented by Roth and Boike (2001) but extends an additional year.
Within the instrumented mud boil the soil mainly consists of silty clay with few interspersed stones.Small coal lenses occur below 0.85 m.The average soil bulk density is 1.60 × 10 3 kg m −3 , porosity ranges from 0.36 to 0.5.

The Tianshuihai site
The Tianshuihai site is located at 35 • 24 N, 79 • 33 E, 4740 m a.s.l., in a vast plain framed by mountains (Fig. 3).Permafrost extends to depth between 75 and 120 m (Li and He, 1989).
The Tibetan Plateau contains the only large permafrost region on Earth that is exposed to strong daily fluctuations in the radiative forcing.This is due to its low latitude  (Li and He, 1989), precipitation is very low in this region.
The plain area around the study site is characterized by loamy ridges of several centimeters in height which alternate with patches consisting of gravel and sand (Fig. 3).At the top of the ridges characteristic cracks occur.They are filled with sand wedges and indicate considerable frost action at this site.There is no vegetation present in the region.
In this study, we use soil temperature and soil water content data from a profile of about 1.75 m depth which has been installed as part of the Tianshuihai climate station and reaches about 0.2 m into the permafrost.It is located in a patch of gravel and sand in the center of a trough which is framed by loamy rims.The upper ten to fifteen centimeters consist of fine sand and some gravel, followed by a layer of sand and gravel that extends down to about one meter depth, which is then underlain by a layer of loam with interspersed patches of gravel.Porosities, densities, heat capacities and hydraulic and thermal properties were estimated from the recorded data or by numerical simulations.
Daily averages of soil temperatures and volumetric water contents, which were measured between 25 March 2008 and 28 December 2009, were used in this study (Fig. 6).The temperature sensors were calibrated with a relative accuracy of 0.02 to 0.03 • C, while for the TDR probes we assume an error of 0.02 to 0.05 in liquid water content.

Model description
To describe the dynamics of the active layer, we make several simplifications: we neglect mechanical processes and concentrate on heat transfer, examining heat conduction and convection by liquid water and water vapor.For this, we assume a one-dimensional, initially uniform soil with an incompressible, rigid matrix.
The model is based on a modified version of Richards' equation coupled with a heat equation based on Fourier's law.To simulate soil freezing, a parameterized freezing curve was included into Richards' equation, and terms for convective heat transport by liquid water as well as a parameterized formulation for water vapor flow were introduced into the equations.

Water movement
We model the one-dimensional hydraulic dynamics based on the conservation of mass in a molar formulation, where c w is the molar concentration of water with respect to the entire volume (mol m −3 ), given by c w = θ l ν l + θ i ν i .j w is the total water flux (mol m −2 s −1 ), including liquid water and water vapor.θ l and θ i are the volumetric contents of liquid water and ice, respectively, ν l = ρ l M w and ν i = ρ i M w the molar densities (mol m −3 ), M w is the molar mass of water (kg mol −1 ), and ρ l and ρ i are the mass densities of liquid water and ice, both set to 1000 kg m −3 .This assumption is necessary, because without it the change of the water density upon freezing would lead to unrealistically large gauge pressures that cannot be converted into an expansion of the soil matrix due to the lack of a mechanical model.
The water flux is given by Buckingham-Darcy's law, where K l is the liquid phase conductivity (m 3 s kg −1 ), T the temperature ( • C), g the gravitational acceleration (m 2 s −1 ) and p w the water pressure (Pa).Following Guymon and Luthin (1974) we neglect the ice pressure, assuming that both liquid water and ice pressure are given by p w .
To relate the total water content θ w to the water pressure p w , we use the soil water retention curve proposed by van Genuchten (1980): where φ is the porosity, θ r the residual water content, p atm the atmospheric pressure (Pa), m=1− 1 n and α (Pa −1 ) and n (-) are scaling and shape parameters, respectively.
Several approaches exist to obtain the latent heat contributions to the apparent heat capacity in dependence of the total water content θ w and the temperature T .Hansson et al. (2004) use the generalized Clapeyron equation to relate the water pressure p w to the temperature T and model the coexistence of the liquid and solid phases by introducing an empirical relationship for the apparent heat capacity, while Hinzman et al. (1998) apply a parameterized exponential function to obtain the latent heat content of the soil in dependence of the temperature.For the study sites that are discussed in this paper we use the following empirical relationship to obtain the ice content from the total water content θ w and the temperature T : The parameters a, b, c and d can be used to fit this curve into the saturated branches of the soil freezing characteristic, the function that relates the liquid water content θ l to the temperature T .In order to ensure that no ice is present at 0 • C we require b= a φ−d , which leaves three free parameters for the ice content.
To obtain the liquid phase conductivity we use the model proposed by Mualem (1976) combined with the van Genuchten parameterization (Eq.3), where κ is the hydraulic permeability (m 2 ) and µ l the dynamic viscosity of liquid water (Pa s).
In the above formulation we assume that the decrease of the liquid water content due to freezing is equivalent to a reduction due to drainage, which implies that in the unsaturated zone ice will form preferentially at the gas-liquid interface.

Heat movement
The heat equations are based on the principle of energy conservation, which reads with the thermal energy density E (J m −3 ) and the energy flux j e (W m −2 ).
The Cryosphere, 5, 741-757, 2011 Sensible heat of liquid water, ice, gas and soil matrix as well as latent heat contribute to the total thermal energy, so if we define the unfrozen soil at 0 • C as reference state, we obtain C i , C l and C g are the molar heat capacities of ice, liquid water and gas (J mol −1 K −1 ), C soil is the volumetric heat capacity of the soil matrix (J m −3 K −1 ), L f is the latent heat of freezing (J mol −1 ) and θ g is the volume fraction occupied by gas, given by θ g = φ − θ w .The molar density of the gas ν g (mol m −3 ) is given by the ideal gas law as ν g = p atm R(T +273.15K) , with the universal gas constant R (J mol −1 K −1 ).
In Eq. ( 7) we included the latent heat of vaporization in the constant effective heat capacity C g of the gas, assuming that the amount of water vapor is proportional to the temperature T .As described in the next section, this is not correct, but due to the small amount of water vapor that is present at typical temperatures in permafrost soils, the approximation seems reasonable at this place.The integral in Eq. ( 7) can be evaluated using Eq. ( 4).
When energy transport by water vapor is not considered, the energy flux can be written as Following de Vries (1975), the heat conductivity K h (W m −1 K −1 ) is parameterized as with the form factors k ι (-), and the sums ranging over all constituents of the soil, namely liquid water, ice, gas and soil matrix.The k ι describe the ratio of the average temperature gradient in particles of type ι to the average temperature gradient in the soil.If the particles are assumed to be independent spheres, the form factors are given by de Vries ( 1975) as assuming that the soil matrix forms a continuous background medium.As Ippisch et al. ( 1998) have shown, this approximation gives good results when compared to simulations that consider the structure of a soil sample explicitly.

Vapor movement
We implemented water vapor transport in a model following Bittelli et al. (2008) that describes water movement and latent heat transport within the soil by water vapor diffusion, including flow that is driven by vapor gradients as well as thermally driven flow.Also, we only consider thermal energy transport by latent heat and neglect sensible contributions, because the latent heat is two to three orders of magnitude larger than the sensible heat for the relevant temperatures.
When we include vapor transport, Eq. ( 2) is changed to while Eq. ( 8) becomes with the latent heat of vaporization L v (J mol −1 ), and the molar vapor flux j g (mol m −2 s −1 ) given by accounting for vapor flow by concentration gradients (first term) and thermally driven flow (second term).
The molar concentration of water vapor c v (mol m −3 ) is given by the Magnus formula (Murray, 1967) as with the relative humidity h r given by Philip ( 1957) Following Bittelli et al. (2008), the water vapor diffusivity is parameterized as with D 0 = 2.12 × 10 −5 m 2 s −1 .Using the equation of Buck (1981) for the saturation vapor pressure, the slope of the saturation vapor concentration s (mol m −3 K −1 ) is given by Bittelli et al. (2008) as

Numerical solution
We solved the set of two coupled, parabolic partial differential equations for the state variables p w and T in one dimension using the finite-element software package COM-SOL Multiphysics (2008).The equations are discretized with a Galerkin method on a regular grid with distances between the nodes of 0.01 m and Lagrange polynomials of second order as test and basis functions for both state variables.For the time stepping we use IDA, which is an implicit solver that uses variable-order variable-step-size backward differentiation formulas (BDF), and is described in detail in www.the-cryosphere.net/5/741/2011/The Cryosphere, 5, 741-757, 2011 Hindmarsh (2005).We work with interpolating polynomials of order one to five with a maximum time step of 5 min.
As nonlinear solver COMSOL Multiphysics ( 2008) uses an affine invariant form of the damped Newton method, which is described in Deuflhard (1974).With a maximum of four iterations and an update of the Jacobian on every iteration, this solver calls the linear solver UMFPACK (Davis, 2004), which is a direct solver for linear systems of equations.It uses the nonsymmetric-pattern multifrontal method and direct factorization of the sparse matrix A into lower and upper triangular matrices to solve the the equation Ax = b for the vector of unknowns x, which in our case contains the coefficients that describe T and p w in dependence of the basis functions.
Solutions of the Richards' equation were compared with solutions done by SWMS ( Šimůnek et al., 1995), and solutions of the heat equation including freezing were compared to analytical solutions to Stefan's problem (Carslaw and Jaeger, 1959), yielding good agreement.As the finiteelement method does not inherently fulfill the conservation laws, we calculated the numerical mass and energy defects by comparing the temporally integrated boundary fluxes with the water content and energy of the system.These defects were then divided by cumulative rain and net radiation, as these meteorological variables give an intuitive estimate for the boundary fluxes.The numerical mean daily mass defect stays below 0.7 % of the average daily rain per square meter for the Bayelva scenarios, and the mean daily energy defect does not exceed 0.3/0.1 % of the average daily absolute net radiation for the Bayelva/Tianshuihai simulations.As no hydrological simulations were run for the Tianshuihai site, no mass balance is given there.

Model set up
The formulation of the physical basis of our model representation was given in the previous section.An equally important part of a model encompasses the specific values for the chosen parameterizations on the one hand, and the model's initial conditions and subsequent forcing on the other.

Material properties
Specific values for the chosen parameterizations are necessarily more heuristic than the general forms of the parameterizations, which in turn are more heuristic than the differential equations that represent the dynamics.The reason for this is that the former depend on details of the material that are in general experimentally inaccessible, e.g. the geometry of the pore space or of the ice-and water-phase, while the latter are firmly rooted in physical principles.A bruteforce approach in such a situation is to numerically invert the available data for the desired parameters.Due to the large number of parameters and the notoriously ill-posed nature of the problem, we chose an alternative way that relies heavily on expert judgment.For any parameter, or set of parameters, we choose some initial guess from previous analyses of the data, from available information from comparable sites or from some general insight.These chosen values are then used for a numerical simulation of the experiment and adjusted manually, sometimes in conjunction with partial inversion, until the agreement between simulated and measured quantities was deemed sufficient.Such an exercise may lead to quite extensive multi-parameter analyses when structural deviations between simulation and data are encountered.An example for this are the different 0 • C-contour lines shown in Figs. 9 and 10 and discussed in Sect.5.2.
Before describing the steps to estimate some of the parameters in detail, we emphasize that such a non-objective approach always yields a possible set of parameters, but hardly ever an optimal one.Hence, the best achievable agreement between model and data is at least as good as the one obtained through our heuristic approach, probably better.Thus, the quantification of the impact of processes that are not represented by the model is an upper limit in the sense that in a model with optimal parameters, the neglected processes would be less important than estimated here.
For the Bayelva site, initial guesses of thermal and hydraulic parameters are based on previous analyses by Roth and Boike (2001), the resulting parameters for both sites are given in Table 1.To obtain the parameters of the freezing characteristics, the measured freezing curves were plotted (Fig. 4), and a manual fit with θ w = φ was put into measurements from probes situated in the water saturated section of the profile.The values were then varied by some percent such that the simulated temperatures and water contents within the whole profile fit the measurements in a better way.The solid line in Fig. 4 shows the optimized function.
For the Bayelva site we used homogeneous thermal properties for the complete soil profile.The heat conductivity of the soil material K h soil was obtained by manual optimization, while the heat capacity C soil was taken from Roth and Boike (2001).The quantity that determines heat conduction is the thermal diffusivity D h (m 2 s −1 ), which is defined as the ratio K h /C eff , where C eff (J m −3 K −1 ) is the effective soil heat capacity that can be obtained by averaging the heat capacities of the components (solid matrix, air, water and ice), weighed by their respective volume fractions.
At the Tianshuihai site, test runs of the simulation with a homogeneous medium provided unsatisfactory results, with temperature deviations between measured and simulated data of more than 3 • C. As the soil shows a quite pronounced layering and its texture strongly changes with depth, the modeled domain was separated into three layers with boundaries at depths of 0.13 m and 0.97 m, all with different heat conductivities K h soil .All other parameters were kept constant over the spatial domain.
For this site, parameters were adapted using an approach that alternates automatic estimation for adjustment of the   three different heat conductivities K h soil in the different layers and manual optimization of all other parameters.The automatic estimation was done by rastering the parameter space, running forward simulations for a discrete set of values for each parameter, with an increased maximum time step of 3 h for computational speed.The optimal parameter sets for both sites were determined by the minimal sum of squared differences between the measured and the simulated values.This method accounts for the observation that a homogeneous soil did not provide satisfactory results, while an inversion of all parameters was not feasible computationally.The heat capacity will only influence the values we obtain for K h soil , but not the simulated temperature distribution, because the heat equation only depends on the thermal diffusivity D h .Resulting values are given in Table 1.

Simulation scenarios
To examine the influence of heat convection by liquid water and water vapor compared to heat conduction, three scenarios with different complexity of the model of the Bayelva site were run as given in Table 2.The most complex scenario (i) includes conductive as well as convective heat transport by liquid water and water vapor.Scenario (ii) neglects water vapor transport and scenario (iii) considers conductive heat transport only.Each run was performed with the same set of parameters, which was manually optimized for the most complex scenario (i).
Unlike at the Bayelva site, the data from the Tianshuihai site has not been analyzed in previous studies.Also, as no texture analysis was performed, most material properties had www.the-cryosphere.net/5/741/2011/The Cryosphere, 5, 741-757, 2011 (iii) Conductive heat transport only, assuming variable water saturation in space but not in time.
to be estimated from textural information based on field observations.As we will discuss in Sect.5.4, hydraulic processes were switched off, so only scenario (iii) was analyzed for the Tianshuihai site.

Boundary and initial conditions
The model is driven at the upper boundary by imposing the values of soil temperature and water content measured by the respective uppermost sensors.This choice means that all surface processes that transport water and sensible heat across the soil boundary are captured, as long as they do not bypass the sensors, e.g. as flux through macropores.This choice also does not permit to represent transport of latent heat in conjunction with vapor flux across the upper boundary.
For soil temperatures below 0 • C we used the last soil water content from temperature measurements above freezing as total water content for the upper forcing.This was necessary because we could not use the TDR measurements for the total soil water content in winter, as they only return the liquid water content, while we needed the total (frozen and unfrozen) water content as upper boundary condition.Due to the low hydraulic conductivity and the presumed low vapor fluxes in frozen ground, this assumption should be valid, only when rain occurs shortly after freezing starts, water infiltration might be underestimated.
Because of the permafrost below the modeled domain, we assumed that no water moves through the lower boundary.For modeling the thermal lower boundary condition, in natural soils there is no physically obvious choice anywhere close to the modeled domain.As we do have temperature measurements within the profile, we used these as the lower boundary condition.While it would be desirable to use a lower boundary in greater depth where temperatures stay constant, our approach enabled us to calibrate the model by adjusting the soil hydraulic and thermal parameters, and allowed to identify structural deviations between the simulation and the measured data.
As initial conditions for the heat equation we applied a linear interpolation between the measured temperature values.For the water contents, we apply a temporal average of measured soil water contents in summer, because the model is initialized when at least parts of the soil are frozen.At greater depth, where the soil does not thaw completely in summer and thus total water contents could not be measured, we used measurements from the TDR probe right above the permafrost table, assuming that the total soil water content stays constant below this depth.Between the sensors, linearly interpolated values were used.As the water content in general is not a continuous quantity, these values might not actually represent the stationary problem, but should be close to a solution if the model represents reality well, and thus equilibrate quickly.

Results
The results for the most complex scenario (i) of the simulation of the Bayelva site with conductive and convective heat transport through liquid water and water vapor is shown in Fig. 5, characteristic features of all scenarios of the Bayelva simulation are provided in Fig. 9. Results of the Tianshuihai simulation are provided in Figs. 6 and 10.Deviations between observed and simulated temperatures are quantified in Table 3, excerpts of the results for the first summer are shown as line plots in Figs.7 and 8.We first get an overview of the results and then focus on some particular aspects.

Overview
We first notice that the qualitative agreement between measured and simulated quantities is very good for all simulations.Typical features found in permafrost soils are represented correctly: the low temperatures during the cold periods in winter, the warming periods that end with sharp thawing fronts propagating into the soils in late spring, the thawed periods, and finally the zero curtains, isothermal plateaus and the subsequent cooling periods.There is furthermore a good quantitative agreement in that corresponding contourlines for temperature and liquid water content for the data and the simulations are near to each other.This is also testified by the difference plots in Figs. 5 and 6.
As a backdrop for assessing the agreement between simulation and observed reality, we recall the results of some other studies.For a site near Toolik Lake, Alaska, Kane et al.
The Cryosphere, 5, 741-757, 2011  (1991) and Hinzman et al. (1998) reported deviations of up to 1 • C and Wang et al. (2010) obtained deviations of up to almost 2 • C for a site in eastern Tibet.Both report the largest deviations during spring where a small error in the timing of the thawing front leads to large discrepancies in tempera-tures as well as liquid water content.On the other hand, using a semi-analytical pure-conduction model for the Bayelva site, Roth and Boike (2001) reported deviations as low as ±0.2 • C. Note however, that their analysis was restricted to periods and depth intervals without freezing or thawing.Looking at the results for the Bayelva site obtained in the current study, we first recognize that the simulations tend to underestimate temperatures during the summer and to slightly overestimate them during the winter (Table 3 and Fig. 5).The somewhat worse simulation for the second winter is attributed to the thinner snow cover which led to lower soil temperatures and thus, assuming a similar relative error, to a larger absolute deviation.The deviations during the summer will be discussed in Sect.6.3.
Turning to the Tianshuihai site, we first comment that a satisfactory simulation of the temperatures was not possible with a uniform soil profile.Hence a layered model was set up, guided by the profile description available for the site.This led to the much better performance of the model that is illustrated by Fig. 6.The corresponding values for the heat conductivity K h soil are shown in Table 1.The colors represent: black: measured data, blue: conduction with convection of liquid water and water vapor, scenario (i), green: conduction with convection of liquid water, scenario (ii), red: conduction with constant water distribution, scenario (iii).Blue (i) and green (ii) overlap almost entirely, so only the green line can be seen clearly.The deviations between simulations and measurements at the Tianshuihai site appear structurally different from those at the Bayelva site.During the summer, there are short periods for which simulated temperatures are too high and others for which they are too low.We hypothesize that this is caused by episodic exchange of water vapor across the soil surface and to depths below the topmost temperature sensor.During the winter, simulated temperatures tend to be lower than measured ones, which is in contrast to the Bayelva site.We emphasize, however, that deviations are low throughout, some 0.28 • C on average.
In closing this overview, we notice that our choice of boundary conditions -driving the processes with state variables -results in boundary fluxes that are highly dependent on the soil hydraulic and thermal properties.Since we obtain good agreements between measured and simulated data, this indicates that the simulated fluxes and the chosen material properties are good representations of the corresponding reality.We comment, however, that a direct experimental verification of this statement cannot be obtained with current technology.

Heat entry in early summer
The thaw depth estimate is very good in all scenarios, which is to be expected due to the nearby temperature forcing at the lower boundary.Still, we observe in Figs. 9 and 10 that the simulated thawing front at greater depths lags behind the measured one by up to two weeks, yet propagates about five centimeters deeper than what is deduced from the measurements.Since this occurs at both sites and for all three simulated scenarios, possible explanations for this effect will be discussed in Sect.6.3.
These two observations -lagging thawing front yet greater thawing depth -fit together well: the thawing front may be considered as a fixed-temperature boundary for heat conduction in the thawed regime.If this boundary is not as deep in the simulation as in reality, the temperatures above it will be underestimated.
This effect may be observed in Figs.7 and 8: when the thawing front passes the sensor, there is a fast increase in liquid water content, which at greater depths occurs later in the simulation than in the measured data (Fig. 7).At the Bayelva  site, the chosen parameters lead to an overcompensation in 0.41 m depth, where the thawing front arrives too early after moving too fast through the dryer upper soil.Some weeks later, when the thawing front has moved deeper, the underestimation of the temperature becomes apparent.

Heat transport processes
In order to assess the relative importance of heat conduction and convection, we consider the mean heat fluxes in the observed profile for the Bayelva site (Fig. 11).It reveals that heat transport is dominated by conduction in all simulated scenarios and for all times.Indeed, convection of liquid water carries just 3.2 % of the total heat flux during the unfrozen period, and is obviously negligible in the frozen state.Contributions by water vapor can be observed in summer as well as in winter, but they are very low with 0.01 % during summer and 0.03 % throughout the year.In the upper 0.05 m, the corresponding numbers are larger, 0.03 % during summer and 0.2 % throughout the year.They are not reliable, however, since our model formulation prevents vapor exchange across the upper boundary.
A similar assessment for the Tianshuihai site is not possible, since we found that convective transport is not required to simulate the observed dynamics.

Hydraulic effects
In the scenarios that include hydraulic processes (Bayelva, scenarios (i) and (ii), Fig. 5), the water movement in summer is represented fairly well, the liquid soil water contents deviate by less than 0.05 during most of the time.For short periods during the phase transitions, the error increases to at most 0.25, corresponding to energy defects of up to 10 MJ m −2 , because the time of freezing or thawing does not fit exactly.In contrast, at the Tianshuihai site (Fig. 6) simulation of water flow by forcing the model with the water content measured at the top-most TDR probe was not possible for two main reasons: (i) the surface is almost always very dry, hence the functions K l (θ w ) and θ w (p w ) are very steep and, as a consequence, the simulated water flux is highly sensitive to small errors in the parameterization of these functions and in the measurements of θ w .Indeed, the observed higher water contents in the top-most layer as a result of small precipitation events led to massive infiltration events in the numerical simulation, which were not observed in the measurements.(ii) Infiltrating water, even though it may pass the top-most TDR-probe, will evaporate quickly, since potential evaporation is about 20 to 50 times higher than precipitation (Gasse et al., 1991).However, vapor transport across the upper boundary is not included in the numerical model due to a lack of relevant data.
An alternative approach to force the hydraulic model would be to use the flux as boundary condition, incorporating rain measurements and estimates of evapotranspiration.However, especially the latter is a topic of active research and so far no reliable method is available to measure it accurately over long periods.
Except for two rain events in the summer of 2008, the assumption that convective transport does not play an important role at the Tianshuihai site is supported by the TDR measurements, which show an almost constant water content within the soil with a water table that varies only by a few centimeters during the unfrozen period.Hence, the scenarios (i) and (ii) of the Tianshuihai simulation were not analyzed any further.

Discussion
We modeled processes in the active layer at two sites that show rather different surface as well as subsurface structures and that are subject to very different atmospheric forcing.We were able to reproduce the thermal dynamics at both sites.Relevant processes are sufficiently similar to be represented by the same mathematical model, although the inclusion of site-specific knowledge such as atmospheric conditions, surface structure and subsurface texture is crucial.Including this knowledge together with measured temperature and TDR data all required parameters could be estimated.Differences in thaw depth were well below 0.01 m.Hence, while the estimated values of D h differ significantly, the impact on simulated temperatures and thaw depths is negligible.
In the following, we will discuss some specific requirements for the good performance of the model.

Hydraulics in the frozen regime
Water movement in frozen soils is a very complicated process which is not yet fully understood.Issues at the macroscopic scale involve the modifications of the soil water characteristic and of the hydraulic conductivity function by ice and the appropriate thermodynamic potentials for the different phases.For operational models, they are typically resolved in a pragmatic way by introducing a heuristic impedance factor (e.g.Hansson et al., 2004) and by postulating relations between the chemical potentials in liquid and frozen water (e.g.Christoffersen and Tulaczyk, 2003).In a first approach, we implemented a rather crude representation of the complicated physics and found that convective transport is of minor importance for the thermal dynamics at the sites considered.Hence we did not attempt to improve the chosen representation.

Calibrated thermal conductivities
From Table 1 we notice that the values estimated for the thermal conductivity K h soil are at extreme ends of what is expected for the given materials.This can be understood as a consequence of our choice for the parameterization of the thermal conductivity.In the de Vries model, there are two different classes of adjustments: (i) the conductivity of the pure constituents and (ii) their geometry.We chose to only adjust the conductivities and to absorb into them all effects of geometries that deviate from the implicitly presumed model of spheres widely dispersed in a uniform background.The low value of K h soil in the upper layer at the Tianshuihai site is found in a sandy soil with very low water contents of 0.05 to 0.08 with peaks reaching 0.16 during rain events.There is almost no water present at the contact interfaces of the soil particles, so the assumption of the de Vries model about the background medium does not hold.Also the high value in the lower layer of the Tianshuihai site and in the whole profile at the Bayelva site may be attributed to geometry.In effect we observe that it is important to include strongly differing thermal diffusivities even in a small layer.

Temperature and thaw depth deviations in summer
At both sites, the most significant deviation between the measured and the simulated data is the underestimated heat input in early summer.This behaviour has also been observed in previous studies for different field sites, e.g.Kane et al. (1991) for an Alaskan site and Wang et al. (2010) for a site in eastern Tibet.However, no detailed discussion was given there.
Possible explanations are that the responsible process is not detected by the uppermost temperature sensor (such as transport through macropores) or some other non-diffusive process.In addition, we need explanations for a very similar effect at both sites, although the deviations are about twice as strong at Bayelva than at Tianshuihai.
In the fall of 1999, strong rain events occurred at the Bayelva site shortly after the onset of freezing.This led to infiltration into an already frozen ground in the simulation, and thus to an overestimation of the ice content in winter.This might be one of the reasons for the underestimation of the thaw depth in the following summer in scenarios (i) and (ii) and the resulting underestimation of temperatures in the unfrozen regime.However, as the temperature is also underestimated in the heat conduction scenario (iii) at both sites in all summers, there must be another effect involved in causing the observed deviations.In the following sections, we will discuss possible explanations.
While we could not simulate water migration for the Tianshuihai site successfully, some effects of the rain can be observed in the temperature estimates: In the summer of 2008, two distinct infiltration events can be identified in the TDR data (Fig. 6, 8 July and 5 to 16 August), and during these times the error in the temperature estimation increases significantly, such that temperatures are underestimated by the simulation.This hints at either an increase of the effective thermal conductivity by the additional water or convection of rain water that is warmer than the soil.

Soil-freezing characteristic
As mentioned in Sect.4, the soil freezing characteristic was optimized manually to match the observed temperature and water content measurements rather than the observed freezing curve itself.However, the overall influence of the steepness of the freezing characteristic is not very large: on both sites, changes in the steepness -parameters a and c in Eq. ( 4) -of 50 % modify the temperature distribution by less than 0.15 • C. The sensitivity towards changes in the fraction of the total water content subject to freezing and thawing (parameter d) is higher: at the Tianshuihai/Bayelva site, an increase by 50 % in d (corresponding to a decrease in freezing water) increases the thaw depth in early summer by about 7/5 cm.At the same time this shortens the duration of the isothermal plateau by about 6/3 days.As the length of the isothermal plateau is already underestimated by up to several days (the exact number depending on the site and the scenario, see Figs. 9 and 10), neither a different choice of the shape of the freezing characteristic nor the total amount of water in the profile can improve the results significantly.
The linear scaling of the freezing characteristic with the total water content as implemented in Eq. ( 4) is an approximation that does not consider the layer boundaries.In Fig. 8 for example one can see that the difference in the measured liquid water content in winter differs by a factor of about three between the two lowermost probes, which are only 10 cm apart.The simulated curves on the other hand overlap almost entirely, because the almost identical total water content implies a very similar freezing characteristic.This could only be avoided by using different freezing characteristics for each layer.However, we find that our approximation yields sufficiently good results in the thermal regime, where the different latent heat contributions of the layers average out over the depth.

Snow melt infiltration through macropores
At the Tianshuihai site, precipitation is very low and the infrequently occuring thin snow covers usually disappear on the same day, while most of the snow likely sublimates.For this reason, snow melt infiltration is only significant at the Bayelva site.
Snow melt infiltration is captured by the upper boundary if it does not pass below the uppermost sensors through macropores.Notice however, that also snow melt that infiltrates through macropores cannot alter the thaw depth in spring: the melted snow has a temperature of 0 • C, so when it reaches the thawing front and deposits latent heat there, it can only do so by freezing itself.Therefore, the total amount of water that has to be thawed by other processes and thus the thaw depth does not change.

Dry cracking and water vapor condensation
One possible explanation for the underestimation of the heat entry in early summer would be dry cracking, which was observed at the Bayelva site on top of the mud boils.Within these cracks, vapor convection could transport large amounts of heat from the surface below the upper temperature sensor.In contrast, at the Tianshuihai site the only observed potentially relevant structures were larger sand-filled wedges that reached to a depth of about 0.3 m.Additional smaller cracks are unlikely, considering the coarse-grained surface sediments.
As we have shown for the Bayelva site, the amount of heat transported by convection of water vapor within the soil is two orders of magnitude smaller than the heat convected by liquid water.However, the way vapor transport was implemented in this model, vapor flow across the upper boundary is not included.Hence, vapor that diffuses from the atmosphere through cracks past the uppermost temperature probe into the soil and deposits its latent heat there constitutes a possible heat input that is not captured by the model.Due to the temperature gradient between the atmosphere and the soil surface, this effect would be strongest during early summer, when the temperature difference is high.This is also the period for which the thaw depth is underestimated the most.

Mechanical soil expansion
In natural soils, the different densities of liquid and frozen water lead to an expansion of the soil upon freezing.As these mechanical effects are not implemented in the model, setting the densities to different values leads to unphysical effects.Exploring this numerically, we observe a huge liquid phase gauge pressure that causes a rise of about 0.2 m in the water table during the isothermal phase and some small fluctuations during the frozen period.However, the water content distributions in summer as well as the temperature distributions do not change significantly in either scenario, because the total amount of water does not change.
In real soils, the expansion of the soil matrix may lead to different distances between the probes in winter, because the probes have been installed in summer, and thus the recorded depths describe positions in the unfrozen soil.If we now adjust the parameters such that the temperatures in winter are represented accurately, we will underestimate the thermal diffusivity and thus also the heat conductivity of the soil K h soil in summer, because the distances in the simulated soil are smaller than in reality.This would lead to an underestimation of the temperatures as soon as they start to drop.
However, this effect does not seem to be strong enough.Assuming that the volume expansion of ice, some 9 %, is directly transformed into a corresponding expansion of the soil and in vertical direction only, an average water content of 0.3 leads to an expansion of vertical distances by some 3 %.Roth and Boike (2001) state that if we want to project a given temperature distribution into increasing depths we would have to increase the effective thermal diffusivity with the square of the depth to obtain the same distribution.In our case we would therefore have to increase the diffusivity by 6 % to represent the temperature distributions in summer.While test runs of the simulations showed that with this value the error of the temperature in summer improves, but only by about 0.03 • C.This is not sufficient to explain the observed differences.

Three-dimensional effects
The restriction of our model to one dimension is an approximation that neglects larger-scale processes and especially lateral fluxes.However, the topography of the selected sites suggests that these fluxes do not play an important role: at the Tianshuihai site, the surroundding area is flat and we presume similar thermal properties throughout the plain.Largescale ground water flux may be present at the site, but due www.the-cryosphere.net/5/741/2011/The Cryosphere, 5, 741-757, 2011

Fig. 3 .
Fig. 3. Overview over the sedimentary plain at the Tianshuihai station.

Fig. 3 .
Fig. 3. Overview over the sedimentary plain at the Tianshuihai station.
, which is equal to that of Crete in Europe or of Albuquerque in the USA.Measurements recorded at the Tianshuihai climate station, which was installed in September 2007, ranged from −100 W m −2 at night to more than 800 W m −2 in June 2008, while daily averages mostly stay between 0 and 150 W m −2 .Measured air temperatures ranged from −8 to 21 • C in July 2008 and from −38 to −5 • C in January 2009.With about 24 mm yr −1 of rain

Fig. 5 .
Fig. 5. Measured data and simulation results for the Bayelva simulation with conductive and convective transport by liquid w vapor (scenario (i)).Top two: Measured and simulated soil temperatures, with solid lines every 2 • C, gray lines every 0.2 • C b and 2 • C and a magenta line for 0 • C. 3rd from top: Simulated minus measured temperatures, with the measured temperature for cross referencing.4th and 5th from top: Measured and simulated liquid soil water contents, with the lines from the tem repeated.Bottom: Simulated minus measured water contents, with the measured temperature lines repeated.Horizontal line sensor positions.

Fig. 5 .
Fig.5.Measured data and simulation results for the Bayelva simulation with conductive and convective transport by liquid water and water vapor (scenario (i)).Top two: measured and simulated soil temperatures, with solid lines every 2 • C, gray lines every 0.2 • C between −2 • C and 2 • C and a magenta line for 0 • C. 3rd from top: simulated minus measured temperatures, with the measured temperature lines repeated for cross referencing.4th and 5th from top: measured and simulated liquid soil water contents, with the lines from the temperature plots repeated.Bottom: simulated minus measured water contents, with the measured temperature lines repeated.Horizontal lines represent the sensor positions.

Fig. 6 .
Fig.6.Measured data and simulation results for the conduction scenario (iii) of the Tianshuihai simulation.Graphical representa to that of Fig.5.Arrows mark the rain events discussed in Sect.6.3.

Fig. 6 .
Fig.6.Measured data and simulation results for the conduction scenario (iii) of the Tianshuihai simulation.Graphical representation identical to that of Fig.5.Arrowheads mark the rain events discussed in Sect.6.3.

Fig. 7 .
Fig.7.Measured (solid)  and simulated (dashed) temperature and liquid water content for the Bayelva simulation with conductive and convective transport by liquid water and water vapor (scenario (i)).The data show the period from 22 May 1999 to 27 January 2000.The probes at 1.14 m depth are situated right at the edge of the thaw depth, here a slight misestimate of the thaw depth can be seen clearly in the liquid water content for days 340 to 380.

Fig. 7 .Fig. 8 .
Fig.7.Measured (solid)  and simulated (dashed) temperature and liquid water content for the Bayelva simulation with conductive and convective transport by liquid water and water vapor (scenario (i)).The data show the period from 22 May 1999 to 27 January 2000.The probes at 1.14 m depth are situated right at the edge of the thaw depth, here a slight misestimate of the thaw depth can be seen clearly in the liquid water content for days 340 to 380.

Fig. 8 .Fig. 9 .Fig. 10 .
Fig. 8. Measured (solid) and simulated (dashed) temperature and liquid water content for the conduction scenario (iii) of the Tianshuihai simulation.The data show the period from 25 March 2008 to 19 January 2009.

Fig. 9 .
Fig. 9. Simulation results for the different scenarios at the Bayelva site in the summer of 1999.Shown are the 0 • C and the −1 • C isotherms.The colors represent: black: measured data, blue: conduction with convection of liquid water and water vapor, scenario (i), green: conduction with convection of liquid water, scenario (ii), red: conduction with constant water distribution, scenario (iii).Blue (i) and green (ii) overlap almost entirely, so only the green line can be seen clearly.

Fig. 9 .Fig. 10 .
Fig.9.Simulation results for the different scenarios at the Bayelva site in the summer of 1999.Shown are the 0 • C and the −1 • C isotherms.The colors represent: Black: Measured data, blue: Conduction with convection of liquid water and water vapor, scenario (i), green: Conduction with convection of liquid water, scenario (ii), red: Conduction with constant water distribution, scenario (iii).Blue (i) and green (ii) overlap almost entirely, so only the green line can be seen clearly.

Fig. 10 .
Fig. 10.Simulation results for the Tianshuihai site in the summer of 2008.Shown are the 0 • C and the −1 • C isotherms.The colors represent: black: measured data, red: conduction with a water distribution that is constant in time, scenario (iii).

Fig. 11 .
Fig. 11.Mean energy flux at the Bayelva site averaged over the depth of the profile.Blue: conductive, green: convective by liquid water, red: convective by water vapor.

Fig. 11 .
Fig. 11.Mean energy flux at the Bayelva site averaged over the depth of the profile.Blue: conductive, green: convective by liquid water, red: convective by water vapor.
For www.the-cryosphere.net/5/741/2011/The Cryosphere, 5, 741-757, 2011 J. Weismüller et al.: Modeling thermal dynamics of active layers the thermal diffusivity D h in the frozen regime in 1998/1999 at the Bayelva site, Roth and Boike (2001) obtained 0.8 × 10 −6 m 2 s −1 .Values used in this work change over time and range from 1 × 10 −6 m 2 s −1 to 1.4 × 10 −6 m 2 s −1 , depending on the soil composition.For comparison with the previously obtained value, we artificially decreased D h by a factor of 1.5, leading to values between 0.67 × 10 −6 m 2 s −1 and 0.93×10 −6 m 2 s −1 .This increased the error by up to 0.06 • C in the first and third winter, and by up to 0.2 • C in the second.

Table 1 .
Soil parameters used in the simulations.As no convective flow has been modeled for the Tianshuihai site, no Mualem-van Genuchten parameters were necessary.

Table 2 .
Simulation scenarios.(i)Conductiveaswellasconvectiveheat transport by liquid water and water vapor.Using all processes that are described by Eqs.(1) to (17); water vapor only for transport within soil profile, no soil-atmosphere coupling is considered.