Drying and Wetting Trends and Vegetation Covariations in the Drylands of China

: The semi-arid and arid drylands of China, which are located in the inland region of Eurasia, have experienced rapid climate change. Some regions in particular, have shown upward trends in the observational records of precipitation. However, there is more to drying and wetting than just changes in precipitation which still have large uncertainties. Coherent results, however, can be obtained, at the regional scale, with the use of multiple indices as shown in the recent literature. We divided the drylands of China into three sub-regions, i.e., a semi-arid (SA), an eastern-arid (EA) and a western-arid (WA) region. Precipitation from the China Meteorological Administration (CMA) and Climatic Research Unit (CRU), statistical and physical drought indices, including the Standardized Precipitation Evapotranspiration Index (SPEI), the Palmer Drought Severity Index (PDSI), self-calibrating PDSI (sc_PDSI), Root zone soil moisture (Root_sm) and Surface soil moisture (Surf_sm) from Global Land Evaporation Amsterdam Model (GLEAM), and Normalized Di ﬀ erence Vegetation Index (NDVI) were used to identify temporal and spatial patterns in drying and wetting. Data were selected from 1982–2012, in line with the availability of the remotely sensed vegetation data. Results show that the drylands of China exhibits a pattern of wetting in the west and drying in the east. The semi-arid region in the east is becoming drier and the drought area is increasing, with the values of CMA_P, CRU_P, PDSI, sc_PDSI, SPEI-01,SPEI-06, SPEI-12, Root_sm, Surf_sm at − 1.064 mm yr − 1 , − 0.834 mm yr − 1 , − 0.050 yr − 1 ( p < 0.1), − 0.174 yr − 1 ( p < 0.1), − 0.014 yr − 1 , − 0.06, − 0.021 ( p < 0.1), − 0.257 × 10 − 3 m 3 m − 3 yr − 1 , − 0.024 × 10 − 3 m 3 m − 3 yr − 1 , respectively. The arid region generally exhibits a wetting trend, while the area in drought declines only in the western arid region, but not in the eastern arid part. In the semi-arid region, growing season (May to September) NDVI is signiﬁcantly correlated ( p < 0.1) with eight out of nine indicators. We show in this study that the semi-arid region needs more study to protect the vegetation ecosystem and the water resources.


Introduction
Drying and wetting patterns due to climate change have an important impact on social development and ecosystems. Therefore, they have always been one of the key issues of climate change. Dozens of drought indices, besides trends in precipitation, have been developed and extensively used in

Land Cover Data
The land use data was obtained from China's 1:100,000-scale land use status remote sensing monitoring database, including six land cover types and 25 sub-categories. We selected the year 2000 as a reference. The dataset is based on Landsat TM and ETM images, using human-computer interaction interpretation to interpret images that are not covered by Landsat TM images or that cover poor image quality, supplemented with Ambient 1 satellite CCD multispectral data. Unified quality

Land Cover Data
The land use data was obtained from China's 1:100,000-scale land use status remote sensing monitoring database, including six land cover types and 25 sub-categories. We selected the year 2000 as a reference. The dataset is based on Landsat TM and ETM images, using human-computer interaction interpretation to interpret images that are not covered by Landsat TM images or that cover poor image quality, supplemented with Ambient 1 satellite CCD multispectral data. Unified quality control and integration checks were performed on each database to ensure high quality and consistent interpretation [23][24][25]. Figure 1b and Table 1 show the land cover types. Grasslands and unused land are the main covers in the arid areas, while in the semi-arid areas the land was mainly covered by grassland and farmland. This study area is shaped with basins, mountains and plains, etc.

Meteorological Data
Precipitation and temperature from the China Meteorological Administration (CMA) and Climatic Research Unit (CRU) are selected. The CMA-precipitation is generated based on more than 2000 observatories in China [26]. Comparatively, the CRU TS3.22 dataset is more widely used [27].

