Rock fragment and spatial variation of soil hydraulic parameters are necessary on soil water simulation on the stony-soil hillslope

Journal of Hydrology 565 (2018) 354–364

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Rock fragment and spatial variation of soil hydraulic parameters are necessary on soil water simulation on the stony-soil hillslope

T

⁎

Xiaoming Laia,b, Qing Zhua,b, , Zhiwen Zhoua,b, Kaihua Liaoa a b

Key Laboratory of Watershed Geographic Sciences, Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences, Nanjing 210008, China University of Chinese Academy of Sciences, Beijing 100049, China

A R T I C LE I N FO

A B S T R A C T

This manuscript was handled by Corradini, Editor-in-Chief, with the assistance of Wei Hu, Associate Editor

Whether the spatial variation and influence of rock fragment contents (RFCs) on soil hydraulic parameters (SHPs) should be jointly considered in the soil water simulation have been less investigated in previous studies. In this study, we tested whether considering these factors were necessary in the soil water simulation on a representative stony-soil hillslope located in Taihu Lake Basin, China. Five schemes of SHPs in HYDRUS-3D model were established. They were (i) MultiRFHet: spatially varied SHPs extracted from Durner’s multimodal retention function based on the observed soil water retention data; (ii) RosHom: spatially uniformed SHPs derived from ROSETTA; (iii) RosHet: spatially varied SHPs derived from ROSETTA; (iv) RosRFCHom: spatially uniformed SHPs derived from ROSETTA and adjusted by RFCs; (v) RosRFCHet: spatially varied SHPs derived from ROSETTA and adjusted by RFCs. Results indicated when the spatial variation and influence of RFCs on SHPs were both considered (MultiRFHet and RosRFCHet), acceptable accuracies (Nash-Sutcliffe efficiency or NSE ≥ 0.58 for MultiRFHet and > 0.17 for RosRFCHet) were achieved in simulating the soil water storage (SWS) variation. Since the MultiRFHet required the calibration by the observed soil water retention data, RosRFCHet was a good alternative in the simulation. On the contrary, the SHPs acquired without considering neither the RFCs nor the spatial variation yielded unacceptable simulation results (NSE general < 0). This demonstrated that spatial variation and influence of RFCs on SHPs should be included in the SWS simulation on this stony-soil hillslope, although it had limited effect in subsurface flow simulation.

Keywords: Rock fragment Stony soils Dual-porosity HYDRUS-3D model

1. Introduction With superiority in time and economic costs and in deriving continuous hydrological processes over space and time, three-dimensional (3D) physically-based modeling has been widely adopted in hydrological researches (e.g., Herbst et al., 2006; Vereecken et al., 2015; Baroni et al., 2017). Numerous models have been developed and applied in different areas over the last decades (e.g., Fiori and Russo, 2007; Fang et al., 2016; Baroni et al., 2017). Sources of uncertainties affecting the model simulation (e.g., model structure, input parameters and dependent variables) have also been investigated and determined (Fu et al., 2015; Vereecken et al., 2015). Among these sources, uncertainty of soil hydraulic parameters (SHPs) was recognized as one of the major obstacles for achieving the optimal simulations (Chirico et al., 2007; Baroni et al., 2017). Traditional approaches in acquiring the SHPs include direct measurement and indirect estimation. Direct laboratory or field measurements of SHPs are costly and time-consuming, thus restrict the

⁎

availability of high spatial resolution SHP data for model simulation (Chirico et al., 2007; Vereecken et al., 2010). Thereby, indirect methods especially the pedotransfer functions (PTFs) have been developed to estimate the SHPs from easily measurable soil physical and chemical properties (Wösten et al., 2001; Vereecken et al., 2010), for example the HYPRES (Wösten et al., 1999) and Vereecken (Vereecken et al., 1989). However, previous studies have demonstrated that applying existing PTFs always introduced substantial uncertainty when they were used outside the datasets they developed (Chirico et al., 2007; Liao et al., 2014). Reliable approaches are needed to cost-effectively estimate the SHPs in areas with different geological, pedological and meteorological backgrounds, for example the stony soil areas. Stony soils are soils containing over 30% in volume of rock fragments (soil particles > 2 mm) as defined by Tetegan et al. (2011). Due to soil development processes and actions of mankind as well as soil erosion, stony soils are worldwide distributed. For example, in Western Europe, stony soils account for 30% of soil resources, and account for 60% in the Mediterranean region (Poesen and Lavee, 1994). In China,

Corresponding author at: Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences, Nanjing 210008, China. E-mail address: [email protected] (Q. Zhu).

https://doi.org/10.1016/j.jhydrol.2018.08.039 Received 16 April 2018; Received in revised form 13 August 2018; Accepted 18 August 2018 Available online 20 August 2018 0022-1694/ © 2018 Elsevier B.V. All rights reserved.

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

the tea and bamboo hillslopes were used to assess the relationships between the SHPs and soil/terrain properties. However, soil water dynamics were only simulated on the tea hillslope since a trench were opened along the boundary of tea and bamboo hillslopes to separate them into two hydrological units. The downslope boundary of the tea hillslope was near a pond. Detailed descriptions of this study hillslope can be found in Liao et al. (2016).

mountains composed of stony soils account for 18% of the national area (Ma and Shao, 2008). Therefore, comprehensively and accurately revealing the hydrological processes of upstream stony-soil hillslope is critical for the protection of the downstream ecosystems. Rock fragment content (RFC) affected the SHPs and consequently affected the hydrological processes in stony soils (Cousin et al., 2003; Ma et al., 2010; Coppola et al., 2013). As rock fragments could reduce pore volume available for containing water and increase the tortuosity of water flow, the water content and hydraulic conductivities of stony soils generally decreased with increasing RFC (Baetens et al., 2009; Hlaváčiková et al., 2016). Rock fragments have an obvious influence on the stony soils’ water retention, especially in the low suction region (Ravina and Magier, 1984; Baetens et al., 2009). The widely used van Genuchten-Mualem model (van Genuchten, 1980; Mualem, 1976) cannot characterize the low suction region of the soil water retention characteristic of stony soils (Vereecken et al., 2010; Assouline and Or, 2013). Recently, several dual-porosity functions have been developed for structured soils (Šimůnek et al., 2003; Arora et al., 2011) and can be used to derive SHPs of stony soils (Ma and Shao, 2008). The one proposed by Durner (1994) (multimodal retention function, abbreviated as MultiRF in this study) is promising since it requires less parameters that are relatively easy to be obtained (Šimůnek et al., 2003; Parajuli et al., 2017). Besides, as existing PTFs did not consider the RFC, empirical functions to modify the PTFs have been built to predict the SHPs in stony soils based on the SHPs of fine earth and RFCs (Bouwer and Rice, 1984; Brakensiek et al., 1986). However, doubts emerged with the applications of these empirical functions (e.g., Ma et al., 2010; Novák et al., 2011). Investigations are still needed to verify performances of the dual-porosity functions and modified PTFs in the soil water simulations for stony soils. The spatial variations of SHPs also affect the simulation of soil water dynamics. The SHPs show spatial variations at different scales and detailed characterizing these variations is a crucial challenge over large areas (Herbst et al., 2006; Baroni et al., 2017). As the discrepancy between scales at which SHPs were measured or predicted and at which hydrological model ran, the upscaling procedure was inevitable, and thus the uncertainty of spatial variations of SHPs generated (Vereecken et al., 2007). To date, many scientists have focused on revealing the uncertainty of spatial variation of SHPs on simulating hydrological process (e.g., Herbst et al., 2006; Fiori and Russo, 2007). For example, Sciuto and Diekkrüger (2010) investigated the impacts of the SHPs’ spatial variations on the water balance and found the spatial resolution has a larger impact on the soil moisture than the aggregation of the soil properties. Baroni et al. (2017) quantified the impact of different types of uncertainties that arose from the unresolved soil spatial variability on simulated hydrological states and fluxes. However, previous studies were generally conducted through virtual experiment, only a few of them based on exact characterizing the actual spatial variation of SHPs of the study regions (e.g., Jin et al., 2015; Fang et al., 2016). The objectives of this study were to test the hypotheses that considering (i) the influence of rock fragments on SHPs and (ii) the spatial variation of SHPs improve the soil water simulation on the stony-soil hillslope. Specifically, the MultiRF was adopted to extract the SHPs by fitting observed soil water retention data. We also investigated whether adjusting the ROSETTA derived SHPs by RFC and adopting spatial variation of SHPs would be necessary.

