Interactions between soil thermal and hydrological dynamics in the response of Alaska ecosystems to fire disturbance

] Soil temperature and moisture are important factors that control many ecosystem processes. However, interactions between soil thermal and hydrological processes are not adequately understood in cold regions, where the frozen soil, fire disturbance, and soil drainage play important roles in controlling interactions among these processes. These interactions were investigated with a new ecosystem model framework, the dynamic organic soil version of the Terrestrial Ecosystem Model, that incorporates an efficient and stable numerical scheme for simulating soil thermal and hydrological dynamics within soil profiles that contain a live moss horizon, fibrous and amorphous organic horizons, and mineral soil horizons. The performance of the model was evaluated for a tundra bum site that had both preburn and postburn measurements, two black spruce fire chronosequences (representing space-for-time substitutions in well and intermediately drained conditions), and a poorly drained black spruce site. Although space-for-time substitutions present challenges in model-data comparison, the model demonstrates substantial ability in simulating the dynamics of evapotranspiration, soil temperature, active layer depth, soil moisture, and water table depth in response to both climate variability and fire disturbance. Several differences between model simulations and field measurements identified key challenges for evaluating/improving model performance that include ( 1 ) proper representation of discrepancies between air temperature and ground surface temperature; (2) minimization of precipitation biases in the driving data sets; (3) improvement of the measurement accuracy of soil moisture in surface organic horizons; and (4) proper specification of organic horizon depth/properties, and soil thermal conductivity.


