LGM permafrost distribution : how well can the latest PMIP multi-model ensembles perform reconstruction ?

Here, global-scale frozen ground distribution from the Last Glacial Maximum (LGM) has been reconstructed using multi-model ensembles of global climate models, and then compared with evidence-based knowledge and earlier numerical results. Modeled soil temperatures, taken from Paleoclimate Modelling Intercomparison Project phase III (PMIP3) simulations, were used to diagnose the subsurface thermal regime and determine underlying frozen ground types for the present day (pre-industrial; 0 kya) and the LGM (21 kya). This direct method was then compared to an earlier indirect method, which categorizes underlying frozen ground type from surface air temperature, applying to both the PMIP2 (phase II) and PMIP3 products. Both direct and indirect diagnoses for 0 kya showed strong agreement with the present-day observation-based map. The soil temperature ensemble showed a higher diversity around the border between permafrost and seasonally frozen ground among the models, partly due to varying subsurface processes, implementation, and settings. The area of continuous permafrost estimated by the PMIP3 multi-model analysis through the direct (indirect) method was 26.0 (17.7) million km 2 for LGM, in contrast to 15.1 (11.2) million km 2 for the pre-industrial control, whereas seasonally frozen ground decreased from 34.5 (26.6) million km2 to 18.1 (16.0) million km2. These changes in area resulted mainly from a cooler climate at LGM, but from other factors as well, such as the presence of huge land ice sheets and the consequent expansion of total land area due to se -level chang . LGM permafrost boundaries modeled by the PMIP3 ensemble – improved over those of the PMIP2 due to higher spatial resolutions and improved climatology – also compared better to previous knowledge derived from geomorphological and geocryological evidence. Combinatorial applications of coupled climate models and detailed stand-alone physical-e ological models for the cold-region terrestrial, paleo-, and modern climates will advance our understanding of the functionality and variability of the frozen ground subsystem in the global ecoclimate system.


