Assessing the Risk of Permafrost Degradation with Physics-Informed Machine Learning

Polina Pilyugina, Timofey Chernikov, Alexey Zaytsev,
Alexander Bulkin, Evgeny Burnaev, Ilya Belalov,
Nazar Sotiriadi, Yury Maximov, Oleg Anisimov

Global warming accelerates permafrost degradation, impacting the reliability of critical infrastructure used by more than five million people daily. Furthermore, permafrost thaw produces substantial methane emissions, further accelerating global warming and climate change and putting more than eight billion people at additional risk. To mitigate the upcoming risk, policymakers and stakeholders must be given an accurate prediction of the thaw development. Unfortunately, comprehensive physics-based permafrost models require location-specific fine-tuning that is challenging in practice. Models of intermediate complexity require few input parameters but have relatively low accuracy. The performance of pure data-driven models is low as well as the observational data is sparse and limited. In this work, we designed a physics-informed machine-learning approach for permafrost thaw prediction. The method uses a heat equation to regularize data-driven approach trained over permafrost monitoring data and climate projections. The latter leads to higher precision and better numerical stability allowing for reliable decision-making for construction and maintenance in the areas endangered by permafrost thaw with a time horizon of decades.

Introduction

A distinctive feature of the Arctic environment is perennially frozen ground, or permafrost, which is defined as any subsurface material that remains below $0^{\circ}$ C for more than two consecutive years. Permafrost can be found on land, in high mountains and even on the shelf of some Arctic seas. Collectively, these three types of permafrost currently occupy about a quarter of the Northern Hemisphere land surface, including about 16.7 million km² in Eurasia, 10.2 million km² in North America (?); see. Fig. 1 for details.

The Arctic Circumpolar Permafrost Region encompasses 1162 permafrost settlements (?), accommodating approximately 5 million inhabitants.Many objects of critical infrastructure are located on permafrost (?), including freeways, railroads, oil and gas pipelines, nuclear stations, which have been under the effect of permafrost thaw (?, ?, ?, ?). For example, according to the Circum-Arctic Map of Permafrost and Ground Ice Conditions by the International Permafrost Association (IPA Map) (?), permafrost comprises almost 65% of Russian territory, 80% of Alaska (?) and about 50% of Canada (?). Permafrost severely affects many aspects of daily life and industrial processes, e.g., mining, in these regions (?, ?, ?) and provide a substantial risk to the global supply chain (?). Furthermore, most of the critical infrastructure and houses in high-latitudes are built on permafrost (?), which makes it crucial to control and monitor the thaw.

The permafrost is characterized by the two major parameters: active layer thickness and ground temperature. The active layer is the uppermost soil above the permafrost, which thaws in the summer and re-freezes in the winter. The active layer is the soil above the permafrost, which thaws in the summer and freezes in the winter. The active layer thickness (ALT) is defined as the maximal seasonal thaw depth and measured at the end of a summer season. The mean annual ground temperature (MAGT) is defined as the temperature at the depth where the annual temperature amplitude is less than 0.1°C. Conventional permafrost degradation models connect the climate parameters, particularly air temperature and precipitation, with the ALT and MAGT. Soil and vegetation properties largely govern this relationship, and most permafrost models explicitly include them. Both ALT and MAGT are vulnerable to climate change. Global warming results in a higher ground temperature, which causes permafrost thaw and increases the active layer thickness. These changes in the permafrost state negatively affect the stability of buildings and infrastructures and decrease the bearing capacity of the building foundations (?, ?, ?, ?).

There is growing observational evidence that such effects are already taking place. Various studies (?, ?, ?, ?) report the effect of global warming on the destruction of buildings and infrastructure built on permafrost. Permafrost thaw also forms gas emission craters, which release large amounts of methane and other gases into the atmosphere (?, ?). Permafrost carbon feedback, a phenomenon of significant concern, involves the release of substantial soil carbon deposits through permafrost thaw, amplifying the impacts on ecosystems (?, ?, ?, ?, ?). Additionally, it is worth noting that methane emissions, with a considerably stronger greenhouse effect compared to an equal amount of CO₂, pose particular concerns in this context. Other studies in this field outline economic consequences of climate change related to permafrost. In (?), the authors approach this problem from the point of its influence on the macroeconomic variables. Studies (?, ?) have estimated that the cost of critical infrastructure damage due to permafrost thaw in North America could reach about $5 billion by the end of the 21^st century.A similar study for Northern Europe and Asia is given in (?, ?, ?, ?, ?).

Mitigating the consequences of permafrost degradation requires a method to estimate the permafrost characteristics under various climate change scenarios. This paper focuses on predicting permafrost changes related to climate change in permafrost-rich areas and combines physics-based mathematical models with data-driven methods aiming at improving their quality over limited data.

Contribution

This paper presents a novel physics-informed machine-learning approach for permafrost thaw prediction to address the challenges posed by global warming and its impact on critical infrastructure and climate change. The key contributions of this research are as follows:

•

a novel robust and accurate prediction framework combining the strengths of physics-based permafrost degradation models and machine learning. The solution is properly regularized and constrained by fundamental physical principles;
•

