Modeling Freezing and Thawing of Subsurface Soil Layers for Global Climate Models

Freezing ground, including permafrost and seasonally-frozen ground, is found in more than half of the Northern Hemisphere (NH). The changes in soil freezing/thawing, or the subsurface thermal conditions, are a local phenomenon, but can have large impacts on climate and life in cold regions. Permafrost (soil is at or below 0°C for more than two consecutive years), occupying about 20% of the exposed NH land, is an essential component in the hydro-eco-climate system in the Arctic. It defines and controls biogeochemical activities, types of dominant vegetation, hydrological states and conditions (e.g., water holding and storage capacity, river and channel routing, limnology) in the regions. However, the influences of local subsurface hydrological-thermal regime changes are not confined to the high-latitude regions. Multiple pathways, which are physical, ecological and/or biogeochemical, can transfer the subsurface changes to other regions, and affect the regional to global climate with different consequences. The carbon and nitrogen cycles are closely controlled by the presence of permafrost, degradation of which may cause an irreversible change in the cycle and accelerate the emission of greenhouse gases, such as carbon dioxide and methane. Some types of vegetation may not survive when the underlying permafrost degrades, leading to migration of other vegetation types. Consequent changes in albedo and evapotranspiration functionality will shift the budget of energy and water on large scales. Precipitation in the high-latitude, continental summer owes strongly to water recycling capability. Water availability in summer, therefore, may be altered following such changes in the surface and subsurface hydrological conditions.

Development and examination of permafrost dynamics in numerical models started in the early days of climate modeling. The performance and the sensitivity of land process models with the freeze/thaw process have been examined in the context of climate studies since the 1990s within the global climate modeling community (e.g., Riseborough et al., 2008), through both offline and in coupled simulations. Saito (2008b) examined, in a series of onedimensional offline experiments for a tundra site, the sensitivity of the subsurface thermal and hydrological regimes to the complexity of the resolved physics under present-day conditions and under a warming scenario. The importance of the physical refinement was demonstrated in terms of the effect of freezing on the evaluation of hydrological-thermal property values, presence of the top organic layers, and presence of unfrozen water under the freezing point. Studies by coupled simulations evaluated the comprehensive performance of the climate models, as a whole, that incorporate the freeze/thaw process (e.g., Lawrence et al., 2008 and Saito et al., in preparation). Such evaluations, however, would not make sense if the used land process models are unable to reproduce the historical subsurface regimes when forced by historical data.
In this chapter, the author aims to evaluate the performance of land process models with different complexities of soil freeze/thaw dynamics, regarding reproducibility of surface and subsurface hydrological-thermal conditions during the latter half of the last century, and its sensitivity to the implemented physics. However, there are not so much observational data for the subsurface hydrological or thermal states that cover the spatial and temporal extent matching large-scale climate model outputs. Therefore, the comparison and validation tend to be limited in time or space, or indirect.

Model, data and experimental design
Two sets of global-scale offline simulations were conducted using different types of forcing data. One set of simulations is a hindcast where the models were forced by a compiled dataset of the historical, observed surface meteorology. Another set is a climatology forced by the climatological values derived from historical values. The land process model employed in this study is called MATSIRO (originally developed by Takata et al. 2003, and www.intechopen.com improved in terms of freeze/thaw dynamics by Saito 2008aSaito , 2008b, a component of the global climate model MIROC (K-1 developers, 2004) developed in Japan. The refinement of MATSIRO for surface and subsurface hydrological-thermal regimes was documented by Saito (2008). The design of the sensitivity examinations is described in 2.1. The forcing data and other data used for evaluations and validations are explained in 2.2.

Experimental design and boundary factors
Four different experiments were conducted to examine the difference and sensitivity of the relevant configurations and factors (Table 1), i.e., a) total depth, total number of layers and thickness of layers in the soil column, b) presence or absence of top organic layers (TOL; e.g. a litter layer and a humus layer), c) basic value of the soil characteristics, such as thermal conductivity and heat capacity for different soil types, and d) freeze/thaw dynamics parameterizations. The last factor includes evaluation of thermal and hydraulic properties such as thermal and hydraulic conductivity, heat capacity, and presence of unfrozen water under the freezing point.

