The current research uses an integrated geoinformation approach, whereby we applied a fuzzy-object-based image analysis (FOBIA) and spectral analysis to obtain a time series of LUCC and soil salinization monitoring. We used the SAR satellite images to detect the land subsidence from 2015 to 2020. The GIScience spatiotemporal analysis was also applied for trend analysis of climate indicators, soil characteristics, and water resources to develop the sustainable food production map. The central aspect of the research methodology is explained in the following:

## Land Use/cover Mapping And Trend Assessment

We processed Landsat time series satellite images for LUCC and soil salinization monitoring. To detect changes, we obtained the satellite images for July of 1990, 1995, 2000, 2005, 2010, 2015 and 2020. To obtain the most accurate results, we applied the integrated FOBIA approach to LUCC monitoring and mapping. Therefore, we applied a multi-resolution segmentation to obtain image objects. For a FOBIA classification, segmentation parameters are critical since they directly impact the size of image objects and the final classification results. Therefore, we estimated the scale parameter 31 and performed the segmentation with the scale of 20, a shape index of 0.6, and a compactness value of 0.4. We identified the primary land use/cover pattern of the LUB based on a comprehensive discussion with the authorities and decision-makers and through field work. We employed a FOBIA classification to derive the characteristics of the LUCC subclasses. The training data were collocated in field operations and from existing and historical LUCC and cadaster maps and aerial photography, which were available in both of the agricultural resource organizations of West- and East Azerbaijan Provinces (ARO-W/EAP). We used a total of 21000 points (3000 points for each study year) to identify the characteristics of each LUCC subclass and as training data for the FOBIA classification. We were able to identify the relevant object features based on the spectral and spatial characteristics of each LUCC subclass and evaluated their effectiveness for driving each LULC subclass by comparing the training data with the fuzzy membership value of each object feature based on the following equations:

o Spectral attributes/ Brightness

\(\varvec{B}=\frac{1}{{\text{n}}_{\text{v}\text{i}\text{s}}} \sum _{\text{i}=1}^{{\text{n}}_{\text{v}\text{i}\text{s}}}{\stackrel{-}{\text{c}}}_{\text{i}\left(\text{v}\text{i}\text{s}\right)}\) (1)

Where *B* is the mean brightness of an object and\({\stackrel{-}{\text{C}}}_{\text{i}\left(\text{v}\text{i}\text{s}\right)}\) is the sum of all the mean brightnesses in the visible bands divided by the corresponding number of bands \({\text{n}}_{\text{v}\text{i}\text{s}}\)

o Normalized Difference Vegetation Index (NDVI)

Tv = mean NDVI

f (Object)= \(\left\{\begin{array}{c} LC if f\left(\text{O}\text{b}\text{j}\text{e}\text{c}\text{t}\right)\le {\text{T}}_{\text{v}})\\ VA if f\left(Object\right)>{\text{T}}_{\text{v}})\end{array}\right.\) \({\text{T}}_{\text{v}}^{{\prime }}=\frac{{\text{m}\text{e}\text{a}\text{n}\text{N}\text{D}\text{V}\text{I}}_{\text{L}\text{C} }+{\text{m}\text{e}\text{a}\text{n}\text{N}\text{D}\text{V}\text{I}}_{\text{V}\text{A}}}{2}\) (2)

\({T}_{v}^{{\prime }}\) is an average of the mean NDVI values for the candidate object of LUCC class and vegetated areas (VA). The NDVI, which has a value between − 1.0 and + 1.0

Green Normalized Difference Vegetation Index (GNDVI)

GNDVI = 100*(1+(([Mean band 4]-[Mean band 2])/([Mean band 4]+[Mean Layer 2]))) (3)

o Modified Normalized Differenced Water Index (MNDWI)

MNDWI = 100*(1+(([Mean Band4]-[Mean Band5])/([Mean Band4] +[Mean Band5]))) (4)

o Normalized Built-up Index (NDBI)

NDBI= ([Mean band 5]-[Mean band 4])/([Mean band 5]+[Mean band 4]) (5)

