Investigating permafrost carbon dynamics in Alaska with artificial intelligence

Frozen soil and carbon-rich permafrost characterizes nearly 14 million square kilometers of the global terrestrial surface, with total soil organic carbon stock estimates near 1307 ± 170 PgC (Hugelius et al 2014). Across the Circumarctic, quantifying the persistent irregularities and impacts attributed to permafrost degradation remains a scientific challenge. The transitional state of permafrost and spatiotemporal ALT heterogeneity drives abrupt changes emerging from rapid, nonlinear carbon-climate feedback mechanisms. These processes are correlated with several biotic and abiotic factors throughout the tundra and boreal, including tundra shrub encroachment, boreal forest migration, caribou migration patterns, topography, precipitation, solar radiation, land surface temperature, and subsurface hydrologic flow (Lloyd et al 2003, Evans et al 2020, Aguirre et al 2021, Joly et al 2021). Carbon release originating from the permafrost-carbon feedback is a climate change catalyst that amplifies localized warming patterns, disrupts carbon cycle partitioning, and destabilizes land–atmosphere coupling mechanisms and feedback thresholds (Sannel and Kuhry 2011, Koven et al 2015, Schuur et al 2015, 2022).

Abrupt permafrost thaw triggers ground subsidence, thermokarst formation, and the proliferation of wetlands, ponds, and methane emission hotspots near littoral zones (Jorgenson et al 2006, 2013, Walter et al 2007, Klapstein et al 2014, Turetsky et al 2014, Olefeldt et al 2016). Across the northern latitudes, increased warming and precipitation trends continue to foster an exponential increase in permafrost carbon release, i.e. permafrost carbon loss is a function of global mean temperature change: ${\text{p}}{{\text{C}}_{{\text{loss}}}} = f\left( {{\text{GMT}}} \right)$ (Pastick et al 2015). Although most model projections arbitrarily terminate between 2100 and 2500 CE, it is purported that anthropogenic-induced warming—coupled with permafrost degradation and soil carbon release—may trigger a positive carbon-climate feedback and sustain a warmer, inhospitable equilibrium for hundreds of thousands of years (Koven et al 2011, McGuire et al 2018, Steffen et al 2018, Wang et al 2023).

Permafrost degradation remains a substantial contributor to the uncertainty of soil carbon release and may promote a positive feedback on the climate system (Koven et al 2015, Burke et al 2017). Synergistic relationships exist among permafrost degradation, thaw subsidence, and terrestrial carbon cycling dynamics. Despite these revelations, the physical processes governing the rate, magnitude, and extent of permafrost degradation and soil carbon remobilization remain uncertain (Koven et al 2015). The mercurial nature of thaw-induced carbon release presents many challenges to quantifying the extent and variability of this feedback (Miner et al 2021). Therefore, it is necessary to constrain drivers of change in model projections to disentangle and interpret the PCF signal (Lawrence et al 2015, Song et al 2020). Rectifying these response patterns requires the ability to scale, simulate, and project the physical dynamics of moisture and energy at the land–atmosphere interface. However, the mechanics and interactions that currently govern each earth system model vary widely (Randall et al 2007, Li et al 2017). Furthermore, few earth system models simulate and account for the mechanisms and long-term effects of the PCF (e.g. UKESM, CanESM, CESM).

Spatiotemporal limitations and instrument constraints inherent to remote sensing technology present considerable obstacles for extracting information from the subsurface profile and restrict the ability to quantify and/or infer thaw variability with high precision and confidence (Gabarró et al 2023). In the high latitudes, the synthesis of optical remote sensing products remains a challenge due to data availability limitations from mechanical restrictions and natural phenomena, such as frequent cloud cover, short summer periods, and low illumination angles (Esau et al 2023, Gay 2023b). These acquisition challenges, compounded with instrument constraints and limited representation in process-based models, prevent accurate quantification of the PCF (i.e. frost heave-thaw dynamics, subsurface freeze- and break-up events, and nutrient cycling).