state-of-the-art prediction accuracy substantially improving existing algorithms and supporting risk assessment and decision-making;
•

extensive empirical validation and uncertainty estimation over Circumpolar Active Layer Monitoring (CALM) observational data.

Overall, our physics-informed machine-learning approach offers a powerful tool for predicting permafrost thaw, mitigating risks to critical infrastructure, and understanding the implications of permafrost degradation on global warming and climate change. This research contributes to the growing field of climate science and offers valuable insights for infrastructure development in high-latitude regions facing the challenges of permafrost thaw.

Algorithms and Models

A compendium of methodologies in permafrost modeling is available in (?). A simplest semi-empirical approach considers two key parameters, ALT and MAGT, which are linked statistically to the climatic parameters using historical observational data. Permafrost observations are spatially and temporally sparse and irregular, which poses additional challenges to a researcher. The lack of a unified standard for the observations further complicates the problem.

Another type of permafrost models (?) takes into account the climatic parameters and variability of snow, vegetation, and soil properties. The results demonstrate a good fit with the observational data. Such a model is used in (?) to estimate the geocryological risks of permafrost thawing. A similar model was used to simulate the temperature of the permafrost and depth of a seasonally thawed layer in (?). The authors used kriging to distribute the values over a spatial grid with a resolution of $0.25^{\circ}$ . Predictive calculations were made with the scenario of climate change derived from six CMIP 5 (Coupled Model Intercomparison Project Phase) Earth System Models under the assumption of the continuing growth of greenhouse gas emissions (so-called Representative Concentration Pathway (RCP) 8.5 scenario). Streletskiy et al. in (?) compared two periods: the present (2005-2015) and the middle of the 21st century (2050-2059).

A different numerical approach based on statistical models was described in (?). Hjort et al. computed characteristics such as the temperature, precipitation, organic carbon content in the soil, soil type, water bodies, solar radiation, and topographic properties for each point of the grid. These factors were predicted with the models from CMIP 5 climate scenario. The variables depend on greenhouse gas concentrations. Three scenarios adopted by the IPCC (RCP 2.6, RCP 4.5 and RCP 8.5) resulted in different values for the variables. In further analysis, these factors served as predictors for each scenario of greenhouse gas concentrations. The analysis included four statistical models (generalized linear model (GLM), generalized additive model (GAM), random forest (RF), and generalized boosting model (GBM)) and their ensemble. The analysis resulted in a map of soil the temperature and thickness of the seasonally thawed layer. The prediction uncertainty of the model based on the ensemble of four methods was estimated as $\pm 0.77^{\circ}$ C for soil temperature and $\pm$ 37 cm for the thickness of the seasonally thawed layer.

An analogous approach was considered in (?) with more focus on the modeling procedure. The predictive error (root mean squared error) of this approach, which was evaluated on the hindcast data for past periods, was reported to be 53 cm for ALT and 1.58°C for MAGT. To date, to the best of our knowledge, this is the best precision achieved in permafrost thaw prediction. However, a model with a prediction error of 53 cm is a poor guide for maintaining a structure built over permafrost.

Ultimately, the result of (?) imposes the inability of machine learning methods to predict environmental processes, including the permafrost thaw, at high latitude carbon cycle balances. The crucial contribution of this paper is regularizing the classical machine learning approaches with a special version of the heat equation (?). The latter substantially improves the prediction quality of permafrost thaw and opens the door to efficient yet accurate environmental risk assessment models in high latitudes.

There are many approaches for permafrost modeling and prediction, which may or may not include the effects of climate change, but the quality of the models limits their use in applications related to the maintenance of buildings and construction design for the long-term sustainable development of the northern regions.

In this study, we focus on spatial models and their interaction with climatic projections for the future. Kudryavtsev equilibrium model (?, ?, ?) is particularly useful. Equilibrium models suggest that the MAGT is in balance with the atmospheric parameters. They have relatively low data requirements; they use the mean monthly temperature and precipitation data as climate forcing with a few edaphic parameters that characterize the soil thermal properties, snow, and vegetation. The equilibrium model developed by V. Kudryavtsev is one of the most successful examples (?). With a slight modification, the model has been used in many subsequent studies (?, ?, ?, ?, ?, ?).

Despite the recent improvement in the availability of observational MAGT and ALT data through dedicated web-portal (?), permafrost modeling remains data-limited. Thus, more powerful transient models cannot be effectively used because all necessary forcing data and edaphic parameters are often not available, and the model cannot be appropriately calibrated. Given these limitations, in this study we used the equilibrium permafrost model.

The model is based on the numerical solution of the nonlinear parabolic equation by Kudryavtsev (?). Kudryavtsev model accounts for the snow, vegetation, organic and mineral substrate, and variable thermal properties (?, ?). The model performs well in extensive validations by empirical observations (?).

In the Kudryavtsev model, every grid element (a pixel on the map) is associated with observable variables such as the snow cover, vegetation, and soil characteristics. The model outputs the ground temperature and active layer thickness. In this study, we employed a modification of the Kudryavtsev model in (?). For a detailed description of the model, see (?, ?).

