Assessing a multi-method approach for dryland soil salinization with respect to climate change and global warming – The example of the Bajestan region (NE Iran)

strongly affects soil properties, with severe consequences for regional ecology, agri-culture and the aeolian dust dynamics. Given its climate-sensitivity it forms a serious environmental hazard


Introduction
Drylands (aridity index < 0.65) cover > 40% of the global land surface today, and host > 38% of the world's population (UNEP-WCMC, 2007;Huang et al., 2017).Intensified greenhouse gas emissions, leading to global warming and climate change, augmented the frequency and severity of extreme climate events and environmental challenges worldwide.This holds especially true for ecologically sensitive drylands, making them especially vulnerable towards these changes (Kweku et al., 2018;Aalijahan et al., 2019;Aalijahan et al., 2022;Nama et al., 2022).Hence, drylands will increase and often become drier during the near future (Feng and Fu, 2013;Huang et al., 2017;Cicek et al., 2022), what will significantly change their hydro-geomorphic and ecological processes with severe consequences for local societies (Ravi et al., 2010;Burrell et al., 2020;Khosravichenar et al., 2020).One fast environmental process in drylands is the soil salinity dynamics, i.e. soil salinization due to surficial accumulation of water-soluble salts.This is caused by an annual rainfall deficit that hinders salt leaching into the subsoil, and strongly affects the physical and chemical soil properties leading to severely changing ecological site conditions (Stockle, 2002;Metternicht and Zinck, 2008).As a consequence, agricultural land use in affected areas is complicated or even rendered impossible (Rengasamy, 2006).Furthermore, such areas are sources of saline aeolian dust affecting larger downwind areas (Ge et al., 2019, Zucca et al., 2021).Hence, soil salinization forms a serious global environmental hazard in the context of the current global climate change (Dehhan and Taylor, 2002;Metternicht and Zinck, 2003;Nawar et al, 2015).Therefore, to plan effective strategies in combatting dryland degradation and to forecast the aeolian dust dynamics, knowledge about regional soil salinization is necessary (Hassani et al., 2020).However, former multi-annual soil salinity studies generally used large and irregular time steps or only short sample periods (e.g.Wu et al., 2014;Scudiero et al., 2015;Dasgupta et al., 2015;Bannari and Al-Ali, 2020;Hassani et al., 2020;Rafik et al., 2022;Ma et al., 2023), hindering a deeper understanding of this environmental process.Furthermore, apart from singular studies suffering from the above-mentioned shortfalls (Bannari and Al-Ali, 2020;Rafik et al., 2022), regional soil salinization was not compared with global climate change.
Playas are intracontinental endorheic dryland basins with a negative water balance for more than half of the year, being characterized by evaporitic sediments.They often contain lakes that are neither dry nor wet > 75% of the time, showing seasonally and interannually strongly varying lake levels (Briere, 2000).Playas deliver ecosystem services such as groundwater recharge, surface water storage and wetland A. Khosravichenar et al. habitats (Bowen and Johnson, 2017).Given the constant accumulation of soluble salts delivered by their tributaries, playas form hotspots of soil salinization (Metternicht and Zinck, 2008).Furthermore, seasonally dry playa surfaces are one main source of saline aeolian dust (Ge et al., 2019, Zucca et al., 2021, Khosravichenar et al., 2021).Caused by the current global climate change and human activities such as irrigation, soil salinization and dust mobilization increased in many playas during the last years, strongly affecting their ecosystem services and regional dust storm intensity (Liu et al., 2015;Bowen and Johnson, 2017;Zucca et al., 2021).Hence, to understand and forecast such processes causing severe consequences for larger areas, soil salinization studies especially in playa regions are necessary.
The extensive Iranian drylands, partly intensively populated today (Guilmoto and Oliveau, 2018), are characterized by an ongoing aridification trend (Sayari et al., 2013;Emadodin et al., 2019).Globally, these drylands are between those being most affected by soil salinization (Hassani et al., 2020).The Bajestan Playa in NE Iran is the second largest playa in this region (Fig. 1a).The playa and its catchment host several protected areas with diverse desert vegetation and wildlife (Sokhanvar et al, 2013; Department of Environment of Iran: https://en.doe.ir) (Fig. 1b).During the dry summer season the playa is influenced by the strong northeasterly "Wind of 120 Days", being linked with dust mobilization.Hence, it forms an important aeolian dust source for inner and southern Iran (Middleton, 1986;Alizadeh-Choobari et al., 2014;Thiam et al, 2021), where the number of climatic refugees increased during the last years (Ebrahimzadeh and Esmaeil Negad, 2017;Mianabadi et al., 2022).However, the soil salinity dynamics of the Bajestan region, significantly influencing local and supra-regional ecosystems and societies, was not systematically studied so far.
We monitored the spatio-temporal soil salinity dynamics in the Bajestan region using a comprehensive approach.Our methodology included a continuous annual remote sensing-based salinity time series spanning from 1992 to 2021, on-site field and laboratory measurements and systematic correlations with regional and global climate parameters.This continuous multi-dimensional approach sets our work apart from previous studies on soil salinity, which predominantly used irregular and/or short time series and hardly linked the soil salinity dynamics with regional and global climate.The goals of this pioneer study were to understand (i) spatio-temporal soil salinity changes in the Bajestan region during the last three decades, and (ii) the impact of regional and global climate changes on these processes, and (iii) the future potential of this comprehensive way to analyse soil salinization with respect to climate change and global warming also for other drylands.

Study area
The Bajestan Playa in the drylands of NE Iran, extending between 58 • 33′ − 58 • 46′N and 34 • 58′ − 35 • 1′E, has a size of ~ 3.725 km 2 and a catchment of ~ 34.170 km 2 .The playa surface is very flat, showing an altitude of ~ 760 m a.s.l in the northern central part where a saline playa lake is developed (Fig. 1b).
The playa basin had formed in the northern part of the Lut Block in the eastern part of the Central Eastern Iranian Microplate (Ahmadirouhani et al., 2017).In synchronicity with Quaternary climate changes, lakes had regularly formed within the Bajestan basin and other large topographical basins of Iran, the last phase occurring during the deglaciation around 14 ka.Subsequently, these lakes gradually dried out due to increased Holocene temperatures, forming the current playas (Torshizian and Moosavi Harami, 1999;Ashoori et al., 2005).The southern catchment is dominated by Cretaceous limestones and Mesozoic to Cenozoic igneous rocks, and the northern and western parts by Paleozoic to Cenozoic sedimentary and metamorphic rocks such as slates and phyllites (National Iranian Oil Company, 1957;Ahmadirouhani et al., 2017).The playa basin hosts different kinds of Quaternary deposits, i.e. mostly alluvial fans, dunes and paleolake deposits.The latter, especially found in the central part, are formed by up to several hundred meters thick salt deposits (Behrouzi, 1987;Ashoori et al., 2005).Large parts of the playa surface show extensive salt flats with salt polygons of different sizes.The size of the saline playa lake strongly varies seasonally and inter-annually.Due to the asymmetric catchment area most intermittent streams enter the playa from the northeast and southwest.Only small alluvial fans cover the zone between the northern foot-slopes and the playa basin, whereas a mega-fan extends from the southern shore towards the playa lake (Krinsley, 1970;Ullmann et al., 2019).Next to surficial hard rocks, the region is covered by bare and partly saline desert soils (Dewan and Famouri, 1961).
According to the Köppen-Geiger climate classification the region shows a semiarid to arid steppe climate today (Raziei, 2021), characterized by dry and extremely hot summers and cool winters (Aalijahan and Khosravichenar, 2021).Mean annual precipitation is 77.3 mm, mean annual temperature is 23 • C, and mean annual evaporation is 290.7 mm (https://www.irimo.ir)(Fig. 1c).During summer, the playa is influenced by the northeastern "Wind of 120 Days" linked with regional dust storm activity (Alizadeh-Choobari et al., 2014).During winter westerly disturbances dominate, bringing most of the annual precipitation (Alijani and Harman, 1985).
The region belongs to the Irano-Turanian phytogeographical region (Khosravi et al., 2014;Noroozi et al., 2019).The patchy vegetation is characterized by halophytic species such as Salsola, Salicornicum and Atriplex.Parts of the region host the Eftekhari, Mozarafi, Helali and Doruneh protected areas (Fig. 1b; unpublished map Department of Environment of Iran).Small agricultural fields (pomegranate, saffron, pistachio) are concentrated in the surroundings of the regional settlements.

Analysis and calculation of salinity indices
To analyze regional soil salinity in the study area encompassing parts of the Bajestan region (17909.5 km 2 ), we used remote sensing (RS) data between 1992 and 2021 (in the following "RS-study period").These were taken from Landsat 5 Thematic Mapper (1992-2001), Landsat 7 Enhanced Thematic Mapper Plus (2002-2012) and Landsat 8 OLI sensors (2013-2021) (https://earthexplorer.usgs.gov;row 36, path 160).In accordance to our field campaign in 2021, we only used images from the first week of July.The dry summer season is ideal to detect soil salinity by RS, since low soil moisture causes surficial salt precipitation being easily detectable by well-developed topographical features, such as efflorescence, crusts, puffy structures and cracks (Metternicht and Zinck, 2003;Stavi et al, 2021).The data were pre-processed by radiometric, atmospheric and spectral normalization using ENVI 5.3 software.Soil salinity was calculated by applying ten different salinity indices using combinations of the blue, green, red and/or the near infrared RS-bands (subsequently called "EC-Pred" [= predicted]) (see Supplementary Table ST-1), applying the TerrSet geospatial monitoring and modeling software (https://clarklabs.org/terrset).

Calculation of wet-/water-covered and vegetated areas
We used widely applied spectral RS-indices and coding in Google Earth Engine to identify fluctuating wet/water-covered and vegetated areas that are known to influence soil salinization processes (Hou and Rusuli, 2022;Perri et al., 2018): Wet/water-covered surfaces during February and March were calculated using the modified Normalized Difference Water Index (mNDWI) (equation (1) (Feyisa et al., 2014;Szabó et al., 2016;Wang et al., 2020).The vegetation cover for July was calculated using the Normalized Difference Vegetation Index (NDVI) (equation 2), being closely related to green and healthy vegetation showing values > 0.2 (Tucker, 1979;Benedict and Jaelani, 2021;Cetin et al., 2022;Zeren Cetin et al., 2023).

Field work and laboratory analysis
We collected 70 surficial soil samples from 10 * 10 m plots during the first half of July 2021.To obtain representative samples not affected by random effects, five sub-samples were collected from the corners and the center of each plot and mixed into one sample.The geographic positions were determined from the plot center using a hand-held GPS.To representatively encompass the geomorphic-environmental diversity of the study area, the sampling locations were based on remote sensingderived soil salinity (see chapter 3.3) as well as on topographic and geomorphological features (Fig. 1b).
Electrical Conductivity (EC), an indicator of soil salinity, was measured in the laboratory of the Ferdowsi University of Mashhad (Iran) following Shaw (1994) (subsequently called "EC-Lab" [= measured in laboratory]).

Statistical soil salinity analysis
To construct a robust soil salinity detection model, the ten RSderived salinity indices (EC-Pred) (chapter 3.1.1)were calibrated with the field-based soil salinity (EC-Lab) (chapter 3.2) using a linear regression model.Root mean squares were used to indicate the inequality of variance and correlation of that model.Given that salinity index SI (Khan et al. 2005) yielded the best agreement between EC-Pred and EC-Lab (R 2 = 0.896; see supplementary Table ST-1), this index was used to categorize six different soil salinity classes according to Jenks Natural Break Algorithm that is based on natural groupings inherent in the data and minimizes the variance within and maximizes the deviations between the classes (Chen et al, 2013;Stavi et al., 2021, and citations therein).Based on these soil salinity classes, we used the Land Change Modeler for Ecological Sustainability (LCM) to detect spatiotemporal soil salinity changes within the study area over the RS-study period (https://clarklabs.org;Eastman and He, 2020).

Climatic analysis
To analyse possible climatic influences on the regional soil salinity dynamics, we used annual averages of daily resolved local temperature and precipitation data from the nearby Tabas synoptic station (Fig. 1a) between 1961 and 2021 (hereafter referred to as "climatic time scale"; see supplementary Table ST-2).The longer time scale of the climatic analyses (61 years) compared with the shorter RS-study period (30 years) allows to better evaluate the robustness of observed climatic trends.Furthermore, to understand the impact of global climatic factors on the climate of the Bajestan region, we compared 37 global climate indices with regional climatic parameters by applying Pearson correlation analyses (source: https://psl.noaa.gov/data/climateindices/list;for the MJO [Madden-Julian Oscillation] Index see: https://www.cpc.ncep.noaa.gov/products/precip/CWlink/MJO/mjo.shtml).The statistical methods are described below.

Regression analysis
Simple linear regression was applied to establish the relationship between two variables using a straight line being closest to the data, by determining slope and intercept and minimizing the regression error.
Furthermore, we applied sixth-order polynomial regression to analyze nonlinear relationships between an independent variable × and a dependent variable y in the form of E(y|x) (Montgomery et al., 2021).

Cross correlation
Cross correlation was used to examine possible relationships between regional climate parameters, the water-covered/wet area during February/March, and NDVI and soil salinity during July for the RS-study period 1992-2021.This allowed to determine the time lags of their best agreement for periods of up to 4 years (Menke and Menke, 2016;Begum et al., 2021).The basic equations are given below in equation ( 3) and ( 4) (Oktaviani and Setiawan, 2021): (mx = average rows of data on variable ×, my = average rows of data on variable y, r = cross-correlation value, d = time of highest correlation between the two variables, L = time lag).

Mann-Kendall test
Nonparametric statistics, being much less affected by outliers and other abnormalities (Lanzante, 1996), provides a measure of monotonic linear dependence (Davis & Sampson, 1986;Rossi et al., 1992).The most commonly used nonparametric test for trends in environmental variables is the Mann-Kendall test (MK) (Yadav et al., 2014;Pohlert, 2016;Aalijahan et al., 2023).This test does not require a specific data distribution, and can be used to analyze change trends, break points and abrupt changes within time series (Ye et al., 2013;Sung et al., 2017;Zhang et al., 2018).Hence, it was used for both the long climatic time scale 1961-2021 and the short RS-study period 1992-2021 to examine significant points of regional climate change.
For an assumed data series K (k 1 ,k 2 ,⋯k n ), n is the length of the data series.First, the cumulative statistic parameter S k , being the cumulative number of values at time i being greater than at time j, is calculated (equations ( 5) and ( 6): Subsequently, assuming random independence of the time series, the statistic parameter U'I is defined as follows (equation ( 7): U'I = 0, and E(S k ) and Var(S k ) represent expected value and variance of the cumulative value S k , respectively: U'I is a standard normal distribution, i.e. a statistical sequence computed from the sequence of time series X.At a significance level α, | U'I|>U α indicates an apparent trend change within the time series.
A. Khosravichenar et al.Using the inverse time series, we calculated U'I again applying the calculation procedure described above, where UI = -U'I and k = n,n − 1,⋯,1 (Liu et al., 2020).

Soil salinity dynamics
The calibration of RS-derived salinity index SI (EC-pred) with fieldderived soil salinity (EC-Lab) resulted in six different soil salinity classes between non-saline to extreme salinity (Table 1; exemplary photos with typical features of the six classes please see in supplementary Figure SF-1).During the studied month of July no open water areas existed in the study area, hence no separate class was designated for it.
Figs. 2 and 3 illustrate the spatio-temporal dynamics of the soil salinity classes.The very high and high classes showed continuously increasing trends throughout the RS-study period, whereas the low class continuously decreased.In contrast, the moderate class hardly showed any systematic change.Also the non-saline and extreme classes did not systematically change during most of the time, however, their areas strongly increased in 2021.Fig. 4 shows area changes between these classes during the three studied decades.In particular, outside the playa basin moderate saline areas went over to the high and very high classes during the first and third studied decades, whereas within the playa basin extreme saline areas passed over to the high and very high classes throughout the complete study period (Fig. 4a).The changes of different classes to the extreme class are shown in Fig. 4b.Obviously, this happens primarily for areas within the playa basin that were originally classified as the very high, high and moderate classes.The latter two contributed particularly during the third studied decade.The extent of the wet/ water-covered area within the playa basin in February/March showed erratic fluctuations during the RS-study period, and was mainly confined to the areas of extreme and very high (summer) soil salinity.
The NDVI was used to analyse vegetation changes during the RS-study period.Fig. 5a shows the vegetation distribution exemplarily for 1992.Most vegetation is found on areas of low soil salinity (34.2%), closely followed by the non-saline (33.6%), moderate (30.0%), high (7.7%),very high (0.6%) and extreme (0.0%) classes.For the entire RSstudy period fluctuations in vegetation density over several years are visible, and extraordinary peaks were reached in 2019 and 2020 (Fig. 5b).However, the average vegetation distribution remained the same over almost the entire study period.

Temporal evolution of regional annual precipitation and temperature
Mean annual precipitation did not show a clear trend for the long climatic time scale 1961-2020.However, linear and 6th order polynomial regression for the shorter RS-study period 1992 -2021 show a decrease especially from 2004 until 2016 (Fig. 6).At the end of the RSstudy period a short three-year wet period occurred 2018-2020, but in 2021 precipitation dropped again.
Mean, maximum and minimum annual temperature linearly increased since 1961 (Fig. 7).However, on the long climatic time scale the increase of minimum annual temperature, being especially strong during the 1970 s and 1980 s, was much higher (linear regression slope a = 0.094 • C) than for mean and maximum annual temperature (each a = 0.018 • C).In contrast, for the shorter RS-study period the slope of the increase of maximum annual temperature (a = 0.069 • C) was higher than for average (a = 0.060 • C) and especially minimum (a = 0.037 • C) annual temperature.Also 6th order polynomial regression, recognizing periodic trends of the temperature parameters, confirms the increasing trends for both time scales.

Mann Kendall trend analysis
The Mann-Kendall trend analysis of mean annual precipitation on the longer climatic time scale did not show a significant trend, since the significance level of ± 1.96% was hardly exceeded despite some intersections between the UI and U'I lines.However, on the shorter time scale of the RS-study period the 2009 intersection is significant, forming a sudden downward change point of precipitation (Fig. 8a).
Maximum and mean annual temperature indicate long-term downward trends during the first years of the longer climatic time scale (Fig. 8b, c), followed by strong upward trends lasting until today.However, the intersections of these temperature variables during their upward trends are not significant (Fig. 8b, c).In contrast, the upward trend of minimum annual temperature started much earlier, and the intersection in 1988 forms a sudden upward change point of this parameter (Fig. 8d).On the shorter time scale of the RS-study period all three temperature variables show upwards trends, and consecutive intersections occurred in the time series of all these variables.For mean and maximum annual temperature the intersections in 2007 and 2012 suggest significant trend changes (Fig. 8b, c), as well as the intersections of minimum annual temperature in 1996 and 2008 (Fig. 8d).
The regional temperature variables correspond well with the global temperature increase (Fig. 8e), what holds especially true for the shorter RS-study period.The global temperature increase intensified since ca.1992, possibly causing the strong upward trends of maximum and mean annual temperature since that time.Interestingly, minimum annual temperature highly corresponded to the global temperature rise already since 1975, i.e. much earlier than the other parameters.

Influence of regional climate changes and the wet/water-covered area in the playa basin on soil salinity, and the role of regional vegetation
To depict the influence of regional climate changes on soil salinity in the Bajestan region, we cross-correlated regional climate parameters with soil salinity (see supplementary Table ST-4; Fig. 9).Furthermore, to understand a possible influence of the wet/water-covered area in the playa basin during February/March on summer soil salinity, we crosscorrelated its extent with the extreme soil salinity class being dominant within the basin.These cross-correlations were calculated for time lags (delays) of up to 4 years.Negative lags indicate a lead of the climatic parameters and of the wet/water-covered area extent, whereas zero time lags indicate a simultaneous reaction of soil salinity to these factors.Furthermore, soil salinity was also correlated with the NDVI, and given the strong influence of regional climate on dryland vegetation density the NDVI was also correlated with mean annual precipitation and temperature.ST-3).
We observe a significant positive correlation (r = 0.604) between soil salinity and regional annual precipitation for a time lag of two years, and a significant negative correlation (r = -0.473)for a time lag of one year (Fig. 9a).Maximum and mean annual temperature showed the largest significant positive correlations with soil salinity for a time lag of zero, with r-values of 0.609 and of 0.482, respectively.In contrast, minimum annual temperature showed a significant negative correlation (r = -0.433)with a time lag of one year (Fig. 9b-d).The correlation between the wet/water-covered area in the playa basin during February/March and extreme soil salinity in July did not show any significant correlation (Fig. 9e).The correlation of the NDVI with mean annual precipitation and temperature showed a very slight significant positive correlation of r = 0.373 and a significant negative correlation of r = -0.493,respectively, both with time lags of zero.The positive correlation of NDVI with mean annual precipitation was possibly caused by the strong vegetation increase 2019 -2020 during an extraordinarily wet period.Finally, the correlation of NDVI with the low soil salinity class showed a significant positive correlation of r = 0.619 with a time lag of zero.No significant correlation was found between NDVI and the non-saline class (Fig. 9f-i).

Teleconnection analysis
The correlations between global climatic indices and climatic parameters in the study area are shown in Table 2.
Out of the 35 regarded teleconnection indices only some sub-indices of the Madden-Julian Oscillation (MJO) Index show significant negative correlations with mean annual precipitation in the Bajestan region, i.e. sub-indices 9, 10, 1 and 2 (the MJO indices at 20, 70, 80 and 100 • eastern longitude) (Table 2).Between those MJO -E70 Index 10 shows the highest negative correlation (r = -0.482,confidence level 0.99%).Also on a short multi-annual periodic time scale, expressed by 6th order polynomial regression, this negative correlation is well recognizable despite a missing long-term trend (Fig. 10a).
Nine out of the 35 global teleconnection indices, and the Global Land-Ocean Temperature Anomaly Index (GLOTI), show significant correlations with regional temperature (6 positive and 4 negative, respectively) (Table 2).Mean annual temperature shows the highest number of significant correlations, followed by maximum and minimum annual temperature.All regional temperature parameters, and especially mean and maximum annual temperature, show significant positive correlations with GLOTI that all show long-term upward trends.They also positively correlate on a multi-annual time scale, as is demonstrated by the 6th order polynomial regression of mean annual temperature and GLOTI (Fig. 10b).Hence, the increasing GLOTI Index obviously strongly controls the temperature increase in the Bajestan region.
Among the correlated teleconnection indices, POL and TSA show the most negative and positive effects on regional temperature.The annual temperature parameters and POL show negative correlations (Table 2), and according to 6th order polynomial regression these opposite trends also hold on a short multi-annual periodic time scale (see Fig. 10c for mean annual temperature).In contrast, regional temperature is positively correlated with TSA, and maximum annual temperature shows the highest correlation (Table 2).Similar with regional temperature TSA shows a long-term linear upward trend, but based on 6th order polynomial regression their correlation is not very clear on a multi-annual periodic time scale (see Fig. 10d for mean annual temperature).

Soil salinity dynamics in the Bajestan region and possible regional causes
Soil salinity in the Bajestan region strongly fluctuated between 1992 and 2021 (Figs. 2, 3), and areas with high and very high soil salinity generally increased whereas low saline areas decreased (Fig. 3).Outside the playa basin mainly moderate saline areas, not showing a total net change, were converted to the high and very high salinity classes (Fig. 4a).Hence, this change must have been paralleled by a simultaneous conversion of low saline areas to the moderate class.In contrast, increases of high and very high saline areas within the playa basin occurred mostly at the costs of the extreme class.Conversely, changes towards extreme soil salinity almost exclusively occurred within the playa basin mostly at the costs of very high and high saline areas (Fig. 4b).Hence, our data demonstrate different soil salinity dynamics within and outside the playa basin: Outside the basin, areas of high and very high soil salinity continuously increased over the RS-study period mainly at the costs of low saline areas.Within the basin, areas with high, very high and extreme soil salinity regularly changed into each other (Figs. 2, 4), but without a clear long-term trend (Fig. 3).
The Bajestan region generally became drier (Fig. 6) and warmer (Fig. 7) during the last three decades.However, whereas the drying trend is only recognizable since 1992 but not on the longer climatic time scale since 1961 (Fig. 6), all regional temperature parameters and especially minimum annual temperature distinctly increased since 1961 (Fig. 7).The close relationship between continuously increasing soil salinity, decreasing precipitation and increasing temperature since 1992 is also confirmed by the significant correlations between the merged highest soil salinity classes and regional climatic parameters (Fig. 9, supporting Table ST-4).Hence, similar with former studies (Bannari and Al-Ali, 2020;Rafik et al., 2022) increasing high and very high saline areas outside the playa basin (Fig. 3) are obviously linked with increasing aridity and temperature.Furthermore, there are different time lags between the considered parameters (Fig. 9): Mean annual precipitation positively correlates with a two-year time lag and negatively correlates with a one-year time lag with soil salinity, whereas mean and maximum annual temperature significantly positively correlate with zero time lags.Therefore, high soil salinity in the measurement year can be traced back to higher precipitation two years and low precipitation one year earlier.In contrast, high temperatures increase soil salinity during the same year.The significant positive correlation between the NDVI (vegetation density) and low soil salinity with a zero time lag means that during warmer years both vegetation density and low saline areas decrease, what is also mirrored by the significant negative correlation between mean annual temperature and NDVI (Fig. 9).Since most extreme saline areas (during summer) are found within the playa basin (Fig. 2), we cross-correlated this class with the extent of the wet/water-covered area during spring.However, the not significant correlation indicates only a minor effect of wet/watercovered areas on soil salinization in the playa basin.
Summing up, for the areas outside the playa basin showing the largest surfaces with high and very high soil salinity (Fig. 2) we suggest the following causal chain of soil salinization: (i) During wet years easily soluble salts are leached and mobilized from salt-bearing rocks and sediments into the Bajestan catchment (e.g. from saline lacustrine sediments originating from larger paleolake phases; Behrouzi, 1987;Ashoori et al., 2005;Etikala et al., 2021).While most of these salts are transported by ephemeral tributaries and saline groundwater along the hydraulic gradient into the playa basin, another part still remains in the catchment outside the basin.(ii) A subsequent dry year then promotes high soil salinity in the catchment, as indicated by the significant negative correlation between soil salinity and precipitation with a time lag of one year.Presumably, salt mobilization and transport still continue with the remaining water from the previous wet year even in this second dry year.However, the moisture is not sufficient to transport large salt amounts further into the playa basin.Therefore, most of the salts remain in the catchment area outside the basin.(iii) High temperatures during the year of soil salinity measurement lead to increasing evaporation, causing strong surficial salt accumulation due to capillary rise from ground or soil water (Forkutsa et al., 2009;Stavi et al., 2021).In this context, the lower correlation of minimum annual temperature with soil salinity compared with mean and maximum temperature can be explained by the exponential relationship between temperature and evaporation (Gabler et al., 2008).This demonstrates a comparatively lower importance of minimum temperaturemostly measured during night-time -for evaporation.In parallel, vegetation density decreases during such warmer years.Whether this decrease is only caused by increasing evaporation and related lower available soil moisture, or also by increasing soil salinity, cannot be resolved by our data.
The extent of extreme saline areas within the playa basin remained largely continuous since 1992 (Fig. 3), and regular increases and decreases of this class were at the costs of very high and high saline areas and vice versa (Fig. 4b).However, there was an unprecedented strong increase of the extreme saline areas in 2021 (Fig. 3), following an unprecedented three-years wet period 2018-2020 (Fig. 6).Hence, the strong increase of extreme saline areas in 2021 can possibly be explained by an unprecedented intensive transport of easily soluble salts from the catchment into the playa basin during the wet period.This, however, obviously represents an extreme event rather than a trend.The long-term stability of the non-saline areas since 1992 opposes the general trend of increasing regional soil salinity outside the playa basin, indicating that these areas are little affected by the processes described above.Their unprecedented strong increase in 2021 (Fig. 3) especially in the SE part of our study area (Fig. 2) can thus not be explained so far.However, we hypothesize that strong leaching and transport of easily soluble salts into the playa basin during the preceding 2018 -2020 wet period could have happened in particular from this SE part of our study area, causing the simultaneous increase of extreme saline areas within the basin.Another explanation could be an enlargement of agricultural activities in the non-saline areas after the wet period, as agricultural use is generally tied to these areas (see Table 1).In addition, a combination of both processes is also conceivable.
Unlike former soil salinity studies that generally used large and irregular time steps or only short sample periods (e.g.Wu et al., 2014;Scudiero et al., 2015;Dasgupta et al., 2015;Bannari and Al-Ali, 2020;Hassani et al., 2020;Rafik et al., 2022;Ma et al., 2023), we analysed soil salinity in the Bajestan region annually over a continuous period of 30 years.Therefore, with this respect we could not only identify long-term soil salinity changes, but also the causal process relationships in detail.Next to the long-term fostering effect of higher temperatures and lower precipitation for increasing soil salinity (Bannari and Al-Ali, 2020;Rafik et al., 2022), we also identified short-term time lags linked with these main controlling factors.Furthermore, unlike many former study areas that were located near coastlines (Dasgupta et al., 2015;Bannari and Al-Ali, 2020) or in intensively irrigated agricultural areas (Wu et al., 2014;Scudiero et al., 2015;Rafik et al., 2022;Ma et al., 2023), our study area is located far away from any coastline and only shows minor human influence.Hence, our detailed results should be representative for the soil salinity dynamics in the vast inland drylands of SW Asia without strong human influence, helping to understand the current and future soil salinity dynamics of these areas.Nevertheless, further factors probably also influenced soil salinity in the Bajestan region or have at least biased our results, such as aeolian activity leading to aeolian removal of salt crusts (Wang et al., 2012) or the regional groundwater dynamics (Wang et al., 2008).Therefore, to quantify these factors further investigations are necessary.

Global climatic controls on the soil salinity dynamics in the Bajestan region
The soil salinity dynamics in the Bajestan region between 1992 and 2021 was mainly controlled by regional precipitation and temperature.Correlation of the latter with different global climate indices showed that regional annual precipitation best correlated with different subindices of the global MJO Teleconnection Index (see Table 2).MJO represents a dominant intra-seasonal tropical oscillation, characterized by an eastward zonal shift of the whole pattern of 12,000-20,000 km along the equator within 30-60 days (Madden andJulian, 1971, 1972;Zhang, 2005).It is one of the main factors controlling precipitation in tropical and subtropical regions, widely impacting wind and precipitation patterns in SW-Asia (Iran, Afghanistan, Pakistan) (Barlow et al., 2005;Madhura et al., 2015;Ahmadi et al., 2019).Accordingly, we see a negative relationship between several MJO sub-indices and mean annual precipitation in the Bajestan region, what is also displayed by a completely inverse behavior of mean annual precipitation and the MJO -E70 Index 10 on a multi-annual time scale, i.e. when applying 6th order polynomial regression (Fig. 10a).All regional temperature parameters showed the highest correlations with the global warming GLOTI Index.GLOTI displays global anomalies of Land Surface Air Temperature (LSAT) over land and sea ice and of Sea Surface Temperatures (SSTs) over ice-free water (NASA GISS, 2022).Furthermore, among the considered teleconnection indices POL and TSA show the highest correlations with regional temperature.The negatively correlated POL index, appearing during the cold season, consists of one main polar anomaly center and separate opposite centers over Europe and NE-China.Thus, it reflects major changes in the strength of the circumpolar circulation by revealing the accompanying systematic midlatitude circulation changes over large parts of Europe and Asia (Gao et al., 2019).The positively correlated TSA index indicates the SSTs in the Gulf of Guinea between 30 • W − 10 • E and 20 • S − 0 • (Enfield et al., 1999).
Summing up, we could identify clear influences of global climatic phenomena on the climate of the Bajestan region.Compared with former studies (Pankova and Konyushkova, 2013;Bannari and Al-Ali, 2020;Okur and Örcen, 2020;Khadimov et al., 2022), this allowed a much deeper look into the linkages between global climate change and regional soil salinization.Whereas mean annual precipitation in the Bajestan region did not show significant systematic changes on the longer time scale since 1961, regional annual temperature was clearly impacted by the global climate change since that time.Due to the location of the Bajestan region in a climatic transitional region, both high and mid-latitude synoptic systems (TSA, POL) play a role in this context.Similar to other regions (Dai, 2011;Ghavidel Rahimi et al., 2017;Sarkar and Maity 2021), we observe the first effects of global warming on the regional climate since the 1970 s.This is particularly evident in the strongly increasing annual minimum temperatures.Mean and maximum annual temperature, on the other hand, were influenced by global warming only during the RS-study period since 1992.This indicates an intensified regional impact of global warming during the recent years (Figs. 6,7,8).Given some influence also of minimum annual temperature on regional soil salinity (Fig. 9), despite missing RS data before 1992 we suspect an influence of global climate change on soil salinity already since the 1970 s.Due to ongoing global warming (Pörtner et al., 2022) probably causing further increasing regional annual temperatures, we also expect a further increase of soil salinization outside the Bajestan playa basin.This will lead to increasing salt dust emissions and increasingly threaten the diverse desert vegetation and wildlife in the regional protected areas.Furthermore, assuming repeated extreme soil salinity conditions within the playa basin such as in 2021, we also expect increasing salt dust emissions from the basin over the coming decades.

Conclusions
During this study, we applied a multi-disciplinary research approach encompassing remote sensing, field and laboratory analyses as well as statistical analyses of regional and global climatic parameters to study soil salinity changes and their relation with global climate change in the Bajestan drylands of NE Iran 1992-2021.Thanks to our so far not applied annual time resolution over three decades, we gained much deeper insights into the causal chain of regional soil salinization than during previous studies.Our study shows that both regional annual precipitation and temperature changes control regional soil salinity, but that soil salinity responds to precipitation changes with time lags of up to two years.This probably reflects the time that is needed for the transport of leached soluble salts from their source areas following humid years.In contrast, the reaction of soil salinity to higher temperatures due to increasing capillary rise does not show any time lag.When closer looking into the study area, no systematic changes of extreme soil salinity within the playa basin but just an extreme salinization event in 2021 were observed.In contrast, the areas outside the basin showed a systematic long-term increase of high and very high soil salinity, whereas low saline areas decreased.This was possibly caused by conversions of the moderate to high and very high saline areas, and a parallel increase of the moderate at the costs of low saline areas.Whereas regional annual precipitation did not systematically change during the last six decades, minimum annual temperature increased since the 1970s and mean and maximum annual temperature during the last years.Accordingly, when correlating these regional climate parameters with global climatic indices, the temperature parameters showed the highest correlations with the global warming GLOTI index as well as with the POL and TSA indices being linked with both high-and midlatitude climatic phenomena, all showing increasing trends during the last years.Given ongoing global climate change, a further increase of  a) -(g), as soil salinity does neither influence regional climate nor the extension of the wet/water-covered area, and patchy dryland vegetation does not significantly influence regional climate.soil salinity outside the playa basin should hence be expected.This will have serious consequences for saline dust emissions and the diverse desert vegetation and wildlife in the regional protected areas, and saline dust emissions could furthermore be reinforced by extreme salinity events within the playa depression as that in 2021.This will possibly also increase the number of climatic refugees from the affected regions during the coming decades.Altogether, our multi-disciplinary approach shows a high potential for future applications also in other salinityaffected drylands, making it a pioneer study for research on soil salinization in the context of global warming as a base for the planning of effective strategies in combatting dryland degradation and to forecast the aeolian dust dynamics.

Fig. 2 .
Fig. 2. Maps of the annual distribution of the six RS-derived soil salinity classes and the wet/water-covered area in the study area between 1992 and 2021.The corresponding values of each class and year can be found in Supplementary TableST-3).

Fig. 3 .
Fig. 3. Total area changes of the six soil salinity classes during the RS-study period 1992-2021.

Fig. 4 .
Fig. 4. Maps showing changes between the soil salinity classes for the three studied decades by comparing 1992 with 2001, 2002 with 2011, and 2011 with 2021 (calculated by Land Change Modeler for Ecological Sustainability).(a) Changes from different classes to the merged high and very classes, (b) Changes from different classes to the extreme class.

Fig. 5 .
Fig. 5. NDVI-based vegetation density in the study area (vegetated areas show NDVI-values > 0.2) (see Supplementary Table ST-3 for the NDVI values).(a) Map exemplarily showing the extension of the vegetated areas and the six soil salinity classes for May 1992 (this earlier month was chosen for this figure due to a better visibility of the vegetation coverage compared with the NDVI-study period in July).(b) Temporal evolution of the vegetated areas during July.

Fig. 6 .
Fig. 6.Mean annual precipitation 1961-2021.Linear and 6th order polynomial regression are separately shown for the long climatic time scale 1961-2021 and the RS-study period 1992-2021.

Fig. 7 .
Fig. 7. Temperature changes between 1961 and 2021.Linear and 6th order polynomial regression are separately shown for the long climatic time scale 1961-2021 and the RS-study period 1992-2021.(a) Maximum, (b) average, (c) minimum annual temperature.

Fig. 9 .
Fig. 9. Graphs of the different cross correlations taking into account possible time lags between soil salinity (merged high, very high and extreme soil salinity classes) and (a) mean annual precipitation, (b) mean annual temperature, (c) maximum annual temperature and (d) minimum annual temperature.Cross correlation between the wet/water-covered area within the playa basin during February/March and extreme soil salinity during July (e).Cross correlation between the NDVI and (f) mean annual precipitation, (g) mean annual temperature, and (h) the non-saline as well as (i) the low soil salinity class.Please note, positive time lags are not possible in reality for the cross-correlations shown in figures (a) -(g), as soil salinity does neither influence regional climate nor the extension of the wet/water-covered area, and patchy dryland vegetation does not significantly influence regional climate.

Fig. 10 .
Fig. 10.6th order polynomial relationships between annual climate parameters of the Bajestan region and global climatic indices showing the highest correlations for the RS-study period 1992 -2021: (a) MJO -E70 Index 10 and mean annual precipitation, (b) Global Land-Ocean Temperature Anomaly Index (GLOTI Index) and mean annual temperature, c) POL Index and mean annual temperature, d) TSA Index and mean annual temperature.

Table 1
Characteristics and surface area of the six soil salinity classes.