Introduction
[2] Soil temperature is considered one of the most important environmental factors affecting soil organic mat- [3] There are two primary ways in which soil moisture can influence soil thermal dynamics: (1) the thermal conductivity of a dry organic soil (i.e., organic soil horizon composed of greater than 18% organic C) is substantially lower than that of a wet organic soil, which makes dry organic soil a good heat insulator [Yi et al., 2007]; (2) the seasonal amplitude of soil temperature is damped through the release and absorption of latent heat. Conversely, the thermal state of soil can also influence hydrology: (1) frozen soils have limited infiltration capacity, which results in a large runoff during spring snowmelt [Shanley and Chalmers, 1999]; and (2) the base flow depends on the extent of unfrozen soil in the hydrologically active zone. For example, deeper unfrozen soil layers are thought to contribute to an increase in winter discharge from rivers flowing into the Arctic Ocean [Oelke et al., 2004). The interactions between soil thermal and hydrological dynamics are generally implemented in third generation land surface models, which are used in general circulation models to simulate the lower boundary water, heat and momentum fluxes [Verseghy, 1991;Bonan, 1996;Oleson et al., 2004]. However, these interactions are generally neglected in most large-scale ecosystem models where one or two soil layers are used to simulate soil moisture dynamics in ecosystem models [e.g., Sitch et al., 2003], and analytical or empirical functions are used for simulating soil thermal dynamics [e.g., Bond-Lamberty et al., 2005). Large-scale ecosystem models have considered vertical soil thermal dynamics (e.g., the Terrestrial Ecosystem Model, TEM [Zhuang et al., 2001[Zhuang et al., ,2002[Zhuang et al., ,2003Euskirchen et aI., 2006]), or the effects of soil moisture on biogeochemical dynamics [e.g., Zhuang et al., 2004], but have not fully included a coupling between soil thermal and vertical soil moisture regimes.
[4] The active layer is the top portion of the soil that thaws during summer and freezes again during autumn in permafrost regions. Most physical, chemical and biological processes happen in the active layer. A change in active layer depth (ALD) may have substantial implications for ecosystem carbon balance [Goulden et al., 1998). The position of water table is an important indicator of soil wetness, and water table depth (WTD) is an important control on soil carbon decomposition [Dunn et al., 2007]. Wildfire disturbances and drainage are important factors that influence soil thermal and hydrological regimes to affect ALD and WTD [Harden et al., 2000[Harden et al., ,2006. Wildfire not only reduces the thickness of organic horizons, but also alters the surface properties [Liu and Randerson, 2008). Soil drainage not only affects the hydrological dynamics of soil directly, but also influences the soil thermal dynamics indirectly through effects on the thickness of organic horizons that result from the responses of ecological processes and fire disturbance to drainage [Harden et al., 2006]. For example, poorly drained ecosystems usually have thicker organic horizons than well drained ecosystems because they experience less frequent and less severe fires [Harden et al., 2000]. However, the effects of drainage and wildfire are seldom considered in the simulation of soil temperature, active layer, soil moisture, and water table dynamics by ecosystem models [see Zhang et al., 2002;Bond-Lamberty et al., 2007;Ju and Chen, 2008).
[s] In this study, we use a process-based modeling approach to investigate interactions between soil thermal and hydrological dynamics in the response of Alaska ecosystems to fire disturbance. Such an approach requires a modeling framework that incorporates the interactions among organic horizon thickness and properties, soil temperature, and soil moisture into an ecosystem model. It is the main goal of this paper to build on the progress of previous versions of the Terrestrial Ecosystem Model (TEM) as represented by Zhuang et at. [2002,2004] to (I) incorporate different types of organic horizons, i.c., live moss, and fibrous and amorphous organic horizons, into the soil structure; (2) differentiate the effects of drainage by incorporating two broad drainage classes, i.e., moderately well drained and poorly drained; and (3) develop and evaluate an efficient and stable numeric scheme for simulating soil thermal and hydrological dynamics within heterogeneous soil. We used the new version of TEM to investigate the following questions: (1) What is the effect of organic horizon thickness on soil thermal and hydrological dynamics? and (2) How does soil drainage intluence interactions between soil thermal and hydrological dynamics? To address these questions we conducted a sensitivity analysis of how modeled active layer depth and water table depths responded to changes in important parameters and atmospheric driving data.

Background and Overview
[6] The Terrestrial Ecosystem Model (TEM) is a processbased ecosystem model designed to simulate the carbon and nitrogen pools of vegetation and soil, and carbon and nitrogen fluxes among vegetation, soil, and atmosphere [Raich et al., 1991). While previous model development efforts have improved the soil thermal and hydrological processes in TEM for application in high-latitude regions [Zhuang et al., 2001[Zhuang et al., , 2002[Zhuang et al., , 2003[Zhuang et al., , 2004Euskirchen et {II., 2006], soil thermal and hydrological processes are not comprehensively coup led, and fire disturbance reduced the amount of soil carbon without affecting organic soil thickness and associated changes in the thermal and hydrological properties of organic soil [e.g., see Balshi et al., 2007J. Zhuang et al. [2002 conducted model experiments that demonstrated that changes in organic matter horizons during and after fire potentially have important influences 011 soil temperature and moisture, hut subsequent modeling efforts have not dealt with the issue of dynamic changes in organic horizons. Therefore, the model development research reported here is focused on the explicit coupling of soil thermal and hydrological processes in the context of a changing organic horizon, which is a necessary step toward dynamically simulating how changes in organic matter horizons during and after fire influence the interactions among soil thermal, hydrologic, and biogeochemical processes.
[7] In this study, we first describe the ncw environmental module (hereafter EnvM) that is responsible for simulating soil thermal and moisture dynamics in the dynamic organic soil framework of TEM (DOS-TEM). We then describe the field sites that were used to evaluate the EnvM. We also describe the information required for the operation of the 20f20 G()2()15 YI ET AL.: SOIL THERMAL AND HYDROLOGICAL DYNAMICS G02015 EnyM including site-specific parameters and atmospheric driving data. We then describe the validation data used to evaluate the model and the sensitivity analyses we conducted to better understand interactions between soil thermal and hydrological dynamics that are influenced by fire disturbance and soil drainage.

Model Development
[8] The detailed descriptions of the water and energy fluxes among the atmosphere, canopy, snow and soil, and within the soil are provided in Appendices A-E. In the EnvM, the ground surface is represented by a snow horizon, three soil organic horizons, two mineral soil horizons, and a rock horizon. Each horizon is further divided into layers that are explicitly treated with respect to energy and moisture exchange. For example, the snow horizon can consist of up to five snow layers.
[9] The organic soil horizons include live moss, fibrous organic soil, and amorphous organic soil. For accurate simulation of soil temperature and moisture, the soil horizons near the surface are divided into thin layers (e.g., layers are a few centimeters thick in the live moss horizon), and layers become thicker as the distance from the surface gets deeper (e.g., layers are approximately 10 m thick in the rock horizon). The number of layers is variable among the three organic soil horizons, which can consist up to 7 layers.
[10] Following the method used in land surface models, e.g., the Canadian Land Surface Scheme [Verseghy, 1991 J, the EnvM considers upper and lower mineral soil horizons. Each mineral soil horizon is 1 of 11 mineral soil types as defined by Beringer et al. [2001]. The upper mineral soil horizon consists of 2 layers with thicknesses of 10 and 20 ern, The lower mineral soil horizon consists of 3 layers with thickness of 50, 100, and 200 cm. Thus, the total thickness of the mineral soil is 3.8 m in EnvM. The rock horizon consists of 5 layers. The total thickness of the soil and rock horizons is around 50 m.
[11] During model simulations, a freeze-thaw front, using Stefan's algorithm, is introduced into each layer to overcome computational problems associated with heterogeneous thicknesses of the layers. Each layer can contain an unlimited number of fronts. The temperature of each layer is updated daily, and the phase changes can only happen in layers of the snow and soil horizons. The water content can only be updated in layers for which the temperature is above O°c. Each model run spans 1901Each model run spans -2006 with the years before 1999 used for model spin-up.

Description of Sites Used for Model Evaluation
[12] Data from seven sites in Alaska were used to evaluate the EnvM in this study (Table 1): a tussock (Eriophorum vaginatum ) tundra site located at Kougarok (K2) on the Seward Peninsula; a poorly drained black spruce (Picea mariana) site (FBKS) located near Fairbanks, Alaska; and two black spruce fire chronosequences located in Donnelly Flats near Delta Junction, Alaska (DFTC, DFT87, and DFT99; DFCC and DFC99), which represent well drained (DFTC, DFT87, and DFT99 sites) and intermediately drained (DFCC and DFC99 sites) conditions. The K2 site, whieh is located in an area of transition between continuous and discontinuous permafrost, experienced a severe to moderate bum in 2002, which removed 7 -9 cm of organic soil between the Eriophorum vaginatum tussocks [Liljedahl et al., 2007]. The K2 site is unique in that it has soil temperature and moisture measurements obtained from the same location both before and after the fire [Liljedahl et al., 2007].
[13] The FBKS site, which is located on the campus of the University of Alaska Fairbanks, is a poorly drained black spruce forest. The dominant overstory vegetation is Picea mariana with a forest floor that is a mixture of tussocks, vascular plants, shrubs, Spaghnum moss, feather mosses, and lichen. These plants include Betula glandulosa, Ledum palustre, Vaccinium vitis-ideae, Carex lugens, Sphagnum spp., Hylocomium splendens, Thuidium abietinum, and Cladina stellaris [Kim et al., 2007]. There is about 50 cm of organic soil overlying loess, with 8 cm of feather and Sphagnum spp. mosses on top of the organic soil.
[14] The DFTC, DFT87, and DFT99 sites are part of a well drained fire chronosequence located in Donnelly Flats, which is near Delta Junction, Alaska; with the most recent stand replacing fires in~1920, 1987, and 1999, respectively [Liu and Randerson, 2008. The tower control site in Donnelly Flats, DFTC, has an overstory of mature 80-year old black spruce trees and an understory dominated by feathermoss (Hylocomium splendens, Pleurozium schreberi, and Aulacomnium spp.) and lichen (genera: Cetraria, Cladonia, Cladina, and Peltigera [15] The DFCC and DFC99 sites are part of the intermediately drained fire chronosequence located in Donnelly Flats along Donnelly Creek. We treat the Donnelly Flats Creek Control (DFCC) site, which last burned in 1885, as the control site for DFC99, which burned in 1999. The ground cover is dominated by feather moss tHylocomium splendens, Pleurozium schreberi, and Aulacomnium spp.) at DFCC and by the early successional postfire mosses Ceratodon purpureus and Polytrichum spp. at DFC99. These two sites were underlain by permafrost within 1 m of the soil surface at the time of the fire in 1999.

Site-Specific Input Parameters
[16] Site-specific parameters, including monthly projected leaf area index (LAI, Table 2) and organic soil thicknesses (Table 3), were used by EnvM in this study. For LAI of the K2 site, the value of 0.52 was used for all seasons based on Beringer et al. [2005]. The simulated monthly LAI from a previous version TEM for mature black spruce, which ranged from 1.1 to 2.0, were used for the DFTC, DFCC, and FBKS sites. For the DFT99 and DFC99 sites, we assumed that the LAI was 0.05 after the fire in 1999. The monthly LAI values of DFT87 were estimated based on satellite data [Liu and Randerson, 2008].
[17] The organic soil thicknesses were specified from field measurements of live moss, fibrous organic soil that includes dead moss or Oi horizons, and amorphous organic soil that includes mesic Oe or amorphous Oa horizons. For the K2 site, the soil organic horizons were 14 cm thick in year 2000, half of which was burned in 2002[Liljedahl et al., 2007. At the FBKS site, the soil organic horizons of 50 cm thick. The mean thicknesses of organic horizons of the DFCC and DFTC sites arc 20 and 10 cm, respectively. After the fires in 1999, only 10 and 3.5 cm organic soil remained at the DFC99 and DFT99 sites, respectively [Manies et af., 2001[Manies et af., , 2004. The values in Table 3 for the Donnelly Flats sites represent the organic horizon thicknesses of the soil cores collected at the locations of soil moisture and temperature measurements.
[18] We used the measured fraction of fine root production at a black spruce site of the Bonanza Creek Long-term Ecological Research Program [Ruess et al., 2006] to represent the root distribution used for estimating transpiration. Since we could not identify measurements of root distribution for aspen and tussock tundra in Alaska, the rooting distributions of balsam poplar and black spruce from Ruess et al. [2006] were used to specify root distributions for the DFT87 and K2 sites, respectively (   [22] The N factor (/1) is used to estimate temperature at the ground surface, either the soil surface or the snow surface, from atmospheric surface temperature [KIene et al., 2001]. To our knowledge, snow surface temperature has never been measured for this purpose. Therefore, we defined as n to be 1.00 when snow is present. When snow is not present, we defined 11 as the ratio between thawing degree day sums of approximately 2 m air temperature and the ground surface. We found that n was 0.66, 0.94, and 1.10 for DFTC, DFT87 and DFT99, respectively, based on measured air temperature and surface soil temperatures. For the K2 site, n was assumed to be 1.00 since there were no measurements of ground surface temperature. For DFCC and DFC99, n was assumed to be 0.66 and 1.10, respectively, to be consistent with the stand age based estimates for DFTC and DFT99. Because the FBKS site is a mature black spruce site, we assumed that n as 0.66 to be consistent with the estimates for the other mature black spruce sites in the study.

Validation Data
[23] Eddy covariance measurements of water fluxes are available for DFT99, DFT87, and DFTC during the years 2002-2004. Seasonal variations and site differences in the closure of the surface energy balance were summarized by Liu and Randerson [2008]. In general, seasonal mean values of the closure ranged from 0.69 to 0.81 for the DFT99 site, from 0.75 to 0.86 for the DFT87 site, and from 0.71 to 0.87 for the DFTC site. The closure estimates are within the range of these reported by FLUXNET cornmunity [Wilson et al., 2002]. These data were summed to daily resolution if there were more than 30 valid half-hour measurements in one day, otherwise, the daily water flux was considered as missing. Similarly, daily estimates were aggregated to make monthly estimates. The simulated monthly evapotranspiration, canopy evaporation, sublimation, transpiration, snow sublimation, and soil evaporation were summed and compared with the field-based estimates of water fluxes [Liu and Randerson, 2008] to evaluate the performance of the EnvM in estimating these water fluxes.
[24] Soil temperature and moisture measurements were available for all seven sites at different depths (2-100 em) and for different periods (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006). The soil temperature and moisture estimates were output from the EnvM at the same depth as each measurement for purposes of comparison. The thaw depth was measured approximately monthly using a permafrost probe at DFTC, DFT99, DFCC, and DFC99.

Sensitivity Analysis
[25] A sensitivity analysis was performed to investigate the effects of changes in parameters and driving data on model estimates for both a moderately well drained and a poorly drained black spruce forest. The baseline thicknesses of the live moss, fibrous organic and amorphous organic 50f20 G02015 Yl ET AL.: SOlL THERMAL AND HYDROLOG1CAL DYNAM1CS G02015 horizons were set to 2, 7, and 5 cm for the moderately well drained forest, and to 3, 17, and 14 cm for the poorly drained forest based on the mean organic horizon thicknesses of black spruce stands in Manitoba, Canada [Manies et al., 2006;n et al., 2009] and Alaska [Manies et al., 2003]. For the baseline simulation, the 2006 atmospheric driving data of the Donnelly Flats sites were used to run EnvM without changing any parameters. Changes to individual driving data inputs or parameters were then applied for each simulation of the sensitivity analysis. The input data and parameters considered in this sensitivity analysis are: atmospheric warm season temperature (i.e., monthly temperature > O°C), snow, rain, surface solar radiation, organic soil thickness, LAI, topography factor, maximum drainage rate, minimum soil thermal conductivity, maximum snow density, snow albedo, and minimum soil evaporation ratio. The results of the equilibrium state after 100 years of simulation were analyzed in the sensitivity analysis.
[26] Because of the importance of active layer depth (ALD) and water table depth (WTD) on biogeochemistry of high-latitude terrestrial ecosystems, we focused the sensitivity analysis on evaluating the sensitivity of ALD and WTD to changes in the chosen input data and parameters. We evaluated the sensitivity of ALD and WTD using a sensitivity index similar to that used by Friend et al. [1993]. As an example, the sensitivity of ALD to air temperature (SAW) is defined as: driver. Here the absolute value of the relative change in air temperature is used so that the sensitivity index provides information on whether ALD increased or decreased.

Evapotranspiration
[27] The monthly field-based estimates of evapotranspiration (ET), which were available for sites of the tower chronosequence from 2002 to 2004, were substantially different among the DFTC, DFT87, and DFT99 sites ( Figure 2). The simulated ET explained 80, 89, and 83% of the variability of the field-based ET estimates for the DFTC, DFT87 and DFT99 sites, respectively. The slope of the regressions between simulated and observed was less than 1 for both the DFTC and DF87 sites because of underestimation of ET in summer ( Figure 2). The intercepts of the regressions between simulated and observed werc not different from 0 at any of the three sites. See Table S1 in thc auxiliary material for more detail.

Freezing and Thawing Fronts
[28] In the EnvM, the freeze-thaw fronts were explicitly simulated using the Stefan algorithm. At the K2 site, the simulated prefire ALD [29] The simulated maximum ALDs of DFTC and DFT99 in 2001 were around 0.65 and 1.8 rn, respectively (Figures 3c and 3d). Similarly, the simulated maximum ALD at the DFCC and DFC99 sites are about 0.6 m and 1.4 m, respectively (Figures 3e and 3f). The simulated and seasonal ALDs agree well at DFTC and DFCC in all months except September (Figures 3c and 3e) The ALD in September is underestimated at DFTC and overestimated at DFCC, however, the simulation is within 1 standard deviation of measurement, 37 and 11 cm for DFTC and D FCC, respectively, At the recently burned sites, DFT99 and DFC99, the simulated thaw depth is similar to the observed thaw depth in April and May, but is deeper than observed in June, July, and September (Figures 3d and 3f), The accurate measurement of thaw depth at these sites in July, August, and September was constrained by the existence of shallow gravel layer that prevented the permafrost probe from penetrating deeper. Soil pits excavated shortly after the fire in 1999 indicated that the permafrost table at DFT99 and DFC99 were deeper than 2 m and I m, respectively [Harden et al., 2006].
[30] Additional simulations were conducted for DFTC, DFT99, DFCC, and DFC99 with the assumption that n factor equals 1.00. The simulated ALDs are about 0.8, 1.5, 0.8, and .1.3 m for DFTC, DFT99, DFCC, and DFC99, respectively. Differences in organic soil depths between where ALD hl is the simulated baseline ALD which used Ta bl as the baseline air temperature driver and ALD a is the simulated ALD using the Ta a as the altered air temperature G02()15 YI ET At.: SOIL THERMAL AND HYDROLOGICAL DYNAMICS G02015 DFTC and DF99 and between DFCC and DF99 were responsible for the ALD differences of 0.7 and 0.5 m between each pair of sites, respectively. Comparison of the differences in ALD between the mature and burned sites with an n factor of 1 versus the n factors we used revealed that approximately 60% of the differences can be explained by organic soil thickness. Thus, fire appears to be an important factor responsible for differences in ALD between mature and burned sites.

Soil Temperatures
[31J After updating the freeze-thaw fronts, the soil temperature is updated in the EnvM. We analyzed scatterplots of simulated and measured soil temperature at (1) Figure 4d). On average, the model-based estimate of soil temperature explained more than 78% of variation in measured soil temperatures. At the K2 site, the slope of linear regressions between simulated and measured soil temperatures were not different from 1 (95% confidence interval) at cach of the three depth comparisons; the intercepts were slightly negative but within 1.5°C of 0 o C At the FBKS site, while the slopes of the regressions are approximately 0.85, they arc technically not significantly different from a slope of 1. However, the intercepts of the regressions at the FBKS site, which are approximately -3 o C, are significantly less than 0 o C. In general, the slopes 70f20 G02015 YI ET AL SOIL THERMAL AND HYDROLOGICAL DYNAMICS G02U15 of the regressions for the Donnelly Flats sites are not significantly different from I. However, the range of the intercepts at the Donnelly Flats sites is between -0.5 o C and -3.5 o C, with most of the intercepts significantly less than 0°c. The negative intercepts are associated with underestimates of soil temperature that largely occur during winter. See Table S2 in the auxiliary material for more detail.
[32] The EnvM performed well in simulating soil temperature at the K2 site before the fire event in 2002, but underestimated the summer soil temperature (Figure 5a). It is possible that fire alters the surface characteristics, and that the n factor following fire may be greater than 1. At the Donnelly Flats sites, the EnvM underestimated soil temperatures during the winter of 2003-2004. The underestimation was about 6-10°C at the DFCC site (Figure 5b). In comparison to other winters reported in Figure 5b, the winter of 2003 -2004 was characterized by substantially lower precipitation and colder atmospheric temperature in our driving data sets (Figures lc and 1f). The simulated colder soil temperatures in the winter 2003-2004 are consistent with the climate data used to drive the simulations. In contrast, the measured soil temperature implies that the snowpack was thicker than in the other years, which is inconsistent with the precipitation in our driving data sets. Therefore, we suspect that the precipitation data are biased in 2003 -2004 or that wind redistribution of snow results in similar snow depths at the sites in each year. Table Depth [33] As with the analysis of the simulated freeze-thaw fronts, we analyzed the simulated WTDs fix the K2 site both before (2001) (Figure 3b). At the Donnelly Flats sites, the simulated WTDs decrease in April because of snowmelt infiltration. During the growing season, the simulated WTDs generally increase because of evapotranspirational water losses, but there are occasional decreases in WTD because of rainfall events. In general, the simulated WTDs of the mature forests (DFTC and DFCC) are shallower than those of the burned sites (DFC99 and DFT99) as deeper ALDs increase base tlow at the burned sites.

Soil Moisture
[34] The soil moisture of each soil layer in the EnvM is updated immediately after soil temperature is updated. We analyzed scatterplots of simulated and measured soil temperature at (1) all depths (Figure 6a), (2) shallow depths (~2 cm; Figure 6b (~25 cm: Figure 6c), and (4) deep depths (>38 cm; Figure 6d). In general, the correspondence between simulated and measured soil moisture is not as tight as the comparison between simulated and measured soil temperature. At the K2 site, the slope of linear regressions between simulated and measured soil moistures was significantly less than 1, with the slope highest (0.84) for the intermediate depth comparison, which was the only comparison in which the intercept was not significantly different from ODe. At the FBKS site, the slopes of the regressions were not different than 1 and the intercepts were not significantly different from 0 o C. At Donnelly Flats sites, the slopes of the regressions were generally not significantly different from 1 and the intercepts were generally not significantly different from 0 o C. See Table S3 in the auxiliary material for more detail.
[15J The most successful simulations at all sites where those of the intermediate depths (~25 cm) (Figure 6c), which captured the seasonal variations (Figure 7). It is difficult to determine if the mismatches at shallow depths indicate problems with the model at shallow depths or are associated with biases in the measurements as it is very difficulty to accurately measure soil moisture in very porous organic soils ( Yoshikawa et al., 2004;Overduin et al., 2005]. We suspect that the mismatches in the deeper soil 90f20 layers are associated with uncertainties of representing lateral drainage in our simulations as the model assumes that lateral drainage only occurs in mineral soils at depths deeper than 80 cm. For the only very poorly drained site, FBKS, the simulated soil moistures are saturated at intermediate depth (30 cm), which are overestimated compared with observations ( Figure 7). Several additional simulations were conducted for FBKS in which several parameters and driving variables were varied, but the results were not improved. The assumption of no lateral drainage in the organic soil horizons appears to not be valid for this site, which usually has ALD occurring in organic soil horizons. For moderately well drained sites, at intermediate depths, there are two notable types of differences between simulated and measured soil moisture. One is an artifact associated with the thermal timing of phase change in which the simulated increase in soil moisture lags the measured increase in soil moisture because the model assumes 5% volumetric unfrozen water content when soil is frozen, e.g., in the springs 2001 and 2002 at the DFTC site (Figure 7). The other difference is associated with improper representation of soil texture in the model, e.g., the porosity of the mineral soil at 25 cm of the DFTC site is specified to be less than 50%, but the measured volumetric water content for this site exceeded 60% in summer 2002. [37] Warm season air temperature (Ta), organic soil thickness (OS), and minimum soil thermal conductivity (TcS1Min) directly affect the thermal regime of soil ( Figure 9). An increase of Ta or TcS1Min or a decrease of OS caused an increase of ALD and vice versa (Figures 8a  and 8b). Increases of ALD cause increases of WTD for moderately well drained black spruce forest, duc to the increases of base flow/drainage. But the change of WTO for poorly drained black spruce forest is small, since the base flow is set to O. These differences are demonstrated in Figures 8a and 8b; the symbols in Figure 8b are more clustered near the vertical axis than those in Figure gao [38] The role of subsurface drainage is also seen in contrasting responses-of WTD to decreases in rain in which WTD of the poorly drained black spruce forest increases (see Figure 8b, region 11), but WTD of the moderately well drained forest increases (see Figure 8a, region ][1). The decreased soil moisture associated with rain initially decreases the thermal conductivity of the soil, which causes a decrease in ALD (see both Figures 8a and Sb). In the moderately well drained forest the decrease in ALD results in less subsurface drainage (Figure 9), which decreases WTD. Because drainage is not affected by ALO in the poorly drained forest, the decrease in rain increases WTO.
[39] The estimates of ALD and WTD simulated by the model were not responsive to changes in the other driving variables and parameters considered in the sensitivity analysis. One would expect that soil moisture would be affected by the surface evaporation ratio parameter (EvS1Min) and the topography parameter (Fsat). This insensitivity is caused by short-term self-adjustment of model, i.e., if infiltration is reduced by a decrease in EvSIMin, then the surface soil will become drier and there will be less runoff and more infiltration with the next rain event (Figure 9). Other factors we considered in the sensitivity analysis (radiation, LAI, and maximum snow albedo, AlbSnMax) have complex effects on the dynamics of the EnvM and involve negative feedback loops that result in low sensitivity of ALD and WTD (Figure 9). For example, as LAI increases ET will increase, which will tend to result in a drier soil. However, when the soil is drier, the surface infiltration becomes greater, which tends to result in a wetter soil. These negative feedbacks limit the effect of increasing LAI on WTD. Thus, while the model does not demonstrate substantial sensitivity to radiation, LAI, and maximum snow albedo, and AlbSnMax, these factors are still quite important in soil thermal and hydrological dynamics of the EnvM.

Discussion
[40] With respect to evaluating model performance in simulating responses to fire, the tundra (K2) site provided both preburn and postburn measurements to evaluate model performance. In contrast, the black spruce chronosequences represent space-for-time substitutions. While space-for-time substitutions have been a powerful means for studying progressive changes in ecosystems after a major disturbance [Rastetter, 1996], it is important to recognize that differences between sites of a chronosequence, for example differences in mineral soil chemistry and disturbance history, may result in comparisons that are not completely homologous with sites that have predisturbance and postdisturbance measurements. This aspect of space-for-time substitutions complicates model-data comparisons in which it is difficult to attribute whether model-data mismatches are associated with model performance issues or data quality issues. Even with these limitations to model-data comparison, the EnvM demonstrates substantial ability in simulating the dynamics of evapotranspiration, soil temperature, active layer depth, soil moisture, and water table depth in response to both climate variability and fire disturbance. The model is capable of simulating changes in active layer depth that affect water table depth through changes of subsurface drainage, changes that have the potential to influence biogeochemical dynamics. Our analyses also identified several differences between model simulations and field measurements that warrant attention. Below we discuss what we have learned from evaluation of the EnvM with respect to issues involving (1) atmospheric driving variables and (2) model parameters. We also briefly discuss the issue of soil moisture measurements in organic soils.

Issues Involving Atmospheric Driving Variables
[41] In our simulations, we used the estimates of daily air temperature and n factor to estimate ground surface temperature for purposes of estimating heat transferred between the ground and the atmosphere. The n factor, which is the ratio of temperature at the soil surface to that in the air, is used by some studies to estimate ground surface temperature from air temperature [Kiene et al., 2001;Karunaratne and Burn, 2004]. However, the n factor is affected by ground thermal properties and subsurface conditions [ Karunaratne and Bum, 2004], and by vegetation cover [KIene et al., 2001]. The n factor values reported in the literature range from 0.63 to 1. 25 [KIene et al., 2001]. Most ecosystem models do not use the n factor approach to estimate ground surface temperature and assume that ground surface temperature is the same as air temperature that is generally measured at 2 m height. We have used the n factor derived from the measurements at the sites analyzed in this study. Among the Donnelly Flats sites, the n factor is much greater at the recently bumed sites than at the mature sites. At the K2 site, we assumed that the 11 factor of the burned site was equal to 1.00, which might represent an underestimate that could be responsible for the underestimation of summer soil temperatures following fire in 2002. Land surface models calculate ground surface temperatures through calculating ground surface energy balance based on subhour resolution of atmospheric driving data, (e.g., CLM3 [Oleson et al., 2004]). Because ecosystem models generally use daily or monthly resolution atmospheric driving data, they cannot adequately make use of the methodology used in land surface models to calculate ground surface temperature. Given the sensitivity of ALD to air temperature identified in this study, some attention should be given to estimating more realistic ground surface temperature.
[42] The accuracy of precipitation data sets has been an important concern for applications of large-scale hydrological and ecosystem models in northern high-latitude regions [McGuire et al. 2008] because of the low density of observations and numerous issues involved in accurately measuring snow fall [Yang et al., 2005]. We believe that biases in precipitation or issues of redistribution of snow are responsible for some of the differences between simulated and measured soil temperature, e.g., the snow fall at the beginning of 2004 at Delta Junction is lower than snow fall at the beginning of the other years of our simulations.

Issues Involving Model Parameters
[43] The factors analyzed in the study can be broadly divided into three categories: (I) factors that directly affect ALD; (2) factors that directly affect WTD; and (3) factors that indirectly affect ALD and WTD through effects on various processes (Figure 9). Factors that directly affect ALD included warm season air temperature, minimum soil thermal conductivity, and organic soil thickness. An important dynamic identified in our analysis of the model simulations is that changes in ALD often affect WTD through changes of subsurface drainage. However, the strength of this effect depends on whether a site is well or poorly drained.
[44] In the sensitivity analysis we conducted, we did not evaluate all of the parameters in the model that affect hydrology. The parameters defining the hydraulic properties of moss, organic soil, and mineral soil can substantially affect the hydrological dynamics of the model. With respect to 1110SS, there are two main types of bryophyte groups in boreal black spruce forests, fcathermoss (e.g., Hvlocomium G020!:; YI ET AL.: SOIL THERMAL AND HYDROLOGICAL DYNAMICS G02005 splendens). which generally occurs in well to intermediately drained conditions, and Sphagnum, which generally occurs in more poorly drained conditions [Manies et al., 2006).
These moss types have distinct thermal, hydrological, and ecophysiological characteristics [e.g., see Bisbee et al., 2001]. Feathermoss depends heavily on the precipitation and dew formation for water, while Sphagnum can wick water from deeper soil through capillary action. The soil moisture measurement at the FBKS site used in this study provides two data sets, one for soil pits with feathermoss cover, and the other with Sphagnum [Kim et al., 2007). The soil moisture at 30 ern of a "Sphagnum soil" is always greater than that of a "feathermoss soil." In this study we did not implement the different hydrologic properties of these different types of moss. It will be important to consider the distinct hydrological and thermal properties of different types of moss that occur in moderately well drained and poorly drained situations in future applications of thc model.
[45] There are significant differences in hydraulic properties among different types of organic soil [Letts et al., 2000]. Beringer et al. [2001] did not implement differences among different organic soils. We therefore defined porosities for moss, fibrous organic soil, and amorphous organic soil based on field measurements for black spruce forests in Manitoba, Canada [Yi et al., 2009]. In addition to porosity, pore size distribution is another important factor affecting water movement in organic soils [Clapp and Hornberger, 1978;Cosby et al., 1984]. If the same saturated hydraulic conductivity and pore size distribution are used for fibrous and amorphous organic soil, the fibrous organic soil with larger porosity would tend to suck water from the amorphous organic soil horizons below. However, amorphous organic soil usually has much larger values of the pore size distribution parameter than fibrous organic soil, which means that amorphous organic soil can hold water more tightly than fibrous organic soil. In this study, a larger value of pore size distribution (i.e., 6) has been used for the amorphous organic horizon than for fibrous organic horizon (4 based on Beringer et al. [200l]). The proper parameterization of different types of organic soil continues to be a key challenge in the simulation of soil water dynamics.

Issues of Soil Moisture Measurement
[46] A key challenge in evaluating the simulation of soil moisture is confidence in the accuracy and precision of the measurements of soil moisture. We found that the EnvM is better in simulating soil moisture at depth than in simulating soil moisture near the surface. It is currently very challenging to measure soil moisture in SUlfate organic soils because the high porosity creates substantial difficulties in calibrating instruments used to measure soil moisture [Yoshikawa et al., 2004;Overduin et al., 2005]. The potential for substantial measurement error in field-based soil moisture estimates makes it difficult to ascertain what discrepancies between model-and field-based estimates require attention in the model,

Conclusion
[47] This paper presented and evaluated a stable and efficient scheme for incorporating the interactions of soil thermal and hydrological processes in a terrestrial ecosys-tem modeling framework. Our analysis using this framework reveal cd that changes in active layer depth often affect water table depth through changes of subsurface drainage. Because these changes have the potential to influence biogeochemical dynamics, we suggest that it is important for models of carbon dynamics in cold regions to more comprehensively represent interactions between soil thermal and hydrological processes.
[48] Key challenges identified in this study for evaluating! improving model performance include (1) proper representation of discrepancies between air temperature and ground surface temperature; (2) minimization of precipitation biases in the driving data sets; (3) improvement of the measurement accuracy of soil moisture in surface organic horizons; and (4) proper specification of organic horizon depth, hydrological properties, and soil thermal conductivity. Resolution of these observational issues will help in reducing model uncertainties. As changes in organic horizons caused by wildfire have substantial effects on soil thermal and hydrological regimes, an important next step is to investigate how the dynamics of organic soil thickness associated with wildfire disturbance and ecological succession affect the dynamics of permafrost and soil carbon in high latitudes.
Appendix A [49] This part contains the detailed descriptions of the water and energy fluxes among the atmosphere, canopy, snow and soil, and within soil layers. In the environmental module (EnvM; Figure S1 in the auxiliary material), the ground layers include snow layers, soil layers, and rock layers. The soil layers include live moss, fibrous organic soil, amorphous organic soil, and mineral soil. Eleven mineral soil types have been considered in this study following Beringer et al. [2001]. The parameters of the EnvM are provided in Table A1.

Incident Radiation
[50] The monthly short-wave irradiance incident on land surface (NIRR) is provided as input (W m -2). Daily radiation is estimated from monthly values through linear interpolation.

A1.2. Air Temperature
[51) The daily air temperature (T a ) is derived from monthly values using linear interpolation.

A1.3. Precipitation
[52] The daily precipitation (mm/day) is downscaled from monthly precipitation (mill/month) based on the algorithm of Liu [1996], as used by Zhuang et al. [2004]. The daily precipitation is identified as rain when the daily air temperature greater than O°C, otherwise it is identified as snow.  and radiation that is passed through the canopy and becomes incident on the ground surface.

B1.1. Radiation Reflectance to the Atmosphere (Rv2a)
where (xv is vegetation albedo. The albedos of tundra, deciduous, and conifer forest were assigned to be 0.19, 0.19, and 0.10, respectively [Beringer et al., 2005] [55] Beer's law is used to calculate the extinction of radiation in canopy: where E R ( = 0.5) is a dimensionless extinction coefficient of radiation in the canopy. LAl is projected leaf area index. B 1.3. Radiation Intercepted by the Canopy (R,,) [56] B2. Water B2. (1) evaporation of canopy liquid water (Qv,vep) (2) sublimation of canopy snow (Qv,sub), and (3) transpiration (Qvtrs) If there is intercepted rain or snow, the intercepted radiation (Rv) is first used to evaporate rain or sublimate snow. The remaining energy is then used to drive the transpiration formulation. The detailed equations for these fluxes arc described in Appendix D of Zhuang et al. [2004]. Several modifications have been made for the calculation of stomatal conductance: where gs and gs,max, are the stomatal conductance and maximum stomatal conductance (m . s-1), and f(Tmin), f(PPFD), and f( ) are the effects of minimum air temperature. radiation, and leaf water potential on stomatal conductance, respectively. The effect of f(PPFD), which was not considered by Zhuang et al. [2004], was implemented in the EnvM following Waring and Running [ 1998]: where Wi is plant wilting factor for layer i, and r, is fraction of roots in layer i. The plant wilting factor Wi is calculated as: where max is a specified constant value for the wilting point potential of leaves (-1.5 x 105 mm), and sat,i are soil water matric potential and saturated matric potential (mm) of layer i (see Appendix D), and T, is the soil temperature in layer i ( 0 C).

H2.3. Water Fluxes From the Canopy to the Ground
[60] The water fluxes from the canopy to the ground consist of four components including throughfall of snow or rain, and the "dripping" of snow or rain from the canopy after throughfall. The throughfall fluxes of snow or rain arc calculated as the precipitation input of snow or rain minus the interception of snow or rain. The "dripping" fluxes of snow or rain are calculated as the interception of snow or rain minus the sublimation or evaporation of snow or rain.

Appendix C: Effects of Snowpack on Radiation and Water Fluxes
[61] The snowpack is divided into a maximum of five layers, similar to the method used in CLM3 (Community Land Model, V3 [Oleson et al., 2004]). The snowpack accumulates from the fluxes of snow throughfall and drip from the canopy, and it is subject to ablation from sublimation and melt. The dynamics of the snow layers are implemented using a double linked list structure, so that when a snow layer is too thin it is combined with an adjacent snow layer or it is removed if there is no adjacent snow layer. When a snow layer is too thick, it is divided into two snow layers. The criteria for the thickness range of snow layers are similar to those criteria used in CLM3.

Cl. Radiation
[62] The radiation from the snowpack to the atmosphere is calculated as: where Osn is snow albedo, which is calculated following Roesh et al. [2001].

of 20
where PPFD is the absorbed solar radiation by canopy (J x m2 x s-1 and PPFD50 is the value of PPFD at which caused f(PPFD) is 50% saturated. In this study, PPFD50 is set to 75 mol· m-2 . S-1 or 16.5 J x m 2 x S-1.
[59] With the vertical distribution of soil water content calculated, f( ) is modified to include the effects of root distribution on transpiration, following Oleson et al. [2004]: [64] If air temperature is less than -1°C, then snowpack is sublimated using available radiation: where SA = 0.6 is snow absorptivity [Zhuang et al., 2002], Lvap) = 2.51 x 10 6 J . kg-1 is the latent heat of vaporization, and Lfus = 3.337 x 10 J. kg -1 is the latent heat of fusion.

C2.3. Water From Snowpack to the Soil
[65] If air temperature is greater than -1°C, snowpack can be melted by radiation or by heat conduction. The radiation-driven melt is calculated as: The heat conduction-driven melt is described below in E. where Osl is soil albedo, which is calculated following the NCAR LSM [Bonan, 1996] using color classes related to soil types. The color class of moss is assigned to be 3 following Beringer et at. [2001], and the color classes of the fibrous and amorphous organic layers are assigned to be 8, which can be exposed by fire, the darkest class. The overall reflectivity of soil is considered to be the average of reflectivity of visible and near infrared radiation.

D2. Water
[67] The methods of the CLM3 are used to calculate the surface runoff and base flow [Oleson et al., 2004]. The infiltration is estimated as the difference between input, i.e., rain and/or snowmelt, and runoff If the topsoil layers (more than 2) are unfrozen, the finite difference equations are used to solve for the movement of water between layers, with the infiltration as top boundary condition, and zero drainage as the bottom boundary condition.
D2.1. Surface Runoff (Qover) and Infiltration (Qinf l) [68] Surface runoff (Qover mm . S -1) is calculated following method of CLM3 [Oleson et al., 2004]. The rest of input of water is considered as infiltration. Surface runoff is calculated based on the water table depth and fraction of saturation, which is determined by a topography parameter, Fsat. In this study, the default value from CLM3, 0.3, is used.

D2.2. Transpirational Loss of Soil Water
[69] On the basis of the fine root production fraction over soil profile, each soil layer loses soil water through transpiration. The sum of soil water lost to transpiration is equal (0 the transpirational flux calculated for the canory (see Appendix B). (Qsl,evap) of Soil Water

D2.3. Evaporation
[70] The calculation of potential evaporation (Qsl,pe) from the soil surface is same as that described in Appendix A of Zhuang et al, [2002], which is based on the Penman equation, If the rain and snowmelt are greater than Qsl,pe, the actual evaporation is calculated as: where Evr is defined as Evr = 0.3/DSR2, and DSR is the number of days since the last rain. However, experiments have shown that the evaporation from the soil surface is generally underestimated when downscaling precipitation from monthly values to daily values. Therefore, in this study a minimum evaporation ratio (0.08) is assigned, based on the observed water fluxes from DFT99 site, which was burned in 1999.

Water Movement Between Soil Layers
[71] On the basis of the calculations of evaporation.
infiltration, and transpiration, the movement of soil water between soil layers is calculated by solving Richards' equation. Because there arc cases for which the updated soil water content of a soil layer can become negative or more than porosity at a daily time step, an iterative method is used to solve Richards' equation. The initial time step at beginning of each day is a half day. If after a time step any soil water content is not reasonable (negative or greater than porosity), then the time step is halved, and this continues until the solution is reasonable. The solution continues until a 1-day period is fully covered.
where is saturated matric potential (mm), b is the Clapp and Hornberger [1978] constant, g is gravitational acceleration (m . S-2), Tsl is soil temperature (K), and Tf; is freezing point of liquid water (K). where Ksat is saturated hydraulic conductivity (mm . S·l).

Sat,
and b are specified for each type of soil following Beringer et al. [2001]. (Qbf) [73] After the soil water contents of all layers are updated, the base flow is calculated following method of CLM3 [Oleson et al., 2004]. Drainmax is an important parameter determining the maximum rate of drainage from unsaturated fraction, in this study, the default value of 0.04 mm/s, is used. In the CLM3, only the deep soil, 6-10 layers, has base flow. In this study, if the mineral soil layers, deeper than 80 cm, are unfrozen, then there can be base flow from these layers, otherwise the base flow is O. For poorly drained soils, the base flow is assumed to be O.

D2.S. Base flow
D2.6. Water Table Depth [74] The WTD (depth from ground surface to groundwater surface, m) is calculated following the method of Oleson et al. [2004]: where zhn is the distance between soil surface and the bottom of soil layer n, .1'; is soil wetness of soil layer n, and Zi is the thickness of soil layer n. [77] If a layer is frozen, then phase change will happen in this layer. TDSA calculates where ddneed,i is the degree day (Oe day) needed to completely thaw layer i, Oi is the volumetric water content of layer i, which is calculated following the method of Johansen [1975]: where Ke is Kersten number, which is related to soil water content. The parameters ksat,i and kdray,i are saturated and dry thermal conductivities of a soil layer, which are specified for each soil layer type.
[78] If ddneed,i is less than ddleft the frozen state is changed to unfrozen, and a thawing front is moved to the top of next layer, and ddleft is recalculated: Rsum,i is then updated by adding the thermal resistivity of the current layer to the old value of Rsum,i, and the TDSA then proceeds to the next layer if it is not rock.
[79] If ddneed,i is greater than ddleft, where dpart is partial depth that can be thawed, and kfrzi,i is the frozen thermal conductivity of layer i (J/Kms). A thawing front will be created at a depth dl'ort relative to the top of layer i. ddleft is set to zero and the iteration stops, [80] If a layer is unfrozen, the ddleft is kept unchanged, Rsum,i; is updated and the TDSA then proceeds to the next layer if the layer is not rock. After the movement of thawing front downward, the same procedure is used to adjust the deepest front upward.

E2. Temperatures of All Layers
[81] After the positions of FTFs are determined, the temperature of each layer will be updated. If there is no front in whole ground column, the temperature of each layer G02015 YI ET AL.: SOIL THERMAL AND HYDROLOGICAL DYNAMICS G()2()15 will be updated by solving finite difference equations of all layers, with the derived ground surface temperature from air temperature as the top boundary condition. Because the heat flux around 50-100 m can neglected for a time period of centuries [Nicolsky et al., 2007], we assumed a zero heat t1ux as the lower boundary condition. Ifthere is one front in whole ground column, the layers above and below the front will be updated separately by solving two different sets of equations assuming no phase change. If there are two or more fronts in the whole soil column, the temperatures of layers above first front will be updated by solving the finite difference equation of layers above first front, and a similar method will be used to update temperatures below the last front. For layers between first and last front, the temperatures are assumed to be 0 o c.
[82] The Crank-Nicholson scheme is used to solve the finite difference equations of ground temperatures. To keep the calculation stable, an adaptive step-size integration approach is used. The initial time step is a half day. After advancing one time step, if the change of a layer's temperature is greater than a specified threshold (0.1 o C in this study), the time step will be halved. Iteration continues until the calculation covers a full day.
where Psn is the snow density in unit of kg x m -3, and the value of Psn ranges from new snow density Psn,new to maximum snow density Psn,max following Verseghy [1991], with a time step of 1 day: where T = 0.24 corresponds to an e-folding time of about 4 days. PSn,new is specified to be 100 kg x m -3, while Psn,max is the maximum snow density, which is specified to be 250 kg x m-3 for forest ecosystems [Pomeroy et al., 1998] and 362 kg x m-3 for tundra ecosystems [Ling and Zhang, 2006].
where K; is Kersten number, and S, is soil wetness. The parameters ksl,sat, which is the saturated thermal conductivity, and ksl,dry, which is the dry thermal conductivity, are calculated as:

R2.2. Soil Thermal Conductivity
[84] Soil thermal conductivity ksl (W x m-1 x K-I) is calculated following Farouki [1986]: where kwlid, ksolids, and kliq arc thermal conductivity of solid material, liquid, ice, and air, respectively. The parameters 0sats, 0 liqa and 0ice are porosity, volumetric liquid content, and volumetric icc content, respectively. The parameters ksolid and 0sat are specified for each type of soil following Beringer et al. [2001]. In some cases, the simulated ksl is too small for topsoil layers under very dry conditions. Therefore, the minimum thermal conductivity has been set 0.1 W m-I K-I following Hinzman et al. [1991] to account for nonconductive heat exchange.

E2.3. Rock Thermal Conductivity
[85] The thermal conductivity of rock may vary by as much as a factor of two or three for the same type of rock. Because sedimentary rock underlies 72.9'% of Alaska [Peuck.er-Ehrenbrink and Miller, 2003], we used a mean value of 2 W x m-I x K-I , the thermal conductivity of bedrock following Clauser lind Huenges [1995].

E3. Volumetric Thermal Capacity
[X6] The volumetric thermal capacity (J x m -3 X K-1) is calculated based on the specific heat capacity (J x kg -I X K -I) and bulk density (kg x m -3). The specific heat capacities of soil are related to soil type, that of snow is 2117.27 J x kg-1 X K-1 , and that of rock is 1000 J X kg-1 X K-1. The bulk densities of soil arc related to soil type, the bulk density of snow is calculated during simulation, and the bulk density of rock is set to 2700 kg x m -3 [87]