Saito et al. Progress in Earth and Planetary Science
https://doi.org/10.1186/s40645-020-00345-z(2020) 7:31

Progress in Earth andPlanetary Science

RESEARCH ARTICLE

Open Access

Mapping simulated circum-Arctic organiccarbon, ground ice, and vulnerability of
ice-rich permafrost to degradation
Kazuyuki Saito1* , Hirokazu Machiya1, Go Iwahana2, Hiroshi Ohno3 and Tokuta Yokohata4

Abstract
Permafrost is a large reservoir of soil organic carbon, accounting for about half of all the terrestrial storage, almost
equivalent to twice the atmospheric carbon storage. Hence, permafrost degradation under global warming may
induce a release of a substantial amount of additional greenhouse gases, leading to further warming. In addition togradual degradation through heat conduction, the importance of abrupt thawing or erosion of ice-rich permafrost
has recently been recognized. Such ice-rich permafrost has evolved over a long timescale (i.e., tens to hundreds ofthousands of years). Although important, knowledge on the distribution of vulnerability to degradation, i.e., locationand stored amount of ground ice and soil carbon in ice-rich permafrost, is still limited largely due to the scarcity ofaccessible in situ data. Improving the future projections for the Arctic using the Earth System Models will lead to abetter understanding of the current vulnerability distribution, which is a prerequisite for conducting climatic andbiogeochemical assessment that currently constitutes a large source of uncertainty. In this study, present-day
circum-Arctic distributions (north of 50° N) in ground ice and organic soil carbon content are produced by a new
approach to combine a newly developed conceptual carbon-ice balance model, and a downscaling technique withthe topographical and hydrological information derived from a high-resolution digital elevation model (ETOPO1).
The model simulated the evolution of ground ice and carbon for the recent 125 thousand years (from the Last
Interglacial to the present) at 1° resolution. The 0.2° high-resolution circum-Arctic maps of the present-day groundice and soil organic carbon, downscaled from the 1° simulations, were reasonable compared to the observationbased previous maps. These data, together with a map of vulnerability of ice-rich permafrost to degradation served
as initial and boundary condition data for model improvement and the future projection of additional greenhousegas release potentially caused by permafrost degradation.
Keywords: Circum-Arctic, Soil organic carbon, Distribution of carbon and ice in permafrost, Permafrost degradation,
Conceptual modeling, Topographical and hydrological downscaling

Introduction
Permafrost is a large reservoir of soil organic carbon(SOC), accounting for about half of all the terrestrial
storage, almost equivalent to two times of the atmospheric carbon storage (IPCC 2013; Hugelius et al. 2013a,
2013b, 2014). Therefore, its degradation under global* Correspondence: ksaito@jamstec.go.jp
1
Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 3173-25,
Syowamachi, Kanazawaku, Yokohama, Kanagawa 2360001, Japan
Full list of author information is available at the end of the article

warming may induce a release of a substantial amountof additional greenhouse gases, leading to further warming (Schuur et al. 2015). In addition to gradual degradation through heat conduction, the importance of abrupt
thawing or erosion of ice-rich permafrost has recently
been recognized (Schuur et al. 2015; Strauss et al. 2017;
Turetsky et al. 2020). Such ice-rich permafrost has
evolved over a long timescale (i.e., tens to hundreds ofthousands of years. Murton et al. 2015). Since this icerich permafrost contains large amount of organic matter© The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License,
which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give
appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate ifchanges were made. The images or other third party material in this article are included in the article's Creative Commonslicence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commonslicence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain
permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

in the sediments and greenhouse gases (i.e., carbon dioxide and methane) entrapped in massive ground ice, it isimportant to include the effect of additional impacts induced by degradation of the ice-rich permafrost in future projections (Strauss et al. 2017; Plaza et al. 2019).
Furthermore, the degradation of the ice-rich permafrost
causes irreversible changes in local infrastructures andlivelihoods through surface subsidence and/or lateralerosion (AMAP 2011, 2017; Turetsky et al. 2020) andmay act as a tipping point in the Arctic (Lenton 2012).
Although important, knowledge on the location andstored amount of ground ice and soil carbon in ice-richpermafrost is still limited, largely due to the scarcity ofaccessible in situ data (Hugelius et al. 2014; Olefeldt
et al. 2016; Morris et al. 2018). As for the ground ice(ICE) distribution, the circum-Arctic permafrost map
produced by the effort of the International Permafrost Association (IPA) is the only available data (Brown et al.
1998). Understanding of the current geographical distribution of SOC and ground ice as well as their vulnerability
to degradation is necessary for the enhancement of Earth
System Models as the permafrost-carbon process remains
one of the large sources of uncertainty in climatic and biogeochemical assessment and projections of the Arctic(Schuur et al. 2015; IPCC 2013; MacDougall and Knutti
2016; AMAP 2017; Dean et al. 2018b).
In this study, a simple two-box numerical model developed to simulate the long-term evolution of SOC and ICE
and a downscaling technique with topographichydrological information derived from a high-resolution
digital relief model (ETOPO1, Amante and Eakins 2009)were used to produce the present-day distribution mapsof SOC and ICE in the circumarctic region (north of 50°
N). In contrast to the current datasets produced from synchronic aggregations of contemporary distribution of soil
carbon contents or frozen ground, the maps produced inthis study incorporate diachronically local characteristicsof the soil carbon and subsurface hydrology in thecircum-Arctic region. Hence, though not presented in thispaper, we can readily reconstruct the distribution of SOC
and ICE at any period in the recent 125 thousand years.
The outline of the model structure and integration settings as well as the downscaling methodology are described in the next section. Examples of model
simulations and the resulting distribution maps are presented and examined in the “Results and discussions”
section, together with the implications of the findings
and future improvements.

Methods
Model description

The conceptual model (SOC-ICE-v1.0) used in thisstudy consists of two boxes, i.e., the “above-surface” andthe “subsurface” (Fig. 1). The above-surface box

Page 2 of 15

determines the seasonality, frozen ground state of either
continuous or discontinuous permafrost, seasonal freezing or no freezing (see Saito et al. 2014 for detail), andthe amount of carbon passing down to the subsurface
box. These values and states are calculated from the climate forcing data, namely, mean annual air temperature(MAAT, Ta) and the annual total precipitation (Precip,Pr), modified by the boundary conditions such as latitude, longitude, altitude, distance from the closest ocean
body, presence or absence of ice sheet cover, and background carbon dioxide concentration (CO2). Based onthese input values from the above-ground box and climatic conditions, the subsurface box computes thebudget of SOC, soil moisture, and ICE amount throughthe sets of simple equations. The outlines of the methodology relevant to this study are described below while
the details of the model are explained in the description
paper (Saito et al. 2020).
The carbon accumulation/decay dynamics of the SOCICE-v1.0 model largely followed the idea of the bog
growth model (Clymo 1984, 1992; Belyea and Baird
2006; Yu et al. 2008). The soil carbon content at each location was updated annually (Δt = 1 year) according to
the following equations.
ΔSOCn¼ Litter FallðT a ; P r ; CO2 Þ−κ n ∙SOCn−1
Δt

ð1aÞ





ðT a −T 0 ÞcT
Pr þ P0∙ log
LitterFallðT a ; P r ; CO2 Þ ¼ a1 ðCO2 Þ exp −
bT
bP

ð1bÞ

κ n ¼ κ n−1 þκ crit;n −κ n−1∙Δtτ.

κ crit;n ¼ SOCn

LitterFalln

ð2Þ
ð3Þ

SOCn: soil organic carbon content [kg m−2]
LitterFall: annual carbon input from the above-surface
box [kg m−2 year−1]
κn: decomposition rate of the soil carbon content[year−1]κcrit, n: critical equilibrium rate of decomposition[year−1]
τ: relaxation time constant representing the speed toreach the equilibrium [year]
Litterfall was modeled using temperature and precipitation as predictors (e.g., Liu et al. 2004) under
different atmospheric carbon dioxide concentration(CO2), where its parameters in Eq. 1b were determined by regression to the outputs taken from thebiogeochemical models that participated in the
GRENE-TEA
model
intercomparison
project(GTMIP) for circum-Arctic regions (Miyazaki et al.
2015) and described in the model description paper.

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