The input parameters of the model include characteristics of the atmosphere (monthly air temperatures and monthly precipitation) for climatic forcing. Additionally, the model configuration includes various soil specifications. By varying the main parameter — the type of soil (i.e., clay, loan, peat) — we obtained an ensemble of models. Other parameters, i.e.,, the snow cover and organic layer depth, were held constant. The organic layer depth is defined based on the persistent vegetation type in the area.

Materials and Methods

The key idea of the physics-informed machine learning approach for permafrost thaw degradation is to boost the prediction quality by combining state-of-the-art data-driven methods with physical-based models. The latter allows for building a high-fidelity thaw predictor over limited observations data.

The Circumpolar Active Layer Monitoring (CALM) program includes data from 265 sites in the Northern Hemisphere in 15 countries (?). Figure 1 demonstrates the locations of the CALM sites. As a part of the approach to assessing the risks of permafrost degradation, the seasonal thaw depth is modeled for all points of the grid.

We trained a machine learning model based on gradient boosting using the data from the CALM sites. The input data included dynamically measured parameters for climate, static data on soil and vegetation, and results for an ensemble of the Kudryavtsev model for different initializations of the parameters. Figure 2 shows the model flowchart. The first step consists of processing the input data and submitting the data to Kudryavtsev model (?).

Kudryavtsev model is a solution to the equation of heat flow theory (?), described by Eq. 1.

C\frac{\delta T}{\delta t}=k\frac{\delta^{2}t}{\delta z^{2}}

(1)

where $C$ is volumetric head capacity in Jm^-3, $T$ is temperature, $t$ is time and $z$ is depth in meters. The mathematical setup of the model is the following. The depth of seasonal thawing, $Z_{thaw}$ , is calculated using the semi-empirical equation:

Z_{thaw}=\frac{2(A_{s}-T_{z})\cdot\left[\frac{\lambda\cdot P_{sn}\cdot C}{\pi}\right]^{1/2}+\frac{(2A_{z}\cdot C\cdot Z_{c}+Q_{ph}\cdot Z_{c})\cdot Q_{ph}\cdot\left[\frac{\lambda\cdot P_{sn}}{\pi\cdot C}\right]^{1/2}}{2A_{z}\cdot C\cdot Z_{c}+Q_{ph}\cdot Z_{thaw}+(2A_{z}\cdot C+Q_{ph})\cdot\left[\frac{\lambda\cdot P_{sn}}{\pi\cdot C}\right]^{1/2}}}{2A_{z}\cdot C+Q_{ph}}

(2)

where

A_{z}=\frac{A_{s}-T_{z}}{\ln\left[\frac{A_{s}+Q_{ph}/2C}{T_{z}+Q_{ph}/2C}\right]}-\frac{Q_{ph}}{2C}\qquad\textrm{and}\qquad Z_{c}=\frac{2(A_{s}-T_{z})\cdot\sqrt{\frac{\lambda\cdot P_{sn}\cdot C}{\pi}}}{2A_{z}\cdot C+Q_{p}h}.

(3)

Here, $A_{s}$ is the annual amplitude of the soil-surface temperature (in ^∘ C), $T_{z}$ is the mean annual temperature at the depth of seasonal thawing (in ^∘ C), $\lambda$ is the thermal conductivity of soil in the thawed state ( $Wm^{-1}K^{-1}$ ), $C$ is the volumetric heat capacity of soil in the thawed state ( $Jm^{-3}K^{-1}$ ), $P_{sn}$ is the period of the temperature wave (sec), and $Q_{ph}$ is the volumetric latent heat of phase changes ( $Jm^{-3}$ ).

Then, we ran Kudryavtsev model with various initializations of parameters to obtain the ALT and MAGT values for each initialization. These values were forwarded as inputs to a Machine Learning (ML) model with other data, as described in the Data section. The ML model was subsequently trained and evaluated using available historical observations from the CALM and TSP stations. We used gridded climate data from the ensemble of CMIP6 Earth System Models to obtain predictions for the periods of 2010 – 2015 and 2040 – 2060. We considered two greenhouse-gas emission scenarios defined by so-called Shared socio-economical pathways (SSP). SSP-scenarios consider the impact of climate change on the development of societies and economies (?).

Refer to caption — Figure 1: Distribution of permafrost types according to the IPA Permafrost Map with CALM sites and TSP boreholes.

Data

Input data	Parameters	Units	Period
IPA Permafrost Map	Permafrost type	Categorical	N/A
swamp	Proportion of swamps in the area	%	N/A
CEDA	Temperatures, precipitation	^∘C, mm	1901-2020
WorldClim (historical)	Temperature, precipitation	^∘C, mm	1960-2018
CALM	Thaw depth	cm	1969-2021
			(irregular data)
WorldClim (CMIP6)	Temperature, precipitation	^∘C, mm	2006-2100
vegetation	Type of vegetation	Categorical	N/A
GTNP-TSP	Temperature, zero-amplitude depth	^∘C, cm	1901-2020

Table 1: Input data

The spatial resolution for historical data was 0.5 angle degrees, except for the GTNP-ALT dataset (see the description below). We used high-resolution temperature and precipitation data from WorldClim (?, ?) for future predictions and historical assessment. The temporal resolution for the time series data varied from 1 month to 1 year. We aggregated values of interest for each pixel with the datasets detailed below:

