Temporal Variation in Preferential Water Flow during Natural Vegetation Restoration on Abandoned Farmland in the Loess Plateau of China

Abstract

The restoration of vegetation in abandoned farmlands is an effective approach to control soil erosion on the Chinese Loess Plateau. However, few studies have investigated the effect of natural restoration age on the infiltration patterns and preferential flow in soil layers. This study examined the effect of the temporal variations in the degree of preferential flow and their contribution on the total infiltration in abandoned farmlands restored with different vegetation communities. In total, four sites were investigated, namely site AF3, AF12, AF18 and AF25 of restoration age of 3, 12, 18 and 25 years, respectively. A traditional dye tracer experiment was adopted to visualize the water flow behavior in the four sites. A semi-variance analysis was conducted to interpret the spatial variability of the soil water induced by infiltration. The results showed that both the soil properties and the root systems improved with the increased natural restoration ages combined with the effect of the plant species. The degree of preferential flow from the lowest to the highest occurred in AF3, AF12, AF18 and AF25, which was mainly affected by the development of root systems. Also, the abundant lateral root and steeper slope may have caused the presence of lateral flow in AF18 and AF25. Meanwhile, the spatial variability of the soil water in the infiltration patterns through vertical soil profiles increased from AF3 to AF25. Compared to AF3, the infiltration volume through preferential pathways in AF12, AF18 and AF25 increased by 1.73 times, 1.94 times and 4.09 times, respectively. Therefore, the preferential flow was improved in abandoned farmland in the presence of natural vegetation restoration that can increase water storage in the deep soil layer and reduce surface runoff and erosion.

1. Introduction

The Chinese Loess Plateau has suffered from severe soil erosion for several decades. The average annual erosion rates generally range between 5000 and 10,000 t/km²/a because of the intensive human disturbance and the presence of a thick erodible soil [1]. The impact of human activities in the plateau can be attributed to continuous and widespread stress, such as large-scale monocultures (wheat and maize) and over-grazing [2]. To control soil erosion and ecosystem degradation, the Chinese government initiated the “Grain-for-Green” project across the country in 1998. As important management of the project, vast areas of farmland in the Loess Plateau have been abandoned and converted into grasslands, shrublands or forests with natural vegetation [3]. Consequently, the soil loss of the plateau decreased by 57% and the vegetation cover increased by 13 times since 1998 [4,5]. The dense vegetation restoration activities and the increasing vegetation cover require a greater amount of water to survive [6]. However, water is usually scarce in the Loess Plateau, where the precipitation is low, the evapotranspiration is high, and the groundwater is deep [7]. Therefore, understanding the responses of the hydrological process to vegetation restoration in abandoned farmlands is essential for estimating their sustainable development.

Infiltration is a key component for storing water in soil layers within the hydrological process [8]. Soil water infiltration involves matrix flow and preferential flow [9]. Matrix flow represents the slow and even movement of water and solutes through the bulk soil [10], while preferential flow refers to the quick flow of water through the soil macrospores that bypasses most of the soil matrix [11]. At a soil profile, cracks, fissures, voids, fauna activities and root channels form a complex network that offers the ideal preferential pathways for water movement [12]. Through the preferential pathways, water and solutes can move to far greater depths than those predicted through the Richards’ equation [13]. As a result, the preferential flow could noticeably increase soil water conservation and groundwater recharge [14]. Dye tracer experiments are valuable techniques for visualising preferential water flow pathways [15]. Quantitative studies on preferential infiltration volume and its contribution to total infiltration and water storage remain nonetheless very limited.

Experiential evidence from field studies suggested that preferential flow processes varied considerably from one site to another [16]. Preferential water flow behaviour is primarily influenced by the soil properties, the rock fragments, the faunal activities, and the plant root systems [12,17,18]. With increasing abandonment time, the soil physical and chemical characteristics, the vegetation species and the plant root system in farmland can change significantly. Previous studies have confirmed that farmlands undergo significant increase in soil total porosity, soil organic matter content [19], soil clay content [20], total nitrogen [21] and aggregate stability [2] after a period of abandonment. Once farmland is abandoned, the vegetation recovers into tall herb, shrub, or forest ecosystems [3]. The plant species diversities and the plant communities vary greatly during the succession process [6]. The root system is the fundamental part of vegetation, and its mass density increases with the abandonment time [20]. Also, roots of different plant species, with diverse morphology, traits and distribution in soil layers tend to be vitally important for shaping the soil structure and the development from the soil surface to subsurface [22]. Those changes stimulated by the abandonment time certainly influence the soil hydrological properties and the water flow behavior. However, few studies have been conducted to investigate the temporal variations in the preferential flow of different vegetation communities after farmland abandonment.