Page 3 of 15

Fig. 1 Schematic of the conceptual model—SOC-ICE-v1.0

The parameter τ in Eq. 2 is a hypothetical variable
that represents a timescale for the subsurface carbon
dynamics system to reach an equilibrium state at agiven climate condition. τ serves as an indicator forthe model to evaluate the impact of the allogenic(i.e., climatic and topographic-hydrological) control
on the soil carbon dynamics at the location (Belyea
and Baird 2006; Loisel et al. 2017). The value of thishypothetical parameter may reflect, directly or indirectly, various relevant processes indicating the timescale from years to thousands of years such as
turnover of soil and carbon (Perruchoud et al. 1999;
Luo et al. 2019), or ecological secondary succession
after disturbances (Yoshikawa et al. 2002 and Narita
et al. 2015 under permafrost condition, or Svenningand Sandel 2013 under temperate conditions). In thisstudy, we examined four different baseline values forτ in a geometric series with a factor of 5 to cover
the above range, namely τ = 4, 20, 100, and 500
years. Regardless of the assigned baseline value for a
125 thousand-year integration, however, the instantaneous value of τ was changed to 2500 (or 1500)
when the ground was underlain by continuous (ordiscontinuous) permafrost (as diagnosed by the“above-ground” box) in the simulations. This condition represents the suppression of the carbon dynamics,
especially its decomposition phase due to inactivity of microbial activities, under the frozen conditions (Luo et al.
2019). The interpretation and sensitivity of simulated carbon pool for the allogenic control are discussed in the“Results and discussion” section.
The hydrological components were updated by thefollowing equations.

ΔW n¼ γP r −ξW n−1 −φnΔt

ð4ÞICEn ¼ φn þ ICEn−1

ð5Þ



φn ¼ f k nf ðSOCn ; W n ; ICEn Þ; k tn ðSOCn ; W n Þð6Þ

Wn: soil moisture content [mm = 10−3 m3 m−2]
ICEn: ground ice content [mm = 10−3 m3 m−2]
γ: effective infiltration rate after evapotranspirationand surface runoff [-]ξ: base runoff rate [-]φn: amount of water experiencing the phase change[m3 m−2]
k nf : thermal conductivity for frozen soil [W m−1 K−1]
k tn : thermal conductivity for thawed soil [W m−1 K−1]
Although the model does not explicitly resolve thesubterraneous stratigraphy and structures, we assumed
the depth of hydrological column (ubiquitously set to3000 mm in this study). The physical characteristics andproperties such as soil types (sand, silt, and clay), porosity, and the hydrological parameters in Eq. 4 were determined by the preliminary analyses on the Global Land
Data Assimilation System (GLDAS, Rodell et al. 2004)
datasets. ξ depends on soil type (e.g., 0.09 for sand andclay, 0.045 for silt), and γ varies in the range of 0.61 ~0.99 depending on soil type and temperature (freezing
conditions). No lateral connection between adjacentpoints for hydrology exists. Amount of melting or freezing of ice φn was determined by the energy balance between thawing in summer and freezing in winter. The
ICE budget in Eqs. 5 and 6 was supposed to express the
maximum possible amount of ICE accumulated underSaito et al. Progress in Earth and Planetary Science(2020) 7:31

the given climate history. Furthermore, the absolutevalues of simulated ICE do not necessarily correspond
literally to the observed amount in the field. Rather, they
are meant to express relative distribution among different locations or changes in temporal development.
The driving and boundary conditions

The climate forcing data of annual temperature andprecipitation was adopted from the SeaRISE (Sea-level
Response to Ice Sheet Evolution; Bindshadler et al.
2013) project data. The time series, provided in termsof deviation (temperature) from and ratio (precipitation) to the present-day climatology value, are reconstructed from the Greenland Ice Core Project (GRIP)
as shown in Fig. 2. We chose this dataset because it
was the only time series continuous for the 125
thousand-year period running up to the present day(i.e., from the Last Interglacial to the year 1950) available at the time of the model integration. In order toPage 4 of 15

extend and improve the temporal resolutions in the
recent millennium, simulation results from the PMIP3(Braconnot et al. 2012; especially, from the past millennium runs and historical runs) and the reanalysis
data, namely, the University of Delaware reconstruction product (UDel_AirT_Precip, Willmott and Matsuura 2001) and the ERA-Interim reanalysis data (Deeet al. 2011), were used. These two reanalysis datasets
were also used to calculate the present-day climatology in surface air temperature and precipitation.
For the boundary data, the ICE-6G_C ice sheets data
(Argus et al. 2014; Peltier et al. 2015) was used forthe ice sheet distribution, ice thickness, and sea/landmask for the integration period. Since the climate forcing data were produced from a single Greenland ice
core (the SeaRISE project), they were virtually synchronous for all the locations north of 50° N except
for those regions and periods that were covered by
ice sheets or under water.
The model was run for the 125 thousand years forthe circum-Arctic region, at the 1° resolution for thenorth of 50° N with individual four different baseline
values of τ, considering the presence or absence ofice sheets, and expansion or retreat of the coastlines.
Details of the construction of the driving and boundary condition data are also described in the model description paper (Saito et al. 2020). After a series ofpreliminary sensitivity examinations, the initial SOC
and groundwater were ubiquitously determined as 25
kgC m−2 and 500 mm, respectively. The initial ICE
was set to 0 mm, and the decomposition rate was set
to its critical value. The spin-up run was performedfor 5 thousand years under the initial boundary anddriving conditions to set up subsurface carbon andhydrological (including ice) states before the 125
thousand-year integration.
Downscaling method

Fig. 2 125-k time series of the baseline driving forces fortemperature and precipitation. The SeaRISE time series from the
GRIP ice core are provided in terms of a temperature deviation andb precipitation ratio relative to the present-day values (set as 1950)

The 1 Arc-Minute Global Relief Model (ETOPO1,
Amante and Eakins 2009) was used to determine thehigh-resolution topographic-hydrological information(Fig. 3) for downscaling the present-day SOC and ICE
distribution, which were originally simulated by the
SOC-ICE-v1.0 at the 1° low resolution. The
topographic-hydrological information consisted of thealtitude above the sea level for the topographic feature (Fig. 3c), the specific catchment area (SCA; Fig.
3d), and topographic wetness index (TWI; Fig. 3e) forthe hydrological features. SCA is defined as the contributing area (upstream catchment area) per unit
contour length in the unit of meter (= m2 m−1), while
TWI is defined as the ratio of the natural logarithmof SCA to the tangent of surface slope, indicating thedepth to water table (Tarboton 1989). To calculate

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

Page 5 of 15

Fig. 3 Maps of topography classes and topographic-hydrological variables. Circum-Arctic distributions of topography classes used for a SOC andb ICE. Maps for the Alaskan area of c altitudes, d specific catchment area (SCA), and e topographic wetness index (TWI) calculated from a digital
elevation model ETOPO1 and by the Terrain Analysis Using Digital Elevation Models (TauDEM) software developed by David Tarboton, Utah State
University. Details of classification criteria are described in text

