SOIL WATER CONTENT TEMPORAL-SPATIAL VARIABILITY OF THE SURFACE LAYER OF A LOESS PLATEAU HILLSIDE IN CHINA

Surface soil moisture exhibits an important variability in terms of spatial and temporal domains, which may result in critical uncertainties for agricultural water management. The purposes of this study were (i) to characterize the temporal dynamics and stability of the spatial variability of the surface 0-6 cm soil water content θ on a hill-slope; (ii) to investigate issues related to soil moisture conditions including dominating factors on soil moisture and to the estimation of the mean θ. During a period of more than one month θ was measured on thirteen days by Frequency Domain Reflectometry using a 10 × 10 m grid of measurement points covering a 60 × 280 m domain within a hill-slope of the Loess Plateau in China. Soil water content exhibited a moderate variability for each measurement date, and the correlation length (λ) for θ ranged from 8.4 to 27.7 m. With the soil becoming drier, λ decreased, the CV% and the sampling number for accurate mean θ estimation increased. Aspect, elevation, organic matter content, clay content, and bulk density were the main influencing factors, whose extent of influence weakened with decreasing θ. Based on time stability analysis and on the correlation of mean relative difference of θ with the relative difference of dominating factors, mean θ values were well estimated, with a better accuracy under wetter conditions.


INTRODUCTION
Shallow surface soil moisture is a key status variable in hydrologic processes on the land surface, which tends to be variable in both time and space.Therefore, the knowledge of the characteristics of soil moisture variability is important for understanding and predicting many hydrologic processes (Western et al., 2004).Variability of soil moisture has been studied extensively in the past half century using either classical statistics (Hawley et al., 1982), geo-statistics (Western et al., 2002), or time series analysis (Timm et al., 2006).Relatively less emphasis has however been given to the temporal domain as compared to the spatial scale.
Top layer soil moisture variability can be controlled by a large number of factors, such as soil properties (Hawley et al., 1983), topography (Western & Blöschl, 1999), and vegetation (Reynolds, 1970).It is not straightforward to recognize the different dominating processes or factors influencing soil moisture distributions along different study areas since the dominating factors may also be different under different soil moisture conditions (Grayson et al., 1997;Famiglietti et al., 1998;Leij et al., 2004).
Since the contribution of Vachaud et al. (1985), there has been an increasing interest in temporal stability of soil moisture (Grayson & Western, 1998;Starks et al., 2006).Based on the time stability analysis, Grayson & Western (1998) pointed out that if there are definable features of the terrain and of the soils that could be used to define locations a priori, then the mean soil moisture could also be well determined only from these selected points.
The purpose of present study was, therefore, to investigate the temporal and spatial variability of the soil water content of the 6 cm top layer on a hillslope of the Loess Plateau.Specific objectives were (i) to characterize the temporal dynamics of the magnitude and pattern of the soil moisture variability; (ii) to understand the relative roles of soil and topographic features under different soil moisture conditions; (iii) to obtain the number of samplings needed for the estimation of the mean soil moisture within a chosen significance level based on classical statistics; and (iv) to estimate the mean soil moisture based on time stability analysis and the relationship of soil moisture with its dominating factors.

Field site description
The study was carried out at the Beigeba hillslope of the Liudaogou watershed located in the Shenmu County, Shaanxi Province, China (110º21' to 110º23' E and 38º46' to 38º51' N) (Figure 1), with a maximum elevation of 1,256 m above sea level.The average slope is 19° with a total slope length of 280 m.The soil (cultivated loessial soil), is an Ust-Sandic Entisol of sandy loam texture, and is mainly covered by Purple medic (Medicago sativa L.) and scattered shrub (Caragana korshinskii).The area is characterized by steep gullies and hills, and the soil is loess-derived and easily eroded (Tang et al., 1993).Mean annual precipitation is 437 mm, nearly half of which is received from June to September.The potential evapotranspiration is 785 mm, with a mean aridity index of 1.8 and annual temperature of 8.4ºC, belonging to moderate -temperate and semi-arid zones.