CEDA: (?)
Various climate variables: the cloud cover (average fraction of the sky obscured by clouds), daily temperature range, proportion of time with negative temperature in a 24-hour period, precipitation, monthly average daily minimum temperature (the minimum temperature was taken daily and averaged over the month), monthly average maximum daily temperature, monthly average temperature, proportion of time during which precipitation occurs. The temporal resolution for the data is 1 month.
CALM: (?)
Statistics of the thaw depth collected at CALM sites (Figure 1). Each data point has eight entries. Seven of those describe the characteristics of the thaw depth: average, median, minimum, maximum, 25th percentile, 75th percentile, and standard deviation. For each point, measurements were performed at the vertices of a rectangular grid of 1 km $\times$ 1 km with 100-m steps; there were 121 measurements in total. In some cases, several measurements were missing. Data from some points were irregularly collected or during short time intervals. To bring the data to a common format with a spatial resolution of 0.5 angle degrees, we considered the nearest CALM site for each point on the map. This distance is the eighth entry for each data point. The maximum distance between the CALM site and the corresponding point on the map was 200 km (4 cases). The temporal resolution for this dataset is 1 year.
WorldClim (historical): (?)
High-resolution historical climatic data projected from CRU-TS-4.03 using WorldClim 2.1 for bias correction (Climatic Research Unit, University of East Anglia). The data cover the period from 1960 to 2018. The spatial resolution of the original dataset is 2.5 minutes ( $\sim$ 21 km²). We projected it to the 0.5-angular-degree grid. The time resolution is 1 month.
IPA: (?)
Map of the distribution of permafrost types. We rasterized the original map in the vector format to bring it to a resolution of 0.5 angle degrees to fit the common format for other data (Figure 1). These data are constant in time.
WorldClim (CMIP6): (?)
The dataset from the CanESM5 model contains data on the air temperature (mean, minimum, and maximum) in ^∘C and precipitation (mm). The time resolution is 1 month.
swamp: These data were kindly provided by Professor Oleg Anisimov. The data show the percentage of wetlands for each point on the map (see Supplementary Figure 1). The data are constant in time.
vegetation: These data were kindly provided by Professor Oleg Anisimov. The data describe the types of landscape or biotopes (see Fig. S2). These data are constant in time.
GTNP-TSP: Zero Annual Amplitude – the depth at which the temperature change throughout the year was less than $0.1^{\circ}$ C (?). These data describe a thermal state of the permafrost, which was measured at boreholes (TSP boreholes in Figure 1). Since the exact depth of the ZAA is not explicitly reported in the dataset, we manually determine the ZAA depth for each weather station and for each year using the GTNP dataset (?). Stations with insufficient data to define the ZAA were omitted. Then, we projected the resulting dataset onto a regular two-dimensional grid to create a map. The final dataset includes both temperature at the ZAA level and ZAA depth value.

Data for modeling

Let $\mathcal{P}$ be the domain. This is the part of the Earth’s surface underlaid by permafrost. We project this two-dimensional locally continuous surface to the grid $\mathbf{P}$ parameterized by two angles $i$ and $j$ . The steps for both angular parameters are 0.5 degrees. For each pixel on the map with coordinates $i$ and $j$ , we have observations of different features $\mathbf{p}$ in different moments of time $\tau$ . The time step is 1 month. Thus, $\mathcal{P}\supset\mathbf{P}=\{\mathbf{p}_{ij}^{\tau}\}$ , where each observation has spatial and time indices.

Input data $\mathbf{p}$ for our model include the dynamic data on the temperature and precipitation and static information on the swampiness, biotopes, and permafrost type:

\mathbf{p}_{ij}^{\tau}=\{i,j,\tau,X_{ij}^{\tau},X_{ij},Y_{ij}^{\tau}\}

(4)

where

$i$ , $j$ - latitude and longitude of the center of a pixel, for which data are observed
$\tau$ - year of the observation
$X_{ij}^{\tau}$ - 12 monthly values of climatic variables (min and max temperatures, monthly precipitation)
$X_{ij}$ - stationary parameters, which do not depend on time: swampiness, biotopes, permafrost type
$Y_{ij}^{\tau}=\{alt_{ij}^{\tau},magt_{ij}^{\tau}\}$ - observed ALT and MAGT values

Our model takes $X_{ij}$ and $X_{ij}^{\tau}$ as the inputs. The key variables for this study, $alt_{ij}^{\tau},magt_{ij}^{\tau}$ , were only recorded for a subset of the entire grid $\mathbf{P}$ . Furthermore, we used this subset for training and validation.

Implementation

Our model considers both observed meteorological data and output of Kudryavtsev model (see Fig. 2). The input data consist of 12 monthly measurements for temperature and precipitation $X_{ij}^{\tau}$ , the thickness of an organic layer, and soil type $X_{ij}$ .