2.2. Data collection Automatic monitoring systems were installed at four sites on the tea hillslope (TG01, TG02, TG03 and TG04) and four sites on the bamboo hillslope (BF01, BF02, BF03 and BF04) (Fig. 1). At each site and each depth (10- and 30-cm), three EC-5 paired with three MPS-6 sensors (Decagon Devices Inc., Pullman WA, USA) were circled installed to monitor the soil water content and matrix potential respectively. The averages of the measured values of the three EC-5 and MPS-6 sensors were respectively used as the final soil water content and matrix potential data at each site and depth. The rock fragments on this hillslope are usually small (< 1 cm in diameter) and comparatively homogenously distributed in the soil. Therefore, although these sensors were installed in the fine-textured soils, their measuring volumes included both the fine-textured soils and rock fragments. The final measured data by these sensors could reflect the soil water content and matrix potential of the stony soils. Sensors were calibrated for each soil depth and before installation. Before installing these sensors, soil profiles were opened and three soil horizons (A, B and C) were observed within each profile (Liao et al., 2016). After all sensors were installed, the profiles were carefully backfilled. A weather station (Decagon Devices Inc., Pullman WA, USA) was established to record the meteorological data (Fig. 1). Two rain gauges were equipped respectively above and below the canopy of the tea plant to monitor the gross rainfall and throughfall on the tea hillslope (Fig. 5). As the ratio (about 75%) of the monitored throughfall to the gross rainfall (1821 mm) during from 29 March 2016 to 28 March 2017 was in the ranges of previous studies, which with approximate leaf area index (e.g., Siles et al., 2010; Zheng et al., 2018). Thus, we think the monitored throughfall by one rain gauge was acceptable. All measurements (soil water content, matrix potentials, and meteorological data) were collected every 5 min. The daily averaged of these measured data were used in the following analyses. Soil samples at 0- to 20-cm and 20- to 40-cm depth intervals were collected using a hand auger at 39 sites on the tea hillslope and at the four automatic monitoring sites on the bamboo hillslope (Fig. 1). Three subsamples were collected for each depth interval at each site and then fully mixed. These samples were air dried, weighted, ground and sieved through a 2 mm polyethylene sieve. Particles > 2 mm were weighed to determine the RFC. Soils passed through the 2 mm polyethylene sieve were used to determine the contents of clay (< 0.002 mm), silt (0.002–0.05 mm) and sand (0.05–2 mm), as well as the contents of soil organic matter (SOM) (Liao et al., 2016). In addition, the depths to bedrock (DB) of all 43 sites (39 on the tea hillslope, 4 on the bamboo hillslope) were determined when taking soil samples using a hand auger (Fig. 1). A high-resolution (1 m) digital elevation model (DEM) was derived from a 1:1000 contour map. Terrain attributes including elevation, slope, plane curvature (PLC), profile curvature (PRC), and topographic wetness index (TWI) were determined from this DEM in ArcGIS 10.0 (ESRI, Redlands, CA). The statistics of terrain attributes and soil properties of the tea hillslope are shown in Table 1.

2. Materials and methods 2.1. Study hillslope

2.3. The predictions of SHPs This study was conducted on a hillslope (31°21′N, 119°03′E) (has an area of about 0.4 ha) in the hilly area of Taihu Lake Basin, China (Fig. 1). Green tea (Camellia sinensis (L.) O. Kuntze) is the dominated plant on the study hillslope. The adjacent hillslope is dominated by Moso bamboo (Phyllostachysedulis (Carr.) H. de Lehaie). The data of both

Scheme I: spatially varied SHPs extracted from the MultiRF based on the observed soil water retention data (MultiRFHet) Twelve groups of soil water retention data were derived in this study (Fig. 2). They were TG01-10 (10 indicated the 10 cm depth), 355

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

Fig. 1. (a) Locations of the study area and (b) the sampling and monitoring sites on the tea hillslope and adjacent bamboo hillslope, as well as (c) the topography of the study hillslope, the boundary conditions and the monitoring sites (red dot) specified in HYDRUS-3D model.

Table 1 Statistical summaries of terrain attributes and soil properties on the study hillslope. DB: depth to bedrock; PLC: plane curvature; PRC: profile curvature; RFC: rock fragment content; SOM: soil organic matter; SP: slope percent; TWI: topographic wetness index; Min: minimum; Max: maximum; SD: standard deviation. Properties

10 cm

Se (h) =

Min

Max

SD

Terrain attributes SP (%) Elevation (m) TWI PLC PRC DB (cm)

11.53 84.10 4.77 0.14 −0.12 48.74

5.43 80.71 1.95 −22.43 −20.70 33.00

19.49 87.17 7.33 34.37 12.63 77.00

3.11 1.51 1.49 8.88 5.29 9.47

Soil properties RFC (kg kg−1) Sand (%) Silt (%) Clay (%) SOM (%)

0.50 12.80 72.25 14.96 2.14

0.36 4.60 59.62 11.95 1.31

0.62 28.35 78.08 19.90 2.95

0.08 4.23 3.46 1.95 0.43

Mean

Min

Max

SD

0.51 13.54 71.62 14.84 1.74

0.32 5.44 61.96 11.91 0.94

0.78 26.13 77.75 20.64 2.87

0.11 4.43 3.60 1.87 0.43

k

1

∑ wi ⎡⎢ 1 + (α i=1

⎣

i

mi

⎤ |h|)ni ⎥ ⎦

(1)

where Se is the effective water content, h (m) is the pressure head, and θr (m3 m−3) and θs (m3 m−3) is the residual and saturated water contents, respectively. k is the number of subsystems that form the total pore-size distribution, and wi is empirical weighing coefficients for the subsystem i, and αi (m−1), ni, mi (1–1/ni) are the van Genuchten curveshape parameters of the subsystem i. In this study, the complicated pore structures of stony soils can be generally characterized by two pore systems: fine earth and rock fragments. Seven SHPs were then obtained: the residual and saturated soil water content (θr and θs), the van Genuchten curve-shape parameters of the fine earth of soil (αm and nm), the volumetric fraction of the rock fragments (wf), and the van Genuchten curve-shape parameters of the rock fragments (αf and nf). The wf was derived by transforming the RFC from gravimetric to volumetric. As the αm, nm and αf were found approximate among these twelve observation data, they were set identically as the averaged values, with αm = 3.00 m−1, nm = 1.48 and αf = 1.00 m−1, respectively, to reduce the prediction uncertainty and complexity. This disposal was reasonable and feasibility as the soil texture were homogenous on the tea hillslope, which meant the soil water retention characteristics of the fine-textured soils were approximate, to some degree. The stepwise multiple linear regression was used in Statistical Product and Service Solutions 19.0 (SPSS 19.0, SPSS Inc., IL, USA) to construct the regression model between the SHPs (θr, θs, nf) and the soil properties (RFC, sand, silt, clay, OM) and topographic attributes (slope, elevation, TWI, PLC, PRC, DB). The significance entry level for an independent variable was set as P < 0.15. The soil properties on different sampling sites were interpolated