Exp.
Soil column depth and thickness  Figure 1 illustrates the distribution of prescribed vegetation (left) and soil characteristics (right). Both distributions were basically the same as the ones used in Saito et al. (2007), but are adapted to the 1948 cropland condition (Ramankutty & Foley 1999) and fixed through the simulations. They were common to all the experiments, except that top organic layers were neglected for CNV and noTOL14 runs. Consult Takata et al. (2003) for details of each category of vegetation and soil characteristics. Figure 2 shows the annual range of soil temperature for the bottom layer (i.e., for 4-14 m for the 6-layer runs and 20-100 m for the 12layer run). It demonstrates that the soil column needs at least 20 m or more to be consistent with the specified lower boundary condition of a constant (often set to zero) heat flux.
Vertical profiles of soil thermal conductivity (lines with symbols) and soil moisture (simple lines) under different frozen ground conditions are shown in Figure 3. Figures 3a and 3c show results for grid points close to the southern permafrost region in North America and Eurasia, respectively. Figure 3b is for the northern part of permafrost region, and Figure 3d is for a seasonally-frozen region. The figure illustrates that the improved parameterization (i.e. TOL14 and TOL100 runs shown in green and blue, respectively) capture the feature well, for example, the summer thermal conductivity values are generally lower than those for winter because solid ice has higher conductivity than liquid water. The conventional model, the CNV run (shown in orange), shows variations following the total moisture but is www.intechopen.com Climate Models 212 less sensitive to the ratio of soil moisture and ice (difference between values in summer and in winter). The vertical structure of thermal conductivity is not clear for the cases of mineral layers only (i.e., CNV and noTOL14). The cases with top organic layers (i.e., TOL14 and TOL100), however, reproduce the vertical structure of the soil property that is comparable to the observations (shown in black in Figure 3a).

Data used for forcing and evaluations
The historical atmospheric data covering different periods in the 20th century were collected from different data sources and used to evaluate statistical uncertainties resulting from the forcing data. A global dataset of meteorological forcing that was developed, compiled and provided by the Land Surface Hydrology Research Group, Princeton University was used (Sheffield et al. 2006). The 3-hourly data at 1.0-degree resolution was interpolated into 1-hourly, T42 horizontal grids (c. 2.85x2.85 degree grids in the direction of longitude and latitude) for the experiments for the simulation period from 1948 to 2000. The climatology of the forcing data was derived from the first thirty years of the period.
The permafrost map, compiled by the International Permafrost Association (IPA) (Brown et al. 1997) and referred to as IPA map hereafter, is consulted to validate the modeled distribution of permafrost and the ground ice volume. River discharge of the selected large arctic basins was compared with the R-Arctic Net dataset collected and compiled by the University of New Hampshire (R-ArcticNET v4.0). Ground temperature distribution in the former Soviet Union was also used; it was interpolated from the observations at more than 200 stations onto the T42 grid at the National Snow and Ice Data Center (NSIDC) (described in Saito et al. 2007). Surface air temperature data was taken from the 0.5 degree latitude/longitude dataset of monthly surface temperature CRU05 produced by the Climate Research Unit, East Anglia University (New et al. 2000). Observational values of active layer depth were taken from the Circumpolar Active Layer Monitoring (CALM) dataset (Brown et al. 2003). The vertical profile of soil thermal conductivity and soil moisture was measured at a coastal tundra site on the Seward Peninsula in summer of 2010 by the chapter author. The Decagon KD2 Pro thermal analyzer and Campbell Scientific Hydrosense were used to measure the respective aforementioned quantities.

Results and discussions
In 3.1, the hindcast runs, which used historical values of the forcing data compiled by Sheffield et al. (2006), were analyzed to highlight how the modeled surface heat budget (3.1.1), subsurface thermal regimes (3.1.2) and hydrological budget (3.1.3) in the cold regions are dependent on the different modeling configurations and different implementations of physics. In addition, 3.2 shows results of the climatological runs, forced by 1948-1978 climatology derived from the above historical data. The climatological runs produced a repeating two-year cycle in which snow accumulation is markedly larger in one year and smaller in the other. One arbitrary two year sequence (after sufficient spin-up integrations) was used for analysis to demonstrate the effects of snow on the subsurface thermal and hydrological regimes, and how sensitive those effects are to the complexity of the resolved land processes.

