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 warming may induce a release of a substantial amount of 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 of thousands of years. Murton et al. 2015). Since this ice-rich permafrost contains large amount of organic matter in the sediments and greenhouse gases (i.e., carbon dioxide and methane) entrapped in massive ground ice, it is important 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 and livelihoods through surface subsidence and/or lateral erosion (AMAP 2011, 2017; Turetsky et al. 2020) and may act as a tipping point in the Arctic (Lenton 2012). Although important, knowledge on the location and stored amount of ground ice and soil carbon in ice-rich permafrost is still limited, largely due to the scarcity of accessible 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 topographic-hydrological information derived from a high-resolution digital relief model (ETOPO1, Amante and Eakins 2009) were used to produce the present-day distribution maps of 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 in this study incorporate diachronically local characteristics of the soil carbon and subsurface hydrology in the circum-Arctic region. Hence, though not presented in this paper, 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 this study consists of two boxes, i.e., the “above-surface” and the “subsurface” (Fig. 1). The above-surface box determines the seasonality, frozen ground state of either continuous or discontinuous permafrost, seasonal freezing or no freezing (see Saito et al. 2014 for detail), and the 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 on these input values from the above-ground box and climatic conditions, the subsurface box computes the budget of SOC, soil moisture, and ICE amount through the 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).

Fig. 1
figure 1

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

Full size image

The carbon accumulation/decay dynamics of the SOC-ICE-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.

$$ \frac{\Delta {\mathrm{SOC}}_n}{\Delta t}=\mathrm{LitterFall}\left({T}_a,{P}_r,{\mathrm{CO}}_2\right)-{\kappa}_n\bullet {\mathrm{SOC}}_{n-1} $$
(1a)
$$ \mathrm{LitterFall}\left({T}_a,{P}_r,{CO}_2\right)={a}_1\left({\mathrm{CO}}_2\right)\ \exp \left(-\frac{{\left({T}_a-{T}_0\right)}^{c_T}}{b_T}\right)\bullet \log \left(\frac{P_r+{P}_0}{b_P}\right) $$
(1b)
$$ {\kappa}_n={\kappa}_{n-1}+\frac{\left({\kappa}_{\mathrm{crit},n}-{\kappa}_{n-1}\right)}{\tau}\bullet \Delta t $$
(2)
$$ {\kappa}_{\mathrm{crit},n}=\raisebox{1ex}{${\mathrm{SOC}}_n$}\!\left/ \!\raisebox{-1ex}{${\mathrm{LitterFall}}_n$}\right. $$
(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 to reach 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 the biogeochemical 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. 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 a given climate condition. τ serves as an indicator for the 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 this hypothetical 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 Svenning and Sandel 2013 under temperate conditions). In this study, 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 (or discontinuous) 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 the following equations.