We used variations of the latter to parametrize Kudryavtsev model. The initialization of Kudryavtsev model contains parameters for four soil types: sand, loam, clay, and peat.The dominant soil type in the permafrost areas is clay soil, which has similar physicochemical parameters to loam. Therefore, we obtained the model predictions for both soil types. Moreover, we obtained the model predictions for dry and wet soils. For swamp-associated locations, we used the average for the initializations with the regular soil type and peat, which were weighted according to the swamp dataset (see Supplementary Figure). Thus, for each pixel, we obtained four pairs of ALT and MAGT values to initialize Kudryavtsev model.

Our model predicts the values for $alt_{ij}^{\tau}$ and $magt_{ij}^{\tau}$ using $\mathbf{p}_{ij}^{\tau}\in\mathbf{P}_{alt}$ and initializations of Kudryavtsev model from the previous step.

We conducted experiments with various machine Learning models to evaluate their performance and select the most suitable model. These models are: Linear regression (?), Multi-layer Perceptron regression (?), Random Forest regression (?), Elastic Net regression (?) and CatBoost Regressor (?). We selected CatBoost Regressor (?) as the most reliable method. This model is a supervised learning meta-algorithm that constructs an ensemble of decision trees. The CatBoost Regressor reduces the bias and variance in predictions to produce a high-quality and diverse ensemble.

For uncertainty estimation, we used the algorithms implemented in Catboost. We trained the model using a specific loss function, which enables uncertainty estimation via virtual ensembles, according to the original study (?). We also provided rejection curves for the estimated uncertainty values based on datasets that were depleted of points with the highest predictive uncertainty (Fig. 3).

Predictions of permafrost degradation in the future, which is the purpose of our model development, uses data on the temperature and precipitation from the ensemble of CMIP6 earth system models instead of historical data. Therefore, results for Kudryavtsev model were updated according to CMIP6 data.

Results

In total, we used 2 729 (ALT) and 961 (MAGT) data points, which were collected in 1990 – 2010, from the areas of 364 000 km² and 322 000 km², respectively. These pixels on the map contain or are neighbors to the CALM sites and/or TSP boreholes. To estimate the performance of our model and its ability to predict future values, we divided our dataset into training and test subsets. The training subsets preceded 2013, while the years for the test were 2013 – 2020. Figure 4 demonstrates an example of the training-testing split that we used. Furthermore, we used the K-fold algorithm to divide the training set into five folds of 80% randomly selected data points.

Tables 2 and 3 show the prediction errors for different models that were constructed using four sets of climatic factors with or without the outputs of Kudryavtsev models. We focus our evaluation of model performance on two metrics: the root mean square error (RMSE) (?) and the coefficient of determination ( $R^{2}$ ) (?).

The RMSE for Kudryavtsev model alone without using climate data was 32.39 $\pm$ 0.17 cm for the ALT and 1.26 $\pm$ 0.02 for the MAGT. This RMSE is a significant improvement over the existing approaches. When we accounted for all data in the data section, the test RMSE improved to 25.53 $\pm$ 0.40 cm for the ALT and 1.08 $\pm$ 0.02 for the MAGT.

Model	Factors	Train	Test	Train	Test
		RMSE		R² ( $\times 100$ )
CatBoost	All data	0.72 $\pm$ .11	1.08 $\pm$ .02	89 $\pm$ 3	53 $\pm$ 1
	Kudr. (ALT + MAGT)	1.04 $\pm$ .02	1.26 $\pm$ .02	77 $\pm$ 2	37 $\pm$ 2
	Kudr. (ALT)	1.27 $\pm$ .03	1.32 $\pm$ .01	67 $\pm$ 2	30 $\pm$ 1
	Kudr. (MAGT)	1.06 $\pm$ .08	1.22 $\pm$ .04	76 $\pm$ 3	41 $\pm$ 4
	Only climate data	0.64 $\pm$ .06	1.14 $\pm$ .02	91 $\pm$ 2	48 $\pm$ 2
Elastic Net	All data	1.24 $\pm$ .02	1.27 $\pm$ .03	68 $\pm$ 1	35 $\pm$ 3
	Kudr. (ALT + MAGT)	1.50 $\pm$ .02	1.25 $\pm$ .01	53 $\pm$ 1	38 $\pm$ 1
	Kudr. (ALT)	1.80 $\pm$ .04	1.50 $\pm$ .02	32 $\pm$ 1	10 $\pm$ 3
	Kudr. (MAGT)	1.61 $\pm$ .02	1.40 $\pm$ .01	46 $\pm$ 1	21 $\pm$ 1
	Only climate data	1.33 $\pm$ .02	1.55 $\pm$ .04	63 $\pm$ 1	4 $\pm$ 6
LinearRegression	All data	1.18 $\pm$ .02	1.41 $\pm$ .07	71 $\pm$ 1	20 $\pm$ 8
	Kudr. (ALT + MAGT)	1.42 $\pm$ .02	1.18 $\pm$ .02	58 $\pm$ 2	42 $\pm$ 2
	Kudr. (ALT)	1.55 $\pm$ .03	1.19 $\pm$ .01	50 $\pm$ 2	44 $\pm$ 1
	Kudr. (MAGT)	1.52 $\pm$ .02	1.17 $\pm$ .01	52 $\pm$ 1	45 $\pm$ 1
	Only climate data	1.29 $\pm$ .02	1.64 $\pm$ .05	65 $\pm$ 1	-7 $\pm$ 7
