Climate Change Impact on Agricultural Land Suitability: An Interpretable Machine Learning-Based Eurasia Case Study

This study addresses the critical global issue of food security, particularly under the influence of climate change on agricultural land suitability. The primary objective of our research is to predict the risks associated with land suitability degradation and changes in irrigation patterns, directly impacting food security. This research aligns with the United Nations’ sustainable development goals to reduce hunger and malnutrition. Central Eurasia, a region facing unique economic and social challenges, serves as the focal point of our investigation, providing a pertinent example for analyzing the effects of climate change on food security. In our approach, we employ interpretable machine learning techniques to analyze the impact of climate change on agricultural land suitability under different carbon emission scenarios. The developed model demonstrates strong performance, evidenced by an accuracy of 86% and a mean average precision of 72% in a multi-class land suitability classification task. Focusing on the most vulnerable regions in Eastern Europe and Northern Asia, our research provides crucial insights for policymakers. These insights are instrumental for strategic planning, including the allocation of critical resources like water and fertilizers, to prevent humanitarian crises. The results demonstrate that machine learning can be a powerful tool in predicting and managing the impacts of climate change on food security.


I. INTRODUCTION
The impact of global climate change extends across various spheres of human activity, significantly influencing global pandemics, food insecurity, and political as well as social instability.As the Earth's land surface contends with The associate editor coordinating the review of this manuscript and approving it for publication was Li He .
concerning temperature increases and declines in humidity, nearly 40% of the planet's croplands and pastures are under threat [1].Projections indicate that global food demand may surge by approximately 110% by the year 2050 [2], [3].Compounding this issue, the ascending mean temperature, snow-water equivalent, and greenhouse gas concentrations pose substantial challenges to maintaining optimal conditions for crop cultivation [4], [5], [6].
This study examines the impact of climate change on agricultural sustainability following the Shared Socioeconomic Pathways (SSP) [7] representing probable paths of carbon emissions.The Coupled Model Intercomparison Project (CMIP) affords a broad set of models describing evaluation of key climate indicators such as mean temperature, humidity, and atmospheric pressure.We have utilized three CMIP6 ensembles-CMCC-ESM2 (Centro Euro-Mediterraneo sui Cambiamenti Climatici, Italy), CNRM-CM6-1 (Centre National de Recherches Météorologiques, France), and MRI-ESM2-0 (Meteorological Research Institute Earth System Model, Japan) [8], [9], [10].Each ensemble was analyzed under three distinct SSP scenarios: SSP1-2.6, which envisions a sustainable, low-emission green energy future; SSP2-4.5, portraying a 'business-as-usual' trajectory with moderate emissions; and SSP5-8.5, reflecting a high fossil fuel dependency scenario with significantly increased greenhouse gas emissions.Additionally, we have integrated the Global Food Security-support Analysis Data at a 1 km x 1 km resolution (GFSAD1km) to scrutinize cropland irrigation conditions [11].
Guided by a high demand from policymakers, academic and industrial stakeholders our study seeks to determine: (1) the potential effects of climate change under various SSP scenarios on the geographical distribution of croplands; (2) the regions with significant food security risks; and (3) the key factors influencing the suitability of cropland, particularly in relation to the type of irrigation employed.
In pursuit of these objectives, we have adopted a stacked ensemble methodology, combining eight state-of-the-art machine learning models to derive a cohesive meta-classifier.Such an amalgamation is generally observed to yield superior performance relative to any solitary model, as delineated by Mohandes et al. [12].This technique enables the precise classification of agricultural land into four distinct categories-major irrigation, minor irrigation, rainfed, and non-cropland-reflective of prevailing climatic conditions.Our research forecasts a marked enlargement in territories apt for agriculture by 2050, with a predilection for rainfed agriculture, potentially mitigating food security apprehensions in susceptible regions.Additionally, we predict a northward shift in the zones of agricultural suitability, alongside an elevated jeopardy for currently cultivated lands.
Employing the meta-classifier, we delineate the pivotal factors influencing the suitability of agricultural land, encompassing alterations in irrigation and land-use patterns.Our investigation offers invaluable insights into the prospective transformations in land suitability over the forthcoming decades, pinpointing the dominant elements that will shape this progression.In particular, our research: • Elucidates the intrinsic linkages between the climate parameters and the associated risks to agricultural land utilization; • Explores the progression of agricultural land, spotlighting shifts in irrigation modalities and the hazards of suitability diminution; • Identifies salient factors that affect the land's agricultural potential through the application of machine and deep learning paradigms.
Our conclusions will be instrumental for policymakers and land-use strategists dedicated to developing efficacious management plans that take into account the complex interconnections between climate change, land suitability, and agricultural output over an extended period.
Section II provides an overview of the research related to our study.The methodology, data employed, and the particulars of the implemented algorithms are explicated in Section III.Section IV elucidates the assessment of our models, along with projections and analyses for prospective scenarios.The limitations of our current study and potential areas for future research are discussed in Section V. We conclude the paper with a summary of our main findings in Section VI.