these values, the Terrain Analysis Using Digital Elevation Models (TauDEM) software developed by David
Tarboton was used (Tarboton 1989, 1997).
The downscaling was done in the following manner.
Four categories in the topography classes for SOC(Fig. 3a) were designated to examine the topographical and hydrological conditions to soil carbon dynamics as follows. The “upland” is likely drier andwith thin overburden so that the adaptation is faster,while the “lowland” and “riverine” topography are expected to sustain wetness if climatic condition allows,
and likely to have thicker overburden for SOC to be
decomposed more slowly. The “interior” is of anintermediate character in both wetness and sedimentations. These categories were defined by the altitudeand SCA, i.e., upland (altitude is greater than Altup,

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

and SCA is less than SCAth), interior (altitude is between Altlow and Altup, and SCA is less than SCAth),riverine (SCA is greater than SCAth), and lowland(altitude is less than Altlow, and SCA is less than
SCAth). In this study, Altup was set to 800 m, Altlow
to 100 m, and SCAth to 5 m so that the categories
other than “upland” were largely corresponding to theareas classified as the thermokarst landscapes by Olefeldt et al. (2016), that is, wetland, lake, and hillslope.
A different value of τ was assigned to each of the four classes in consideration to the controls on soil carbon dynamics with respect to wetness and topography: upland (τ = 4),interior (τ = 20), riverine (τ = 100), and lowland (τ = 500).
Then, the value of SOC at a high-resolution grid point
was taken from the results of the 1° simulation that was
executed with the value of τ corresponding to the topography class of the point.
ICE content was downscaled in the following way. As
described in the “Methods” section, ICE in Eqs. 5 and 6
were supposed to express the maximum possibleamount of ground ice accumulation under the given climate history at a location. Downscaled ICE is expected
to reflect the local surficial hydrological conditions. We
used the TWI value to express the localized
topographic-hydrological control on re-distribution of
the water to the ground and to model the high.
resolution ground ice content as ICE, where1:2−TWIICE denotes the ground ice content at the parent 1° grid
point, and TWI is the topographic wetness index calculated at the 1-arc minute resolution. The map was thenaggregated to a 0.2° resolution to remove high-frequency
noise and to highlight the landscape-scale characteristics
on the circum-Arctic maps.

Results and discussions
Simulated time series of SOC

Examples of the 125 thousand-year integration resultsby the SOC-ICE-v1.0 model at two neighboring
points in the lower central Alaska with contrastingbackground climate are shown in Fig. 4. The southernpoint (60° N, 147° W; left-hand column) is warmer
and wetter (Ta = 6.0 ° C, P = 1369 mm) and contains
lowland and upland classes. The northern point (62°N, 147° W; right-hand column) is cooler and drier(Ta = − 2.3 ° C, P = 381 mm) and contains interior andriverine classes. The upper two rows of the panels in
Fig. 4 illustrate carbon and water input to the subsurface box, i.e., LitterFall and infiltration (γP), respectively. Although the input values and resulting SOC
show distinctive difference in the level reflecting thedifferences in the background climate between thetwo points, temporal changes in the resulted SOC
outputs show a common feature: the SOC

Page 6 of 15

accumulated during relatively warmer periods (i.e.,interstadial and interglacial periods), while it mostly
decayed during the colder periods (i.e., stadial andglacial periods). It is also found that the recent accumulation started around 14 thousand years ago during
the deglaciation period after the last glaciation (Vittet al. 2000; MacDonald et al. 2006; Yu et al. 2008;
Morris et al. 2018). The behavior of the modeled
SOC evolution can be interpreted as below. The frozen ground condition physically suppresses the decomposition due to cold temperature (i.e., small κvalue). Under cold temperature, carbon input (litterfall) also decreases. The decrease in litterfall is more
rapid than the adapted subsurface carbon dynamics(due to large τ value) because litterfall responds simultaneously to the change in climate conditions (Eq.
1). These responses lead to a limited amount of SOC
accumulation during the cold periods. In the thawing
phase (e.g., post-glaciation and Holocene, 10–20 ka),
the climate favors more carbon input (i.e., ecosystem
activities are enhanced). In a case of smaller τ (e.g., 4years) the decomposition rate κ catches up faster with
the new warm climate and reaches an equilibrium,
while for larger τ (e.g., 500 years), it takes a longer
time for κ to increase and catch up to a new warmer
state. Therefore, decomposition is slower compared tothe increased carbon input, and the accumulation ofcarbon is more apparent for the case of the larger τ.
These results are consistent with the current knowledge and core-based reconstruction of carbon accumulation (Vitt et al. 2000; MacDonald et al. 2006; Yu
et al. 2008). The scatter between the curves in different colors in SOC in Fig. 4 manifests the impact of τon the SOC evolution.
The value of τ is independent of the (autogenic)
microbial and/or internal conditions, but it reflects
the allogenic factors on the determination of the carbon accumulation dynamics (Clymo 1984, 1992; Boudreau and Ruddick 1991; Belyea and Baird 2006;
Loisel et al. 2017). Interpretation suggests that the
combination of the value of the τ and the litterfall
function determines the eco-climatic conditions of thelocation. For example, equilibrium will be reached in
the order of years for tropical rainforests (except forthe swamps or inundated areas) where the input, i.e.,
growth and replacement of the plants, is large andoutput, i.e., decompositions, is fast, whereas the required time scale can be in the order of 103 years at
cold and wet taiga and tundra regions in high latitudes or altitudes because of slow growth and limited
litterfall and slow decomposition (Harden et al. 1992;
Luo et al. 2019). In another interpretation of the results, the value of τ illustrates the topographichydrological control. The small value corresponds toSaito et al. Progress in Earth and Planetary Science(2020) 7:31

Page 7 of 15

Fig. 4 Example simulation results of litterfall, infiltration, SOC, and ICE for the areas (60° N, 147° W, left column) and (62° N, 147° W, right column)
for the 125 thousand years. The colors in SOC denote the value of τ (years): 5 (red), 20 (orange), 100 (yellow), and 500 (light blue)

dry locations or thin sedimentations owing either tohydrology and soil textures conditions (e.g., fastdrainage of water), climate conditions (e.g., small net
precipitation), and/or topography (e.g., steep slopes totransfer sediments), while large values are for wet or
frozen locations where the changes in the carbon dynamics system are substantially suppressed. Nonetheless, the adequate range and possibly theparameterization scheme of the value of τ as well as
the decomposition rate κ in Eq. 2 should be examinedfurther when the SOC-ICE-v1.0 model is extendedly
applied to the temperate to tropical regions.
Circum-Arctic SOC distribution

Figure 5 shows the circum-Arctic maps of the presentday SOC distribution derived from the SOC-ICE-v1.0
model simulations and observation-based datasets. The
variations in magnitude of the simulated SOC amountand their geographical distribution presented in Fig. 5a–d depict such sensitivity on τ on SOC amount that theSaito et al. Progress in Earth and Planetary Science(2020) 7:31

Page 8 of 15

Fig. 5 Present-day circum-Arctic SOC distribution. The simulated distribution of SOC for the circum-Arctic region north of 50° N at 1° resolution
for a τ = 4, b τ = 20, c τ = 100, and d τ = 500. e The integrated SOC map at 0.2° resolution. Observation-based map for f the NCSCDv1 and g
NCSCDv2 datasets. Difference maps between the integrated map and h the NCSCDv1 and i NCSCDv2 datasets. The common color bar for a to g
is shown under g, and that for h and i is shown under hpresent-day SOC amount decreases with a smaller valueof τ and away from the permafrost regions. We shouldnote that Fig. 5a–d panels include the subsea distribution of SOC, as a result of accumulation during the glacial period when those areas were above ground but

eventually submerged after the deglaciation. The amountand distribution of subsea SOC are currently not known
well.
The integrated map (Fig. 5e) at 0.2° resolution produced by the downscaling procedures described in theSaito et al. Progress in Earth and Planetary Science(2020) 7:31