Meanwhile, water flow moves into the soil layers through preferential pathways, which lead to the heterogeneity in the soil water through the soil profile [23]. Geostatistical methods were not only used to reveal the spatial variation of the soil water at slope and watershed scale but also at centimeter scale [24]. In this study, we investigate the vertical water flow behavior and the preferential infiltration based on a dye tracer experiment at four farmlands with different abandonment time and vegetation community. The aims of this study were threefold: (1) to identify the temporal variations in the soil properties and the root systems, and their relationship with the preferential flow with natural restoration age; (2) to contrast the water flow behavior in the soil profiles of the four abandoned farmland; (3) to assess the contribution of the preferential infiltration to the soil water storage during vegetation restoration.

2. Materials and Methods

2.1. Study Area

The study site was located in the Fangta catchment (36°46′13″–36°49′18″ N, 109°15′22″–109°17′01″ E) in the An’sai County of the Shaanxi Province, in North China (Figure 1). This region is representative of the loess hill and the gully region of the Loess Plateau. It has a semi-arid continental climate with a mean annual temperature of 8.8 °C, and mean annual precipitation of 505 mm, over 70% of which occurs from July to September [20]. Soils in this watershed mainly are red Loess (the Early or Mid-Pleistocene Loess) and yellow loess (the Late Plesitocene Loess). According to the World Reference Base for soil Resource (WRB), yellow Loess soil was classified as Calcaric Regosols and red Loess soil belonged to Calcaric/Chromic Cambisols. The dominated plant communities are composed of artificial forest and shrubs (Robinia pseudoacacia, Caragana korshinski and Hippophae rhamnoides), natural shrubs (Rosa xanthine and Sophora davidii), natural grasses (A. gmelinii, Stipa bungeana, Bothriochloa ischaemum, A. giraldii and Lespedeza davurica). Following the “Grain to Green” project, a significant portion of the steep farmland area in this region was converted to forestland or grassland.

2.2. Experimental Design

After a detailed field survey and several interviews with the residents, four abandoned farmlands with respective abandonment times of 3, 12, 18 and 25 years were selected in this study (

Table 1. Elementary information on the experimental sites.

Site	Altitude (m)	Slope Aspect (°)	Slope Gradient (°)	Slope Length(m)	Vegetation Coverage (%)	Dominant Communities
AF3	1157	NE80	11	43	35	Artemisia scoparia-Green bristle grass
AF12	1233	NE55	18	55	44	Artemisia sacrorum-Artemisia argyi
AF18	1254	NE65	20	66	37	Bothriochloa ischaemum-Artemisia sacrorum
AF25	1287	NE30	16	58	47	Periploca sepium Bunge-Artemisia sacrorum

Note: AF means the abandoned farmland and the number means the abandonment time.

). The four sites were covered by different plant communities of Artemisia scoparia, Artemisia sacrorum, Bothriochloa ischaemum and Periploca sepium Bunge, representing the succession sequence of the natural vegetation in this area [6]. The terrain, the soil type and the previous tillage practices of these four abandoned farmlands are similar except for their plant communities and their respective abandonment time. In each test site, three plots at the middle slope were selected as replicates. The selected plots were located at least 3 m apart to avoid interference. From July to August 2018, all the experiments were conducted on these four sites. Soil and plant root sampling, as well as dye tracer experiments, were completed in the fields where rainfall had not occurred for five days.

2.3. Dye-Tracer Experiment

Twelve dye-tracer experiments were conducted in the four test sites (four sites × three plot replications). In each plot, one hollow stainless-steel cube (length: 0.4 m; width: 0.4 m; height: 0.4 m) was inserted into the soil at 10–15 cm depth by using a rubber hammer. The aboveground plants were clipped to ground level and litter were removed carefully before the application of the dye tracer. The 12 plots were all irrigated with 10 L tap water which was dyed with Brilliant Blue FCF at a concentration of 4 g/L [24]. The falling head method was used to measure the soil infiltration process [5]. After irrigation, the surface of the stainless-steel cube was covered by plastic sheeting to prevent rainfall and evaporation. Twenty-four hours after the end of the infiltration, the plastic sheeting and the cube were removed carefully. The vertical soil profiles (width: 0.6 m; depth: 0.6 m) were excavated at the center of the cube. The dyed soil profiles were photographed using a digital camera (AF-S D5600, Nikon, Tokyo, Japan) under daylight conditions.

Soil sampling was conducted after the images were obtained through subdividing the profile into small grids of 2 × 2 cm (width × height) to extract soil samples with a volume of 4 cm³ (width × height × thickness: 2 × 2 × 1 cm), following the procedures thoroughly detailed by Mei et al. [24]. A total of 400 soil samples for each profile (width: 40 cm; depth: 40 cm) were dried at 105 °C for 24 h to measure the gravimetric water content and the infiltration volume.

The images were analyzed using the ERDAS IMAGINE v9.2 software (Erdas Inc., Atlanta, GA, USA) to identify the dyed coverage. According to the image processing procedures described by Cey and Rudolph [25], the images were geometrically corrected (Figure 2a) and subsequently classified into dyed and non-dyed regions. After that, the images were resampled to a 4000 × 4000 pixel resolution and stored as bitmap format for further analysis (Figure 2b). To quantify the vertical soil water flow path, the following parameters were determined as preferential flow indices [16] from the classified images of the soil profile using the Matlab 2015a software (Figure 2c):