II. RELATED WORK
Relevant studies can be divided into the following categories: climate and remote sensing data related to land usage, artificial neural networks, and machine learning paradigms for geospatial data, the latter applied to cropland analysis, feature importance, and, finally, climate projections.

A. CLIMATE AND SATELLITE DATA
Climate data related to agriculture are typically represented in the form of geospatial datasets.For example, historical global reanalysis models such as ERA5 [13] and TerraClimate [14] are represented using a spatio-temporal grid that covers the Earth uniformly.ERA5, with a spatial resolution of 31 × 31 kilometers, has extensive temporal resolution support for 3-hour, 1-day, and 1-month resolutions and is continually refined in physical models, data assimilation, and core dynamics.TerraClimate has a spatial resolution of 4 × 4 kilometers and limited to monthly temporal resolution.Also, remote sensing data is of frequent use in cropland analysis [15].The specific examples of such data are the Normalized Difference Vegetation Index (NDVI), different types of Normalized Difference Indexes (NDI47, NDI45), C-band radar backscatter data from satellites, and many others [16].

B. CLIMATE PROJECTIONS
Climate Model Intercomparison Project phase 6 (CMIP6) that began in 2013 [17] sets a framework that allows for the numerical modeling climate projections extending into the future.Such projections allow one to describe the trends and average behavior of the Earth's atmosphere, land, and ocean interconnections from 2015 up to the year 2100, based on the assumed patterns of human behavior.Using this data, Shoaib et al. [18] investigated crop yield change in China under different representative concentration pathways, namely RCP 4.5, RCP 6.0, and RCP 8.5 using temperature and rainfall data.They employed a linear regression model to analyze the crop yield based on the World Bank dataset [19].Another study by Müller et al. [20] analyzed a set of 79 CMIP5 and CMIP6 projections of pressure and temperature, their statistics, and the impact on the crop yield.It is worth mentioning the original IPCC special report [21] that provides a general overview of the climate projections regarding cropland suitability.This report outlines recommendations for policymakers on, ideally, preventing global average temperature increase.Finally, there is a common practice of assembling multi-model ensemble projections by combining several CMIP projections [22], [23], [24], [25].Essentially, such approaches aim to construct a distribution of specific variables of the global climate models, with the distribution is estimated empirically from the variability between multi-model ensembles.The value of a climate variable is defined by the mean of the ensemble members and its variance can be estimated as well.

C. CROPLAND DATA
In the open source, there are several prominent examples of datasets that either contain or address data on the suitability of land usage.One such example is the Global Food Security Support Analysis Data with 1km resolution (GFSAD1km) [11], which describes landcover type for each grid cell in a single 2010 year.There are non-cropland and cropland of 5 types, depending on the irrigation pattern.The previously mentioned ERA5 [13] partially addresses the potential land usage by having a soil type parameter, which can indicate pixels that have organic or tropical organic soil within.Moreover, there is a sub-type of CMIP projections: Land Usage MIP (LUMIP) [26] with yearly temporal resolution.A portion of such projection describes land cover type in fracLut variable -the fraction of the grid cell for each landuse tile.In other words, it represents the portion of a pixel covered with urban, cropland, pasture, and the other types.A large portion of projections either do not allow these values to vary or do not describe the irrigation required.

