Spatial Pattern and Temporal Stability of Root-Zone Soil Moisture during Growing Season on a Larch Plantation Hillslope in Northwest China

Soil moisture plays a decisive role for tree growth and forest ecosystems services supply in dryland regions. Hence, it is necessary to clarify the spatio-temporal variation of soil moisture under field conditions. This study selected a hillslope in the Liupan Mountains covered by the plantation of Larix principis-rupprechtii Mayr (larch), a main afforestation tree species in north and northwest China. The volumetric soil moisture (VSM) in root zone layers was monitored with a time interval of about 15 days during the growing season (from May to October) in 2016 at 48 points on this hillslope. The aim was to evaluate the spatial pattern and temporal stability of soil moisture at slope scale. The results showed a moderate spatial variability of VSM in each soil layer, with the variation coefficients range of 17.12–22.63%. The spatial variability of VSM showed a dependence on the soil wetness and a threshold effect, it increased with rising VSM until the VSM reached a threshold of about 15%, but thereafter decreased. The mean relative difference (MRD) among the 48 points ranged from −30.56% to 27.20%, −29.89% to 39.58%, and −28.13% to 33.71% for the soil layers of 0–20, 20–40, and 40–60 cm, respectively. The associated standard deviation (SDRD) (and range) was 11.38% (5.20–26.06%), 8.28% (4.64–15.63%), and 6.51% (2.00–14.16%) for the soil layers of 0–20, 20–40, and 40–60 cm, respectively. The high Spearman’s rank coefficients (p < 0.05) among the measuring dates at each soil layer indicated that the spatial distribution of VSM in the root zone had strong temporal stability. The decrease of Spearman’ rank correlation coefficient and mean SDRD with rising soil depth indicated an increasing temporal stability of VSM with rising soil depth. The mean VSM of the three soil layers on the entire hillslope can be estimated by the direct method (using representative points determined by the index of temporal stability (ITS)) successfully, and these representative points determined by ITS were mainly located at the points with a ratio of field capacity to leaf area index (LAI) close to the slope mean. Moreover, the mean VSM of the three soil layers on the entire hillslope can also be estimated by indirect method (using the time-stable points determined by mean absolute bias error (MABE) and considering the offset between slope mean VSM and observed VSM at time-stable points), and the prediction accuracy of the indirect method was better than the direct method. Significant correlation between MRD and soil bulk density, field capacity, capillary porosity, and LAI were observed for all soil layers, indicating that both the water-retention ability in root zone soil (expressed mainly by field capacity) and water-consumption ability of trees (expressed mainly by canopy LAI) are the main factors controlling the spatial pattern of root-zone VSM on the larch plantation hillslope studied.