Drought Index
Drought indices are widely used in drying and wetting studies. In this study, besides precipitation, three drought indices were used: SPEI, PDSI, sc_PDSI. SPEI is a statistical drought index which is based on a similar algorithm to that of the original SPI. Therefore, SPEI has a multi-scale characteristic. Unlike SPI, temperature effect is considered in SPEI by calculating the potential evapotranspiration [28]. The first calculation of the times series of the differences between precipitation (P) and potential evapotranspiration (PET) is given as: where k is the time scale, D k j is the P − PET difference of month j at the time scale k. After that, a three-parameter log-logistic distribution was chosen to model the D series. Finally, SPEI was generated by standardization [16]. In this study, the SPEI v2.3 data (available from http://sac.csic.es/ spei/database.html) was used. Similar to previous studies [29,30], the 1, 6, 12 months scale SPEI were used to analyze the drying and wetting trends in this study.
PDSI is a physical drought index calculated based on a two-layer soil water balance [2]. PDSI combines effects of the water deficit and impacts of the duration to the severity of drought. Firstly, the difference (d) between climatically appropriate for existing conditions (CAFEC) precipitation (p ) and the actual precipitation (p) is calculated. Then, the moisture anomaly index (Z index) is generated by applying a climatic characteristic coefficient K on (d). Finally, the PDSI value is caculated based on Z index and duration coefficient. The calculation formula is as follows: Water 2020, 12, 933 5 of 17 PET is the mean potential evapotranspiration, R is the mean recharge moisture, RO is the mean runoff, P is the mean precipitation, and L is the mean moisture loss. D is the mean absolute value of d. k is a climatic characteristic coefficient. Z is the moisture anomaly index which is generated by applying k on d. The 0.897 and 1/3 are constant coefficients which reflect the influence of previous moisture and month Z index. Meanwhile, the sc_PDSI changes the climatic characteristic coefficient k and the constant coefficient 0.897 and 1/3 in the PDSI based on local historical climatic data for each station or grid. Therefore, sc_PDSI presents better spatial compatibility than PDSI [3,31].
Previous studies have shown that the PET calculation scheme may have a substantial impact on the estimation of PDSI trends. The Penman-Monteith method is generally more accurate than the Thornthwaite method in PET [6,32,33]. The full Penman-Monteith equation can be expressed as follows: where λ is the latent heat of vaporization, ET is the evapotranspiration flux described as depth per unit time, ∆ is the slope of the saturation vapor pressure minus temperature, R n is the net radiation flux density at the surface, G is the soil flux, ρ a is the air density, c p is the specific heat of moist air under constant pressure, e s is the saturation vapor pressure at air temperature, e a is the actual vapor pressure of the air, r a is the aerodynamic resistance to turbulent heat and/or vapor transfer from the surface to some Z height above the surface, γ is the psychrometric constant, r s is the canopy conductance that describes the water vapor flowing from inside the leaf, vegetation canopy or soil to outside the surface. Therefore, this paper adopted PDSI and its deformation sc_PDSI with the Penman-Monteith equation for PET estimations [7]. The PDSI and sc_PDSI data used in this study is produced by Sheffield [34].

The GLEAM Soil Moisture
The GLEAM (http://www.gleam.eu/) model is used to calculate the evaporative flux over land with maximize the use of satellite-derived observations [35]. The root-zone soil moisture (Root_sm) and surface soil moisture (Surf_sm) can also be obtained from GLEAM. Because the use of constant parameters in GLEAM, the bare soil, short vegetation and vegetation with a tall canopy was chosen. The global model contains four modules including a separate module of the water budget (rain and snow) over the root-zone. The water balance of three land surface types was calculated separately and each type has a different number of layers. The bare soil has one layer (0-0.05 m) and short vegetation has the second layer (0.05-1.00 m). For the vegetation with a tall canopy, two layers (0.05-1.00 m and 1.00-2.50 m) were considered. The soil moisture content (ω) of each layer (l) on a given day (i) is calculated as follows: where F (l−1) and F (l) is the downward flux from the above layer, E l denotes the loss of soil water caused by evaporation, ∆z (l) is the thickness of the layer. F (l) is estimated by formula as follows: where ω f c is the field capacity. Simultaneously, microwave remote sensing datasets of surface soil moisture were used to constrain the uncertainty of the modelled results. Additionally, the Kalman filter assimilation approach was used to correct running water balance estimates at the daily time step.

Normalized Difference Vegetation Index (NDVI)
NDVI data is widely used to describe vegetation conditions; for example, the inter annual changes and spatial patterns [36,37]. In this study, the NDVI data was obtained from the third-generation data of the Global Inventory Monitoring and Modeling System (GIMMS). This NDVI data is aggregated by different AVHRR sensors and considers various detrimental effects such as calibration losses, orbital drift, and volcanic eruptions. The latest version of the GIMMS NDVI dataset covers the time range from July 1981 to December 2015 [38,39]. The equation used is: where NIR and VIS indicate spectral reflectance at the near-infrared and visible (red) spectral band range.

Statistical Methods
The least squares linear trend analysis was chosen to calculate the linear trends of precipitation and the drought indices. Simultaneously, the significance of these linear trends was examined by Mann-Kendall trend analysis [40,41]. The Mann-Kendall trend analysis is a nonparametric statistical method. The advantage of this method is that it does not require the data to have a normal distribution. Therefore, this method is widely applied to evaluate changes of trends of meteorological and hydrological time series. In addition, the Pearson correlation coefficients between NDVI and precipitation drought indices were calculated.

Spatial Patterns of Regional Climate and Vegetation
The annual average patterns of precipitation, temperature and vegetation in the study area from 1982 to 2012 are shown in Figure 2. Generally, both the CMA and CRU datasets show that precipitation decreased from the east to the west, with the lowest precipitation occurring in the desert areas of the western arid region. The western Tianshan Mountains and the Junggar Basin (upper left corner area) are relatively humid in the arid regions. However, CRU precipitation clearly underestimates the wetness of this area (see the circle in Figure 2) compared to CMA records. Specifically, in the Tianshan area, the annual average of the gridded CMA precipitation can reach 400-500 mm or even more. In contrast, most of the gridded CRU precipitation was less than 300 mm. Moreover, the vegetation condition, as reflected by the NDVI data, further confirmed the above conclusion that the precipitation was underestimated by the CRU dataset. In particular, this can be found in the Tianshan Mountains and its western regions.  In contrast, the annual averaged temperature of CMA and CRU datasets show a better consistency than that recorded for precipitation. The main difference was also found in the Tianshan area. The annual averaged CRU temperature was a litter higher than that of CMA. For the CRU data, a limited number of observations were available in this area for interpolation which makes it difficult to generate accurate temperature reanalysis data. The annual averaged NDVI decreased from the east to the west. The annual averaged NDVI in semi-arid areas is mostly 0.2-0.4 (75.58% of the semi-arid area), whereas in the arid region, the annual averaged NDVI values were reduced to about 0-0.2 (90.71% of the arid area).
However, in the Tianshan Mountains and its western areas, the NDVI values were higher, with a magnitude greater than 0.3 (1.71% of the arid area).

Regional Averaged Wetting and Drying Trends Calculated by Nine Dryness Indicators
The drying and wetting trends of dryness indicators were further investigated in the study area as shown in Table 2 and Figure 3. For the whole dryland area, no consistent drying or wetting trend was found among the nine dryness indicators. Specifically, an increasing trend was found among six dryness indicators, i.e., SPEI-01, SPEI-06, Root_sm, Surf_sm, CMA_P and CRU_P. However, a decreasing trend was detected among another 3 dryness indicators, i.e., SPEI-12, PDSI and sc_PDSI.  In contrast, the annual averaged temperature of CMA and CRU datasets show a better consistency than that recorded for precipitation. The main difference was also found in the Tianshan area. The annual averaged CRU temperature was a litter higher than that of CMA. For the CRU data, a limited number of observations were available in this area for interpolation which makes it difficult to generate accurate temperature reanalysis data. The annual averaged NDVI decreased from the east to the west. The annual averaged NDVI in semi-arid areas is mostly 0.2-0.4 (75.58% of the semi-arid area), whereas in the arid region, the annual averaged NDVI values were reduced to about 0-0.2 (90.71% of the arid area).
However, in the Tianshan Mountains and its western areas, the NDVI values were higher, with a magnitude greater than 0.3 (1.71% of the arid area).

Semi-Arid Eastern-Arid Western-Arid
Similar trends of three sub-regions were found by comparing the annual results with the values during the growing season. However, for the entire dryland, both CMA_P and CRU_P exhibited a decreasing trend (Table 2).

Spatial Patterns of Drying and Wetting Trends by Means of nine dryness Indicators
The spatial distribution of the trends based on the indices at the grid scale is shown in Figure 4. These indices exhibit generally coherent spatial patterns. Drying trends dominate the semi-arid region, and the drying range indicated by CMA_P is slightly smaller than CRU_P (Figure 4 a-b). The Surf_sm also represents a smaller drying range than Root_sm which is significant to the ecosystem. The region with the most severe dryness is in the eastern corner. This area is a desert grasslandshrubland transition zone. Evidently, all the indices indicated a larger drying trend than other regions. In the eastern arid area, wetting trends contribute most to areas indicated by these indices except sc_PDSI. Over the western arid region, both drying and wetting areas are found. Referring to the land use and elevation distribution indicated in Figure 1, it can be found that the drying trends were mainly located in the desert area, such as the Taklimakan desert and western Gobi desert and the wetting trends were close to the Tianshan Mountains and Altai Mountains.

Spatial Patterns of Drying and Wetting Trends by Means of Nine Dryness Indicators
The spatial distribution of the trends based on the indices at the grid scale is shown in Figure 4. These indices exhibit generally coherent spatial patterns. Drying trends dominate the semi-arid region, and the drying range indicated by CMA_P is slightly smaller than CRU_P (Figure 4a,b). The Surf_sm also represents a smaller drying range than Root_sm which is significant to the ecosystem. The region with the most severe dryness is in the eastern corner. This area is a desert grassland-shrubland transition zone. Evidently, all the indices indicated a larger drying trend than other regions. In the eastern arid area, wetting trends contribute most to areas indicated by these indices except sc_PDSI. Over the western arid region, both drying and wetting areas are found. Referring to the land use and elevation distribution indicated in Figure 1, it can be found that the drying trends were mainly located in the desert area, such as the Taklimakan desert and western Gobi desert and the wetting trends were close to the Tianshan Mountains and Altai Mountains.

Drought Area Variations
Figure 5 demonstrates the regional average time series of the area in drought across each region based on SPEI-01, PDSI, and sc_PDSI. The interannual variation of the area in drought demonstrated by these three indices is generally coherent over the whole dryland and three sub-regions. A feature in the temporal patterns is that SPEI-01 changed ahead a few months of the PDSI and sc_PDSI. Based on the maximum correlation coefficient analysis, the typical time of lag is three months. Although the drylands did not show a very typical increase or decrease in drought, the drought changes more obviously in the three sub-regions. Especially in the semi-arid region, all three indicators showed significant increases in the area in drought during the second half of the study period. Specifically,

Drought Area Variations
Figure 5 demonstrates the regional average time series of the area in drought across each region based on SPEI-01, PDSI, and sc_PDSI. The interannual variation of the area in drought demonstrated by these three indices is generally coherent over the whole dryland and three sub-regions. A feature in the temporal patterns is that SPEI-01 changed ahead a few months of the PDSI and sc_PDSI. Based on the maximum correlation coefficient analysis, the typical time of lag is three months. Although the drylands did not show a very typical increase or decrease in drought, the drought changes more obviously in the three sub-regions. Especially in the semi-arid region, all three indicators showed significant increases in the area in drought during the second half of the study period. Specifically, from 1982-1996 to 1997-2012, the average area in drought is 5.73% to 16 There is no particular trend change for the eastern arid region, although there are some years of frequent droughts like the year of 2001. For the western arid region, the area of drought has slightly declined. From 1982From -1998From and 1999From -2012, the average area in drought is 8.95% to 5.20% for SPEI-01, 14.28% to 13.52% for PDSI, and 14.76% to 12.79% for sc_PDSI, respectively.

NDVI Co-variations
Usually, water availability is the primary constraint on vegetation growth over drylands and precipitation; drought index and soil moisture are widely used in the studies of water related vegetation dynamics. Here we investigated the covariations of vegetation and water by examining the correlation coefficients (r) between M-S NDVI and these dryness indicators based on time series of O-S and M-S. The regional average and gridded results are shown in Table 3 and Figure 6. There is no particular trend change for the eastern arid region, although there are some years of frequent droughts like the year of 2001. For the western arid region, the area of drought has slightly declined. From 1982-1998 and 1999-2012, the average area in drought is 8.95% to 5.20% for SPEI-01, 14.28% to 13.52% for PDSI, and 14.76% to 12.79% for sc_PDSI, respectively.

NDVI Co-Variations
Usually, water availability is the primary constraint on vegetation growth over drylands and precipitation; drought index and soil moisture are widely used in the studies of water related vegetation dynamics. Here we investigated the covariations of vegetation and water by examining the correlation coefficients (r) between M-S NDVI and these dryness indicators based on time series of O-S and M-S. The regional average and gridded results are shown in Table 3 and Figure 6. High correlations between NDVI and the multiple dryness indicators including precipitation, drought indices and soil moisture, are expected. It can be seen that for the regional averages, the highest positive correlations was found between NDVI and precipitation as well as the soil moisture, that is, the CMA_P, Root_sm and Surf_sm (Table 3) during M-S and O-S periods. However, in the western arid region, the correlations between NDVI and the dryness indicators, except the Root_sm and Surf_sm, during the growing season were lower than during O-S. Consistently, the CMA_P of growing season was also lower than the time range of O-S for the eastern-arid region. The SPEI-06 correlated higher with NDVI than SPEI-01 and SPEI-12 in the semi-arid region. Seeing the spatial pattern, 86.731% grids exhibited significant positive correlations between CMA_P and NDVI (p < 0.1). However, the correlation between CRU_P and NDVI is smaller than that between CMA_P and NDVI (Table 2 and Figure 6b), particularly over the western arid region. Although the PDSI and sc_PDSI all used the CRU precipitation as input, the correlation between NDVI and sc_PDSI was much lower than that between PDSI and NDVI over both region average and spatial patterns (Table 3, Figure 6), especially over the semi-arid region. The grids of positive correlation were only 39.735% for sc_PDSI over the semi-arid regions. As to the SPEI-01, SPEI-06 and SPEI-12, they had the similar spatial distribution (Figure 6e-g) with 56.796%, 70.227% and 42.071% grids indicating a positive correlation, which was significant (p < 0.1). The Root_sm had a higher positive correlation with NDVI than CMA_P over the eastern and western arid regions during the growing season (  High correlations between NDVI and the multiple dryness indicators including precipitation, drought indices and soil moisture, are expected. It can be seen that for the regional averages, the highest positive correlations was found between NDVI and precipitation as well as the soil moisture, that is, the CMA_P, Root_sm and Surf_sm (Table 3) during M-S and O-S periods. However, in the western arid region, the correlations between NDVI and the dryness indicators, except the Root_sm and Surf_sm, during the growing season were lower than during O-S. Consistently, the CMA_P of growing season was also lower than the time range of O-S for the eastern-arid region. The SPEI-06 correlated higher with NDVI than SPEI-01 and SPEI-12 in the semi-arid region. Seeing the spatial pattern, 86.731% grids exhibited significant positive correlations between CMA_P and NDVI (p < 0.1). However, the correlation between CRU_P and NDVI is smaller than that between CMA_P and NDVI (Table 2 and Figure 6 b), particularly over the western arid region. Although the PDSI and sc_PDSI all used the CRU precipitation as input, the correlation between NDVI and sc_PDSI was much lower

Discussion
In this study we evaluated drying and wetting trends over China's drylands based on multiple sets of indicators. Two sets of precipitation were used, the CMA and CRU precipitation. For the calculation of SPEI and PDSI drought indices, the CRU precipitation was used rather than the CMA precipitation because of computational limitations and because it is widely used. The correlations are shown in Table 4. The precipitation of CMA and CRU had similar correlations with other dryness indicators in the eastern-arid and western-arid areas. However, in semi-arid regions, the correlation between CRU precipitation and other dryness indicators was better than that of CMA precipitation with drought indices. This is largely a result of the fact that the drought indices were calculated with the CRU precipitation rather than the CMA set. Here, significant differences between CMA and CRU in precipitation and temperature data from annual average trends and spatial patterns were detected. One reason is that the CMA dataset is more detailed in the Tianshan Mountains, and the other reason is that the mean estimates are significantly larger than the CRU in the multi-year trend. Thus, in this study, all the adopted drought indices including SPEI, PDSI and sc_PDSI were calculated based on CRU precipitation. However, as the CMA precipitation product is more accurate than CRU in terms of the number of observations, it is necessary to apply the CMA precipitation for drought indices calculation in future studies within the same topic.
Based on the annual trends of various indices, spatial distribution and changes in the dry areas, we have obtained consistent conclusions to the aggravation of drought in semi-arid areas. In contrast, the eastern-arid and western-arid areas became more wet with dry areas decreasing.
Another interesting finding is that the NDVI time series correlated more with the precipitation (both CMA and CRU) than with the drought indices such as SPEI-01, SPEI-06, SPEI-12, PDSI and sc_PDSI in both regional averaged and gridded results. This is unexpected as one would expect the NDVI to be better correlated with drought indices that express something closer to water availability than just precipitation. Therefore, drought indices are generally regarded as the more appropriate indicator for describing the water conditions rather than precipitation. This is in contrast to previous studies that showed that antecedent precipitation was a significant factor to the vegetation growth in the growing season [42,43]. Above all, the application of drought indices in the study of vegetation dynamics may need further evaluations in this region. For instance, the use of the effective drought index (EDI) and standardized evapotranspiration deficit (SEDI) [44,45], which can evaluate the available water resources and detect the response of vegetation to drought at a biological time scale, can be explored in further research.

Conclusions
In this study, the drying and wetting trends of the drylands of China and its three sub-regions, i.e., semi-arid, east arid and west arid areas, were investigated by means of nine dryness indicators, including two precipitation indicators as CMA_P and CRU_P, and seven statistical and physical drought indices as SPEI_01, SPEI_06, SPEI-12, PDSI and sc_PDSI, as well as two soil moisture indicators, i.e., root zone soil moisture (Root_sm) and surface soil moisture (Surf_sm), during 1982-2012. These dryness indicators were selected based on Chen's hypothesis that robust conclusions can be obtained by using multiple and sophisticated indices at the regional scales [14]. On the whole, the drylands show a pattern of wetting in the west and drying in the east. Locally, in the sub-regions, drying trends are found in the semi-arid region in the east with the drought area significantly expanded, whereas a wetting trend was detected in the arid region, with the drought area decreasing in the west of this region but not in the eastern part.
The covariations between vegetation and water condition were also evaluated by the correlation coefficients between NDVI and each of the dryness indicators. Precipitation (both CMA and CRU sets) and soil moisture, showed significant correlation values at both regional and grid scales. This implies that vegetation pattern is mainly determined by precipitation and soil moisture. Due to the fragile vegetation ecosystem in the arid regions of northern China, it is important to select appropriate indicators to reflect the water condition for the vegetation. Furthermore, further evaluation should be performed on the rationale of applying drought indices as dryness indicators in the study of vegetation dynamics. Critically important, the results of this paper emphasized that more attention should be paid to the semi-arid region to protect the vegetation ecosystem and the water resources.
The difference between CMA and CRU records of precipitation was also quite large in the western arid region. Regarding to the vegetation spatial pattern related, it was shown that the CMA precipitation was more reliable than CRU precipitation in terms of its spatial distribution. This also implies that careful investigation of using drought indices, calculated based on CMA precipitation, for the vegetation dynamics study, is required in further research.