D. MACHINE LEARNING AND NEURAL NETWORKS FOR GEOSPATIAL DATA
Research in the field of geospatial data analysis spans a broad range of topics, from developing innovative architectures to focusing on specific modeling techniques.Studies [27], [28] explore a variety of artificial neural networks, including Long Short-Term Memory (LSTM) [29] and Multi-Layer Perceptron (MLP) [30], for temperature forecasting applications.Dikshit et al. [31] evaluate several machine learning approaches, including Artificial Neural Networks (ANN), Support Vector Machines (SVM), Adaptive Neuro-Fuzzy Inference Systems (ANFIS), Extreme Learning Machines (ELM), decision trees, and random forests, for predicting droughts in various continents.A review by Dharani et al. [32] delves into the use of deep learning models for predicting crop yields, emphasizing classification and regression models tailored for this sector.Diaconu et al. [33] apply the ConvLSTM network architecture to predict NDVI and RGB values in satellite imagery for land cover classification.In soil science, Yadav et al. [34] investigate the effectiveness of traditional machine learning models like SVM and random forests in assessing soil fertility for croplands.Finally, Hounkpatin et al. [35] utilize classic machine learning techniques to study the chemical properties of soil in Benin.

E. FEATURE IMPORTANCE
Understanding which features significantly influence the predictions of machine learning models is crucial for both analyzing subjects using these models and evaluating model performance.Various methods exist for estimating feature importance, either on an object-specific basis or using statistical approaches.Sundararajan et al. [36] delve into the 'attribution' concept, which elucidates the significance of features for a specific object, introducing the Integrated Gradients (IG) approach.The attribution concept is rooted in two essential axioms: sensitivity and completeness.Authors in [37] explore different methods for determining feature importance across eight diverse environmental datasets, focusing on ensembles of machine learning models.They classify these methods into three main categories: filter, wrapper, and embedded techniques.Their findings highlight Permutation Importance and Shapley Additive Explanations (SHAP) as particularly effective, underscoring their adequacy in explaining feature relevance.Moreover, the authors emphasize permuation feature importance technique due to its universality and simplicity among machine learning models.

F. RESEARCH GAP
Despite comprehensive studies on climate and satellite data, along with machine learning applications in cropland analysis, gaps in research persist.One area with room for exploration is the use of machine learning and neural networks to differentiate and analyze croplands, specifically regarding irrigation practices in future climate scenarios projected by CMIP.Moreover, the interpretability of machine learning models in cropland analysis is still nascent, with limited methodologies for comprehensively understanding feature importance in model predictions.Lastly, the regional specificity in the behavior of these models under changing climate conditions remains insufficiently examined.Addressing these gaps is crucial for enhancing crop yield predictions and developing region-specific agricultural strategies in response to global climate change.

III. APPROACH
In this study, we develop a machine learning pipeline to interconnect overall climatic conditions with cropland classification by irrigation type.Our pipeline utilizes climatic and morphological data as predictive variables to estimate the state of agricultural land.The prediction is based on the variability of the present combinations of land use, climate, and morphological variations.We assume that the link between the possibility of agricultural land use and climate conditions is invariant, hence, we can forecast future land use based on climate information if sufficient present variability is captured as a learning sample.However, this approach requires a significant amount of data, encompassing all possible states of climate and agriculture, including the absence of agriculture, related to a possible change in climate conditions for the area of interest.
Our research specifically focuses on cropland types within Eastern Europe and Northern Asia.This region, which includes Eastern Europe and Northern Asia, provides a rich diversity for an agricultural land suitability analysis.We homogenize and align all data to the exact spatial and temporal resolution, then proceed with feature engineering and data processing (see Section III-A).Subsequently, we split the dataset into train and test subsets.
The next stage of our pipeline involves training and evaluation of the machine learning models.Finally, we forecast the distribution of cropland types based on various climate models, Shared Socioeconomic Pathways (SSP) scenarios, and analysis.The workflow also considers morphological data to constrain conditions where agricultural cultivation is impossible.
In summary, our methodology consists of a three-stage workflow: data acquisition and preprocessing, training of machine learning models, and result evaluation through the prediction of cropland distribution based on different climate models and SSP scenarios.This workflow yields robust results, leverages historical data, and is capable of simulating the distribution of cropland types under future climate projections.The detailed workflow methodology is illustrated in Figure 1.