“Methods” section showed better results than any of thecoarse-resolution individual maps shown in Fig. 5a–d in
comparison with the observation-based maps produced
from the Northern Circumpolar Soil Carbon Datasets(NCSCD). Figure 5f shows the map of the version 1(NCSCDv1. Hugelius et al. 2013a) which covers theupper 100 cm of soil layer, and Fig. 5g for the version 2(NCSCDv2. Hugelius et al. 2013b) which covers theupper 300 cm. Comparison in regional distributions between Fig. 5e–g shows the following features. The higher
SOC distribution in Siberia, northern and western
coastal Alaska, and Canadian High Arctic underlain bypermafrost is well represented. High SOC at the southof Hudson Bay is also well reproduced. These are consistent with the current knowledge and observations(e.g., Beilman et al. 2009; Klein et al. 2013; Hugelius
et al. 2013a, 2013b, 2014; Olefeldt et al. 2016; Morris
et al. 2018; Dean et al. 2018a; Charman et al. 2019).
However, there are some regions that show disparity between the simulated and observed. One such disparity is
the West Siberia area between 60° E and 90° E south of
63° N, and the continental area of Northwest Territories
in Canada that show underestimation, and another is theapparent overestimation west of the Hudson Bay. The
reasons for these discrepancies are discussed together
with the ground ice results in the next section.
The total amount of simulated and observation-based
SOC are summarized in Table 1 and compared for theentire and quartile regions in the circum-Arctic. The
total circum-Arctic SOC for the integrated map was
about 1032 PgC. This value is only 60% of the value
shown in IPCC (2013), but compares well to the more
recent estimates, 1119 PgC for NCSCDv2, or 1035 ± 150
PgC (Hugelius et al. 2014). Except for the 0–90° E quarter, the simulated SOC amount lies between the valuesof the two NCSCD datasets. Difference maps of the integrated map from the NCSCD observation-based maps
are shown in Fig. 5h and i, respectively.
Table 2 summarizes the variations in the simulated
SOC amounts among the different overburden types andcontained ground ice content according to the IPA classifications (Brown et al. 1998. Fig. 6c). Types ofPage 9 of 15

overburden cover (sedimentation) are classified either
thick (lowlands, highlands, and depressions) or thin(mountains, ridges, plateaus, and exposed bedrock), andground ice content are categorized by the volumetricmeasure as high (more than 20% for the thick overburden or more than 10% for the thin overburden), medium(10–20%), and low (0–10%). Observational data
(NCSCD) commonly show that SOC is higher for thickoverburden areas and lower in thin overburden areas is
proportionate to the ground ice content within each
overburden category. This proportional tendency appears also for the simulated SOC for the thick overburden areas, but the contrasting difference between thetwo overburden categories is found not only in the integrated map, but it is also common for different values ofτ. This suggests that the model should incorporate processes regarding sedimentation conditions into the soil
carbon dynamics.

Ground ice evolution and distribution

The bottom panels in Fig. 4 show an example of ICE
evolution at the two locations. Simulated ICE at thenorthern point (right-hand column; interior Alaska) accumulated through the glacial period, c. from 110
through to 15 thousand years ago, and melted away during the Holocene. In fields, massive ground ice (icewedges and pool ice) of the Late Pleistocene origins wasfound near Fairbanks (64.8° N, 147.8° W) in Interior Alaska (Hamilton et al. 1988, Plug 2003, Shur et al. 2004,
Kanevskiy et al. 2008). Determination of the age and rateof ground ice accumulation based on field observations
is still a challenge because some of the ice has melted
away during the past warm periods (including the mostrecent one) and age detection using samples within andoutside the ground ice does not necessarily reflect theage of ice formation.
Ground ice evolution only has negligible dependenceon τ, although it could have exerted an impact throughchanges in thermal conductivity of the ground in association with carbon content with smaller conductivity (cf.
Eq. 6).

Table 1 Total amount of simulated SOC and ground ice in the circum-Arctic region
Total SOC [PgC]:

Circum-Arctic

0–90° E

90–180° E

180–270° E

270–360° E

τ=4

887.1

149.2

409.3

203.9

124.7

τ = 20

992.9

198.0

427.6

229.4

137.9

τ = 100

1215.3

304.6

460.3

285.0

165.3

τ = 500

1456.8

429.2

487.6

345.2

194.9

Integrated

1032.2

234.0

431.4

241.9

124.9

NCSCDv1 (0–100 cm)

489.6

79.9

211.9

125.2

72.6

NCSCDv2 (0–300 cm)

1118.3

158.0

448.9

342.3

169.1

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

Page 10 of 15

Table 2 Comparison of simulated SOC in the circum-Arctic region for different physiographical subsections. Averaged values ineach area are shown with the standard deviation in the parentheses
SOC (kgC m^-2)

Area (10^6 km^2)

NCSCDv1

NCSCDv2

Integrated

τ: 4

τ: 20

τ: 100tau500Thick (High)

2.0

37.6 (13.7)

78.2 (24.4)

41.7 (9.8)

37.9 (14.4)

39.1 (12.6)

41.5 (9.4)

43.6 (7.1)

Thick (Med)

2.5

30.0 (14.5)

73.6 (36.1)

36.1 (13.5)

29.5 (17.8)

32.2 (15.4)

37.6 (10.8)

42.5 (7.1)

Thick (Low)1.8

30.7 (16.1)

76.7 (42.9)

29.7 (14.1)

22.8 (16.0)

25.9 (14.1)

31.8 (10.8)

37.2 (8.1)

Thin (High)

2.7

26.0 (13.3)

55.3 (22.7)

37.3 (13.4)

36.2 (14.4)

36.9 (13.4)

38.7 (11.2)

40.4 (9.9)

Thin (Low)

5.6

21.7 (11.6)

49.8 (21.3)

41.4 (9.9)

40.9 (10.6)

41.3 (9.8)

42.2 (8.5)

43.0 (7.9)

Circum-polar distribution of the simulated ICE at 1°
resolution is shown in Fig. 6a, and the integrated
map at the 0.2° resolution is shown in Fig. 6b. Bothmaps show a monotonous increase in ICE along thelatitudes, starting from 60° N to 65° N up to morethan 16 m near the arctic coasts (Fig. 6b). Evidencesof high ground ice content in Siberia (e.g., Murton
et al. 2015; Strauss et al. 2017), Alaska (e.g., Hamilton
et al. 1988, Plug 2003, Shur et al. 2004, Kanevskiy
et al. 2008, 2011), and Canadian High Arctic Archipelago (e.g., Kokelj and Burn 2013, O’Neill et al.
2019) are available and some of them extend muchgreater than 8 m below the present-day ground surface.
Those thick massive ice contents were often thought to besyngenetic ground ice (French and Shur 2010) and wereformed during the Last Glacial periods (MIS2-4. e.g.,French 2007). Following the categories by French andShur (2010) (cf. O’Neill et al. 2019), the simulated amountof ice (solid-state water) in this model conceptually includes pore ice (frozen moisture in the capillary spaces),
segregated ice (ice lenses), vein (wedge) ice, and pool ice(freezing water from pond or lake). Structures of groundice vary according to many factors, such as sizes and typesof matrix material, history of freezing and thawing, andhydrology. One of the advantages in our conceptual modelcompared to conventional depth-resolving land surface
models is the ease in introducing the structure-free consideration of ice budget.
For the circum-Arctic ground ice distribution, the only
currently available dataset is expressed in terms of categorical classification based on volumetric ice content andtypes of overburden (Brown et al. 1998. Fig. 6c), while thesimulated map (Fig. 6a–b) is shown in terms of the thickness of ground ice. Comprehensive data on the amountand distribution of ground ice contents is currently unavailable for comparison with the results. Use of remotesensing techniques, airborne or satellite, may add new information in the future. However, this will require a substantial breakthrough especially where immobile parts ofground ice underlie (Muskett and Romanovsky 2011; Iwahana et al. 2016; Gido et al. 2019).
The integrated map (Fig. 6b) exhibits regional reproducibility as it shows reasonable correspondencewith respect to the latitudinal and altitudinal dependence in terms of ground ice distribution, i.e., itis thicker in higher latitudes and lower altitudes,

Fig. 6 Circum-Arctic maps of ICE distribution. a The simulated ICE at 1° resolution. b The integrated ICE map at 0.2° resolution. c Map of groundice content taken from the IPA map (Brown et al. 1998)

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