Sampling and measurement
A total of 203 sampling points were setup every 10 m along two perpendicular directions (Figure 1).Since 26 points were located deep in the gullies, 177 remained available for measurements.At each sampling point the Frequency Domain Reflectometry (FDR) technique was used to determine the volumetric soil water content (θ, cm 3 cm -3 ) of the top soil layer (0-6 cm).On 13 days, from May 17 to July 21, 2005, θ was measured along the 280 m of the toposequence slope.In order to reduce the possible uncertainties caused by the timing of samplings, the measurements started about 90 minutes before darkness everyday, and each measurement finished within 70 minutes.Furthermore, measurements with the FDR probe of 6 cm length were taken within a soil surface radius of 25 cm at each sampling point for different dates in order to decrease the influences brought by possible microscale variability.
After θ measurements were completed, soil bulk density was determined gravimetrically and other three soil samples were collected near to the moisture measurement points with a 5 cm diameter hand auger, down to the depth of 5 cm, which were mixed to obtain composite samples.After air-dried, samples were passed through a 1 mm sieve for laboratory analysis.Soil organic matter content was evaluated using the potassium dichromate method.Soil particle size analysis was performed with a MaterSizer2000 laser particle size analyzer manufactured by Malvern, and then the proportion of physical clay (< 0.01 mm) content was calculated according to the Soviet standard (Shao et al., 2006).
Elevation at each sampling location and several additional hillslope points was measured with a differential kinematic GPS, to construct a digital elevation model (DEM).Several terrain-based attributes of the sampling points were also observed with the Mapinfo software, including elevation, slope (tan b, where b is the angle of the maximum slope direction, the N/S line taken as reference with b=0 degrees), cos (aspect) (aspect refers to the direction to which a mountain slope faces, usually expressed as an angle within 0 to 360 degrees), roughness, profile curvature, planform curvature, specific contributing area (a), and wetness index (lna /tan b).

Methods of statistical analysis
The coefficient of variation CV% and the correlation length (λ) were chosen to describe the magnitude and spatial pattern of θ variability.For CV% ≤ 10%, heterogeneity was considered weak, when 10% < CV% < 100%, moderate, and when CV% ≥ 100%, the heterogeneity was taken as strong.According to Nielsen & Wendroth (2003), λ can be obtained from equation (1) and available data of autocorrelation coefficients r(h) for different separation distances h: where r 0 is the autocorrelation coefficient for a separation distance of 0, so that r 0 = 1.When h=λ, r(h) = e -1 .Being concerned only about the spatial pattern in the direction along the slope, the correlation length for this direction was obtained averaging the auto-corre-lation coefficients from all the seven columns of the grid (Figure 1).With the aim of determining which and to which extent factors dominate the soil moisture distribution, Pearson correlation coefficients and crosscorrelation coefficients were calculated.Pearson correlation coefficients can be derived using the SPSS statistical software directly, while cross-correlation coefficients r c (h) for variables A and B are derived from: where, cov and var refer to covariance and variance, respectively, and x is the position coordinate along the slope.
The significance of the cross-correlation coefficient r c (h) is often assessed by the critical values of t, which can be written as follows: where n is the number of pairs used for the calculation of r c (h).
For the purpose of testing the existence of time stability of θ, Spearman's rank correlation coefficient and the relative difference method were used.Generally, the closer the Spearman rank correlation coeffi- cient is to 1, the more stable the data will be.Detailed procedure of calculation for these two methods can be found in Vachaud et al. (1985), Grayson &Western (1998), andStarks et al. (2006).
Under the assumption that the θ data are normally distributed for all the measurement dates, the number n of samples required for their estimation at a given accuracy level k, is given by: where µ a is a critical µ value corresponding to a certain confidence limit p 0 , say, 95% or 90% and k can be set to 5%, 10%, 15%, 20% depending on the demanded accuracy.
Theoretically, within the study area, there may exist sampling locations by which the mean θ of the area can be well estimated.Two methods were applied to determine the sampling locations for estimating mean θ on the hillslope scale.First, these locations can be chosen as those having mean relative differences in relation to the mean close to 0 and also with lowest standard deviation, which can be called as the relative difference (RD) method.Second, these locations can be specified based on the relationship between the mean relative differences of θ and the relative difference of the dominating factors, denoted as dominating factors (DF) method.The second method can be processed as follows: first, the relative difference values of the dominating factors for each location are calculated, then their correlation coefficients with the mean relative difference of θ are also calculated.After this, the products of relative differences of the dominating factors by their correlation coefficients with the mean relative difference of θ at each location are summed.Ranking the sum of products for each location it is possible to select the locations with the sum of products closest to 0, taking them as predictive points.Since many locations can meet the condition of mean relative difference of θ or of the sum of products closer to 0, several locations were considered to investigate the possible effects of sampling number on their predictive ability.

