Modelling the spatial pattern of ground thaw in a small basin in the arctic tundra

of the frost table, which here acts as a relatively impermeable surface, are therefore crucial in determining the hillslope drainage rate. This work aims at understanding how the topography of the ground surface affects the spatial variability of the depth of thaw in a 1 km low-elevation arctic tundra basin with a fine resolution model that fully couples energy and water flow processes. The simulations indicate that the spa10


Introduction
Arctic tundra is a treeless terrain of continuous permafrost, with a relatively thin (i.e.0.05-1.0m) organic soil layer, overlying mineral sediment.This configuration has strong implications to the runoff and thaw processes (Quinton et al., 2005).The organic soil layer commonly includes peat (Tarnocai, 2004), and is characterized by low thermal conductivity, very high porosity of up to 0.95 (The Canadian System of Soil Classification, 1998), and, in particular, very high hydraulic conductivity of up to 10 mm s −1 (Quinton et al., 2008).This is opposed to the relatively high thermal conductivity and low hydraulic conductivity (between 10 ) of the underlying mineral sediment (Quinton et al., 2000).As a consequence, hillslope drainage is normally minimal on the surface, and occurs predominantly within the organic layer (Quinton and Marsh, 1999).
In these environments the ground below the frost table, defined as the zero-degree isotherm, is normally saturated with ice and a small fraction of unfrozen water.As a result, the frost table forms a relatively impermeable surface, which acts as the lower boundary of subsurface flow (Quinton et al., 2009;Hayashi et al., 2003).Its position, therefore, affects the hillslope drainage rate.This becomes even more important if it is considered that the hydraulic conductivity in the organic soil dramatically decreases with depth due to increasing degree of decomposition.In particular, it has been observed that the saturated hydraulic conductivity of peat decreases by 2 orders of magnitude with increasing depth between 0.1 and 0.2 m, regardless of the total thickness of the peat (Quinton et al., 2008).Therefore, the rate of subsurface flow is strongly dependent on the degree of thaw.Early in the thaw season, lateral flow is relatively high, since the thawing-front (i.e.frost table) is relatively close to the ground surface, where the hydraulic conductivity is relatively high.Later in the thaw season the lateral flow rate is substantially lower since the thawing-front occupies a deeper position, and, therefore, horizontal drainage occurs through the portion of organic soil where the hydraulic conductivity is lower (Quinton and Marsh, 1999).Spatial and temporal variations of ground thaw rates result in continuous changes to the frost table topography.As a result, the patterns of subsurface drainage, including the volume, timing and direction of flow, are highly variable over space and time.Knowledge of the frost table topography in relation to the surface topography is therefore critical to predicting runoff from organic-covered permafrost terrains.Abbey et al. (1978) found that the depth of thaw is closely related to the cumulated heat flux transferred into the ground from the atmosphere.This implies that the spatial variability of the degree of thaw is strongly affected by the spatial distribution of solar and terrestrial radiation, and turbulent fluxes of sensible and latent heat.Carey and Woo (2001) showed that north-facing slopes in a sub-arctic mountain site were characterized by greater near-surface wetness as a result of shallow thaw depths, while deeper subsurface pathways due to deeper thaw were more evident on south-facing slopes.However, the pattern of soil thaw is also influenced by the spatial pattern of snow disappearance (Quinton and Carey, 2008), which, in turn, is determined by the spatial distribution of the snow mass and melt energy at the snow surface.Given the highly insulating properties of the snow cover, it is commonly observed that appreciable soil thaw only starts once snow is completely removed (Quinton et al., 2009).
Subsurface water flow can also affect the thermal regime of the seasonally thawed (i.e.active) layer, and, therefore, affects the spatial variability of thaw.In areas of converging flow where water accumulates, the soil thermal conductivity, which increases with soil water content (Kane et al., 2001), is raised, producing greater thawing, even if the rate of energy transfer into the ground is offset by the energy loss due to enhanced evaporation of the relatively wet ground surface.
The purpose of this work is to understand how the topography of the ground surface affects the rate and spatial patterns of ground thaw in a catchment in the continuous permafrost region of arctic tundra using a fine resolution digital elevation model, limited field observations, and a model that fully couples the soil heat and water flow equations, represents the energy exchange with the atmosphere, and describes the subsurface water flow in its three-dimensional nature.By using a model that couples energy and mass flows it is possible to investigate the linkages between the spatial patterns of energy fluxes and runoff.Since the frost table is relatively impermeable, and the lateral hydraulic conductivity decreases with depth below the ground surface in a manner that we can describe with a simple analytical formula (Quinton et al., 2008), a map of the frost table topography is in effect one that indicates the horizontal hydraulic conductivity at the bottom of the thawed layer as well as the magnitude and direction of the hydraulic gradients at a given moment in time.
This work represents one of the first attempts to describe the temporal and spatial evolution of the depth of thaw with a fully coupled three-dimensional distributed model at a small scale.Existing models applied in permafrost environments are normally: (i) models applied at local, regional, and continental scale that integrate a one-dimensional form of the heat equation with phase change to predict the evolution of the depth of thaw, but neglect, or address with simple methods, the coupling with subsurface lateral water flow (e.g.Hinzman et al., 1998;Lawrence and Slater, 2005;Marchenko et al., 2008); (ii) hydrological models that are commonly applied at a large scale to predict river discharge, but do not accurately describe the coupling with the soil energy balance and the temporal evolution of the depth of thaw (e.g.Kuchment et al., 2000;Zhang et al., 2000); and (iii) models that very accurately describe and couple water and energy subsurface flow in frozen soil, but do not consider the heat flux exchanged with the atmosphere, given by the sum of net radiation and turbulent fluxes (McKenzie et al., 2006;Hanson et al., 2004;Daanen et al., 2007).