30 cm

Mean

θ (h)−θr = θs−θr

TG02-10, TG02-30, TG03-10, TG04-10, TG04-30, BF01-10, BF02-10, BF02-30, BF03-30, BF04-10, BF04-30. The sensors at TG01-30, TG0330, BF01-30, and BF03-10 were damaged. In this study, the MultiRF proposed by Durner (1994) was adopted to fit the observed soil water retention data using the nonlinear least-squares optimization program RETC (van Genuchten et al., 1991). The MultiRF was constructed by a linear superposition of subcurves of the van Genuchten (van Genuchten, 1980) type:

356

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

Fig. 2. The observed soil water retention data and the multimodal retention function fitted (MultiRF), ROSETTA predicted (ROSETTA), and rock fragment content adjusted (ROSETTA + RFC) soil water retention curves.

(1979) algorithm was used to run the k-means clustering, and the elbow criterion (Twarakavi et al., 2010) was used to estimate an optimal number of clusters.In our study, 12 clusters were sufficient to present the spatial variance of the SHPs at both 10- and 30-cm depths, as adding one additional cluster did not add sufficient information based on the elbow criterion (Twarakavi et al., 2010). The attributes at the cluster centroid were determined as the SHPs on the corresponding cluster. In addition, the saturated hydraulic conductivities (Ks) were measured on each cluster at 10- and 30-cm depths by tension infiltrometer.

with a resolution of 1 m using ordinary kriging in ArcGIS 10.0. The spatial distributions of SHPs at 10- and 30-cm depths were obtained based on the regression models and the maps of soil properties and terrain attributes. The k-means clustering algorithm (Jain, 2010) was applied with the maps of SHPs as inputs to classify the SHPs at 10- and 30-cm depths. In the k-means clustering algorithm, cluster centers are initially chosen randomly from the set of input observation vectors (i.e. the vectors of SHPs). An iterative approach was applied by minimizing the objective function ξ(k) between the input vector and the centroid vector: k

ξ (k ) =

ni

Scheme II: Spatially uniformed SHPs derived from ROSETTA (RosHom) The ROSETTA pedotransfer function (Schaap et al., 2001) was adopted to predict the SHPs at the 39 soil sampling sites based on the contents of sand, silt and clay. The maps of ROSETTA-derived SHPs were generated by spatial interpolating the SHPs at these 39 sites using ordinary kriging in ArcGIS 10.0. The arithmetic means of these maps of SHPs were used as the spatially uniformed SHPs in this scheme.

∑ ∑ ||Xij −ui ||2 i=1 j=1

(2)

where k is the number of clusters determined a priori, ni is the number of data objects belonging to the cluster i, Xij is the vector of attributes of the j data object of the cluster i, ui is the mean value of the attributes for the data objects in the cluster i, and || || describes the Euclidean distance between two data objects. In this study, Hartigan and Wong’s 357

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

domain were treated as zero flux boundaries since the upslope was the top of the mountain and the left side had a trench opened to separate the hillslope with the surrounding area (Fig. 1). A seepage face boundary condition was assigned to the downslope boundary (near a pond) and the right side of the domain (isolated by a trench) (Fig. 1). This allowed water to leave the domain through sections of the boundary that were above the field capacity (pressure head = −0.18 m). A free drainage boundary condition was specified for the bottom (the bedrock), as the groundwater table was deep (> 8 m) and the bedrock had high water permeability (Fig. 1). Initial conditions were defined as a soil water content of 0.25 m3 m−3. The ETp was determined as crop reference evapotranspiration (ET0) multiplied by a crop coefficient of the tea (assumed to a constant equal to 1 according to Stephens and Carr, 1991, and Kigalu, 2007). The Hargreaves method (Hargreaves and Samani, 1985) was applied to calculate ET0 based on the monitored daily maximum, minimum and mean temperature. The ETp was partitioned into potential evaporation (Ep) and potential transpiration (Tp) based on Beer’s law (Ritchie, 1972):

Scheme III: Spatially varied SHPs derived from ROSETTA (RosHet) The maps of ROSETTA-derived SHPs in Scheme II were used to obtain the 12 clusters using the k-means clustering for 10- and 30-cm depths. As the SHPs predicted by ROSETTA were more spatially homogenous than those fitted by MultiRF, 12 clusters were adequate to capture the spatial variation. The attributes at the cluster centroid were determined as the SHPs on the corresponding cluster. Scheme IV: Spatial uniformed SHPs derived from ROSETTA and adjusted by RFCs (RosRFCHom) The influences of RFC were considered by adjusting the maps of ROSETTA-predicted SHPs (θr, θs and Ks) in Scheme II by empirical functions proposed by Ravina and Magier (1984) and Brakensiek et al. (1986) (denoted as “ROSETTA + RFC”).

θT = (1−wf ) × θm + wf × θf

(3)

KsT = (1−wf )×Ksm

(4)

where θT (m3 m−3) and KsT (m d−1) are the adjusted residual/saturated soil water contents and saturated hydraulic conductivity of the stony soils, respectively; wf is the volumetric fraction of the rock fragments, m3 m−3; θm (m3 m−3) and Ksm (m d−1) are the ROSETTA predicted residual/saturated soil water contents and saturated hydraulic conductivity of the matrix part, respectively; θf (m3 m−3) is the residual/ saturated soil water contents of the rock fragments (set as 0.00- and 0.25-m3 m−3, respectively).Then the maps of SHPs derived from ROSETTA and adjusted with RFC were generated for different depths on the study hillslope. The arithmetic means of these maps were used as the spatially uniformed SHPs in this scheme.

Ep = ETp × e−k × LAI

(5)

Tp = ETp−Ep

(6)

where k is the extinction coefficient, and LAI is the leaf area index. In this study, k was set to 0.55 according to White et al. (2000). The LAI was assumed to be a constant equal to 6 through the little varied monthly observed data by LAI-2200 plant canopy analyzer (LI-COR Inc., Lincoln, NE, USA). The Tp was employed for simulating the root water uptake according to the S-Shaped water stress response function (van Genuchten, 1987). The soil water potential at which root water uptake was reduced by 50% (P50) in the water stress response function was usually set from −10 m to −50 m. However, as high RFC in the soils, a lower P50 value (-8 m) was suggested by Zhu et al. (2009). The exponent in the water stress response function (p) was set to the recommended value of −3. The root length density was measured using the standard auger method (Smit et al., 2013). Total 600 days of observed data (from 7 August 2015 to 28 March 2017) were used to run the model. The first 234 days were set as spinup period, and the simulation results of the latter 365 days (from 29 March 2016 to 28 March 2017) were validated by the observed data. As no finite element node was exactly corresponded to the depth (10- or 30-cm) where soil water content was measured, the soil water storage (SWS) calculated by soil water contents at 10- and 30-cm depths at these four sites were used to validate the simulated SWS. The root mean squared error (RMSE) and Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) were used to evaluate the simulation accuracies.

