Temporal variation of soil moisture over the Wuding River basin assessed with an eco-hydrological model , in-situ observations and remote sensing

The change pattern and trend of soil moisture (SM) in the Wuding River basin, Loess Plateau, China is explored based on the simulated long-term SM data from 1956 to 2004 using an eco-hydrological process-based model, Vegetation Interface Processes model, VIP. In-situ SM observations together with a remotely sensed SM dataset retrieved by the Vienna University of Technology are used to validate the model. In the VIP model, climate-eco-hydrological (CEH) variables such as precipitation, air temperature and runoff observations and also simulated evapotranspiration (ET ), leaf area index (LAI), and vegetation production are used to analyze the soil moisture evolution mechanism. The results show that the model is able to capture seasonal SM variations. The seasonal pattern, multi-year variation, standard deviation and coefficient of variation ( CV ) of SM at the daily, monthly and annual scale are well explained by CEH variables. The annual and inter-annual variability of SM is the lowest compared with that of other CEH variables. The trend analysis shows that SM is in decreasing tendency at α=0.01 level of significance, confirming the Northern Drying phenomenon. This trend can be well explained by the decreasing tendency of precipitation ( α=0.1) and increasing tendency of temperature ( α=0.01). The decreasing tendency of runoff has higher significance level ( α=0.001). Because of SM’s decreasing tendency, soil evaporation ( ES) is also Correspondence to: S. Liu (liusx@igsnrr.ac.cn) decreasing ( α=0.05). The tendency of net radiation ( Rn), evapotranspiration ( ET ), transpiration ( EC), canopy intercept (EI ) is not obvious. Net primary productivity (NPP), of which the significance level is lower than α=0.1, and gross primary productivity (GPP) at α=0.01 are in increasing tenency.


Introduction
The Wuding River basin, in the Loess Plateau, China has been suffered severe soil erosion damages.Understanding long-term change of soil moisture (SM) in the basin, located in a transition zone from farmland and grassland to desert in the Plateau, is very important as it provides useful information for making scientific suggestions for early warning of land desertification, ecosystem recovery and also efficient soil and water management.
Unfortunately, all over the world, including the study basin, long-term SM data are very scarce.Hence, the models and remote sensing technique are used to extend the SM data both in space and time.
Remote sensing data can be effectively used to retrieve SM information at regional and global scale and promote the SM study.However, the length of the available remotely sensed SM data has not been long enough for long-term trend studies yet.
S. Liu et al.: Temporal variation of soil moisture over Wuding River basin Physically process-based models are thus used to obtain reliable prediction of long-term SM.So far several SM datasets have been produced using different land surface process models.However most of the models have simple SM schemes (e.g., Liu et al., 2003) or are forced with monthly average, which made the simulated SM results not agree well with the observations (Chen et al., 1997;Entin et al., 1999;Schlosser et al., 2000).Soil moisture is an important element in hydrological cycle which is closely related to water and energy transfer between soil, vegetation, and atmosphere.It has been shown that much of the global warming so far was at night (Karl et al., 1991(Karl et al., , 1993;;Folland et al., 1992;Stenchikov and Robock, 1995;Robock et al., 2000).So at least models with diurnal cycle are needed to correctly simulate SM.The SM may also be characterized by autocorrelation in time, which means that the lagged effects in inputs or losses can be important as much as those occurring at the time that the impacts are actually observed (Hamlet et al., 2007).A physically process-based model with detailed information about soil-vegetation-atmosphere water and energy transfer is needed to correctly simulate long term SM time series for trend analysis.
There are at least three prevailing methods to do the trend analysis.The first is to draw the linearly regressive trend line over the time series of data to see if it is in upward (increasing) or downward (decreasing) trend (Robock et al., 2000).The second is performing the nonparametric Mann-Kendall (M-K) trend test (Mann, 1945;Kendall, 1975;Hirsch and Slack, 1984;Gilbert, 1987) and judging the trend by the value of the M-K statistic (e.g., Sheffield and Wood, 2008).The last not the least is to compare the SM averaged among each of the decades, which is broadly used in the related studies in China (e.g., Yan et al., 1999).Within the longterm trends, there are noticeable inter-annual and decadal variations in SM for most regions, which weaken the robustness of the trends and have to be considered in trend studies (Sheffield and Wood, 2008).
The M-K trend test has received the greatest attention in hydrology due to its free-distribution nature, which makes it more suitable than parametric alternatives for testing hydrological data.However, the use of M-K trend test is limited to the assumption that the data must be independent and identically distributed (IID).One of the solutions to overcome the problem is transforming the real data into a time series, in which its autocorrelation coefficient is zero.Zero autocorrelation coefficient indicates that the data are independent and identically distributed.So the solution turns to be how to make data series with non-zero autocorrelation coefficient into data series with autocorrelation coefficient being zero.One of such solutions is called "prewhitening" (von Storch, 1995) -von Storch method.It is assumed that data series x i with a non-zero autocorrelation coefficient can be modelled by an autoregressive first order model (Box et al., 1994): where ρ 1 is the lag-1 autocorrelation coefficient (Kulkarni and von Storch, 1995;Douglas et al., 2000;Zhang et al., 2001).However, it is found that the trend computed from such a prewhitened time series is smaller, since the trend in the prewhitened data has a slope (Wang and Swail, 2001) as where µ pw and µ are the trend for the prewhitened and the original data respectively.At least two methods are proposed to avoid this problem.One of them is using the iterative procedure proposed by Zhang et al. (2000), refined by Wang and Swail (2001) and Zhang and Zwiers (2003).In this method, called Zhang method, the autocorrelation is computed after removing the significant trends from the series.This process is continued until the differences in the estimates of the slope and the ρ 1 in two consecutive iterations become negligible.The M-K test for trend is run on the resulted time series and then Sen approach (Sen, 1968) is used to compute the slope.
Another method is Trend Free Pre-whitening procedure proposed by Yue and Pilon (Yue et al., 2002) -called Yue method.In this method, the slope of a trend in sample data is estimated using the approach proposed by Sen (1968).Assuming the trend is linear, the sample data are detrended by subtracting the trend from the sample data.The lag-1 serial correlation coefficient (ρ 1 ) of the detrended series is then computed.If ρ 1 is not significantly different from zero, the sample data are considered to be serially independent and the M-K test is directly applied to the original sample data.Otherwise, it is considered to be serially correlated and prewhitening is used to remove the first order autoregressive process from the detrended series.The residual series after applying this procedure should be an independent series.The identified trend and the residual are combined as the blended series, which just includes a trend and a noise and is no longer influenced by serial correlation.Then the M-K test is applied to the blended series to assess the significance of the trend.
Besides the research shown in the above methods, there have been a long and wide discussion and yet different sounds on the prewhitening (Bayazit andOnoz, 2004, 2007;Hamed, 2007Hamed, , 2008Hamed, , 2009)).Many of such discussions are based on the synthetic data generated by Monte-Carlo method.The prewhitening procedure relies basically on theory.In reality the observations are much more complex than the pure synthetic data.It might not possible to describe the real data with a first order autoregressive model.In this way, discussion of the accuracy of the trend estimation based on the pure synthetic data as described above may be not fully applied.
In this paper, a distributed and physically process-based hydrological model, Vegetation Interface Processes (VIP) model (Mo and Liu, 2001;Mo et al., 2004aMo et al., , b, 2005)), describing the detailed information about water and energy transfer in soil-vegetation-atmosphere system, is used to simulate long term SM data for trend analysis in the Wuding River basin.Observed SM data by both in-situ and remote sensing from 1992 to 2004 are used to validate the model.The validated model is then used to obtain a long term SM series.The M-K scheme is employed to identify the SM trend.The trend results based on original M-K and M-K with three prewhitening techniques to remove the possible influence of auto-correlation on trend analysis are compared.
In the following section the methodology is briefly introduced.The study region and datasets are described in Sect.3. The results of the analysis are given in Sect.4, and the discussions are in Sect. 5. Finally a short conclusion is made in Sect.6.