Dye coverage (DC, %) referred to the percentage ratio of the dyed region area to the sum areas of the dyed and non-dyed regions, and was expressed as follows:

$$\mathrm{DC}=100·\left(\frac{\mathrm{D}}{\mathrm{D}+\mathrm{ND}}\right),$$

(1)

where D was the area of the dyed region (cm²), and ND was the area of the non-dyed region (cm²).

The Preferential Flow Proportion (PFP), expressed as a percentage, was the fraction of the total infiltration that flows through the preferential flow paths and was expressed as follows:

$$\mathrm{PFP}=100·\left(1-\frac{\mathrm{Unifr}·40}{\mathrm{TDA}}\right),$$

(2)

where PFP was the preferential flow fraction, Unifr was the uniform infiltration depths (cm) as the depth at which the dye coverage decreases below 80%, and was multiplied by the 40 cm wide profile in this study, TDA was the total dyed area (cm²).

The length index (LI) was the sum of the absolute differences between the consequent DC values with depth in the soil profile as follows:

$$\mathrm{LI}={{\displaystyle \sum }}_{i=1}^{100}\left|{\mathrm{DC}}_{i+1}-{\mathrm{DC}}_i\right|,$$

(3)

where i represents a given depth interval (or zone) of the 100 for which DC was calculated.

To partition preferential flow infiltration from the total infiltration, the total infiltration volume (TIV, mm), the preferential infiltration volume (PIV, mm) and the contribution of PIV to TIV (Con, %) was calculated as follows [24]:

$$\mathrm{Con}=100·(\frac{\mathrm{PIV}}{\mathrm{TIV}}),$$

(4)

$$\mathrm{TIV}={{\displaystyle \sum }}_{i=0}^{i=\mathrm{TID}}{\mathrm{Q}}_i,$$

(5)

$$\mathrm{PIV}=\mathrm{TIV} - {{\displaystyle \sum }}_{j=0}^{j=\mathrm{Unifr}}{\mathrm{Q}}_j,$$

(6)

$$\mathrm{Q}=10{{\displaystyle \sum }}_{k=1}^{k=20}\left({\mathrm{W}}_{dk} - {\mathrm{W}}_{ik}\right)/40,$$

(7)

$${\mathrm{W}}_k = {\mu }_k·{\mathrm{BD}}_k·v,$$

(8)

where Q (mm) represented the amount of the infiltration volume at a soil layer of i and j; TID (cm) was the total infiltration depth, Unifr (cm) was the uniform infiltration depths (cm); the W_dk and W_ik (mm) represented the soil water content after and before irrigation; k was the horizontal soil sample sequence, which had a maximum value of 20 because the soil was sampled every 2 cm within each horizontal soil profile of 40 cm wide; μ_k (%) was the gravimetric soil water content of the k sample; BD_k (g/cm³) was the soil bulk density of the layer; v (cm³) was the volume of the soil sample (cm³), which was 4 cm³.

2.4. Soil and Root Sampling

Prior to dye tracer experiment, soil properties including soil water content (SWC, %), the soil bulk density (BD, g/cm³), the particle size and the soil organic matter (SOM, g/kg) were measured in adjacent area to the dye plot. Soil samples were collected with three replications for four layers of 0–10 cm, 10–20 cm, 20–30 cm and 30–40 cm in each soil profile. The samples were transported to the laboratory for further analyzing the soil water content (SWC, %), the soil bulk density (BD, g/cm³), the particle size and the soil organic matter (SOM, g/kg). The SWC was measured with the oven-drying method (105℃ for 24 h). The BD was determined from the SWC and was calculated according to the weight of the soil in the cutting cylinders (inner diameter: 50.46 mm; height: 50 mm; Volume: 100 cm³). The particle size (clay: <0.002 mm; silt: 0.002–0.05 mm; sand: 0.05–2 mm) was measured using a Mastersizer 2000 laser diffraction device (Malvern, UK). The distribution of the soil aggregates was measured by the wet sieve method, and the water-stable aggregate content > 0.25 mm (WR_0.25) was calculated using the weight proportion. The SOM was determined by the potassium dichromate colourimetric method.

Additionally, the intact soil with root was sampled in each plot from 0 to 40 cm depth of the soil layers using a sampler that was built with a metal frame (length × width × height: 20 cm × 10 cm × 10 cm). The samples were sealed with a plastic bag and stored in a refrigerator for 12 h at 4 °C. The roots were carefully separated from the attached soil particles. The root images were taken at a resolution of 300 dpi by using a V700 scanner (Epson, Tokyo, Japan). The root length (cm), the root mean diameter (cm) and root surface area (cm²) and root volume (cm³) was measured using the Regent WinRhizo software (2009a, Regent Instruments, Quebec, QC, Canada). All the roots were subsequently oven-dried for 24 h at 65 °C to obtain the root biomass. The root mass density (RMD, mg/cm³), the root length density (RLD, cm/cm³), and the root volume density (RVD, cm³/cm³) were calculated as described by Yu et al. [26].

