Quantitative Soil Wind Erosion Potential Mapping for Central Asia Using the Google Earth Engine Platform

: A lack of long-term soil wind erosion data impedes sustainable land management in developing regions, especially in Central Asia (CA). Compared with large-scale ﬁeld measurements, wind erosion modeling based on geospatial data is an e ﬃ cient and e ﬀ ective method for quantitative soil wind erosion mapping. However, conventional local-based wind erosion modeling is time-consuming and labor-intensive, especially when processing large amounts of geospatial data. To address this issue, we developed a Google Earth Engine-based Revised Wind Erosion Equation (RWEQ) model, named GEE-RWEQ, to delineate the Soil Wind Erosion Potential (SWEP). Based on the GEE-RWEQ model, terabytes of Remote Sensing (RS) data, climate assimilation data, and some other geospatial data were applied to produce monthly SWEP with a high spatial resolution (500 m) across CA between 2000 and 2019. The results show that the mean SWEP is in good agreement with the ground observation-based dust storm index (DSI), satellite-based Aerosol Optical Depth (AOD), and Absorbing Aerosol Index (AAI), conﬁrming that GEE-RWEQ is a robust wind erosion prediction model. Wind speed factors primarily determined the wind erosion in CA ( r = 0.7, p < 0.001), and the SWEP has signiﬁcantly increased since 2011 because of the reversal of global terrestrial stilling in recent years. The Aral Sea Dry Lakebed (ASDLB), formed by shrinkage of the Aral Sea, is the most severe wind erosion area in CA (47.29 kg / m 2 / y). Temporally, the wind erosion dominated by wind speed has the largest spatial extent of wind erosion in Spring (MAM). Meanwhile, a ﬀ ected by the spatial di ﬀ erence of the snowmelt period in CA, the wind erosion hazard center moved from the southwest (Karakum Desert) to the middle of CA (Kyzylkum Desert and Muyunkum Desert) during spring. According to the impacts of land cover change on the spatial dynamic of wind erosion, the SWEP of bareland was the highest, while that of forestland was the lowest.


Introduction
During the past few decades, global climate change and human disturbance have meant that land degradation has become one of the most serious environmental problems of the 21st century [1]. Despite the lack of strong political will, the land degradation problem has attracted much attention

Study Area
The most common definition of CA is the official one of the Soviet Union, which includes the five former Soviet republics of Kazakhstan (KZ), Uzbekistan (UZ), Turkmenistan (TK), Kyrgyzstan (KG), and Tajikistan (TJ). The total area of CA is nearly 4 × 10 5 km 2 , which is mainly covered with bareland and sparse vegetation. The landform types of CA are mainly plains and hills. Additionally, the mountains (Tianshan Mountain, Pamir Mountain, and Altai Mountain), which are known as the "Water Tower of Central Asia", are mainly distributed in the southeast [57]. As the Tianshan Mountain and Pamir Mountain block rain clouds that should enter CA from the east and south, CA is one of the largest land-locked arid regions in the world [58]. Most of CA lies in an arid climatic zone, which has low annual precipitation (less than 300 mm), a high air temperature, and strong evaporation. Five large temperate deserts (Karakum Desert, Kyzykum Desert, Muyunkum Desert, Sarresi-Atyray Desert, and Aralkum Desert) are distributed from the southwest to middle east ( Figure 1). Additionally, desertification caused by large-scale agriculture practices has been an issue since 1960 and enhanced climate change presents many economic, social, and environmental problems in CA [15,24,59,60]. The most notorious example is the Aral Sea Crisis, which has been considered to be one of the planet's worst environmental disasters of the 21st century [59]. The large-scale construction of irrigation canals has reduced runoff from Syr Darya river and Amu Darya river into the Aral Sea, which in turn reduced the Aral Sea surface area from 68,000 square kilometers in 1960 to less than 7000 square kilometers in 2016 [60]. Meanwhile, a new anthropogenic desert known as Aralkum Desert in the eastern dry basin appeared in 1960. Salt and dust storms, which are caused by wind erosion occurring in Aralkum Desert, represent one of the most serious problems for human health and agricultural activities in CA [16].
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 27 of empirical models based on batch geospatial data and high-performance computing. The main conclusions could be beneficial for desertification control and land resource management in CA.

Study Area
The most common definition of CA is the official one of the Soviet Union, which includes the five former Soviet republics of Kazakhstan (KZ), Uzbekistan (UZ), Turkmenistan (TK), Kyrgyzstan (KG), and Tajikistan (TJ). The total area of CA is nearly 4 × 10 5 km 2 , which is mainly covered with bareland and sparse vegetation. The landform types of CA are mainly plains and hills. Additionally, the mountains (Tianshan Mountain, Pamir Mountain, and Altai Mountain), which are known as the "Water Tower of Central Asia", are mainly distributed in the southeast [57]. As the Tianshan Mountain and Pamir Mountain block rain clouds that should enter CA from the east and south, CA is one of the largest land-locked arid regions in the world [58]. Most of CA lies in an arid climatic zone, which has low annual precipitation (less than 300 mm), a high air temperature, and strong evaporation. Five large temperate deserts (Karakum Desert, Kyzykum Desert, Muyunkum Desert, Sarresi-Atyray Desert, and Aralkum Desert) are distributed from the southwest to middle east ( Figure 1). Additionally, desertification caused by large-scale agriculture practices has been an issue since 1960 and enhanced climate change presents many economic, social, and environmental problems in CA [15,24,59,60]. The most notorious example is the Aral Sea Crisis, which has been considered to be one of the planet's worst environmental disasters of the 21st century [59]. The largescale construction of irrigation canals has reduced runoff from Syr Darya river and Amu Darya river into the Aral Sea, which in turn reduced the Aral Sea surface area from 68,000 square kilometers in 1960 to less than 7000 square kilometers in 2016 [60]. Meanwhile, a new anthropogenic desert known as Aralkum Desert in the eastern dry basin appeared in 1960. Salt and dust storms, which are caused by wind erosion occurring in Aralkum Desert, represent one of the most serious problems for human health and agricultural activities in CA [16].