Introduction
Frozen ground (i.e., permafrost and seasonally frozen ground) constitutes a critical environmental subsystem of the Arctic eco-climate, closely linked with snow and vegetation (Nelson, 2003;Saito et al., 2013a).Frozen ground can affect local and remote regions through interactions and feedbacks between energy, momentum, water, and materials (carbon, nitrogen, etc.).This subsystem has shown an enhanced and wide response to the recent global warming trend (Romanovsky et al., 2010a, b;Smith et al., 2010;Christiansen et al., 2010).Understanding its functionality, variations, and stability under a different (e.g., glacial) and/or changing (e.g., present-day) climate is a crucial element in Earth science research, as well as in social sciences and policy making (UNEP, 2012).The areal extent of permafrost during the Last Glacial Maximum (LGM) has been mapped using regional to local evidence (e.g., Baulin and Danilova, 1984;Velichko, 1984;Frenzel et al., 1992;Baulin et al., 1992;Petit-Maire et al., 2000;French, 2007), and the northern and western parts of Eurasia and the northern parts of North America and Greenland have been defined relatively well (Vandenberghe et al., 2008;Hubberten et al., 2004;Rozenbaum and Shpolyanskaya, 1998;Washburn, 1980), especially when compared with other areas such as northeast Asia (Vandenberghe et al., 2004(Vandenberghe et al., , 2012;;Saito et al., 2012Saito et al., , 2013b) ) and the Southern Hemisphere (Trombotto, 2008;Saito et al., 2012).
Numerical reconstructions of subsurface temperature regimes from past glacial periods to the present day have been investigated in Alaska, Europe, and Siberia, with comparison to borehole temperature or permafrost thickness data (Lachenbruch et al., 1982;Osterkamp and Gosink, 1991;Delisle, 1998;Delisle et al., 2003;Sueyoshi and Hamano, 2003;Kitover et al., 2012).While these analyses have provided minute information about permafrost evolution at sites under climatic changes, they are not easily extended to wider areas such as continents or hemispheres, partly owing to scarce networks of deep borehole data and to high heterogeneity of the terrain.
Large-scale climate models, from generations of general circulation models to advanced global climate models (hereafter, GCMs) and Earth system models (ESMs), have been used for global-scale paleoclimate studies since the 1970s (Gates, 1976).However, less attention has been paid to terrestrial dynamics, including those of high-latitude regions, until recently.At the time of the Paleoclimate Modelling Intercomparison Project phase II (PMIP2; Braconnot et al., 2007), subsurface variables such as soil temperature and soil moisture had not been archived within the common database.It was in the 2000s that the subsurface regime modeled by GCMs and/or ESMs began to undergo intensive analysis (Renssen et al., 2000;Renssen and Vandenberghe, 2003;Lawrence and Slater, 2005;Roche et al., 2007;Saito et al., 2007).Subsequently, the shortcomings and problems of the processes implemented in the earlier-generation GCMs have been widely addressed and improved (e.g., Nicholsky et al., 2007;Alexeev et al., 2007;Saito, 2008;Lawrence and Slater, 2010).During the third phase of PMIP (PMIP3 2011), conducted in tandem with the fifth phase of the Climate Modelling Intercomparison Project (CMIP5; Taylor et al., 2012), a number of subsurface variables (for soil columns ranging from 3 m to 44 m in depth) have now been archived in the common database, made publicly available, and used for analysis (e.g., Koven et al., 2013).
When and where subsurface temperature is not available, mean annual surface air temperature (MAAT) is a surrogate variable commonly used to interpret and distinguish types of underlying permafrost (Washburn, 1980;French, 2007;Matsuoka, 2011;Gruber, 2012;Boeckli et al., 2012).Van Vliet-Lanoë (2009) provides the rough criteria: MAAT < −7 • C for permafrost zones covering more than 80 % of the area, −7 • C < MAAT < −3 • C for zones with permafrost coverage of less than 80 %, and −3 • C < MAAT for patchy permafrost.However, threshold values are different for different application areas, periods, and types of permafrost (for example, Washburn (1980) cites −2 • C as the upward limit for permafrost, and uses the range of −5 • C to −10 • C for the presence of continuous permafrost for different regions).This is partly a manifestation of subsurface thermal regime not being a function of air temperature alone, but it also depends on other factors such as snow cover, vegetation, soils, and micro-topography (French, 2007;Saito et al., 2013a).Still dependent on atmospheric thermal information alone, however, the use of freezing and thawing indices (defined as cumulative degree-days below and above 0 • C, respectively) have been attempted in order to improve the classification of frozen ground types (Harris, 1981(Harris, , 1982;;Anisimov and Nelson, 1997;Saito et al., 2009).Saito et al. (2009) utilized the PMIP2 surface air temperature outputs to derive freezing and thawing indices, in order to estimate frozen ground distribution at 0 kya and LGM, and assessed these estimations using evidence-based maps.Partly due to coarseness of the horizontal resolution, these comparisons produced mixed results consistent with evidence in some regions but not others.Unsuccessful regions include north of the Alps, high mountainous regions such as the Tibetan Plateau in the Northern Hemisphere, and the Andes and Tierra del Fuego in the Southern Hemisphere.These discrepancies, however, also resulted from known biases in PMIP2 products (e.g., warm bias in the western Europe winter; Ramstein et al., 2007;Braconnot et al., 2007).Upon the simulation outputs of PMIP3 experiments becoming publicly available, this is a good occasion for assessing the models' ability to reconstruct LGM frozen ground.There have been efforts from the observational side for an Action Group of the International Permafrost Association (IPA) -"Permafrost Extension during the Last Permafrost Maximum (LPM) in the Northern Hemisphere", in particular -to compile an evidence-based map of maximum permafrost extent during the last glaciations period (J.Vandenberghe, personal communication, 2013).
The issues addressed in this study are the following: -How good are explicitly computed cold-region, nearsurface ground states in the PMIP3 models that employ newer generation GCMs/ESMs, with improved physics and finer spatial resolutions?And how do modeled ground temperature distributions compare with previous knowledge from the evidence-based maps?
-How do frozen ground distributions based on air temperature (indirect method) and on ground temperature (direct method) compare?
-How well do the derived distributions reconstruct permafrost boundary locations (especially in the Northern Hemisphere)?
This paper provides a descriptive analysis of the multi-model PMIP3 ensemble, with respect to distribution of modeled surface and subsurface thermal states, and the geographical extent of reconstructed near-surface frozen ground.Progress in modeling frozen ground characteristics after PMIP2 simulations is assessed through comparisons with the evidencebased knowledge and maps currently available.

Experimental design and boundary conditions
Simulations by global paleoclimate models that contributed to PMIP2 (Braconnot et al., 2007) andPMIP3 (PMIP3 2012) were used in this study for the pre-industrial (piControl; 0 kya) and Last Glacial Maximum (LGM; 21 kya) periods.A summary of the used models, institutes or groups, and simulations is shown in Table 1.In PMIP3 simulations, orbital parameters (e.g., eccentricity and obliquity) and trace gas concentrations (such as carbon dioxide, methane, N 2 O, CFC, and ozone) were prescribed to common values for the respective periods of 0 kya and 21 kya (cf.PMIP3 https://wiki.lsce.ipsl.fr/pmip3/doku.php/pmip3:design:21k:final{#}ka experimental design).The pre-industrial control run experiment of PMIP3 is common to and follows design specifications and boundary conditions determined by the CMIP5 experiment protocol (consult Taylor et al., 2012, for details).Ice sheet extent and land/sea distribution, as well as their altitude (topography) at the LGM period, were also prescribed commonly by participating models.The used ice sheet extent was a blend product of ICE-6G v2.0 (Argus and Peltier, 2010;Peltier and Drummond, 2008;Peltier, 2009), MOCA (Tarasov andPeltier, 2002, 2003), and ANU (Lambeck et al., 2010).Detailed information is available at https://wiki.lsce.ipsl.fr/pmip3/doku.php/pmip3:design:pi:final:icesheet.Vegetation was either computed by the model or prescribed to the pre-industrial conditions.We did not stratify this analysis using the difference in vegetation distribution or its mechanisms; these tasks are left for future research.
In PMIP2 experiments, adherence to common boundary conditions was not strict, and, therefore, the orography, land/sea mask, and extent and height of ice sheets varied between the models.Braconnot et al. (2007) offers details regarding the specifications of the participating models.
All boundary conditions and monthly temperature data in PMIP2 and PMIP3 were provided in a netCDF format.The PMIP3 data files were taken from the Earth System Grid Federation (ESGF) system, common to the other CMIP5 data (for example, at http://pcmdi9.llnl.gov/esgf-web-fe/).PMIP2 data files were taken from the PMIP2 database at http:// pmip2.lsce.ipsl.fr/database.The climatology of the variables (i.e., near-surface air and ground temperatures) was computed using the equilibrated thirty years of the simulations.
Since the horizontal resolutions (i.e., the number of grid boxes in the longitude and latitude directions) are different from one model to another (Table 1), spatial interpolation was performed in advance of analysis for a respective common grid system for PMIP2 and PMIP3, by choosing the finest resolution among the models used in the analysis.For PMIP2, this resolution was 128 segments in longitude and 64 segments in latitude (approximately 2.8 degrees by 2.8 degrees resolution), while it was 288 by 192 boxes (approximately 1.25 degrees by 0.94 degrees) for PMIP3.A common mask for the PMIP2 ice sheet and ocean for analysis and display were determined for each grid box by a majority principle among the models used.

Frozen ground zonation
Permafrost is defined as "ground (soil or rock and included ice and organic material) that remains at or below 0 • C for at least two consecutive years" by the IPA (van Everdingen, 1998).According to this definition, permafrost regions can be divided into continuous and discontinuous permafrost zones -in the former, "permafrost occurs everywhere beneath the exposed land surface", or practically more than 90 % of the area underlain by permafrost, while in the latter "permafrost occurs in some areas beneath the exposed land surface throughout a geographic region where other areas are free of permafrost".In the "circum-Arctic map of permafrost and ground-ice conditions" (the IPA map; Brown et al., 1998), however, a part of a discontinuous permafrost zone is further subdivided to sporadic and isolated permafrost zones, in which 10-50 % and 0-10 % of the area is underlain by permafrost, respectively.In this study, we use only two subcategories for permafrost zones, unless otherwise noted: "continuous permafrost" and "discontinuous permafrost".
For the PMIP3 results, soil temperature output (T sl ) at a grid point can determine explicitly the frozen ground type of the grid box between permafrost, seasonally frozen ground, or no freezing.In this study, frozen ground types at a grid point classified from monthly soil temperature include the following three categories, using the thirty-year data examined for climatology (T sl -based criteria): Permafrost(P sl ) : temperature remains at or below 0 • C at any soil layer (1) Seasonally frozen ground(SFG sl ) : the top soil layer(s) freeze and thaw annually, but the layers below remain unfrozen (2) No freezing (NF sl ) : no freezing occurs at any layer (3) However, the determination of frozen ground distribution (i.e., an areal zonation or coverage percentage of the ground thermal state) from a one-point value of the GCM output for an entire area represented by the grid point is not a trivial issue.This will be discussed in Sect. 4. At the time of PMIP2 no soil temperature was archived, so it was not possible to compute the ground thermal regime directly.Saito et al. (2009) developed an indirect methodology for using near-surface air temperature at 2 m (T as ) for estimating the frozen ground type underneath.They employed freezing and thawing indices -originally based on daily values, but computed from the modeled monthly mean outputs (Frauenfeld et al., 2007); thereby a 365 day cycle pro year was assumed for all simulations to enable mutual comparison, although some PMIP2 simulations used a 360-day cycle (Table 1).They derived a four-category classification by comparing the 0.5 • -grid IPA map (Brown et al., 1998) with the 1981-2000 climatology of the Equal-Area 25 km-grid freezing (I f ) and thawing (I t ) indexes interpolated and compiled by Zhang et al. (2005) from the 0.5 • -grid Climate Research Unit (CRU) data (Mitchell and Jones, 2005).Saito et al. (2013b) refined the classification by selecting IPA map subsets of higher certainty with the presence of borehole observations.The resulting T as -based criteria are as follows: Climate-driven permafrost (CP) : No freezing (Nf) : Examination of the T as -based map produced from the observed air freezing and thawing indices showed the following category correspondence: climate-driven permafrost (CP) to "continuous permafrost", environmentally conditional (EP) to "discontinuous permafrost", seasonal frost (Sf) to deep and/or long-lasting seasonal frost, intermittent frost (Im) to short-time frost (i.e., less than two weeks), and Nf to "no freezing" (Saito et al., 2013b).Note that the above classification criteria derived by air indices would be different if derived by surface indices (e.g., snow pack effect).This difference is illustrated in Fig. 2 of Saito et al. (2013b) In the analyses of this study, the statistics are shown by multi-model median (50th percentile value) rather than arithmetic mean, as the number of samples used in the analysis was not large enough to evaluate the Gaussian assumption.

Simulated surface temperature climatology
Climatological differences in the atmospheric thermal conditions at the surface for the LGM, relative to the pre-industrial period, were evaluated for PMIP2 and PMIP3 (Fig. 1).Freezing indices measure the severity of the cold season, while thawing indices denote warm-season characteristics.An increase in the freezing index in the middle to high latitudes was commonly found in the PMIP3 and PMIP2 in the Northern Hemisphere.This increase was enhanced for the PMIP3, reaching up to 7000 • C days along the southern boundary of the continental ice sheets (for the Laurentide and Fennoscandian ice sheets, especially), and up to 4000 • C days in the northern areas across Eurasia and North America, including Beringia (Fig. 1a and b).This increase indicates either additional cooling during the sub-zero season or an extension of the season, or both.Cooling was also noted in the mountain ranges from Caucasus to the Pamir, Tibetan and Ordus plateaus, by up to 2500 • C days along the southern limit of the frozen zones.In the Southern Hemisphere, the highest part of the Andes showed a similar cooling at a magnitude of up to 2000 • C days (not shown).An overall decrease in the global thawing index is found in both experiments except for some warming in northeast Eurasia (Fig. 1c and d).Similar to the increase in the freezing index, this trend indicates either a decrease in the warm-season temperature or a shortening of the season, or both.These results are in line.
MAAT is an amalgamation of the two indices (i.e., MAAT = (I t − I f )/365).The maps of MAAT show clearly that the overall characteristics in MAAT difference are governed by the cold-season trends in both PMIP2 and PMIP3 (Fig. 1e and f).One important implication is that the degrees of cooling at 21 kya between the cold and warm seasons (i.e., changes in I f and I t ) were different from region to region.The geographical distribution of this difference in seasonality, or asymmetry, is essential to interpreting the spatial difference of frozen ground response, resulting in the sometimes insufficient or misleading use of MAAT (Harris, 1981(Harris, , 1982)).Previous experiences have indicated that this is especially applicable in cases of "coarse" spatial scales of 10 000 km 2 or larger, where the effects of other small-scale factors are apt to be smoothed out by spatial averaging (Saito et al., 2009(Saito et al., , 2013b)).

Present-day frozen ground distributions
Reconstructed frozen ground distributions are shown as maps in Fig. 2 for the Northern Hemisphere (NH).The areal extent for the entire NH and its four sub-regions is tabulated in Table 2. Results for 0 kya from the direct method using soil temperature (T sl -criteria; Fig. 2a, and PF3 in Table 2) and the    -71.8 - * Total of "discontinuous", "sporadic", and "isolated" permafrost.
indirect method using surface air temperature (T as -criteria; Fig. 2c, and PG3 in Table 2) show that PMIP3 reconstructions compare well in general to the observation-based evidence for the entire NH, in terms of geographic distribution (Fig. 2g) and areal extent (column 4 in Table 2a).The observed continuous and discontinuous permafrost zones occupy 10.1 % and 12.1 % of the NH land area, respectively; the multi-model results lie in the range of 11.2 % to 15.1 % for the former, and 11.7 % to 11.9 % for the latter (cf.Table 2 of Koven et al., 2013).Further examination of the sub-regions to the north of 40 2c), and North America (170 Table 2d) -also shows good correspondence to observationbased distribution.Presence of permafrost was also reproduced to the south of 40 • N by all diagnostics (Table 2e; extending to low-to mid-latitudes, though the presence of frozen ground is confined to the Tibetan Plateau), despite variations in areal extent and types of permafrost.Underestimation of total permafrost extent in this area is primarily due to smoothed orography in the models, leading to warmer temperatures.The T sl -based map (Fig. 2a) showed dark blue areas (P sl ) that lie between "continuous permafrost" (cold permafrost zone) and "discontinuous" permafrost (warm permafrost zone).Similarly, the SFG sl area (green in Fig. 2a) does not coincide with "seasonally" and "intermittently" frozen ground (green and pale green in Fig. 2g, respectively).This is one example of the gap between a point soil temperature profile and the determination of an areal zonation.These discords with observations have also resulted from the time difference between the pre-industrial period (circa 1850; Taylor et al., 2012) and the observations for the IPA map (late 20th century).This likely owes to a recovery from the Little Ice Age and a long time constant of frozen ground dynamics.In fact, Koven et al. (2013) showed that the modeled permafrost area for 1850 had changed by 2005 by a value ranging from 2 % (increase) to −47 % (decrease) in the CMIP5 multi-model simulations.Another factor, not irrelevant to the previous point, is that the piControl experiment is an equilibrium experiment, while the IPA map includes transient states of the permafrost system.The T as -based method, however, reproduced a result closer to the IPA map, as shown below.Considering the shorter timescale of atmospheric temperature relative to its subsurface counterpart, the discrepancies between the T sl -based map and observations is rather due to this issue of associating point value with a zone of areal extent (including limitations in the classification construction given by Eqs.1-3), and/or implemented subsurface physics.LGM periods, from PMIP3 results.Multimodel mode (the most frequent categories in each grid box) was used.Colors denote permafrost zone (Pf; dark blue), transitional zone (Tr; blue), seasonally freezing zone (Sf; dark green), intermittently freezing zone (If; pale green), no-freezing zone (Nf; orange), and ice sheets (LI; yellow).White solid (dashed) lines denote the southern boundary of continuous (discontinuous) permafrost reconstructed by Vandenberghe et al. (2004Vandenberghe et al. ( , 2012)), while dotted line denotes the reconstruction by Xu et al. (1988).Cross symbols (x) denote the presence of alpine permafrost by Xu et al. (1988).(Brown et al., 1998), and the seasonally freezing ground is determined by mean monthly air temperature at the surface (Saito et al., 2007).(Brown et al., 1998), and the seasonally freezing ground is determined by mean monthly air temperature at the surface (Saito et al., 2007).
The PMIP3 T as -based classification is more closely aligned with observation-based reconstruction.CP and EP are in good correspondence with "continuous" and "discontinuous" permafrost, while Sf and Im are consistent with seasonally and intermittently frozen ground (Fig. 2c and g).The PMIP2 product (Fig. 2e, PF2 in Table 2) is also consistent with observations but failed to represent regional features, particularly in areas with complex relief -e.g., the Altai and Rocky mountains.The finer resolution of PMIP3 definitely contributed to the improvement of T as -based maps with local details.The PMIP3 map shows less zonal features with clearer orographic effects.One of the advantages of T asbased methodology, as compared to the T sl -based method, is its flexibility to downscale.With use of fine-resolution digital elevation models (DEMs) and the assumed lapse rate of atmospheric temperature with height, surface air temperature data from the GCM outputs or observation-derived maps (for example, reanalysis data) can be used to produce a downscaled frozen ground distribution map.This application has been conducted for northeast Asia (Saito et al., 2013b), Japan (Saito et al., 2012), and South America (Saito et al., 2012), using PMIP2 outputs.It is expected that application of PMIP3 data, which produced more realistic climate representation than PMIP2, will yield improved downscaled distribution.
Our frozen ground maps from the PMIP products were consistent with one another (Fig. 2b, d, and f).The more zonally symmetric feature of the LGM permafrost boundary that earlier studies (Velichko and Nechaev, 1984;Vandenberghe et al., 2008Vandenberghe et al., , 2012) ) had also identified was successfully reproduced in the PMIP3 products as well.Previous model simulation studies largely failed to reconstruct permafrost in western Europe due to insufficient winter cooling in the area (Levavassuer et al., 2011;Saito et al., 2009).Vandenberghe et al. (2012) demonstrated that the simulated southern limit of permafrost is in linear relationship with winter North Atlantic sea ice extent.In the PMIP3 results, the T sl -based (direct) method successfully reconstructed the presence of permafrost in western Europe, to the north of the Alps reaching to the Atlantic coast.This is an obvious improvement from the PMIP2 products.However, the T as -based (indirect) method still failed to reconstruct permafrost in the region.This suggests a further need for improvement of the T as classification method.
Differences among these methods are more apparent near the borders between seasonal and non-freezing zones (Fig. 2), and in the Southern Hemisphere (Fig. 1 in the Supplement).However, there is not enough evidence to validate the boundaries between seasonal and non-freezing zones, even for the present day.For the Southern Hemisphere, geomorphological evidence was reported for mountain permafrost in the Andes and Tierra del Fuego during the last glacial period (Trombotto, 2002).The T sl -based estimate produced areas with the presence of LGM permafrost more clearly than the T as -based estimate did, and also appeared closer to the field observations.The PMIP2 product was generally too coarse for mountain permafrost reconstruction, despite some hint of permafrost presence reproduced at the highest area of the Andes.
The statistics on the areal extent of frozen ground regions at 21 kya (column 3) and their changes from 0 kya (column 5) are summarized in Table 2 for the entire Northern Hemisphere and the four sub-regions.Continuous permafrost increased at 21 kya relative to 0 kya in all cases except for North America, which is affected by the substantial decrease of land by the Laurentide ice sheet.The differences in percentage varied between the experiments and methods from 1 to 11 million km 2 (in the case of NH), to which the warm bias of PMIP2 21 kya climate in the European sector greatly contributed.Seasonally frozen ground showed unanimous decrease from the LGM north of 40 • N (Table 2b-d), while it increased in the more equatorial region (Table 2e), reflecting the overall equatorward advancement of frozen area at 21 kya.

Inter-model diversity
Figure 3 shows the variations, or diversity, of the estimated frozen ground types among the models for each method and era, measured by the quantity H = − p i ln(p i ), known as the Shannon-Wiener index, defined by the classification probability p i for the ith category (Shannon, 1948), and plotted after standardization to the range [0, 1].Areas of high diversity (shown in reddish colors) denote regions of dispute among the models, while zero diversity indicates the resulting type is unanimous.For 0 kya, the T sl -based classification showed high diversity between permafrost and seasonally frozen ground, while the T as -based version varied more between seasonal and no frost (Fig. 3a and c).The areas of high diversity shifted southward for 21 kya for both methods, following the overall cooling and the equatorward shift of frozen ground.The values of diversity and areal extent were largely greater for 21 kya, reflecting the general diversity of the modeled atmospheric climate among the models.Correlation analysis on the reconstructed categories showed that the higher diversity found in the PMIP2 map (Fig. 3e and  f), in contrast to the PMIP3 map (Fig. 3c and d), resulted primarily from temperature differences between the models (not shown).The T sl -based methods, rather, showed decreases in the areas of high diversity.

Modeled subsurface thermal regime
Seasonal change in the top soil layer of the frozen ground zone is a direct concern to the eco-climate system, and in modern times to the anthroposphere (including the socioeconomy) as well.It is in the active (seasonally thawed) layer of the permafrost zones that most of the hydrological and biogeochemical activities take place.This layer's thickness defines the volume at which water is held and drained and Table 3. Global areas (in million km 2 ) for maximum active layer thickness (ALT) and maximum seasonal frost depth (SFD).Percentage for ALT relative to the total permafrost regions, and that of SFD relative to the total seasonally frozen ground regions are given in parentheses.

ALT SFD
LGM piControl LGM piControl 0-1 m 7.1 (25.0 %) 0.9 (5.5 %) 9.2 (62.1 %) 20.5 (65.2 %) 1-2 m 12.6 (44.6 %) 9.1 (57.4 %) 4.6 (31.2 %) 9.7 (30.7 %) > 2 m 8.6 (30.5 %) 5.8 (37.1%) 1.0 (6.7%) 1.3 (4.1%) the realm in which the root systems can grow and other biotic matter can take up nutrients.In the seasonally frozen ground zone, the length and depth of the frozen layer restrict surface infiltration and biotic activities, including agriculture.The seasonality of the subsurface thermal regime is closely related to the phenology in the regions.The thickness of the active layer is a measure of change in the permafrost (i.e., aggradation and degradation).The maximum active layer thickness (ALT) is known to have close connection to summertime temperature, as well as local factors such as vegetation cover and previous snow cover, while maximum seasonal frost depth (SFD) is related to winter temperature (and concurrent snow depth).Figure 4 shows ALT (in red) for permafrost regions, and SFD (in blue) in seasonally frozen ground regions, computed from monthly soil temperature profiles.Table 3 summarizes the areal extent of different ALT and SFD for the two eras.The results show a shallower active layer and deeper seasonal frost for the LGM, which is consistent with the changes in freezing and thawing indices (Fig. 1). Figure 5a and b extracted the boundaries for CP and EP, respectively, from Fig. 2 (PF3 in black,FG3 in blue,and FG2 in green).Present-day boundaries taken from observations (IPA map) are shown in red.Present-day results are shown by thin lines, while the LGM results are in thick lines.Modeled boundaries for the LGM are generally in good correspondence with the evidence-based knowledge described above.That the LGM is more "latitudinal" (i.e., parallel to the latitude lines) than the present day is clearly seen, as indicated by Velichko and Nechaev (1984) and Vandenberghe et al. (2008Vandenberghe et al. ( , 2012)).
The first implication of these figures is that the T as -based CP boundary agrees well with the observation-based continuous permafrost boundary, in both Eurasia and North America.Secondly, the direct (T sl -base) method does not necessarily produce a better estimate than the indirect (T as -based) method (e.g., in eastern Siberia).The southern boundary of latitudinal permafrost is determined primarily by the zonal climate, by definition.Vandenberghe et al. (2012) demonstrated the sensitivity of permafrost boundaries to overlying climate, which is impacted by sea-ice distribution on a large scale.However, it is also highly dependent upon local modifying agents, such as snow, soil, vegetation, and topography.The subsurface physics and dynamics implemented in the current generations of the GCMs/ESMs may be quite divergent in the way they handle these factors, as described by Koven et al. (2013), or in terms of total soil column depth (column 7 of Table 1).In fact, the formation, aggradation, maintenance, and degradation of permafrost are controlled not only by atmospheric temperature, but also by other factors such as wetness (precipitation, snow, soil moisture, drainage), soil (the composition of soil with matrix and organic materials, layering, and grain size distribution), vegetation, and micro-topography (aspect, face of slopes, roughness) (Shur and Jorgenson, 2007;Saito et al., 2013a).
Due to this complexity of the subsystem, the subsurface isotherms do not necessarily coincide with the permafrost zone boundaries, similar to the MAAT-permafrost zonation correspondence described in the introduction.; 2007-2009) in the Northern Hemisphere (IPY-TSP; Romanovsky et al., 2010;IPA, 2010).In Eurasia, there are areas of overestimation (cooling biases) between the −5 and −10 • C isotherms in eastern Siberia.In North America, to the contrary, the southwestern part of Alaska has shown underestimation (warm biases).These biases related to the exaggerated southward advance in eastern Siberia of the 0 kya T sl -based permafrost boundaries (thin black line) are shown in Fig. 5a and b.Despite these regional discrepancies, overall correspondence between the modeled pre-industrial isotherms and the observations is surprisingly good, considering possible warming since the "pre-industrial" period (the mid-1800s), as argued in Sect.4.1.In the Asian part of Russia, between 60 • E and 120 • E, permafrost boundaries coincide well with the −5 • C isotherm, while in northern Canada the boundaries are slightly warmer: between 0 • C and − 5 • C isotherms.
Similar to the LGM permafrost boundaries, LGM isotherms showed "latitudinal" characteristics in the middle of Eurasia.To the east, the orographic effects bent the isotherms to the south, while to the west the cooling influence of the ice sheet is apparent.Overall cooling of MAGT at 21 kya relative to 0 kya is about 5 • C in both Eurasia and Alaska.Note that MAGT cooling at LGM may have been greater in some regions, as implied by the large MAAT decrease from the field studies in Europe (below −8 • C in northwestern and central Europe by Huijzer and Vandenberghe, 1998) and China (7-8 • C decrease in northeastern China by Jin et al., 2000;13 • C decrease in Inner Mongolia by Vandenberghe et al., 2004).
The results of this study have demonstrated that direct and indirect methods can be used in combination by compensating mutually.Reconstruction of frozen ground distribution using GCMs/ESMs can be done through three approaches: statistical, physical offline (standalone), and physical online (coupled) approaches.The T as -based method is a statistical approach, while the T sl -based is a physical coupled one.These three approaches all have advantages and disadvantages, but are maybe not mutually exclusive.Simplicity and ease in computation and application (in terms of time, space, and coding) are advantages of the statistical method -though it may lack physical justification when the assumptions or background conditions upon which the relationships were built have changed.In contrast, GCMs/ESMs can compute the subsurface temperature and other variables directly in a physically consistent way.Disadvantages of this method include a high load of computational resources and the general inflexibility to change experimental settings (for example, spatial resolutions and choice of processes to be included in the experiments).Multi-layered or multi-tiered use of these approaches, including offline (one-way) simulations using the physical subsurface models (e.g., Fig. 4 of Saito et al., 2013a, by Marchenko and Romanovsky with their spatially distributed permafrost model GIPL-1.3l),will aid us more effectively to widen and deepen our understanding regarding the roles and contributions of the frozen ground subsystem, and the impact of its change for the global eco-climate system.