especially in the Beringian areas from eastern Siberia to
Alaska where no large ice sheets had covered during thelast glacial period. The simulated map overestimated ICE
content in Eastern to Central Canada, covering most of
the Canadian shield (Wheeler et al. 1996). This area also
exhibited overestimation for SOC (Fig. 5e). Primarily,
these overestimations were likely caused by the disagreement between the limited soil layer over the bedrock in
the area and the current configuration and settings of themodel (i.e., universally common soil thickness of 3000 mm
for subsurface hydrological processes and allowance of
unlimited accumulation of carbon). Secondly, it is possibly
because the translated climate forcing data derived fromthe Greenland ice core may have inadequately reconstructed for this region although the evolution of theLaurentide Ice Sheet, especially its presence until around8 thousand years ago in this region (Dyke 2005), wasalready incorporated in integration. In contrast, the discrepancy does not appear to be substantial in the Scandinavian region probably because of difference in geology,
the earlier retreat of the Fennoscandia Ice Sheet, and its
generally more maritime climate (Stroeven et al. 2016).
On the contrary, the integrated map underestimated
ice in the meridional bands below ca. 65° N, flat lowland in southern West Siberia to Central Asia, and in
the western coastal regions in Alaska. This may bepartly due to the limitations in the employed frozen
ground classification (Saito et al. 2014). It has been
noted that the classification method tends to underestimate the extent of continuous permafrost in thoseregions that are mountainous or have less or little
surface insulation such as snow cover or surface organic layer (Saito and Machiya 2018; Bertran P, personal communications). An effort to improve theclassification is underway (Saito 2019).

Page 11 of 15

Table 3 summarizes the variations in the simulated
ICE amount among the areas of different overburden
and ground ice content as defined by Brown et al.
(1998. Fig. 6c) (similar to Table 2). It also compares
the simulated results with another observation-based
data for Canada (the three right-most columns;
O’Neill et al. 2019). O’Neill et al. produced a new distribution map of ground ice abundance in northern
Canada, using a geographic information system (GIS)enhanced paleogeographical modeling, in which three
types of ground ice of different geocryological characteristics were distinguished (i.e., relict ice, segregatedice, and wedge ice). Relict ice is buried or injectedmassive ice, segregated ice is formed by thin ice
lenses, and wedge ice is refreezing of infilled meltwater in the cracks formed under harsh winter environment. They categorized the abundance into fiveordered classes for each of the ground ice types, i.e.,“high,” “medium,” “low,” “negligible,” and “none,” andassigned numbers from 4 to 0, respectively. The average and standard deviation of the categorical numbers
for each overburden–ground ice content area ofBrown et al., similar to the ICE amount in the left
columns, are presented in the three right columns ofTable 3.
Marked differences were found between the areas ofdifferent overburden types, similar to the SOC results.
For the thick-overburden areas, the amount of simulated ICE follows the gradient of the observational
ground ice content well. In contrast, it failed to differentiate the ice amount for the thin-overburden
areas. This strongly implies that consideration of sedimentation thickness is also needed for downscaling
and localization in addition to the topographichydrological consideration using TWI.

Table 3 Comparison of simulated ground ice in the circum-Arctic and Canadian region for different physiographical subsections.
Averaged values in each area are shown with the standard deviation in the parentheses
ICE

Area
Circumpolar(10^6Integrated τ: 4km^2)

Canada
τ: 20

τ: 100τ: 500

Amount [m^2 m^-3]

CanadaIntegrated τ: 4

τ: 20

τ: 100τ: 500

Amount [m^2 m^-3]

Wedge
Ice

Segregatedice

Relictice[0: none, 1: negligible 2: low,
3: med, 4: high]Thick(High)

2.0

6.3 (5.1)

19.7(12.4)

19.7(12.4)

19.7(12.4)

19.7(12.4)

6.4 (4.2)

23.0
(9.1)

23.0
(9.1)

23.0
(9.1)

23.0
(9.1)

1.9 (0.9) 2.3 (0.8)

1.2
(1.2)

Thick
(Med)

2.5

3.9 (4.5)

13.8(12.5)

13.8(12.5)

13.8(12.5)

13.8(12.5)

5.5 (5.6)

17.2(11.4)

17.2(11.4)

17.2(11.4)

17.1(11.4)

1.1 (0.5) 2.2 (1.0)

0.1
(0.5)

Thick
(Low)

1.8

2.9 (4.5)

9.5(12.8)

9.4(12.8)

9.4(12.8)

9.4(12.8)

1.4 (1.9)

4.9
(5.9)

4.9
(5.9)

4.9
(5.9)

4.9
(5.9)

0.8 (0.4) 1.2 (0.6)

0.0
(0.0)

Thin
(High)

2.7

4.0 (3.4)

18.2(12.3)

18.2(12.3)

18.2(12.3)

18.2(12.3)

6.9 (3.6)

27.6(10.4)

27.6(10.4)

27.6(10.4)

27.6(10.4)

1.9 (1.0) 2.1 (0.9)

0.9
(1.3)

Thin
(Low)

5.6

4.0 (3.3)

18.1(12.3)

18.1(12.3)

18.1(12.3)

18.1(12.3)

5.5 (3.2)

21.0
(9.8)

21.0
(9.8)

20.9
(9.8)

20.9
(9.8)

1.1 (0.6) 1.4 (0.8)

0.4
(0.8)

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

Application and improvements

Under warming conditions, permafrost containing alarge amount of ground ice can induce a large hazard
due to the melting of the massive ice potentially leading
to erosion and subsidence. Similarly, the larger theamount of SOC, the greater is the potential releaseamount of additional GHGs upon degradation of permafrost. The integrated SOC and ICE maps (Figs. 5e and6b) can be combined to produce a distribution map ofpotential hazards due to ice-rich permafrost degradation.
The vulnerability index V(x) of a location x to permafrost degradation was defined as a function of the standardized SOC(x) and ICE(x) values. Standardization was
performed with respect to the respective maximum
values in the entire circum-Arctic domain A, i.e., max fx∈ASOCðxÞg and max fICEðxÞg . Therefore, it is expressedx∈Aby the following equation:

Page 12 of 15

8
<

98
9
=<
=
SOCðxÞ
ICEðxÞ
V ðxÞ ¼∙
: max fSOCðxÞg; : max fICEðxÞg;
x∈A

ð7Þ

x∈A

The overall feature of the resulting vulnerability map(Fig. 7) is largely consistent with the previous studies
showing high vulnerability in riverine and/or tundra areas
in regions such as Eastern Siberia and North Slope, Alaska, where large ice sheets did not cover during the lastglacial periods (Murton et al. 2015; Olefeldt et al. 2016;
Strauss et al. 2017). We observed two marked differencesfrom the previous data: apparent overestimation of vulnerability in Central Canada, west of the Hudson Bay, which
resulted from the overestimation of SOC and ICE as discussed in the previous subsections, and similarly, underestimation of vulnerability in Southern Siberia to Central
Asia and in the western coastal regions in Alaska.
The circum-Arctic distribution maps of SOC, ICE, andthe vulnerability of ice-rich permafrost to degradation inthis study are potentially applied either to the initial or

Fig. 7 High-resolution circum-Arctic map of vulnerability of ice-rich permafrost to degradation. The vulnerability is determined by the values of
SOC and ground ice shown in Figs. 5e and 6b

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

boundary conditions for the Earth System Models forthe present reconstruction and future projection simulations and to evaluate the impacts of additional greenhouse gas release on global warming through various
pathways resulting from permafrost degradation (Yokohata et al. 2020a, 2020b).
This study shows the potential of the new approach inproducing contemporary maps of SOC and ICE distributions. It also indicates needs for improvement in the implemented SOC and ICE dynamics on local and regional
scales. The current model does not resolve the subsurface structure. The thickness, profile and compositionsof materials (e.g., sediments, water, and carbon), andfunctionality of the subterraneous vertical movement
should be considered in future research.
Although the model forcing data were constructedbased on a single ice core record, the simulated timeseries produced reasonable results for the constructionof the present-day distributions. However, climatic evolutions and variations in reality are known to exhibit distinct differences among regions due to various local
conditions such as latitude, topography (e.g., coastline
changes), orography (loss of ice sheets and isostatic rebounds), land surface conditions, and the resulting atmospheric circulations (Bradley 1999; Braconnot et al.
2012). Adjustment and localization of climate forcing
data are also necessary for future improvement.