Study area
This study focuses on Siksik Creek (Fig. 1), a subcatchment of Trail Valley Creek (68 • N, 133 The ground surface is characterized by the presence of mineral earth hummocks, covering ∼ 50% of the ground surface, with diameters of between 0.4 to 1.0 m, and crests between 0.1 and 0.4 m above the surface of the inter-hummock area.The latter is composed of peat ranging in thickness between 0.2 m and 0.5 m (Quinton and Marsh, 1999), with an abrupt transition between the organic material and underlying mineral sediment.The upper 0.1 m of the peat is comprised of living vegetation rooted in weakly decomposed plant material.Below this, the peat is composed of moderately to very strongly decomposed plant material with plant structures much less distinct.
Other studies have been conducted at Siksik Creek.Quinton and Marsh (1998) and Quinton et al. (2000) studied the effect of hummocks on hillslope drainage.Quinton and Marsh (1999) performed systematic observations of the water table in several points and developed a conceptual framework of runoff generation.They showed that the peat thickness in the inter-hummock area is normally higher close to the stream channel.As such, the "near-stream area" buffers streamflow due to its relatively large water storage capacity.However, they defined only empirically the evolution of the frost table topography.This study, focusing in detail on the complex mechanisms controlling the spatial variability of the depth of thaw, represents a step forward in this limitation.

Methodology
The GEOtop model (Rigon et al., 2006) searches a coupled solution of the heat and water flow equations in the soil.As such, it is an appropriate tool to examine how the spatial pattern of energy relates to that of runoff in organic-covered permafrost terrains.Moreover, GEOtop has recently been improved with respect to the representation of surface energy balance (Endrizzi and Marsh, 2010), the soil freezing and thawing processes (Dall'Amico et al., 2010), and the numerical algorithms.The snow constitutes a further element of complexity that we will not consider here.Instead, we ran GEOtop over the summer and assumed an initial condition of uniform frost table at the ground surface at the end of snow melt, implicitly assuming that snow is removed instantaneously across the entire basin.The analysis will in particular focus on the end-of-summer frost table topography, as this is the cumulative product of all ground thawing mechanisms active throughout the thaw season.It is also expected that the influence on the frost table topography of the variable snow cover duration on the landscape is greatest immediately after snow removal, but diminishes with time, and is least by the end of summer.
Another element of complexity is given by the heterogeneity of the vegetation cover, which, for the above-mentioned reasons, will not be addressed in this work.The ground surface is assumed uncovered by shrub canopy.The living and lightly decomposed vegetation on top of the peat is described as a particular type of soil of low heat capacity and high porosity (Quinton and Gray, 2003).
Measurements of the thaw depth commenced once the ground became snow-free in May 1993, and continued until the end of August 1993.Measurements were focused in the inter-hummock zone at 3 plots in proximity to the mainstream channel (Fig. 1).

GEOtop
GEOtop is a grid-based model that seeks a numerical solution of both the heat and subsurface water flow equations in the soil.This characteristic is common among several Land Surface Models (e.g.CLM, Oleson et al., 2004;CLASS, Verseghy, 1991), but, differently from them, the equations are solved at small grid sizes (metres to tens of metres).The model domain consists of a soil volume of user-specified depth (typically of a few meters) from the ground surface, which is, in turn, defined by a Digital Elevation Model (DEM).The domain is discretized in a finite regular grid, given by the intersection of the DEM cells.The heat and subsurface water flow equations are then solved with finite differences schemes.While the heat equation is solved in a one-dimensional form in the grid cell columns normal to the surface, the equation describing the water flow in the soil (Richards equation) is solved in a fully three-dimensional way, which allows a complete coupled description of vertical and lateral flows.The surface water flow over land is also described, with a solution of a simplified version of the shallow water equations that neglects the inertia terms.The integration time step typically ranges from a few minutes to a few hours.The determination of the heat flux exchanged from the atmosphere to the ground surface (hereafter referred to as ground heat flux) is crucial, since it constitutes the upper boundary condition of the heat equation.The ground heat flux is given by the algebraic sum of net shortwave (solar) radiation, net longwave (terrestrial) radiation, and turbulent fluxes of sensible and latent heat.
The resolution methods of the heat and subsurface water flow equations as well as the parameterization of the surface energy flux are briefly illustrated below.

Heat equation
The heat equation is derived from the conservation of energy, and its 1-D form is written as ] the soil thermal conductivity.Equation ( 1) is integrated assuming that the boundary condition at the surface is given by the surface heat flux, and the boundary condition at the bottom is a fixed temperature at a certain depth (namely, the zero-amplitude annual soil temperature).The internal energy includes the effects of sensible and latent heat, and can be written as where C is the soil thermal capacity [J m −2 K −1 ], T f is the freezing temperature (273.15K), L f the thermal heat of fusion (3.34 × 10 6 J kg −1 ), ρ w the liquid water density (1000 kg m −3 ), and θ w the volumetric water content [-].It is assumed that a certain amount of unfrozen water is always present, even at several degrees below freezing temperature (Quinton et al., 2005).A freezing soil characteristic curve defines the unfrozen water content present for a given temperature.This curve is derived assuming that the negative soil pressure resulting from the freezing of a certain amount of liquid water is the same that would result from the drying process leading to the removal of the same amount of water (Miller, 1963;Spaans and Baker, 1996).The negative pressure is then related to a degree of depression of freezing temperature through the generalized Clausius-Clapeyron relation (Kay and Groenevelt, 1974).Most of freezing occurs at a temperature very close to and slightly lower than 273.15K.As a consequence, the freezing soil characteristic curve is very steep and highly non-linear at these temperatures.Further details on the implementation of this theory can be found in Dall 'Amico et al. (2010).
The thermal capacity and conductivity are a combination of the corresponding properties of each component of the soil multiphase mixture (soil solid, ice, liquid water, and air in the free spaces of the pores).While a simple additive mixing law is exact for the thermal capacity, the behavior of a multiphase mixture concerning the thermal conductivity is much more complex.Several non-linear mixing laws have been proposed (De Vries, 1963;Johansen, 1975;Balland and Arp, 2005), and they all behave differently (Zhang et al., 2008).While most of the formulations use different parameters according to the soil type, the formulation of Cosenza et al. (2003) is attractive for its simplicity and generality.On the analogy with the behavior of dielectric permittivity in electrical lattices, they proposed the following quadratic parallel mixing law where i corresponds to the different phases (soil solid, liquid water, ice, air), x i and k i are respectively the volumetric content and the thermal conductivity of the phase i .Equation ( 1) is discretized at the finite differences in time and in the vertical coordinate, introducing the concept of layers.A system of non-linear equations with temperature at the next time step as unknown is obtained, while temperature at the previous time step is know.A solution is sought with the Newton method (Kelley, 2003).

Surface heat flux
The surface energy balance is the forcing term of the heat equation.An accurate description of the surface fluxes is essential for modelling energy-based processes.
Incoming shortwave radiation is normally available as a measurement, and the albedo of the surface is given as a parameter (0.20 for open tundra).Incoming longwave radiation in clear sky conditions is calculated with empirical formulations (e.g.Brutsaert, 1975), which in general apply the Stefan-Boltzmann law using the air temperature measured at the surface with an effective atmosphere emissivity.In the presence of a cloudy sky, the emissivity of the atmosphere is significantly increased.The more practical way to infer cloudiness is from the ratio of the measured to clear-sky incoming shortwave radiation (Crawford and Duchon, 1999).At night, an interpolation between the values at two consecutive days is used.
The turbulent fluxes of sensible (H) and latent heat (LE ) are calculated with the flux-gradient relationship: where ρ is the air density, c p the specific heat at constant pressure, u the wind speed, T s the temperature of the surface, T a the temperature of the air, L e the specific heat of vaporization, Q * s the saturated specific humidity at the surface, Q a the specific humidity of the air, r a the aerodynamic resistance, and α and β are coefficients that take into account the soil resistance to evaporation.These coefficients are calculated according to the parameterisation of Ye and Pielke (1993), which considers evaporation as the sum of the proper evaporation from the surface and water vapour diffusion from the soil pores.The aerodynamical resistance is obtained using the theory of Monin-Obukhov.A topographically dependent wind field is obtained using topographic parameters (Liston and Elder, 2006).
GEOtop also considers the energy balance of vegetation, and the portioning of the surface energy balance into soil surface and vegetation components, as is described in Endrizzi and Marsh (2010).However, in this work, given the purpose stated above, the effect of vegetation is not considered.

Subsurface water flow equation
The Richards equation is solved in a three-dimensional form.This permits a full coupling of the normal and lateral flows.The equation is solved using the hydraulic head (H) as the unknown variable, given by the sum of water pressure head and height over a reference level.This allows writing the Richards equation as where K [m s −1 ] is the hydraulic conductivity, and S [s −1 ] a source term.Equation ( 5) is non-linear since both θ w and K depend on the water pressure, and, in turn, on the hydraulic head.The Van Genuchten model (Van Genuchten, 1980) is used to relate the volumetric water content with the pressure head when the soil is unsaturated, while for saturated soils water content is linked to the pressure through the specific storativity.The Mualem model (Mualem, 1976) provides a relation between pressure head and hydraulic conductivity in unsaturated soil, while a constant value is used for saturated soils.In partially frozen soil the hydraulic conductivity is further decreases with increasing ice content (Hanson et al., 2004).The Richards equation is discretized in time and space at the finite differences, using the same grid used for the heat equation, and solved with the Newton method.

Modelling assumptions of soil properties
Since the earth hummocks have horizontal diameters (0.4 to 1 m) smaller than the smallest grid dimension at which GEOtop can be run at Siksik creek in a reasonable computational time (10 m), they must be represented as sub-grid variability.We addressed this problem considering in the model the peat thickness uniform in the grids and equal to the values observed in the inter-hummock area, but an effective value for the hydraulic conductivity that takes into account the reduction of the overall conductivity due to the presence of the hummocks has been used.This is also motivated because the attention is turned to the spatial variability of the thawed soil depth in the portion of the soil where hillslope runoff occurs, which is the inter-hummock area.Quinton (1997) measured peat thicknesses in the inter-hummock zone along 2 transects running perpendicular to the stream channel on its west side.Despite their high degree of spread, the measurements reveal a marked decrease of the peat thickness in a 50 m distance from the main channel from around 50 cm to 25-28 cm, while, at farther distances no particular trend can be recognized (Fig. 2a).A quadratic regression curve (parabola branch) has been derived in the region where the most variation occurs, and a constant value, equal to the value at the semi-parabola vertex, has been assumed representative of the peat thickness at farther distances from the main stream.This information has been used to estimate the peat thickness distribution for the whole Siksik creek.The peat layer thickness has been considered to be dependent only on the linear distance from the main stream (Fig. 2b).
The Richards equation is solved assuming the validity of the Van Genuchten model, which describes the relation between water pressure and water content (also referred to as retention curve).However, a large uncertainty exists on the value of the Van Genuchten parameters (i.e. the α and n coefficients, the saturated and residual water contents) that are most appropriate for the peat.An estimation of α and n has been performed by Carey et al. (2007) for Granger Creek (60 • N, 135 • W), Yukon Territory, Canada, and by Quinton et al. (2005) for Scotty Creek (61 • N, 121 • W), Northwest Territories, Canada, which, in general, present similar, but thicker organic soil layers.In both cases a small decrease of water content close to saturation conditions corresponds to a large drop in water pressure.Larger pressure drops observed at the surface, showing little effect of capillary forces, while relatively smaller drops have been found at higher depths as a result of the increasing degree of decomposition, and, in turn, of the decrease in the size of particles and interparticle pores.This corresponds to high values of α and n (α ∼ 40 m −1 , n ∼ 2) at the surface, and then decreasing with depth (α ∼ 8 m −1 , n ∼ 1.35 at 1 m depth at Scotty Creek).
The concept of porosity in the peat is more complex than in other types of soil since peat contains an "active" porosity (Hoag and Price, 1997) with well-interconnected pores, and an "inactive" porosity composed of closed and dead-end pores formed by the remains of plant cells.Quinton and Gray (2003) measured the total and actual porosity at Siksik, and found that, while the total porosity is roughly constant (0.90-0.95) with depth, the active porosity is very close to the total porosity at the surface, but significantly decreases with depth.This has been represented in the model setting the residual water content equal to the inactive porosity, and the saturated water content equal to the total porosity.This implicitly assumes that the inactive porosity is always filled with liquid water that cannot drain, evaporate or freeze.
The Richards equation is solved using the Mualem model for the hydraulic conductivity, which requires the estimation of the saturated hydraulic conductivity.Quinton et al. (2008) compared data of saturated hydraulic conductivities in the direction of the hillslope drainage through peat at Scotty Creek, Granger Creek and Siksik Creek.Despite their different climates and organic soil thicknesses, they found similar relationships between depth and saturated hydraulic conductivity.In particular, values around 4 mm s −1 are typical at the ground surface, while values around 0.016 mm s −1 are common at lower peat levels, with the transition occurring mostly between 0.1 and 0.2 m depths.2008) proposed an analytical expression fitting their results.The vertical hydraulic conductivity of the peat may be different from the lateral hydraulic conductivity, since there is evidence of anisotropy in the literature (e.g.Kruse et al., 2008).In this study peat hydraulic conductivity is assumed isotropic (i) for simplicity, because as far as we know, no specific analysis has been performed for the peat at the present study site, and (ii) because vertical water flow in the peat occurs so rapidly and in very short distances, being the water table commonly very shallow, that its sensitivity to the value of the vertical conductivity is expected to be very low.The grid-averaged hydraulic conductivity is described as a weighted average of the hydraulic conductivity in the hummock and inter-hummock areas.Since the hydraulic conductivity in the hummocks is several orders of magnitudes than in the inter-hummock zone, it is reasonable to exclude the former from the calculation, and approximate the weighted average with the product of the peat hydraulic conductivity with the area fraction of the inter-hummock zone (50%).
The simulation has been carried out for the period of time spanning from 24 May to 13 September 1993 assuming an initial condition given by soil frozen and saturated, with an uniform soil temperature of −6 • C and saturation front at the surface.The soil has been discretized with 30 layers with variable thicknesses covering 3.5 m of soil depth.Soil layers are thinner in proximity of the surface (around 2 cm), become slightly thicker (up to 5 cm) until a depth of 80 cm, and significantly thicker (10-100 cm) below (a fine discretization is actually not needed at these depths, because the soil will be permanently frozen).The bottom boundary condition has been provided assigning a constant temperature of −6 • C at a depth of 8 m, which has been observed to be representative of the annual zero-amplitude temperature in proximity of the Trail Valley Creek area (Kokelj et al., 2009).