Model introduction
The VIP model is used to simulate hydrological cycle over the Wuding River basin.The ecological and hydrological processes in the study basin are modelled by implementing a coupled one-dimensional soil-vegetation-atmosphere transfer (SVAT) scheme, a distributed runoff routing scheme (kinematic wave) and a vegetation dynamic scheme.The digital elevation model (DEM) is employed to create the channel flow directions corrected with the channel pattern.Geographical information of vegetation type and land use is incorporated to assign the land surface attributes spatially (Mo and Liu, 2001;Mo et al., 2004aMo et al., , b, 2005)).The estimation of evapotranspiration E T , SM, runoff and vegetation productivity are briefly introduced as follows.
Energy fluxes are described with a two-source scheme discerning canopy and soil surface separately (Shuttleworth and Wallace, 1985).The total canopy E T (kg m −2 s −1 ) composed of soil evaporation E S , canopy transpiration E C and its intercept water evaporation E I is expressed in the form of Penman-Monteith equation (Monteith and Unsworth, 1990).The SM budget is estimated with a six-layer scheme (0-2, 2-10, 10-30, 30-70, 70-130, 130-200 cm).The depth of the root zone is 130 cm.SM simulated by the VIP model represents volumetric soil moisture (cm 3 cm −3 ).The depth of soil layer is 2.0 m.Soil hydraulic parameters are estimated using the Clapp and Hornberger (1978) empirical formula.
The magnitude and timing of overland runoff are affected mainly by rainfall intensity, SM condition and land use/cover.The overland runoff is treated as the ratio of the square of the precipitation, deducted the canopy interception, to the summation of this precipitation and the SM deficit in the root zone (Choudhury and Digirolamo, 1998).The photosynthetic production is input into the vegetation growth module for biomass and leaf area estimation.In the vegetation growth module, vegetation phenological stages are expressed with air temperature degree-day which determines the fractions of assimilation partitioned to vegetation components (leaf, stem, root and grain), and leaf area is estimated by leaf biomass with specific leaf area.Besides SM, such Climate-Eco-Hydrological (CEH) variables as observed precipitation, air temperature and runoff and simulated evapotranspiration E T , leaf area index (LAI), and production by the VIP model are used for SM trend analysis.

Trend analysis method
The M-K test for trend and Sen's slope estimates (by partly using Excel template MAKESENS, developed by Salmi et al., 2002) is used for detecting and estimating trends in the time series of annual mean of SM.The M-K test for a monotonic trend testing in the data and the Sen's method for estimating slope of the trend are introduced as follows, respectively.
The data values x i of the time series can be assumed to obey the model where f (t) is a continuous monotonic increasing or decreasing function of time and the residuals ε i is assumed to be from the same distribution with zero mean.
To test the null hypothesis of no trend, H 0 , i.e. the observations x i are randomly ordered in time, against the alternative hypothesis, H 1 , where there is an increasing or decreasing monotonic trend, a statistic is constructed as follows where x j and x k are the annual values in year j and k, j >k, respectively, and The number of annual values in the data series is denoted by n.
Under the IID assumption (Kendall, 1955), the expectance of the statistic S is zero, and the variance of S is where q is the number of tied groups and t p is the number of data values in the p-th group.By using the statistic S and the above variance VAR(S), a new statistic Z, which is a standardized test statistic, is constructed to be used in the test for n being larger than 10, i.e.,   A positive (negative) value of Z indicates an upward (downward) trend.The standardized M-K statistic Z has a normal distribution with mean of zero and variance of one.In this case, Z 1−α/2 is obtained from the standard normal cumulative distribution tables.As shown in Fig. 1, if the absolute value of Z equals or exceeds the specific value Z α/2 , it means that in this case at the significant level of α, the probability of accepting H 0 is very low, in other words, the existence of a monotonic trend (H 1 ) is very probable.For the significance level of α=0.01 as shown in Fig. 1, it means that there is a 1% probability that the values x i are from the normal distribution and with that probability we make a mistake when rejecting H 0 of no trend.
For an existing trend (as change per year), f (t) in Eq. ( 3) is equal to: where µ is the slope and B is a constant.
To get the slope estimate of µ, the slope values of all data pairs are calculated as: wherej >k .If there are n values of x j in the time series we get as many as N=C 2 n =n(n−1)/2 slope estimates µ i .In Sen's method (cited by Salmi et al., 2002), the slope estimate µ is the median of these N values of µ i which are ranked from the least to the top.The Sen's estimator is In order to give a two-sided confidence interval about the slope estimate, a working variable is defined as Next, M 1 =(N−C α )/2 and M 2 =(N+C α )/2 are computed.
The lower and upper limits of the confidence interval, µ min and µ max are defined as the M 1 th largest and the (M 2 +1)th largest of the N ordered slope estimates of µ i .If M 1 is not a whole number, the lower limit is interpolated.Correspondingly, if M 2 is not a whole number, the upper limit is interpolated.
With the slope estimate µ, it is easy to obtain the estimate of B from Eq. (8) as x i −µt i .The median of these values gives an estimate of B. The estimates for the constant B for the confidence interval are calculated by a similar procedure of µ.
Considering the possible influence of auto-correlation on the M-K method, comparisons are made among the trend results based on original M-K and M-K with three prewhitening techniques, i.e., von Storch, Zhang and Yue method described in the introduction section.In order to compare the trend between SM and the other CEH variables, data standardization (subtracting the mean and dividing by the standard deviation, Liu et al., 2001) is conducted before the trend analysis.
In order to compare and illustrate the effects of prewhitening, the linear trend analysis with significance level determination is also carried out.We propose the following simple method to assess the significance of the linear trend: based on the fact that the correlation coefficient between the time and standardized variable is equal to the linear trend of the standardized data on the standardized time scale, the method to assess the significance of the correlation can be used to assess the significance of the linear trend.