Scheme V: Spatial varied SHPs derived from ROSETTA and adjusted by RFCs (RosRFCHet) The maps of SHPs derived from ROSETTA + RFC in Scheme IV were used to obtain the 12 clusters using the k-means clustering for 10- and 30-cm depths. Twelve clusters were also adequate to capture the spatial variation of the SHPs derived from ROSETTA + RFC. The attributes at the cluster centroid were determined as the SHPs on the corresponding cluster. 2.4. The simulation by HYDRUS-3D model The HYDRUS-3D model (Šimůnek et al., 2012) was used in this study. The model domain was generated by importing the x, y, z values (respectively convert from longitude, altitude and elevation) of the surface, 20-cm depths and bedrock layer on the study hillslope. The z values of the 20-cm depth and bedrock layers were derived by subtracting 20 cm and the DB from the surface DEM, respectively. Two sublayers were then obtained, the topsoil (0–20 cm) representing horizon A and the subsoil (from 20 cm to the bedrock) representing the horizon B. The average spacing between nodes in horizontal direction was set as 0.4 m, and the total finite elements mesh contained 86,850 nodes, arranged in ten mesh layers and resulting in 151,524 3-D elements in the form of triangular prisms. The topsoil consisted of the top six mesh layers whereas the remaining lower four mesh layers belonged to the subsoil and the lowest one mesh layer represented the bedrock. The spatial distributions of SHPs at 10- and 30-cm depths in these five schemes were used to represent the SHPs in the topsoil and subsoil. In the scheme of MultiRFHet, the composite hydraulic conductivity function proposed by Durner et al. (1999) was used to describe the soil hydraulic conductivities in HYDRUS-3D model. The bedrock SHPs were assumed as θr = 0.05 m3 m−3, θs = 0.30 m3 m−3, α = 1.02 m−1, n = 6.05, Ks = 0.08 m d−1 and l = 0.5 referred to Heilweil and Solomon (2004) and Katsura et al. (2006) with the similar lithology. An atmospheric boundary condition with daily throughfall and potential evapotranspiration (ETp) was imposed on the surface of the model domain (Fig. 1). The upslope boundary and the left side of the

3. Results 3.1. The fitted and predicted SHPs The MultiRF well fitted the observed soil water retention data at different depths and sites on the study tea hillslope and the adjacent bamboo hillslope (Table 2, Fig. 2). The RMSE values were generally below 0.02 m3 m−3, and the NSE values were generally above 0.90 (Table 2). Therefore, the SHPs derived from the MultiRF were treated as the observed SHPs. The fitted θr, θs, wf, and nf varied at different depths and sites (Table 2), implying that large spatial variances of SHPs on this study hillslope. For example, the largest wf was found at TG02-30 (0.60 m3 m−3), while the smallest wf was found at BF02-10 (0.18 m3 m−3); the largest nf value was 22.29 at BF02-10, while the smallest nf was 1.48 and at BF04-30. The SHPs derived from ROSETTA and ROSETTA + RFC had relatively poor performances in fitting the observed soil water retention data (Fig. 2). The RMSE values were generally above 0.04 m3 m−3, and the NSE values were generally below zero. However, SHPs derived from 358

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

Table 2 Fitted SHPs by Durner’s (1994) multimodal retention function for the twelve observed soil water retention datasets and the fitting accuracies (RMSE and NSE). The subscripts of m and f refer to the fine earth and rock fragments of the stony soils, respectively. As the αm, nm and αf were found approximate among these twelve datasets, they were set as the identical, with αm = 3.00 m−1, nm = 1.48 and αf = 1.00 m−1, respectively. Site-depth

θr (m3 m−3)

θs (m3 m−3)

αm(m−1)

nm

wf (m3 m−3)

αf (m−1)

nf

RMSE (m3 m−3)

NSE

TG01-10 TG02-10 TG02-30 TG03-10 TG04-10 TG04-30 BF01-10 BF02-10 BF02-30 BF03-30 BF04-10 BF04-30

0.09 0.11 0.12 0.16 0.09 0.12 0.10 0.13 0.08 0.22 0.17 0.18

0.32 0.37 0.41 0.52 0.38 0.44 0.31 0.54 0.55 0.59 0.60 0.43

3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00

1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48

0.43 0.47 0.60 0.23 0.31 0.27 0.33 0.18 0.28 0.29 0.20 0.31

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

10.95 5.29 1.78 1.62 3.09 4.71 6.47 22.29 13.81 2.96 2.18 1.48

0.018 0.016 0.015 0.017 0.029 0.015 0.011 0.015 0.022 0.010 0.016 0.006

0.76 0.87 0.91 0.94 0.80 0.93 0.93 0.93 0.92 0.98 0.97 0.99

ROSETTA + RFC were generally better than those derived from ROSETTA (Fig. 2). The θr values derived from ROSETTA and ROSETTA + RFC were similar (ranging from 0.06- to 0.08-m3 m−3) and generally smaller than these derived from the MultiRF (ranging from 0.08- to 0.22-m3 m−3). However, the θs predicted by ROSETTA + RFC were closer to those derived from the MultiRF than those predicted by ROSETTA (Fig. 2). Smaller curvatures were found for the soil water retention curves derived from ROSETTA and ROSETTA + RFC than those derived from the MultiRF (Fig. 2). The mean α was 0.53 m−1 for both the ROSETTA and ROSETTA + RFC and ranged from 0.44- to 0.69-m−1, and the mean n was 1.66 and ranged from 1.62 to 1.70.

depths were observed for RosHet (Table 5). The θr, θs and n were about 0.07 m3 m−3, 0.45 m3 m−3 and 1.67, respectively, in different clusters and at these two depths. The mean α was 0.52 m−1 and ranged from 0.46- to 0.61- m−1 at 10-cm depth and was 0.50 m−1 and ranged from 0.45- to 0.59- m−1 at 30-cm depth. At both depths, the mean Ks was 0.21 m d−1 and ranged from 0.18- to 0.24- m d−1. However, for RosRFCHet, obvious differences of some SHPs were observed among different clusters and at different depths (Table 5). The θs ranged from 0.33- to 0.40-m3 m−3 and from 0.30- to 0.40-m3 m−3 at 10- and 30-cm depths, respectively, while the Ks ranged from 0.10- to 0.20-m d−1 and from 0.13- to 0.15-m d−1.

3.2. The spatial distributions of SHPs

3.3. The simulated SWS

Local PTFs were developed by constructing the regression models between the SHPs derived from the MultiRF (θr, θs and nf) and the soil properties and topographic attributes (Table 3). These regression equations well predicted the θr, θs and nf with both the determination coefficient (R2) and NSE above 0.68. The independent variables for predicting the θr were silt, clay and elevation, for predicting the θs were sand, clay, SOM, RFC and DB, and for predicting the nf were TWI, clay and SOM. Based on these regression models, maps of SHPs derived from the MultiRF were generated (Fig. 3). Great spatial variation can be observed for these SHPs. The spatial mean θr values were both 0.11 m3 m−3 at 10- and 30-cm depths, while the coefficients of variations (CVs) of the θr were 28.1% and 18.6% at these two depths, respectively. The spatial means of θs was 0.43 m3 m−3and 0.42 m3 m−3 at 10- and 30-cm depths, respectively, while the CVs for these two depths were 15.3% and 17.5%. The spatial means of wf at 10- and 30-cm depths were 0.34 m3 m−3 and 0.36 m3 m−3, respectively, while the CVs for these two depths were 18.15% and 24.77%. The nf value was generally smaller than 8.00. The CVs of nf were 46.0% and 32.1% at 10and 30-cm depths, respectively. The SHPs of MultiRFHet in different clusters were shown in Table 4 and the spatial distributions of the clusters were shown in Fig. 4a and b. The 12 clustered SHPs for RosHet and RosRFCHet and their spatial distributions at different depths were shown in Table 5 and Fig. 4c–f. Negligible distinctions of the SHPs among clusters and at different