NeuralNetwork	All data	1.06 $\pm$ .07	1.56 $\pm$ .13	77 $\pm$ 3	2 $\pm$ 16
	Kudr. (ALT + MAGT)	1.41 $\pm$ .03	1.39 $\pm$ .05	59 $\pm$ 2	23 $\pm$ 5
	Kudr. (ALT)	1.50 $\pm$ .09	1.34 $\pm$ .01	53 $\pm$ 4	28 $\pm$ 1
	Kudr. (MAGT)	1.45 $\pm$ .05	1.35 $\pm$ .11	57 $\pm$ 2	27 $\pm$ 11
	Only climate data	1.05 $\pm$ .10	1.40 $\pm$ .23	77 $\pm$ 4	19 $\pm$ 26
RandomForest	All data	0.38 $\pm$ .01	1.16 $\pm$ .03	97 $\pm$ 1	47 $\pm$ 2
	Kudr. (ALT + MAGT)	0.45 $\pm$ .02	1.24 $\pm$ .04	96 $\pm$ 1	38 $\pm$ 4
	Kudr. (ALT)	0.51 $\pm$ .02	1.41 $\pm$ .04	94 $\pm$ 1	21 $\pm$ 5
	Kudr. (MAGT)	0.48 $\pm$ .01	1.27 $\pm$ .04	95 $\pm$ 1	35 $\pm$ 4
	Only climate data	0.37 $\pm$ .01	1.28 $\pm$ .08	97 $\pm$ 1	35 $\pm$ 8

Table 2: Mean annual ground temperatures (MAGT) prediction errors on the proposed training-testing splits. For each configuration, we trained the model with 80% of the training data and tested it on the test data. The values show the average of five estimations for each model configuration. Here, Kudr. denotes the Kudryavtsev model results. The best result is shown in bold.

Model	Factors	Train	Test	Train	Test
		RMSE		R² ( $\times 100$ )
CatBoost	All data	5.65 $\pm$ .90	25.53 $\pm$ .40	98 $\pm$ 1	62 $\pm$ 1
	Kudr. (ALT + MAGT)	27.32 $\pm$ .55	32.39 $\pm$ .17	61 $\pm$ 2	38 $\pm$ 1
	Kudr. (ALT)	25.71 $\pm$ 1.69	36.71 $\pm$ .13	65 $\pm$ 5	21 $\pm$ 1
	Kudr. (MAGT)	27.41 $\pm$ 1.48	33.62 $\pm$ .54	60 $\pm$ 4	33 $\pm$ 2
	Only climate data	9.96 $\pm$ 1.63	28.57 $\pm$ .47	95 $\pm$ 2	52 $\pm$ 2
Elastic Net	All data	29.80 $\pm$ .31	37.79 $\pm$ .18	53 $\pm$ 1	16 $\pm$ 1
	Kudr. (ALT + MAGT)	38.93 $\pm$ .39	40.17 $\pm$ .18	20 $\pm$ 1	5 $\pm$ 1
	Kudr. (ALT)	39.79 $\pm$ .35	40.80 $\pm$ .39	17 $\pm$ 1	2 $\pm$ 2
	Kudr. (MAGT)	39.48 $\pm$ .40	39.96 $\pm$ .18	18 $\pm$ 1	6 $\pm$ 1
	Only climate data	31.14 $\pm$ .28	38.48 $\pm$ .27	49 $\pm$ 1	13 $\pm$ 1
LinearRegression	All data	28.48 $\pm$ .23	36.02 $\pm$ .23	57 $\pm$ 1	24 $\pm$ 1
	Kudr. (ALT + MAGT)	36.43 $\pm$ .27	39.76 $\pm$ .38	30 $\pm$ 1	7 $\pm$ 2
	Kudr. (ALT)	39.13 $\pm$ .38	40.76 $\pm$ .15	20 $\pm$ 1	2 $\pm$ 1
	Kudr. (MAGT)	38.23 $\pm$ .36	39.62 $\pm$ .33	23 $\pm$ 1	7 $\pm$ 2
	Only climate data	31.07 $\pm$ .28	38.93 $\pm$ .28	49 $\pm$ 1	11 $\pm$ 1
NeuralNetwork	All data	19.21 $\pm$ 1.30	29.53 $\pm$ .60	81 $\pm$ 3	49 $\pm$ 2
	Kudr. (ALT + MAGT)	33.50 $\pm$ .62	39.31 $\pm$ 1.88	41 $\pm$ 3	9 $\pm$ 9
	Kudr. (ALT)	38.65 $\pm$ .48	38.35 $\pm$ .31	21 $\pm$ 1	13 $\pm$ 1
	Kudr. (MAGT)	31.13 $\pm$ .94	36.73 $\pm$ .28	49 $\pm$ 3	21 $\pm$ 1
	Only climate data	17.81 $\pm$ 2.89	32.30 $\pm$ 1.14	83 $\pm$ 5	38 $\pm$ 4