Conclusions and implications
Frozen ground distribution at LGM was reconstructed from a multi-experiment analysis using multi-model ensembles from PMIP2 and PMIP3 GCM simulations.Direct (soiltemperature-based) and indirect (surface-air-temperaturebased) estimation approaches were employed, and both showed results more consistent with the present day and glacial periods in comparison with previously known distributions compiled from multiple observational studies.The PMIP3 soil temperature output from the models with finer horizontal resolutions and improved model climatology produced a refined reconstruction with regional details and quantitative information for the subsurface thermal regime for the LGM period.Both the direct and indirect approaches proved to be successful in many aspects, each showing different strength and shortcomings, which also prove somewhat mutually complementary.
The fact that larger inter-model diversity was obtained from the direct method (i.e., based on soil temperature) frozen ground distribution than from the indirect method (based on near-surface air temperature) has implied that the implemented subsurface regime of the GCMs is still at the development phase in comparison to the implemented atmospheric processes, and further suggests the substantial possibility of improvement, in terms of both resolved processes (parameterizations) and specification of subsurface characteristics (settings and parameters; e.g., total soil column depth, thickness of soil layers, choice of appropriate physical property values).Coupling, such as between thermal and hydrological processes and with biogeochemical processes at and below the surface (including snow-vegetation dynamics), is definitely a necessary step in the next generation.Interactions and feedback with vegetation and ecological processes   are very important to understanding the stability and variability of the Arctic terrestrial subsystem in the changing global eco-climate system.To this end, transitional integration, across the late Pleistocene through the Holocene deglaciation to the present day, with differing complexity of implemented mechanisms and resolutions, will provide multi-tiered information regarding different aspects of the frozen ground subsystem in the past and present climate system, and throughout its future projections.
Acknowledgements.This study was supported by the National Science Foundation of the US under ARC-0327664.Portions of this research were conducted under the JAMSTEC-IARC Collaboration Study, with funding provided by the Japan Agency for Marine-Earth Science and Technology (JAMSTEC) under a grant to the International Arctic Research Center (IARC).This paper is also partly supported by the GRENE Arctic Climate Change Research Project.The authors acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups listed in Table 1 of this paper for producing and making available their model output.For CMIP, the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and has led the development of software infrastructure in partnership with the Global Organization for Earth System Science Portals.The authors also acknowledge the international modeling groups for providing their data for analysis, and the Laboratoire des Sciences du Climat et de l'Environnement (LSCE) for collecting and archiving model data.The PMIP2/MOTIF Data Archive is supported by CEA, CNRS, the EU project MOTIF (EVK2-CT-2002-00153) and the Programme National d'Etude de la Dynamique du Climat (PNEDC).The analyses here were performed using version 5 September 2008 of the database.More information is available at http://pmip2.lsce.ipsl.fr/and http://motif.lsce.ipsl.fr/.The authors thank Nate Bauer (International Arctic Research Center, University of Alaska Fairbanks) for language proofreading and editing of the manuscript.
Edited by: G. Ramstein