Variability of soil moisture
Soil moisture status at different times -Distributions of θ along the toposequence for each measurement time (Figure 2) varied greatly in terms of quantity during the observational period (minimum and maximum of 2.73% for June 21 and 23.63% for May 29, respectively).For all the data, however, the surface soil moisture presented the same spatial pattern along the hill-slope, i.e. first decreasing and then in-creasing along the toposequence, which is in agreement with Hu et al. (2006) for the same hillslope.Despite rainfall during the study period, the surface soil moisture on each of the 13 days show a similar pattern, which is possibly related to geomorphic differences or some intrinsic differences such as soil texture at the study site.
Dynamics of the magnitude and pattern of the spatial variability -During the observational period, coefficients of variation ranged from 14.9% to 48.2% (Table 1), indicating a moderate spatial variability of θ on the hillslope scale.The relationship between heterogeneity of the variability and the soil moisture condition is shown through the plot of the coefficient of variation versus the mean θ (Figure 3).The coefficient of variation decreased exponentially with the mean θ, which indicates that the 0-6 cm soil moisture on the hillslope tends to be more homogeneous when the soil becomes wetter.The distribution of surface θ is influenced by several factors, thus the change of the magnitude of the spatial variability of θ may reflect and also be the result of the evolution of dominating factors in different soil moisture conditions.Soil moisture variability reflects to some extent the variability of soil porosity, thus soil particle distribution.Likewise, the variation of different particle sizes might explain the θ variability.For the particle classes ranging from 0.05-0.1 mm, the variability was larger for the finer particle class (Table 2).For the θ data of this study, the maximum mean is 23.63%, the soil sorptivity is between 10-20 kPa according to the retention curve of the local soil (Yang & Shao, 2000), and most of water is kept in pores sized within 0.015-0.03mm of equivalent diameter.Generally, small pores tend to be formed by fine particles; therefore, it is not difficult to understand that the CV% increases when the soil dried out .In contrast, beyond the particle class of 0.05-0.1 mm, the variability is likely to be increased.Hence, if θ increases from a low value, as observed in this study, to a high value as saturation, then the variability of θ might decrease first and then increase again.
Correlation lengths λ for all the measurements are also listed in Table 1.The relationship between λ and θ (Figure 4) shows that λ increased linearly with the mean θ (P = 0.033), ranging from 8.4 to 27.7 m.This may imply that the extent over which correlation exists becomes larger as the soil becomes wetter.In other words, the surface θ tends to be much more in-   dependent as the soil dries out.This is in accordance with the investigations of Grayson et al. (1997) who found dry catchment conditions being characterized by a stochastic moisture pattern, while wet catchment conditions showed a much more organized pattern.The evolution of λ along time can be explained by the prevalence of different soil moisture dynamics.Under wetter conditions, variable lateral water flow occurred on the larger area along slope, which can lead to a more organized pattern of soil moisture.As the soil dries out, however, upward evapotranspiration plays an increasingly important role in the soil surface wa- ter movement, which, obviously, results in a greater independence from neighbor areas.