Performances of HYDRUS-3D based on SHPs derived from these five schemes varied on this study hillslope (Fig. 5). The simulated SWSs by MultiRFHet were close to the observations at these four sites, with the NSE ≥ 0.58, and the RMSE < 16.50 mm. Both of the RosHom and RosHet generally overestimated the SWSs, especially during the wet period. The NSE values of these two schemes at sites TG01, TG02 and TG04 were < 0. The SWSs at sites TG01 and TG04 were overestimated by the RosRFCHom and RosRFCHet, while were well simulated at sites TG02 and TG03. Low NSE values (< 0.30) were found for both RosRFCHom and RosRFCHet at sites TG01 and TG04. At site TG03, the NSE values of these two schemes were > 0.80, while at TG02, the NSE value was > 0.60 for RosRFCHet but < 0 for RosRFCHom. 3.4. The accumulative boundary flux Slight differences were also found among the accumulative boundary fluxes derived from different schemes (Fig. 6). For all five schemes, the accumulative free drainage fluxes and evapotranspiration rates were two major sink terms of the water balance, which individually accounted for about 50% of the accumulative precipitation (1821.00 mm). The smallest accumulative drainage flux (896.00 mm) and the largest accumulative evapotranspiration (860.00 mm) were found for the RosHom and RosHet, while the accumulative drainage flux (around 928.00 mm) and evapotranspiration (around 840.00 mm) were approximate for RosRFCHom, RosRFCHet and MultiRFHet.

Table 3 The best regression models that predicted the Durner’s multimodal retention function SHPs of θr, θs, nf by soil properties and topographic attributes. DB: depth to bedrock; RFC: rock fragment content; SOM: soil organic matter content; TWI: topographic wetness index. parameter

Regression equations

R2

NSE

θr (m3 m−3) θs (m3 m−3) nf

0.014 + 0.006*Silt + 0.011*Clay-0.006*Elevation 0.374 + 0.006*Clay + 0.026*SOM-0.009*Sand + 0.003*DB-0.225*RFC 22.946–1.997*TWI-1.025*Clay + 5.026*SOM

0.722 0.724 0.719

0.721 0.686 0.719

359

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

Fig. 3. The spatial distributions of soil hydraulic parameters (θr, θs, wf, nf) derived from Durner’s multimodal retention function fitted soil water retention curves at 10- and 30-m depths, respectively.

4. Discussion

observed soil water retention data, performances in simulating the SWS were obviously better than the other schemes (Fig. 5). However, since the soil at TG03 was fine-textured (wf < 0.25 m3 m−3, sand < 9.0%), the NSE of MultiRFHet in simulating SWS was lower than those of RosRFCHom and RosRFCHet (Fig. 5). This indicated that the MultiRF may overestimate the influence of RFC for the fine-textured soils. When the ROSETTA predicted SHPs were adjusted by RFC, improved accuracies were observed in fitting the soil water retention data and simulating the SWS variations (Figs. 2 and 5). For all the depths and sites, RMSE values in fitting the soil water retention data decreased and NSE values increased when adopting the SHPs derived from

4.1. The influences of RFC The SHPs were well extracted by the MultiRF (Fig. 2) and the observed SWSs were well simulated (Fig. 5) by adopting MultiRFHet in HYDRUS-3D model. As comprehensively considering the pore systems of stony soils (the fine earth and the rock fragments), the MultiRF can better characterize the SHPs in low suction region than the traditional single-porosity models (Durner, 1994; Šimůnek et al., 2003; Dettmann et al., 2014) (Fig. 2). As SHPs in MultiRFHet were calibrated by the

Table 4 The Durner’s multimodal retention function derived SHPs on the 12 clusters at 10- and 30-cm depths of the tea hillslope. The subscripts of m and f refer to the fine earth and rock fragments of the stony soils, respectively. Depth (cm)

cluster

θr (m3 m−3)

θs (m3 m−3)

αm (m−1)

nm

Ks (m d−1)

l

wf (m3 m−3)

αf (m−1)

nf

1

1 2 3 4 5 6 7 8 9 10 11 12

0.09 0.13 0.13 0.11 0.11 0.07 0.04 0.08 0.03 0.10 0.09 0.11

0.34 0.53 0.46 0.46 0.49 0.36 0.32 0.39 0.32 0.43 0.40 0.41

3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00

1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48

0.12 0.11 0.10 0.13 0.07 0.12 0.10 0.11 0.11 0.11 0.17 0.09

0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50

0.30 0.27 0.35 0.29 0.26 0.28 0.42 0.32 0.30 0.39 0.40 0.43

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

4.29 9.44 6.95 12.98 1.29 16.42 4.00 8.67 10.69 8.95 15.10 2.97

30

1 2 3 4 5 6 7 8 9 10 11 12

0.11 0.13 0.07 0.12 0.13 0.12 0.11 0.09 0.08 0.12 0.08 0.08

0.50 0.43 0.31 0.46 0.56 0.34 0.39 0.40 0.42 0.42 0.35 0.34

3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00 3.00

1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48 1.48

0.06 0.11 0.12 0.10 0.11 0.11 0.06 0.10 0.10 0.10 0.16 0.13

0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50

0.24 0.28 0.41 0.38 0.28 0.32 0.56 0.45 0.29 0.34 0.41 0.44

1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00

3.74 3.47 3.00 5.89 7.30 6.82 2.35 3.88 7.87 8.95 11.00 6.32

360

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

Fig. 4. The spatial patterns of the clusters of SHPs in schemes of (a, b) MultiRFHet, (c, d) RosHet and (e, f) RosRFCHet at topsoil and subsoil in the HYDRUS-3D model.

Fig. 5. Observed and simulated soil water storage (SWS) and the corresponding simulation accuracies (NSE and RMSE) at sites (a) TG01, (b) TG02, (c) TG03 and (d) TG04 in schemes of RosHom, RosHet, RosRFCHom, RosRFCHet and MultiRFHet, and (e) time series of the gross rainfall and throughfall. As sensors were damaged, no observed SWSs were derived during certain period.

361

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

Table 5 The spatial uniformed and clustered SHPs predicted by ROSETTA (RosHom and RosHet) and adjusted with rock fragment contents (RosRFCHom and RosRFCHet) at 10- and 30-cm depths on the tea hillslope. Scheme

Cluster

θr (m3 m−3)

θs (m3 m−3)

Α (m−1)

Ks (m d−1)

n

l

10

30

10

30

10

30

10

30

10

30

10

RosHom RosHet

/ 1 2 3 4 5 6 7 8 9 10 11 12

0.07 0.07 0.07 0.07 0.06 0.06 0.07 0.07 0.06 0.07 0.07 0.07 0.07

0.07 0.07 0.07 0.06 0.07 0.06 0.07 0.07 0.06 0.06 0.06 0.06 0.07

0.45 0.46 0.46 0.45 0.46 0.45 0.46 0.44 0.44 0.45 0.45 0.45 0.45

0.45 0.45 0.45 0.45 0.46 0.45 0.44 0.45 0.44 0.45 0.44 0.44 0.46

0.51 0.61 0.57 0.49 0.52 0.49 0.56 0.46 0.46 0.49 0.53 0.52 0.52

0.50 0.53 0.54 0.51 0.59 0.51 0.47 0.50 0.45 0.49 0.47 0.47 0.55

1.67 1.64 1.65 1.67 1.67 1.67 1.66 1.68 1.69 1.66 1.66 1.66 1.67

1.67 1.66 1.65 1.68 1.64 1.67 1.67 1.67 1.69 1.68 1.68 1.68 1.65

0.21 0.18 0.18 0.22 0.23 0.22 0.20 0.24 0.23 0.19 0.20 0.19 0.22

0.21 0.21 0.20 0.23 0.18 0.21 0.20 0.20 0.24 0.22 0.23 0.22 0.19

0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50