The judgement of the significance of serial autocorrelation
From Yue et al. (2002), to judge if observed sample data are serially correlated, the significance of the lag-1 serial correlation ρ 1 at the significance level of α=0.01, for example, of the two-tailed test is assessed using the following approximation (Anderson, 1942;Yevjevich, 1972;Salas et al., 1980): 3 Materials

Basin description
The Wuding River basin (108 • 18 ∼110 • 45 E, 37 • 14 ∼ 39 • 15 N) with an area of 30 260 km 2 is located in a transition zone from farmland and grassland to desert over the Loess Plateau (Fig. 2).The elevation in the basin ranges from 600 m to 1800 m.The northwestern part of the region is composed of sand desert with a gentle undulation landscape and the southern part is characterized as steep hillslopes with incised channels.The soil texture in the basin  is principally categorized as sandy loam, sandy silt, sand, silt and coarse sand.Since the basin locates in a temperate and semi-arid monsoon climate zone, the annual mean air temperature ranges from 7.9 • C to 11.2 • C and the annual precipitation from 300 mm to 550 mm decreasing gradually from the southeastern to northwestern part.Precipitation events mainly occur in summer monsoon season (June to September) along with occasional heavy storms.The average annual runoff depth over the past thirty years is about 35 mm.

Data
The data can be categorized as: (1) in-situ and remotely sensed SM data; (2) meteorological and hydrological data; and (3) land surface characteristic data.

SM data
In-situ data at two SM measurement sites, Suide (110.21 • E, 37.5 • N) and Yulin (109.7 • E, 38.23 • N), are available in the Wuding River basin.The locations of the sites are indicated by cross symbols in Fig. 2. SM measurements (g g −1 ) in these sites were made with gravimetric method from the top layer down to the 50 cm at the 10 cm interval on the 8th, 18th and 28th each month from 1992 to 2005.The in-situ SM data are used to validate the model.In addition to in-situ SM measurements, a global remotely sensed SM dataset is employed to validate the model at regional scale, which is retrieved from long term scatterometer measurements using a change detection method developed in the Vienna University of Technology (TUW) (Wagner et al., 1999).In the TUW retrieval algorithm, SM time series are extracted from backscatter measurements of the scatterometers onboard the European Remote Sensing satellites ERS-1 and ERS-2.Due to the penetration depth of scatterometer signal decreasing with increasing soil moisture, the TUW data represent the water content in 2-5 cm soil depth in rela-   tive units between the driest and wettest conditions in the period of 1 August 1991 to 31 May 2007 (data acquisition has been done irregularly with on average two measurements per week), resampled into a 12.5×12.5km 2 equally-spaced discrete global grid, which is indicated by solid dots in Fig. 2. In the TUW algorithm, the microwave backscatter measurements normalized at 40 • incidence angle (σ 0 (40)) are used to extract soil moisture dynamics.Eventually the σ 0 (40) measurements are scaled between the lowest and highest values ever observed within the long-term observations representing the driest and wettest conditions.In this way, the σ 0 (40) corresponds to the relative soil moisture values at topmost 2-5 cm soil surface ranging between 0 and 1 (0% and 100%) (Naeimi et al., 2009).

Meteorological and hydrological data from 1956 to 2004
The climatic data at 15 stations in and around the basin from 1956 to 2004 are used to drive the model.The atmospheric forcing variables are daily maximum and minimum air temperature, humidity, wind speed, precipitation and sunshine duration.The daily climatic data are interpolated to the whole basin with the inverse distance square method.Discharge data from the 9 sub-catchments from 1956 to 2004 are used for model validation.