A. DATA ACQUISITION AND PROCESSING
This subsection outlines the data sources and preprocessing methods utilized in our study.
Our study area spans from 20 • E to 152 • E and from 42 • N to 64 • N.This region encompasses a variety of climate zones, each with different levels of agricultural development.The datasets leveraged in our analysis are itemized in Table 1.
To identify cropland areas, we use the Global Food Security Analysis Data (GFSAD1km) which provides 1-km resolution data for the nominal year 2010.This dataset, compiled from remote sensing, secondary data, and field-plot information from 2009 to 2013 [11], originally categorizes croplands into five classes: major irrigation (class 1), minor irrigation (class 2), rainfed (class 3), rainfed with minor fragments (class 4), and rainfed with very minor fragments (class 5).The first three categories represent approximately 75% of the world's total cropland area, while the latter two make up about 20%.
In our study, we note that classes 4 and 5 have sparse cropland presence.Consequently, we reclassify them into Class 0, indicating areas with minimal or no cropland.This adjustment helps to focus on significant cropland types and accommodating potential shifts in land usability due to changing climatic conditions, where previously non-arable lands may become suitable for agriculture, or conversely, current croplands may lose their agricultural value.
The modified cropland mask now categorizes regions into one of four cropland classes, as demonstrated in Table 2.
For terrain data, we incorporate not only elevation but also ten morphological variables.These variables are calculated based on four hierarchical levels-3, 11, 33, and 47 kilometers-following the approach suggested by Turcotte [39].This methodology is designed to capture the landscape's variations across various spatial scales.Through spectral analysis, we decompose the absolute relief height into distinct frequencies and characterize each frequency's degree of expression through amplitude.Hierarchical levels correspond to the period.The fractal analysis reveals model deviations and generates morphometric characteristics listed in the second part of Table 3.These features describe elevation gradient, surface shape, and shaded relief, determining the risks of natural and anthropogenic activities while regulating heat and vapor distribution.
For historical climate data, our strategy involves collecting monthly climate statistics spanning 2000 to 2010.
To minimize the effects of random fluctuations or ''noise'' in these observations, we compute the ten-year average for each climate indicator.This process yields a set of 12 average monthly values per climate attribute.An exception is the Standardized Precipitation Index (SPI), for which we compute a single average value.We adopt the same methodology of data averaging for future climate projections for two time periods: 2022-2032 and 2040-2050.
The input for our models, referred to as X , is determined by the dimensions x • × y • and includes 121 climate-related features and 41 terrain features.The models aim to predict the cropland class label, which could be any of the four values: 0, 1, 2, or 3.
For the LSTM, ConvLSTM, and Transformer models, we specifically address the time series aspect of our  classification problem.We structure our input features into sequences spanning 12 months to reflect the annual cycle.Consequently, the shape of the input is expressed as (N , L, F), where N is the number of data samples, L is the sequence length (which for our study is 12 months), and F indicates the total number of features.In our case, out of 52 features, 42 are static, meaning that their values remain constant for each month.
The ConvLSTM model, initially developed for analyzing sequential multi-channel images, processes data through a convolutional layer before it is introduced to the memory cell.In our research, we apply the same ConvLSTM architecture but adapt it for one-dimensional data, ensuring the temporal sequence of monthly features is maintained and effectively utilized.
We split the historical data into training, validation, and test subsets, stratified by the target variable.The allocation is 80% for training, 10% for validation, and 10% for testing.For normalization, we use the Min-Max Scaler, a preprocessing method that scales features to a range between 0 and 1.We calculate the scaling parameters using the training set and apply these to the validation, test sets, and future data to prevent data leakage.
Given the significant class imbalance in the target variable, we adopt random undersampling of the dominant class in the training subset to achieve a more balanced class distribution.This approach helps to enhance model performance by preventing the overrepresentation of the predominant class from biasing the predictions.
Hyperparameter optimization is crucial for achieving optimal performance of machine learning algorithms.We optimize hyperparameters of logistic regression, CatBoost, XGBoost, and LightGBM by random search (RS) during a five-fold cross-validation (CV) process, with the test set reserved for final evaluation.F1-score serves as our primary performance metric.The ranges for the hyperparameters are provided in Appendix A.
We denote [P c M ] ij as a machine learning (ML) model prediction for a Class c for a pixel ij from the climate projection/reanalysis model M. We minimize the cross-entropy loss during training between the model output class probability and the current land use type from GFSAD1km.
Furthermore, we construct a meta-classifier (stacking ensemble) on top of all models.This second-level classifier is built based on the probability outputs of our primary models.The goal of the meta-classifier is to determine the best combination of the outputs of the base classifiers to produce the final prediction [48].This strategy synergizes insights from diverse models, yielding a unified probability for each class in every sample.In our study, CatBoost is selected as the meta-classifier model.Details on the inference procedure can be found in Section IV.