2.5. Statistical Analysis

The semivariance (γ(h)), which indicate the variance of several samples in a certain space interval, was adopted to characterise the spatial heterogeneity of the vertical soil water through infiltration in this study. It was expressed as [27]:

$$\mathsf{\gamma }\left(h\right) = \frac{1}{2N\left(h\right)}{{\displaystyle \sum }}_{i=1}^{N\left(h\right)}{\left[Z\left(x_i\right) - Z\left(x_i + h\right)\right]}^2,$$

(9)

where γ(h) was the semivariance; N(h) was the number of samples at each distance interval h, Z(x_i) and Z(x_i + h) were the values of the soil water at the site x_i and x_i + h, respectively. In this study, three geostatistical models were used to fit the experimental value of γ(h), including the Exponential, the Spherical and the Gaussian models. The best-fit model was selected according to the coefficient of determination (R²) and the root-sum square (RSS) of the cross-validation. The nugget variance (C₀), the sill (C₀ + C) and the range (A) were calculated based on the best-fit model. C₀ indicated the existence of measurement error or the variance at a shorter distance. (C₀ + C) represented the total variance of the variables. A represented the autocorrelation similarity of the variables. On account of the semivariance analysis, the ordinary kriging was adopted to map the spatial distribution of the soil water content. The above geostatistical analyses and the mapping were conducted using the ArcGIS 10.5 (ESRI, Redlands, CA, USA) and Surfer 15 software (Golden Software Inc., Golden, CO, USA).

The one-way analysis of variance (ANOVA) was used to compare the effect of different vegetation type on the preferential flow indices, the soil and the root properties. Pearson correlation coefficient was conducted to determine the relationship between the preferential flow characteristics and the soil and plant root properties. All the statistical analyses were conducted using IBM SPSS 19.0 (IBM, Armonk, NY, USA).

3. Results

3.1. Soil and Root Characteristics

According to the one-way analysis of variance, the soil properties differed substantially among the vegetation communities with different abandonment years following vegetation restoration (

Table 2. Mean soil properties (0–40 cm soil layer depth) of different sites.

Site	BD (g/cm3)	SWC (%)	SOM (g/kg)	WR0.25 (%)	Soil Particle Size Distribution
					Clay (%)	Silt (%)	Sand (%)
AF3	1.31 ± 0.03 a	10.07 ± 0.05 b	4.67 ± 0.03 b	29.45 ± 1.46 c	10.73 ± 0.01 a	23.15 ± 0.23 a	66.12 ± 0.73 a
AF12	1.27 ± 0.01 b	11.02 ± 0.08 a	4.60 ± 0.05 b	39.44 ± 1.93 b	10.80 ± 0.01 a	23.55 ± 0.09 a	65.66 ± 0.16 a
AF18	1.22 ± 0.01 c	9.34 ± 0.06 c	4.87 ± 0.05 b	46.64 ± 1.62 a	11.25 ± 0.01 a	24.62 ± 0.54 a	65.13 ± 0.51 a
AF25	1.14 ± 0.01 d	10.68 ± 0.01 a	5.28 ± 0.08 a	48.26 ± 2.16 a	11.54 ± 0.12 a	25.23 ± 1.09 a	64.23 ± 1.58 a

Note: Data represent means ± standard deviation (n = 3 replications × 4 soil layers = 12). Different letters in one column represent significant difference at p < 0.05 among different sites.

). The BD significantly decreased with increasing abandonment time. Compared with AF3, the BD of AF25 decreased by 12.98%. By contrast, the WR_0.25 in AF12, AF18 and AF25 increased by 33.92%, 58.37% and 63.87%, respectively. The changes in the other soil characteristics over time were more complex. Within 0–40 cm of the soil layer, the mean SWC was the lowest in AF18 and the highest in AF12. The SOM in AF3, which was abandoned rapidly, was a little higher than in AF12, but significantly lower than in AF25.

The spatial distribution of the RMD, RVD, RD and RLD in the four sites differed among the soil layers and the vegetation communities after the farmland abandonment (Figure 3a–d)). The four root parameters in all the sites decreased with increasing soil depth. Both the RVD and RMD increased with increasing abandonment time. The maximum values of RVD (0.0253 cm³/cm³) and RMD (3.30 mg/cm³) were found in AF25, which were 5.62 times and 1.76 times greater than the minimum value, which was found in AF3. For the RD and RLD, the AF3 and AF25 exhibited the minimum and maximum values, respectively; the same applied to the RVD and RMD. However, the RD in site AF12 was greater than in AF18, maybe because of the differences between the dominant plant species and their root systems. Bothriochloa ischaemum, the dominant species in AF18, has a fibrous root system; while Artemisia sacrorum in AF12 has a tap root system. The RLD of AF12 in 0–10 cm was greater than that of AF18 and lower than that of AF18 in deeper soil layers.

3.2. Infiltration Patterns