Fig. 1 .
Fig. 1.Difference in freeze index at the LGM (21 kya) from the pre-industrial period (0 kya), computed from (a) the PMIP3 simulations and (b) PMIP2.(c) and (d) same as (a) and (b) except for the freeze index.(e) and (f) same as (a) and (b) except for mean annual air temperature at the surface (MAAT).

Fig. 2 .
Fig. 2.Frozen ground distribution diagnosed from the Tsl for (a) the pre-industrial and (b) the LGM periods, from PMIP3 results.Multimodel mode (the most frequent categories in each grid box) was used.Colors denote permafrost zone (Pf; dark blue), transitional zone (Tr; blue), seasonally freezing zone (Sf; dark green), intermittently freezing zone (If; pale green), no-freezing zone (Nf; orange), and ice sheets (LI; yellow).White solid (dashed) lines denote the southern boundary of continuous (discontinuous) permafrost reconstructed byVandenberghe et al. (2004Vandenberghe et al. ( , 2012)), while dotted line denotes the reconstruction byXu et al. (1988).Cross symbols (x) denote the presence of alpine permafrost byXu et al. (1988).(c) and (d) same as (a) and (b) except that the diagnoses are T as -based.(e) and (f) same as (c) and (d) except for the PMIP2 results.(g) represents observation-based frozen ground distribution at the present day.Permafrost distribution was taken from the IPA map(Brown et al., 1998), and the seasonally freezing ground is determined by mean monthly air temperature at the surface(Saito et al., 2007).
Fig. 2.Frozen ground distribution diagnosed from the Tsl for (a) the pre-industrial and (b) the LGM periods, from PMIP3 results.Multimodel mode (the most frequent categories in each grid box) was used.Colors denote permafrost zone (Pf; dark blue), transitional zone (Tr; blue), seasonally freezing zone (Sf; dark green), intermittently freezing zone (If; pale green), no-freezing zone (Nf; orange), and ice sheets (LI; yellow).White solid (dashed) lines denote the southern boundary of continuous (discontinuous) permafrost reconstructed byVandenberghe et al. (2004Vandenberghe et al. ( , 2012)), while dotted line denotes the reconstruction byXu et al. (1988).Cross symbols (x) denote the presence of alpine permafrost byXu et al. (1988).(c) and (d) same as (a) and (b) except that the diagnoses are T as -based.(e) and (f) same as (c) and (d) except for the PMIP2 results.(g) represents observation-based frozen ground distribution at the present day.Permafrost distribution was taken from the IPA map(Brown et al., 1998), and the seasonally freezing ground is determined by mean monthly air temperature at the surface(Saito et al., 2007).