Historical hindcast from 1948 to 2000
A historical hindcast from 1948 to 2000 (53 years) was integrated after a 500 year spin-up under 1949 conditions, which was used instead of 1948 to avoid the leap year. All soil layers reached equilibrium in all the runs after the spin-up. No ensemble was made for this study. www.intechopen.com

Influences on energy budget at the surface
From the viewpoint of large-scale climate model development, if the subsurface model performance has little or no influence on the atmosphere, the complexity of the physics implemented for the land process models does not matter. Impacts on surface budgets of either energy or water, or both, are evidence that subsurface models do matter to the total climate model. The marked difference and sensitivity found in winter between the experiments is in the albedo and, consequently, the shortwave radiation budget due to different snow accumulations. Figure 4 shows the surface albedo (%) and snow water equivalent (kg m -2 ) in February. In most of the NH land snow accumulation (or for the sake of radiation budget, snow cover extent) reaches its maximum in January to February, and the solar radiation is stronger and spread wider in February than i n J a n u a r y . I n o r d e r t o e m p h a s i z e t h e comparison and due to space limitations, the figures in this subsection show only the results for the TOL14 run and the deviation of the conventional model (CNV-TOL14), but the presence of the top organic layers was the major factor of the difference (not shown). Snow accumulation and, hence, diagnosed snow cover area were larger in the northern highlatitude land for the TOL runs, due to cooler soil temperature below the surface at the beginning of snow accumulation season. A larger snow area led to a higher albedo in winter that potentially cools the near-surface atmosphere if the atmosphere is interactive through reduced absorption of the solar energy at the surface. The areas of large difference for the TOL runs are northern Europe, the European p a r t o f R u s s i a , K a m c h a t k a , t h e R o c k y Mountains and northeastern Canada for the higher albedo, and southern Siberia and British Columbia for the lower albedo. In warmer seasons, the eddy heat fluxes (i.e., sensible and latent heat flux) become more important and enhanced than in cold seasons. Figure 5 shows a result for the summer conditions in which the surface heat balance was different between the experiments because the amount and the partition of surface available energy was changed: evapotranspiration (latent heat flux) was enhanced in permafrost regions, i.e., northeastern Eurasia and northern North America (Figure 5a-b), and sensible heat flux increased in large areas in the northern high-latitude land areas for the refined models (Figure 5c-d). One of the important factors that caused this change was the reduction of geothermal heat to the ground. In spring to summer the geothermal heat flux i s n e g a t i v e ( d o w n w a r d , i . e . , f r o m t h e atmosphere to ground; Figure 6a). The lower thermal conductivity value for the top organic layers hinders heat conduction to the ground in comparison to the only mineral soils ( Figure  6b). Convergent heat at the surface can increase sensible heat fluxes. Another reason is availability of water in upper soil layers in permafrost regions (see profiles of thermal conductivity and soil moisture for summer shown by solid lines, in Figures 3a-c). Layers of mosses or dead organics can hold water more effectively to make it available for evaporation from the ground surface, or for plants to transpire.
Geothermal heat flux showed a marked difference in two transitional seasons, early summer and early winter ( Figure 6). Early in the cold season, the geothermal heat is generally positive (upward, i.e., from the ground to the atmosphere; Figure 6c). The conventional scheme transferred more heat from the atmosphere in summer and to the atmosphere in winter so that the annual amplitude of ground temperature was large. The refined parameterizations take the amount of ground ice into account in computing the working thermal conductivity, which makes wintertime thermal conductivity larger even when the total soil moisture (both in solid and liquid) is the same between summer and winter. This makes seasonal asymmetry in conduction of ground heat leading to different vertical profiles and seasonality of modeled annual ground temperature, as discussed in Saito (2008b). The changes in the surface heat balance and heat conduction influence the subsurface thermal regimes differently in different regions, and modify the distribution of frozen ground (permafrost, seasonally frozen, and no freeze regions) as shown in 3.1.2.