IV. PERFORMANCE EVALUATION AND DISCUSSION
Climatic shifts have the potential to transform regional landscapes, turning previously non-arable lands into croplands or vice versa.Our machine learning (ML) model, informed by climate data from 2000-2010, terrain information, and cropland masks from 2010 as a baseline and target variable, is engineered to forecast cropland types and their extents in Eastern Europe and Northern Asia.We operate under the assumption that terrain characteristics remain substantially unchanged.Following the methodologies described in [23] and [24], we compute an average of the climatic variables from three climate projections to form an ensemble that mitigates territorial-specific inaccuracies of individual models.
We examine the potential agricultural transformations under three Shared Socioeconomic Pathways (SSP) scenarios over two time intervals: the near future (2022-2032) and the midcentury (2040-2050).

A. CLASSIFICATION METRICS
When evaluating classification algorithms applied to imbalanced datasets, the selection of metrics that truly capture model performance is crucial.Our study employs a comprehensive set of metrics, including Accuracy, Precision, Recall, F1-score, and the mean Average Precision (AP).These metrics collectively provide a holistic view of model performance.
Recall, or the true positive rate, is critical as it quantifies the correctly identified positive instances out of all actual positives; a higher recall implies fewer false negatives for a given class.Precision, or the positive predictive value, reflects the accuracy of positive predictions -higher precision means fewer false positives.The F1-score serves as the harmonic mean of precision and recall, providing a single score that balances the two-a higher F1-score indicates better precision and recall balance.The AP, particularly critical in multi-class classification contexts, assesses whether a model can accurately distinguish positive from negative examples across different decision thresholds for each class.It does this by averaging the area under the precision-recall curve for each class, yielding a singular score.
Employing macro averages for these metrics ensures a balanced assessment across variably distributed classes, mitigating the influence of class overrepresentation.As indicated in Table 4, although no single model consistently outperforms others across all metrics, CatBoost and LSTM show slightly better overall performance.
Additionally, the meta-classifier, which combines all individual models, exhibits a modest but consistent enhancement in performance metrics.This improvement is attributable to the meta-classifier's integrated approach, which leverages the distinct advantages of each individual model, capturing a broader array of data patterns and thereby enhancing overall predictive accuracy, as discussed in [49].Further insights into the meta-classifier's performance across different classes are presented in Table 5.Here we observe that the precision for classes 0 and 1 is commendable, with class 0 achieving the highest precision at 0.99.However, recall for class 1 is comparatively lower, suggesting challenges in consistently identifying true positives within this category.The F1-score is the highest for class 0, with a value of 0.99.The average precision score, measuring the area under the precision-recall curve, is high for all classes, signifying the model's effectiveness in ranking positive instances higher than negative ones.