Introduction
Soil moisture is a key factor affecting a series of eco-and hydrological processes in dryland regions, because it controls the rainfall infiltration, evapotranspiration, runoff generation, vegetation growth, and photosynthesis [1][2][3].Meanwhile, it is also an important input for hydrologic, tree growth, and climate models [4][5][6][7].However, soil moisture is highly affected by many factors (e.g., land use, soil properties, terrain, meteorology, and vegetation dynamic), and presents great spatio-temporal variation and scale-dependence [8][9][10].This will certainly affect the associated ecoand hydrological processes and related model applications.Therefore, determining the soil moisture and its spatio-temporal patterns becomes particularly important.Additionally, owing to the high spatio-temporal variability of soil moisture, a large number of sample points are required to derive sufficient information on soil moisture, which is costly in terms of both time and money.An approach that can optimize the sampling points to precisely assess the mean soil moisture is thus desirable [11].
Although there are several methods to investigate the space-time dynamics of soil moisture, such as the empirical orthogonal function (EOF) methodology [12,13] and modeling [14], the temporal stability analysis, which can judge whether the spatial pattern of soil moisture is persistent and help to determine the mean soil moisture as the average of soil moisture at some sites with lowered observation points [15][16][17], has attracted much attention and has been widely used for analyzing the temporal stability of soil moisture patterns [15], providing missing data [2], and making spatio-temporal validations [18].The temporal stability concept has been proved to be an effective way to analyze the spatial pattern and temporal stability of soil moisture [17,[19][20][21][22].In past decades, many such studies have been carried out under different conditions of land uses (e.g., grassland, orchard, shrub, farmland, forestland, and gravel-mulched field), climate (e.g., semi-arid, semi-humid, and humid), and spatial scales (e.g., plot, hillslope, and watershed).For example, Pan et al. [23] reported a stronger temporal stability of soil moisture under humid conditions than drought conditions in an artificially re-vegetated desert area.Gao et al. [24] found higher temporal stability of soil moisture at the points with higher clay content for a jujube orchard on a slope.Zhao et al. [21] used the temporal stability analysis to determine a representative location for estimating the mean soil moisture of a gravel-mulched field.These studies provided a theoretical basis for soil moisture prediction and management under different conditions.However, studies on the temporal stability of soil moisture patterns on forested hillslopes (especially completely forested hillslopes) remain rare.The soil water redistribution and the high heterogeneous canopy structure on forested hillslope will strengthen the variability of soil moisture, and meanwhile, the absorptive capacity of the roots of forest/vegetation vary with soil depth, which might affect the temporal stability of soil moisture across different soil layers.Therefore, how do the spatial patterns and temporal stability of soil moisture vary with soil depth on completely forested hillslopes?Can we use several representative points to accurately estimate the mean soil moisture of forested hillslopes?Moreover, do the representative points have a specific distribution rule?All these questions still need to be further investigated.
The representative points are typically defined as the points where the measured soil moisture values either are close to the mean soil moisture or can be easily converted to obtain such mean values [15].The representative points tend to have good temporal stability features.Several methods have been proposed to determine the representative points by temporal stability analysis.The simplest method is to use the point with the mean relative difference (MRD) closest to zero [16].However, this method ignores the fact that MRD is actually a form of statistics with inherent errors (standard deviation, SDRD) [15].This means that the point with MRD close to zero cannot guarantee its high temporal stability.Therefore, a modification of the above method is to use the points with MRD close to zero and with smaller SDRD values (<5%) to represent the mean soil moisture [25].The soil moisture of the representative points determined by this method can accurately represent the mean soil moisture.However, there might be no such points satisfying both requirements in some areas with large spatio-temporal variation of soil moisture [6].An alternative method is to use the points with the smallest index of temporal stability (ITS) (a combination of MRD and SDRD) [17].
Forests 2018, 9, 68 3 of 22 Hu et al. [20] also proposed an indirect method to estimate the mean soil moisture by introducing a constant offset (MRD) on the point with minimum value of the mean absolute bias error (MABE).However, each method has its own advantages; no individual method can always perform better than another [10], but the performance of different methods varies with land use, climate, and spatial scale.For forested hillslopes with randomly distributed soil moisture, the criteria of several methods might not be satisfied, and therefore the selection of representative points should synthetically consider the performance of different methods.
Several studies found that the representative points tended to equal or were close to the average characteristics of the study area (e.g., the average soil particle content, the mid-slope location) [26,27].Several studies reported that the most time-stable points were related to the soil properties [17,28], and several studies also noted that the time-stable points were poorly related to topography and soil properties [29].These various findings indicate that the controlling factors of the time-stable points (the representative points) are very complicated and dependent on the characteristics of the study area (e.g., climate, soil, and vegetation).The variability of soil moisture on forested hillslopes might be much larger than in the aforementioned studies because of the soil water redistribution and the heterogeneous canopy structures.However, knowledge about the controlling factors of the representative points on forested hillslopes has been limited.
Larix principis-rupprechtii (larch) is one of the major afforestation tree species in the mountain areas of dryland regions in north and northwest China.For example, its plantation area amounts to 31.3% in the Liupan Mountains (LPM) area in northwest China, thus it plays an important role in the supply of forest ecosystem services (FES, e.g., water regulation, soil erosion control, carbon sequestration, climate buffering, and timber production).However, soil drought during the growing season in these dryland regions is often a critical limit to tree growth and FES supply.Because of the soil water redistribution due to the unevenly distributed trees and other driving factors, the spatio-temporal difference in soil moisture is often significantly strengthened on hillslopes [30].However, little is known about the temporal stability of soil moisture patterns on such forested hillslopes.This has greatly limited the accurate understanding of the effect of spatio-temporal variation of soil moisture on the vegetation distribution pattern and tree growth on hillslopes.Several studies have reported that the controlling factors of the temporal stability of soil moisture varied with wetness and seasons [31], thus, the temporal stability of volumetric soil moisture (VSM) pattern identified with the data from an entire observation period (e.g., one year) might not be suitable to explain the temporal stability of soil moisture patterns within a certain time (e.g., the growing season) [32].Our previous studies have shown that the soil moisture during the growing season played an important role for the stem radial growth of larch in this area [33].Therefore, clarifying the temporal stability of soil moisture patterns during the growing season is vitally important for the integrated forest-soil-water management.
This study was carried out on a representative larch plantation on a slope in the LPM area, through the continuous monitoring of soil moisture at different root-zone layers at 48 points on a slope in the growing season, with the aims of: (1) clarifying the spatial variability of root-zone soil moisture; (2) analyzing the temporal stability of soil moisture at different soil depths; (3) assessing the methods for estimating the mean soil moisture of the entire slope; and (4) identifying the major factors influencing the spatial patterns of soil moisture.

Study Area
The study was conducted in the small watershed of Xiangshuihe (XSH, 106 • 12 -106 • 16 E, 35 • 27 -35 • 33 N) in the LPM area, Ningxia Hui Autonomous Region, northwest China (Figure 1).XSH has an area of 43.7 km 2 , an elevation range of 2010-2942 m above sea level (a.s.l.), a temperate semi-humid climate with mean air temperature of 6.0 • C and mean annual precipitation of 632 mm.The main vegetation types are natural secondary forests and plantations.The natural secondary forests mainly consist of Pinus armandii Franch, Betula platyphylla Suk., Betula albosinensis Burk., and Quercus liaotungensis Koidz.The plantations consist mainly of Larix principis-rupprechtii, which accounts for 90% of the total plantation area, and a small portion of Pinus tabuliformis Carrière.
A southeast-facing and semi-sunny hillslope, with a slope length of 480.6 m and covered by a 35-year-old pure larch plantation was selected in XSH (Figure 1).The mean slope gradient is 27.8 • , the elevation ranges from 2259 to 2498 m a.s.l., and the soil thickness is 1-1.2 m.The mean tree height (H) on the slope was 17.1 m, and the mean diameter at breast height (DBH) was 19.83 cm.The tree H and DBH presented a remarkable difference along slope positions, with their maximum at mid-slope.The shrub layer and herb layer had low coverage under forest canopy, with coverage values of 15% and about 40% respectively, because of the high canopy density of 0.74.The main shrub species were Viburnum mongolicum (Pall.)Rehd., Prunus salicina Lindl., and Cotoneaster zabelii C.K. Schneid.The main herb species were Fragaria orientalis Losinsk and Carex hancockiana (Maxim.).
Forests 2018, 9, x FOR PEER REVIEW 4 of 22 liaotungensis Koidz.The plantations consist mainly of Larix principis-rupprechtii, which accounts for 90% of the total plantation area, and a small portion of Pinus tabuliformis Carrière.A southeast-facing and semi-sunny hillslope, with a slope length of 480.6 m and covered by a 35-year-old pure larch plantation was selected in XSH (Figure 1).The mean slope gradient is 27.8°, the elevation ranges from 2259 to 2498 m a.s.l., and the soil thickness is 1-1.2 m.The mean tree height (H) on the slope was 17.1 m, and the mean diameter at breast height (DBH) was 19.83 cm.The tree H and DBH presented a remarkable difference along slope positions, with their maximum at mid-slope.The shrub layer and herb layer had low coverage under forest canopy, with coverage values of 15% and about 40% respectively, because of the high canopy density of 0.74.The main shrub species were Viburnum mongolicum (Pall.)Rehd., Prunus salicina Lindl., and Cotoneaster zabelii C.K. Schneid.The main herb species were Fragaria orientalis Losinsk and Carex hancockiana (Maxim.).

Sampling Points and Measurement of Soil and Vegetation Properties
Polycarbonate tubes were installed at 48 points to monitor soil moisture along the hillslope.The distance between the upper and lower adjacent tubes was about 30 m, and the distance between the neighboring tubes at the same level was about 15 m (Figure 1).The volumetric soil moisture (VSM) within the root zone (0-60 cm) at each point was monitored at a space interval of 20 cm and a time interval of about 15 days, using a time domain reflectometry (TDR) portable soil moisture measuring instrument (TRIME-PICO-IPH, Ettlingen, Germany).The VSM data of totally 13 observation days were collected from May to October in 2016.According to Brocca et al. [34], generally 12 sampling occasions are required to determine the most time-stable locations (MTSLs).Thus, the 13 sampling

Sampling Points and Measurement of Soil and Vegetation Properties
Polycarbonate tubes were installed at 48 points to monitor soil moisture along the hillslope.The distance between the upper and lower adjacent tubes was about 30 m, and the distance between the neighboring tubes at the same level was about 15 m (Figure 1).The volumetric soil moisture (VSM) within the root zone (0-60 cm) at each point was monitored at a space interval of 20 cm and a time interval of about 15 days, using a time domain reflectometry (TDR) portable soil moisture measuring instrument (TRIME-PICO-IPH, Ettlingen, Germany).The VSM data of totally 13 observation days were collected from May to October in 2016.According to Brocca et al. [34], generally 12 sampling occasions are required to determine the most time-stable locations (MTSLs).Thus, the 13 sampling times in our study should be considered sufficient for analyzing the temporal stability of VSM.The weather conditions were monitored using an automatic weather station (WeatherHawk 232, WeatherHawk, Logan, Utah) about 100 m away from the slope foot at the relatively open site of the valley.The daily precipitation distribution during the study period is shown in Figure 2. times in our study should be considered sufficient for analyzing the temporal stability of VSM.The weather conditions were monitored using an automatic weather station (WeatherHawk 232, WeatherHawk, Logan, Utah) about 100 m away from the slope foot at the relatively open site of the valley.The daily precipitation distribution during the study period is shown in Figure 2. Soil samples were collected near each tube, by using a ring knife with a volume of 200 cm 3 , at the 0-20, 20-40, and 40-60 cm soil layers, to measure the bulk density, field capacity, and capillary, noncapillary, and total porosity.The elevation at each tube was measured using GPS (GPSMAP 631SC, Garmin, Olathe, KS, USA).The leaf area index (LAI) above each point was measured every 15 days using a plant canopy analyzer (LAI-2200C, LI-COR, Lincoln, NE, USA).The descriptive statistics for soil properties, LAI, slope gradient, and elevation across 48 observation points are shown in Table 1.Soil samples were collected near each tube, by using a ring knife with a volume of 200 cm 3 , at the 0-20, 20-40, and 40-60 cm soil layers, to measure the bulk density, field capacity, and capillary, non-capillary, and total porosity.The elevation at each tube was measured using GPS (GPSMAP 631SC, Garmin, Olathe, KS, USA).The leaf area index (LAI) above each point was measured every 15 days using a plant canopy analyzer (LAI-2200C, LI-COR, Lincoln, NE, USA).The descriptive statistics for soil properties, LAI, slope gradient, and elevation across 48 observation points are shown in Table 1.

Spearman's Rank Correlation Coefficient
Spearman's rank correlation coefficient (r s ) was used to analyze the similarity of the spatial patterns of VSM among sampling dates [16], which is defined as: where n is the number of sampling points, R ij is the rank of VSM at point i on date j, and R ik is the rank of VSM at point i on date k.The r s of 1, for any point, represents the VSM having the same rank between sampling date j and k, and indicates the perfect temporal stability.Thus, higher values of r s , or values closer to 1, represent higher temporal stability of spatial patterns of VSM between sampling dates.

Frequency Distribution
The frequency distribution of VSM for all sampling points was computed to determine whether or not a given point kept its rank in the frequency distribution for different sampling dates [35].The points maintaining the same rank in the frequency distribution for different sampling dates can be considered to have high temporal stability of VSM.If the point with a probability of 0.5 maintained the same rank in different sampling dates, the VSM of this point can be used to represent the mean VSM [16].

Relative Difference Analysis
The temporal stability of VSM at the point scale can also be evaluated by a relative difference analysis [16].The relative difference is calculated by where θ ij is the VSM at point i on date j, θ j is the mean VSM of all sampling points at sampling date j.
The mean relative difference (MRD i ) of VSM for each point and the associated standard deviation (SDRD i ) are calculated by where m is the number of measuring occasions.A point with lower SDRD has higher temporal stability.Typically, the points with MRD values close to zero and with smaller SDRD values (<5%) can be considered as the most time-stable locations (MTSLs) to represent the slope mean VSM.

The Estimation of Slope Mean VSM by Direct and Indirect Methods
According to Jacobs et al. [17] and Zhao et al. [36], the index of temporal stability (ITS), which is a combination of MRD and associated SDRD, can be used to identify the time-stable points.A point with lower ITS reflects higher temporal stability and can directly represent the slope mean VSM (direct method).The ITS can be expressed as Forests 2018, 9, 68 The non-zero MRD points with the minimum mean absolute bias error (MABE) can be considered as time-stable points that can be used to indirectly estimate the mean hillslope VSM (indirect method) [20].The MABE can be expressed as: and the mean slope VSM (θ) can be calculated using the following equation: where θ i is the VSM of the time-stable points with non-zero MRD i .

Statistical Analysis
Two statistical indicators (Nash-Sutcliffe coefficient of efficiency, NSCE; and root mean square error, RMSE) [37] were introduced to estimate the prediction accuracy of mean slope VSM by direct and indirect methods.
A Pearson correlation analysis performed by SPSS 19.0 software (IBM SPSS., Chicago, IL, USA) was used to determine the relationship between MRD and topographic factors (elevation and slope gradient), soil properties (bulk density, field capacity, and capillary, non-capillary, and total porosity), and vegetation structure (LAI).

SDs and CVs versus Slope Mean VSM
The relationship between CVs and slope mean VSM varied with soil depth.CVs increased firstly and then decreased at the depth of 0-20 cm, and decreased at the depths of 20-40 cm and 40-60 cm as the slope mean VSM increased (Figure 4A).The relationship between CVs and slope mean VSM depended on the slope mean VSM value as a whole; CVs increased as the spatial mean VSM increased when the slope mean VSM was lower than 15%, and decreased as the slope mean VSM increased when the slope mean VSM was higher than 15%.The SDs were positively and linearly correlated with the slope mean VSM for all soil layers (Figure 4B), but the correlation became weaker with rising soil depth (with Pearson's correlation coefficients of 0.957, 0.948, and 0.810 for the 0-20, 20-40, and 40-60 cm soil layers, respectively).

At Point Scale Frequency Distribution
Frequency distribution analysis can determine if a particular point maintains its rank at different sampling dates.The maximum spatial mean VSM for the 0-60 cm soil layer was on the 1 May, with the value of 26.91%; while the minimum (15.41%) was on 30 September.The cumulative probability functions for these two extreme days are given in Figure 6.It shows that only a few points from the 48 points maintained the same rank.They are points 28, 36, and 41 for the 0-20 cm soil layer; points 24, 29, and 38 for the 20-40 cm soil layer; and points the 12, 14, 16, 18, and 24 for 40-60 cm soil layer.No points with the probability of 0.5 maintained the same rank under the two extreme soil moisture conditions for all three soil layers.Frequency distribution analysis can determine if a particular point maintains its rank at different sampling dates.The maximum spatial mean VSM for the 0-60 cm soil layer was on the 1 May, with the value of 26.91%; while the minimum (15.41%) was on 30 September.The cumulative probability functions for these two extreme days are given in Figure 6.It shows that only a few points from the 48 points maintained the same rank.They are points 28, 36, and 41 for the 0-20 cm soil layer; points 24, 29, and 38 for the 20-40 cm soil layer; and points the 12, 14, 16, 18, and 24 for 40-60 cm soil layer.No points with the probability of 0.5 maintained the same rank under the two extreme soil moisture conditions for all three soil layers.
The number of points with a VSM close to the slope mean varied with soil depth.For example, the number of the points with the MRD range from −5% to 5% was 9, 7, and 10 for the 0-20, 20-40, and 40-60 cm soil layers, respectively.The number of time-stable points was also different in diverse soil depths, the number of points with SDRD values less than 5% was 0, 2, and 15 for the 0-20, 20-40, and 40-60 cm soil layers, respectively.
The number of points with a VSM close to the slope mean varied with soil depth.For example, the number of the points with the MRD range from −5% to 5% was 9, 7, and 10 for the 0-20, 20-40, and 40-60 cm soil layers, respectively.The number of time-stable points was also different in diverse soil depths, the number of points with SDRD values less than 5% was 0, 2, and 15 for the 0-20, 20-40, and 40-60 cm soil layers, respectively.

Identification of Representative Point for Acquisition of Mean VSM
According to the relative difference method (MRD close to zero), the points that were closest to the slope mean VSM were point 18 (MRD = 0.15% ± 13.69%) for the 0-20 cm soil layer with the SDRD of 13.69%, point 14 (MRD = 1.12% ± 5.64%) for the 20-40 cm soil layer with the SDRD of 5.64%, and point 29 (MRD = 0.04% ± 14.16%) for the 40-60 cm soil layer with the SDRD of 14.16%, respectively.This indicates that the points with MRDs close to zero cannot guarantee that their SDRDs will be minimum values (<5%) and might not directly represent the slope mean VSM due to their lower temporal stability.Thus, the direct method was used to select the average VSM of three points with the lowest ITS values for representing the slope mean VSM.These points were point 10, 16, and 27 for the soil layer of 0-20 cm (Figure 8A), 9, 14, and 22 for the soil layer of 20-40 cm (Figure 8B), and 5, 14, and 40 for the soil layer of 40-60 cm (Figure 7C), respectively.The estimated accuracy of mean VSM using these points (Figure 9A-C) increased with soil depth; the NSCE and RMSE were 0.95 and 0.91% for 0-20 cm soil layer; 0.98 and 0.63% for 20-40 cm soil layer; and 0.97 and 0.53% for 40-60 cm soil layer, respectively.An estimation with a RMSE less than 2% can be considered to be reliable [39], thus the average VSM at these points can be used to represent the mean VSM for the different soil layers for the studied larch plantation on the hillslope.

Identification of Representative Point for Acquisition of Mean VSM
According to the relative difference method (MRD close to zero), the points that were closest to the slope mean VSM were point 18 (MRD = 0.15% ± 13.69%) for the 0-20 cm soil layer with the SDRD of 13.69%, point 14 (MRD = 1.12% ± 5.64%) for the 20-40 cm soil layer with the SDRD of 5.64%, and point 29 (MRD = 0.04% ± 14.16%) for the 40-60 cm soil layer with the SDRD of 14.16%, respectively.This indicates that the points with MRDs close to zero cannot guarantee that their SDRDs will be minimum values (<5%) and might not directly represent the slope mean VSM due to their lower temporal stability.Thus, the direct method was used to select the average VSM of three points with the lowest ITS values for representing the slope mean VSM.These points were point 10, 16, and 27 for the soil layer of 0-20 cm (Figure 8A), 9, 14, and 22 for the soil layer of 20-40 cm (Figure 8B), and 5, 14, and 40 for the soil layer of 40-60 cm (Figure 7C), respectively.The estimated accuracy of mean VSM using these points (Figure 9A-C) increased with soil depth; the NSCE and RMSE were 0.95 and 0.91% for 0-20 cm soil layer; 0.98 and 0.63% for 20-40 cm soil layer; and 0.97 and 0.53% for 40-60 cm soil layer, respectively.An estimation with a RMSE less than 2% can be considered to be reliable [39], thus the average VSM at these points can be used to represent the mean VSM for the different soil layers for the studied larch plantation on the hillslope.
In order to further clarify the distribution characteristics of time-stable points characterized by ITS on our study hillslope, the relationship between ITS and topographic, vegetation, and soil factors for the soil layers of 0-20 cm, 20-40 cm, and 40-60 cm are shown in Figure 9, with the three points with the lowest ITS values labelled.As compared to the slope means, the slope gradient, elevation, LAI, field capacity, and soil bulk density of the three time-stable points, which can represent the mean VSM, fell in a broad range (Figure 10A-O), but the ratios of field capacity to LAI of these time-stable points were close to the slope means (Figure 10P-R).
We then estimated the slope mean VSM at different soil layers based on Equation (7), relying on the time-stable points by selecting three points having the lowest MABE but with non-zero MRD (indirect method).The selected points were point 20, 21, and 27 for the soil layer of 0-20 cm (Figure 8A), 16, 20, and 26 for the soil layer of 20-40 cm (Figure 8B), and 9, 16, and 48 for the soil layer of 40-60 cm (Figure 8C), respectively.Figure 9D-F showed that the indirect method can provide a good estimate of the slope mean VSM.The NSCE and RMSE were 0.99 and 0.49%, 0.98 and 0.56%, and 0.99 and 0.34%, for the soil layers of 0-20, 20-40, and 40-60 cm, respectively.These NSCE values were higher and these RMSE values were lower than that by the direct method.Thus, the estimation by indirect method was better than the direct method.
In order to further clarify the distribution characteristics of time-stable points characterized by ITS on our study hillslope, the relationship between ITS and topographic, vegetation, and soil factors for the soil layers of 0-20 cm, 20-40 cm, and 40-60 cm are shown in Figure 9, with the three points with the lowest ITS values labelled.As compared to the slope means, the slope gradient, elevation, LAI, field capacity, and soil bulk density of the three time-stable points, which can represent the mean VSM, fell in a broad range (Figure 10A-O), but the ratios of field capacity to LAI of these time-stable points were close to the slope means (Figure 10P-R).
We then estimated the slope mean VSM at different soil layers based on Equation ( 7), relying on the time-stable points by selecting three points having the lowest MABE but with non-zero MRD (indirect method).The selected points were point 20, 21, and 27 for the soil layer of 0-20 cm (Figure 8A), 16, 20, and 26 for the soil layer of 20-40 cm (Figure 8B), and 9, 16, and 48 for the soil layer of 40-60 cm (Figure 8C), respectively.Figure 9D-F showed that the indirect method can provide a good estimate of the slope mean VSM.The NSCE and RMSE were 0.99 and 0.49%, 0.98 and 0.56%, and 0.99 and 0.34%, for the soil layers of 0-20, 20-40, and 40-60 cm, respectively.These NSCE values were higher and these RMSE values were lower than that by the direct method.Thus, the estimation by indirect method was better than the direct method.

Factors Affecting the Spatial Pattern of VSM
To identify which factor (elevation, slope gradient, soil bulk density, field capacity, capillary, non-capillary, and total porosity, and LAI) controlled the spatial pattern of VSM, the relationships between MRD and these factors were determined (Table 3).A significant correlation between MRD and elevation was observed only in the 20-40 cm soil layer, but not in the other soil layers.No significant correlation between MRD and slope gradient was found for all soil layers.The soil bulk density, field capacity, and capillary porosity were significantly correlated with MRD for all soil layers.A significant correlation between MRD and soil total porosity existed only in the 0-20 cm soil layer, but not the other layers.Soil non-capillary porosity was only significantly negatively correlated with MRD in the 40-60 cm soil layer, but not the other layers.LAI was significantly negatively correlated with MRD for all soil depths.

Spatial and Temporal Variability of VSM
The obvious spatial variation of VSM in each soil layer at any given date showed moderate variability; this is similar to the study results on forested slopes by others [22,40], and can be attributed to the spatial heterogeneity of influencing factors, including soil properties, topography, and canopy structure [11].The temporal variability of the slope mean VSM decreased with rising soil depth; this is similar to that reported under other land use types, such as grassland [41], gravel-mulched field [21], and artificial re-vegetated desert [42].Higher temporal variation in spatial variability was also observed in the shallow layer.The larger variation in the shallow layer might be ascribed to a greater vulnerability to the effects of environment variation (i.e., evaporation, precipitation).
The relationship between CVs and mean VSM varied in different studies: both decreases [21,43,44] and increases of CVs [45,46] with rising mean VSM were reported.In our study, we observed a VSM threshold below which the CVs correlated with the mean VSM positively but above which they correlated negatively.This difference was mainly caused by the variation range of VSM in different studies.In fact, the relationship between CVs and mean VSM is dominated by soil texture because of its influence on soil hydraulic properties (i.e., evaporation and drainage) [47].The CVs are mainly controlled by soil hydraulic conductivity and porosity under (extreme) wet conditions.When the mean VSM starts to decrease due to the drainage and evapotranspiration, the variability of saturated hydraulic conductivity, particle size distribution, and air entry pressure will result in an increase in CVs, thus CVs are negatively correlated with the mean VSM in relatively wet conditions [47,48].When the mean VSM decreases to a threshold, the dominant process controlling soil drying is switched from drainage to evapotranspiration.With further reduction of mean VSM, CVs are decreased by the effect of evapotranspiration, showing a positive correlation with mean VSM in dry soils [48].The VSM threshold depends on soil, topography, climate, and vegetation, and is typically between 15% and 25% [48,49].The VSM threshold (about 15%) in our study was within this range.
The positive correlation between SDs and mean VSM in all soil layers in this study is in accordance with the findings of others [21,50,51], but in contrast with some results [52,53].This difference might be ascribed to the different soil water statuses among those studies [54].Besides, we found that the correlation of SDs with mean VSM became weaker with increasing soil depth.The lower dependency with soil depth indicated that the role of soil moisture condition in determining the spatial heterogeneity of VSM becomes weak with increasing soil depth.Gao et al. [50], however, found this correlation increased with soil depth.The different results might be related to the differences of vegetation cover.The soil water consumption in the grassland in the study of Gao et al. [50] was mainly due to the soil surface evaporation, which typically had a decreasing impact on VSM with rising soil depth; while the soil water consumption in the forestland in our study was mainly due to tree transpiration, which typically had an increasing impact on VSM with rising soil depth within the root zone.

Temporal Stability of VSM
Spearman's rank correlation coefficient can reflect the similarity of spatial patterns of VSM among observation dates [50].The significant Spearman's rank correlation coefficients for all soil layers in this study indicate a strong temporal stability of the VSM on the larch planation on the hillslope, similar to other forested hillslopes, such as hillslopes forested with Pinus tabulaeformis Carr.[22,40] and Robinia pseudoacacia L. [11].The decrease of Spearman's rank correlation coefficient with rising soil depth indicates a temporally more stable VSM in deep soil layers than in shallow layers.This result is consistent with most previous studies [21,42,50].The reasons for the stronger temporal stability with soil depth can be explained by three factors.Firstly, the effects of terrain and water uptake by vegetation roots on the soil moisture in deeper soil layers remained relatively stable across time [55].Secondly, the dynamics of water retention and soil structure are much more pronounced at greater soil depths, thus enhancing temporal stability with increasing soil depth [56].Thirdly, the seasonal dynamics of understory vegetation caused more variability of VSM for the shallower layer.
SDRD can reflect the degree of temporal stability for a given point.The decrease of mean SDRD of 48 points with rising soil depth indicates a stronger temporal stability of VSM in deep soil layers than in shallow soil layers.The relationship between the temporal stability and soil depth in our study was consistently indicated by the SDRD and Spearman's rank correlation coefficient.This is consistent with the result reported by Zhao et al. [21].However, in the watershed scale study, Hu et al. [20] found a contrasting result in that the temporal stability of the points in the shallow soil layer (20 cm) was weaker than in the deeper soil layers (60 and 80 cm), but the similarity of spatial patterns in the shallow soil layer was stronger than in the deeper soil layers.This difference might be related to the different study scales, as there should be more factors (e.g., land use cover, climate, and topography) influencing the soil water movement and soil moisture variation at the watershed scale than at the slope scale.

Representative Points for Mean VSM
Vachaud et al. [16] reported that the points associated with a probability of 0.5 could be used to represent the mean VSM for an entire field, but none of the points with a probability of 0.5 maintained the same rank under the two extreme moisture conditions for all soil layers in our study hillslope, indicating that estimating the mean VSM based on the analysis of frequency distribution was not a workable solution for our larch hillslope.The rank of a given point, especially the time-stable point, underwent few variations when soil wetness among different sampling dates had little fluctuation (e.g., dry or wet soil conditions), but, its rank changed a lot when soil wetness transitioned from dry to wet [24].The fact that none of the points with a probability of 0.5 demonstrated this characteristic (the same rank) in our larch hillslope might be related to the larger fluctuation of VSM during the study period (Figure 3A).Similar results were also reported in other land use types (e.g., gravel-mulched field, re-vegetated abandon farmland) with larger fluctuation of soil moisture during the experimental period [10,21].
The points with the MRD close to zero and minimal SDRD (<5%) can be selected as the MTSLs to estimate the slope mean VSM.However, in this study, no such point was found due to the high spatial heterogeneity on the larch plantation slope [57][58][59].A similar result was found in a study on a shrub slope in a karst area [37].Thus, in some studies [17,36], the time-stable points with the minimum ITS value, which considers both the MRD and SDRD, was selected.Such selection with the direct method yielded a good accuracy of mean VSM estimation in this study, but also in other studies involving mulched field [21], karst slopes with shrub and grass [37], and desert area [60].This indicates that the slope mean VSM can be directly estimated by the mean VSM value of several representative points.Several studies have found that the representative points (time-stable points) that can directly estimate the mean soil moisture had some particular distribution rules.For example, Jacobs et al. [17] reported that the points with 27% clay content had the best time-stable feature at the watershed scale; Joshi et al. [61] found that the time-stable points were located at higher elevations at the airborne remote sensing footprint scale (800 m × 800 m).In our study at hillslope scale, the representative points were located at the points with a ratio of field capacity to LAI close to the slope mean.This implies the representative points might be jointly determined by the soil's water retention capacity (field capacity) and the forest's water consumption ability (LAI), and thus these factors should be considered when identifying the representative locations in this area.The fact that determining factors for representative points (time-stable points) were different with the various scales might be attributed to the complexity and diversity of the factors determining the spatio-temporal variation of soil moisture at different scales.
Several studies have reported that the points with the minimum MABE with non-zero MRD can indirectly estimate the mean VSM, and have a higher accuracy than the direct method [6,62].A similar result was found in our study, since the prediction errors by indirect method were lower than that by direct method.The reason may be that the mean MRD of time-stable points based on the direct method deviated slightly from zero (1.41%, 0.66%, and −0.66% for the 0-20, 20-40, and 40-60 cm soil layers), leading to a deviation of the VSM of time-stable points from the real slope mean VSM; whereas the indirect method can eliminate this deviation by introducing a constant offset (MRD) [37].Therefore, if no point with the MRD close to zero and a minimal SDRD (<5%) can be selected, the estimation of slope mean VSM based on the indirect method could be considered.
It was noted that no single point could represent the mean VSM for all three soil layers.This agreed with the findings of Guber et al. [63] and Gao et al. [50] for agricultural field and grassland areas, but disagreed with the findings of Zhao et al. [21] for a gravel-mulched field.The reasons for different representative points at different soil depths can be explained by two factors.Firstly, the vertical variability of soil properties can have a critical impact on the temporal dynamics of soil moisture [36], while the soil properties (e.g., bulk density and field capacity) on our forested hillslope had strong vertical variability (Table 1).Secondly, the temporal stability of VSM was related to the combined effects of terrain, soil, and vegetation, however, these factors, especially the roots of vegetation in forestland, had diverse effects on VSM for different soil depths [50].
Based on the representative points, the slope mean VSM can be derived with a lower time cost.Thus, it might be more beneficial to estimate the slope mean VSM during shorter campaigns.The representative points can also be used to estimate the VSM for other growing seasons with similar environmental conditions.However, the temporal stability of VSM depends on the mean soil wetness, so the representative points determined based on the observation of one growing season of a dry hydrographic year in this study cannot be used to estimate the slope mean VSM in other hydrographic years (e.g., wet hydrographic years).Therefore, for long-term monitoring and estimation, it should be better to determine the representative points based on the data including all of the three hydrographic (normal, wet, and dry) years.

Factors Affecting the Spatial Pattern of VSM
Topography plays an important role in the spatial pattern of VSM.Elevation is one topographic factor that may affect the spatial variation of VSM, as reported in several studies [11,64], especially under complex terrains conditions [27,35,37,41].However, in our study, a significant correlation between MRD and elevation was observed only in the 20-40 cm soil layer, similar to the result observed in a study on terraced land [22].In fact, elevation is linked with many other factors (e.g., potential evapotranspiration, rainfall, soil properties, and tree growth characteristics) and many hydrological processes (e.g., soil evaporation, interception, and tree transpiration).The complex interaction of these factors and processes might be one of the reasons leading to the weak correlation with elevation in our study.In addition, the vegetation cover could weaken the effect of topography (e.g., elevation) on the VSM, because the unevenly distributed trees and roots can influence soil water redistribution, and then complicate this correlation [65,66].Slope gradient is also one topographic factor that can have a great influence on the spatio-temporal variation of VSM through its effect on runoff and soil water movement [28].However, no significant correlation between MRD and slope gradient was found in this study; this may be explained by the low surface runoff and interflow in our forested hillslope, which give less chance to form the effect of slope gradient.
Soil properties are important factors influencing the spatial variation of VSM [17,36,60].For example, Cosh et al. [39] reported that soil bulk density could explain 31.8% of MRD variability at the watershed scale.Wang et al. [37] found that the MRD was significantly correlated with soil bulk density and saturated hydraulic conductivity on karst slopes.In our study, significant correlations between MRD and soil bulk density and field capacity were observed.Soil bulk density and field capacity can directly affect soil hydraulic properties such as infiltration and evaporation, and then the spatial variation of VSM [67,68].This indicates that the soil water retention capacity may have played an important role in the spatial variation of VSM on our forested hillslope.Capillary porosity is closely related to the soil water retention capacity (p < 0.01), and thus was significantly correlated with MRD.Soil porosity was significantly correlated with MRD only in the 0-20 cm soil layer, but not in the other soil layers.This may be caused by the effect of non-capillary porosity, which does not affect the water retention capacity of soil (there was no significant correlation between non-capillary porosity and field capacity).
Canopy structure has a strong influence on the spatial variation of VSM at both plot and slope scales [10,36,69].For example, Duan et al. [11] reported that MRD was significantly negatively correlated with LAI in Robinia pseudoacacia forestland on the Loess Plateau of China.As the most important canopy structure, higher LAI leads to higher evapotranspiration rates and water consumption, and then lower soil moisture, especially in the case of insufficient soil moisture, thus showing a negative correlation with MRD [31,70,71].Similarly, a significant negative correlation between MRD and LAI was observed for all soil layers in this study; the remarkable spatial LAI difference on our forested hillslope and lower soil moisture during the study period should be the major causes [59].
The above results show that the spatial variation of VSM was mainly controlled by the combined effects of soil and vegetation.Teuling et al. [31] also reported that soil and vegetation controlled the spatial variation of soil moisture, and the main discriminating factor between both factors depends on the soil wetness.Vegetation (e.g., LAI) might have a more pronounced effect on the spatial pattern of soil moisture in the "dry" domain, and soil (e.g., soil properties) in the "wet" domain [71].Our unpublished paper also found that soil properties, not LAI, were the major factors affecting the spatial variation of soil moisture in a year with sufficient precipitation (2015).In this study, the correlation between MRD and LAI is higher than that between MRD and soil properties (e.g., bulk density, field capacity).The precipitation during the growing season of 2016 was 413 mm, which was significantly lower than the long-term mean of 550 mm in the period of 1970-2010; thus, the soil moisture during the growing season of the study year was relatively lower.Therefore, the LAI played a more important role in the spatial variation of VSM compared with soil properties, and this might be ascribed to the lower soil moisture during the study period of 2016.

Conclusions
The spatial variation of VSM for each soil layer (0-20, 20-40, and 40-60 cm) of the root zone showed moderate variability.The spatial variability of VSM increased with rising VSM until a threshold of about 15%, and decreased thereafter.The spatial distribution of VSM in all three soil layers had a strong temporal stability, but increased with rising soil depth.The slope mean VSM estimation using the VSM of representative points selected by ITS (direct method) was successful for the larch plantation slope, and these representative points were mainly located at the points with a ratio of field capacity to LAI close to the slope mean.Moreover, the slope mean VSM can also be estimated by indirect method (using the time-stable points selected by MABE and introducing a constant offset (MRD)), and the prediction accuracy of the indirect method was higher than the direct method.If no representative point with a MRD close to zero and a minimal SDRD (<5%) can be selected, the slope mean VSM can be estimated by the indirect method (on the premise that the temporal stability of soil moisture is known).The water retention capacity of soil (mainly field capacity) and the water consumption ability of forest (mainly canopy LAI) were the main processes (factors) controlling the spatio-temporal variation of the root-zone VSM in the larch plantation slope.

Figure 1 .
Figure 1.Location of study area and sampling points across the larch plantation on the slope.

Figure 1 .
Figure 1.Location of study area and sampling points across the larch plantation on the slope.

Figure 2 .
Figure 2. The distribution of daily precipitation during the study period in 2016.

Figure 2 .
Figure 2. The distribution of daily precipitation during the study period in 2016.

Figure 3 .
Figure 3.The temporal dynamics of slope means of VSM (A) and associated coefficients of variation (B) and standard deviations (C) in different soil layers.

Figure 3 .
Figure 3.The temporal dynamics of slope means of VSM (A) and associated coefficients of variation (B) and standard deviations (C) in different soil layers.

Forests 2018, 9 , 68 9 of 22 22 Figure 4 .
Figure 4.The relationships between the slope mean VSM and the coefficient of variation (A), and the standard deviation (B).

3. 3 .
Temporal Stability of VSM 3.3.1.At Slope Scale The temporal stability of the spatial patterns of VSM was determined based on the Spearman's rank correlation coefficients (rs) between different sampling dates.The spatial patterns of VSM at different sampling dates were all significantly correlated (p < 0.05), and all the average rank correlation coefficients of VSM corresponding to different time lags at three different depths are greater than 0.5 (Figure 5), indicating a strong temporal persistence of the spatial distribution of VSM in various depths.The temporal stability of the spatial distribution of VSM increased with rising soil depth, with the mean rs values over the study period of 0.65, 0.78, and 0.84 for the 0-20, 20-40, and 40-60 cm soil layers, respectively.A one-way ANOVA showed that the mean rs among the three soil depths were significantly different (p < 0.01).

Figure 5 .
Figure 5. Mean Spearman's rank correlation coefficients (rs) of VSM corresponding to different time lags in different depths.

Figure 4 .
Figure 4.The relationships between the slope mean VSM and the coefficient of variation (A), and the standard deviation (B).

3. 3 .
Temporal Stability of VSM3.3.1.At Slope ScaleThe temporal stability of the spatial patterns of VSM was determined based on the Spearman's rank correlation coefficients (r s ) between different sampling dates.The spatial patterns of VSM at different sampling dates were all significantly correlated (p < 0.05), and all the average rank correlation coefficients of VSM corresponding to different time lags at three different depths are greater than 0.5 (Figure5), indicating a strong temporal persistence of the spatial distribution of VSM in various depths.The temporal stability of the spatial distribution of VSM increased with rising soil depth, with the mean r s values over the study period of 0.65, 0.78, and 0.84 for the 0-20, 20-40, and 40-60 cm soil layers, respectively.A one-way ANOVA showed that the mean r s among the three soil depths were significantly different (p < 0.01).

Forests 2018, 9 , 22 Figure 4 .
Figure 4.The relationships between the slope mean VSM and the coefficient of variation (A), and the standard deviation (B).

3. 3 .
Temporal Stability of VSM 3.3.1.At Slope Scale The temporal stability of the spatial patterns of VSM was determined based on the Spearman's rank correlation coefficients (rs) between different sampling dates.The spatial patterns of VSM at different sampling dates were all significantly correlated (p < 0.05), and all the average rank correlation coefficients of VSM corresponding to different time lags at three different depths are greater than 0.5 (Figure 5), indicating a strong temporal persistence of the spatial distribution of VSM in various depths.The temporal stability of the spatial distribution of VSM increased with rising soil depth, with the mean rs values over the study period of 0.65, 0.78, and 0.84 for the 0-20, 20-40, and 40-60 cm soil layers, respectively.A one-way ANOVA showed that the mean rs among the three soil depths were significantly different (p < 0.01).

Figure 5 .
Figure 5. Mean Spearman's rank correlation coefficients (rs) of VSM corresponding to different time lags in different depths.Figure 5. Mean Spearman's rank correlation coefficients (r s ) of VSM corresponding to different time lags in different depths.

Figure 5 .
Figure 5. Mean Spearman's rank correlation coefficients (rs) of VSM corresponding to different time lags in different depths.Figure 5. Mean Spearman's rank correlation coefficients (r s ) of VSM corresponding to different time lags in different depths.

Figure 6 .
Figure 6.The cumulative frequency of 48 points at the two days with maximum and minimum VSM in the soil layers of 0-20 cm (A), 20-40 cm (B), and 40-60 cm (C).

Figure 6 .
Figure 6.The cumulative frequency of 48 points at the two days with maximum and minimum VSM in the soil layers of 0-20 cm (A), 20-40 cm (B), and 40-60 cm (C).

ForestsFigure 9 .
Figure 9.The measured and predicted slope mean VSM of each of the sampling dates for different soil layers using direct (A-C) and indirect (D-F) methods.

Figure 9 .
Figure 9.The measured and predicted slope mean VSM of each of the sampling dates for different soil layers using direct (A-C) and indirect (D-F) methods.

Figure 10 .
Figure 10.Relationships between ITS and slope gradient (A-C), elevation (D-F), LAI (G-I), soil field capacity (J-L), soil bulk density (M-O), and the ratio of field capacity to leaf area index (LAI) (P-R) in different soil layers (0-20, 20-40, and 40-60 cm).The dash lines denote the mean values of affecting factors.

Figure 10 .
Figure 10.Relationships between ITS and slope gradient (A-C), elevation (D-F), LAI (G-I), soil field capacity (J-L), soil bulk density (M-O), and the ratio of field capacity to leaf area index (LAI) (P-R) in different soil layers (0-20, 20-40, and 40-60 cm).The dash lines denote the mean values of affecting factors.

Table 1 .
The statistics for soil properties, leaf area index (LAI), slope gradient, and elevation across 48 observation points.
Statistic Variables Soil Depth Mean SD CVSoil bulk density (cm 3 •cm −3 )Note: Mean is the mean value of each variable for 48 points; SD is the standard deviation of each variable for 48 points; CV is the variation coefficient of each variable for 48 points (%).

Table 1 .
The statistics for soil properties, leaf area index (LAI), slope gradient, and elevation across 48 observation points.Note: Mean is the mean value of each variable for 48 points; SD is the standard deviation of each variable for 48 points; CV is the variation coefficient of each variable for 48 points (%).

Table 2 .
The temporal statistics of slope (spatial) means of volumetric soil moisture (VSM) and associated coefficients of variation (CVs) and standard deviations (SDs) for three soil layers.

Table 2 .
The temporal statistics of slope (spatial) means of volumetric soil moisture (VSM) and associated coefficients of variation (CVs) and standard deviations (SDs) for three soil layers.

Table 3 .
The correlations between mean relative difference (MRD) and topographic, soil physical, and vegetation structure factors for three soil layers.