Conclusions
Integrations from a newly developed simple conceptualmodel for soil organic carbon and ground ice and a downscaling method with topographic-hydrological information
were applied to construct a present-day circum-Arctic distribution of SOC and ICE at a high spatial resolution. The
results compare reasonably well with the previously obtained observation-based knowledge and data. However,improvement in terms of adjustment in the climate forcing
data for localities, re-evaluation of the decomposition rate,
and refinement of the model’s vertical structure for SOC
and ground ice are needed. The newly obtained datasets
are suitable for use as either initial or boundary conditions
for global-scale future projection simulations on the impacts of additional greenhouse gas release due to permafrost degradation. The prepared vulnerability map can serve
as a guiding tool for policymakers and the local (mostly indigenous) population to evaluate the hazard of subsidenceor other types of land devastations and help in mitigating
the impacts of climate change in the high latitudes. Since
the maps were produced by diachronic incorporation of
local evolutions of the soil carbon and subsurface hydrology, similar maps of the circum-Arctic distributions on
SOC, frozen ground, and ICE contents can be reconstructed for any period in the recent 125 thousand years.
Such a reconstruction can provide useful information forPage 13 of 15

the understanding of terrestrial paleo–eco–climate evolutions in the Arctic region.
Abbreviations
AMAP: The Arctic Monitoring and Assessment Programme; GHG: Greenhousegas; ICE: Ground ice; IPA: The International Permafrost Association; IPCC: The
Intergovernmental Panel on Climate Change; MAAT: Mean annual air
temperature; MIS: Marine Isotope Stage; NCSCD: Northern Circumpolar Soil
Carbon Dataset; PMIP3: Paleoclimate Intercomparison Project, Initiative 3;
SCA: Specific catchment area; SeaRISE: Sea-level Response to Ice Sheet
Evolution; SOC: Soil carbon content; TauDEM: The Terrain Analysis Using
Digital Elevation Models; TWI: Topographic wetness index
Acknowledgements
We thank Drs. Atsushi Sato and Hideki Miura for their advice in conducting our
project. We also thank Dr. Jun’ichi Okuno for his comments and analysis ininterpreting the ICE-6G_C data. The UDel_AirT_Precip data was provided by the
NOAA/OAR/ESRL PSD, Boulder, CO, USA, from their Web site at https://www.
esrl.noaa.gov/psd/. We would like to thank Editage (www.editage.com) for the
English language editing.
Authors’ contributions
KS proposed and managed the project and designed both the study and themodel. HM helped model development and carried out the numericalexperiments. GI, HO, and TY helped in the results interpretation and validations.
All authors read and approved the final manuscript.
Funding
This study was conducted as a part of the Environment Research andTechnology Development Fund project (2-1605) “Assessing and Projecting
Greenhouse Gas Release from Large-scale Permafrost Degradation,” supported by the Ministry of Environment and the Environmental Restoration
and Conservation Agency.
Availability of data and materials
The datasets generated during the current study are available from thecorresponding author at reasonable request. The model source codes andrelevant datasets are included in the accompanying model description
article (Saito et al. 2020)
Competing interests
The authors declare that they have no competing interest.
Author details
1
Japan Agency for Marine-Earth Science and Technology (JAMSTEC), 3173-25,
Syowamachi, Kanazawaku, Yokohama, Kanagawa 2360001, Japan.
2International Arctic Research Center, University of Alaska Fairbanks, 2160
Koyukuk Drive, Fairbanks, AK 99775, USA. 3Kitami Institute of Technology, 165
Koencho, Kitami 0908507, Japan. 4National Institute for Environmental
Studies, 16-2 Onogawa, Tsukuba 3050053, Japan.
Received: 27 December 2019 Accepted: 22 June 2020

References
Amante C, Eakins BW (2009) ETOPO1 1 Arc-Minute Global Relief Model:
procedures, data sources and analysis. In: NOAA Technical Memorandum
NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:https://doi.
org/10.7289/V5C8276M. Accessed 23 Dec 2019
AMAP (2011) Snow, water, ice and permafrost in the Arctic (SWIPA): climatechange and the cryosphere. Arctic Monitoring and Assessment Programme(AMAP), Oslo, Norway.
AMAP (2017) Snow, water, ice and permafrost in the Arctic (SWIPA) 2017. Arctic
Monitoring and Assessment Programme (AMAP), Oslo, Norway
Argus DF, Peltier WR, Drummond R, Moore AW (2014) The Antarctica componentof postglacial rebound model ICE-6G_C (VM5a) based upon GPS positioning,
exposure age dating of ice thicknesses, and relative sea level histories.
Geophys J Int 198(1):537–563. https://doi.org/10.1093/gji/ggu140
Beilman DW, MacDonald GM, Smith LC, Reimer PJ (2009) Carbon accumulation inpeatlands of West Siberia over the last 2000 years. Gloval Biogeochemical
Cycles 23:GB1012. https://doi.org/10.1029/2007GB003112

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