The twelve dye stained images of the vertical soil profiles in the four sites are shown in Figure 4. For all images, the non-uniform spreading of the stained region and the uneven wetting front implied that water flowed into the soil along the vertical direction through the preferential pathways and bypassed parts of the soil matrix. However, the stained dye patterns differed in detail among the four sites. In the four sites, the stained region was more evenly distributed, and the wetting front was more stable in AF3 (first row in Figure 4).

With the natural restoration ages, a less uniform distribution of the stained region and a sharper wetting front were observed at AF12 (second row in Figure 4), AF18 (third row in Figure 4) and AF25 (fourth row in Figure 4). Interestingly, parts of the stained region were distributed separately at AF18 (Figure 4g,i), and AF25 (Figure 4j,k), which indicated the infiltration of water into the soil in the two sites also followed other pathways from other directions.

The measured infiltration time and the preferential flow characteristics based on the dye stained images are shown in

Table 3. Dye stain characteristics in the vertical soil profiles for the twelve plots of four sites.

Site	IT	MID	Unifr	DC	FFP	LI	PIV	Con
	(min)	(cm)	(cm)	(%)	(%)		(mm)	(%)
AF3
Plot 1	29.60	23.2	18.9	50.44	6.32	273.75	5.03	13.85
Plot 2	32.50	24.6	20.8	55.28	5.93	327.50	3.92	11.46
Plot 3	29.50	19.9	14	35.12	0.34	280.00	6.82	26.34
Mean	30.31 ± 0.98 a	22.57 ± 1.19 a	17.9 ± 3.50 a	46.95 ± 10.52 a	4.2 ± 3.35 b	293.75 ± 16.97 a	5.26 ± 1.46 b	17.22 ± 7.99 b
AF12
Plot 1	20.60	29.5	15	43.86	14.5	436.75	16.19	41.58
Plot 2	29.60	24	8.1	37.21	45.58	242.00	17.06	53.2
Plot 3	23.70	20.4	11.4	37.15	23.28	226.75	9.79	35.59
Mean	24.63 ± 2.64 b	24.63 ± 4.58 a	11.5 ± 3.54 a	39.41 ± 3.86 a	27.79 ± 16.02 a	301.83 ± 69.14 a	14.35 ± 3.97 a	43.45 ± 8.95 a
AF18
Plot 1	43.20	25.9	13	44.26	26.57	402.00	13.38	40.71
Plot 2	35.50	32.8	21.2	55.44	4.4	361.50	12.79	31.2
Plot 3	36.60	22	4.5	23.88	52.89	185.25	20.27	72.21
Mean	38.01 ± 2.40 a	26.9 ± 5.46 a	12.9 ± 8.35 a	41.19 ± 16.00 a	27.95 ± 24.27 a	316.25 ± 79.41 a	15.48 ± 4.16 a	48.04 ± 21.46 a
AF25
Plot 1	25.70	18	9	28.12	19.99	386.00	9.2	41.99
Plot 2	27.60	25.3	5.9	35.51	41.54	331.25	30.16	72.37
Plot 3	17.90	22	2.8	22.65	30.91	291.75	25.07	79.33
Mean	23.65 ± 2.96 b	21.77 ± 3.66 a	5.9 ± 3.10 c	28.76 ± 6.45 a	30.81 ± 10.78 a	336.33 ± 27.32 a	21.48 ± 10.93 a	64.56 ± 19.85 a

Note: IT means the infiltration time, MID means the maximum infiltration depth, Unifr means the uniform infiltration depth, DC means the dye coverage to the total soil profile region, FFP means the preferential flow proportion, LI means the length index, PIV means the preferential infiltration volume, Con means the contribution of the preferential infiltration to the total infiltration. The mean value of each site represents mean ± SD, n = 3. The different letters for mean value in one column represent a significant difference at p < 0.05 among different sites.

. For infiltration time with the same volume of water (62.5 mm), it took about 23.65 min and 24.63 min to infiltrate in AF25 and AF12, respectively, which was significantly less than the time required for AF3 and AF18 to infiltrate, which was about 30.31 min and 38.01 min, respectively. Although the maximum infiltration depth (MID) was 18.00 cm at plot 1 in site AF25 and 32.8 cm at plot 2 in site AF18, the difference of the mean MID was not significant among the four sites.

The least DC was observed in the site AF25, whereas the largest one was observed in AF3. Meanwhile, the longest mean Unifr was found in AF3, indicating the water flowed uniformly through the soil matrix down to the deep depths in this site. However, the FFP, LI, PIV and Con exhibited increasing trends with increasing natural restoration ages. The mean FFP, PIV, LI and Con of the site AF25 was 7.34 times, 4.09 times, 1.17 times and 3.75 times greater than that of the site AF3.

3.3. Contribution of Preferential Flow