o Soil Water Content Index (InfraRed Index)/ Soil Color Index

SWCI (IR)= (NIR- ETM7)/ (NIR + ETM7) (6)

SCI = R – G / R + G

o Normalized Built-up Index (NDBI)

NDBI= ([Mean band 5]-[Mean band 4])/([Mean band 5]+[Mean band 4]) (7)

Salinity Index (SI)

o Normalized Difference Salinity Index

Salinity Index (SI)= ([Mean Layer 1]*[Mean Layer 3])^0.5 (8)

NDSI=([Mean band 3]-[Mean band 4])/([Mean band 3]+[Mean band 4]) (9)

o Shape Rectangular

\(1-\sqrt{\frac{{{\lambda }}_{\text{m}\text{i}\text{n}}}{{{\lambda }}_{\text{m}\text{a}\text{x}}}}\) (10)

\({{\lambda }}_{\text{m}\text{i}\text{n}}\) is the minimal eigenvalue

\({{\lambda }}_{\text{m}\text{a}\text{x}}\) is the maximal eigenvalue

\(\frac{\{\#\left(x,y\right)\epsilon {P}_{v}:\rho v(x,y\rho )\le 1\}}{\#{P}_{v}}{P}_{v}=\)(x,y) is the elliptic distance at a pixel (x,y) (11)

[0,1]; where 1 is a perfect rectangle

o Shape Index

\(\frac{{\text{B}}_{\text{v}}}{{4\sqrt{{\#}_{\text{V}}}}_{\text{v}}}\) (12)

\({\text{B}}_{\text{v}}\) is the image object border length

\({4\sqrt{{\#}_{\text{p}}}}_{\text{v}}\) is the border of square with area of \(\#{P}_{v}\) (13)

After identifying the relevant object-based features for LUCC classification, we exported the computed membership values for all objects in LUB to Python for the deep learning classification and to produce the final LUCC map. To validate the results, we applied an accuracy assessment based on the fuzzy synthetic evaluation (FSE), which was proposed by Feizizadeh 32 and acknowledged for its effectiveness in FOBIA classifications by several studies 3, 33. Technically, the FSE makes use of two sets of data, namely reference data (ground truth data) and the results of the OBIA-based classification map. The FSE is used to compute the interpretation of confidence ratings (ICR), which is the classification confidence and observed level of error for each class 33 (See supplementary Table 2). Therefore, 30% of all ground control points (6300 out of 21000) were used as references data, and the overall accuracy for the study period of 1990–2020 with 5-year intervals was computed to be 0.94, 0.93, 0.98, 0.95, 0.96, 0.93 and 0.97, respectively (See supplementary Table 3–4 for results of accuracy assessment based on the FSE).

## Land Degradation Monitoring (Soil Salinity And Land Subsidence)

We used the obtained Landsat images to carry out a LUCC classification to monitor and map soil salinization. According to earlier studies, the spectral properties of Landsat allow us to compute the soil salinization rate based on the soil salinity indices 34–35. The soil characteristics are impacted by various factors (e.g., vegetation cover, moister, texture, parent material, etc.), which must be considered in soil salinity mapping 36–37. Therefore, we applied the Combined Spectral Response Index (CSRI) technique to compute the soil salinity ratio for each study year. The CSRI method was developed by Fernandez-Buces et. al 38 and has been confirmed to be an efficient method for determining soil salinity by the remote sensing community 37,39−40. The CSRI is based on the spectral information of RGB and the near-infrared bands, and it is computed as follows:

\(\text{C}\text{S}\text{R}\text{I}=\frac{\text{B}+\text{G}}{\text{R}+\text{N}\text{I}\text{R}}\times \text{N}\text{D}\text{V}\text{I}\) and (B+G)/(R+NIR) × NDVI (14)

NIR: Near infrared and NDVI: Normalized Differences Vegetation Index

We used the field observation data together with the soil electrical conductivity (SEC) data to validate the results. Time series SEC data were obtained from the ARO-W/EAPs. In order to validate the results, a linear regression analysis was applied to compute the spatial correlation between the reference data and the results of the CSRI method for each study year. The result of the linear regression analysis from 1990–2020 with 5-year intervals was computed to be 0.93, 0.92, 0.90, 0.92, 0.96, 0.97 and 0.95, respectively, which represents a very high spatial correlation between the reference data and the obtained soil salinization maps.

As indicated in the context of land degradation, we also aimed to monitor the land subsidence in the LUB. We applied differential interferometric synthetic aperture radar analysis (DInSAR) for the 2015–2020 period to create a land subsidence map. The DInSAR technique is a reliable and fast approach to derive long- and short-term deformations. The technique mainly relies on the ‘master’ and ‘slave’ SAR image processing of the same part of the earth from the same satellite orbit. In the repeated pass, the difference of phase value correlation of two SAR images can be used to estimate the ground subsidence. Generally, the phase correlation or coherence () of two acquisitions in DInSAR analysis contains various sources. After reducing non-deformation phases, such as the atmospheric phase, topographic phase, and noises, the multi-temporal DInSAR can be implemented. In the multi-temporal analysis, images are taken at different times (T1, T2, …, Tn), and the phase value of the image is a simple way to measure the changes in satellite-to-target direction (R1, R2, …,Rn) across the same study area 41−43. For the DInSAR analysis, we used 6 single look complex (SLC) images of descending orbits (T 79) of Sentinel-1 to cover the study area. The images were taken on the 2015.11.11, 2017.11.12, and 2020.11.02, and cover the upper and lower parts of the Lake Urmia basin (6 images in total) in interferometric wide (IW) mode. The images are in VV (vertical - vertical) and VH (Vertical - horizontal) polarizations, and the incidence angle of the images are 39.08°. The incidence angle difference of the first pair (2015.11.11 and 2017.11.12) and the second pair (2017.11.12 and 2020.11.02) are 0.009° and 0.01°, respectively. Due to the semi-arid to arid climate in the study area, the large temporal baseline (~ 2 years) is not a serious problem for the interferometric analysis, which has also been acknowledged in similar studies 21−24.

## Water Salinization Spatiotemporal Analysis

The significance of freshwater for agricultural-, industrial-, and residential purposes is beyond debate. Due to the semi-arid climate in the LUB, agricultural activities were traditionally based on surface run-off water and groundwater. The LUCC results reveal a substantial increase in croplands and farmland over the past 30 years, which predominantly depend on groundwater. The regional water organizations in both West and East Azerbaijan Provinces (RWO-E/WAPs) have monitored the LUB water quality parameters annually since 2000 and quarterly since. We made use of observations from deep- and semi-deep wells, springs, and Qanats during the wet and rainy seasons. Since the number of observation wells has increased every year since the Lake Urmia drought, the groundwater observation data used in the spatiotemporal modelling of the groundwater quality assessment was limited to 630 in 2005, 684 in 2010, 751 in 2015, and 859 in 2020. We obtained and analyzed the following water quality parameters for the studied time period (2000–2020): electrical conductivity (EC), power of hydrogen (pH), potassium (K+), total hardness (TH), sulfate (SO42-), magnesium (Mg2+), chloride (Cl-), sodium (Na+), total dissolved solids (TDS), calcium (Ca2+), carbonate (CO32-), and bicarbonate (HCO3-). We obtained groundwater physicochemical data and evaluated the water status using the water quality index (WQI), spearman correlation, principal component analysis (PCA), and agglomerative hierarchical clustering (ACA) based on the GIS spatial analysis for determining the hydro-geochemical attributes. The statistical evaluation of the physicochemical groundwater parameters is essential to comprehend the main factors controlling water quality variations over time 44.

Multivariate statistical techniques such as the PCA are effective approaches for data visualizing and mining. A PCA technique was used to identify correlations between physical and chemical characteristics in the aquifer and to elucidate complicated patterns in data matrices 45–46. Squared cosines with absolute values of more than or equal to 0.4 were used for the exploration of the observed ions 46–47. We used the Spearman correlation coefficient to identify the sources of different elements in the groundwater samples. High coefficients reveal the significance of the relationship between two parameters. A positive coefficient indicates that the associated parameters are similar and harmonious, while a negative coefficient indicates that the variables are moving in opposite directions. We carried out an agglomerative hierarchical cluster analysis (ACA) of the groundwater data to identify a specific pattern of similar observations within the studied variables 46. The ACA allowed us to identify suitable structures in chaotic and complex data and thus simplify the explanation of observations 48. Using a GIS spatial analysis, we obtained the distribution characteristics and created the groundwater quality maps. While many interpolation methods are available, we choose to use the inverse distance weighted (IDW) interpolation technique, which is the most frequently used deterministic modelling tool 49. The model is based on the hypothesis that a node shares more similarities with nearby points and that it is influenced more strongly by the surrounding data values. We used 856 points and derived the node values by averaging the weighted total of all of them.

## Environmentally sustainable food production mapping

We applied a GIS-based multi criteria decision analysis (GIS-MCDA) to multiple factors of food production in the LUB to examine sustainable food production under the impact of the Lake Urmia drought. Sustainable food production in a fragile ecosystem such as the LUB must consider the interaction of several causal spatial indicators. The GIS-MCDA provides an effective approach for dealing with the spatial decision problems and patterns based on the concept of ‘spatial thinking’ 50–51. The GIS-MCDA allows us to consider the relevant spatial indicators and their characteristics as attributes in combination with decision maker’s preferences to analyze the spatial problems 51. In this research, we analyzed sustainable food production in the LUB by considering the relevant causal climate indicators affecting agricultural production, namely annual average precipitation, temperature, sunshine hours, and humidity. In the context of soil characteristics, the soil degradation, fertility, texture, and depth were also used as causal indicators. From the land characteristics perspective, the salt scattering spots and land use/cover were also considered as relevant indicators for SFP (See Supplementary Fig. 6). The selection of these causal indicators is based on expert opinions, data availability, and relevant research literature 29,52−53.

After finalizing the causal environmental indicators, the initial data were collected from the relevant governmental departments, satellite images, and the Spatial Data Infrastructure (SDI) of LUB. All data were gathered, relevant geometric edits were applied, and the indicators were developed into a spatial GIS dataset. Since the spatial GIS data were collected from heterogeneous resources with different scales, we used the fuzzy standardization process to standardize the indicators in the same scale suited for criterion weighting. This approach considers the degree of fuzzy membership values on a scale of 0–1 based on the benefit/ cost context of the indicators 49. In the GIS-MCDA analysis, the significance of each criterion is computed as part of the decision analysis. We employed the Fuzzy Analytical Network Analysis (FANP) to determine the significance of each indicator. The FANP is an efficient GIS-MCDA weighting tool 50. We included the knowledge of 35 experts from different departments of the University of Tabriz and the University of Urmia. We asked these experts from agricultural-, natural resource-, and food security departments to initially rank the provided indicators. Then, based on the FANP method, we yielded the following criteria weights: precipitation (0.098), temperature (0.096), sunshine hours (0.075), humidity (0.058), groundwater depth (0.071), water quality (0.095), soil degradation (0.069), soil fertility (0.098), soil texture (0.099), soil depth (0.098), salt scattering spots (0.086), and land use/cover (0.098). For the criteria weighting, it is necessary to compute the consistency ratio (CR) to validate the obtained weights. According to Saaty 54, a CR < 0.1 indicates an acceptable level of consistency among the experts and that the results can be employed for a spatial aggregation. Our study yielded an acceptable CR value of 0.065.

Criteria weighting in GIS-MCDA significantly influences the uncertainty and reliability of the results 55. Therefore, we applied a global sensitivity analysis (GSA) to determine the validity of the computed weights using the FANP method. The GSA approach is used to compute the two critical indexes of *S* (first order) and *St* (total effect). The letters *S* and *St* denote that the FANP's weights were assigned in a semantic manner. Any discrepancy in the value and order of the S and St indices, as well as the reference weights (e.g. FANP's weights), can be regarded as uncertainty connected with the criteria weights 56. Accordingly, we used a ‘weighted overlay’ spatial aggregation function to produce the food production capability map.

## Environmental Prediction Based On The Ca-markov

We used a combination of a Markov model and a Cellular Automata model (CA-Markov) to predict the land degradation and water salinity based on the trend observed in the past years (1990–2020). Both techniques have been identified as being effective techniques for GIScience prediction problems 57. A typical CA-Markov model separates the discrete cellular, finite state, neighbor, and rule features into four categories while analyzing the trend as follows:

$${Z}_{(I+1)}={Z}_{\left(I\right)}\times Q$$

15

Z(I) = the first map in year *A l, Z(l + 1)* = the second map year *B, l + 1, Q* = state transition matrix, and Z can be described as the following matrix:

Z = [Z1 Z2 Z3] (16)

*Zi (i = 1, 2, 3)* represents the changes between classes of map *A* and *B, Q* can be described as an *[n, n]* matrix in the following, where *n =* total number of changes between map *A* and *B, Qij =* transition probability for any changes between map *A* and *B* from the time line of *i* to *j*, and the sum of each row of the matrix should be equal to 1.

$$Q=\left[\begin{array}{cccc}{Q}_{11}& {Q}_{12}& \dots & {Q}_{1n}\\ {Q}_{21}& {Q}_{22}& \dots & {Q}_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ {Q}_{n1}& {Q}_{n2}& \dots & {Q}_{nn}\end{array}\right], \sum _{(j=1)}^{n}{Q}_{ij}=1, 1 \le i, j \le n, n=4$$

17

According to a transformation function, the next state cell is decided by the current state and its surroundings 57–58. This method is based on generating a transition probability matrix between two maps in different timelines. The transition probability matrix enables an assessment that indicates the likelihood of each pixel in class A of the first map class converting to another class (e.g. B, C, D, ...) or remaining in class A in the second map class 59. The Markov chain, which is technically a separate random process, uses transition probability to forecast the next state and all future states based on the current state 60. Except for some of the after-effect occurrences, it is a viable method for predicting regional characteristics. We used CA-Markov to predict the land use/cover, soil, and water salinization maps for 2030, 2040, and 2050 based on the trend obtained for each environmental indicator.

## Spatial Uncertainty Analysis

Understanding the spatial uncertainty is critical in geospatial analysis and modelling. Data quality, correctness, inaccuracy, vagueness, fuzziness, and imprecision, are all referred to as spatial uncertainty in GIScience 32. Uncertainty in GIS-based modelling is inevitable due to the variety caused by a heterogeneous dataset, expert opinions, and model error, which has forced the GIScience community to develop spatial and statistical approaches to understand and quantify uncertainty 60. The Dempster Shafer Theory (DST) has been identified as an efficient technique for spatial uncertainty analysis. It is based on mathematical operations employing Bayesian probability theory 3. The DST's features make it an effective tool for modeling the imprecision and computing the spatial uncertainty 60. The DST can compute the epistemic uncertainty that affects expert knowledge of the probability P(M_) within the alternative model M_, = 1,...,n. This is also known as the ‘theory of evidence’, which aims to compute the BBA m (A) on sets A (the focal sets) of the power set P(Z) of the event space Z, i.e., on sets of outcomes rather than single elementary events 60. The belief function computes the lower limit value for a (known) probability to determine the spatial uncertainty. The plausibility function additionally predicts the upper bound value for an (unknown) probability. The geographic uncertainty can be calculated using the difference between the plausibility (*Pl*) and belief (*Bel*) functions. We employed the *Belief* function in Idrisi software to calculate DST and to define the spatial uncertainty of the predicted soil and water salinization maps for 2030, 2040, 2050 as well as the computed sustainable food production map (See Supplementary Fig. 7).

$$m: {P\left(Z\right)}^{ }\to \left[\text{0,1}\right], m\left(0\right)=0;\sum _{A∊\left(\text{Z}\right)}^{ }m\left(A\right)=1 \left(18\right)$$