Fig. 2 .
Fig. 2. (g) represents observation-based frozen ground distribution at the present day.Permafrost distribution was taken from the IPA map(Brown et al., 1998), and the seasonally freezing ground is determined by mean monthly air temperature at the surface(Saito et al., 2007).

Fig. 3 .
Fig. 3. (a) to (f) same as Fig. 2 except for the inter-model diversity of the diagnosed frozen ground categories.For definition of the diversity index see text.
Figure 5c compares the distribution of the multi-model median of mean annual ground temperature (MAGT) for the permafrost regions for 0 kya (red) and 21 kya (blue).The 0 kya MAGT distribution was compared to observations obtained by the project "Permafrost Observatory Network: a Contribution to the Thermal State of Permafrost (TSP)", conducted during the International Polar Year (IPY

Fig. 4 .
Fig. 4. (a) Multi-model median of maximum active layer depth (in red; for permafrost zones) and seasonal frost depth (in blue; for seasonally frozen ground zones) for the pre-industrial period.(b) same as (a) except for the LGM period. 50

Figure 5
Figure 5. a) Southern boundary of the latitudinal permafrost zones in the Northern 1 Hemisphere for continuous permafrost for pre-industrial (0k; thin lines) and LGM (21k; thick 2 lines) periods, as defined in Figures 2a and 2b, respectively.The boundaries for the T as -based 3 diagnosis (CP) are shown in color (blue for PMIP3 and green for PMIP2), while those for the 4 T sl -based diagnosis (P sl ) are shown in black.Observed boundaries for the present-day 5 continuous permafrost are shown by red lines.b) represents the southern boundary of overall 6 permafrost for 0k (thin lines) and 21k (thick lines).The boundaries for the T as -based diagnosis 7 (colored) are taken as environmentally-conditional (EP), while T sl -based boundaries (black) 8 are defined by areas with seasonal frost depth of less than 1.4 m. c) is the multi-model median 9 of mean annual ground temperature (MAGT) for the pre-industrial (red) and LGM (blue) 10

