Drought level definition based on PAI.
It is important to highlight energy-water balance and ecosystem response to climate changes. The change of water-energy balance and ecosystem due to climate change will affect the regional ecological and human living significantly, especially in Southwest China which is an ecologically fragile area. This chapter presents the retrieval methodology of parameters (reconstruction of vegetation index, land cover semi-automatic classification, a time series reconstruction of land surface temperature based on Kalman filter and precipitation interpolation based on thin plate smoothing splines), time-series analysis methodology (land cover change, vegetation succession and drought index) and correlate analysis methodology (correlation coefficient and principal component analysis). Then, based on the above method, remote sensing data were integrated, a time series analysis on a 30-year data was used to illustrate the water-energy balance and ecosystem variability in Southwest China. The result showed that energy-water balance and ecosystem (ecosystem structures, vegetation and droughts) have severe response to climate change.
- energy-water balance
- climate change
- droughts estimation
- vegetation index
- time series analysis
- land cover change
Southwest China consists of the municipality of Chongqing and the four provinces of Guangxi, Yunnan, Guizhou and Sichuan, 21°07′N to 34°14′N, 97°30′E to 112°05′E (Figure 1). Southwest China, which was considered to be rich in water resources, has suffered climate changes in the last decades. It is particularly important to highlight energy-water balance and ecosystem response to climate changes in Southwest China using the remote sensing techniques and Geographic Information System(GIS). This chapter introduced the methodologies for retrieving parameters of land cover, vegetation index, gridded precipitation, land surface temperature, which can be used as the factors of climate change and ecosystem. Furthermore, the climate change response of droughts estimation methods was presented. The time series analysis of water-energy balance and ecosystem variability of Southwest China in recent 30 years showed that ecosystem structures, vegetation and droughts were affected by the change of climate. The following paragraphs describe the recent achievements in more detail.
1.1. Retrieval of parameters of energy-water balance and ecosystem
1.1.1. Ecosystem (vegetation index): reconstruction of vegetation index based on Savitzky-Golay filter
Terrestrial ecosystem is extremely sensitive to climate change, especially the corresponding change of the surface vegetation is the most significant. Vegetation is the important feature of the land surface and is the core and function parts of the biosphere and ecological system. It is widely accepted that energy and water together drive diversity and form of vegetation [1, 2]. In order to show the vegetation information, many kinds of vegetation index (VI), such as Normalized Difference Vegetation Index(NDVI), Soil-Adjusted Vegetation Index(SAVI), Transformed Soil-Adjusted Vegetation Index(TSAVI), Modified Soil-Adjusted Vegetation Index(MSVAI), Difference Vegetation Index(DVI), Green Vegetation Index (GVI), Perpendicular Vegetation Index(PVI), Enhanced Vegetation Index(EVI), etc., have been developed . And most of these attempts focused on the differences of the absorptive and reflective electromagnetic spectrum properties between the visible (VIS) and near-infrared (NIR) portions. According to some researchers, the time series Vis can derived from NOAA/AVHRR, SPOT/VEGETATION, TERRA or AQUA/MODIS, and the time series Vis can detect long-term development of vegetation, evapotranspiration, drought, plants phenology, corn yield, land-use/cover changes, and terrestrial ecosystems at different spatial and temporal resolutions globally [4–7]. Theoretically, because of the subtle vegetation canopy changes with respect to time, a generalized VI temporal profile is continuous and smooth . However, the time series VI data always fluctuate with remarkable rises and falls which is the result of disturbances possibly caused by cloud contamination, atmospheric variability and bi-directional effects. Besides, there are a lot of no-data pixels in some VI data products because of hardware or human factors . For instance, the VI datasets of MOD13A2 of MODIS need cooperate with its quality assurance (QA) data layer for application. The reliability and quality of VI values can be checked from the QA data layer. A series of approaches have been developed to reduce the errors of the noises in VI data and to reconstruct high quality time series datasets. Among them, a simple algorithm based on Savitzky-Golay filter explored by Chen was thought to be most efficient .
To reconstruct the high-quality MODIS EVI time series data, the Savitzky-Golay filter was used based on two hypotheses proposed by Chen: (1) The EVI data from a satellite sensor are primarily related to vegetation changes. Such as, an EVI time series follows annual cycle of vegetation growth and decline; (2) Clouds and poor atmospheric conditions made an EVI time series incompatible with the gradual process of vegetation change, because these usually depress EVI values and cause sudden drops in EVI which are regarded as noise and removed .
The Savitzky-Golay filter process requires continuous data, while the MOD13A2-EVI is not continuous as it contains poor quality pixels (for example, cloudy, not being processed). Based on the quality assurance (QA) dataset generated in preprocessing of MODIS EVI data in spatial, we produce a continuous dataset by interpolating the poor quality pixels of EVI data. Firstly, we searched eight pixels with good QA nearest a given pixel with poor quality and recorded the value of those eight pixels and their distance to the given poor quality pixel. Figure 2 describes the eight pixels by eight different directions searching.
After that, the inverse distance weighted interpolation method (IDW) is applied to compute the values of the given pixel. According to the definition of IDW, we can define EVI as:
where EVI is EVI value of poor quality pixel after interpolation and EVIi is the EVI value of a good pixel. The continuous initial EVI time series is recorded as EVI0.
Each pixel in the 163 images in seven years was studied using time series analyzing methods. Based on the EVI0 data, we firstly generated a time series curve of each pixel and dealt with noisy signal for each curve using Savitzky-Golay filter, and then extracted a long-term change trend (EVItr). After that, the weight of each point, which was in the 163 samples time series () data, was computed and the EVI time series were traversed.
1.1.2. Ecosystem (land cover): land cover semi-automatic classification from multispectral remote sensing imagery
Through the conversion of forests and grasslands to croplands and pastures, humans have affected the exchange of energy, water and carbon between the atmosphere and the land surface. As the key input to ecosystem researches, various studies focused on land cover mapping . Whereas, land cover data at large scale is hard to approached except remote sensing . A variety of approaches was used to map land cover based on remote sensing data [13–17], including visual interpretation classification, unsupervised clustering coupled with extensive ancillary data and manual labeling of clusters, supervised classification, expert system classification, artificial intelligence neural network classification, and decision tree classification. However, the accuracy and the efficiency of land cover classification is not guaranteed and the land cover classifications are arbitrary using these method. For example, supervised classification methods require selecting training samples which rely on substantial expertise and human participation, so the result of land cover classification is influenced greatly, and it is impossible to classify land cover automatically with these methods. The algorithms such as neural network classification and fuzzy logic classification are difficult to understand and apply widely because of highly complicated in their algorithm basis. The construction of the decision tree and the assignment of thresholds for each sub-nodes are the key problem of decision tree classification, and they heavily depends on human experience and varies spatially and temporally. In order to solve these problems to improve the accuracy and efficiency of classification, Jiang et al. proposed an efficient automatic landscape classification approach using prior accurate land-cover data as the background experience [18, 19]. This method consists of two steps: (1) semi-automatically detecting land cover changed pixels from satellite images compared with prior land cover map; (2) semi-automatically classifying the land cover of changed pixels based on pattern recognition and changed rules.
Automatic collection of training samples: In this method, pure pixels of land cover were extracted automatically with an accurate previous land cover dataset as prior knowledge. The interiors of individual land cover areas and larger patches are considered to be more ecologically stabile areas. Based on accumulation area threshold (Pa), the samples of different land covers sorted in descending order of their area. The accumulation area threshold is calculated based on the percent of largest patches of specific land cover categories occupied in its total area.
After spatial buffer analysis, the joint region of different land cover were discarded. The buffer area are obtained using the different distance to all patches because of patches of vary areas. The distance of buffer analysis is
where is area threshold for buffer analysis, is the buffer area of the patch with a distance of d, with d negative, and A is the area of the patch.
Establishment of three-dimensional feature space: The data in all spectral bands of each land cover class, which were extracted from the region of interest, were processed by principal component statistical analysis. The first three principal components for orthogonal decomposition was selected to construct the three-dimensional feature space of different land cover classes.
where P is the principal component (PC), MP and are mean and standard deviation of PC, respectively, c is the radius of the three-dimensional feature space.
Change detection and classification of changed pixels: Combining the satellite images and early land cover maps, the spectral data of the images were extracted according to accurate early land cover maps. Based on spatial intersect analysis with corresponding three-dimensional feature space, the pixels with spectral data outside the ellipsoid were detected as changed pixels. After obtaining the changed land cover pixels, the satellites images and three-dimensional feature space were used to classify the changed area of land cover by calculating the minimum spectral distance. For each changed pixel, the spectral data of all bands were input to the formula for all land cover classes in three-dimensional feature space to calculate the minimum spectral distance .
To express the rules of changed land cover, we proposed a drag coefficient of changed land cover (r). Combining and r, we determined the final land cover classification of changed pixels. The minimum distance of the land cover classification based on changed rules is defined as:
where is the minimum distance of the mth pixel to the ith land cover class based on the changed rules and is the drag coefficient of the ith land cover class that changed to the jth land cover class.
The class information would be contaminated by adjacent class codes, so the post-classification results was modified to solve the problem. Any group of pixels, which was <3 pixels in size, was identified as noise in dilation operation.
1.1.3. Energy parameter: land surface temperature
The parameter of temperature of surface energy balances has be important index for studies of ecosystem response to climate change. As land cover mapping, the land surface temperature retrieval based on remote sensing has been the primary method for gaining temperature at large scale. This is due to the accessibility and convenience of satellite remote sensed information. However, the surface energy parameters retrieved from remote-sensing data are often interrupted on the spatial and temporal scales by clouds, aerosols, solar elevation angle and bidirectional reflection. so, it is always difficult to obtain complete dataset for a large region. In addition, the accuracy of surface parameter retrieval was affected with various degrees by an indirect retrieval method and the instantaneous features of monitoring [2, 20]. To reduce such impacts, a time-composite method is generally adopted. For surface energy parameters, the time series fitting and noise-removal methods commonly used to include the mean diurnal variation and nonlinear regression methods. Compared with nonlinear regression method, in-situ measurement data are necessary for the mean diurnal variation method [10, 21]. However, the modeled results often cannot represent the actual situation under significant changes of environmental conditions [20, 22]. In recent years, data-assimilation methods have been adopted in the reconstruction of time series data to solve these problems. In this paper, a time series reconstruction method based on an ensemble Kalman filter was established, which is can be applied in various studies[19, 23–26]. It focused on evaluating the state of a discrete-time controlled process. The state Xk can be expressed as :
where Wk−1 are noise in normal probability distributions; and are the noise covariance of process and measurement, respectively.
Further assumption was defined as:
where is the nonnegative covariance matrix of and is the positive covariance matrix of ; is the function of Kronecker-δ.
We define N as the number of days of measurement, then:
where is defined to be our a priori state estimate at step k, given knowledge of the process prior to the step k, and is our a posteriori state estimate at step k, given measurement . Then, the a priori state is defined as follows:
and the a posteriori state as follows:
The calculation process of a Kalman filter is a constant “forecast correction” process. Based on time-update and observation-update, the values were reconstructed with minimum variance by compared to in situ values.
1.1.4. Water balance: precipitation interpolation based on thin plate smoothing splines
Precipitation data are one of the important input elements of ecological mechanism model, and it play an extremely important role in simulating and researching ecosystem changes at regional or global scale. It also a key index in global change researches, which is a direct parameter for drought estimation . The precision of precipitation data may affect the simulation and prediction precision. At present, the precipitation data acquisition mainly rely on the long-term observation of the weather stations, but scarcity and uneven distribution of weather stations reduce the accurate of spatial precipitation data.
The precipitation measurements from the WMO stations are point measurements. Spatial interpolation technique for meteorological elements is often used to obtain the meteorological data of every position in the scope which is covered by the meteorological stations. Several methods have been developed to interpolate these point data to a real estimation of rainfall. The software ANUSPLIN, developed by Hutchinson, Australian National University, was applied to generate gridded precipitation data . ANUSPLIN is a kind of analysis and interpolation tool for multivariable data using local thin plate spline function. The model required to enter the location of the meteorological site, elevation and other ancillary data. Its statistical analyzes and diagnose multivariable data to realize spatial interpolation function.
In order to fit to datasets distributed across an unlimited number of climate station locations, interpolating methods of ANUSPLIN were applied, which use thin plate smoothing splines for spatial interpolation with a third parameter of elevation [30–36]. The main advantage of thin plate splines is that splines do not require prior estimation of spatial auto-covariance structure, which was difficult to estimate and validate. The partial spline observational model for n data values zi at positions x is given by setting :
where is the thin smoothing spline, are jth unknown parameter and are jth known function, which all have to be estimated. is spatial position with elevation, are errors with covariance structure given by
where , is positive definite matrix and may be known or unknown, the errors are uncorrelated if is diagonal and correlated otherwise. The function and the parameters are estimated by minimizing
where and with
where is a measure of the roughness of the spline function defined in terms of mth order derivatives of , and is a positive number called the smoothing parameter
1.2. Time series analysis
1.2.1. Land cover change
Land cover change is one of the main methods by which the human activities have effect on the land surface environment. Research on dynamic models of land use change process is an important approach and means to deeply understand land use change process and its causes, and the research also can reveal the response of the land cover to anthropogenic activities. Land cover change mainly shows in land cover change speed and transfer direction, and comprehensive land use dynamic degree, single land use dynamic degree and transition matrix were used to express them .
2. Comprehensive land cover dynamic degree
Comprehensive land use dynamic degree is used to describe regional difference of land use types change speed and reflect the influence on change of land cover types by human activities. The mathematical model are as follows:
where is the comprehensive land cover dynamic degree in t period, is the total area of land cover change of ith type converted into jth type from monitor beginning to end. is the area of ith type when the monitor started.
3. Single land cover dynamic degree
Single land use dynamic degree is used to describe the speed and amplitude of different land cover types change in a certain period. It reflects the influence on change of single land cover type by human activities. The mathematical model is as follows:
where is the single land cover dynamic degree of ith type from t1 to t2 time, and were areas of ith type in and time separately.
4. Transition probability matrix
Transition probability matrix was proposed by the Russian mathematician Markov. At the beginning of the twentieth century, Markov found that the nth result affected by the n-1th result in the transfer of some factors of a system. In Markov’s analysis, the quantitative descriptions of the system state and state transition were reflected in the transform process of a metastable system from time T to time T + 1 in a certain time interval, thus revealing the land cover pattern time and space evolution process. Transition matrix of land cover depicts comprehensively and specifically structural characteristics of land cover change and reflects the change direction led by human activities. Transition matrix are as follows:
4.1. Vegetation succession
Net primary production (NPP) represents the accumulated organic matter by plants per unit area and time. From an ecological perspective, it measures the rate at which solar energy is stored by plants as organic matter. NPP is influenced by climate, soil, vegetation type and human activities, for various ecological monitoring activities, and is generally regarded as an important factor that provides a comprehensive evaluation of ecosystem status and services, including productivity capability, habitat, and wildlife, and ecological footprint [39, 40]. NPP is not a directly observable ecosystem characteristic, and it is difficult to measure accurately over large areas due to the spatial variability of environmental conditions. A number of NPP models for different ecosystems have been developed. These models are broadly classified into regression-based and process-based. Regression-based models are established by empirically derived relationships between climate values and NPP, such as Miami . Although regression-based models, with the advantages of simplicity and fewer parameter requirements, can be extrapolated for most land ecosystems, uncertainties are also involved when considering heterogeneous vegetation, standard errors of measurements and novel climatic conditions, which may not be appropriate for the regressions [42, 43]. Process-based models, ranging from simple models based on light use efficiency (LUE) to more mechanistic models based on “soil-vegetation-atmospheric-transfer” (SVAT) schemes, are based on physiological and ecological processes such as photosynthesis, evapotranspiration, respiration and nutrient cycling [44, 45]. These models have more parameter requirements and complexities; however, they better describe mechanisms and have the potential to estimate NPP more accurately when compared with regression-based models. The models based on LUE are called production efficiency models (PEMs), which use LUE for the conversion of absorbed photosynthetically active radiation (APAR) to biomass . They are widely acceptable to map NPP at different scales as it follows the basic principles of the photosynthesis process and is easily amenable to remote sensing data . The satellite data-driven PEMs, such as CASA  and GLO-PEM , have been used to analyze the spatiotemporal patterns of NPP over continents and global land surfaces [50–53].
The CASA model simulates NPP directly thus avoiding a Ra (autotrophic plant respiration) calculation and taking environmental conditions (temperature, rainfall/soil moisture) and vegetation characteristics into consideration [54, 55].
where x represents the grid cell and t represents the period in which NPP is accumulated, for example, a month. APAR is determined by the fraction of photo synthetically active radiation (FPAR) and the total solar surface radiation (SOL) (MJ m−2)  as
where the constant 0.5 represents the ratio of the total solar radiation (with a wavelength range of 0.4–0.7 μm) used by the vegetation .
LUE is calculated as the product of maximum light use efficiency, and its temperature and moisture stressors  as
where represents the actual light use efficiency, the maximum light use efficiency, and the value for grass (0.604 g/MJ), simulated by Running based on BIOME-BGC model , was used here; and are temperature scalars, and is the moisture stress coefficient. , and were computed at every location at each time step. and are calculated as [56, 57].
where is an optimal temperature, defined as the mean temperature in the month of maximum normal differential vegetation index (NDVI). T is the monthly mean temperature; , when ; it decreases to 0.5 when T is 10°C above or 13°C below .
reflects the effect of water condition, and it generally increases when available water increases. Atmospheric vapor pressure deficit reflects air humidity, which affects transpiration and then the LUE . Therefore, there are currently studies using vapor pressure deficit (D in kPa) to calculate the moisture stress coefficient [61, 62], computed as 
where is dew point temperature (K) and is surface temperature (K). When . Td was derived from Guo Jie’s regression model for Sichuan province based on Yang Jingmei’s findings of a significant linear relationship between dew point temperature and the logarithm of total perceptible water [63, 64] as follows:
where U is total perceptible water (mm).
4.2. Estimation of droughts: drought index, drought level definition, index of drought frequency
Drought is an kind of extreme water deficit processes, which can be used as indicators of ecosystem deteriorate and climate change. Drought can be classified to four categories of Meteorological Drought, Hydrological Drought, Agricultural Drought and Socioeconomic Drought [65–68]. There are numerous drought indices were formulated by integrating variables to identify and quantify the duration, magnitude, intensity and spatial extent of a drought, such as precipitation, evapotranspiration, temperature, terrestrial water storage (TWS), the TWS anomaly index (TWSI), vegetation, etc.[69–71]. Whereas, precipitation is the most direct parameter for evaluating meteorological droughts, which also can be applied in estimating agricultural drought, hydrological drought and socioeconomic drought [72–73]. To be a worldwide natural hazard, meteorological droughts can be measured by various indices such as Precipitation Anomaly Index (PAI) [69, 74, 75], Palmer Drought Severity Index(PDSI) , Z-score or standardized rainfall anomalies , Standardized Precipitation Index(SPI), Standardized Precipitation Evapotranspiration Index(SPEI) , et al. However, compare with complex indices, a simple measure may be applied more easily to evaluate drought disaster at large scale . For example, PDSI requires, in addition to precipitation, soil moisture, runoff, evapotranspiration, potential evapotranspiration and other factors of plant growth to assess droughts. Furthermore, PDSI cannot be used to identify drought at short time scales, e.g., less than nine months . Compared to complex indices, such as PDSI, SPI maybe a better index based on precipitation alone, as it also compares drought conditions among different time periods and regions. Among these indices, precipitation anomaly index (PAI) that uses precipitation alone is the simplest index; it is a dimensionless number in which negative/positive values indicate dry/wet conditions. It is precisely because of these advantages of simple computation, spatiotemporal consistency, and easy comparison to historical records, PAI is an important meteorological drought index for large area drought assessment in China .
The meteorological droughts index of the precipitation anomaly index (PAI) was used for drought analysis. PAI was calculated from the monthly precipitation data from the China meteorological data sharing service system (CMDSSS) of the China meteorological administration (CMA). The PAI of SC from 1961 to 2012 is calculated as:
where and are precipitation and mean value.
|Droughts level||PAI (%)|
|None||[−40, )||[−25, )||[−15, )|
|Mild||[−60, −40)||[−50, −25)||[−30, −15)|
|Moderate||[−80, −60)||[−70, −50)||[−40, −30)|
|Severe||[−95, −80)||[−80, −70)||[−45, −40)|
|Extreme||(− , −95)||(− , −80)||(− , −45)|
Because of difficulty to evaluate absolute droughts, the meteorological drought level was used to estimate the drought severity. Table 1 shows the relationship between the PAI and meteorological drought level in three temporal scales of month season and year .
The drought frequency (DF) were defined as:
where n is the number of years of droughts, N is number of study years .
4.3. Correlate analysis between energy-water balance and ecosystem parameters
4.3.1. Correlation coefficient
Correlation coefficient method is used for studying the closeness relations of variables and is described by a quantitative index which is called correlation coefficient. The calculation process is simple and clear and the result is intuitional and easy to interpret, so the method is considered to be the best for analyzing long-term vegetation trends. When we study the interrelation of a plurality of geographic features, and study the impact of certain factor on the other feature without taking into account of other features, the correlation coefficient can be used. We chose correlation coefficient as the quantitative indicator of evaluation relevant. The formula of correlation coefficient is as follows :
where is the correlation coefficient of variable X and variable Y; n is the number of sample; is the mean of variable X; is the mean of variable Y.
The significance test of correlation coefficient generally uses the t-test. The statistic calculation formula is:
where r is the correlation coefficient, n is the number of samples.
4.3.2. Principal component analysis
Principal component analysis (PCA) is a kind of multivariate statistical method. Through orthogonal transformation converts, a set of possible correlation between variables convert to a set of linear irrelevant variables which is called principal components.
Principal component analysis is to investigate correlation among multiple variables. It was used to study that how to reveal the internal structure among multiple variables by a few principal components. Namely, a few principal components were derived from the original variables and make them as much as possible to retain the information of the original variables and unrelated to one another. Usually, in mathematical treatment, the linear combination of the original index was regard as a new composite indicator. The classical approach is to use the variance of F1 (first linear combination) to express. The bigger Var(F1) is, the more information you gather F1 contains. Thus, the variance of F1 should be the biggest of all linear combinations and F1 is the first principal component. The number of principal components is decided by the quantity information that principal components represented .
5. Result and discussion
5.1. Climate change pattern of Southwest China
5.1.1. Land surface temperature
Taking Sichuan province as an example, the temperature increased from 1982 to 2010 (Figure 3). In regional terms, the temperature increased in 90.2% of the region in Sichuan province, in addition to Batang County, Derong County and Xiangcheng County. In western Sichuan region, Annual average temperature showed a rising trend, and it is very obvious in Jiuzhaigou County, Li County, Mili Tibetan Autonomous County, Shiqu County and Leibo County, the biggest increase among them is 1.89°C every 10 years. In eastern Sichuan, including Guangyuan, Bazhong, Nanchong, Ziyang, Leshan, Luzhong and so on, temperature increased by over 0.2°C every 10 years. The temperature rising trend in 76% of the region of the Sichuan province are extremely significant or significant. The temperatures have a decreased of as much as 1.25°C every 10 years in Batang County, Derong County and Xiangcheng County.
Figure 4 shows the change of mean temperature of Sichuan Province in four seasons from 1982 to 2010. It turns out the spatial distribution pattern changes of mean temperature at different seasons are not significant. The changes of mean temperature in summer ranged from −6.59 to 22.06°C, the change ranged from −6.59 to 22.06°C in spring, from −6.59 to 22.06°C in autumn, and from −6.59 to 22.06°C in winter. The mean temperature variations in winner is biggest (30.746°C), and it in autumn is smallest (25.27°C).
Figure 5 shows the monthly mean temperature changes of time series from 1982 to 2010. As can be seen from the Figure 5, the monthly mean temperature display cyclic variations under this time series. It offered upgrade firstly than descending latter tendency. In June, July, August and September, the temperature was higher and reaches a maximum in July and August. In January, the monthly mean temperature is lower. The changes of every year were basically similar, and they have certain periodicity with an obvious fluctuation. The fluctuation of the highest temperature of each year is very small. A maximum of the highest temperature appears in the June 2006 was 23.65°CA minimum of the highest temperature has an obvious fluctuation ranged from 1.23(1984) to 4.28(1987).
Figure 6 shows annual changes of the mean temperature from 1982 to 2010. The mean temperature fluctuated between 12.48 and 13.90°C. The minimum of the mean temperature occurred in 1992 and the maximum occurred in 2006. From 1982 to 2010, the annual mean temperature was rising obviously which was consistent with global warming trends.
Taking Sichuan province as an example, we selected the monthly mean precipitation data of nine weather stations from 1982 to 2010. As can be seen from Figure 7, the spatial distribution of annual average precipitation was uneven in 29 years with descending from the east to the west. Taking Songpan Country-Li Country-Kangding City-Jiulong Country-Yanyuan Country as a boundary, the annual precipitation of the east of the boundary was abundant, and it is more than 900 mm/a. Especially, the annual mean precipitation of Ya’an is more than 1200 mm/a and ranked first in the whole province. The annual mean precipitation of most areas in the west is not more than 800 mm. Among them, there is a very little rain in Shiqu country, less than 600 m/a.
As shown in Figure 8, the annual average precipitation was concentrated mostly in summer, and there are obvious differences in the spatial distributions in different seasons. The spring precipitation mainly concentrated in the central and eastern regions in Sichuan, and the maximum precipitation was 311.55 mm. The precipitation in summer moves westward concentrating in the central regions with the maximum precipitation of 859.68 mm. The scope of autumn rainfall mainly concentrated in the central and southern Sichuan with uneven spatial distribution. The maximum precipitation was 340.02 mm. The precipitation of most regions in winter was low, and decreased from the west to the east.
Monthly mean precipitation changing trends in every years were basically the same, namely, it offered upgrade firstly than descending latter tendency (Figure 9). In June, July, August and September, the precipitation was higher and reaches a maximum in July and August. In January, the monthly mean precipitation is lower. The precipitation change has certain periodicity with an obvious fluctuation from 1982 to 2010. The maximum of annual highest precipitation was 294.49 mm in 1984, and the minimum is 154.03 mm in 2006. The annual lowest precipitation has no obvious fluctuation ranged from 2.86 to 15.71 mm. The monthly mean precipitation had not a visible downtrend by 0.11 mm per decade.
Figure 10 shows the changes of the annual precipitation from 1982 to 2010. The annual precipitation showed obvious decreasing trend. The maximum of the annual precipitation occurred in 1989 (1131.17 mm), and the minimum occurred in 2006 (812.84 mm). The precipitation changed in a wide range and had obvious fluctuations.
5.2. Variability of ecosystem in Southwest China
5.2.1. Land cover change
We integrated remote sensing data and other correlative material to get the land use maps for 1980, 1990, 1995, 2000 and 2005 years (Figure 11). The maps included 18 land cover types: dry land, paddy fields, low coverage grassland, medium coverage grassland, high coverage grassland, sparse woodland, shrub land, woodland, other wood land, lakes, graff, reservoir and pits, coast, beach, urban land, rural residential areas, construction land, waste land. As can be seen from Figure 11, the main land use types were cropland, grassland, woodland. Among them, paddy field and dry land were highly concentrated in Sichuan and Chongqing and dry land accounted for a large proportion. There are only small numbers of them scattered across the other region. The grassland, including low coverage grassland, medium coverage grassland and high coverage grassland, were distributed mainly in Sichuan and Guizhou. The woodland were mainly distributed at Yunnan and Guangxi.
The time series analysis on land cover indicated the whole change is not obvious with certain changes on land use types during 20 years. High coverage grassland, woodland, urban land and construction land increase, others decrease.
5.2.2. Vegetation destruction and recovery
The normal differential vegetation index (NDVI) products MOD13A2 (d209-d225) were obtained from the NASA website (ftp://e4ftl01.cr.usgs.gov/) at the same time from 2004 to 2010. These data were used to analyze the vegetation destruction and recovery. As can be seen from Figure 12, NDVI is smaller in the northwest of Sichuan Province than other area in Southwest China. Comparing 2006 with 2004, NDVI in southeast of Yunnan and Chongqing Province changed significantly with a sudden drop, and it in other areas showed no obvious change. In 2008, the conditions had improved, NDVI in southeast of Yunnan and Chongqing Province increased, but it in the northwest of Sichuan Province decreased. In addition, there is a large scale of the decline of NDVI in the south-central part of Southwest China. NDVI in some regions of northwest of Sichuan Province increased in 2010, and it decreased in the region near the boundary of Sichuan and Chongqing Province. In general, NDVI in Southwest China decreased from 2004 to 2010, especially in Sichuan Province.
5.3. Energy-water balance and ecosystem response to climate change
5.3.1. Energy-water balance response: droughts analysis based on precipitation
The spatial distribution variability of the index of drought frequency (DF) of 12 months from 1961 to 2012 in SC is shown as Figure 13. Monthly DF has clear change in different months. SC suffered from droughts in a large area in January, February, March, October, November and December. As shown in Figure 13, Yunnan Province and Guangxi Province are severely drought-afflicted areas. A drought pattern also appeared in the monthly variability of DF over time. From January to March, central and eastern Yunnan Province, southwestern Sichuan Province, and southern Guangxi Province experienced droughts (DF >40%) over large area. In May, the drought area rapidly narrowed, and the area of DF >40% was located in a limited extension of Yunnan Province. From June to September, the DF was low with a slightly increase. The distribution of droughts was spread from east to west. In October, droughts began to occur again in a large area, and the DF of eastern Guangxi Province and the DF of northwestern Yunnan Province exceeded 40 and 30%, respectively. From November, the drought area increased rapidly. Half of study area is at a high drought risk. The DF of the entire Yunnan Province in western SC is >40%. In December, the DF of most of the area of Yunnan Province exceeded 50%, except in the limited extension of the central portion, which meaning that these regions suffered droughts biyearly.
To illustrate the spatial distribution of drought variability, the annual DF of SC from 1961 to 2012 was calculated. As shown in Figure 14, the annual droughts DF were scattered, in contrast to the monthly scale. More than 61% of the area in SC had a relatively high DF (>15%), that is, all of SC has been a region of high drought risk for more than half century. The southern and eastern mountain zones have high drought frequency. In contrast, Sichuan Basin, which occupies a large portion of the study area with relatively flat and fertile grounds, suffered fewer droughts events.
5.3.2. Ecosystem response: droughts analysis based on vegetation index
The vegetation index data were 1 km, 16 days composited MODIS EVI (MOD13A12), which were downloaded from the NASA EOS Data Gateway (EDG) (http://modis.gsfc.nasa.gov/index.php). To be consistent with monthly GRACE estimates, the days EVI product was recomposited using the maximum value composite (MVC). The composited monthly EVI data were then linearly interpolated to a 1 × 1 grid.
Droughts are common in Southwest China, and several episodes of potential severe droughts were detected using temporal variability analysis of the GRACE TWS change. To evaluate drought events in Southwest China from 2003 to 2013, three indicators were taken into account: the PAI, the AVI and the annual cycle removed TWS (TWSI). To compute the TWSI, monthly averaged GRACE TWS change data for 10 years were removed at each pixel. Using the spatiotemporal variability analysis described above for the TWS, the precipitation, and the EVI, we chose representative pixels from three regions for drought analysis. Figure 15 shows the time series of the three droughts indicators between 2006 and 2012 when droughts were likely to occur in Southwest China.
The correlation coefficients between EVI and precipitation were all more than 0.84 from 2003 to 2013, whereas they were 0.05, 0.07, 0.08 and 0.16 between PAI and AVI, respectively. The low correlation coefficients between PAI and AVI imply that these two drought indicators predict different water resources deficit conditions that accompany droughts. However, correlation analysis showed that TWSI present low correlation with EVI, precipitation and corresponding droughts indicators. Figure 15 shows that the occurrence and release of drought AVI lagged behind PAI for 1–3 months and droughts of AVI were more severe than PAI. This is because of the delay in recharge of surface and soil water from rainfall and the vegetation growth. Both indicators mainly reflected water depletion in surface water and shallow soil water. PAI and AVI were both invalid for detecting droughts. This was most obvious during the period from the end of 2011 to the beginning of 2012 in Figure 15, when the time series for PAI and AVI had no extreme droughts and a normal and gentle amplitude fluctuation. However, TWSI showed that Southwest China experienced great water resource decreases during this period, which may have caused an extreme drought in 2006. Considering the annual cycle of precipitation and vegetation growth and their relation to shallow water, the TWS change mainly contributed to the discharge of groundwater to surface water, which implied a drought risk.
Southwest China is an ecologically fragile area with more than 242 million populations. The change of water-energy balance and ecosystem due to climate change will affect the regional ecological and human living significantly. Our study of Southwest China in past 30 years shows that ecosystem structures destroyed by shrinking of water ecosystem, forest ecosystem and grass ecosystem, together with the precipitation reduction and temperature rise. In addition, the climate changes also affected the artificial ecosystems such as crops land, which resulting in food production mainly caused by frequently occurred severe droughts. Our research showed that, even in the region with abundant water resources of Southwest China, energy-water balance and ecosystem have severe response to climate change, which is significant to human productions and activities of daily livings.