RosRFCHom RosRFCHet

/ 1 2 3 4 5 6 7 8 9 10 11 12

0.06 0.06 0.06 0.06 0.06 0.05 0.06 0.06 0.06 0.06 0.06 0.06 0.06

0.06 0.06 0.06 0.05 0.06 0.05 0.06 0.06 0.06 0.06 0.06 0.06 0.06

0.38 0.38 0.38 0.39 0.35 0.33 0.40 0.38 0.39 0.37 0.40 0.33 0.38

0.38 0.36 0.40 0.33 0.39 0.30 0.33 0.39 0.36 0.39 0.37 0.39 0.38

0.51 0.51 0.46 0.56 0.47 0.53 0.60 0.46 0.53 0.51 0.57 0.49 0.49

0.50 0.50 0.59 0.49 0.49 0.46 0.50 0.48 0.45 0.54 0.52 0.55 0.51

1.67 1.68 1.69 1.65 1.68 1.67 1.64 1.68 1.66 1.66 1.66 1.67 1.67

1.67 1.67 1.64 1.68 1.67 1.68 1.67 1.68 1.69 1.66 1.70 1.65 1.67

0.14 0.15 0.17 0.11 0.15 0.13 0.11 0.20 0.12 0.10 0.14 0.14 0.16

0.13 0.13 0.13 0.14 0.13 0.13 0.13 0.14 0.14 0.14 0.15 0.13 0.14

0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50 0.50

30

especially in hilly areas (Baetens et al., 2009). Even so, imperfect performances were still observed when using the ROSETTA predicted and RFC adjusted SHPs as inputs to fit the soil water retention data and simulate the SWS variations (Figs. 2 and 5). For example, the NSE values in simulating the SWS variations were below 0.25 at TG01 and TG04 for RosRFCHom and RosRFCHet (Fig. 5). This indicated that the empirical functions proposed by Ravina and Magier (1984) and Brakensiek et al. (1986) may not be effective in all cases, which has also been reported by Novák et al. (2011) and Coppola et al. (2013). When using the ROSETTA + RFC to fit the soil water retention curves, influences of RFC on the curve parameters α and n were not considered. This can be a reason for the imperfect performances of ROSETTA + RFC (Baetens et al., 2009; Ma et al., 2010).

ROSETTA + RFC. The changed RMSE and NSE values were significantly correlated with wf (P < 0.05), which indicated the improved performance by considering RFC were more obvious for the coarse textured soils. This is due to that ROSETTA was constructed without considering the rock fragments, thus SHPs of stony soils were overestimated when predicted by the ROSETTA. After adjusted by RFC based on the empirical equations, the predicted SHPs (θr, θs and Ks) decreased and thus the accuracies in fitting the soil water retention data were improved. Additionally, comparing the RosHom with RosRFCHom and RosHet with RosRFCHet, the NSE values all improved and the RMSE values all reduced (Fig. 5). However, for TG03 with relative finetextured soils, the improvements of prediction accuracies were less obvious. Coppola et al. (2013) also demonstrated that better simulations of the average evolution of water content were achieved when the PTFs predicted SHPs were scaled for stoniness. Therefore, it is necessary to incorporate the RFC as a predictor for the construction of PTFs, and the measurements of RFC should be included in the soil survey,

4.2. The influence of spatial variation of SHPs Adopting the spatial varied SHPs can improve the simulation of the Fig. 6. (a) The accumulative boundary fluxes and the ratios to the accumulative throughfall in schemes of RosHom, RosHet, RosRFCHom, RosRFCHet, and MultiRFHet in the HYDRUS-3D model, and (b) time series of the throughfall on the tea hillslope. Cum_drainage: the accumulative free drainage and seepage face fluxes; Cum_evapotranspiration: the accumulative evapotranspiration flux; Cum_throughfall: the accumulative throughfall amounts; Cum_δSWS: the accumulative daily changed soil water storage.

362

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

evapotranspiration but high spatial variations of SHPs, the changed SWS may account for large portion in the water balance and thus spatial variation and RFC may show strong effects on subsurface flow simulation. To improve the simulation of subsurface flow, accurate inputs of precipitation and evapotranspiration can be more important. Precipitation is usually monitored by raingauge. For the evapotranspiration, numerous studies have reported its significant influences on the hydrological modelling (e.g., Ivanov et al., 2010; Vereecken et al., 2015; Feng et al., 2017). To acquire the accurate evapotranspiration parameters, we should either equip with more advanced monitoring systems (e.g., eddy-covariance flux tower and large aperture scintillometer), or use more advanced algorithms/data to calculate/estimate the evapotranspiration rate (Wang and Dickinson, 2012).

SWS only when the influence of RFC was considered (RosRFCHom and RosRFCHet) (Fig. 5). Since the spatial variations of SHPs predicted by ROSETTA was low (Table 5), this improvement was not presented between RosHet and RosHom (Fig. 5). However, comparing the simulation accuracies by RosRFCHet with those by RosRFCHom, obvious improvements were observed at sites TG01 and TG02 (reduce the RMSE by 47.9% and 48.6%, respectively). This can be attributed to distinct differences between the SHPs derived from RosRFCHet and RosRFCHom at these two sites and the SHPs derived from RosRFCHet may be closer to the actual SHPs. For TG03, as the SHPs derived from RosRFCHom and RosRFCHet may be both close to the real values, the simulation accuracies were both high. For TG04, as the SHPs derived from RosRFCHom and RosRFCHet both deviated from the actual SHPs, the simulation accuracies were both low. Numerous previous studies have emphasized the influence of the spatial variation of SHPs on simulated hydrological states and fluxes (e.g., Herbst et al., 2006; Fiori and Russo, 2007; Fang et al., 2016; Baroni et al., 2017). Many approaches have been proposed to precisely scale and characterize the spatial variation of SHPs (e.g., Vereecken et al., 2007; Fang et al., 2016). However, only a few studies were devoted to improving the simulation accuracy by exactly characterizing the actual spatial variation of SHPs based on the existing scaling approaches. For example, Jin et al. (2015) found the simulation by Soil and Water Assessment Tool with spatially clustered SHPs had more accurate performance than using the average soil properties; Fang et al. (2016) employed the information entropy concept for the scale dependent parameterization of soil hydraulic conductivity and found it can improve model performance both in terms of soil water content and runoff simulation. The optimal performances of MultiRFHet (Fig. 5) demonstrated that the k-means clustering was appropriate to cluster and simplify the SHPs’ spatial distribution. This could significantly reduce the simulation times that were needed to run the model with finer resolution of SHPs. The k-means clustering has been also used in previous studies. For example, Twarakavi et al. (2010) used k-means clustering algorithm to determine the optimal number of soil hydraulic classes; Zhang et al. (2016) proposed a combination of k-means clustering of kriged soil texture data, PTFs, and model inversion to parameterize the vadose zone hydraulic properties.

5. Conclusions In this study, whether considering RFC and spatial variation of SHPs were necessary in simulating soil water dynamics were investigated on a representative stony-soil hillslope. Results shown adopting the MultiRFHet which was calibrated by the observed soil water retention data and considered the spatial variation of SHPs achieved the best accuracies in simulating the SWS for this high RFC soil. However, RosRFCHet which considered the influence of RFC and the spatial variation of SHPs could be the most proper choice in SWS simulation since it did not need the observed soil water retention data and calibration. The RosHom, RosHet and RosRFCHom which did not both consider the influence of RFC on SHPs and the spatial variation of SHPs yield unacceptable simulations. This demonstrated that for stony soils, RFC impacts and spatial variation of SHPs should be both considered in simulating soil water dynamics. However, this consideration had limited effects in subsurface flow simulation. Findings of this study will facilitate the hydrological modeling in stony soils and supplement the knowledge of the influences of RFC. Meanwhile, the results of this study will also be beneficial for the agricultural development and management in the hilly areas, in which large amount of rock fragments were always contained in the soils. Acknowledgments This study was financially supported by the National Natural Science Foundation of China (41622102), Key Research Plans of Frontier Sciences, Chinese Academy of Sciences (QYZDB-SSW-DQC038) and Jiangsu Natural Science Foundation (BK20177777).