Despite these knowledge gaps, modeling and reanalysis products are realistic solutions that adopt simple sinusoidal annual temperature gap-filling algorithms. Additionally, cross-sensor probability density functions and probabilistic modeling approaches are used to compensate for crude spatial resolution (Westermann et al 2011). Benchmarking and coupled model intercomparisons (CMIP6) have reduced the broad swaths of uncertainty. Model frameworks continue to improve as a breadth of observational data provide intelligence for calibration and validation (CALVAL) (Collier et al 2018, Beauchamp et al 2020). Concurrently, multimodal applications with AI and observational data quantify, emulate, resolve, and improve upon the physical representation of the climate system with increased flexibility, efficiency, and precision (Koppe et al 2019, Xu et al 2021). Limitations inherent to tractable, modular-limited frameworks and sensitive, computationally intensive AI approaches present opportunities for unifying data harmonization, information sharing, and bias correction among physics-based and data-driven methods (Cuomo et al 2022, Slater et al 2022, Fernandes et al 2023). In response to these sustained and co-emerging limitations, we capture abrupt and persistent changes in subsurface conditions and disentangle control factors driving the PCF with a process-informed, data-driven formulation.

DL is a collection of machine learning algorithms that integrate artificial neural networks (ANN) in multiple layers to discover and extract high-level features and patterns. Convolutional neural networks (CNN) reduce the computational burden of parameterization during training without sacrificing performance. Alternatively, recurrent neural networks (RNN) are complex nonlinear systems that produce feature maps based on impulse sets at discrete time intervals (You et al 2017, Körner and Rußwurm 2021). The equation below (E1) contextualizes the feedback modulation critical for long short-term memory networks (LSTM) wherein learned representation (i.e. recurrent output) at the current time step combines hidden states and loss iteratively (S6):

$$\begin{equation}{y_{\left( t \right)}} = \emptyset \left( {W_X^T{x_{\left( t \right)}} + W_y^T{y_{\left( {t - 1} \right)}} + b} \right).\end{equation} \tag{ E1 }$$

The ability to handle complex multidimensional earth system data while improving data assimilation methods, capturing non-linear relationships, quantifying uncertainty with probabilistic measures, and enhancing modeling precision, efficiency, and forecasting power present important advantages and motivations for implementing AI into current ESM frameworks. Furthermore, encoding previous hidden states is crucial for preserving the internal dynamics and learned representation of memory. This allows the model to reconfigure and assess future cell states based on performance metrics (Peterson 2002). A classic example is soil moisture, a complex and sensitive covariate governed by abiotics (i.e. climatology, topography, soil hydraulics) functioning as a key explanatory variable responsible for encoding the ecological memory of the earth system (Ogle et al 2015). Ecological memory refers to the impact of previous events and drivers of change on the dynamics and complexities of ecosystem response patterns (Khalighi et al 2022). These elusive, interconnected signals illustrate the mercurial nature of memory and change and demonstrate why capturing these effects is critical for simulating real-world dynamics.

The study region is composed of four regional subdomains delineated with borough and census tracts in Alaska (The North Slope [TNS], Interior Boreal [INT], Seward Peninsula [SEW], Yukon-Kuskokwim Delta [YKD]) and over 13.1 million unique observations across 790 field sites. Figure 1 illustrates the distribution of the flux tower network and field campaigns across Alaska. These subdomains are characterized by zone-defined climatological variability and indicate dynamic subsurface conditions and aboveground geomorphology.

Zoom In Zoom Out Reset image size