Results and discussion
The modelled end-of-summer water table depth (Fig. 3a) ranges between 10-20 cm, and always lies within the peat layer.Greater water table depths occur in well-drained areas including those with higher slope and diverging curvature.The water table is relatively close to the surface where the topography is flat or in an area of flow convergence.In particular, it is at or slightly below the ground surface in the middle of the near-stream area.
The modelled end-of-summer thaw depth (Fig. 3b) mirrors the spatial pattern of the water table depth in that deeper thaw depths are found to correspond with shallower water tables, and shallower thaw depths always occur where the water table is deep.
Markedly deeper thaw depths (up to 80 cm) occur in the internal part of the near-stream area, where the water table is very close to the ground surface.This result suggests that the spatial variability of the depth of thaw is mostly controlled by subsurface water flow, which is, in turn, largely determined by topographic variations.Other topographic factors that could affect the spatial variability of the depth of thaw, such as (i) slope and aspect, affecting shortwave radiation, and, in turn, the heat flux exchanged with the atmosphere, (ii) the spatial variability of the highly insulating organic layer thickness were found to be insignificant.It is therefore predominantly through its control on subsurface water flow that the surface topography appears to control the spatial distribution of the soil thaw.Similar patterns are encountered in the spatial distribution of the depth of thaw in other moments during the summer, with increasing thaw as summer progresses (Fig. 4).
Subsurface water flow can affect soil thaw process for the following reasons: 1. Subsurface water flow controls the liquid water content in the soil.This, in turn, affects the bulk soil thermal conductivity.The thermal conductivity of the thawed portion of the active layer controls the downwards propagation of energy to the thawing front.Given the high porosity of the peat (0.95), the low thermal conductivity of the peat solids (0.21 W m −1 K −1 ), and that the thermal conductivity of ) is significantly higher than that of the soil solid fraction, the water content becomes critical in raising the overall thermal conductivity, and in augmenting the degree of thaw (Kane et al., 2001).A spatially-variable heat conduction therefore may explain the spatial variability of the depth of thaw.However, higher liquid water content also results in a higher thermal capacity, which decreases the degree of thaw.
2. Subsurface flow in partially frozen soil also results in an advective heat transfer to the downslope areas, which may provide some energy contribution.GEOtop represents the advection of latent heat that occurs when a water particle moving into a mixture of water and ice freezes and releases heat, but it does not describe the advection of sensible heat, resulting from the temperature gradient between the moving liquid water particles and the mixture of water and ice.However, this thermal gradient is likely very small since the temperature of the water flowing over the frost table is very close to 0 • C.