4.3. Determining proper approaches When the observed soil water retention data were available, the MultiRFHet which calibrated by the observed data and considered the spatial variation of SHPs yielded the best simulation of SWS variations (Fig. 5). When the observed soil water retention data were unavailable, the RosRFCHet which both considered the influence of RFC on SHPs and the spatial variation of SHPs also had acceptable accuracies in SWS simulation and could be an alternative (Fig. 5). Relative to the MultiRFHet, the RosRFCHet is more time-efficient and economical as it is not necessary to be calibrated by observed soil water retention data. As the influence of RFC on SHPs and the spatial variation of SHPs were neither considered in RosHom, RosHet and RosRFCHom, these schemes were incapable to simulate the SWS variations on the stony-soil hillslope (NSE values generally < 0) (Fig. 5). Considering the spatial variation of SHPs and effects of RFC on SHPs had limited influences on the subsurface flow simulation, although they had critical influences on SWS simulation. This was consistent with the study by Coppola et al. (2013). The reasons for comparatively small differences of the boundary fluxes among different schemes maybe: (i) the SWS differences of different schemes only accounted for a small portion of the soil water balance comparing with the precipitation and evapotranspiration on this study hillslope, and (ii) the spatial variations of SHPs were not strong enough to show obvious effect on the subsurface flow simulation since the study hillslope was comparatively small (0.4 ha). Therefore, in other areas with low precipitation or

References Arora, B., Mohanty, B.P., Mcguire, J.T., 2011. Inverse estimation of parameters for multidomain flow models in soil columns with different macropore densities. Water Resour. Res. 47 (4) 2010WR009451 10.1029/2010WR009451. Assouline, S., Or, D., 2013. Conceptual and parametric representation of soil hydraulic properties: a review. Vadose Zone J. 12 (4) https://doi.org/10.2136/ vzj2013.07.0121. Baetens, J.M., Verbist, K., Cornelis, W.M., Gabriels, D., Soto, G., 2009. On the influence of coarse fragments on soil water retention. Water Resour. Res. 45 (7). https://doi.org/ 10.1029/2008wr007402. Baroni, G., Zink, M., Kumar, R., Samaniego, L., Attinger, S., 2017. Effects of uncertainty in soil properties on simulated hydrological states and fluxes at different spatio-temporal scales. Hydrol. Earth Syst. Sci. 21 (5), 2301–2320. https://doi.org/10.5194/ hess-21-2301-2017. Bouwer, H., Rice, R., 1984. Hydraulic properties of stony vadose zones. Groundwater 22 (6), 696–705. Brakensiek, D., Rawls, W., Stephenson, G., 1986. Determining the saturated hydraulic conductivity of a soil containing rock fragments. Soil Sci. Soc. Am. J. 50 (3), 834–835. Chirico, G.B., Medina, H., Romano, N., 2007. Uncertainty in predicting soil hydraulic properties at the hillslope scale with indirect methods. J. Hydrol. 334 (3–4), 405–422. https://doi.org/10.1016/j.jhydrol.2006.10.024. Coppola, A., Dragonetti, G., Comegna, A., Lamaddalena, N., Caushi, B., Haikal, M.A., Basile, A., 2013. Measuring and modeling water content in stony soils. Soil Till. Res. 128, 9–22. https://doi.org/10.1016/j.still.2012.10.006. Cousin, I., Nicoullaud, B., Coutadeur, C., 2003. Influence of rock fragments on the water

363

Journal of Hydrology 565 (2018) 354–364

X. Lai et al.

stony-soil water retention. Agr. Forest Meteorol. 244–245, 1–8. Poesen, J., Lavee, H., 1994. Rock fragments in top soils: significance and processes. Catena 23 (1), 1–28. Ravina, I., Magier, J., 1984. Hydraulic conductivity and water retention of clay soils containing coarse fragments. Soil Sci. Soc. Am. J. 48 (4), 736–740. Ritchie, J.T., 1972. Model for predicting evaporation from a row crop with incomplete cover. Water Resour. Res. 8 (5), 1204–1213. Schaap, M.G., Leij, F.J., Van Genuchten, M.T., 2001. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. 251 (3), 163–176. Sciuto, G., Diekkrüger, B., 2010. Influence of soil heterogeneity and spatial discretization on catchment water balance modeling. Vadose Zone J. 9 (4). https://doi.org/10. 2136/vzj2009.0166. Siles, P., Vaast, P., Dreyer, E., Harmand, J.M., 2010. Rainfall partitioning into throughfall, stemflow and interception loss in a coffee (Coffea arabica L.) monoculture compared to an agroforestry system with Inga densiflora. J. Hydrol. 395 (1–2), 39–48. Šimůnek, J., Jarvis, N.J., Van Genuchten, M.T., Gärdenäs, A., 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. 272 (1), 14–35. Šimůnek, J., Van Genuchten, M.T., Šejna, M., 2012. The HYDRUS software package for simulating two-and three-dimensional movement of water, heat, and multiple solutes in variably-saturated media. Technical manual, version 2.0, 1: 260. Smit, A.L., Bengough, A.G., Engels, C., van Noordwijk, M., Pellerin, S., van de Geijn, S., 2013. Root methods: a handbook. Springer Science & Business Media. Stephens, W., Carr, M., 1991. Responses of tea (Camellia sinensis) to irrigation and fertilizer. II. Water use. Exp. Agr. 27 (2), 193–210. Tetegan, M., Nicoullaud, B., Baize, D., Bouthier, A., Cousin, I., 2011. The contribution of rock fragments to the available water content of stony soils: proposition of new pedotransfer functions. Geoderma 165 (1), 40–49. Twarakavi, N.K.C., Šimunek, J., Schaap, M.G., 2010. Can texture-based classification optimally classify soils with respect to soil hydraulics? Water Resour. Res. 46 (1). https://doi.org/10.1029/2009wr007939. Van Genuchten, M.T., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44 (5), 892–898. Van Genuchten, M.T., 1987. A numerical model for water and solute movement in and below the root zone. United States Department of Agriculture Agricultural Research Service US Salinity Laboratory. Van Genuchten, M.T., Leij, F., Yates, S., 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils. Vereecken, H., Huisman, J.A., Hendricks Franssen, H.J., Brüggemann, N., Bogena, H.R., Kollet, S., Javaux, M., van der Kruk, J., Vanderborght, J., 2015. Soil hydrology: recent methodological advances, challenges, and perspectives. Water Resour. Res. 51 (4), 2616–2633. Vereecken, H., Kasteel, R., Vanderborght, J., Harter, T., 2007. Upscaling hydraulic properties and soil water flow processes in heterogeneous soils: a review. Vadose Zone J. 6 (1), 1–28. https://doi.org/10.2136/vzj2006.0055. Vereecken, H., Maes, J., Feyen, J., Darius, P., 1989. Estimating the soil moisture retention characteristic from texture, bulk density, and carbon content. Soil Sci. 148 (6), 389–403. Vereecken, H., Weynants, M., Javaux, M., Pachepsky, Y., Schaap, M.G., van Genuchten, M.T., 2010. Using pedotransfer functions to estimate the van genuchten-mualem soil hydraulic properties: a review. Vadose Zone J. 9 (4), 795. https://doi.org/10.2136/ vzj2010.0045. Wang, K., Dickinson, R.E., 2012. A review of global terrestrial evapotranspiration: observation, modeling, climatology, and climatic variability. Rev. Geophys. 50 (2). White, M.A., Thornton, P.E., Running, S.W., Nemani, R.R., 2000. Parameterization and sensitivity analysis of the BIOME–BGC terrestrial ecosystem model: net primary production controls. Earth Interact. 4 (3), 1–85. Wösten, J., Lilly, A., Nemes, A., Le Bas, C., 1999. Development and use of a database of hydraulic properties of European soils. Geoderma 90 (3), 169–185. Wösten, J., Pachepsky, Y.A., Rawls, W., 2001. Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics. J. Hydrol. 251 (3), 123–150. Zhang, Y., Schaap, M.G., Guadagnini, A., Neuman, S.P., 2016. Inverse modeling of unsaturated flow using clusters of soil texture and pedotransfer functions. Water Resour. Res. 52 (10), 7631–7644. Zheng, J., Fan, J.L., Zhang, F.C., Yan, S.C., Xiang, Y.Z., 2018. Rainfall partitioning into throughfall, stemflow and interception loss by maize canopy on the semi-arid Loess Plateau of China. Agr. Water Manage. 195, 25–36. https://doi.org/10.1016/j.agwat. 2017.09.013. Zhu, Y., Ren, L., Sakggs, T.H., Lü, H., Yu, Z., Wu, Y., Fang, X., 2009. Simulation of Populus euphratica root uptake of groundwater in an arid woodland of the Ejina Basin. China. Hydrol. Process. 23 (17), 2460–2469.