Factors influencing soil moisture distribution
Correlation analysis between soil moisture and topographic and soil properties -Under the assumption that variations in vegetation, meteorological factors and other significant influences on soil moisture variability are negligible at the hillslope scale, focus is given only on soil and topographic properties to identify the dominating factors for soil moisture.Table 3 shows the Pearson's correlation coefficients among topographic and soil characteristics as well as their correlation lengths, which indicate the existence of different extents of correlation between soil and topographic properties and their spatial structure.
Pearson's correlation coefficients for soil and topographic properties in relation to soil moisture observed on different dates (Table 4), indicate that for the same factors their correlation extent with θ can be quite different in terms of their significance level as well as the correlation sign.Generally, when disregarding their significance level, soil organic matter content, clay content, cos (aspect) were positively correlated with θ, while soil bulk density, elevation, slope (tan b), profile curvature, planform curvature, roughness, specifically contributing area, and wetness indexness (lna/tan b) were negatively correlated.When considering significance levels, however, profile curvature and planform curvature were negatively correlated with θ only for observation dates 3 and 5 out of the 13 dates, which agrees with Leij et al. (2004) who virtually found no correlation between retention parameters and curvature of the land surface.For slope and roughness, although negatively correlation existed under all soil moistures of interest, significance can hardly be discovered.For the specific contributing area and wetness index, their negative correlations with θ, although not significant under most of the soil moisture conditions, may not be expected.The wetness index and the specific contributing area are both negatively correlated with elevation (Table 3) and θ under all conditions were also negatively correlated with elevation (Table 4).So, θ should be positively correlated with the wetness index and the specific contributing area, which is also an essentially normal phenomenon in θ estimation.However, the unexpected negative relation obtained in this study may imply in the fact that the wetness index and the specific contributing area may not be good estimators for the soil moisture distribution at the semiarid hillslope scale.
Considering other factors such as aspect, elevation, organic matter content, clay content, and soil bulk density it was found that for most soil water conditions (9 out of 13 at least), their correlations with soil moisture were significant.Among those, organic matter content, clay content, and cos (aspect) were in general positively correlated with θ, while soil bulk density and elevation were negatively linked to θ.Therefore, aspect, elevation, organic matter content, clay content, and soil bulk density are the main influencing factors on the θ distribution.The relative importance of dominating factors in determining soil moisture distribution -The distinct correlation coefficients for aspect, elevation, organic matter content, clay content, and soil bulk density with θ under different observational data can be found in Table 4. Figure 5 characterizes the relationship of their correlation coefficients with mean θ, indicating that their magnitude of correlation increased with increasing values of θ (all passed the significant level test), no matter if negative or positive.Therefore, under higher soil moisture conditions, surface θ was much more controlled by the five main influencing factors.Under lower soil moisture conditions, however, effects of these five main influencing factors were weaker, presenting no or adverse correlations when especially low soil moistures are considered, such as the average θ values of 3.85% on May 15 and of 2.73% on June 21.Therefore, with the drying process of the soil, θ may tend to be much more erratic rather than organized, echoed by the increasingly shorter correlation lengths (Figure 4).
Immediately after a rainfall event, it can be expected that the soil bulk density (through porosity) would dominate the distribution of θ, which can be reflected by the larger negative correlation at higher soil moistures as compared with lower moisture conditions.Meanwhile, because of the high evapotranspi-ration rates in this area, soil evaporation and plant transpiration would be the main active process soon after precipitation, and then organic matter and clay would control θ through their water retention ability.With the soil becoming drier, effects of organic matter and clay on θ became gradually negligible.In fact, since sampling was temporally lagged in relation to precipitation, the correlations of soil bulk density and soil moisture are far weaker than expected because of the dominating process of soil evaporation that occurred soon after that.
Aspect can influence the θ distribution by solar irradiance redistribution.At the beginning soil evaporation tends to be stronger and the influence of the aspect would be largest.As the soil becomes drier, evaporation would decrease until minimum rates, and obviously, the effects of the aspect on θ distribution would also be decreased.This is in agreement with Hébrard et al. (2006) who did not find correlation between the variation of insolation and the evolution of topsoil θ during the second drying sequence.Chanzy & Bruckler (1993) also pointed out that soil insolation was not the main factor controlling the dominant influence when the soil was rather dry.
Elevation can exert large influence on soil moisture by its effects on the possible soil lateral flow along the slope.However, on this research area, possible lat-  eral flow may soon be replaced by the strong evaporation, and therefore, no lateral flow at a dry condition can explain the lower correlation with θ under drier conditions.In addition, the increasingly weaker correlation between elevation and θ may also be the outcome of the combined effects of soil properties.Many researchers have attempted to find out the real driving forces for θ distributions and, for this study, it seems unlikely to have success when considering the complicated correlations among the various topographic and soil properties shown in Table 3. Taking elevation as an example, it is difficult to conclude whether the significant correlation with θ is representative of the effects of its own or other soil properties on θ.However, this may not be so important when choosing surrogates for θ distribution from the statistical point of view.
Cross-correlation analysis between soil moisture and dominating factors -With the aim of understanding the effects of the five dominating factors on the θ distribution with respect to its influencing area, crosscorrelation coefficients were calculated for the five dominating factors with all θ observations.Figure 6 shows only the cross-correlograms of θ for May 17, which did not show symmetry with lag 0 (especially for elevation and aspect versus θ ), implying in a different extent of influence in different directions.The absolute value of the cross-correlation coefficients of θ with all the dominating factors decreased with the lag distance increasing in both directions, which implies that the longer the distance separated by two variables, the weaker the influence imposed by dominating factors.As θ decreased, the cross-correlation coefficients also exhibited a decreasing tendency for a given lag distance as occurred with the Pearson's correlation coefficient (data not shown).
The lag distance coverage within which the cross-correlation coefficients have confidence limits was considered to be the influence distance.After estimating the influence distance with crosscorrelograms, the coefficients were plotted for all the dominating factors against mean θ (Figure 7), and as it can be seen, the influence distance for elevation, soil bulk density, and organic matter content increased obviously with mean θ, while those for clay content and aspect showed a slightly increasing tendency with θ.Therefore, in the process of drying, we can conclude that the factors influencing the decreasing effects on θ can be explained in two ways: (i.) their Pearson's correlation coefficients as well as their cross-correlation coefficients with θ decreased with soil moisture; and (ii.) their influence distance on θ would be shortened as well.Temporal stability analysis of soil moisture Spearman's rank correlation coefficient analysis -The matrix of the rank correlation coefficients of θ for different dates during the entire measurement period is listed in Table 5.In general, the closer the rank correlation coefficient is to 1, the more stable is the process (Vachaud et al., 1985).All the values are highly significant (P = 0.01), indicating that a very strong time stability in the ranks of the locations can be assumed, as did Vachaud et al. (1985) and Pelt & Wierenga (2001).Aspect, elevation, organic matter content, clay content, and bulk density were generally the main influencing factors for the θ distribution for the different measurement dates, although their influences weakened under dry conditions.Since these factors kept relatively invariable for a long time, the strong time stability for θ for different dates becomes very clear.Based on soil water content measurements made with neutron probes, Reichardt et al. (1997) suggested that part of the time stability is due to the use of incorrect calibration curves.Generally, the same calibration curve is used over a whole field and due to the soil matrix variability, each sampling point is affected by a systematic error, or a systematic deviation from the mean.This might well also be true for FDR measurements, and has to be better explored in the future.
Relative difference analysis -After calculating the time average of the relative differences and the temporal standard deviations for each location, the mean relative differences for θ were ranked from the smallest to the largest (Figure 8).Locations R19 and R21, with values close to 0 and with temporal standard deviations below 15%, were considered representative of mean θ at any observational day.Similarly certain locations systematically either overestimate (B2, R22, and B3) or underestimate (L10, E13, and L11) the hillslope average θ, regardless of the observation time.