Evaluations of subsurface thermal regimes
The current thermal state of permafrost from the last century was widely investigated and documented during the International Polar Year (IPY 2007(IPY -2009. Permafrost in Eurasia  and North America (Smith et al., 2010) is stable in the northern cold regions, but shows a warming trend in the southern warmer regions.
Thermal state means a seasonal cycle and distribution of temperature. Figure 7 shows an example of modeled subsurface temperature (Tg) distribution and a comparison with an observation-based estimate for the regions of the former Soviet Union. Notice that the distribution of modeled Tg is very close to that of the forcing data (surface air temperature at 2m, T2; Figure 7c). The difference between the prescribed T2 and the modeled Tg was mostly in the range between 2 to 4 o C, and Tg was warmer. This difference is due to several agents, such as vegetation, snow cover, and seasonality of water storage. Estimates of subsurface temperature distribution, such as the one shown in Figure 7b, are neither common nor easy to produce. It is a valuable map, but has weaknesses too. As discussed in Saito et al. (2007) the regions of permafrost, e.g., eastern Siberia, feature a less dense network of observation stations so that the estimated values have large uncertainty. In the seasonally frozen regions (the southwestern part of the map) the correspondence between the modeled and estimates are good, but the difference is large in the northeastern part. Judging together with other observational evidence (e.g. Figure 9f in this chapter) the estimates are too warm in the region. The active layer is "the layer of ground that is subject to annual thawing and freezing in areas underlain by permafrost" (van Everdingen (ed), 1998 ). The maximum depth of the active layer is usually reached in the late summer, and referred hereafter as maximum ALD or mALD. Change of mALD occurred due to several reasons, both natural and anthropogenic. Temperature change, as well as surface disturbances such as wildfire, construction, and cultivation, can modify mALD. Figure 8 shows examples of the evolution of mALD, both modeled and observed, at two different sites characterized as permafrost.

www.intechopen.com
Although the average and the interannual variations of the forcing T2 (thin black lines) do not differ from the amplitude of the interannual variations in the observed (thick black lines. only shown as anomaly from the average) and modeled mALD (colored lines) are very different. In Figure 8a the CNV run showed a larger interannual amplitude than the TOL14 run. This is because the large and rapid conduction of heat to ground for the CNV run reached to a deeper layer in some years but not in other years. In the TOL runs, the heat conduction is slower and damped, and the depth of mALD is more stable. In Figure 8b the maximum active layer depth is deeper than in Figure 8a, and, since the CNV run allows no unfrozen water, once a layer is frozen it stays at 0°C until all water freezes to ice or ice melts away completely, and the table of permafrost stays at the same depth. The refined parameterization, however, allows unfrozen water and the ground temperature can be below 0°C while some water stays in liquid phase, so that the depth of the maximum ALD can fluctuate over years. At both points, the mALD was shallower in the 1960s and the early 1970s, and deepened after the mid 1980s. All the runs with the refined physics showed this similar trend, but the TOL14 run computed an abrupt deepening of ALD in Figure 8b Scale is shown on the right axis) at a) 82E, 72N, and b) 130E, 60N in the continuous permafrost region. Also shown is the forcing T2 at the points (Scales shown on the right axis. Note the values need to be multiplied by ten). Figure 9 shows the geographical distribution of mean maximum ALD for the integration period. The CNV run computed the deepest ALD of the four, ranging from 1m along the northern coast and more than 3m in some southern regions. The noTOL14 run simulated somewhat shallower ALD, but had a wide zone of ALD deeper than 3m compared to the CNV run. This is again because of the presence of unfrozen water (or continuous transition at and below the freezing point), which can allow a gradual change of ALD. The experiments with top organic layers modeled much shallower ALD (Figures 9c-d). Figure 9e is derived from the NSIDC soil temperature data for the former Soviet Union, and 9f is taken from the CALM dataset. Since the CALM dataset is directly taken from the observations (no spatial or temporal interpolation) it is one of the most trustable sources. The result from CALM, in comparison with Figure  Those differences in the simulated thermal conditions (temperature and heat flow) led to a difference in the modeled distribution of the thermal regimes, and frozen ground ( Figure 10 and Table 2). The regions of modeled permafrost were defined by whether any of the subsurface layer stays frozen for two or more consecutive years. Regions were classified as seasonally-frozen if some soil layer freezes at any time of a year but excluded permafrost.   (Brown et al., 1997) for permafrost, and derived from observed monthly surface air temperature for seasonally frozen ground (cf. Saito et al., 2007), and mapped on to the T42 grids.
An observation-based map (Figure 10e) was derived from the IPA map for the permafrost region by choosing the continuous (90-100%) and discontinuous (50-90%) permafrost regions. Distribution of seasonally-frozen ground was diagnosed from the CRU temperature data following the methodology of Zhang et al. (2003) as described in Saito et al. (2007). In general, the modeled permafrost extent is larger than the values from the observation-based one, while the area of seasonally frozen ground is smaller for the simulations. Since the surface air temperature was prescribed (as forcing), the southern limit of seasonally frozen ground was mostly fixed for the modeled case in these offline simulations. It is of interest to investigate whether the area of seasonally frozen ground would change when the land process models are coupled to the atmospheric (and even to the oceanic) climate models. The results of the coupled simulations are compiled in another paper.

