**3. Materials and methods**

#### **3.1. Statistical methods**

#### *3.1.1. Trend analysis*

The NDVI trend from 1982 to 2013 at pixel scale was estimated using the ordinary least squares (OLS) based on the ArcGIS 10.1 platform:

$$\Theta\_{slope} = \frac{n \times \sum\_{i=1}^{n} i \times NDVI - \sum\_{i=1}^{n} i \sum\_{i=1}^{n} NDVI\_i}{n \times \sum\_{i=1}^{n} i^2 - \left(\sum\_{i=1}^{n} i\right)^2} \tag{1}$$

where *θ* is the regression slope and *n* represents the study year during the research period. The positive value of *θ* means increasing NDVI.

#### *3.1.2. Mann-Kendall (MK) test*

Mann-Kendall analysis, applied as a nonparametric, rank-based method for evaluating trends in time-series data [11], was used to detect the changing trend because it is known as more resilient to outliers. A rank sequence (*Sk*) for time series was built:

$$LS\_K = \sum\_{i=1}^k r\_{i(k-2,3,\dots,n)}\tag{2}$$

where *k* is the dataset record length over years, and *ri* is the altered data series for original dataset:

$$\mathbf{r}\_{i} = \begin{pmatrix} 1 & \mathbf{x}\_{i} > \mathbf{x}\_{j} \\ \mathbf{0} & \mathbf{x}\_{i} < \mathbf{x}\_{j} \end{pmatrix}\_{(j-1,2,\ldots,l)} \tag{3}$$

Under the assumption of random and independent time series, the statistic *Z* is defined:

$$Z\_k = \frac{[S\_k - E(S\_k)]}{\sqrt{Var(S\_k)}} (k = 1, 2, \dots, n) \tag{4}$$

Moreover, *Z*1 = 0, *E* (*Sk*) and Var (*Sk*) is the mathematical expectation and variance, respectively:

$$E(S\_k) = \frac{n(n-1)}{4} \tag{5}$$

$$Var(S\_k) = \frac{n(n-1)(2n+5)}{72} \tag{6}$$

The positive *Zk* value means the trend is increasing. Compared *Zk* with *Zα*, the result of |*Zk*| > *Zα* (*Z*0.05 = 1.96) means the trend is statistically significant.

#### *3.1.3. Ordinary linear square*

**3. Materials and methods**

32 Land Degradation and Desertification - a Global Crisis

(OLS) based on the ArcGIS 10.1 platform:

q

The positive value of *θ* means increasing NDVI.

*3.1.2. Mann-Kendall (MK) test*

dataset:

=

The NDVI trend from 1982 to 2013 at pixel scale was estimated using the ordinary least squares

1 1 1

*n i NDVI i NDVI*

å åå

*i i i*

*n n n*

<sup>=</sup> = =

*slope n n*

resilient to outliers. A rank sequence (*Sk*) for time series was built:

where *k* is the dataset record length over years, and *ri*

´´ -

2

*ni i*

1 1

= =

*i i*

where *θ* is the regression slope and *n* represents the study year during the research period.

Mann-Kendall analysis, applied as a nonparametric, rank-based method for evaluating trends in time-series data [11], was used to detect the changing trend because it is known as more

( 2,3, , )

( 1,2, , )

*i j j i*

1 *k K ik n i S r* = ¼ =

*i j*

Under the assumption of random and independent time series, the statistic *Z* is defined:

[ ( )]( 1, 2, , ) ( )

*x x <sup>r</sup> x x* = ¼ æ ö ç ÷ è ø

1 0

*k k*

*Var S*

*k S ES Z kn*

*i*

*k*

´ - ç ÷

å å

2

*i*

<sup>=</sup> å (2)

<sup>&</sup>gt; <sup>=</sup> <sup>&</sup>lt; (3)


is the altered data series for original

(1)

æ ö

è ø

**3.1. Statistical methods**

*3.1.1. Trend analysis*

In order to compare the relative importance of temperature and precipitation for NDVI, the multivariate regression and the standardized coefficients were applied together. The higher standardized values mean important roles. The MATLAB 8.1 was used to establish multivariate linear model:

$$NDVI = b\_0 + b\_1 \times Temperature + b\_2 \times Prexipation + \varepsilon \tag{7}$$

where *b*0, *b*1, and *b*<sup>2</sup> are the regression parameters, while *ε* is the regression residual. Because of the different range for values of temperature and precipitation, it required normalization to compare the relative importance of climatic factors in the NDVI variations:

$$b\_i^\cdot = b\_i \times \sqrt{\sum\_{\iota=1}^n \left(\chi\_\iota - \overline{\chi}\right)}\tag{8}$$

#### *3.1.4. Geographically weighted regression (GWR)*

The GWR analysis, coupled in ArcGIS 10.1, was conducted to reveal the spatial variations in relationships between NDVI and climatic variables. Both the spatial distribution and the dynamics of NDVI were considered by the GWR model. GWR extends the traditional OLS to consider the spatial heterogeneity in climate-vegetation correlations by assigning weight values [12]:

$$\mathbf{y}\_i = \beta\_0(\boldsymbol{\mu}\_i, \boldsymbol{\nu}\_i) + \sum\_{k=1}^p \beta\_k(\boldsymbol{\mu}\_i, \boldsymbol{\nu}\_i)\mathbf{x}\_{ik} + \boldsymbol{\varepsilon}\_i \tag{9}$$

where *yi* , *xik*, and *εi* , represent the dependent variable, the independent variables, and the random error term at location *i*, respectively. Note that (*µi* , *νi* ) expresses the coordinate location of the *i*th point, *k* denotes the independent variable number. *β*0 and *βk* are the regression parameters at location *i*.

The regression coefficients were estimated by:

$$\beta(\mu\_i, \nu\_i) = (X^\top W(\mu\_i, \nu\_i)X)^{-1} X^\top W(\mu\_i, \nu\_i)Y \tag{10}$$

*β* is the unbiased estimate of the regression coefficient. *W* is the weighting matrix, and *X* and *Y* are matrices for independent and dependent variables, respectively. The kernel function, used to determine the weight, was performed as the exponential distance decay:

$$\alpha o\_{\boldsymbol{\upbeta}} = \exp \left( \frac{d\_{\boldsymbol{\upbeta}}^2}{b^2} \right) \tag{11}$$

*ωij* expresses the weight of observation *j* for location *i*, *dij* represents the Euclidean distance between points *i* and *j*, and *b* is the kernel bandwidth.

#### **3.2. WRF climate model and experimental design**

The WRF-ARW was developed as the next generation for regional climate model. It includes different parameterization schemes for longwave and shortwave radiation, cloud microphysics, cumulus, and land surface processes. The simplified simple biosphere model (SSiB), coupled with WRF model, was selected to simulate land surface energy balance. According to the SSiB model description, there are 12 types of vegetation cover, while the vegetation and soil parameters were set for every types. Defining different vegetation cover types in this study enabled investigation of the impact of land degradation and Karst rocky desertification using the WRF-SSiB model. The domain for WRF model was set as follows: dimensions of 196 × 154 horizontal grid points with center at 35°N and 110°E. In this domain, the influencing factors for East Asian summer monsoon can be included, for example, the upper level westerly jet (ULJ) and low-level jet (LLJ), the Bay of Bengal and the southeast trade wind, and so on [13]. The WRF downscaling ability was assessed by comparing the simulations with different physical schemes (**Table 1**), and the optimal combination was concluded from the assessment. For the execution of the WRF, we used the NCEP DOE Reanalysis-2 [14], hereafter NCEP R-2, at 6-h intervals to provide initial conditions and lateral boundary conditions.

Two experiments were done. One was the Case C, using the original SSiB vegetation map (as shown in **Figure 2a**), the other was Case D with the degraded land cover types (**Figure 2b**). The degraded types were decided based on the spatial pattern of different rocky desertification degrees [15]. For example, if the deserted areas accounted more than 30% of the corresponding counties, the SSiB vegetation was modified to bare soil (type 11 in SSiB model). The type 9 (broadleaf shrubs with bare soil) was used to replace original vegetation types in areas described as desert and potential desert areas larger than 45% of the counties and smaller than 30% of the counties, respectively. Based on the reset of vegetation cover types, two vegetation maps were used in WRF model, and was further used to conduct Case C and Case D.


R, correlation coefficient; RMSE: root mean square error.

where *yi*

, *xik*, and *εi*

34 Land Degradation and Desertification - a Global Crisis

parameters at location *i*.

random error term at location *i*, respectively. Note that (*µi*

bmn

between points *i* and *j*, and *b* is the kernel bandwidth.

**3.2. WRF climate model and experimental design**

The regression coefficients were estimated by:

, represent the dependent variable, the independent variables, and the

 mn


è ø (11)

of the *i*th point, *k* denotes the independent variable number. *β*0 and *βk* are the regression

<sup>1</sup> (,) ( (,)) (,) *T T*

*i i XW X XW Y i i i i*

*β* is the unbiased estimate of the regression coefficient. *W* is the weighting matrix, and *X* and *Y* are matrices for independent and dependent variables, respectively. The kernel function,

2

*d b*

*ωij* expresses the weight of observation *j* for location *i*, *dij* represents the Euclidean distance

The WRF-ARW was developed as the next generation for regional climate model. It includes different parameterization schemes for longwave and shortwave radiation, cloud microphysics, cumulus, and land surface processes. The simplified simple biosphere model (SSiB), coupled with WRF model, was selected to simulate land surface energy balance. According to the SSiB model description, there are 12 types of vegetation cover, while the vegetation and soil parameters were set for every types. Defining different vegetation cover types in this study enabled investigation of the impact of land degradation and Karst rocky desertification using the WRF-SSiB model. The domain for WRF model was set as follows: dimensions of 196 × 154 horizontal grid points with center at 35°N and 110°E. In this domain, the influencing factors for East Asian summer monsoon can be included, for example, the upper level westerly jet (ULJ) and low-level jet (LLJ), the Bay of Bengal and the southeast trade wind, and so on [13]. The WRF downscaling ability was assessed by comparing the simulations with different physical schemes (**Table 1**), and the optimal combination was concluded from the assessment. For the execution of the WRF, we used the NCEP DOE Reanalysis-2 [14], hereafter NCEP R-2,

<sup>2</sup> exp *ij*

æ ö = = ç ÷ ç ÷

 mn

used to determine the weight, was performed as the exponential distance decay:

*ij*

at 6-h intervals to provide initial conditions and lateral boundary conditions.

Two experiments were done. One was the Case C, using the original SSiB vegetation map (as shown in **Figure 2a**), the other was Case D with the degraded land cover types (**Figure 2b**). The degraded types were decided based on the spatial pattern of different rocky desertification degrees [15]. For example, if the deserted areas accounted more than 30% of the corresponding counties, the SSiB vegetation was modified to bare soil (type 11 in SSiB model). The type 9

w

, *νi*

) expresses the coordinate location

**Table 1.** Descriptive statistics of precipitation and temperature from WRF/SSiB with different microphysics and radiation schemes for June 2000 over 18°-52°N, 86°-136°E.

**Figure 2.** Potential LCC based on the spatial pattern of KRD in GKP. (a) The percentage of areas with KRD for counties. (b) SSiB vegetation map for GKP and vegetation cover conversion.