Figure 5 shows the vertical distribution of the difference value of SWC before and after irrigation among the four sites with different abandonment time and vegetation communities. With increasing soil depth, the SWC generally decreased for all plots in the four sites. However, the distribution pattern of the SWC was various due to the different water flow pathways network in every soil profile. The higher SWC appeared as the vertical belt was in the contour map, serving as the preferential flow. The belts occurred in the shallow layer of 0–8 cm in site AF3, and in the deep soil layer of about 0–16 cm in AF12, AF18 and AF25. The vertical distribution of the water infiltration volume induced by the dye-stained experiment is shown in Figure 6. The data shows that some water infiltration volume remained mostly below the Unifr depth than above. Besides, compared to AF3, the preferential infiltration volume in AF12, AF18 and AF25 increased by 1.73 times, 1.94 times and 4.09 times, respectively (Table 3).

Table 4. Geostatistical analysis of the main characteristics of the SWC in the four sites.

Site	Model	Nugget (C0)	Sill (C0 + C)	A(cm)	C0/(C + C0) (%)	R2	RSS
AF3
Plot 1	Sp	0.0404	2.60	27.22	1.55	0.894	0.0019
Plot 2	Ga	0.4526	3.56	28.83	12.71	0.999	0.0011
Plot 3	Sp	0.0305	0.57	20.85	5.36	0.934	0.0017
Mean		0.1745 ± 0.1460 a	2.24 ± 1.23 b	25.63 ± 4.22 b	6.54 ± 2.39 a	0.943	0.0016
AF12
Test 1	Ex	0.1545	3.64	42.21	4.25	0.934	0.0010
Test 2	Ex	0.1577	1.52	19.81	10.36	0.865	0.0021
Test 3	Ex	0.0651	2.27	54.57	2.87	0.891	0.0016
Mean		0.1258 ± 0.0303 a	2.48 ± 1.07 b	38.86 ± 9.70 a	5.82 ± 2.30 a	0.897	0.0016
AF18
Test 1	Ex	0.0585	3.65	49.57	1.60	0.914	0.0013
Test 2	Ex	0.0463	4.76	46.37	0.97	0.924	0.0011
Test 3	Ex	0.0646	1.24	26.62	5.23	0.854	0.0023
Mean		0.0565 ± 0.0054 b	3.22 ± 1.04 a	40.85 ± 7.98 a	2.60 ± 1.32 b	0.897	0.0016
AF25
Test1	Ex	0.0064	5.42	56.57	0.12	0.903	0.0016
Test 2	Sp	0.0074	1.93	36.93	0.38	0.910	0.0016
Test 3	Ga	0.1260	3.13	32.13	4.03	0.864	0.0056
Mean		0.0466 ± 0.0039 b	3.49 ± 1.02 a	41.88 ± 11.65 a	1.51 ± 1.26 b	0.892	0.0029

Note: Sp means the spherical model; Ga means the Gaussian model; Ex means the exponential model; RSS means the residual sum of squares. The mean value of each site represents mean ± SD, n = 3. The mean value of each site represents mean ± SD, n = 3. The different letters for mean value in one column represent significant difference at p < 0.05 among different sites.

further reveals the spatial distribution of the SWC by the parameters of the semivariogram analysis. For each plot, the best-fit model was performed as it exhibited the greatest R² and the lowest RSS. The small nugget variances (C₀) were observed in all the sites, suggesting that the random variations in the SWC were relatively small. The nugget percentage (C₀/(C + C₀)) values among the 12 plots were less than 25%, which indicates that the SWC of all the sites exhibited strong spatial dependence. The sill (C₀ + C) values and the mean spatial autocorrelation range (A) values increased from AF3 (2.24 and 25.63 cm, respectively) to AF25 (3.49 and 41.88 cm, respectively), and the C₀/(C + C₀) values decreased from AF3 (6.54) to AF25 (1.51). This suggested that the SWC was more variable as the abandonment time increased.

4. Discussion

4.1. Relationship between the Preferential Water Flow and the Soil and Root Features