RandomForest	All data	7.36 $\pm$ .12	29.57 $\pm$ .32	97 $\pm$ 1	48 $\pm$ 1
	Kudr. (ALT + MAGT)	10.43 $\pm$ .16	34.53 $\pm$ .26	94 $\pm$ 1	30 $\pm$ 1
	Kudr. (ALT)	12.89 $\pm$ .12	37.05 $\pm$ .16	91 $\pm$ 1	19 $\pm$ 1
	Kudr. (MAGT)	11.58 $\pm$ .28	35.03 $\pm$ .44	93 $\pm$ 1	28 $\pm$ 2
	Only climate data	7.62 $\pm$ .12	31.70 $\pm$ .28	97 $\pm$ 1	41 $\pm$ 1

Table 3: Active layer thickness (ALT) prediction errors on the training-testing splits. For each configuration, we trained the model with 80% of the training data and tested it on the test data. The values show the average of five estimations for each model configuration. Here, Kudr. denotes the Kudryavtsev model results. The best result is shown in bold.

Figure 5 shows the predicted values by the CatBoost model, which was trained with all factors included versus the actual active layer thickness computed using all climatic data in the Data section. Although our model demonstrates excellent predictive power, it slightly underestimates the ALT.

We applied our model to predict the extent of permafrost degradation under the optimistic CMIP 6 SSP245 scenario (updated RCP 4.5 (?)). The prediction for the permafrost degradation by 2050 is demonstrated in Figure 6.

Figures 7 and 8 demonstrate the predicted active layer thickness and temperature of the soil at the zero-amplitude level for 2050 (under the CMIP 6 SSP245 scenario) and as a comparison to 2010.

Discussion

Significant progress in machine learning algorithms and methods for natural sciences cannot be ignored while studying climate change and permafrost degradation as a key part of it. While promising in predicting environmental changes, machine learning methods demonstrate mediocre quality over sparse, limited, and incomplete data. The latter is intrinsic for sparsely populated Arctic areas and permafrost melting there.

The major challenge in data-driven permafrost modeling and thaw depth prediction comes from sparse and irregular data on the active layer thickness dynamics. Moreover, the locations of CALM sites do not cover the actual temperature range for all types of landscapes or soils.

As such, classical physics-based models are still dominant in permafrost thaw and environmental modeling in the Arctic. In this paper, we are the first to demonstrate the power of machine learning algorithms properly augmented by physical models for permafrost analysis and long-term thaw prediction. In particular, we show that a coarse-grained physical model, such as a Kudryavtsev solution to the heat equation (?), used as a regularization to a machine learning method can improve the performance of both machine learning and physics-based models by 20% or more.

We believe the latter opens the door to studying a wide range of environmental processes in high latitudes from the physics-informed machine learning viewpoint, as the latter may substantially improve existing algorithms and practices. Furthermore, it will provide stakeholders with a high-fidelity tool to mitigate the climate change impact and protect people and the planet.

Supplementary materials

Figs. S1-S2

References