TNS exhibits continuous permafrost, mesic tundra, acidic/non-acidic loess foothills, hillslope water tracks, colluvial-basin deposits, wet tussock/non-tussock sedges, open/closed erect prostrate and dwarf shrubs, moss and forb communities, and sporadic wet graminoid vegetation at higher altitudes (Walker et al 1989, Raynolds et al 2019). The INT displays extensive land cover heterogeneity, i.e. rolling hills and lowland terrain with coniferous and deciduous forest canopies, continuous and discontinuous frozen soil and permafrost, thermokarst, lowland peat bogs, and wet bottomlands. This region is characterized by continental climate regimes, frequent wildfire activity, high elevation intersections between alpine tundra and arctic lowland tundra, and rapid ecological succession ranging from wet, disturbed regions with shrub and productive graminoid colonization to warm, dry hillsides with spruce replacement. The SEW is characterized by surface/subsurface dynamics ranging from sparse upland tundra flats with low productivity and shallow active layers, well-drained embankments with average productivity and spruce-dominated stands, and densely forested riparian networks with deep active layers (Raynolds et al 2019). The YKD is defined by low-and-upland terrain with extensive wet tussock sedge tundra, lichen-rich peat cover, low vascular plant speciation, and dry tussock graminoid-dwarf shrub communities (Raynolds et al 2019). Plant functional types and species richness contain mesic lichen, sedges, and dwarf shrub vegetation, and isolated instances of erect, prostrate shrubs across volcanic outcrops and drainages (Frost et al 2020). The YKD displays spatially variable, terrain-dependent frozen soil with spatially variable continuous, discontinuous, and sporadic permafrost peppering the upland tundra, flat rolling hills, and eolian surficial deposits. This region exhibits coastal maritime climate patterns with frequent upland wildfires fueled by lichen-rich peat and tussock litter, high shrub abundance, and low seed recruitment (Frost et al 2020).