retention and water percolation in a calcareous soil. Catena 53 (2), 97–114. https:// doi.org/10.1016/s0341-8162(03)00037-7. Dettmann, U., Bechtold, M., Frahm, E., Tiemeyer, B., 2014. On the applicability of unimodal and bimodal van Genuchten-Mualem based models to peat and other organic soils under evaporation conditions. J. Hydrol. 515, 103–115. Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res. 30 (2), 211–223. Durner, W., Priesack, E., Vogel, H.J., Zurmühl, T., 1999. Determination of parameters for flexible hydraulic functions by inverse modeling. In: van Genuchten, M.Th., Leij, F.J., Wu, L. (Eds.), Characterization and Measurement of the Hydraulic Properties of Unsaturated Porous Media. University of California, Riverside, CA, pp. 817–829. Fang, Z., Bogena, H., Kollet, S., Vereecken, H., 2016. Scale dependent parameterization of soil hydraulic conductivity in 3D simulation of hydrological processes in a forested headwater catchment. J. Hydrol. 536 (Supplement C), 365–375. https://doi.org/10. 1016/j.jhydrol.2016.03.020. Feng, H., Zou, B., Luo, J., 2017. Coverage-dependent amplifiers of vegetation change on global water cycle dynamics. J. Hydrol. 550, 220–229. https://doi.org/10.1016/j. jhydrol.2017.04.056. Fiori, A., Russo, D., 2007. Numerical analyses of subsurface flow in a steep hillslope under rainfall: The role of the spatial heterogeneity of the formation hydraulic properties. Water Resour. Res. 43 (7). https://doi.org/10.1029/2006wr005365. Fu, C., James, A.L., Yao, H., 2015. Investigations of uncertainty in SWAT hydrologic simulations: a case study of a Canadian Shield catchment. Hydrol. Process. 29 (18), 4000–4017. https://doi.org/10.1002/hyp.10477. Hargreaves, G.H., Samani, Z.A., 1985. Reference crop evapotranspiration from temperature. Appl. Eng. Agric. 1 (2), 96–99. Hartigan, J.A., Wong, M.A., 1979. Algorithm AS 136: A k-means clustering algorithm. Appl. Stat. 28 (1), 100–108. Heilweil, V.M., Solomon, D.K., 2004. Millimeter‐to Kilometer‐Scale Variations in Vadose – Zone Bedrock Solutes: Implications for Estimating Recharge in Arid Settings. Groundwater Recharge in a Desert Environment: The Southwestern United States: 49–67. https://doi.org/10.1029/009WSA04. Herbst, M., Diekkrüger, B., Vanderborght, J., 2006. Numerical experiments on the sensitivity of runoff generation to the spatial variation of soil hydraulic properties. J. Hydrol. 326 (1–4), 43–58. https://doi.org/10.1016/j.jhydrol.2005.10.036. Hlaváčiková, H., Novák, V., Ŝimůnek, J., 2016. The effects of rock fragment shapes and positions on modeled hydraulic conductivities of stony soils. Geoderma 281, 39–48. Ivanov, V.Y., Fatichi, S., Jenerette, D., Espeleta, J., Troch, P.A., Huxman, T.E., 2010. Hysteresis of soil moisture spatial heterogeneity and the “homogenizing” effect of vegetation. Water Resour. Res. 46 (9). Jain, A.K., 2010. Data clustering: 50 years beyond K-means. Pattern Recogn. Lett. 31 (8), 651–666. https://doi.org/10.1016/j.patrec.2009.09.011. Jin, X., Zhang, L.H., Gu, J., Zhao, C., Tian, J., He, C.S., 2015. Modelling the impacts of spatial heterogeneity in soil hydraulic properties on hydrological process in the upper reach of the Heihe River in the Qilian Mountains, Northwest China. Hydrol. Process. 29 (15), 3318–3327. https://doi.org/10.1002/hyp.10437. Katsura, S.Y., Kosugi, K.I., Yamamoto, N., Mizuyama, T., 2006. Saturated and unsaturated hydraulic conductivities and water retention characteristics of weathered granitic bedrock. Vadose Zone J. 5 (1), 35–47. Kigalu, J.M., 2007. Effects of planting density on the productivity and water use of tea (Camellia sinensis L.) clones: I. Measurement of water use in young tea using sap flow meters with a stem heat balance method. Agr. Water Manage. 90 (3), 224–232. Liao, K.H., Lv, L.G., Yang, G.S., Zhu, Q., 2016. Sensitivity of simulated hillslope subsurface flow to rainfall patterns, soil texture and land use. Soil Use Manage. 32, 422–432. Liao, K.H., Xua, F., Zheng, J.S., Zhu, Q., Yang, G.S., 2014. Using different multimodel ensemble approaches to simulate soil moisture in a forest site with six traditional pedotransfer functions. Environ. Modell. Softw. 57, 27–32. https://doi.org/10.1016/ j.envsoft.2014.03.016. Ma, D., Shao, M., 2008. Simulating infiltration into stony soils with a dual-porosity model. Eur. J. Soil Sci. 59 (5), 950–959. https://doi.org/10.1111/j.1365-2389.2008. 01055.x. Ma, D., Shao, M., Zhang, J., Wang, Q., 2010. Validation of an analytical method for determining soil hydraulic properties of stony soils using experimental data. Geoderma 159 (3–4), 262–269. https://doi.org/10.1016/j.geoderma.2010.08.001. Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12 (3), 513–522. Nash, J.E., Sutcliffe, J.V., 1970. River flow forecasting through conceptual models part I—A discussion of principles. J. Hydrol. 10 (3), 282–290. Novák, V., Kňava, K., Šimůnek, J., 2011. Determining the influence of stones on hydraulic conductivity of saturated soils using numerical method. Geoderma 161 (3–4), 177–181. https://doi.org/10.1016/j.geoderma.2010.12.016. Parajuli, K., Sadeghi, M., Jones, S.B., 2017. A binary mixing model for characterizing

364