1. J. Obu, Journal of Geophysical Research: Earth Surface 126, 1 (2021).
2. J. Ramage, et al., Population and Environment 43, 22 (2021).
3. J. Hjort, et al., Nature Reviews Earth & Environment 2022 3:1 3, 24 (2022).
4. D. Reimchen, G. Doré, D. Fortier, R. Walsh, 2009 Annual Conference and Exhibition of the Transportation Association of Canada-Transportation in a Climate of Change (2009).
5. P. H. Larsen, et al., Global Environmental Change 18, 442 (2008).
6. A. M. Melvin, et al., Proceedings of the National Academy of Sciences 114, E122 (2017).
7. J. Ramage, et al., Population and Environment 43, 22 (2021).
8. J. Brown, O. Ferrians Jr, J. A. Heginbottom, E. Melnikov, Circum-Arctic Map of Permafrost and Ground-Ice Conditions, version 2 uuser guide (2002).
9. M. Jorgenson, et al. (2008).
10. O. J. Ferrians, The Geology of North America 1, 845 (1994).
11. D. Streletskiy, O. Anisimov, A. Vasiliev, Permafrost Degradation, vol. December (Springer, 2015).
12. D. A. Streletskiy, A. B. Sherstiukov, O. W. Frauenfeld, F. E. Nelson, Environmental Research Letters 10 (2015).
13. D. Streletskiy, Snow and Ice-Related Hazards, Risks, and Disasters February, 297 (2021).
14. A. Alexander, C. Blome, M. C. Schleper, S. Roscoe, International Journal of Operations & Production Management 42, 1061 (2022).
15. D. A. Streletskiy, L. J. Suter, N. I. Shiklomanov, B. N. Porfiriev, D. O. Eliseev, Environmental Research Letters 14 (2019).
16. O. A. Anisimov, V. Kokorev, D. A. Streletskiy, N. Shiklomanov, Assessment Report The Main Natural and Socio-economic Consequences of Climate Change in Permafrost Areas: A Forecast Based upon a Synthesis of Observations and Modelling, Assessment Report – Greenpeace USA (2014).
17. O. Anisimov, D. Streletskiy, Arctika XXI Century 2, 60 (2015).
18. V. Bogoyavlensky, Funnels of Yamal and Taymyr. Drill. Oil 10, 4 (2014).
19. A. Khimenkov, et al., Natural Hazards and Risk Research in Russia (Springer, 2019), pp. 305–316.
20. S. M. Natali, et al., Nature Climate Change 9, 852 (2019).
21. E. A. Schuur, et al., Nature 520, 171 (2015).
22. L. Mattsson, O. A. Anisimov, Environmental Research Letters 2, 045016 (2007).
23. T. R. Christensen, et al., Geophysical Research Letters 31 (2004).
24. K. R. Miner, et al., Nature Reviews Earth and Environment 3, 55 (2022).
25. B. Porfiriev, et al., Arctic: Ecology and Economy 4, 4 (2017).
26. P. H. Larsen, et al., Global Environmental Change 18, 442 (2008).
27. A. M. Melvin, et al., Proceedings of the National Academy of Sciences of the United States of America 114, E122 (2017).
28. J. Hjort, et al., Nature Communications 9 (2018).
29. N. I. Shiklomanov, D. A. Streletskiy, T. B. Swales, V. A. Kokorev, Geographical Review 107, 125 (2017).
30. S. V. Badina, Polar Science 29, 100730 (2021).
31. D. Riseborough, N. Shiklomanov, B. Etzelmüller, S. Gruber, S. Marchenko, Permafrost and Periglacial Processes 19, 137 (2008).
32. O. A. Anisimov, Kriosfera Zemli February, 36 (2016).
33. J. Aalto, O. Karjalainen, J. Hjort, M. Luoto, Geophysical Research Letters 45, 4889 (2018).
34. I. A. Shirley, Z. A. Mekonnen, R. F. Grant, B. Dafflon, W. J. Riley, Environmental Research Letters 18, 014026 (2023).
35. V. Kudryavtsev, L. Garagulya, V. Melamed, et al., Fundamentals of frost forecasting in geological engineering investigations (osnovy merzlotnogo prognoza pri inzhenerno-geologicheskikh issledovaniyakh), Tech. rep., COLD REGIONS RESEARCH AND ENGINEERING LAB HANOVER NH (1977).
36. N. Shiklomanov, F. Nelson, Ecological Modelling 123, 105 (1999).
37. E. E. Jafarov, S. S. Marchenko, V. E. Romanovsky, The Cryosphere 6, 613 (2012).
38. T. Sazonova, V. Romanovsky, Permafrost and Periglacial Processes 14, 125 (2003).
39. D. A. Streletskiy, N. I. Shiklomanov, V. I. Grebenets, Earth’s Cryosphere 16, 22 (2012).
40. GTNP - global terrestrial network for permafrost, http://gtnpdatabase.org/activelayers.
41. D. Riseborough, N. Shiklomanov, B. Etzelmüller, S. Gruber, S. Marchenko, Permafrost and Periglacial Processes 19, 137 (2008).
42. N. Shiklomanov, G. Washington, The circumpolar active layer monitoring (calm) program: Data collection, management, and dissemination strategies informal roads: The impact of unofficial transportation routes on remote arctic communities view project arctic pire view project dmitry a streletskiy.
43. K. Riahi, et al., Global environmental change 42, 153 (2017).
44. S. E. Fick, R. J. Hijmans, International journal of climatology 37, 4302 (2017).
45. I. Harris, P. D. Jones, T. J. Osborn, D. H. Lister, International journal of climatology 34, 623 (2014).
46. University of east anglia climatic research unit; harris, i.c.; jones, p.d.; osborn, t. (2021): Cru ts4.05: Climatic research unit (cru) time-series (ts) version 4.05 of high-resolution gridded data of month-by-month variation in climate (jan. 1901- dec. 2020). nerc eds centre for environmental data analysis, https://catalogue.ceda.ac.uk/uuid/c26a65020a5e4b80b20018f148556681.
47. J. Brown, O. Ferrians, J. A. Heginbottom, E. Melnikov, Circum-arctic map of permafrost and ground-ice conditions, version 2 (2002).
48. World climate research programme - cmip6 https://www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6.
49. V. E. Romanovsky, S. L. Smith, H. H. Christiansen, Permafrost and Periglacial processes 21, 106 (2010).
50. F. Galton, The Journal of the Anthropological Institute of Great Britain and Ireland 15, 246 (1886).
51. G. E. Hinton, Machine learning (Elsevier, 1990), pp. 555–610.
52. L. Breiman, Machine learning 45, 5 (2001).
53. H. Zou, T. Hastie, Journal of the royal statistical society: series B (statistical methodology) 67, 301 (2005).
54. L. Prokhorenkova, G. Gusev, A. Vorobev, A. V. Dorogush, A. Gulin, Advances in neural information processing systems 31 (2018).
55. A. Malinin, L. Prokhorenkova, A. Ustimenko, International Conference on Learning Representations (2020).
56. J. S. Armstrong, F. Collopy, International Journal of Forecasting 8, 69 (1992).
57. R. A. Fisher, Metron 1, 3 (1921).

Supplementary materials

The PDF file includes:

Figures S1 to S2