Role of ground heat flux
Several studies (e.g.Abbey et al., 1978) have shown a strong correlation between the depth of thaw and cumulative ground heat flux.At Siksik Creek, the ground heat flux averaged over the simulation time (Fig. 5a) mirrors very well the end-of-summer depth of thaw, and show values ranging from 14 W m −2 where the thaw is shallower to 22 W m −2 where the thaw is deeper.The average net shortwave radiation (Fig. 5b) has instead higher spatial variability (130 W m −2 on north-facing slopes, 155 W m −2 on south-facing slopes, and 145-150 W m −2 on flat surfaces such as valley bottoms and plateaus), but shows no significant correlation with the soil thaw depth and ground heat flux.Therefore, the spatial variability of the ground heat flux is here mostly induced by soil moisture and surface temperature spatial patterns, rather than shortwave radiation spatial patterns (slope and aspect).The water content in the 3 cm soil layer closest to the ground surface, averaged over the simulation time (Fig. 6a), scales very well with the end-of-summer water table depth.The surface temperature averaged over the simulation period (Fig. 6b) shows instead a clear negative correlation with the soil water content at the surface: lower temperatures (∼ 11 • C) are associated to wetter areas, and higher temperatures (∼ 12.5 • C) to drier areas.In the presence of uniform shortwave radiation input, wetter soils lose more energy through evaporation and, therefore, the surface temperature is lower than in drier soils.In addition, drier soils have lower thermal conductivity than wetter soils, which contributes to further increase the surface temperature, although this effect may be offset by their lower thermal capacity.In this application, shortwave radiation is not uniform.Therefore, the modelled surface temperature reflects both the spatial patterns of soil moisture, and those related to the spatial distribution of solar radiation.The former seem much more evident, demonstrating that in this application the strongest control on the surface temperature is exerted by soil moisture.This can be better understood if we look individually at the single components of the surface heat flux averaged over the simulation time.
The maps of sensible and latent heat fluxes (Fig. 7a  The net longwave radiation (Fig. 7d) is controlled by surface temperature, since the outgoing flux scales with the fourth power of the absolute temperature.The timeaveraged net longwave heat flux is always an energy loss, and varies from −78 W m −2 , where the surface temperature is higher and, therefore, the loss is larger, to −68 W m −2 , where the surface temperature is lower.This result has 2 important implications: (i) while the sum of the turbulent fluxes shows higher energy losses in the wetter areas, this effect is offset if we add the longwave radiation and consider the total surface heat flux; (ii) the longwave radiation losses are higher on south-facing than north-facing slopes, following the spatial pattern of surface temperature, and, therefore, contribute to further buffer the effect of spatial patterns of shortwave radiation in the spatial distribution of the ground heat flux.The net result is that wetter areas have relatively low energy losses to the atmosphere, which contributes to increase the rate of ground thaw.Subsurface lateral water flow, therefore, seems to completely control the spatial patterns of the ground heat flux, which, in turn, affects the spatial variability of ground thaw.

Relative importance of factors affecting thaw depth
The relative importance of the factors that were shown to affect the thaw depth spatial distribution was assessed by systematically switching-off each factor individually, and comparing the resulting map showing the spatial pattern of ground thaw with the reference simulation map (Fig. 3b), which accounts for all factors.
When subsurface lateral water flow was turned off by setting the hydraulic conductivity at zero, the spatial pattern of the end-of-summer depth of thaw (Fig. 8a) is dominated by the spatial variability of the peat layer thickness which ranges from 27 cm (50 cm end-of-summer ground thaw) and 60 cm (45 cm end-of-summer thaw).Thus, outside of the near-stream area the spatial variation in thaw depth is low (50-55 cm), however north-facing slopes have slightly greater thaw owing to their relatively high near-surface moisture content and therefore higher thermal conductivity.Evaporation is the only means of water loss with the hydraulic conductivity set to zero, and since evaporation is higher on south facing slopes, their surface are relatively dry and are therefore better able to insulate the underlying frost table.
The advection of latent heat was set to zero by excluding both vertical and lateral subsurface water flow in the partially or completely frozen soil layers.The resulting end-of-summer ground thaw depth map (Fig. 8b) shows only minimal differences to the correspondent map resulting from the reference simulation (Fig. 3b).Therefore, advection of latent heat does not appear to have a significant influence on the spatial pattern of thaw.
The spatial variability of the ground heat flux due to the spatial distribution of soil moisture and temperature at the surface is removed setting the sensible and latent heat fluxes at zero, and deriving the outgoing longwave radiation from the basin-averaged surface temperature rather than the local surface temperature.The resulting end-ofsummer frost table depth (Fig. 8c) shows marked spatial variability (0-90 cm), which appears entirely controlled by the spatial patterns of shortwave radiation and completely independent of the spatial distribution of the soil moisture.
The spatial variability of thermal conductivity was removed by considering the soil thermal conductivity independent of the water and ice contents, and by setting it to a constant value equal to 0.4 W m −1 K −1 .The resulting end-of-summer depth of soil thaw (Fig. 8d) has moderate spatial variability, ranging from 60 to 70 cm, and shows some correspondence to the spatial patterns of solar radiation.However, unlike the reference simulation, Fig. 8b indicates shallower thaw in wet regions.If the surface energy balance is considered, wet and dry areas present larger differences in the net turbulent heat fluxes than in the base simulation.As shown in the previous section, net turbulent fluxes produce a greater energy loss in wet areas.While in the map of ground heat flux resulting from the reference simulation the effect of longwave radiation, which instead generates higher energy losses for the drier areas, prevails, here the energy loss patterns due to net turbulent fluxes are prevalent.This is explained by the smaller differences in the surface temperature between wet and dry areas resulting from a thermal conductivity independent of soil moisture.Wet areas have higher surface temperature, and are subjected to a stronger evaporation than in the base simulation.This increases energy losses in wet areas, and reduces the degree of soil thaw compared to dry areas.
The comparisons with the reference simulation above suggest that the spatial distributions of the ground heat flux and soil hydraulic conductivity are the most important factors controlling the spatial pattern of soil thaw.Since wet surfaces are more thermally conductive, they remain relatively cool, have less evaporative loss of energy, and therefore have greater thaw depths.This process dominates the surface energy balance, and explains why the spatial patterns of the depth of thaw obtained in the model results are largely independent of the spatial distribution of shortwave radiation.

Comparison with field data
The thaw depths measured daily at the three study plots (Fig. 1) were averaged on a weekly basis.The frequency distribution of thaw depth for each week of observations is considered representative of the ∼ 40 m wide intervening stream-bank section (Fig. 1).This procedure allows for comparison of a modeled and observed frequency distributions of thaw depth representing different stages of active layer thaw (Fig. 9).
The measured and modelled modal thaw depths compare reasonably well for each time period, although the simulated thaw progressed at a slightly higher rate so that by the last time period, the modal depth of the simulated thaw was approximately 5 cm greater by the last time period.However, the model overestimates the minima and maxima of thaw in most time periods.In particular, the development of a tail in the above-average side is evident in the frequency distribution of the simulated results, with pronounced maxima slightly below 80 cm, corresponding the wettest portions of near-stream area.However, the range of variability of the observations and simulated results is comparable.It is particularly different only at the beginning of the summer, when the observations, in opposition to the modelled results, are probably significantly affected by the spatial variability the snow cover, and at the end of the summer, as a result of the tail towards the maxima values observed in the simulated results.
Delayed snow melting may be responsible for the overestimation of the minima and maxima of the depth of the spatial distribution.In addition, GEOtop may improperly represent the effect of scale on subsurface flow processes, which may bias the simulation of thaw.Subsurface flow may be controlled by soil and topographic heterogeneities occurring at the sub-grid (i.e.<10 m) scale.As a consequence, the saturated area may extend over a smaller surface, so that the tail towards the maximum values in the frequency distribution of the depth of thaw is less accentuated.Drier areas may also be on average even drier, since water probably accumulates in micro-topographic depressions encompassing a small fraction of the surface.
Transpiration may also reduce the depth of thaw.Since biomass is normally larger in wet areas, this may also attenuate the difference in the degree of thaw between wet and dry areas.Larger biomass corresponds to higher transpiration and therefore greater energy loss.Local areas of thick peat accumulation may further reduce the highest thaw depth values observed in the model results.Other uncertainties can be associated with the estimation of the components of the surface heat flux, i.e. incoming longwave radiation, turbulent fluxes.

Conclusions
The factors controlling the spatial variability of the thaw in a typical arctic tundra basin was investigated with a fine resolution model that couples energy and water flow processes.This work has improved the understanding of how the pattern of thaw is affected by the ground surface, the surface heat flux exchanged with the atmosphere, and the subsurface flow regime.Although the snow cover is know to be an important mediator of energy exchange between the soil surface and the atmosphere, it was not considered in the present study.The thaw simulations indicate that the spatial patterns of ground thaw are not dominated by slope and aspect, but are instead entirely controlled by the spatial distribution of soil moisture, which is determined by subsurface flow patterns.Wet surface being more thermally conductive than dry areas, transmit heat with greater efficiency, while remaining relatively cool.On the other hand, dry surfaces are relatively warm because they are relatively poor thermal conductors.This leads to a spatial distribution of sensible heat flux and outgoing longwave radiation that has lower energy loss in wet areas than in dry areas.These patterns are not offset by the higher energy losses due to evaporation, which are larger in wet areas.Measured thaw depths have a similar range of variability to the simulated values for each stage of active layer development, although the model slightly overestimated the depth of thaw.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

369
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

373
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

375
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

377
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

379
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Quinton et al. ( Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

381
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | liquid water (0.57W m −1 K −1 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | and b) show negative values (i.e.directed towards the atmosphere) and reflect the spatial patterns of the time averaged surface temperature and soil moisture, respectively.The fluxes show a marked spatial variability: the sensible heat flux values range from −15 W m −2 in the colder points to −40 W m −2 in warmer areas, while the latent heat flux values range from −20 W m −2 in the drier areas to −50 W m −2 in the wetter areas.If we add up the turbulent fluxes of sensible and latent heat (Fig. 7c), the sum, which we refer to as net turbulent heat flux, shows a significant negative correlation with solar radiation, varying from −48 W m −2 where shortwave radiation is lower to −63 W m −2 where shortwave radiation is higher, and a moderate variability due to surface soil moisture patterns, with the wetter areas losing more energy (approximately −62 W m −2 for the wet areas and −57 W m −2 for the dry areas).This range of variability has the effect to compensate the effect of shortwave radiation patterns in the spatial distribution of the ground heat flux.383 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

385
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

387
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 5 .
Fig. 5. Ground heat flux [W m −2 ], sum of net all-wave radiation and turbulent fluxes of sensible and latent heat (a), and net shortwave radiation [W m −2 ], incoming minus outgoing (b), both averaged over the simulation time.
(Marsh et al., 2010)imately 55 km northeast of Inuvik (Northwest Territories, Canada).Siksik Creek is biophysically representative of small catchments of the continuous permafrost zone near the northern fringe of the forest-tundra transition zone of northwestern Canada.It is mostly open tundra, with shrub patches more frequent in proximity to the stream channel(Marsh et al., 2010).Its drainage area extends over 1005 km 2 , and its elevation ranges from 47 to 102 m a.s.l.Even if it is located in an area