B. PROBABILISTIC HEATMAPS
In this subsection, our objective is to explore how cropland suitability varies across different locations by the year 2050.To this end, we use a machine learning model trained on land use patterns from 2010 to predict the probabilities of cropland types at a granular level -pixel by pixel.
Given that the meta-classifier provides a probability output P c ij for class c at each pixel location ij, we assemble these probabilities into a comprehensive heatmap for class c, denoted by P c = ∪ ij P c ij .We then extract the predicted class probabilities P c ij for each pixel, constructing probabilistic heatmaps P c .To contextualize changes, we compare these predictions against the 2010 baseline land use for each class, symbolized by GFSAD1km c , where GFSAD1km c ij ∈ {0, 1}.The resulting differential heatmap (P c ) = P c − GFSAD1km c reveals deviations from the baseline.
A differential value close to zero (0) indicates concordance with the historical state of the land.Conversely, a value approaching +1 suggests an increased probability of class c emerging in that locale, while a value near −1 hints at a potential decline or replacement of the current class.For example, within non-cropland areas (class 0), a difference nearing −1 could signify a transition to agricultural land suitability.
Our analysis, through the probabilistic heatmaps (P c ) visualized in Figure 2, particularly under the SSP2-4.5 scenario, points to several implications: regions traditionally used for agriculture are at risk, while previously non-arable northern areas are trending towards suitability for cultivation, even without irrigation enhancements.Furthermore, the Southern Urals region exhibits a pronounced requirement for increased irrigation to compensate for the diminishing rainfed croplands.These findings emphasize the critical need for sustainable agricultural and land-use strategies that are attuned to the ecological capacities and changing conditions of the local environments.Additional heatmaps for the SSP1-2.6 and SSP5-8.5 scenarios are presented in Appendix A, in Figures 7 and 8.The SSP1-2.6 scenario indicates less severe and more gradual changes in cropland types.In contrast, the SSP5-8.5 scenario, which represents the most adverse socio-economic conditions, shows a more stark transformation.
Furthermore, we examine how the number of pixels allocated to each class of croplands changes under various SSP scenarios until the year 2050 and present the distribution of cropland suitability in Figure 3.Our analysis shows that there is a general trend toward increasing the suitability of lands for cropland use across all SSP scenarios.The SSP5-8.5 scenario displays a positive trend despite being environmentally unfriendly.Our findings are consistent with previous studies, including [50].The figures suggest a decrease in the demand for major irrigated lands (class 1), an increase in rainfed (class 3), and a significant relocation and growth of minor irrigated lands (class 2).While this last class may indicate further potential for development, more attention should be given to the lands that transit from rainfed to irrigated, as they are already productive and are becoming more sensitive to water access.The overall trend is similar across all considered SSPs.

C. PERMUTATION FEATURE IMPORTANCE
To enhance the interpretability of our research, we employ the permutation importance algorithm with CatBoost, MLP, and LSTM-our top three models based on the recall score.This model-agnostic method assesses feature importance by 15754 VOLUME 12, 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.We employ the F1-score as the baseline performance metric s.For each feature j, we shuffle the values and compute the perturbed model score s k,j across k iterations.The importance I j for each feature f j is calculated as the mean difference from the baseline: This permutation process is independently replicated 10 times for each feature.Given the sequential nature of input data for the LSTM model, its feature grouping differs; a single importance value is assigned for each climate feature, while for other models, importance is aggregated over 12 monthly values.The importance of terrain features is obtained as their mean value in four hierarchical levels.Subsequently, we normalize these importance scores for each model to adjust for scale differences.We then compute the average of these normalized scores for each feature across all models.Lastly, features are ranked based on their mean importance scores.
This method of assessing feature importance discloses the relationship between the feature and target variable.Consequently, the rate of decrease in the score for a specific feature can serve as an indicator of its importance to the model.
Our analysis, as visualized in Figure 4, uncovers that key indicators include monthly precipitation, along with various temperature-related metrics such as the mean monthly air temperature, the frequencies of days experiencing: temperature rises over 6 • C and temperature transitions through 0 • C. Interestingly, model-specific preferences emerge: MLP prioritizes temperature fluctuations over 6 • C, while both LSTM and CatBoost highlight monthly precipitation as paramount, with the latter showing a marked preference.The frequency of days with temperature rises greater than 6 • C proves particularly significant in northern areas, where agricultural expansion is likely under current climatic trends.
The monthly precipitation feature is anticipated to influence the classification related to irrigation strategy directly.