Belyea LR, Baird AJ (2006) Beyond “The limits to peat bog growth”: cross-scale
feedback in peatland development. Ecol Monogr 76(3):299–322
Bindshadler RA, Nowicki S, Abe-Ouchi A, Aschwanden A, Choi H, Fastook J,
Granzow G, Greve R, Gutowski G, Herzfeld U, Jackson C, Johnson J, Khroulev
C, Levermann A, Lipscomp WH, Martin MA, Morlighem M, Parizek BR, Pollard
D, Price SF, Ren D, Saito F, Sato T, Seddik H, Seroussi H, Takahashi K, Walker
R, Wang WL (2013) Ice-sheet model sensitivities to environmental forcing
and their use in projecting future sea level (the SeaRISE project). J Glaciol
59(214). https://doi.org/10.3189/2013JoG12J125
Boudreau BP, Ruddick BR (1991) On a reactive continuum representation oforganic matter diagenesis. Am J Sci 291:507–538
Braconnot P, Harrison SP, Kageyama M, Bartlein PJ, Masson-Delmotte V, Abe-Ouchi
A, Otto-Bliesner B, Zhao Y (2012) Evaluation of climate models using
palaeoclimatic data. Nat Clim Chang 2:417–424. https://doi.org/10.1038/
NCLIMATE1456
Bradley RS (1999) Paleoclimatology: reconstructing climates of the Quaternary.
Second edition Academic Press, San Diego
Brown J, Ferrians OJ, Heginbottom JA, Melnikov ES (1998, revised 2001) Circumarctic map of permafrost and ground ice conditions. National Snow and Ice
Data Center, Digital media, Boulder, CO
Charman DJ, Amesbury MJ, Hinchliffe W, Hughes PDM, Mallon G, Blake WH,
Daley TJ, Gallego-Sala AV, Mauquoy D (2019) Drivers of Holocene peatland
carbon accumulation across a climate gradient in northeastern North
America. Quat Sci Rev 121:110–119
Clymo RS (1984) The limits to peat bog growth. Philos Trans R Soc Lond B 303:
605–654
Clymo RS (1992) Models of peat growth. Suo 43:127–136
Dean JF, van derVelde Y, Garnett MH, Dinsmore KJ, Baxter R, Lessels JS, Smith P,
Street LE (2018a) Abundant pre-industrial carbon detected in Canadian Arctic
headwaters: implications for the permafrost carbon feedback. Environ Res
Lett 13:034024
Dean JF, Middelburg JJ, Röckmann T, Aerts R, Blauw LG, Egger M, Jetten MSM, deJong AEE, Meisel OH, Rasigraf O, Slomp CP, in’t Zandt MH, Dolman AJ(2018b) Methane feedbacks to the global climate system in a warmer world.
Rev Geophys, 56:207–250. doi:https://doi.org/10.1002/2017RG000559
Dee DP, Uppala SM, Simmons AJ, Berrisford P, Poli P, Kobayashi S, Andrae U,
Balmaseda MA, Balsamo G, Bauer P, Bechtold P, Beljaars ACM, van de Berg L,
Bidlot J, Bormann Delsol C, Dragani R, Fuentes M, Geer AJ, Haimberger L,
Healy SB, Hersbach H, Hólm EV, Isaksen L, Kållberg P, Köhler M, Matricardi M,
McNally AP, Monge-Sanz BM, Morcrette JJ, Park BK, Peubey C, de Rosnay P,
Tavolato C, Thépaut JN, Vitart F (2011) The ERA-Interim reanalysis:
configuration and performance of the data assimilation system. Q J Roy
Meteor Soc 137:553–597
Dyke AS (2005) Late quaternary vegetation history of northern North America
based on pollen macrofossil, and faunal remains. Géog Phys Quatern 59(2-3):
211–262
French HM (2007) The periglacial environment. John Wiley & Sons Ltd.,
Chichester
French HM, Shur Y (2010) The principles of cryostratigraphy. Earth Sci Rev 101:
190–206
Gido NAA, Bagherbandi M, Sjöberg LE, Tenze R (2019) Studying permafrost byintegrating satellite and in situ data in the northern high-latitude regions.
Acta Geophysica 67:721–734. https://doi.org/10.1007/s11600-019-00276-4
Hamilton TD, Craig JL, Sellmann PV (1988) The Fox permafrost tunnel: a late
Quaternary geologic record in central Alaska. Geological Society of America
Bulletin 100:948–969.
Harden JW, Sundquist ET, Stallard RF, Mark RK (1992) Dynamics of soil carbon
during deglaciation of the Laurentide Ice Sheet. Science. 258(5090):1921–
1924
Hugelius G, Tarnocai C, Broll G, Canadell JG, Kuhry P, Swanson DK (2013a) The
Northern Circumpolar Soil Carbon Database: spatially distributed datasets ofsoil coverage and soil carbon storage in the northern permafrost regions.
Earth Syst Sci Data 5:3–13
Hugelius G, Bockheim JG, Camill P, Elberling B, Grosse G, Harden JW, Johnson K,
Jorgenson T, Koven CD, Kuhry P, Michaelson G, Mishra U, Palmtag J, Ping CL, O'Donnell J, Schirrmeister L, Schuur EAG, Sheng Y, Smith LC, Strauss J, Yu
Z (2013b) A new data set for estimating organic carbon storage to 3 m
depth in soils of the northern circumpolar permafrost region. Earth System
Science Data 5:393–402. https://doi.org/10.5194/essd-5-393-2013
Hugelius G, Strauss J, Zubrzycki S, Harden JW, EAG S, Ping CL, Schirrmeister L,
Grosse G, Michaelson GJ, Koven CD, O’Donnell JA, Elberling B, Mishra U,

Page 14 of 15

Camill P, Yu Z, Palmtag, Kuhry P (2014) Estimated stocks of circumpolar
permafrost carbon with quantified uncertainty ranges and identified data
gaps. Biogeosciences 11(23):6573–6593
IPCC (2013) Climate change 2013: the physical science basis. Contribution ofWorking Group I to the Fifth Assessment Report of the Intergovernmental Panelon Climate Change [Stocker, T.F., D. Qin, G.-K. Plattner, M. Tignor, S.K. Allen, J.
Boschung, A. Nauels, Y. Xia, V. Bex and P.M. Midgley (eds.)]. Cambridge
University Press, Cambridge, United Kingdom and New York, NY, USA
Iwahana G, Uchida M, Liu L, Gong W, Meyer F, Guritz R, Yamanokuchi T, Hinzman
L (2016) InSAR detection and field evidence for thermokarst after a tundra
wildfire, using ALOS-PALSAR. Remote Sens 8(3):218
Kanevskiy M, Fortier D, Shur Y, Bray M, Jorgenson T. (2008) Detailed
cryostratigraphic studies of syngenetic permafrost in the winze of the CRREL
permafrost tunnel, Fox, Alaska. In: Kane DL, Hinkel KM (eds) Proceedings ofthe ninth international conference on permafrost, Fairbanks, 889–894
Kanevskiy M, Shur Y, Fortier D, Jorgenson MT, Stephani E (2011) Cryostratigraphyof late Pleistocene syngenetic permafrost (yedoma) in northern Alaska, Itkillik
River exposure. Quat Res 75:584–596. https://doi.org/10.1016/j.yqres.2010.12.
003
Klein ES, Booth RK, Yu Z, Mark BG, Stansell ND (2013) Hydrologymediated differential response of carbon accumulation to late
Holocene climate change at two peatlands in Southcentral Alaska.
Quat Sci Rev 64:61–75
Kokelj SV, Burn CR (2003) Ground ice and soluble cations in near-surface
permafrost, Inuvik, Northwest Territories, Canada. Permafr Periglac Process 14:
275–289
Lenton TM (2012) Arctic climate tipping points. Ambio 41(1):10–22. https://doi.
org/10.1007/s13280-011-0221-x
Liu C, Westman CJ, Berg B, Kutsch W, Wang GZ, Man R, Ilvesniemi H (2004)
Variation in litterfall-climate relationships between coniferous and broadleaf
forests in Eurasia. Glob Ecol Biogeogr 13(2):105–114. https://doi.org/10.1111/j.
1466-882X.2004.00072.x
Loisel J, van Bellen S, Pelletier L, Talbot J, Hugelius G, Karran D, Yu Z, Nichols J,
Holmquist J (2017) Insights and issues with estimating northern peatland
carbon stocks and fluxes since the Last Glacial Maximum. Earth Sci Rev 165:
59–80
Luo Z, Wang G, Wang E (2019) Global subsoil organic carbon turnover times
dominantly controlled by soil properties rather than climate. Nat Commun
10:3688. https://doi.org/10.1038/s41467-019-11597-9
MacDonald GM, Beilman DW, Kremenetski KV, Sheng Y, Smith LC, Velichko AA(2006) Rapid early development of circumarctic peatlands and atmosphericCH4 and CO2 variations. Science 314:285–288. https://doi.org/10.1126/
science.1131722
MacDougall AH, Knutti R (2016) Projecting the release of carbon frompermafrost soils using a perturbed parameter ensemble modelling
approach. Biogeosciences 13:2123–2136. https://doi.org/10.5194/bg-132123-2016
Miyazaki S, Saito K, Mori J, Yamazaki T, Ise T, Arakida H, Hajima T, Iijima Y, Machiya
H, Sueyoshi T, Yabuki H, Burke EJ, Hosaka M, Ichii K, Ikawa H, Ito A, Kotani A,
Matsuura Y, Niwano M, Nitta T, O’ishi R, Ohta T, Park H, Sasai T, Sato A, Sato
H, Sugimoto A, Suzuki R, Tanaka K, Yamaguchi S, Yoshimura K (2015) The
GRENE-TEA model intercomparison project (GTMIP): overview andexperiment protocol for Stage 1. Geosci Model Dev 8:2841–2856. https://doi.
org/10.5194/gmd-8-2841-2015
Morris PJ, Swindles GT, Valdes PJ, Ivanovic RF, Gregoire LJ, Smith MW, Tarasov L,
Haywood AM, Bacon KL (2018) Global peatland initiation driven by regionallyasynchronous warming. PNAS 115(19):4851–4856. https://doi.org/10.1073/
pnas.1717838115
Murton JB, Goslar T, Edwards ME, Bateman MD, Danilov PP, Savvinov GN, Gubin
SV, Ghaleb B, Haile J, Kanevskiy M, Lozhkin AV, Lupachev AV, Murton DK,
Shur Y, Tikhonov A, Vasil'chuk AC, Vasil'chuk YK, Wolfe SA (2015)
Palaeoenvironmental interpretation of Yedoma silt (Ice Complex) deposition
as cold-climate loess, Duvanny Yar, Northeast Siberia. Permafr Periglac
Process 26(3):208–288. https://doi.org/10.1002/ppp.1843
Muskett RR, Romanovsky VE (2011) Alaskan permafrost groundwater storagechanges derived from grace and ground measurements. Remote Sens 3:
378–397. https://doi.org/10.3390/rs3020378
Narita K, Harada K, Saito K, Sawada Y, Fukuda M, Tsuyuzaki S (2015) Vegetationand permafrost thaw depth 10 years after a tundra fire in 2002, Seward
Peninsula, Alaska. Arct Antarct Alp Res 47(3):547–559. https://doi.org/10.1657/
AAAR0013-031