The preferential flow occurs in most soils due to the presence of macropores [17], which is closely related to the soil physical properties. Based on the dye tracer experiment, the calculated parameters such as the FFP, PIV and Con (Table 3) could be adopted to estimate the degree of preferential flow. Generally, relatively loose, fragmented and macropore-rich soil were prone to preferential flow [23]. The results of the correlation analysis (

Table 5. Person coefficients between preferential water flow parameters and soil and plant root characteristics.

Items	Soil							Root
	BD	SWC	SOM	WR0.25	Clay	Silt	Sand	RMD	RVD	RD	RLD
IT	0.357	0.527	−0.401	−0.085	−0.308	−0.212	0.112	−0.358	−0.597 *	−0.493	−0.416
MID	0.199	−0.176	−0.332	−0.565 *	0.268	0.253	−0.237	0.542 *	0.804 **	0.613 *	−0.503 *
Unifr	0.707 **	−0.059	−0.638 *	−0.721	0.267	0.396	−0.306	−0.738 **	−0.780 **	−0.290	−0.674 **
DC	0.306	−0.108	−0.342	−0.424	−0.377	0.298	0.239	−0.342	−0.155	−0.485	−0.341
FFP	−0.715 **	0.082	0.652 **	0.686 **	0.151	0.236	−0.178	0.873 **	0.805 **	0.537 *	0.612 *
LI	−0.601 *	−0.083	0.562 *	0.726 **	0.213	0.144	−0.258	0.783 **	0.648 *	0.572 *	0.589 *
PIV	−0.713 **	−0.076	0.573 *	0.621 *	0.152	0.294	−0.214	0.818 **	0.733 **	0.659 **	0.711 **
Con	−0.675 **	−0.078	0.537	0.656 *	0.155	0.286	−0.192	0.778 **	0.735 **	0.607 *	0.728 **

Note: ** means the significant correlation at 0.01 level (p < 0.01); * means the significant correlation at 0.05 level (p < 0.05).

) showed that BD, WR_0.25 and SOM had significantly positive relationships with the degree of preferential flow, indicating that lower BD, higher WR_0.25 and SOM could enhance the preferential water flow. In general, the reduction of BD was associated with increasing soil porosity [28,29]. Along with decreasing BD, larger aggregates are associated with larger, more widely spaced and more continuous fissures [30]. Kaarakka et al. [31] suggested that SOM decreased the soil BD and increased the soil porosity, thereby inducing to an enhanced infiltration capacity. Hardie et al. reported that high initial SWC could limit the preferential flow [32], but the conclusion derived by Mei et al. [24] was rather the opposite. However, the results provided herein evidenced that the initial SWC had a weak relationship with the preferential flow parameters. Moreover, although Zhang et al.) found that preferential flow was likely to occur in sandy soil [17], soil texture exhibited weak relationship with the preferential flow in our study. The differences was probably ascribed to the different vegetation communities. In their study, the sites were covered by forest, shrub and grass, whereas the abandoned farmland covered by grass or subshrub were selected as our test sites.

Compared with the soil properties, the plant roots and their morphometric features tended to play more important roles on the preferential flow (Table 5), which is in agreement with Perillo et al. [33], who found that both decaying and living root systems were the most pivotal agents for preferential flow pathways. Cui et al. also found that fine roots were the main factor affecting the potential of soil infiltration in semi-arid regions [34]. During plant root growth, the plant roots not only grow into the preferential flow pathway, but they also create the pathways through biopores, cracks and burrows [17]. On the other hand, growing root systems were associated with reducing soil BD and increasing WR_0.25 and SOM in our study (Table 2 and Figure 3). This finding is supported by Benegas et al. [12], who found that the inputs of plant root were critical for improving soil properties during vegetation recovery, thereby motivating preferential flow. In the present study, all the four root features showed positive correlations with the content of the preferential flow. The distribution of the dye stained area corresponded to the distribution regions of the root system (Figure 3), and the preferential flow ended at the layers with fewer roots. The values of RMD, RVD, RD and RLD were the greatest in AF25 in comparison with AF3. Hence, the preferential flow was substantially stronger in AF25 than in AF3, which is consistent with the results reported by Zhang et al. [35], who found that RMD and RLD were greater in preferential pathways than in the soil matrix. Huang et al. [36] also revealed that RMD was the most important factor for the infiltration process. However, Cui et al. [34] found that non-perishable coarse roots could compact soil and block water flow, thereby causing more RMD inhibited soil infiltration process. The differences may be attributed to the variations in the root morphology and the root distribution in the soil profiles amongst the different plant species [22].

4.2. Variation Patterns of Water Flow Behaviors during Natural Vegetation Restoration

In the abandoned farmlands without ploughing and sowing, the soil physicochemical properties and the below-ground biomass changed greatly with vegetation restoration [3]. In the present study, as the abandonment time increased, vegetation restoration not only improved the soil and root properties (Table 2 and Figure 3), but also had different effects on water flow behaviors, as illustrated by the visual dye tracer information (Figure 4) and the preferential flow parameters (Table 3). Although the soil matrix flow was initiated near the soil surface, the preferential flow became the dominant flow behavior after years of abandonment. The degree of preferential flow from lowest to highest occurred in AF3, AF12, AF18 and AF25 (Table 3), following the increasing natural restoration ages. These results confirmed the vital effect of root systems on the formation of preferential flow pathways. Plant root mass density greatly increased with vegetation restoration age both at hillslopes and steep gully slopes in the Loess Plateau [4], which agreed with the present results (Figure 3). The AF3 site, which was dominated by the short-lived grass of Artemisia scoparia with small amount roots, had the lowest extent of preferential flow, suggesting that a significant portion of infiltrating water flowed uniformly through the soil matrix. This finding is in agreement with the conclusion of Jiang [13], who found that the matrix flow was the dominant flow type in the soil profile with no roots. Compared to the site AF3, the dominant vegetation in AF12, AF18 and AF25 were respectively the perennial grass, the subshrub and the shrub with more abundant and dense root systems, resulting in a greater preferential flow. Similarly, Alaoui et al. found that forestlands exhibited greater preferential flow than grasslands [10]. Moreover, old root channels and some wormholes were also found through the soil profile in the field (Figure 4e,h,i,l)), which could form numerous continuous macropores and preferential paths [37]. Old or decaying roots and animal activates are believed to increase with natural restoration ages.