$$ \frac{\Delta {W}_n}{\Delta t}=\gamma {P}_r-\xi {W}_{n-1}-{\varphi}_n $$
(4)
$$ {\mathrm{ICE}}_n={\varphi}_n+{\mathrm{ICE}}_{n-1} $$
(5)
$$ {\varphi}_n=f\left({k}_n^f\left({\mathrm{SOC}}_n,{W}_n,{\mathrm{ICE}}_n\right),{k}_n^t\left({\mathrm{SOC}}_n,{W}_n\right)\right) $$
(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 evapotranspiration and surface runoff [-]

ξ: base runoff rate [-]

φn: amount of water experiencing the phase change [m3 m−2]

\( {k}_n^f \): thermal conductivity for frozen soil [W m−1 K−1]

\( {k}_n^t \): thermal conductivity for thawed soil [W m−1 K−1]

Although the model does not explicitly resolve the subterraneous stratigraphy and structures, we assumed the depth of hydrological column (ubiquitously set to 3000 mm in this study). The physical characteristics and properties 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 and clay, 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 adjacent points 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 the given climate history. Furthermore, the absolute values 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 and precipitation was adopted from the SeaRISE (Sea-level Response to Ice Sheet Evolution; Bindshadler et al. 2013) project data. The time series, provided in terms of 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 to 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 (Dee et 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 for the ice sheet distribution, ice thickness, and sea/land mask 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.

Fig. 2
figure 2

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

Full size image

The model was run for the 125 thousand years for the circum-Arctic region, at the 1° resolution for the north of 50° N with individual four different baseline values of τ, considering the presence or absence of ice 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 of preliminary 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 performed for 5 thousand years under the initial boundary and driving conditions to set up subsurface carbon and hydrological (including ice) states before the 125 thousand-year integration.

Downscaling method

The 1 Arc-Minute Global Relief Model (ETOPO1, Amante and Eakins 2009) was used to determine the high-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 the altitude above the sea level for the topographic feature (Fig. 3c), the specific catchment area (SCA; Fig. 3d), and topographic wetness index (TWI; Fig. 3e) for the 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 logarithm of SCA to the tangent of surface slope, indicating the depth to water table (Tarboton 1989). To calculate these values, the Terrain Analysis Using Digital Elevation Models (TauDEM) software developed by David Tarboton was used (Tarboton 1989, 1997).

Fig. 3
figure 3

Maps of topography classes and topographic-hydrological variables. Circum-Arctic distributions of topography classes used for a SOC and b 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

Full size image

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 and with 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 an intermediate character in both wetness and sedimentations. These categories were defined by the altitude and SCA, i.e., upland (altitude is greater than Altup, 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 the areas 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 possible amount 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 \( \raisebox{1ex}{$\mathrm{ICE}$}\!\left/ \!\raisebox{-1ex}{${1.2}^{-\mathrm{TWI}}$}\right. \), where ICE 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 then aggregated 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 results by the SOC-ICE-v1.0 model at two neighboring points in the lower central Alaska with contrasting background climate are shown in Fig. 4. The southern point (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 and riverine 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 the differences in the background climate between the two points, temporal changes in the resulted SOC outputs show a common feature: the SOC accumulated during relatively warmer periods (i.e., interstadial and interglacial periods), while it mostly decayed during the colder periods (i.e., stadial and glacial periods). It is also found that the recent accumulation started around 14 thousand years ago during the deglaciation period after the last glaciation (Vitt et 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., 4 years) 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 to the increased carbon input, and the accumulation of carbon 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.

Fig. 4
figure 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)

Full size image

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 the location. For example, equilibrium will be reached in the order of years for tropical rainforests (except for the swamps or inundated areas) where the input, i.e., growth and replacement of the plants, is large and output, 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 topographic-hydrological control. The small value corresponds to dry locations or thin sedimentations owing either to hydrology and soil textures conditions (e.g., fast drainage of water), climate conditions (e.g., small net precipitation), and/or topography (e.g., steep slopes to transfer 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 the parameterization scheme of the value of τ as well as the decomposition rate κ in Eq. 2 should be examined further 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 present-day SOC distribution derived from the SOC-ICE-v1.0 model simulations and observation-based datasets. The variations in magnitude of the simulated SOC amount and their geographical distribution presented in Fig. 5a–d depict such sensitivity on τ on SOC amount that the present-day SOC amount decreases with a smaller value of τ and away from the permafrost regions. We should note 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 amount and distribution of subsea SOC are currently not known well.

Fig. 5
figure 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 h

Full size image

The integrated map (Fig. 5e) at 0.2° resolution produced by the downscaling procedures described in the “Methods” section showed better results than any of the coarse-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 the upper 100 cm of soil layer, and Fig. 5g for the version 2 (NCSCDv2. Hugelius et al. 2013b) which covers the upper 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 by permafrost is well represented. High SOC at the south of 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 the apparent 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 the entire 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 values of 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 1 Total amount of simulated SOC and ground ice in the circum-Arctic region
Full size table

Table 2 summarizes the variations in the simulated SOC amounts among the different overburden types and contained ground ice content according to the IPA classifications (Brown et al. 1998. Fig. 6c). Types of overburden cover (sedimentation) are classified either thick (lowlands, highlands, and depressions) or thin (mountains, ridges, plateaus, and exposed bedrock), and ground ice content are categorized by the volumetric measure 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 thick overburden 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 the two 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.

Table 2 Comparison of simulated SOC in the circum-Arctic region for different physiographical subsections. Averaged values in each area are shown with the standard deviation in the parentheses
Full size table
Fig. 6
figure 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 ground ice content taken from the IPA map (Brown et al. 1998)

Full size image

Ground ice evolution and distribution

The bottom panels in Fig. 4 show an example of ICE evolution at the two locations. Simulated ICE at the northern 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 (ice wedges and pool ice) of the Late Pleistocene origins was found 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 rate of 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 most recent one) and age detection using samples within and outside the ground ice does not necessarily reflect the age of ice formation.

Ground ice evolution only has negligible dependence on τ, although it could have exerted an impact through changes in thermal conductivity of the ground in association with carbon content with smaller conductivity (cf. Eq. 6).

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. Both maps show a monotonous increase in ICE along the latitudes, starting from 60° N to 65° N up to more than 16 m near the arctic coasts (Fig. 6b). Evidences of 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 much greater than 8 m below the present-day ground surface. Those thick massive ice contents were often thought to be syngenetic ground ice (French and Shur 2010) and were formed during the Last Glacial periods (MIS2-4. e.g., French 2007). Following the categories by French and Shur (2010) (cf. O’Neill et al. 2019), the simulated amount of 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 ground ice vary according to many factors, such as sizes and types of matrix material, history of freezing and thawing, and hydrology. One of the advantages in our conceptual model compared 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 and types of overburden (Brown et al. 1998. Fig. 6c), while the simulated map (Fig. 6a–b) is shown in terms of the thickness of ground ice. Comprehensive data on the amount and distribution of ground ice contents is currently unavailable for comparison with the results. Use of remote sensing techniques, airborne or satellite, may add new information in the future. However, this will require a substantial breakthrough especially where immobile parts of ground 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 correspondence with respect to the latitudinal and altitudinal dependence in terms of ground ice distribution, i.e., it is thicker in higher latitudes and lower altitudes, especially in the Beringian areas from eastern Siberia to Alaska where no large ice sheets had covered during the last 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 the model (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 from the Greenland ice core may have inadequately reconstructed for this region although the evolution of the Laurentide Ice Sheet, especially its presence until around 8 thousand years ago in this region (Dyke 2005), was already 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 be partly 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 those regions 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 the classification is underway (Saito 2019).

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, segregated ice, and wedge ice). Relict ice is buried or injected massive 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 five ordered classes for each of the ground ice types, i.e., “high,” “medium,” “low,” “negligible,” and “none,” and assigned numbers from 4 to 0, respectively. The average and standard deviation of the categorical numbers for each overburden–ground ice content area of Brown et al., similar to the ICE amount in the left columns, are presented in the three right columns of Table 3.

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
Full size table

Marked differences were found between the areas of different 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 topographic-hydrological consideration using TWI.

Application and improvements

Under warming conditions, permafrost containing a large 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 the amount of SOC, the greater is the potential release amount of additional GHGs upon degradation of permafrost. The integrated SOC and ICE maps (Figs. 5e and 6b) can be combined to produce a distribution map of potential 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 \( \mathcal{A} \), i.e., \( \underset{\boldsymbol{x}\boldsymbol{\in}\mathcal{A}}{\max}\left\{\mathrm{SOC}\left(\boldsymbol{x}\right)\right\} \) and \( \underset{\boldsymbol{x}\boldsymbol{\in}\mathcal{A}}{\max}\left\{\mathrm{ICE}\left(\boldsymbol{x}\right)\right\} \). Therefore, it is expressed by the following equation:

$$ V\left(\boldsymbol{x}\right)=\left\{\frac{\mathrm{SOC}\left(\boldsymbol{x}\right)}{\underset{\boldsymbol{x}\boldsymbol{\in}\mathcal{A}}{\max}\left\{\mathrm{SOC}\left(\boldsymbol{x}\right)\right\}}\right\}\bullet \left\{\frac{\mathrm{ICE}\left(\boldsymbol{x}\right)}{\underset{\boldsymbol{x}\boldsymbol{\in}\mathcal{A}}{\max}\left\{\mathrm{ICE}\left(\boldsymbol{x}\right)\right\}}\right\} $$
(7)

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 last glacial periods (Murton et al. 2015; Olefeldt et al. 2016; Strauss et al. 2017). We observed two marked differences from 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.

Fig. 7
figure 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

Full size image

The circum-Arctic distribution maps of SOC, ICE, and the vulnerability of ice-rich permafrost to degradation in this study are potentially applied either to the initial or boundary conditions for the Earth System Models for the 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 in producing 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 compositions of materials (e.g., sediments, water, and carbon), and functionality of the subterraneous vertical movement should be considered in future research.

Although the model forcing data were constructed based on a single ice core record, the simulated time series produced reasonable results for the construction of 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 conceptual model 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 subsidence or 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 for the understanding of terrestrial paleo–eco–climate evolutions in the Arctic region.