Estimation for mean soil moisture
In this study, the techniques to obtain accu- rate mean of surface θ is considered in two ways.First, the sampling number was obtained with no exact locations required for a reasonable estimation of mean θ, based on classical statistics.Second, the estimate of mean θ was attempted with a limited sampling number, at fixed locations based on relative dif-  the normal distribution (P > 0.05), and therefore equation (4) can be used to calculate the sampling number required for the estimation of the mean θ under different confidence levels and precision, for all measurement dates listed in Table 6.Two obvious trends existed.The resulting sampling number obtained from equation ( 4) was generally larger for the 95% confidence level than for 90%, and the sampling number increased when the precision increased from 20% to 5%.The sampling number was intimately associated to the CV(%) of θ, since the CV(%) differed with the soil moisture condition, so that it is supposed that the sampling number for the estimation would be different according to the soil moisture condition.Sampling numbers for different soil moisture conditions (Figure 9), under a confidence level of 95% and a precision level of 5%, showed that the sampling number decreased exponentially with increasing θ.Therefore, for the estimation of θ, much more samples have to be taken under dry conditions because of their stronger  heterogeneity.As a consequence, under especially dry conditions, such as the mean soil moistures of 5.04%, 3.85%, and 2.73%, the sampling number under a precision of 5% was even much greater than the sampling that was actually taken (177 points).Considering the large difference of sampling numbers required for mean θ estimation, therefore, caution should be taken in relation to the initial θ values when choosing the input parameters for hydrological modeling.
Estimation of mean soil moisture -Using the RD and the DF methods, surface mean θ was estimated using different numbers of locations at the hillslope scale.Because aspect, elevation, organic matter content, clay content, and soil bulk density are the main influencing factors, these factors were chosen for the DF method.Figure 10 shows the relationship of estimated mean θ versus measured mean θ.Most of the points are close to the 1:1 line, indicating that the two methods can well estimate the mean soil water content.
Ordinarily, only one location whose mean relative difference is closest to 0 and with the least standard deviation is sufficient to estimate the mean soil moisture accurately (Grayson & Western, 1998).Considering a location with no exact zero mean relative difference and the existence of a temporal standard deviation, more than one location have to be used to improve the accuracy of estimation by the RD method.Similarly, this is also true for the DF method.Figure 11 shows the absolute values of relative estimated errors for different locations employed for both methods.The more locations are adopted, the more the absolute values of relative errors decrease, indicating an increasingly enhanced predictive ability.Beyond ten lo- cations however, the predictive ability did not improve much as indicated by the extremely slow change of the absolute values of the relative errors.These relative errors were small beyond ten locations, which were less than 5%.
Based on the knowledge of relatively poorer performance using ten locations in terms of estimation accuracy, only more than ten locations were used to investigate the relative predictive ability under different soil moisture conditions.Figure12 shows their relative predictive ability under different soil moisture conditions, and no matter which method used, their predictive ability tended to be better for wetter conditions, as indicated by the lower absolute value of the relative error for wetter conditions, which agrees well with Leij et al. (2004) who found more accurate predictions for wetter conditions using Artificial Neural Networks.This may be related to the stronger spatial structure in wetter conditions and much more erratic distribution in drier conditions, as analyzed before.On this point, with the soil becoming drier, the decreased accuracy by both methods was in accordance with the increased number of samplings required for the estimation of the mean θ analyzed above.
As for the relative predictive ability of these two methods, the absolute value of the relative error for the DF method were less than that for the RD method, for most number of the locations, except for two or four locations as shown in Figure 11, and the absolute value of the relative error for the DF method was less than that for the RD method under lower soil moisture conditions.They were close to each other although a larger error of DF is shown as compared to RD (Figure 12).Results of samples compared by the T test with respect to absolute value of errors show that the predictive ability of the DF method was slightly better than the RD method, as indicated by the lower, not significant level of the absolute value of errors for the former method as compared to the latter method in Figure 11 (P = 0.216) and Figure 12 (P = 0.372).