Fig. 5 .
Fig. 5. (a) Southern boundary of the latitudinal permafrost zones in the Northern Hemisphere for continuous permafrost for pre-industrial (0 kya; thin lines) and LGM (21 kya; thick lines) periods, as defined in Fig. 2a and b, respectively.The boundaries for the T as -based diagnosis (CP) are shown in color (blue for PMIP3 and green for PMIP2), while those for the T sl -based diagnosis (P sl ) are shown in black.Observed boundaries for the present-day continuous permafrost are shown by red lines.(b) represents the southern boundary of overall permafrost for 0 kya (thin lines) and 21 kya (thick lines).The boundaries for the T as -based diagnosis (colored) are taken as environmentally conditional (EP), while T sl -based boundaries (black) are defined by areas with seasonal frost depth of less than 1.4 m.(c) is the multi-model median of mean annual ground temperature (MAGT) for the pre-industrial (red) and LGM (blue) periods.Contour interval is 5 • C.Only temperatures below 0 • C on the exposed land areas are shown.The observed boundaries for the present day are reproduced from (a) and (b) in black thick and thin lines, respectively.Observed MAGT by IPY-TSP are shown by closed circles.

Table 1 .
Summary of model simulations in PMIP2 and PMIP3 used in the study.The horizontal resolution for land is given by the number of grid boxes in the longitudinal and latitudinal directions.Number of soil layers are given for PMIP3 models, followed by the depth of the total soil column.

Table 2 .
Coverage of frozen ground (percentage relative to the land area) in the Northern Hemisphere, simulated by the models in PMIP2 and PMIP3 for the pre-industrial and LGM periods.PF3 denotes T sl -based diagnosis in PMIP3.FG3 and FG2 denote T as -based diagnosis in PMIP3 and PMIP2, respectively.Values are given in median, followed by the minimum and maximum (left and right numbers in the square parenthesis, respectively), and the difference between the two eras.The difference between the two eras is also shown in the rightmost column.Numbers of models used for deriving the statistics for each category are given at the bottom of Table2a.