This research study examines over 13.1 million unique field observations from over 790 sites in Alaska. A considerable proportion of these micrometeorological and biogeochemical data originate from flux towers (AmeriFlux, NEON; appendix C), providing near-continuous measurements (i.e. 30 minute sampling frequency over a 3 to 19 year period). In addition, thaw depth survey measurements were derived via mechanical probing and ground-penetrating radar (GPR) during soil moisture and active layer monitoring campaigns from 1969 to 2022 (e.g. NASA Arctic Boreal Vulnerability Experiment Soil Moisture and Active Layer Thickness (ABoVE SMALT) and Remotely-Sensed Active Layer Thickness (ReSALT (Schaefer et al 2015), Global Terrestrial Network for Permafrost (GTNP), and Circumpolar Active Layer Monitoring program (CALM) (Hinkel and Nelson 2003)). The base dataset consists of 91 variables collected across Alaska over a 53 year period (table 1) and the corresponding Data Sources are provided in appendix C for reference and reproducibility.

Following the workflow in figure 2, the construction of the base dataset and subsequent model development required spatiotemporally constrained preprocessing methods to address latent barriers (i.e. study design), reduce dimensionality, initialize gap-filling and data transformations (i.e. scaling), and preserve information legacies. Dimensionality reduction was used to reduce model complexity and improve generalization, reduce multicollinearity, and enhance computational efficiency. Replicate samples were considered independent features and were not resampled to reduce overgeneralization and bias inflation. The preprocessing pipeline used supervised learning while missing values were imputed by backward-forward gap-filling bilinear interpolation. This gap-filling method served as a padding mask for flux data, reducing sample mismatch and enabling temporal alignment over 53 years. Feature engineering addressed data type conflicts and replicate sampling by removing string objects from the dataframe followed by stationarity tests, i.e., logarithmic transformation and trend differencing (S3).

Zoom In Zoom Out Reset image size

We examined the properties and distribution of the dataset with diagnostics including KDEs and correlograms under various scaling methods to clarify and justify scaling and feature selection (figure 3, table 2). These methods included dimensionality reduction (i.e. Variance Inflation Factorization [VIF]; Principal Component Analysis [PCA]) and spatial autocorrelation metrics (i.e. Wiener–Khinchin theorem; Augmented Dicky Fuller Test [ADF]; Kwiatkowski-Phillips-Schmidt-Shin Test [KPSS]). VIF functions as a diagnostic tool in multivariate forecasting to quantify multicollinearity severity, providing a scalar measurement that encapsulates how much the variance of an estimated regression coefficient increases when the predictors are correlated (Craney Surles 2002, Akinwande et al 2015).

Zoom In Zoom Out Reset image size

PCA compresses data in linear systems; however, for nonlinear relationships, multivariate dimensionality reduction with autoencoding networks was used to explore features' spatiotemporal response patterns via neural-based decomposition. ADF and KPSS methods address deterministic trends and lagging complexities by examining unit root stationarity and non-stationarity in time series data (Lopez 1997). The most significant principal components exhibiting maximum variance were extricated from the 56-feature VIF-reduced multivariate dataset in accordance with dimensionality reduction methods at an 85% cumulative explained variance (CEV) threshold (86.12% [VIF: 61.96%; PCA: 52.17%]). Moreover, pairwise correlation matrices and serial autocorrelations were computed with ranking criteria and Fourier transformations of the power spectra. Based on these results, three PCF features were targeted within the base dataset: ALT, methane flux, and carbon dioxide flux (ALT; CO2_1_1_1; FCH4_1; table 3). These measurements were reframed and prepended with time-lag⁻³ lead⁺² supervision (i.e. learning encouragement) prior to model development.

GeoCryoAI is a state-of-the-art DL architecture composed of bidirectional memory-based recurrent networks (Schuster and Paliwal 1997, Graves and Schmidhuber 2005, Shi et al 2015, Gay et al 2022, Gay 2023b) capable of uncovering hidden response patterns and disturbance signals from in-situ measurements, remote sensing observations, and modeling output (figure 4). This framework presents a novel method to unify spatiotemporal multimodal data. The application quantifies and resolves real-world dynamics and disparities among abrupt and persistent drivers of change. For this study, the in-situ module and corresponding layer architecture is utilized for teacher forcing (S4 and S5). The process-driven component of the network accounts for real-world physical dynamics. The model implements test data and evaluation criteria to simulate ecological memory, highlighting the direct, indirect, concurrent, lagged, internal, and external effects of thaw-induced carbon release on the surrounding environment. The data-driven component is characterized by a CNN-filtered LSTM-encoded RNN that uses teacher forcing from in-situ and flux observations into a network of processing layers with activation and regularization functions (Williams and Zipser 1989). The outputs are trained within an autoencoder framework that encodes and imputes optimized decoding protocol for generative adversarial training and synthetic time series' reconstruction (Yoon et al 2019, Klopries and Schwung 2023). This generative process utilizes encoder-decoder compression to denoise, transform, and learn the latent feature space representation of inputs while using the identity function to reconstruct the inputs from the compressed latent-space vector (Kim et al 2021, Bourland et al 2022, Chen and Guo 2023).

Zoom In Zoom Out Reset image size

The structure of the model configuration is optimized with multimodal inputs and a predefined hyperparameter dictionary generated by a sequential model-based optimization method (i.e. Bayesian optimization search algorithm). These nonlinear programs culminate in a series of logical operations that scan and determine the best-performing model configuration based on regularizing model complexity and optimizing learning efficiency (i.e. early stopping and hyperparameter tuning in latent space). During training, validation loss is used as the objective function. The tuner constructs an a priori distribution model and iterates through a set of hyperparameters (H_p ) to compute the lowest objective score (f(x)) across a distribution of validation metrics (E2). The P(f(x)|H_p ) model is updated every trial with these results, and new criteria are established from a selection function (Wu et al 2019):

$$\begin{equation}{H_p} = \arg \mathop {\min }\limits_{x \in X} f\left( x \right).\end{equation} \tag{ E2 }$$

The GeoCryoAI model emulated hidden dynamics, linking approximately 2.511 million parameters with 1.177 billion interpolated measurements to computationally reflect the state space of the earth system from in-situ measurements. The model design required an efficient, flexible approach to emulate nonlinearities with ground-truth observations with module reconstructions including consolidated tabular time-series layer processing (i.e. LSTM) and sequential time-distributed convolving layers (i.e. Conv1DLSTM). During construction, expanded layer selection was guided by the loss of spatial information with LSTM applications, and the persistent degradation in quality inputs in both data structure and processing resilience. Model performance was optimized via hyperparameter tuning, with the loss function reconfiguring weight matrices while facilitating active learning and convergence, generating output vectors with every batch over each epoch. Furthermore, error metrics were the objective of the NN, with each iteration seeking to reduce the effects of the loss function on weight matrices and vector components. In response to Bayesian optimization hyperparameter tuning and validation loss captured after training, the training and validation subsets were shifted temporally to provide more data for the model to learn (ALT: 1969:2018, 2019–2020; CH₄: 2011:2018, 2019–2020; CO₂: 2003:2018, 2019–2020).

The cleaned 56-feature dataframe yields dynamic relationships; more specifically, thaw depth and subsurface soil respiration demonstrates moderate anticorrelation (r_{soil_[CO2]_15cm}: −0.311942), while thaw depth, CH₄ and CO₂ fluxes illustrate positive correlations (r_{FCH4_1}: 0.334563; r_{FCH4_2}: 0.309447; r_{CO2_1_1_1}: 0.243649). These corroborate the rate-limiting biogeochemical mechanisms that fortify the permafrost-carbon feedback: increased thaw yields higher water availability, halting respiration and instead, fostering warmer anoxic environments suitable for methanogenesis. At the surface, increased thaw and water availability promotes aerobic respiration and evapotranspiration in the soil and surrounding vegetation. Further anticorrelated interactions were noted, i.e. thaw depth relative to GHF (r_{G_1_1_4}: −0.267481) and WT (r_{TW_1_1_1}: −0.219022). In addition to flux variables, a positive correlation was identified between thaw depth and WTD (r_{WTD_2_1_1}: 0.237545). The complexities of the earth system are made apparent wherein mixed relationships were observed as well, with AT exhibiting negative and positive correlations relative to thaw depth (r_{TA_1_2_1}: −0.205626; r_{TA_1_3_1}: −0.180273; r_{TA_1_5_1}: 0.263813). We examined the spatiotemporal distribution (i.e. statewide versus regional, 1969–2022) of permafrost degradation across the four subdomains. TNS exhibits the lowest ALT mean and variance in permafrost degradation (48.460 ± 14.976 cm), INT displays the highest ALT variance (±31.164 cm), and the YKD demonstrates the highest ALT mean (62.415 cm). The temporal period among these features varies substantially (i.e. ALT, 1969–2022; CO₂, 2003–2021; CH₄, 2011–2022). Thaw depth measurements were constrained in alignment with the temporal resolution of carbon flux (CH₄, CO₂) to identify potential nuances governing the dynamics of the PCF. The temporally constrained plots in figure 5 illustrate contrasting annual thaw depth, carbon dioxide, and methane flux variability across Alaska.

Zoom In Zoom Out Reset image size

Abrupt thaw signals in 2014 and 2017 are complemented by persistent increases in CO₂ and CH₄, and increased CH₄, respectively. An interesting anticorrelated relationship is also observed: carbon dioxide and thaw signal relationships are less obvious which may be attributed to increased sequestration of CO₂ resulting in tundra greening and increased primary productivity. Another explanation may be rooted in feedback attribution: carbon dioxide regulates thaw rate, but due to increased warming, methane and thaw experienced runaway effects (i.e. 2014–2017). Error metrics (RMSEs) for ALT, CO₂, and CH₄ were computed from feature variances originating from calibration and measurement error prior to model training (i.e. aleatoric stochasticity): 4.348 cm, 20.117 umolCO₂m⁻² s⁻¹, and 79.185 nmolCH₄m⁻² s⁻¹, respectively.

Initial baseline metrics were explored prior to model construction. The implementation of time-delayed (±1) naive persistence forecasting yielded 1.997 cm, 1.906 µmolCO₂m⁻²s⁻¹, 0.884 nmolCH₄m⁻²s⁻¹ (RMSE). Thereafter, a simple RNN was explored to determine the suitability for emulating the complex interactions characterizing the PCF. The sequential linear regression network was evaluated for comparison, yielding 1.024 cm, 1.114 µmolCO₂m⁻²s⁻¹, 0.295 nmolCH₄m⁻²s⁻¹ (RMSE). These results indicate substantial error reduction by an order of magnitude (RMSE: 1.389) and suggest implementing a more comprehensive and efficient NN that accurately reflects the internal dynamics and ecological memory of the Arctic system. Model simulations were informed in a data-driven sequence of site-level teacher forcing via in-situ and flux observations of ALT (ALT^RMSE: 1.327 cm [1969–2022]; CO₂ ^RMSE: 0.697 µmolCO₂m⁻²s⁻¹ [2003–2021]; CH₄ ^RMSE: 0.715 nmolCH₄m⁻²s⁻¹ [2011–2022]). Evaluation metrics were computed on the training configurations followed by the application of inverse transformations to the test sets prior to feature predictions (figure 6).

Zoom In Zoom Out Reset image size

The evaluation metrics demonstrate the performance and predictive capacity of the model (figure 6(a)). Although the predictions indicate accurate impulses over time (figure 6(b)), the magnitudes of ALT predictions are underestimated while CH₄ and CO₂ flux predictions are overestimated. The contributing factor for these disparities is likely due to data imbalance. Although randomization of data splitting (i.e. training, validation, testing sets) is best practice, data partitioning prioritized sample representation after cleaning while retaining the temporal structure. Moreover, the spatiotemporal sparsity of data samples likely contributed to these disparities, which was both accounted for and mitigated with gap-filling and time lag-lead methods. For example, ALT training data consisted of 72.63% of the 98 753 ALT samples (1969–2022), CO₂ training data consisted of 74.05% of the 547 101 CO₂ samples (2003–2021), and CH₄ training data consisted of 71.23% of the 212 642 CH₄ samples (2011–2022). Another possible factor contributing to underestimation is the introduction of excessive regularization (i.e. dropout) in the modeling framework, leading to an oversimplification of the model and feature variance. Conversely, the potential issues facilitating overestimated feature predictions may originate from model overcomplexity thus capturing spurious dynamics. To counteract these potential discrepancies and misinterpretation of derived results, more sensitive evaluation metrics could be evaluated in addition to further refining of the Bayesian hyperparameter optimization tuning process is necessary to include more operands (i.e. loss functions, activation functions) and continuous arrays (i.e. batch size, epochs) to isolate relevant causes of these magnitude disagreements.

Bakian-Dogaheh K, Chen R H, Moghaddam M, Yi Y and Tabatabaeenejad A 2020 ABoVE: Active Layer Soil Characterization of Permafrost Sites, Northern Alaska, 2018 (ORNL DAAC)

Brown J, Hinkel K M and Nelson F E 2000 The circumpolar active layer monitoring (CALM) program: historical perspectives and initial results Polar Geography 24 165–258

Chen A, Parsekian A, Schaefer K, Jafarov E, Panda S K, Liu L, Zhang T and Zebker H A 2015 Pre-AboVE: Ground-Penetrating Radar Measurements of ALT on the Alaska North Slope (ORNL DAAC)

Clayton L K et al 2021 Active layer thickness as a function of soil water content Environ. Res. Lett. 16 055028

Engram M J, Anthony K W and Meyer F J 2020 AboVE: SAR-based Methane Ebullition Flux from Lakes, Five Regions, Alaska, 2007–2010 (ORNL DAAC)

Nicolsky D J, Romanovsky V E, Kholodov A L, Dolgikh K and Hasson N 2020 ABoVE: Soil Temperature Profiles, USArray Seismic Stations, AK and Canada, 2016–2019 (ORNL DAAC)

Schaefer K et al 2021 Field Measurements of Soil Moisture and Active Layer Thickness in Alaska and Canada (ORNL DAAC)

Schaefer K et al 2021 AboVE: Active Layer Thickness Derived from Airborne L- and P-band SAR, Alaska, 2017 (ORNL DAAC)

Watts J D, Natali S and Minions C 2022 Soil Respiration Maps for the above Domain, 2016–2017 (ORNL DAAC)

Whitley M, Frost G, Jorgenson M T, Macander M, Maio C V and Winder S G 2018 AboVE: Permafrost Measurements and Distribution across the Y-K Delta, Alaska, 2016 (ORNL DAAC)

Yi Y and Kimball J S 2020 AboVE: Active Layer Thickness from Remote Sensing Permafrost Model, Alaska, 2001–2015 (ORNL DAAC)

Rocha A, Shaver G and Hobbie J 2020 AmeriFlux BASE US-An1 Anaktuvuk river severe burn, ver. 2-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246142

Rocha A, Shaver G and Hobbie J 2020 AmeriFlux BASE US-An2 Anaktuvuk river moderate burn, ver. 2-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246143

Rocha A, Shaver G and Hobbie J 2020 AmeriFlux BASE US-An3 Anaktuvuk river unburned, ver. 2-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246144

Billesbach D and Sullivan R 2022 AmeriFlux BASE US-A03 ARM-AMF3-Oliktok, ver. 5-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1498752

Billesbach D and Sullivan R 2021 AmeriFlux BASE US-A10 ARM-NSA-barrow, ver. 4-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1498753

Zona D and Oechel W 2016 AmeriFlux BASE US-Atq Atqasuk, ver. 1-1 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246029

Zona D and Oechel W 2019 AmeriFlux BASE US-Ivo Ivotuk, ver. 4-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246067

Zona D and Oechel W 2016 AmeriFlux BASE US-Brw barrow, ver. 2-1 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246041

Euskirchen E 2022 AmeriFlux BASE US-BZB bonanza creek thermokarst bog, ver. 4-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1773401

Euskirchen E 2022 AmeriFlux BASE US-BZF bonanza creek rich fen, ver. 4-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1756433

Euskirchen E 2022 AmeriFlux BASE US-Bzo bonanza creek old thermokarst bog, ver. 3-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1846662

Euskirchen E 2022 AmeriFlux BASE US-BZS bonanza creek black spruce, ver. 3-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1756434

Euskirchen E, Shaver G and Bret-Harte S 2022 AmeriFlux BASE US-Ich imnavait creek watershed heath tundra, ver. 4-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246133

Euskirchen E, Shaver G and Bret-Harte S 2022 AmeriFlux BASE US-Ics imnavait creek watershed wet sedge tundra, ver. 6-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246130

Euskirchen E, Shaver G and Bret-Harte S 2022 AmeriFlux BASE US-Ict imnavait creek watershed tussock tundra, ver. 4-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246131

Kobayashi H, Ikawa H and Suzuki R 2019 AmeriFlux BASE US-Prr poker flat research range black spruce forest, ver. 3-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246153

Randerson J 2016 AmeriFlux BASE US-Bn1 bonanza creek, 1920 burn site near delta junction, ver. 1-1 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246033

Randerson J 2016 AmeriFlux BASE US-Bn2 bonanza creek, 1987 burn site near delta junction, ver. 1-1 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246034

Randerson J 2016 AmeriFlux BASE US-Bn3 bonanza creek, 1999 burn site near delta junction, ver. 1-1 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246035

Torn M and Dengel S 2020 AmeriFlux BASE US-NGB NGEE Arctic barrow, ver. 3-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1436326

Torn M and Dengel S 2020 AmeriFlux BASE US-NGC NGEE Arctic council, ver. 1-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1634883

Ueyama M, Iwata H and Harazono Y 2019 AmeriFlux BASE US-Fcr cascaden ridge fire scar, ver. 2-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1562388

Ueyama M, Iwata H and Harazono Y 2022 AmeriFlux BASE US-Rpf poker flat research range: succession from fire scar to deciduous forest, ver. 7-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1579540

Ueyama M, Iwata H and Harazono Y 2022 AmeriFlux BASE US-Uaf University of Alaska, Fairbanks, ver. 10-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1480322

Sullivan P 2022 AmeriFlux BASE US-KPL Lily Lake Fen, ver. 1–5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1865478

Schuur T 2021 AmeriFlux BASE US-EML Eight Mile Lake permafrost thaw gradient, Healy Alaska., ver. 4-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1418678

Oechel W 2016 AmeriFlux BASE US-Hva happy Valley, ver. 2-1 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246064

Oechel W 2019 AmeriFlux BASE US-Upa Upad, Ver. 3-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1246108

NEON (National Ecological Observatory Network) 2022 AmeriFlux BASE US-xBA NEON barrow environmental observatory (BARR), ver. 5-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1671892

NEON (National Ecological Observatory Network) 2022 AmeriFlux BASE US-xBN NEON caribou creek—poker flats watershed (BONA), ver. 6-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1617727

NEON (National Ecological Observatory Network) 2022 AmeriFlux BASE US-xDJ NEON delta junction (DEJU), ver. 6-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1634884

NEON (National Ecological Observatory Network) 2022 AmeriFlux BASE US-xHE NEON Healy (HEAL), ver. 6-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1617729

NEON (National Ecological Observatory Network) 2022 AmeriFlux BASE US-xTL NEON Toolik (TOOL), ver. 5-5 AmeriFlux AMP, (Dataset) 10.17190/AMF/1617739

NEON (National Ecological Observatory Network) Soil temperature, RELEASE-2021 (DP1.00041.001) 10.48443/zayg-n802 (Dataset accessed from https://data.neonscience.org) (25 January 2021)

NEON (National Ecological Observatory Network) Triple aspirated air temperature (DP1.00003.001), RELEASE-2022 10.48443/q16j-sn13 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Soil water content and water salinity (DP1.00094.001), RELEASE-2022 10.48443/ghry-qw46 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Soil temperature (DP1.00041.001), RELEASE-2022 10.48443/9e56-pj39 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Soil CO₂ concentration (DP1.00095.001), RELEASE-2022 10.48443/h4xx-yz37 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Precipitation (DP1.00006.001), RELEASE-2022 10.48443/6wkc-1p05 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) 2D wind speed and direction (DP1.00001.001), RELEASE-2022 10.48443/77n6-eh42 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Shortwave and longwave radiation (net radiometer) (DP1.00023.001), RELEASE-2022 10.48443/stbf-bh38 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Groundwater and active layer measurements at permafrost sites (DP1.20099.001), RELEASE-2022 10.48443/1v5y-av98 (Dataset accessed from https://data.neonscience.org) (2 November 2022)

NEON (National Ecological Observatory Network) Soil physical and chemical properties, Megapit (DP1.00096.001), RELEASE-2022 10.48443/10dn-8031 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Albedo—spectrometer—mosaic (DP3.30011.001), RELEASE-2022 10.48443/j8yy-ky49 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Discrete return LiDAR point cloud (DP1.30003.001), RELEASE-2022 10.48443/ca22-1n03 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Elevation—LiDAR (DP3.30024.001), RELEASE-2022 10.48443/ymmp-fr93 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) High-resolution orthorectified camera imagery (DP1.30010.001), RELEASE-2022 10.48443/ksv5-ts70 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Soil physical and chemical properties, periodic (DP1.10086.001), RELEASE-2022 10.48443/fk4j-ax76 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Field spectral data (DP1.30012.001), RELEASE-2022 10.48443/qy50-d561 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Soil temperature (DP1.00041.001), RELEASE-2022 10.48443/9e56-pj39 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Dissolved gases in surface water (DP1.20097.001), RELEASE-2022 10.48443/4661-hh45 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Chemical properties of groundwater (DP1.20092.001), RELEASE-2022 10.48443/h7vt-3v55 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Litterfall and fine woody debris production and chemistry (DP1.10033.001), RELEASE-2022 10.48443/b2qt-7z79 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Soil CO₂ concentration (DP1.00095.001), RELEASE-2022 10.48443/h4xx-yz37 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Soil water content and water salinity (DP1.00094.001), RELEASE-2022 10.48443/ghry-qw46 (Dataset accessed from https://data.neonscience.org) (22 November 2022)

NEON (National Ecological Observatory Network) Bundled data products—eddy covariance (DP4.00200.001), RELEASE-2022 10.48443/7cqp-3j73 (Dataset accessed from https://data.neonscience.org) (22 November 2022)