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 physicsinformed 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 A distinctive feature of the Arctic environment is perennially frozen ground, or permafrost, which is defined as any subsurface material that remains below 0  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 2 in Eurasia, 10.2 million km 2 in North America (1); see. Fig. 1 for details. encompasses 1162 permafrost settlements (2), accommodating approximately 5 million inhabitants.Many objects of critical infrastructure are located on permafrost (3), including freeways, railroads, oil and gas pipelines, nuclear stations, which have been under the effect of permafrost thaw (4)(5)(6)(7). For example, according to the Circum-Arctic Map of Permafrost and Ground Ice Conditions by the International Permafrost Association (IPA Map) (8), permafrost comprises almost 65% of Russian territory, 80% of Alaska (9) and about 50% of Canada (10). Permafrost severely affects many aspects of daily life and industrial processes, e.g., mining, in these regions (11)(12)(13) and provide a substantial risk to the global supply chain (14). Furthermore, most of the critical infrastructure and houses in high-latitudes are built on permafrost (15), 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 (11,12,16,17). There is growing observational evidence that such effects are already taking place. Various studies (11)(12)(13)16) 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 (18,19). Permafrost carbon feedback, a phenomenon of significant concern, involves the release of substantial soil carbon deposits through permafrost thaw, amplifying the impacts on ecosystems (20)(21)(22)(23)(24). Additionally, it is worth noting that methane emissions, with a considerably stronger greenhouse effect compared to an equal amount of CO 2 , pose particular concerns in this context. Other studies in this field outline economic consequences of climate change related to permafrost. In (25), the authors approach this problem from the point of its influence on the macroeconomic variables. Studies (26,27) 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 (13,17,(28)(29)(30). Mitigating the consequences of permafrost degradation requires a method to estimate the permafrost characteristics under various climate change scenarios. This paper focuses on pre-dicting 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. 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 physicsbased 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 (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. A compendium of methodologies in permafrost modeling is available in (31). A simplest semiempirical approach considers two key parameters, ALT and MAGT, which are linked statisti-cally 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 (32) 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 (17) 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 (15). The authors used kriging to distribute the values over a spatial grid with a resolution of 0. A different numerical approach based on statistical models was described in (28). 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 +/-0.77  C for soil temperature and +/- 37 cm for the thickness of the seasonally thawed layer. An analogous approach was considered in (33) 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 (34) 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 (34). 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 (32,35,36) 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 (35). With a slight modification, the model has been used in many subsequent studies (15,17,32,(37)(38)(39). Despite the recent improvement in the availability of observational MAGT and ALT data through dedicated web-portal (40), 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 (35). Kudryavtsev model accounts for the snow, vegetation, organic and mineral substrate, and variable thermal properties (36,38). The model performs well in extensive validations by empirical observations (36). 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 (32). For a detailed description of the model, see (32,41). 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. 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. (CALM) program includes data from 265 sites in the Northern Hemisphere in 15 countries (42). 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 (35). Kudryavtsev model is a solution to the equation of heat flow theory (31), described by Eq. 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: where 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),  is the thermal conductivity of soil in the thawed state (W m 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 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 (44,45) 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: Statistics of the thaw depth collected at CALM sites . 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 x 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. 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 (21 km 2 ). We projected it to the 0.5-angular-degree grid. The time resolution is 1 month. 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 . These data are constant in time. 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. Zero Annual Amplitude -the depth at which the temperature change throughout the year was less than 0.1  C (49). 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 (40). 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. Let 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 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 p in different moments of time  . The time step is 1 month. Thus, P  P = {p  ij }, where each observation has spatial and time indices. Input data p for our model include the dynamic data on the temperature and precipitation and static information on the swampiness, biotopes, and permafrost type: where i, j -latitude and longitude of the center of a pixel, for which data are observed  -year of the observation X  ij -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 Our model takes X ij and X  ij as the inputs. The key variables for this study, alt  ij , magt  ij , were only recorded for a subset of the entire grid P. Furthermore, we used this subset for training and validation. Our model considers both observed meteorological data and output of Kudryavtsev model . The input data consist of 12 monthly measurements for temperature and precipitation X  ij , 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 and magt  ij using p  ij  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 (50), Multilayer Perceptron regression (51), Random Forest regression (52), Elastic Net regression (53) and CatBoost Regressor (54). We selected CatBoost Regressor (54) as the most reliable method. This model is a supervised learning meta-algorithm that constructs an ensemble of decision trees. reduces the bias and variance in predictions to produce a highquality 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 (55). We also provided rejection curves for the estimated uncertainty values based on datasets that were depleted of points with the highest predictive uncertainty . 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. 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 2 and 322 000 km 2 , respectively. These pixels on the map 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) (56) and the coefficient of determination (R 2 ) (57). The RMSE for Kudryavtsev model alone without using climate data was 32.39 +/- 0.17 cm for the ALT and 1.26 +/- 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 We applied our model to predict the extent of permafrost degradation under the optimistic CMIP 6 SSP245 scenario (updated RCP 4.5 (43)). The prediction for the permafrost degradation by 2050 is demonstrated in 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 environ- mental 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 (35), 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. The PDF file includes: Supplementary The PDF file includes: Supplementary Figure 1: Proportion of wetlands in each pixel. Supplementary Figure 2: One-dimensional representation of the vegetation types.