Hydrological evaluations in the cold regions
Discharge of the river basins to the Arctic Ocean shows a significantly increasing trend since the last century (Peterson et al., 2002;White, D., et al., 2007). The reason of the increase is not clear yet, but the possible agents include changes in atmospheric moisture conversion and precipitation, dam construction, permafrost thaw, and wildfire (Yang et al., 2003;McClelland et al., 2004). Possible attribution of the permafrost thaw is from changes in the water content held in permafrost through changes in the area of permafrost and the depth of the active layer. Figure 11 shows the difference of the averaged water content between experiments (Figures 11a, 11c and 11e). The CNV run showed smaller water content in the high latitudes, partly because the specified porosity (volume not occupied by soil) for mineral soil was smaller and because infiltration of water at the beginning of the cooling season was small in the conventional model. Figure 11 also shows the distribution of ground ice. The model difference is most prominent between the conventional and the refined one (11a and 11b), ground ice content tended to larger for the CNV run despite smaller water content as a whole in permafrost. This is because some percent (up to 20% or more, depending on the temperature and soil characteristics) of water remains unfrozen in the refined models.

www.intechopen.com
The meridional profile of the modeled ground ice was compared with the observation-based estimate in Figure 12. The experiments successfully captured the general characteristics of the latitudinal distribution in the northern hemisphere, but the absolute value diverts roughly by a factor of three. Some parts of the discrepancy (i.e., overestimation) arguably resulted from neglect of sub-grid scale heterogeneity of land profile. For example, a grid categorized as grassland may have small desert, lakes or rocky hills that tend to have limited ground-ice content compared to a grassland. In turn, the classification of ground-ice content in the IPA map is susceptive to be more qualitative than quantitative in areas with sparse observational evidence. Historical changes in the modeled and observed annual river discharge, together with those in the forcing precipitation, are shown in Figure 13 for the selected arctic basins, and the result of analyses for the attribution of the changes are summarized in Table 3. The observational data are taken from the R-Arctic Net database (R-ArcticNet v4.0). For those river basins underlain in parts by permafrost, such as the Lena, Mackenzie and Anadyr Kolyma, the correlation of the modeled and observed river discharge and precipitation was significant, although the averaged level of discharge tended to differ between the experiments. Those river basins mostly in the seasonally frozen region showed little or no correspondence between the modeled and the observed.
Since precipitation (water input from the atmosphere) is common to all the experiments and basins, there should be other reasons (or mechanisms), e.g., evapotranspiration, or water storage, that determine the discharge from those basins. The lower row of the analysis results in Table 3 shows the correlation between the modeled and observed discharge after removing precipitation. It serves as an index for possible sources of fluctuation other than precipitation.