Figure 1 -
Figure 1 -Location of the Beigeba hill-slope in the Liudaogou watershed and sampling points.(Sampling marks are shown by numbers and letters).

Figure 2 -Figure 3 -
Figure 2 -Soil water content (θ ) of the 0-6 cm soil layer along the toposequence during the observational period.

Figure 5 -Figure 6 -
Figure5-Relationship between magnitude of the correlation and the mean soil water content (θ) (Abscissas denote the mean 0-6 cm soil water content θ, while ordinates are correlation coefficients).

Figure 7 -
Figure 7 -Influence distances for five dominating factors on soil water content (θ) for different soil moisture conditions.

Figure 8 -
Figure 8 -Ranked mean relative differences for the 0-6 cm soil water content θ (One standard deviation bars and sub-figures referring to minimum, middle, and maximum of the mean relative difference sites numbers are also shown.

Figure 10 -Figure 12 -
Figure 10 -Measured versus predicted mean 0-6 cm soil water content (θ) with (a) dominating factors (DF) method and (b) relative difference (RD) method (Numeral in legend are the number of sites used to predict the mean soil water content).

Table 1 -
Classical statistics of the 0-6 cm soil water content (θ ) data during the observational period.
*included in bracket is adjusted R 2

Table 2 -
Soil particle distribution and coefficients of variation.

Table 3 -
Pearson's correlation coefficients among topographic and soil characteristics and their correlation lengths (λ).

Table 4 -
Pearson's correlation coefficients between 0-6 cm soil water content (θ ) and topographic and soil characteristics during the entire observational period.

Table 5 -
Rank correlation coefficients matrix for the 0-6 cm soil water content (θ ) for different soil moisture conditions (dates).

Table 6 -
Number of samples for accurate measurement of soil water content (θ) for different soil moisture conditions (dates).