D. REGION SPECIFICS AND LOCAL FEATURE ATTRIBUTES
In this subsection, we analyze the areas classified as high-risk or potential for cultivation by our meta-model.These regions -Northern China, Southern Urals, and Eastern Europe -are delineated with coordinates in Table 6 and depicted on the probability heatmaps in Figure 2 for predictions of Classes 1, 2, and 3, respectively.
An overview of monthly precipitation, the top influential feature identified in the permutation importance analysis, is provided in Figure 5.It exhibits considerable variability among specific CMIP projections, justifying the rationale for their averaging.On the other hand, the models and historical data differ regionally.We tackle the areas that are likely to experience serious changes in cropland type, as discussed in Section IV-B.Northern China is projected to transform currently major-irrigated lands, while Southern Urals is expected to require minor irrigation and lose its suitability as a rainfed cropland, moving from Class 3 to Class 2. On the other hand, Eastern Europe is expected to become suitable for arable use under rainfed irrigation -Class 3.
Next, for these regions, we aim to reveal specific features and their influence on the cropland type change using the Integrated Gradients (IG) method [36].This method allows to attribute features and their influence on the neural networks prediction.According to the original paper, if F(x) is a neural network with input x = (x 1 , . . ., x n ), where x i is an individual feature, then A F (x, x ′ ) = (a 1 , . . ., a n ) is called an attribution of the network prediction on x relative to a baseline x ′ if a i describe a contribution of x i to the prediction F(x).This method satisfies completeness and sensitivity axioms.
Essentially, the IG attribution for a network F is defined as In other words, it is an integral of the gradient component corresponding to the feature x i , between the baseline and the current input.For a multi-class problem, we consider F from above as a single logit of the class of interest.For example, if the goal is to obtain attributions that impact either a shift to Class 3 (rainfed irrigation cropland) or away from it, we choose F to be the output of our model for this class.
To interpret the models using the IG attribution method, we estimate the attributions for the top 2 neural network based models (LSTM and MLP) by averaging over the corresponding territories depicted in Figure 2. We choose historical climate data as a baseline (x ′ ) and the 2040-2050 climate projection averaged over the SSP2-4.5 scenarios as an input (x).The estimation results are presented in Figure 6.The key observation is that for each region and each model, there is a general consensus on features semantically.Indeed, the attributions provided by IG consider precipitation and temperature-related features as the most important.However,  some feature attributions contradict each other for different neural networks.We discuss this issue on region-specific scope below.6.
In Liaoning, Northern China, Integrated Gradients (IG) attribution analysis reveals that an increase in snow moisture content (snw), and total monthly precipitation (tp) potentially maintains the prevalence of predominantly irrigated croplands (Class 1), as indicated by positive attribution values.

TABLE 7.
Hyperparameter grids for various models with descriptions.The parameters that are not listed here are assumed to be default as in the corresponding packages.
Conversely, the of days with extreme precipitation (pr_p95) exhibits a considerable negative attribution, gesting that significant shifts in precipitation patterns may lead to changes in the irrigation type.Furthermore, Figure 2 illustrates a pronounced shift toward rainfed croplands, implying that changes in precipitation also influence the type of cropland class.Additionally, a decline in monthly average temperature is likely to conserve the dominance of irrigated croplands, as inferred from the negative attribution scores for both models.Notably, discrepancies in feature attributions between the LSTM and MLP models, such as the number of days with extreme temperature (tasmax), and the 12-month Standardized Precipitation Index (12m_SPI), suggest that model-specific focus on different classes and small regions may result in inconsistent local feature significance.
In the Southern Urals, similar to Northern China, a rise in snow moisture content (snw) correlates with an expansion of minor irrigated croplands.The lack of a consensus on the influence of other features is attributed to the dynamic shifts in land classifications within the region.From Figure 2, one can observe that there is also a significant increase for Class 3, let alone Class 2. Similarly, as in China, the models appear to prioritize different locales and cropland classes.
In Eastern Europe, an increase in total monthly precipitation (tp) and SPI (12m_SPI) tends to facilitate the spread of rainfed croplands, along with an increase in the average temperature at 2m (t2m).The inconsistency between attributions for snow moisture content (snw) and the other features might be caused by the fact that models performed differently across classes.As one can observe from the map in Figure 2, there are also transitions to minor irrigated croplands, which affected attribution.
In summary, our approach provides valuable insights into the local feature importance of developed models, VOLUME 12, 2024 15759 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
2.70GHz, and Quadro RTX 8000 GPU.Training of each model has taken not more than 8 hours.

