Study region

The Barrow Peninsula (~1800 km²) is situated on the northern limit of the Arctic Coastal Plain (Fig. 1). The mean annual temperature, precipitation, and snowfall are −11.2 °C, 115 mm, and 958 mm, respectively (1981–2010)³⁵ and the maximum thaw depth ranges from 30 to 90 cm^36,37. This continuous permafrost region is characterized by meso-scale (tens to hundreds of square kilometers) drained thaw lake basins (DTLBs) and interstitial tundra^9,38, which are composed of a mosaic of fine-scale polygonal tundra landforms (tens to hundreds of square meters). Excluding lakes and rivers, the dominant polygonal tundra landforms in this region includes low-center (LC) polygon, flat-center (FC) polygon, high-center (HC) polygon, coalescent LC polygon, drained slopes (DS), nonpatterned DTLB (nDTLB), and thermokarst ponds, which cover an estimated 34, 24, 16, 11, 11, 3, and 1% of the land surface area, respectively^9,13. Due to the similarity in morphological and physiological characteristics of coalescent LC polygons and thermokarst ponds, they are rarely differentiated in field observations. Therefore, both these landforms are combined and referred to as Ponds in the proceeding analysis. Though multiple vegetation communities may be found on each tundra landform, communities typically assemble along a soil moisture gradient representative of each landform²¹. These community–landform associations are identified as follows: dry Salix heath–DS, dry Luzula heath–HC, moist–wet Carex–Oncophorus meadow–FC, moist–wet Carex–Eriophorum meadow–LC, wet Dupontia meadow–nDTLB, and wet Arctophila pond margin–Pond²¹.

Model parameterization and validation

We synthesized an extensive collection of field data measured on the Barrow Peninsula to parameterize and validate DOS-TEM (Supplementary Table 1 and Fig. 3). The majority of this data was acquired by scientific initiatives: (1) International Biological Research Program during the early 1970s^21,38,39,40, (2) Next Generation Ecosystem Experiments between 2010 and 2016^{41,42,43,44,45,46,47}, and (3) Carbon in Arctic Reservoirs Vulnerability Experiment (CARVE) during 2011–2015⁴⁸. In addition, we leveraged key ancillary datasets including: soil carbon pedons (i.e., 100 cm soil cores)^{22,23,24,30,38,39,49,50,51}, vegetation carbon and nitrogen^{21,38,39,40,52}, eddy covariance measurements⁴⁸, and polygonal tundra landform maps^9,13.

Modeled carbon fluxes were compared to net ecosystem exchange (NEE) measurements from the CARVE tower near Utqiaġvik (71°19′22.72′N, 156°35′47.74′W). The tower footprint (~250 m radius) was located in a heterogeneous tundra site composed of all dominant polygonal tundra landforms (exception of Ponds). Although we identified good correspondence with modeled and measured NEE for most of our observations, DOS-TEM underestimated respiratory losses during the zero-curtain seasonal freeze and thaw isothermal period (e.g., September and October)⁵³, resulting in an underestimate of the 1 to 1 line (R² = 0.46, p < 0.001, Fig. 3a, b). Simulated carbon pools for each tundra landform were compared and validated to (i) model benchmarks (Supplementary Table 1) and (ii) an independent subset of soil carbon pools (i.e., pedons), identifying excellent correspondence with modeled and measured carbon pools (Fig. 3c). Together these results demonstrate the ability of DOS-TEM to capture seasonal and inter-annual patterns of carbon dynamics in tundra ecosystems.

Response of soil carbon to climate change