Water flow bypassed the soil matrix and confined it to isolate patches, thereby exhibiting typical flow behavior in AF18 and AF25, which was known as lateral flow (Figure 4g,i,j,k)). In agreement with Sheng et al. [38], the lateral flow was found to be not obvious in the upper soil layers (<10 cm), while it may have been more predominate with increasing soil depth. Jiang et al. revealed that lateral flow occurred when a wetting front penetrated through the soil layers [23]. The factors that influence lateral flow include the soil matrix properties, the macropores, and the macropore-matrix interface [13]. The Bothriochloa ischaemum that dominated in AF18 possessed abundant fine roots and intricate fibrous roots. The channels of random directions created by the dense root network contributed to the non-uniform water movement and the irregular wetting front, likely causing the phenomenon of lateral flow. The Periploca sepium bunge in AF25, which had developed taproot and lateral root, offered potential lateral directions for water infiltration. Also, the relatively steeper slope may further promote the phenomena of lateral flow [39]. In the present study, Figure 4g shows that the isolated patches that occurred at the bottom of the soil profile was far from the upper macrospore flow in Plot 1 in site AF18 with slope gradient of 20°, which may be attributed the lateral water movement from the upper slope driven by gravity from the upper slope. This finding was distressed from those obtained by Jiang et al. [23], who found that compacted matrix by tillage pan and high SWC in cropland lead to the emergence of lateral flow.

4.3. Contribution of the Preferential Flow on Soil Water Storage

Spatial distribution maps of the soil water content better represented the water infiltration process and the preferential flow pathways than the dye-stain images [24]. As shown in Figure 5, the belt with high moisture content was associated with macropore flow paths. The depths of those belts were prone to be deeper as the degree of preferential flow was greater with increasing abandonment time. This result is consistent with the finding reported by Malvar et al. who found that water may concentrate in macrospores and cause preferential flow paths, thereby leaving the surrounding soil matrix relatively dry [40]. Further, it is worth noting that the infiltration volume below Unfir tended to increase with increasing soil depth (Figure 6), which therefore demonstrates that the preferential pathways in the soil profile could promote the delivery of water flow into the deep soil layer. This result agreed with Hanson et al. who reported that preferential flow paths and large macrospores carried most of the water when the soil became saturated [41]. Previous studies indicated that several scattered areas with high water infiltration might be caused by the rock fragment of soil in forests [42] and tillage practice in croplands [23]. However, this phenomenon was herein related to the lateral flow, which was probably induced by the non-uniform void of the soil created by the plant roots and worm holes.

Based on semivariograms analysis, the nugget percentage values of C₀/(C + C₀) were less than 25% in all sites, indicating that the water flow pattern had strong spatial dependence among the four sites with different natural restoration time. Moreover, the spatial heterogeneity of the SWC increased from AF3 to AF25 with the increasing abandonment time. The stronger spatial heterogeneity of the SWC may be induced by the greater degree of preferential flow with the non-uniform soil properties and root systems [43]. Therefore, the improved preferential flow led to the vertical and lateral redistribution of the soil water, and increased water storage in the soil profile, which was in favor of the plant growth and sustainability of vegetation restoration.

Additionally, dye tracer images do not quantify soil macopores and 3-D soil macropore networks, although it can visually reveal the vertical water flow patterns and preferential flow pathways. Future research should prefer to the methods of X-ray computed tomography (CT) and dye tracer along the different direction.

5. Conclusions

This study quantified the temporal variations in water flow behaviour of different vegetation communities in abandoned farmlands on the Loess Plateau. The results showed that the soil organic matter content (SOM) and the water stable aggregate (>0.25 mm) content (WR_0.25) increased with natural restoration age, while the bulk density (BD) decreased from AF3 to AF25. The natural restoration age combined with the plant species induced a positive effect on the root mass density (RMD), the root volume density (RVD), the root diameter (RD) and the root length density (RLD). These improved soil properties and developed root systems substantially influenced the water flow behavior and enhanced the degree of preferential flow.

The results of the dye tracer experiment and the image analysis indicated that the mean FFP, PIV, LI and Con of AF25 were 7.34 times, 4.09 times, 1.17 times and 3.75 times greater than that of AF3. The correlation analysis implied that the plant roots and their morphometric features exhibited a greater effect on the preferential flow in comparison with the soil properties. The abundant lateral root and the steeper slope may have caused the presence of lateral flow that occurred in AF18 and AF25.

Moreover, the contour map of the soil water content showed a belt of high soil water content spread into the deeper soil layer with longer restoration ages. The spatial variability of the soil water through the vertical soil profiles increased from AF3 to AF25 with increasing degree of preferential flow. Also, the contribution of the preferential flow to the total infiltration increased with abandonment time. Therefore, the improvement of the preferential flow in the abandoned farmland during natural vegetation restoration helped soil water storage in the deep soil layer.