V. CONSTRAINTS OF THIS STUDY
This research, while providing valuable insights into agricultural land suitability under climate change scenarios, is subject to several limitations inherent in the data sources and methodologies used.The GFSAD1km dataset, while invaluable for identifying cropland classes, is constrained by its spatial resolution and classification accuracy, making it a broad indicator rather than a definitive source of land usage patterns.
Similarly, the ERA5 climate reanalysis model, despite its extensive temporal resolution and coverage, has its limitations in capturing specific meteorological phenomena.Such limitations can result in misrepresentations of climatic variables, possibly leading to over-or underestimations of real-world conditions.
Furthermore, the CMIP6 projections, which form the basis of our future climate scenarios, contain inherent uncertainties and biases.These variabilities introduce a degree of unpredictability in our long-term projections, affecting the accuracy of our study's foresight regarding cropland suitability.As a result, our model predictions should be interpreted with caution, especially in scenarios where localized agricultural practices are significant.
These constraints suggest a need for cautious interpretation and application of our study's findings in policy and planning.Future research should focus on developing higher-resolution datasets, exploring advanced modeling techniques, and integrating multiple data sources to mitigate these limitations and enhance the reliability of studies in this domain.Our findings, while indicative, should be complemented with additional data and analyses to support robust and informed agricultural and environmental strategies in the face of climate change.

VI. CONCLUSION
In this study, we provide a comprehensive evaluation of future cropland suitability in the context of varying climate projections up to the year 2050, leveraging a suite of advanced machine learning algorithms.Our findings offer a robust foundation for predicting land suitability changes.Importantly, this study aligns with and complements recommendations from the Intergovernmental Panel on Climate Change (IPCC) [53], emphasizing the need for detailed regional assessments in adapting to climate variability and securing food supplies.By mapping out potential shifts in cropland viability and irrigation needs, our work offers practical insights for policymakers, crucial for effective risk management and adaptation strategies as underscored in the IPCC report, particularly in the context of global efforts to combat hunger and malnutrition.
Monthly precipitation and a set of temperature-related features were found to be critically important for classification.Temperature and morphological variations exhibited a strong relationship with croplands differentiated by the irrigation type.Our localized feature importance analysis, especially for at-risk regions in China, underscores the significance of temperature extremes and precipitation patterns.This highlights the need for region-specific interventions, such as modifying irrigation practices or reallocating agricultural resources, to mitigate the risks posed by climate change.
The total amount of land suitable for agriculture is projected to increase, expanding further to the north.However, some currently exploited agricultural regions may require increased irrigation, posing potential risks.This expansion indicates significant potential for land development and the possibility of meeting the growing food demand.To harness this potential, intelligent land management and substantial investments are essential.Neglecting the impacts of climate change may lead to considerable cropland shrinkage, jeopardizing food security.
Our code is available online on GitHub1 and allows replicating the analysis presented in the paper.

APPENDIX A SUPPLEMENTARY DATA
See Table 7.

ACKNOWLEDGMENT
Author Contributions: Valeriy Shevchenko: analytics, data acquisition, software development, and text writing; Aleksandr Lukashevich: analytics and algorithms, and text writing; Daria Taniushkina: data analysis and text writing; Alexander Bulkin: paper revision; Roland Grinis: paper revision; Kirill Kovalev: critical review; Veronika Narozhnaia: paper revision; Nazar Sotiriadi: conceptualized the study; Alexander Krenke: data interpretation and scientific supervision; and Yury Maximov: paper revision.

FIGURE 2 .
FIGURE 2. Heatmaps of class probabilities from the Ensemble model for arable lands classes for 2050 under the SSP2-4.5 scenario.Blue rectangles stand for areas that are studied in Subsection IV-D.

FIGURE 3 .
FIGURE 3. Pixel-wise cropland classification results forecasted by the meta-classifier under different SSP scenarios, with GFSAD1km cropland mask as the initial point.

FIGURE 4 .TABLE 6 .
FIGURE 4. Models (permutation) top features importance.TABLE 6. Regions with significant class change.Coordinates are given in latitude-longitude order.

FIGURE 6 .
FIGURE 6. Top 10 local specific attributes obtained from IG, zones are defined in Table6.

FIGURE 7 .
FIGURE 7. Heatmaps of the class probabilities from the Ensemble model for arable lands classes for 2050 under the SSP1-2.6 scenario.Blue rectangles stand for areas that are studied in Subsection IV-D.

FIGURE 8 .
FIGURE 8. Heatmaps of class probabilities from the Ensemble model for arable lands classes for 2050 under the SSP5-8.5 scenario.Blue rectangles stand for areas that are studied in Subsection IV-D.

TABLE 1 .
The datasets description.

TABLE 4 .
Comparative performance assessment of ML Models for arable land prediction.

TABLE 5 .
Per class performance of the Meta-classifier.