We simulated the response of soil carbon pools to climate change within frozen and seasonally thawed organic (fibric and humic) and mineral horizons. These soil horizons are vertically stratified and vary in the degree of organic matter decomposition. Due to the known data limitations across the Arctic^54,55,56, the complete representation of tundra heterogeneity has not been possible. We used high (i.e., Canadian Centre for Climate Modeling and Analysis (CCCMA) A2) and low (i.e., ECHAM5 B1) climate and emission scenarios, comparable to the five best-performing CMIP5 (Coupled Model Intercomparison Project phase 5) model mean representative concentration pathways 6.0 and 4.5 for the Barrow Peninsula (Supplementary Figs. 1 and 2). Scenarios CCCMA A2 and ECHAM5 B1 project an increase in air temperature (6.96–5.72 °C, precipitation (182–215 mm), and atmospheric CO₂ (509–214 ppm) between 1970 and 2100.

Model simulations indicate that all landforms will gain soil carbon by the end of the twenty-first century (Fig. 4a). However, the trajectories between wetter versus dryer landforms diverged between climate change scenarios. Dependent on climate scenarios from 1970 to 2100, simulated soil C stocks in wet nDTLBs and LC increased between 6127 and 4012 (11.9–7.8%) and between 5246 and 3715 (8.5–6.0%) g C m⁻², respectively (Fig. 4a). Moist FC and aquatic Ponds gained soil carbon at a slightly lower rate than wet landforms, increasing between 2838 and 2073 (5.3–3.9%) and between 4762 and 3970 (7.2–6.0%) g C m⁻², respectively. While dry HC and DS landforms slowly increased in soil carbon content between 2389 and 1895 (3.6–2.8%) and between 2596 and 1834 (3.2–2.3%) g C m⁻², respectively.

The rates of soil carbon accumulation varied by soil horizon (Supplementary Fig. 3). Although all landforms increased soil carbon at a relatively constant rate in the fibric horizon, the rate of accumulation differed between landforms in the humic horizon. Wet landforms continued to increase in carbon content, while moist and dry landforms either did not change (i.e., FC and HC) or lost carbon (i.e., DS). No notable changes in the mineral horizon were identified (Supplementary Fig. 3), with the exception of FC that slowly increased at a rate of ~3.5 g C m⁻² year⁻¹.

We grouped landforms using a single cluster analysis, employing six key biogeophysical characteristics (i.e., water table depth, thaw depth, vegetation carbon, soil carbon, soil nitrogen, and percentage of clay) and re-parameterized and re-calibrated the model to evaluate the uncertainty associated with incrementally reducing the spatial resolution of the representation of tundra heterogeneity. The ensemble of “grouped” landform simulations (Fig. 4b) followed a similar change trajectory as individual landforms (Fig. 4a), identifying the magnitude of change among landforms to be greater than climate uncertainty (CCCMA A2 versus ECHAM5 B1). Similar to individual landform responses, wet landform groups “FC + LC + nDTLB” and “FC + LC + nDTLB + Pond” increased in soil carbon between 4523 and 3397 (7.7–5.8%) and between 3570 and 2570 (6.2–4.4%) g C m⁻², respectively (Fig. 4b). Landform groups moist “FC + LC” and dry “DS + HC” increased in soil carbon content between 3822 and 3126 (6.6–5.4%) and between 1356 and 907 (1.9–1.3%) g C m⁻², respectively (Fig. 4b). The “tundra-biome” landform group (parameterized with data from all landforms) well represented the mean trajectory across landforms, increasing between 2630 and 2094 (4.5–3.6%) g C m⁻² (Fig. 4b). Generally, the change in soil carbon by horizon for landform groups (Fig. 4b) was comparable to that identified by individual landforms (Fig. 4a and Supplementary Fig. 3).

Impact of model spatial scale on soil carbon

DOS-TEM was regionally extrapolated using all landform parameterizations across nine different scales (0.0009, 0.25, 0.5, 1, 2, 4, 8, 16, and 25 km²), to evaluate the influence of model spatial scale on soil carbon dynamics. Input landform distributions were aggregated using a majority filter algorithm (i.e., resampling approach that reclassifies course resolution pixels using the most abundant underlying fine pixel values) from a 30-m resolution (i.e., 0.0009 km²) tundra landform map⁹. Six non-overlapping 75 km² subregions of the Barrow Peninsula (Fig. 5b) were randomly selected to (i) standardize the area extent of the simulation domain and (ii) estimate model uncertainties associated with landform heterogeneity and model spatial scale, expressed as the bias and random error^57,58.

The coarsening of model spatial scales incrementally magnified uncertainties (Fig. 5a) and were directly linked to the misrepresentation of landform distribution (Fig. 6). The bias error became increasingly negative with scale, decreasing by −0.6% for every 1 km² coarsening of spatial scale. Bias errors ranged from 0.5 to 11.9% associated with relatively fine (i.e., ≤4 km²) to coarse spatial scales (i.e., >4 km²), while random error became increasingly positive with scale, increasing by 1.4% for every 1 km² coarsening of spatial scale. Random errors ranged from 3.9 to 22.8% associated with fine to coarse-scale representation of tundra landforms.

Both the bias error and random error were significantly minimized at fine scales (Fig. 5a), as twenty-first century soil carbon was only overestimated by a maximum of 3.7 and ±7.4%, respectively. This is in contrast to coarser spatial scales as bias and random error sharply increased at 8, 16, and 25 km² by −6.1% (±10.7%), −17.0% (±22.1%), and −12.6% (±35.5%), respectively (Fig. 5a). The increase in spatial scale led to the overestimation in the area of low productivity thermokarst lakes (1.1% for every 1 km²) and underestimated wet productive landforms such as Ponds (−0.5% for every 1 km²) and LC polygons (−0.5% for every 1 km²; Fig. 6). This underestimation of wet landforms was particularly concerning as wet landforms have been regionally identified as those most sensitive to change^59,60,61, while representing a significant proportion of the regional carbon cycle^9,60,62,63.

Influence of tundra heterogeneity and model spatial scale

To evaluate the causes, consequences, and mitigation strategies for twenty-first century errors of prediction (i.e., bias and random error), we examined the combined influence of both tundra heterogeneity and model spatial scale. Correlation matrices clarified the potential causes of variable prediction errors, while hierarchical cluster analysis implemented using Euclidean distance and McQuitty linkage methods were used for grouping tundra heterogeneity and model spatial scales with similar errors of prediction to identify potential mitigation strategies or recommendations for future modeling applications.

Correlation matrices supported our presumption that an overestimation of lakes and underestimation of productive wet landforms altered the quantification of landscape-level soil carbon stocks, as bias error was strongly negatively correlated with lake cover (r = −0.98) and positively correlated with wet landforms (r = 0.94; Fig. 7). We found an inverse correlation in bias error as the prevalence of lake cover increased with spatial scale at the expense of nearly all other landforms, but in particular the landforms in low abundance such as tundra ponds (Figs. 6 and 7). Similar to the identified influence of spatial scale on random error (Fig. 4), correlations were highly positively related with model spatial scale (r = 0.99; Fig. 7), reinforcing the impact of coarsening model scale on uncertainty propagation.

Overall, bias error was linked with the misrepresentation of tundra landforms as spatial scale increased (Fig. 7 and Supplementary Fig. 4). Therefore, we next elucidated the influence of heterogeneity and scale on random error. Though random error was correlated with spatial scale, we explored the variability across tundra heterogeneity and scale. The lowest and highest random errors occurred at the finest (≤4 km²) and coarsest (≥16 km²) spatial scales, respectively (Fig. 8). Landform clusters include one or more landforms and landform groups needed to represent tundra heterogeneity on the Barrow Peninsula. Random error was constrained to ±4.5% by considering 5 or 6 tundra landform groups at fine scales. However, at coarse scales these heterogeneous groups also showcased the greatest errors (±28.9%) due to the high number of landforms parameterized within increasingly uncertain landform distributions as scale increased (Figs. 5 and 8). The lowest error among clusters was identified in landform cluster 2 (i.e., ±3.4%; dry and wet), likely due to biogeophysical similarities (i.e., soil anaerobicity, soil available nitrogen, productivity gradients) between dry versus wet landforms (Supplementary Table 1) and similar responses to climate change (e.g., Fig. 4a). Interestingly, even at coarse scales the error found in cluster 2 remained lower than all other landform clusters. Although the “tundra-biome” cluster 1 had a relatively low random error across spatial scales (Fig. 8), this result would not be directly transferable to other modeling applications as we leveraged (i) a robust dataset for model parameterization and (ii) high-resolution polygonal tundra landform maps, currently unavailable across the Arctic for initializing and weighting model parameterization data. The importance of our data assimilation and landform weighting protocol was confirmed by testing the performance of a single unweighted landform parameterization (i.e., HC polygon) extrapolated across the Barrow Peninsula. We found random error to double (±15%) that of cluster 1 at fine scales (≤4 km²) and nearly triple (±45%) at coarse scales (>8 km²). Therefore, to best simulate dynamically changing carbon pools in permafrost soils, our analysis recommends a minimum of two landform groups (i.e., dry and wet) at a maximum model spatial scale of ≤4 km² (Fig. 8).

Implications for modeling soil carbon dynamics in Arctic tundra

Current uncertainties among Pan-Arctic model projections reflect inadequate spatial and temporal data needed to initialize, parameterize, and validate key Arctic ecosystem processes^55,56,64. This study overcame many of these limitations by leveraging a legacy of data (1973–2016) collected from the data-rich Barrow Peninsula to constrain parameter, climate, and model uncertainties, to improve the representation of Arctic tundra heterogeneity across model spatial scales. We identify a scale-dependent balance between tundra heterogeneity and model spatial scale, linked with the decoupling of actual and simulated tundra landform distributions as spatial scales increased (Figs. 5 and 6). The scale-dependency of model process representation is supported by ground-based assessments, as the drivers of carbon dynamics vary across local (e.g., drainage conditions affecting aerobic/anaerobic processes), regional (e.g., vegetation distribution), and landscape scales (e.g., climate variability). Though we identified relatively minimal differences in carbon accumulation rates between polygonal tundra landforms, this was not necessarily surprising as Arctic coastal tundra landforms are relatively young (<5500 years)²⁴, often forming within drained lake basins underlain by the same initial soil substrates⁶⁵. It is not until ice-wedge aggradation alters surface microtopography¹⁴ that local changes in soil moisture, vegetation community composition, and carbon and nitrogen fluxes incrementally alter soil carbon and nitrogen pools. Therefore, our results are ecologically consistent and computationally relevant, as the two recommended dry and wet tundra landforms are at opposite ends of the geomorphological succession spectrum²¹ and found to not only be those most important to adequately represent tundra heterogeneity but also for minimizing prediction errors (bias and random error) while maximizing computational efficiency in Pan-Arctic modeling applications.

Because all models are subject to imprecision associated with imperfect observations⁶⁶, evaluating the influence of various sources of uncertainty in data-poor high-latitude ecosystems are of utmost importance^56,67,68. The inherent randomness in natural systems over space and time was constrained by iteratively synthesizing nearly all possible data combinations to parameterize and simulate change in soil carbon across landforms and groups using relatively moderate climate scenarios for our study region (Supplementary Table 1 and Supplementary Fig. 2). Measurement error (i.e., imprecision of data) was overcome by capitalizing on the robust data collected across multiple projects, initiatives, and programs, thereby minimizing systematic error associated with sampling bias. Due to the difficulty quantifying the influences of cause-and-effect relationships among parameters, we constrained model structural uncertainty in this application by using a single process-based model, ensuring that ecosystem processes are constant among simulations

Despite constraining multiple sources of uncertainty, simulation limitations persist. For example, we did not explicitly simulate carbon dynamics in lakes on the Barrow Peninsula as regional carbon accumulation has been very low throughout the Holocene (~10 g C m⁻² year⁻¹)^69,70,71. In contrast, lakes within yedoma deposits accumulate carbon at a rate of ~47 g C m⁻² year^−1 72. This disparity highlights the importance of not only representing heterogeneous tundra landforms for minimizing uncertainties in Arctic carbon dynamics (Figs. 6 and 8) but also heterogeneous Arctic lakes. Although our results clearly identify the scale-dependency of simulated carbon dynamics, patterns of uncertainty propagation with scale are likely to persist across most spatially explicit model parameters (i.e., available nitrogen), specifically if models are unable to represent sub-grid-scale spatial variability. In addition, similar to many other models⁷³ DOS-TEM does not represent within-grid lateral flow (i.e., adjacent grid cells do not interact), making us unable to evaluate the influence of spatial scale on the lateral transport of carbon, nutrients, and water across the landscape. Recent catchment-scale modeling applications provide a more realistic representation of regional surface and subsurface tundra hydrology and terrestrial–atmosphere fluxes⁷⁴, with the potential to minimize scale-dependent uncertainties. Among the most conspicuous limitations in nearly all terrestrial and earth system models is in the inadequate representation of periglacial landscape evolutionary pathways (i.e., permafrost degradation) that will control spatial and temporal patterns of tundra wetting and drying⁷⁵. Though DOS-TEM partially captures this process associated with the vertically resolved thawing of permafrost, these periglacial ecosystem dynamics are difficult to represent in gridded process models⁷⁶ because (i) data characterizing short- and long-term thermokarst (i.e., subsidence of the surface with permafrost thaw) dynamics are limited, (ii) permafrost degradation does not uniformly initiate across space and time^75,77, and (iii) degradation processes do not occur at a uniform rate for initiated thermokarst^77,78,79. These small but prevalent thermokarst features may represent hot-spots of biogeomorphic change⁸⁰, yet not currently represented in any regional to Pan-Arctic process-based biogeochemistry model.

Due to these known limitations, this assessment focused on improving the representation of heterogeneous Arctic tundra ecosystems in carbon cycle models. Similar to others^4,81,82, simulations indicate that elevated atmospheric CO₂ was the dominant driving factor for projected soil carbon sequestration (sum of carbon from litter, fibric, humic, and mineral soil horizons; Supplementary Fig. 5). This projected increase in soil carbon storage was in line with twenty-first century simulations in northern Alaska⁸³ and the Pan-Arctic⁴. As compared to Pan-Arctic Permafrost Carbon Network model intercomparison⁴, TEM models consistently project higher rates of soil carbon sequestration than other models that do not include nitrogen dynamics, as net primary production is enhanced by nitrogen availability with soil warming. Compared to simulations from a previous version of TEM that used a single parameterization to represent Arctic tundra wetlands of northern Alaska⁸³, our soil carbon accumulation rates were still overestimated by as much as 75.4%. If these tundra wetland simulations⁸³ were implemented with at least two parameterizations (i.e., dry and wet), our findings estimate a threefold decrease in the error of prediction. Collectively, these results highlight that a slight increase in the representation of tundra heterogeneity in terrestrial and Earth System biogeophysical models may notably decrease uncertainties in projected frozen and seasonally thawed soil carbon dynamics in the Arctic.