Data Collection and Source
The meteorological data included the wind speed, soil moisture, and snow depth, which were derived from Global Land Data Assimilation System 2.1 (GLDAS2.1) integrating satellite and ground-based observational products [61]. Three other sets of climate assimilation data, including The Fifth Generation ECMWF Atmospheric Reanalysis Data (ERA5), NCEP Climate Forecast System Reanalysis (CFSR), and the Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System (FLDAS), were used to investigate the wind speed variability across CA ( Table 1). The soil mechanical composition, soil organic matter, and several other soil properties were obtained from the Harmonized World Soil Database (HWSD) and OpenLandMap (OLM), which are based on machine learning predictions from a global compilation of soil profiles and samples. A total of six standard depths (0, 10, 30, 60, 100, and 200 cm) were divided in the OLM dataset (Table 1), due to the lack of soil calcium carbonate content data in GEE datasets and based on the finding that a nonlinear positive correlation exists between the soil pH and soil calcium carbonate [62,63]. Huang, et al. [62] found that the relationship between the calcium carbonate content and pH value of surface soil in East Central Asia has the highest R 2 when they simulated the factors with an exponential equation. Liu, et al. [63] found that the soil pH and CaCO 3 content have a non-linear positive correlation during a study conducted in China. Based on more than 15,000 different soil mapping units, we proposed an exponential equation (Equation (1)) to quantify the relationship between the soil pH and soil calcium carbonate (CaCO 3 ).
where pH is the soil pH of different soil types in HWSD and CaCO 3 is the soil calcium carbonate content in HWSD (%). NDVI was derived from the NASA Terra Moderate-Resolution Imaging Spectroradiometer (MODIS) Vegetation Indices (MOD13Q1). The National Aeronautics and Space Administration Shuttle Radar Topographic Mission (NASA-SRTM) provided Digital Elevation Model (DEM) data, which were used to calculate the slope data. The land cover data were provided by the European Space Agency (ESA)-based Climate Change Initiative (CCI) 300 m global land cover data products developed using the GlobCover unsupervised classification chain and by merging multiple available Earth observation products.
Ground observation wind speed data (2 m height) were derived from NOAA Global Surface Summary of Day (GSOD), which includes global data obtained from the United States Air Force (USAF) Climatology Center. We collected the daily ground measurement wind speed data from more than 400 weather stations in CA. Due to political or other reasons, weather stations in the former Soviet Union were abandoned and several new weather stations were established between 1990 and 2010 [60]. Most weather stations currently working were established in the 1960s. Therefore, we integrated ground observations of the wind speed based on 204 weather stations from 2000 to 2019 (Figure 1). The average wind speed data of all-weather stations were used to study the temporal variation characteristics of the wind speed in CA. The visibility, data which could be employed to calculate the Dust Storm Index (DSI), were derived from GSOD. Additionally, the RS technique has provided a new perspective on the validation of soil wind erosion. Aerosol data included the Aerosol Optical Depth (AOD) and Absorbing Aerosol Index (AAI), which were used to compare and validate the wind erosion in this research. AOD was derived from MODIS based on the Multi-angle Implementation of Atmospheric Correction (MAIAC). Moreover, AAI was derived from the Sentinel-5 Precursor, which launched on 13 October 2017. Details on the research data are listed in Table 1.