Snow effects on hydrological-thermal regimes
Seasonal snow cover is one of the major players in cryosphere-climate interactions. It influences the surface energy balance between the atmosphere and land through its high albedo. It also influences the state of the subsurface hydrological and thermal regimes, and its seasonal change though its insulation effect and storage of water during winter. The climatological simulations all showed a two-year cycle with a high snow winter and a low snow year. Snow accumulation (denoted by snow water equivalent, SWE) is a prognostic variable in the model, i.e., computed within the model instead of externally prescribed or forced. Reproduction of the amount and distribution of snow accumulation is one of the areas in which current global climate models are limited. The produced two-year cycle is an artificial phenomenon of the land process model, flip-flopping between two stable states, when forced by climatological data (i.e., no interannual variations). Nevertheless, this gives a unique opportunity to investigate the effects of snow reflected on the subsurface hydrology and thermal regimes.  Table 3. Pearson's correlation between the observed and simulated annual discharge (upper row) and the partial correlation with the effect of precipitation removed (lower row) in the selected arctic basins. Statistically significant correlation coefficients at the 99% confidence level are shown by bold fonts, and those at the 95% level are shown by italic fonts. Numbers in parentheses are the percentage of the modeled discharge area relative to the discharge area at the corresponding station. Table 4 and Figure 14 summarize the effects of snow on the subsurface hydrological and thermal regimes by the four experiments summarized in Table 1. The grids 14a, 14b and 14c are from permafrost regions, and the grid 14d is from a seasonally frozen region, at which the refined model showed no large difference among experiments. Note that the scales for snow are different among the figures.
Snow accumulation (snow water equivalent) is largest at grid 14b (northwestern Siberia) with the monthly maximum reaching more than 100 kg m -2 even in a low year, whereas grid 14a (Interior Alaska) has only 8 kg m -2 , even in a high snow year. Variations in snow between the experiments are small at all the points, as well as ground temperature for the top soil layer. In contrast, geothermal fluxes and soil water content are highly variable among the experiments. As for the effect of snow, water content and ground temperature show a large difference in summer, while geothermal fluxes have greater sensitivity in spring and autumn.
An observational fact is that the maximum active layer depth (mALD) is more correlated with the summer temperature than the snow amount from the previous winter. The maximum depth of the active layer was larger and soil temperature was higher in summer following the high snow year at all permafrost grids (Table 4a), so it could not be examined here. mALD is largest by a factor of two or more for the CNV run among the experiments, and the difference in snow conditions was also the largest.
The average amount of annual local runoff was proportional to the average SWE, but interannual variations were not proportional to the snow fluctuation (Table 4b), which is partly explained by the differing amount of soil water content after snowmelt (more water content after high snow winter, Figure 14), suggesting the buffering (or smoothing) functionality of soil.
The conventional model showed the largest amount of cumulative geothermal heat conducted both from and to the atmosphere, while it is least sensitive to the effects of snow (Table 4c). Upward flux is larger after high snow accumulation at all grids, but the downward flux varies among the grids and among the experiments. , and soil moisture content (m m -1 ; thin lines) for the top layer (0-5 cm for the 6-layer runs and 0-1 cm for the 12-layer run). Scales are on left and right axes, respectively. The maximum active layer depth (x10 cm) is shown by an asterisk. A year cycle begins in September and is repeated three times (high snow-low snow-high snow years). The color scheme is the same as in Figure 2.

Conclusion
Two sets of global-scale offline experiments were run with land process models of different complexity of freeze/thaw dynamics and related subsurface conditions. One set is a hindcast with historical forcing data, and the other was a climatological run. The simulated surface and subsurface variables were compared to each other for the sensitivity of model conditions, and with observation-based data for validation. It was demonstrated that subsurface physics can affect the atmosphere through the surface energy budget, that the subsurface regime is simulated more correctly using top organic layers and the refined freeze/thaw dynamics. The hydrological regime and snow effects showed significant sensitivity to the model complexity, although the results remained descriptive and qualitative, rather than quantitative, partly due to an immature observation-model comparison framework.