Saito et al. Progress in Earth and Planetary Science(2020) 7:31

O’Neill HB, Wolfe SA, Duchesne C (2019) New ground ice maps for Canada using
a paleogeographic modelling approach. Cryosphere 13:753–773. https://doi.
org/10.5194/tc-13-753-2019
Olefeldt D, Goswami S, Grosse G, Hayes D, Hugelius G, Kuhry P, McGuire AD,
Romanovsky VE, Sannel ABK, Schuur EAG, Turetsky MR (2016) Circumpolar
distribution and carbon storage of thermokarst landscapes. Nature Comm 7:13043. https://doi.org/10.1038/ncomms13043
Peltier WR, Argus DF, Drummond R (2015) Space geodesy constrains ice-age
terminal deglaciation: the global ICE-6G_C (VM5a) model. Geophys. Res. Solid
Earth 120:450–487. https://doi.org/10.1002/2014JB011176
Perruchoud D, Joos F, Fischlin A, Hajdas I, Bonani G (1999) Evaluating timescalesof carbon turnover in temperate forest soils with radiocarbon data. Glob
Biogeochem Cycles 13(2):555–573
Plaza C, Pegoraro E, Bracho R, Kathryn GC, Crummer G, Hutchings JA, Hicks Pries
CE, Mauritz M, Natali SM, Salmon VG, Schädel C, Webb EE, Schuur EAF (2019)
Direct observation of permafrost degradation and rapid soil carbon loss intundra. Nat Geosci 12:627–631. https://doi.org/10.1038/s41561-019-0387-6
Plug LJ (2003) Ground-ice features and depth of peat across a mire
chronosequence, NW Alaska. In Permafrost, Phillips, Springman & Arenson(eds), Proceedings of the Eighth International Conference on Permafrost,
901-906.
Rodell M, Houser PR, Jambor U, Gottschalck J, Mitchell K, Meng C, Arsenault K,
Cosgrove B, Radakovich J, Bosilovich M, Entin JK, Walker JP, Lohmann D, Toll
D (2004) The global land data assimilation system. Bull Amer Meteor Soc 85:
381–394. https://doi.org/10.1175/BAMS-85-3-381
Saito K (2019) Re-visiting permafrost zonation classification by climate variables –considerations on elevation and snow factors –. Paper presented at the JSSI
& JSSE Joint Conference on Snow and Ice Research, Yamagata University,
Yamagata, 9 September 2019
Saito K, Machiya H (2018) Past 122-thousand-year frozen ground distribution
north of 50°N: reconstructed advance and retreat. Paper presented at the 5th
European Conference on Permafrost, Congress Center Le Majestic, Chamonix,
26 June 2018
Saito K, Marchenko S, Romanovsky V, Hendricks A, Bigelow N, Yoshikawa K, Walsh
J (2014) Evaluation of LPM permafrost distribution in NE Asia reconstructed
and downscaled from GCM simulations. Boreas 43:733–749. https://doi.org/
10.1111/bor.12038
Saito K, Machiya H, Iwahana G, Yokohata T, Ohno H (2020) Conceptual model tosimulate long-term soil organic carbon and ground ice budget withpermafrost and ice sheets (SOC-ICE-v1.0). Geosci. Model Dev. Discuss. doi:
https://doi.org/10.5194/gmd-2020-80
Schuur EAG, Abbott B, Network PC (2015) Climate change and the permafrost
carbon feedback. Nature 520:171–179. https://doi.org/10.1038/nature14338
Shur Y, French HM, Bray MT, Anderson DA (2004) Syngenetic permafrost growth:
cryostratigraphic observations from the CRREL Tunnel near Fairbanks, Alaska.
Permafr Periglac Process 15(4):339–347
Strauss J, Schirrmeister L, Grosse G, Fortier D, Hugelius G, Knoblauch C,
Romanovsky V, Schädel C, Schneider von Deimling T, Schuur EAG, Shmelev
D, Ulrich M, Veremeeva A (2017) Deep Yedoma permafrost: a synthesis ofdepositional characteristics and carbon vulnerability. Earth Sci Rev 172:75–86.
https://doi.org/10.1016/j.earscirev.2017.07.007
Stroeven AP, Hättestrand C, Kleman J, Heyman J, Fabel D, Fredind O, Goodfellow
BW, Harbor JM, Jansen JD, Olsen L, Caffee MW, Fink D, Lundqvist J, Rosqvist
GC, Strömberg B, Jansson KN (2016) Deglaciation of Fennoscandia. Quat Sci
Rev 147:91–121. https://doi.org/10.1016/j.quascirev.2015.09.016
Svenning J-C, Sandel B (2013) Disequilibrium vegetation dynamics under futureclimate change. Am J Bot 100:1266–1286. https://doi.org/10.3732/ajb.1200469
Tarboton DG (1989) Terrain analysis using digital elevation models (TauDEM).
Available via Hydrologi Research Group, Utah State University, http://hydrology.usu.edu/taudem/taudem5/. Accessed 23 Dec 2019
Tarboton DG (1997) A new method for the determination of flow directions andupslope areas in grid digital elevation models. Water Resour Res 33(2):309–319
Turetsky MR, Abbott BW, Jones MC, Anthony KW, Olefeldt D, Schuur EAG, Grosse
G, Kuhry P, Hugelius G, Koven C, Lawrence DM, Gibson C, Sannel ABK,
McGuire AD (2020) Carbon release through abrupt permafrost thaw. Nat
Geosci 13(138):138–143. https://doi.org/10.1038/s41561-019-0526-0
Vitt DH, Halsey LA, Bauer IE, Campbell C (2000) Spatial and temporal trends incarbon storage of peatlands of continental Western Canada through the
Holocene. Can J Earth Sci 37:683–693

Page 15 of 15

Wheeler JO, Hoffman PF, Card KD, Davidson A, Sanford BV, Okulitch AV, Roest WR(1996) Geological map of Canada. Geological Survey of Canada, "A" Series
Map 1860A. doi:10.4095/208175
Willmott CJ, Matsuura K (2001) Terrestrial air temperature and precipitation:monthly and annual time series (1950 - 1999). Available via NOAA/OAR/ESRL
PSD http://climate.geog.udel.edu/~climate/html_pages/README.ghcn_ts2.
html. Accessed 23 Dec 2019
Yokohata T, Saito K, Takata K, Nitta T, Satoh Y, Hajima T, Sueyoshi T, Iwahana G(2020a) Model improvement and Future projection of permafrost processes
in a global climate model. Submitted to the same SPEPS issue of Progress inEarth and Planetary Science.
Yokohata T, Saito K, Ito A, Ohno H, Tanaka K, Hajima T, Iwahana G (2020b) Future
projection of climate change due to permafrost degradation with a simple
numerical scheme. Submitted to the same SPEPS issue of Progress in Earthand Planetary Science.
Yoshikawa K, Bolton WR, Romanovsky VE, Fukuda M. Hinzman LD (2002) Impactsof wildfire on the permafrost in the boreal forests of Interior Alaska. J
Geophys Res 107:8148. doi:https://doi.org/10.1029/2001JD000438
Yu Z, Beilman DW, Jones MC (2008) Sensitivity of northern peatland carbon
dynamics to Holocene climate change. Geophysical Monograph Series 184:
55–69. https://doi.org/10.1029/2008GM000822

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.