GEE-RWEQ
The comprehensive assessment of wind erosion in a large-scale region like CA is complex and challenging. Dozens of parameters are employed to calculate the soil wind erosion modulus by using field-scale models, such as WEPS. GEE is a cloud computing platform specially designed to process raster data, including satellite images, climate assimilation grids, and other geospatial data. The advantage of the GEE platform lies in the instant access, manipulation, visualization, and real-time analysis of large amounts of geospatial data [52]. Therefore, the advent of GEE made it possible to launch global-scale environmental mapping and monitoring programs [53]. This is of great potential for integrating an environmental model on the GEE platform to build a GEE-based production framework. Furthermore, most developing countries where resources are limited have suffered from various environmental problems, including droughts, flooding, deforestation, soil degradation, and dust storms caused by wind erosion [16,53,64,65]. These countries often lack monitoring sites and networks for environmental problems, making these problems more serious [64][65][66][67]. In this study, by using multisource geospatial data, we present a fully automated algorithm for mapping NRT monthly wind erosion dynamics at a global scale using the GEE platform.
Although the computational efficiency should not be a concern in the GEE platform, the limited ground observation data present a big challenge for simulating soil wind erosion. Therefore, it is necessary to build a simplified and more practical model that can estimate the SWEP at a large scale on the GEE platform. Although it has a lower accuracy than mechanistic wind erosion models, this relatively simple model is not limited by the input data, location, and scale of the study area. RWEQ has been extensively tested, and good agreements between model results and field measurements were found in previous studies. In this study, a GEE cloud computing-based RWEQ model was developed to conduct quantitative mapping of soil wind erosion in a ground measurement limited area. As mentioned above, based on the progress of Earth observation and numerical modeling, several parameters that used to be filed, measured, or calculated can be easily acquired on the GEE platform. Although most of the input parameters are retained in GEE-RWEQ, the soil roughness factor (K') is difficult to estimate during farming production on a regional scale. Ouyang, et al. [68] replaced the soil ridge roughness with the roughness caused by topography, and it was calculated by the Smith-Carson equation. Because this equation has been widely used in many regions [68][69][70][71], it can applied when the study area is scaled up from a field to a region. Due to the limitations of wind erosion estimation based on RS on regional scales, the combined crop factor (C) was simplified based on previous findings [48][49][50]68,72].
The GEE-RWEQ involved basic equations, as follows [37]: where SWEP is the amount of soil wind erosion potential per unit area (kg/m 2 ); Q max is the maximum transport capacity (kg/m); x is the distance from the upwind edge of the field (m), set to 55 m for the study area; s is the critical field length (m); WF is the weather factor (kg/m); EF is the erodible fraction (dimensionless); SCF is the soil crust factor (dimensionless); K is the soil roughness factor (dimensionless); and C represents combined crop factors (dimensionless). The weather factor can be calculated as where SW is the soil wetness (dimensionless), SD is the snow cover factor (dimensionless), u 2 is the wind speed at 2 m (m/s), u t is the threshold wind speed at 2 m (assumed 5 m/s), N is the number of wind speed observations (u 2 > u t ) in the period, N d is the number of days in the period, ρ is the air density (kg/m 3 ), and g is the acceleration due to gravity (m/s 2 ). The erodible fraction (EF) and soil crust factor (SCF) can be calculated as where Sa is the sand content (%), Si is the silt content (%), Sa/Cl is the sand to clay ratio (%), OM is the organic matter (%), CaCO 3 is the calcium carbonate (%), and Cl is the clay content (%). The soil roughness factor (K ) can be calculated as [70] where α is the slope gradient (degree), which can be calculated by the Digital Elevation Model (DEM). The combined crop factor (C) can be calculated as [70,73] where SC is the vegetation coverage (%), NDVI soil is the NDVI value of a bare soil pixel, and NDVI max is the maximum NDVI value of the study area.

Model Performance Evaluation
Because of the diverse land cover types and large area, it is extremely difficult to measure wind erosion for a whole region. Additionally, in the past two decades, almost no research has conducted field measurements on wind erosion in CA. Therefore, validation methods that can evaluate the reliability of wind erosion model results need to be proposed. Considering that the ground observed dust storm can indicate the frequency and intensity of wind erosion events, DSI was used to validate the spatial variation of SWEP. DSI was calculated based on the meteorological record-visibility, which can represent the frequency and intensity of wind erosion events. DSI was first proposed by McTainsh [74] in the National Collaborative Project on Indicators for Sustainable Agriculture (NCPISA). Based on the relationship between meteorological records and DSI, O'Loingsigh, et al. [75] used daily visibility data acquired from 180 long-term meteorological stations to investigate a long-term national wind erosion record  in Australia. DSI is a methodology employed for monitoring wind erosion based on long-term daily meteorological observations. At present, it is generally accepted as an indicator of broad-scale wind erosion rates in Australia, Iran, and Northeast Asia [75,76]. Based on weather codes relating to wind erosion or visibility, wind erosion events were divided into three categories: (a) Severe Dust Storms (SDS); (b) Moderate Dust Storms (MDS); (c) Local Dust Events (LDE). The DSI was calculated using the following equation [75]: where i is ith value of n stations for i = 1 to n, SDS is a severe dust storm (visibility < 200 m), MDS is a moderate dust storm (200 m < visibility < 1000 m), and LDE is a local dust event (1000 m < visibility < 20,000 m). Due to the lower population and urban density, soil mineral particles produced by wind erosion are the main source of atmospheric aerosols in CA [24]. Therefore, the satellite-derived AOD data were used to evaluate the reliability of SWEP simulated by the RWEQ model. There are several satellite-based aerosol products, which have different spatial and temporal resolutions, such as CALIPSO Lidar Tropospheric Aerosol Profiles All sky data, VIIRS/SNPP Deep Blue L3 daily aerosol data, OMI/Aqua Multi-wavelength AOD Daily data, MODIS MO(Y)D08_M3 Terra (Aqua) Atmosphere Monthly data, MODIS MCD19A2 Terra & Aqua MAIAC Land Aerosol Optical Depth Daily data, and Sentinel-5P NRTI AER AI. However, most satellite-based aerosol products have a low spatial resolution and cannot meet the requirements of quantitative spatial comparisons [77][78][79]. In this study, we used the MODIS MCD19A2 dataset at the 0.47 µm blue band, along with the parameter Optical_Depth_047, which has a spatial resolution of 1 km [79]. Another aerosol dataset named the Absorbing Aerosol Index (AAI), with a 0.01-degree spatial resolution, was extracted from the Sentinel-5P NRTI AER AI product. Because the Sentinel-5P was launched on 13 October 2017, the aerosol dataset was released in 10 July 2018 [80]. Therefore, we used the 2019 annual average AAI to compare with the 2019 SWEP in this study.

Technical Flowchart of this Study
The Land Cover Change (LCC), which is influenced by both climate change and human activity, usually affects wind erosion on surface roughness and soil physical and chemical characteristics. Therefore, we studied the SWEP of different land cover types and SWEP changes caused by the conversion of different land cover types. In this study, we chose ESA-CCI 300 m global land use land cover data products developed using the GlobCover unsupervised classification chain and by merging multiple available Earth observation products. Based on the United Nation Land Cover Classification System's (LCCS) plant functional types (PET), the CCI-LC map is classified into 22 land types. According to the study area land characteristics, the land cover types are reclassified into nine categories based on the look-up table-conversion of CCI-LC classes to PET in the product user guide [81].
Based on the objective of this study, this manuscript is organized as presented in the technical flowchart ( Figure 2). The research consists of four main steps: First, based on a time-series decomposition model, the wind speed variability of ground measurement data and reanalysis data was explored; second, by using multi-source geospatial data, the monthly SWEP across CA was generated based on GEE-RWEQ, and we explored the spatiotemporal variation of SWEP between 2000 and 2019; third, based on DSI and satellite-based AOD, validation was conducted to test the reliability of annual SWEP; and finally, we investigated the responses of wind erosion to ground measurement wind speed change and land cover change.
was explored; second, by using multi-source geospatial data, the monthly SWEP across CA was generated based on GEE-RWEQ, and we explored the spatiotemporal variation of SWEP between 2000 and 2019; third, based on DSI and satellite-based AOD, validation was conducted to test the reliability of annual SWEP; and finally, we investigated the responses of wind erosion to ground measurement wind speed change and land cover change.
The ArcGIS10.6 software was implemented for land cover type reclassification in this research. The Pearson correlation coefficient (r) was calculated in R 3.6.3. The time-series decomposition model based on the "tseries" package was also run in R 3.6.3. The exponential fitting of the soil pH and soil calcium carbonate was performed in MATLAB 2018a.

Variability of the Daily Average Wind Speed across CA
As the key factor of wind erosion, wind speed variability plays a vital role in wind erosion dynamics. A host of studies have reported that there was a declining trend in the global near-surface wind speed from 1970 to 2010 [14,[17][18][19]. However, a recent study described an increase in the global wind speed during a particular year. Zeng, et al. [22] found that, after several decades of global terrestrial stilling, the wind speed has increased rapidly across the globe since 2010. Although Zeng, et al. [22] have investigated the global temporal variation of the wind speed, further studies are required because of the sensitivity of CA to global climate change [60].
To better understand the temporal variations of wind speed, it is possible to decompose wind speed time series data into sub-components by a time-series decomposition model. In this study, we The ArcGIS10.6 software was implemented for land cover type reclassification in this research. The Pearson correlation coefficient (r) was calculated in R 3.6.3. The time-series decomposition model based on the "tseries" package was also run in R 3.6.3. The exponential fitting of the soil pH and soil calcium carbonate was performed in MATLAB 2018a.

Variability of the Daily Average Wind Speed across CA
As the key factor of wind erosion, wind speed variability plays a vital role in wind erosion dynamics. A host of studies have reported that there was a declining trend in the global near-surface wind speed from 1970 to 2010 [14,[17][18][19]. However, a recent study described an increase in the global wind speed during a particular year. Zeng, et al. [22] found that, after several decades of global terrestrial stilling, the wind speed has increased rapidly across the globe since 2010. Although Zeng, et al. [22] have investigated the global temporal variation of the wind speed, further studies are required because of the sensitivity of CA to global climate change [60].
To better understand the temporal variations of wind speed, it is possible to decompose wind speed time series data into sub-components by a time-series decomposition model. In this study, we used a multiplicative decomposition model, which is more effective when a seasonal value changes over time [82]. The calculation of this model included the following three steps [83]. The trend component was first determined and removed from time series by using the moving averages method. Secondly, the seasonal component was calculated and centered by averaging all periods for each time unit. In this study, a time-series decomposition multiplicative model was applied to Ground Measurement Wind Speed (GMWS) data. Figure  used a multiplicative decomposition model, which is more effective when a seasonal value changes over time [82]. The calculation of this model included the following three steps [83]. The trend component was first determined and removed from time series by using the moving averages method. Secondly, the seasonal component was calculated and centered by averaging all periods for each time unit. In this study, a time-series decomposition multiplicative model was applied to Ground Measurement Wind Speed (GMWS) data. Figure   In order to ensure the consistency of the reanalysis data input by the model and the actual observation data in CA, the trends of four reanalysis data (GLDAS2.1, ERA5, CFSR, and FLDAS) were introduced to conduct comparisons with GMWS. The comparison result showed that the daily average derived from GLDAS2.1 has the highest correlation with the ground measurement wind speed (Supplementary Materials Figure S1). Furthermore, in order to better compare the relationship between the wind speed in trend and seasonal components, we decomposed the GLDAS2.   Figure 5 shows the distribution of the annual mean SWEP in CA over the most recent 20 years (2000-2019). SWEP exhibits significant spatial variation, which has a range of 0-256 kg/m 2 . The SWEP in the southwest CA is higher than that in the southeast and north CA, where the vegetation coverage  Figure 5 shows the distribution of the annual mean SWEP in CA over the most recent 20 years (2000-2019). SWEP exhibits significant spatial variation, which has a range of 0-256 kg/m 2 . The SWEP in the southwest CA is higher than that in the southeast and north CA, where the vegetation coverage and soil moisture are higher. During the past 20 years, the Aral Sea dry lake bed (ASDLB), which is one of the most active dust sources, was the most severe wind erosion area in CA (47.29 kg/m 2 /y), followed by Kyzylkum Desert (10.64 kg/m 2 /y), Karakum Desert (10.58 kg/m 2 /y), and Muyunkum Desert (6.81 kg/m 2 /y). Due to the dramatic shrinkage of Aral Sea from the second half of the 20th century, the ASDLB, also known as the Aral Sea Desert, was covered with the original salts and chemicals of the water [24]. The toxic sediments of the Aral Sea were blown away by strong winds and formed white sandstorms. These toxic particles from the dry Aral Sea lake bed had been found in Japan, Norway, Greenland, and even in the South Pole [84]. From 2000 to 2019, the annual mean value of SWEP was 3.45 kg/m 2 , with the lowest value occurring in 2010 and the highest value occurring in 2015. As mentioned above, the wind speed exhibited a significant increase around 2011, and the area and intensity of wind erosion have also increased significantly since 2011. However, the SWEP gradually decreased by −6.85 kg/m 2 /y, over ASDLB, from 2011 ( Figure 5). This reversal may have been caused by the recovery of the water body of the Aral Sea since 2010 [85]. This means that more dry lake beds are covered by water bodies and less bareland suffers from wind erosion. The Amu Darya River (ADR) and Syr Darya River (SDR), From 2000 to 2019, the annual mean value of SWEP was 3.45 kg/m 2 , with the lowest value occurring in 2010 and the highest value occurring in 2015. As mentioned above, the wind speed exhibited a significant increase around 2011, and the area and intensity of wind erosion have also increased significantly since 2011. However, the SWEP gradually decreased by −6.85 kg/m 2 /y, over ASDLB, from 2011 ( Figure 5). This reversal may have been caused by the recovery of the water body of the Aral Sea since 2010 [85]. This means that more dry lake beds are covered by water bodies and less bareland suffers from wind erosion. The Amu Darya River (ADR) and Syr Darya River (SDR), which are the two largest rivers across CA, are the principal water suppliers of the Aral Sea. Since the 1960s, many irrigation canals have been constructed in the middle and lower reaches of ADR and SDR [60]. Therefore, according to ESA land cover data, more than 43% of the irrigation cropland of CA survived on these canals, especially in Amu Darya Delta (ADD) and Syr Darya Delta (SDD). The main land cover type of these two deltas is irrigated cropland, which accounts for more than 20% of the total irrigated cropland in CA. Because of the high vegetation coverage, according to Figure 5, the ADD and SDD regions have a lower SWEP than the surrounded area. Similarly, on the edge of the Kyzylkum Desert, oasis agriculture that relies on irrigation also greatly reduces SWEP.

The Spatiotemporal Variation of Wind Erosion across CA
The seasonal variation characteristics of SWEP in CA are shown in Figure 6. Although it is usually consistent with the spatial distribution of the annual mean SWEP, the spatial pattern of SWEP varies in different seasons. As we can see in Figure 6, more land suffered from severe wind erosion in CA during the spring than in other seasons (spring: 1.48 kg/m 2 > summer: 0.70 kg/m 2 > winter: 0.69 kg/m 2 > autumn: 0.59 kg/m 2 ). Additionally, obvious wind erosion exists in several famous desert regions, such as the Karakum Desert, Kyzylkum Desert, and Aralkum Desert during spring. Due to the impact of snow cover, less wind erosion exists in north CA during the winter. However, the wind speed of the Aral region in winter markedly exceeds that in other seasons, especially in December. Therefore, the ASDLB region has higher SWEP in the winter. Figure 6 shows that the north Kazakhstan region, Kyrgyz, Tajikistan, has been slightly affected by wind erosion. Additionally, in summer, SWEP shows a significant high value in the middle reaches of Amu Darya ( Figure 6). As shown in Figure 7, the SWEP displays significant monthly temporal variability, especially in ASDLB. March is the most severe month of wind erosion in the ASDLB. Alternatively, due to strong winds and dry surface soil in December, January, and February, sandstorms frequently occurred in ASDLB. We can confirm this from the true color image of MODIS (Supplementary Figure S2). Due to the difference in the solar radiation energy received by different latitudes, the melting time of snow cover varies in different regions. Figure 6 shows that the center of wind erosion moves from southwest to northeast during the spring (March, April, and May). These factors can explain why significant wind erosion in other famous deserts occurs in different months. The most severe month of wind erosion in Karakum is April, and that in Kyzylkum Desert and Muyunkum Desert is May. However, in general, Central Asia has suffered from the most severe wind erosion in April and the most widespread wind erosion in May.
varies in different regions. Figure 6 shows that the center of wind erosion moves from southwest to northeast during the spring (March, April, and May). These factors can explain why significant wind erosion in other famous deserts occurs in different months. The most severe month of wind erosion in Karakum is April, and that in Kyzylkum Desert and Muyunkum Desert is May. However, in general, Central Asia has suffered from the most severe wind erosion in April and the most widespread wind erosion in May.

Impacts of Ground Measurement Wind Speed Changes on the SEWP
Based on the research results mentioned above, the wind speed was found to be the dominant factor of wind erosion. Therefore, in this study, we analyzed the influence of wind speed as a factor of climate change on wind erosion. Based on the ground measurement wind speed data, we investigated the influence of wind speed changes on wind erosion. Figure 8 shows the strong and

Impacts of Ground Measurement Wind Speed Changes on the SEWP
Based on the research results mentioned above, the wind speed was found to be the dominant factor of wind erosion. Therefore, in this study, we analyzed the influence of wind speed as a factor of climate change on wind erosion. Based on the ground measurement wind speed data, we investigated the influence of wind speed changes on wind erosion. Figure 8 shows the strong and significant positive correlation (r = 0.700, p < 0.001) between the average GMWS and average ground SWEP. According to the trend line of these two sets of time-series data, the turning point roughly occurred between 2010 and 2011. The average GMWS showed a slowly decreasing trend (−0.07 ms −1 decade −1 ) before the turning point in 2011, while it displayed a significant increasing trend during the period of 2011-2019 (+0.6 ms −1 decade −1 , p < 0.001). The recent increasing rate is almost tenfold the decreasing rate in the first decade. The average ground SWEP also showed a slightly decreasing trend during the period of 2000-2010 (−0.027 kgm −2 decade −1 ), while it displayed a significant increasing trend (+0.37 kgm −2 decade −1 , p < 0.001). Shao, et al. [86] found that the global monthly mean dust concentration decreased from 2000 to 2012. This shows that the end of the quiet period of dust activities in Central Asia or globally marks the beginning of an active period of dust activities. Based on current research, it seems reasonable to relate the dust trend to the climate trend, especially the reversal in global terrestrial stilling [22]. Additionally, the highest monthly average SWEP, which appeared in May 2014, was more than 2.4 kg/m 2 . The seven months with the strongest wind erosion (SWEP > 1.5 kg/m 2 ) were March (2), April (2), May (2), and December (1). There were nine months (May: 3, March: 3, April: 2, Decembeer: 1) with an average wind speed greater than 3.5 m/s. Moreover, the highest monthly average wind speed, with a value of 3.79 m/s, appeared in Decembr 2015. Overall, a similar tendency of wind speed and wind erosion was observed for the past two decades. Both show a distinguishable declining trend, and then a sudden remarkable increase, and before slowly declining or finally stabilizing. During the study period, the wind speed (+0.38 ms −1 decade −1 , p < 0.001) and SWEP (+0.34 kgm −2 decade −1 , p < 0.001) increased very quickly from 2000 to 2019, indicating a more serious soil degradation and air pollution problem in CA.

Divergence of SWEP from Different Land Cover Types
According to the land characteristics of the study area, the ESA CCI land cover types can be reclassified into nine categories (cropland irrigated, cropland rain-fed, forestland, shrubland, Additionally, the highest monthly average SWEP, which appeared in May 2014, was more than 2.4 kg/m 2 . The seven months with the strongest wind erosion (SWEP > 1.5 kg/m 2 ) were March (2), April (2), May (2), and December (1). There were nine months (May: 3, March: 3, April: 2, Decembeer: 1) with an average wind speed greater than 3.5 m/s. Moreover, the highest monthly average wind speed, with a value of 3.79 m/s, appeared in Decembr 2015. Overall, a similar tendency of wind speed and wind erosion was observed for the past two decades. Both show a distinguishable declining trend, and then a sudden remarkable increase, and before slowly declining or finally stabilizing. During the study period, the wind speed (+0.38 ms −1 decade −1 , p < 0.001) and SWEP (+0.34 kgm −2 decade −1 , p < 0.001) increased very quickly from 2000 to 2019, indicating a more serious soil degradation and air pollution problem in CA.

Divergence of SWEP from Different Land Cover Types
According to the land characteristics of the study area, the ESA CCI land cover types can be reclassified into nine categories (cropland irrigated, cropland rain-fed, forestland, shrubland, grassland, sparse vegetation land, bareland, urbanland, and waterbody). The last subplot of Figure 9 shows the areas of different land cover types across CA in 2018. Figure 9 shows that the monthly average SWEP and its change rate of different land cover types were substantially different. The monthly average SWEP of bareland was more than 0.836 kg/m 2 , followed by shrubland (0.572 kg/m 2 ), sparse vegetation land (0.203 kg/m 2 ), grassland (0.073 kg/m 2 ), irrigated cropland (0.043 kg/m 2 ), rainfed cropland (0.033 kg/m 2 ), and forestland (0.017 kg/m 2 ). This relationship is basically consistent with previous research conducted in CA and surrounding regions [20,87]. We also compared the soil wind erosion modulus with respect to regions with similar conditions or the use of different methodologies (Table 2). Li, et al. [20] assessed the soil wind erosion modulus variation in CA (including Xinjiang, China) between 1986 and 2005. Zhang, et al. [87] investigated the RWEQ-based soil wind erosion, which was validated by 137 Cs in Inner Mongolia (IM), during the time period of 1990-2015. Compared to other arid or semiarid regions, CA has relatively higher rates of soil wind erosion, which may be the result of the widespread distribution of deserts and wind speed increase in the past decade. Grassland has a relatively lower soil wind erosion rate, because it is mainly distributed in northern CA, which has a more humid climate and less erodible underlying surface condition [12]. Although the time periods and dataset sources were different, from the perspective of the wind erosion diversity of land cover, the wind erosion result of our research is reliable.   There was a significant increase in SWEP in CA in the past two decades. However, the change rate of wind erosion varied among regions with various types of land cover. Specifically, the change rate of wind erosion was the highest for bareland (0.0978 kg/m 2 /y), followed by shrubland (0.0730 kg/m 2 /y), sparse vegetation land (0.0206 kg/m 2 /y), grassland (0.0062 kg/m 2 /y), irrigated cropland (0.0045 kg/m 2 /y), rainfed cropland (0.0038 kg/m 2 /y), and forestland (0.0016 kg/m 2 /y). Combined with the areas of different land covers, we could calculate the total amount of soil wind loss for different land covers in the period of 2000-2019. More than 2.8255 × 10 11 t soil was eroded by the wind across CA during the past few decades. The soil wind erosion of bareland (1.838 × 10 11 t) contributed more than 65% soil loss by the wind in CA, followed by shrubland (0.3907 × 10 11 t), sparse vegetation land (0.3380 × 10 11 t), grassland (0.1812 × 10 11 t), rainfed cropland (0.0517 × 10 11 t), irrigated cropland (0.0223 × 10 11 t), and forestland (0.0028 × 10 11 t). According to the continuous changes in the land cover area in CA during the period of 2000-2018, the cropland, forestland, urbanland, and shrubland showed an increased trend, while bareland and grassland showed a decreased trend (Supplementary Table S1). More than 2.67 × 10 km 2 land had undergone land cover change, including from bareland to grassland, sparse vegetation land to grassland, grassland to cropland (rain-fed), sparse vegetation land to cropland (rain-fed), and waterbody to bareland. In order to remove wind speed variability effects on SWEP, we calculated the annual SWEP of 2018 based on the wind speed data of 2000 to compare it with the annual SWEP of 2000. We found that the LCCs with the strongest inhibitions of wind erosion activity were bareland into shrubland (−0.782 kg/m 2 ), sparse vegetation land into forestland (−0.106 kg/m 2 ), and grassland into forestland (−0.073 kg/m 2 ). In comparison, the conversions of land cover which accelerated wind erosion the most were waterbody into bareland (+3.784 kg/m 2 ), sparse vegetation land into bareland (+1.124 kg/m 2 ), and grassland into bareland (+0.490 kg/m 2 ).

Validation of the GEE-RWEQ Model
Due to lack of long time series and wide range of ground-measured wind erosion data in CA, validation of the GEE-RWEQ model is challenging. Furthermore, due to the large area and complex terrain conditions, almost no previous research has conducted wind erosion field measurements in CA. Considering that most of the local dust storms are caused by surface wind erosion in CA [24], the dust storm index (DSI) can be used as a proxy to evaluate the wind erosion model performance [75]. Therefore, we used the DSI based on weather station visibility records to evaluate the reliability of the SWEP spatial distribution. Figure 10 displays the annual mean DSI across CA from 2000 to 2019. This map was interpolated by the annual average DSI of more than 200 weather stations based on Natural Neighbor Interpolation method. From Figure 10, we can see that the southwest desert region of CA has a high DSI, which means very frequent dust storms. However, there are some high values in the southeastern parts of CA, where there is a lower wind erosion risk. According to the research of Liu, et al. [89], affected by the strong southwest winds, the dust particles were transported from the western desert to eastern mountains and valley. Moreover, dust episodes were observed in these regions. Additionally, the southeastern parts of CA are the most densely populated areas in CA. Most of the weather stations across CA are located around densely populated cities. The anthropogenic pollutants will also be recognized as dusty weather due to reduced visibility. The spatial distribution of DSI is generally consistent with the spatial pattern of the annual SWEP. Although visibility records based on weather stations are a valuable and useful data resource for wind erosion monitoring, several limitations still exist. The low spatial density of weather stations is a challenge for conducting highly accurate wind erosion mapping, especially in the southeastern part of CA. Therefore, we need to obtain higher spatial resolution and more continuous SWEP verification data. It should be pointed out that satellite-based atmospheric aerosols, which refer to solid and liquid particles suspended in the atmosphere, have strong spatial correlations with wind erosion.
Although visibility records based on weather stations are a valuable and useful data resource for wind erosion monitoring, several limitations still exist. The low spatial density of weather stations is a challenge for conducting highly accurate wind erosion mapping, especially in the southeastern part of CA. Therefore, we need to obtain higher spatial resolution and more continuous SWEP verification data. It should be pointed out that satellite-based atmospheric aerosols, which refer to solid and liquid particles suspended in the atmosphere, have strong spatial correlations with wind erosion.  Figure 11 shows the spatial pattern of the annual AOD (a) and average AAI in 2019 (b), as well as the comparisons with SWEP in the Aral Sea region (ASR). From Figure 11, we can see that the Aral Sea region and its southwest surrounding area has the highest value in CA. This is because the dust of the Aral Sea is transported to the southwest under the action of the dominant wind-northeast wind [90]. Therefore, we chose ASDLB as a research hotspot area to compare with SWEP. Linear relationships between SWEP and aerosol parameters (AOD and AAI) were found, as shown in Figure 11. The results show that they had moderate positive linear relationships, with r values of 0.5623 (p < 0.001) and 0.5660 (p < 0.001), respectively. Ultimately, although it is difficult to verify the SWEP value, the comparison results obtained from the perspective of spatial and temporal distribution patterns showed that the RWEQ-based SWEP data in CA were reliable.  Figure 11 shows the spatial pattern of the annual AOD (a) and average AAI in 2019 (b), as well as the comparisons with SWEP in the Aral Sea region (ASR). From Figure 11, we can see that the Aral Sea region and its southwest surrounding area has the highest value in CA. This is because the dust of the Aral Sea is transported to the southwest under the action of the dominant wind-northeast wind [90]. Therefore, we chose ASDLB as a research hotspot area to compare with SWEP. Linear relationships between SWEP and aerosol parameters (AOD and AAI) were found, as shown in Figure  11. The results show that they had moderate positive linear relationships, with r values of 0.5623 (p < 0.001) and 0.5660 (p < 0.001), respectively. Ultimately, although it is difficult to verify the SWEP value, the comparison results obtained from the perspective of spatial and temporal distribution patterns showed that the RWEQ-based SWEP data in CA were reliable.

Discussion
In this study, our results show that there are significant spatial and temporal differences in the wind erosion in CA. Controlled by the latitude zonality and vertical zonality, higher SWEP values are primarily distributed in the southwestern part of CA, which has low vegetation coverage and

Discussion
In this study, our results show that there are significant spatial and temporal differences in the wind erosion in CA. Controlled by the latitude zonality and vertical zonality, higher SWEP values are primarily distributed in the southwestern part of CA, which has low vegetation coverage and more fragile surface soil [91], while lower SWEP is mainly concentrated in the northern part of CA. Based on the land cover map of 2018, more than 89% of rain-fed cropland and 78% of forestland are distributed in these regions. Furthermore, due to latitude zonality and vertical zonality, the precipitation in these areas is higher than in other places of CA, and the temperatures in these areas are lower. Therefore, the higher soil moisture caused by lower evapotranspiration will reduce the wind erosion to a certain degree. Affected by the restoration of the Aral Sea in recent years, the vegetation coverage and other underlying surface factors of ASDLB are getting better [85]. Although the SWEP showed a decreasing trend in ASDLB (−6.85 kg/m 2 /y), this region was still the most severe wind erosion area in CA. In addition, our results show that the SWEP has a clear seasonal and monthly variation. The land threatened by wind erosion has the largest range in spring, especially in May. Due to the difference in solar radiation heat at different latitudes and altitudes, the snowmelt period varies in different regions of CA [92]. The higher soil moisture caused by snowmelt and the snow cover both affected the movement of the wind erosion center across CA during spring [37]. Therefore, considering the major deserts of CA are located in different latitudes and altitudes, the severe wind erosion regions in CA will also migrate over time. The severe wind erosion regions move from the southwest CA (Karakum Desert) to the middle CA (Kyzylkum Desert and Muyunkum Desert) during the spring (MAM). However, due to the special meteorological conditions in certain areas of CA, such as the middle reaches of the ADR, severe wind erosion occurs in summer, but not spring, respectively [24]. We calculated the monthly average wind speed of four weather stations in this region. A comparison was made for the monthly average wind speed of all the weather stations in CA and the four weather stations in this region (Supplementary Figure S3). The most extreme value of the surface wind speed in CA appeared in March, while the most extreme value of the wind speed in the middle reaches of the ADR appeared in July. This result demonstrates that wind speed plays the key role in the spatial distribution of SWEP.
As the most dominant factors of wind erosion, wind speed variability, such as wind stilling or wind stilling reversal, account for the majority of the spatial-temporal variation of wind erosion in CA. Although global terrestrial stilling has been confirmed by many pieces of research, most studies have only looked at global or regional wind speed changes from the 1980s to 2010, and few have involved recent (after 2010) wind speed changes [14,[17][18][19]. According to several climate assimilation datasets (GLDAS, ERA5, CFSR, and FLDAS) and a ground measurement dataset (GSOD), we found a turning point of wind speed stilling during the period of 2009-2012 in CA. This finding is supported by other studies that have reported that global terrestrial stilling has rebounded over the past few decades and has increased rapidly since 2010 [21][22][23]. Our research proves that the increase rate of the average wind speed in CA (0.6 m s −1 decade −1 ) is higher than the increase rate of the average global wind speed (0.24 m s −1 decade −1 ) over the same period, which means that stronger wind erosion occurred in CA. Indeed, the result shows a strong and significant positive correlation (r = 0.7) between the average GMWS and average SWEP (p < 0.001). A number of studies have demonstrated that CA is more sensitive to climate change compared to the global average [12,20,57,60,91,93]. While it is widely acknowledged that the global wind speed rebound is beneficial to the wind power industry for the near future [22], this study suggests that more severe wind erosion activity happened in CA.
According to the significant differences in natural conditions, such as the air temperature and precipitation, and the disturbance of human activities such as irrigation, the land cover in CA exhibited strong spatial differentiation. We found that SWEP differs greatly in different land cover types. This result is roughly consistent with a previous study on wind erosion in CA and surrounding regions [20,87]. Most of the shrubland in CA is made up of deserts and xeric shrublands in which Haloxylon ammodendron, Calligonum aphyllum, and Ephedra lomatolepis, as well as grasses such as Agropyron fragile, grow [94][95][96]. These kinds of vegetation have good windproof and sand fixing functions in CA. In addition, LCCs are closely related to wind erosion activity and affect each other. From 2000 to 2018, more than 2.6 × 10 5 km 2 of land has changed land cover types (Supplementary  Table S1). We compared the SWEP for the land with the LCC from 2000 to 2018, in which the effect of the wind speed variability was removed. The conversion of bareland to shrubland helped reduce wind erosion by −0.782 kg/m 2 /y. Furthermore, due to the shrinkage of the Aral Sea, the conversion of waterbody to bareland increased wind erosion by +3.784 kg/m 2 /y. The restoration of the Aral Sea has not only superficially reduced the possibility of wind erosion in ASDLB, but also increased the increased vegetation coverage of ASDLB caused by the higher groundwater level, making the long-distance transport of dust difficult. Although the wind speed has shown an increasing trend in the past 10 years, the wind erosion risk in the Aral Sea area is gradually decreasing due to the continuous recovery of the Aral Sea area. In the past 30 years, a large number of engineering projects have attempted to improve the Aral Sea environment, directly or indirectly. Although Aral sea restoration is the most effective way to restrain wind erosion in ASDLB, the complex political relations among countries in the Aral Sea basin make cross-border water management difficult [84]. Besides, in other parts of CA, more effective measures should be taken for wind erosion artificial control, for example, increasing the grassland area in regions with a suitable temperature and precipitation, developing cropland by using limited water resources, and planting cold-and drought-resistant shrub vegetation in bareland.
In this study, the RWEQ model was adopted as a wind erosion model in the GEE cloud computing platform. Compared with the local computing platform, GEE can process large amounts of geospatial data in a short time, which means that its processing power is completely unconstrained by time and space. Therefore, we do not need to spend a lot of time on downloading, preprocessing, and model running of a large amount geo-spatial data, which can greatly shorten the time required for long-term wind erosion mapping. As mentioned above, GEE-RWEQ provides the possibility of wind erosion monitoring in developing regions lacking on-site monitoring data. Meanwhile, the GEE platform makes it easier for researchers to publish their results for decision makers, and even the public [55]. Therefore, our research has a broader application value for decision makers than previous studies on wind erosion. In the future, we will interactively develop Earth Engine App to explore our result, which can then be used by experts and non-experts alike.
The RWEQ is a process-based, field-scale, empirical model that can quantitatively estimate wind erosion. However, the RWEQ was initially developed for the middle western area of the United States [37]. Therefore, it still presents some limitations in other regions [37,50,72,97]. Although the most important input parameters were retained in this study, the dataset required by some parameters was unavailable on GEE platforms. Therefore, several factors were simplified to simulate the global-scale wind erosion more effectively and more accurately. The soil moisture data were used to simulate how the surface wetness influences the wind speed required to erode the soil. Additionally, the cosine of the slope gradient which was calculated by DEM represented the soil roughness factor. On the other hand, we only used wind erosion-related data such as visibility data and remote sensing data for the verification of wind erosion in CA, but this still has uncertainties on a global scale. Central Asia, which is one of the most severe wind erosion regions, is restricted in terms of wind erosion modeling studies due to the lack of wind erosion measurement data for this region. Therefore, more ground soil loss measurement data on a global scale should be acquired to conduct more verification studies.

Conclusions
In this study, we developed a fully automated algorithm for quantitatively mapping wind erosion based on the Google Earth Engine, processed terabytes of geo-spatial data, and retrieved spatial and temporal patterns of monthly SWEP in CA, over 20 years (2000 to 2019). Several conclusions were reached in our study, as follows: (1) With respect to the conventional methods, GEE-RWEQ does not require any ground measurement data, which need lots of manpower and resources, especially in developing countries or sparsely populated regions. However, based on the Cloud computing platform, GEE-RWEQ uses climate assimilation data, soil property data, vegetation data, terrain data, and other underlying data to automatically generate high spatial resolution NRT soil wind erosion potential products. After verification using ground observation-based DSI and satellite-based AOD, the results still reach an acceptable accuracy and can be used for quantitative wind erosion mapping. This methodology provides new ideas for the construction and use of empirical models based on batch geospatial data and high-performance computing; (2) According to the comparison of GMMS and SWEP, the wind speed is the main driving factor of wind erosion (r = 0.7, p < 0.001). Affected by the wind speed variability, the SWEP decreased first and increased remarkably during 2011. From the perspective of the temporal and spatial distribution, due to the sparse vegetation distribution and special meteorological conditions, the deserts in southwestern Central Asia are most affected by wind erosion, especially in ASDLB (47.29 kg/m 2 /y). The severe wind erosion period of CA occurred in spring (MAM), especially in May. We also found that the SWEP distribution has obvious latitude zonality due to the distribution of snow cover and the start time of snow melt, and the wind erosion hot spot in spring moves from the southwest to central area across CA; (3) Land cover change has strong effects on the soil wind erosion in CA, with the most obvious being the conversion of bareland into the water body in ASDLB. Affected by the restoration of the Aral Sea, the SWEP in this area has shown a downward trend (−6.85 kg/m 2 /y) since 2011. Additionally, the conversion of bareland to shrubland helped reduce wind erosion by −0.782 kg/m 2 /y. According to the SWEP variation based on LCC, more effective measures should be taken to maintain wind erosion artificial control, for example, restoring the Aral Sea water area to prevent more bareland from being exposed to wind erosion, increasing the grassland area in regions with a suitable temperature and precipitation, developing cropland by using limited water resources, and planting cold-and drought-resistant shrub vegetation in bareland.