Land surface characterization data
The model characterizes the land surface as topography, vegetation type/density, soil texture and land use.resolution terrain raster.Leaf Area Index (LAI) mainly features the vegetation density, which is estimated by the VIP model.Soil texture data are retrieved from the map at the scale of 1:14 000 000 (Institute of Soil Science, Chinese Academy of Sciences, 1986).Land use data at the scale of 1:100 000 (http://www.resdc.cn/) in 1980s are used for land-use/cover classification, which is divided into six types, namely, farmland, mixed forest, dwarf shrub, grassland and desert with fractions of 29%, 3%, 4%, 43% and 21%, respectively.Farmland is mainly located in the stream valleys, slopes and terraces in the southern part, on which crops, such as maize, millet, soybean, rice and wheat, are planted.Natural vegetation cover in the basin is generally sparse.Excess reclamation and over-grazing have induced vegetation degradation, soil erosion and desertification, which are the main causes of environmental vulnerability in this basin.

Model implementation and analysis
Water and energy balance components are calculated for each grid point separately, neglecting flux exchanges between the grid points.Since the irrigated field covers only a small amount of the farm land area, the farm land is considered as fully rain-fed land.As the energy fluxes' response to atmospheric driving forces is much faster than the response of hydrological processes in the soil, the energy budget module is run on hourly basis and soil water module on daily intervals.
Geographical and vegetation cover data are sampled respectively into the VIP grid resolution of 8 km.The reason of using 8 km sampling is to be in harmony with the resolution of global NDVI data from NOAA.
The data time period is from 1956 to 2004.The model is validated over the period of 1991-2004, by using SM observations in ten-day intervals at two in-situ stations and remotely sensed SM data grid by grid over the basin with at least two measurements per week on average.
Both the modeled and in-situ SM data are scaled with TUW scatterometer data to have uniform datasets by using the following formulations for the validation.
In this way, we extract the relative (normalized) SM (%) from both VIP simulation and in-situ observed SM data.For the comparison at the regional scale, we first record the latitude and longitude of the point of TUW data within the basin.
Then we output the VIP SM data at those 8-km spacing grids with the latitude and longitude matched with those of TUW.These VIP SM data are used to compare with the corresponding TUW data.
By using the validated model, the SM is simulated from 1956 to 2004.One simulation year of 1956 is repeated to avoid of the possible influence of initial status in the beginning of the modelling.The trend analysis was performed using the M-K and Sen's methods after spatially averaging the long-term SM over the whole basin.The seasonal patterns of all CEH variables are compared with each other.Coefficient of variation (C V , Liu et al., 2001) is used to make further analysis of the temporal characteristics of the SM variation.

Validation of the model
Figure 3 shows the simulated SM data at the Suide and Yulin stations compared with the in-situ SM observations and the TUW scatterometer-derived SM data.The Pearson correlation coefficients (R 2 ) between the VIP simulation and the in-situ SM data, between the VIP simulation and the TUW SM data, and between the in-situ and the TUW SM data are 0.28 and 0.19, 0.20 and 0.18 and 0.52 and 0.56 at Suide and Yulin, respectively.The Nash-Sutcliffe coefficients among the in-situ, TUW and the VIP are also low in all cases.Although Nash-Sutcliffe coefficient is often used to evaluate the model efficiency, it may not work well for evaluating of SM simulation because the in-situ observations are sparse and there are many low values in SM dataset.It has been argued  by some authors that the utility of the Nash-Sutcliffe efficiency as a performance measure may be limited by bias in its evaluation (Garrick et al., 1978;Ma et al., 1998;Weglarczyk, 1998;Sauquet and Leblois, 2001), especially because the square summation between the observation and the simulation is not necessarily smaller than the square summation (σ 2 obs ) between the observation and the mean of the observation, when σ 2 obs is small.More details are shown in Mo et al. (2006a).Although R 2 and Nash-Sutcliffe efficiency are not high in all the cases, from Fig. 3, the simulated SM data by VIP model and the TUW data can catch the variation trend and match with the observed SM at Suide and Yulin basically.
In the regional scale, the relative SM (%) averaged over 6 years (1992)(1993)(1994)(1995)(1996)(1997) for both VIP simulation and TUW data are shown in Fig. 4. Generally, at regional scale the averaged simulation SM data comply pretty well with the averaged TUW data except in winter period in which some obvious discrepancies are observable.In winter the simulated SM values are larger than the values of TUW data.It is worthy to mention that this cannot be treated as real large errors as the TUW data in winter should be used in caution.The reason is that the backscatter from frozen soil surface behaves like the backscatter from dry soil.Depending on the wetness of snow cover, backscatter signal changes sporadically under this kind of condition.It is very difficult to predict the behavior of backscatter from a surface covered with snow and frozen soil (Albergel et al., 2009).
Validation is also done based on observed discharge data indirectly.The estimated annual E T from the VIP model is compared with the observed annual E T , which is the difference between the observed annual precipitation and discharge.Figure 5 shows the comparison in one of the validated years from the 9 sub-catchments with R 2 being 0.90.

The seasonal variation and multi-year variability of SM at daily and monthly scales
Before studying the general trend of the SM over the recorded years, it is meaningful to detect the variation of SM at temporal scales including the daily, monthly, annual and decadal scales.
We first look at the seasonal cycle of SM and its multi-year variability (Fig. 6).It is seen that SM follows the summer monsoon pattern in that the maximum SM both for surface and root zone happens in summer, so do SM for all other CEH variables as shown in Figs.7 and 8a, b.During spring, over the basin, SM is relatively stable.There is a sharp decrease of soil water from February to March (in surface) or from January to February (root zone), because the lowest temperature in the basin is not in December, but in January (Fig. 7).This makes the soil surface or the whole of soil profile frozen, with a little lag over the basin, causing the minimum SM occurred in the earlier spring.The small peak of root-zone SM in May is corresponding to warmer weather (higher temperature in Fig. 7), which makes SM higher.This also causes the spring flood in runoff (Q in Fig. 7), but it happens in March, as the surface water (river) gets warm earlier than the soil water.
Besides atmospheric forcing factors, the plant could also influence SM.The SM increases rapidly to 15% in the surface and 12% in the root zone as the growing season proceeds (toward maximum plant growth, see LAI in Fig. 8b).It becomes stable as the vegetation completes its growth cycle.It is easy to see that the peak of surface SM is observable earlier than the peak of LAI and root zone SM, implicating the plant's reliability on rainfall.to the maximum LAI, water demand of plants is high.Due to intensive precipitation during this period, high SM condition is lasting.Since the amount of precipitation is lower in September than in August, SM is also tending to decrease.
The maximum of E T happens in August, which is matched with the maximum of net radiation (R n ) and precipitation (P ) (Fig. 8a, b).The components of E T from canopy (E C ) and intercepted water by vegetation (E I ) reach the maximum values in August and September, which are corresponding to the high LAI in September due to high water consumption by the vegetation.Although E S also reaches the highest in summer due to high temperature and precipitation, its maximum does not happen in September when the leaf covers fully.The maximum of E S reaches maximum in July when there are much water in soil surface for evaporation.The ratio of E C to E T is much higher than the ratio of E S to E T in the growing season, which makes the pattern of total E T being similar to E C at most.The relationship of the two ratios is opposite in winter, where the values for both E T and the components are small.
The patterns of C V can be categorized into two groups.The first is the C V patterns essentially matching with the patterns of the variables, such as root Zone SM, Q, LAI, and E C .The second is the C V patterns being right opposite to the patterns of the variables, such as surface SM, P ,T , R n , E T , E S and E I .For each pattern, the time for the peak (or trough) appearance of the variable is not necessarily matched with the time for that of the C v .
Generally, the multi-year variability of SM over the 48 years  is the least compared with other CEH variables (Table 1).Over the years, runoff at the crosssection of the basin has been used as an important index for trend analysis because of the relative easy data access (e.g., Zhang et al., 1998;Qian et al., 1999;Xu, 2004;Li et al., 2007).However, strong variation of runoff itself may make it hard to get a trend.Compared with runoff, SM, with lower multi-year variability, may be a better and important hydrological component to explore the changes of hydrological elements in the Wuding River basin.This encourages us using SM to describe the CEH trend in the basin.
It is seen that the multi-year change of the averaged SM over the root soil layer is milder than that of the surface SM (Fig. 6).Over the years, the largest variability of SM in the surface happens at the spring peak and lowest variability at the summer peak.This is matched with the lowest variability of the precipitation peak in summer (Fig. 7).Oppositely, the largest variability of SM in the root zone is observed at the summer peak and the lowest variability at the spring peak, illustrating that the surface SM is closely related to precipitation and net radiation, and SM in root zone is more related  to E T , being influenced more comprehensively by both the atmospheric forcing and vegetation dynamics (Fig. 8).
No matter how high or low the variability of SM in summer is, the seasonal variation of C V for surface and root zone SM is obvious, a feature in rainfed situation as in our case.For irrigation case instead, the C V 's variation is much stable (Mahmood and Hubbard, 2004).Under variable hydroclimatic condition (precipitation and temperature), plant demand of moisture and insufficient water supply result in the alternative peak and minimum pattern of SM variability.
By comparing the seasonal cycle and its variability with the monthly (Figs.6-8 ) and daily SM , it is seen that the seasonal cycles are similar to both the SM and other CEH variables.As expected, changes in SM from one to the next month are not as drastic as those at the daily scale.Overall, the monthly variability is relatively "smooth" compared to that at daily scale, with the means of C V at the monthly scale being less than the means of C V at the daily scale for all the variables (Table 1, and Fig. 12).More fluctuating local hydro-climatic conditions at daily scale, namely, wet or dry, are averaged monthly, reducing SM variability at monthly scales.As some variability has been averaged, some pat-terns are clearer at monthly than at daily scale.For example the double peak of runoff and SM is shown clearly at the monthly scale, which is hid at the daily scale.
From Fig. 12, it is shown that the values of C V of SM at daily, monthly and yearly scale are all smaller than those of other CEH variables.Lower variability in SM can be explained at least in the following two courses.Firstly, the Wuding River basin is located in a semi-arid climate zone where annual precipitation is about 400 mm yr −1 and most of it falls in summer.The low precipitation leads to low surface SM, which is close to the wilting point.Secondly, scientists found that owing to its larger inertia and longer memory to atmospheric driving force than other hydrological components, SM, along with snow cover, is the most important component of meteorological memory for the climate system over the land (Delworth andManabe, 1988, 1993;Robock et al., 2000).All of these support our first impression to the low variability of SM in the Wuding River basin.
It is worthy to point out that the value of the variable, its standard deviation and its C V should be used conjunctively for a more comprehensive analysis.Only using standard deviation or C V for variability analysis is not enough (Liu et    2001; Mahmood and Hubbard, 2004).For example, the standard deviation of R n in summer is larger than that in winter season.However, the C V in summer is the lowest.Usually, the relationship between the standard deviation and the value is positive and the relationship between C V and the variable itself (and the standard deviation) is negative (Fig. 13).There is a minor shortcoming in using C V for variability analysis.For example, when the value is turning from negative to positive or vice versa, such as in the R n and T cases (Figs. 7 and 8,Figs. 10,11), the value of C V is much higher than what is expected.When the value is at the turning point from high to low and when the mean values are very low, such as in the LAI and E C cases, the C V values are also very high.However, this is not really meaning that the multi-year variability at this time is much larger than that in other periods.

The interannual variation of SM and its variability within each year
The interannual variation of SM and its within-year variability are shown in Fig. 14.There is a decreasing tendency for SM both in the surface and the root zone.Year 1999 is the driest among the simulation years.The within-year variability keeps the same over the years, also seen in Fig. 15.This pattern follows that of precipitation.So does for runoff.We can see that although the values of SM, P and runoff (Q) are in decreasing tendency, their within-year variability keeps relatively stable.Comparing with P , Q and all other CEH variables (Figs.16 and 17 and Table 1), the within-year variability of SM is the least.This again encourages us using SM to explore the change signal of CEH processes in the basin.
It is clear that over the years, the basin is getting warmer, but the within-year variation of temperature is getting smaller, showing that the seasonal cycle of temperature is weakened over the years.Over the years, R n , E T and its components remain stable, which shows the compensation effect of precipitation and temperature.
Table 2 shows that the SM in 1950s, 1960s and 1970s is higher than the multi-year average.Since then, it keeps drying in 1980s, 1990s and 2000s.Averagely, root zone SM decreases at the rate of 0.01 cm 3 cm −3 per decade, this is in agreement with the trend analysis by Nie et al. (2008) who found the SM at the top 50 cm in this region decreases at the Hydrol.Earth Syst.Sci., 13,2009 www.hydrol-earth-syst-sci.net/13/1375/2009/

Fig. 12.
The mean of the C V at the daily, monthly and annual scale (q1 and q2 represents θ 1 and θ 2 , respectively).For R n , the four large values (41.47, 381.43,131.1, 456.92) are not included.rate of 0.00072 kg kg −1 based on the observed data from 13 stations from 1981 to 1998.One interesting phenomenon is that the decreasing tendency of precipitation has been shown in 1970s already, so does runoff with an immediate correspondence.However, the decreasing tendency of SM has a delay, which appears only in late of 1980s.This shows that there is a lag in SM's response to precipitation.

The total variation trend of SM over the 48 years
The M-K method is used to explore how significant of the trend is.Although there is a suggestion to use the SM in summer (June, July, August) to test the trends of SM as summer drying will accompany global warming (Robock et al., 2005), the annual mean SM is used in the trend analysis in a more comprehensive way.
Some studies made the variables standardized/normalized at first before they did the trend analysis (Zhao and Yan, 2006).Scaling the axes, or standardizing the variables is a common technique (e.g., Liu et al., 2001;Ouarda et al., 2001;Mo and Beven, 2004;Riad et al., 2004).Application of such  Table 3.The statistics of M-K trend and Sen's slope estimate for the 11 CEH variables by using original M-K, M-K with prewhitening by von Stroch, Zhang and Yue method respectively.(Z and µ are explained in Sect.2.2.ρ 1 is the lag-1 auto-correlation coefficient).4. The final estimate of the trend and its significance for each of all the 11 CEH variables (the significance of linear trend and the lag-1 autocorrelation coefficient ρ 1 , as shown in columns 4 and 6, are assessed by the method as described in Yue et al. (2002).The ratio of column 2 and column 3 is the standard deviation of the years from 1956 to 2004, i.e., 14.The meaning of symbol of the significance is the same as in Table 3).

Linear trend Lag-1 cor
The final estimate of trend Var.
Original time scale Standardized time scale Sig.transformations prior to the analysis attempts to remove the influence of scaling effects from the analysis, with the aim to obtain a simpler expression (Ouarda et al., 2001).Furthermore, when plotting the results for a large number of catchments, variables with the largest means tend to dominate the display.As discussed by Friendly and Kwan (2003), scaling axes can obviously lead to an incoherent display in which no systematic trends or relations can be distinguished.There-fore although as shown in Liu et al. (2008), scaling may bring the uncertainty of the results, we still use the normalized SM for trend analysis, in order to make a better comparison with other CEH variables.
In order to remove the effect of auto-correlation, the data are prewhitened by three techniques proposed by von Storch, Zhang and Yue, respectively, before making the trend analysis with M-K method.In order to compare the results with Hydrol.Earth Syst.Sci., 13,2009 www.hydrol-earth-syst-sci.net/13/1375/2009/  and without prewhitening, a linear trend and its significance are also calculated based on standardized dataset.The results are shown in Table 3.
Generally, in most of the cases, no matter which prewhitening techniques are used, they all lead to potentially inaccurate assessments of the significance of a trend.In most cases for P , R n , E T ,E C , E S , and T , they follow the same pattern as those deducted from numeric experiments (Yue et al., 2002;Zhang et al., 2004).That is, removal of positive lag-1 auto-correlation from the data by prewhitening will remove a portion of trend and hence reduce the possibility of rejecting the null hypothesis while it might be false.Contrarily, removal of negative lag-1 auto-correlation by prewhitening will inflate trend and lead to an increase in the possibility of rejecting the null hypothesis while it might be true (Yue et al., 2003).For those not following the above pattern, it may be due to either the very small trend as in the E I and NPP cases, or the large trend with small ρ 1 as in the cases of GPP and Runoff.
However, for all techniques, as shown in Fig. 18 and Table 3, the differences of the estimates of the trends are not very much.The final estimate of the trend is depended on the significance of the auto-correlation coefficient.The significance is assessed using the same method as in Yue et al. (2003).When the data are not autocorrelated significantly, the trend from original M-K and Sen method is used.When the data are autocorrelated significantly, we choose the trend from the three techniques, which most approaches the original trend and with lowest ρ 1 .The final estimate of the trend and its significance for each of all the 11 CEH variables are shown in Table 4.
Although the least square estimator of linear trend is vulnerable to gross errors and the associated confidence interval is sensitive to non-normality of the parent distribution (Sen, 1968, from Wang andSwail, 2001), it is still used to make a comparison with the final estimate of trend.It is found that after normalizing both the variable and time (year), the correlation between the variable and the time is right equal to the linear trend of the normalized variable at the normalized time scale (column 3 in Table 4).It is thus possible to use the same method of assessment of the significance of correlation to assess the significance of the linear trend, as shown in column 4 in Table 4.
From Table 4, it is shown that the Wuding River basin is in drying tendency at α=0.01 level of significance.Runoff, precipitation and soil evaporation (E S ) are also in a decreasing tendency.The significance of SM is lower than that of runoff (α=0.001),but higher than that of precipitation (α=0.1) and soil evaporation (α=0.05).Temperature is in increasing tendency at the same level of significance as SM.The tendency of net radiation (R n ), Evapotranspiration (E T ), transpiration (E C ), canopy intercept (E I ) is not obvious in that the M-K test indicates a decreasing trend for R n , E T and E C and increasing trend for E I , but all at the less significance level than α=0.1.As the trend for E I is very small, prewhitening techniques are possible to change the direction of the trend.Thus for a small trend, discussion on its direction has no realistic meaning.
It is interesting to see from Table 4 that although soil is in drying tendency, vegetation productivity (NPP (net primary productivity), especially GPP (gross primary productivity)) is in increasing tendency.It is seen that the temperature is in increasing tendency.In North China, it is found that if no CO 2 fertilization effect is considered, crop yield is reduced with an increase in temperature (Liu et al., 2009).However, if considering CO 2 fertilization, crop yield may increase.As our model not only deals with hydrological process, but also with vegetation dynamics, the response of productivity may be affected by the precipitation, CO 2 fertilization and temperature.It is thus understandable why a drier soil may produce higher vegetation productivity.
It is worthy to mention that, although the absolute value of SM trend is small, the normalized trend of SM is large (Fig. 18).This indicates that northern drying is an important signal in the basin.

Northern drying?
The drying tendency in the north of China, Northern Drying, has been paid strong attention to.Tao et al. (2003) produced a long-term SM dataset over China, represented by soil deficit index, from 1946 to 1995 by using a conceptual model.They found there was a trend toward soil drying in the North China Plain and the North-eastern China area during this period.By using the observed SM data from 35 stations over North China from 1981 to 2002, Zhao and Yan (2006) found that the whole trends of annual mean SM storage averaged by different soil layers from soil surface to top 50 cm was decreasing generally.By using SM dataset of the top 50-cm soil layers at 178 SM stations in China from 1981 to 1998, Nie et al. (2008) found that there were increasing trends of SM for the top 10 cm, but decreasing trends for the top 50 cm of soil layers in most regions.Our study further identifies this Northern Drying phenomenon.
The same drying trend is also found in West Africa, as driven by decreasing Sahel precipitation.For most of other areas in the world and in some parts of China, the trends are different.For example in Ukraine, Russia, Mongolia (Robock, 2000;Robock and Li, 2006) and most parts of the US (Robock, 2000;Andreadis and Lettenmaier, 2006;Sheffield and Wood, 2008), soil is found in wetter tendency.In the far North China Qaidam Basin of Qinghai Province in northwestern China, Yin et al. (2008) reconstructed the SM conditions from tree ring from year 566 to 2001 and revealed a general trend toward a wetter condition during the most recent 300 years.Tao et al. (2003) reported that there was a significant increase in SM levels in Southwest China, and a generally insignificant increase or decrease trend in SM levels in Southeast China.Europe, southeast and southern Asia appears to having not experienced significant changes in SM (Sheffield and Wood, 2008).

Man-made Northern drying or nature-made
The definitions of man-made and nature-made change in different disciplines are not the same.For example, from climatology's view, man-made change includes a) local changes such as land use, b) large-scale changes that caused by anthropogenic forcing such as global warming.And the nature made changes include a) natural internal variation of the climate system and b) natural forcing external to earth climate system such as changes in solar output.In this paper, we only pay attention to the regional-local land use change as man-made change and natural internal variation of the climate system as the nature-made change.
The catchment of the Wuding River, in the Loess Plateau of China, is one of the well-known areas suffered most serious soil and water loss and thus becomes one of the main tributaries to produce sediments for the Yellow River.In order Hydrol.Earth Syst.Sci., 13,2009 www.hydrol-earth-syst-sci.net/13/1375/2009/  to conserve soil and water, it has taken many soil water conservation countermeasures, including planting trees, building terraces and constructing checking dams.These countermeasures not only mitigate the sediment yield from the basin but also change the hydrological processes in the basin via land use and cover change (LUCC).Therefore any changing tendency from the observed data may be contributed by human activities or climate change.It is a big challenge to distinct their contributions.There are some methods such as "baseline" method (Lorup et al., 1998;Yang et al., 2005;Li et al., 2007), "regression" (e.g., Xu, 2004) and "optimal fingerprint" method (Hasswlmann, 1997;Allen and Scott, 2003;Zhang et al., 2007) to deal with this problem.In our research the long-term SM data are simulated by actual at-mospheric forcing and land cover, the trend shown by the data should reflect more the influence of nature factors.The decrease in precipitation seems the primary reason for the trend of SM.Logically, with the increase in temperature in Wuding River, E T may increase consequently, and will further reduce the soil water storage.However in our case the total E T and canopy transpiration shows a declining trend at an insignificant level.Soil evaporation shows a downward trend at the significant level of α=0.05.Actually, E T is not linearly related with temperature as described by VIP, which plays a complicate role in temperature-E T -SM relationship.

The representative of SM simulation via VIP
Although the simulated SM from VIP mainly reflects the changes of climate, it is still a good indicator of the wetness status of the basin in the real world.There are at least three reasons.Firstly, from our previous study and documented studies, man-made change is a main, at least half of the contribution to the change of the runoff which is showing the decreasing tendency (Mo et al., 2006;Li et al., 2007).However, for soil moisture, as it has the least coefficient of variation relative to other CEH variables, it has been used as a more effective hydrological indicator of climate change (Robock et al., 2000).Secondly, as VIP model is a distributed model, it can be used to detect the soil moisture under different land use scenarios, which has been a worldwide research topic.From the documented reports, it is found that there are discernible SM differences between land use or crop types planted (Cai and Wang, 2006;Yang and Rong, 2007;Yang et al, 2008).This is one of the important reasons to use VIP to simulate SM for the Wuding River.Last but not the least, from the observed records, it is shown that the crop types haven't changed much during the last twenty years in the Wuding River basin.The land use change over that period is less than 5%.This is the reason that although we do not consider very much the temporal variation of vegetation type for SM simulation, the SM simulation still matches the observed in an acceptable way.
Based on the above analysis, owing to SM's slow temporal and yet obvious spatial variation over the basin under different Land use/covers, the simulated SM time series from Hydrol.Earth Syst.Sci., 13,2009 www.hydrol-earth-syst-sci.net/13/1375/2009/  the VIP model, which reflect the temporal changes attributed more to climate and yet also reflect the spatial variation over the basin, can be assumed to represent the realistic conditions.The agreements between the simulated and observed SM at the stations and the regional remotely sensed SM support this assumption.

Conclusions
The Northern Drying phenomenon is found in the Wuding River basin, China, one of areas suffering most serious water and soil loss.As the basin just has relatively short-term observed soil moisture (SM) data at two stations, an ecohydrological processes-based model (VIP model) is used to simulate the long-term daily SM data from 1956 to 2004 to explore this trend.
The model is at first validated by both the in-situ SM data at two stations and the remotely sensed SM data produced by the Vienna University of Technology (TUW).The averaged TUW SM data over 189 points matched with the simulated SM by VIP well except significant differences in winter.Although there are differences between the simulation and measured/remotely sensed SM data, their variation trends are matched each other.From the first aspect, this encourages us to use the model to simulate the long term SM data from 1956 to 2004 for the trend analysis.
With the long term SM data, the seasonal cycle, summer monsoon pattern, and multi-year variation are analyzed based on the values of SM themselves, their standard deviation and the coefficients of variation at the daily, monthly and annual scale.The temporal variation patterns of SM are well explained by the other CEH variables such as precipitation, temperature, net radiation, E T and LAI.This encourages us to use the data for the trend analysis of SM from the second aspect.
It is found that the inter-annual and seasonal variability of SM is the lowest compared with that of other CEH variables.This, from the third aspect, encourages us using SM data for the trend analysis.For a variable, which has a slight change over years, to identify its trend is more meaningful than to identify the trend for a much varying variable.
The trend analysis is then done by the M-K based on the long-term simulated SM data with the consideration to autocorrelation.It is shown that SM is in decreasing tendency at α=0.01 level of significance, confirming the Northern Drying phenomenon as shown in the documented studies done elsewhere.Its significance is lower than the significance of the decreasing tendency of runoff and increasing tendency of temperature (α=0.001),whereas higher than the decreasing tendency of precipitation (α=0.1).This indicates that the change of SM is not only due to the climate (precipitation) but also other natural factors.

*
From the standard normal cumulative distribution 1.The normal distribution for the two-tailed hypothesis test in Mann-Kendall method.

Fig. 1 .
Fig. 1.The normal distribution for the two-tailed hypothesis test in M-K method.

Fig. 2 .Fig. 2 .
Fig. 2. The SM data distribution (Red dots: TUW grid point (12.5×12.5 km) of the ERS scatterometer-derived SM data; 189 blue square: VIP grid point (8x8 km, matched with TUW grid) of SM (Relative water content of soil surface (2.5 cm, wc/wc sat ) simulated by VIP model Black star: Observed point of SM (relative water content of soil surface (10 cm, wc/wcsat) observed at Suide and Yulin)

Fig. 3 .
Fig. 3.The yearly evolution (upper two panels) of and the correlation (lower three panels) among relative SM simulated by VIP, retrieved from remote sensing TUW data and observation at Suide and Yulin.

S.Fig. 4 .
Fig. 4. The averaged (a) seasonal variation and (b) yearly evolution of relative SM over the 189 grid points simulated by VIP and retrieved from remote sensing by TUW.

Fig. 4 .
Fig. 4. The averaged seasonal variation of relative SM over the 189 grid points simulated by VIP at 10 cm soil layer and retrieved from remote sensing by TUW.

Fig. 5 .
Fig. 5.The simulated averaged annual evapotranspiration (ET) and the observed ET which is the difference between the observed annual precipitation and discharge based on water balance (wb) of the 9 sub-catchments of the Wuding River Basin in 2000.

Fig. 5 .
Fig. 5.The simulated averaged annual and the observed E T which is the difference between the observed annual precipitation and discharge based on water balance (wb) of the 9 sub-catchments of the Wuding River basin in 2000.

Fig. 6 .
Fig.6.Seasonal cycle (left panel) of monthly mean SM in the 0-2 cm surface layer (θ 1 ) and root zone (θ 2 ) averaged over the entire data records and their variability (right panel).Months are indicated with the numbers on the abscissa.Shown in the figure is also the one standard deviation.

Fig. 7 .
Fig. 7. Average seasonal cycle of precipitation (P), air temperature (T) and runoff (Q) at the

Fig. 7 .
Fig. 7.Seasonal cycle (left panel) of precipitation (P ), air temperature (T ) and runoff (Q) at the control station of the basin averaged over the entire data records and their variability (right panel).Months are indicated with the numbers on the abscissa.Shown in the figure is also the one standard deviation.

Fig. 8a .
Fig. 8a.Seasonal cycle (left panel) of net radiation (R n ), total evapotranspiration (E T ), canopy transpiration (E C ), averaged over the entire data records (1957∼2004) and their variability (right panel).Months are indicated with the numbers on the abscissa.Shown in the figure is also the one standard deviation.

Fig. 9 .
Fig. 9. Mean daily SM (Left panel) and SM variability (Right panel) at daily time scale in the

Fig. 11a .
Fig. 11a.Mean daily values (left panel) and their variability (right panel) at daily time scale for net radiation (R n ), total evapotranspiration (E T ) and canopy transpiration (E C ). Day 1=1 January.

Fig. 11b .
Fig. 11b.Same as Fig. 11a but for soil evaporation (E S ), evaporation from the intercepted precipitation by the canopy (E I ) and leaf area index (LAI).

Fig. 13 .
Fig. 13.The comparison among the value, standard deviation and the coefficient o

Fig. 13 .
Fig. 13.The comparison among the mean value, standard deviation and the coefficient of variation for the 11 variables.

Fig. 14 .
Fig. 14.The annual mean of SM in the 0-2 cm surface layer (θ 1 ) and root zone (θ 2 ) (left panel) and their variability (right panel) over the years.

Fig. 16 .
Fig. 16.The annual mean of precipitation (P), air temperature (T) and runoff (Q) at the

Fig. 16 .
Fig. 16.The annual mean of precipitation (P ), air temperature (T ) and runoff (Q) at the control station of the basin (left panel) and their variability (right panel) over the years.

Fig. 17a .
Fig. 17a.The annual mean of net radiation (R n ), total evapotranspiration (E T ), canopy transpiration (E C ) (left panel) and their variability (right panel) over the years.

Fig. 17b .
Fig.17b.The same as Fig.17abut for soil evaporation (E S ), evaporation form the intercepted precipitation by the canopy (E I ) and leaf area index (LAI).

Table 1 .
The mean, range (=max-min) and the maximum of C v for the 11 CEH variables at the daily, monthly and annual scale. al.,

Table 2 .
The decade mean of SM (the first line) and the ratio of the difference between the decade mean and the multi-year average to the multi-year average (the second line) for precipitation (P ), net radiation (R n ), evapotranspiration (E T ) and its components, GPP, NPP, runoff (Q), root zone SM (θ 2 ) and air temperature (T ) of the Wuding River basin.
Fig.12.The mean of the CV at the daily, monthly and annual scale ( q1 and q2 represents θ and θ2 respectively).For Rn, the four large values(41.47,381.43,131.1,456.92)are not included.
if trend at α=0.001 level of significance; b if trend at α=0.01 level of significance; c if trend at α=0.05 level of significance; d if trend at α=0.1 level of significance a if trend/correlation at α=0.001 level of significance, for correlation coefficient,lower bound=-0.49,upper bound=0.45b if trend/correlation at α=0.01 level of significance, for correlation coefficient,lower bound=-0.39,upper bound=0.35c if trend/correlation at α=0.05 level of significance, for correlation coefficient,lower bound=-0.30,upper bound=0.26d if trend/correlation at α=0.1 level of significance, for correlation coefficient,lower bound=-0.26,upper bound=0.22 a In this semi-arid basin, annual E T is mainly regulated by precipitation.S. Liu et al.: Temporal variation of soil moisture over Wuding River basin