The study area, Jianpinggou catchment, is located in southwest China’s Sichuan Province (Figures 1A,B), has an area of about 3.5 km², and ranges in elevation from 1,051 to 2,199 m above sea level. Derived from the 5-m resolution digital elevation model (DEM), the average slope is 32.4° (ranges from 0° to 63.3°). The main channel length is 1.7 km. The climate is subtropical and humid, with average annual precipitation of 1,134 mm, and about 80% of this occurs in the rainy season (from May to September) (Ran et al., 2015). The areas are covered by the perennial trees and shrubs (Figure 1C). The soil in the catchment is classified as loam (sandy loam and silty loam), with typically a 50-cm-thick silty clay and silt loam horizon, a 250-cm-thick sandy loam horizon, and an underlying bedrock. The soil horizons are full with bedrock fragments. The third horizon, bedrock, is composed of the basalt, tuff, and andesite rocks, which are of low permeability.

The main rainfall gauging was located about 2 km (960 m elevation) downstream of the study area with an accuracy of 0.2 mm/min. Velocity was measured at the catchment outlet by the large-scale particle image velocimetry (LSPIV) system at every 45 min from 8:00 am to 20:00 pm in rainy season (Ran et al., 2016; Li et al., 2019). The hinge of LSPIV is to calculate the displacement of the natural tracer particles on the water surface in the paired images to obtain the surface flow field. There were two long-term observation sites of soil volumetric water content in the catchment, one was on the hillslope (HS) close to the outlet, and the other was on the ground (GND) around 1 km upstream of the outlet. The monitoring depths and frequency were 100 cm in total at every 10 cm, and a 5- to 60-min time interval, respectively. Based on the dielectric method, the soil moisture sensor determines the soil moisture according to the change of the soil dielectric constant. In the soil, the dielectric constant is determined by the soil matrix, water, and air, whose dielectric constants are about 3–5, 81, and 1 at room temperature, respectively. Therefore, the soil dielectric constant is mainly influenced by the soil water content. Soil properties were sampled along the main channel by the cutting ring method, mainly measuring the saturated hydraulic conductivity, porosity, and mechanical composition. Although the field measurements were mostly located at mid-downstream in the study area, the model performances were consistent with upstream and downstream. For example, Beville et al. (2010) successfully employed InHM to estimate the spatiotemporal pore pressure distributions (from upstream to downstream along the watershed) for Lerida Court landslide in Portola Valley, CA, United States.

Previous studies had explored that the hydrologic response here was dominated by the Dunne overland flow and subsurface flow (Ran et al., 2015). In our recently field surveys, stones appear commonly in the soil, ranging in diameter from several centimeters to several meters, and there were also some exposed rocks on the surface. The vertical saturated hydraulic conductivity (Ks) and porosity of the first horizon are 2.35*10⁻⁵, 1.44*10⁻⁵, and 9.79*10⁻⁵ m/s and 0.65, 0.61, and 0.33 for the >1,300-, 1,180 to 1,300-, and <1180-m-elevation region, respectively. The corresponding values of the second and the third horizons are 1.33*10⁻⁴ and 1*10⁻⁸ m/s, and 0.45 and 0.30 for the whole region, respectively. The hydraulic conductivity of the second horizon is almost 1 order of magnitude larger than that of the first horizon, promoting the lateral flow in the soil.

The installed soil moisture sensors interrogated about 0.32 m³ of soil to measure volumetric water content (m³/m³) and soil temperature (°C). The data loggers recorded sensor measurements hourly or more frequently when solar energy was more efficient. We replaced volumetric water content values equal to 0 with ‘NaN’ values (Gasch et al., 2017), and confirmed that all records had the correct date and time stamps.

The outlet discharge was estimated by the velocity–area method, and the three-dimensional topography (mainly to get the cross section) was developed by the stereo imaging method in which the paired images were also obtained from the LSPIV system (Ran et al., 2016). However, in case of flood, the stereoscopic imaging effect was poor, and the cross section was calculated by measuring riverbed elevation and the section water table. Those interested in a detailed description of LSPIV system are referred to Fujita et al. (1998).

The study by Mirus and Loague (2013) examined that when the soil water passes through the preferential channel, the change of water content in the lower layer may be faster than that in the upper layer. Based on the above analysis, as shown in Figure 4, we supposed two overlapping triangular prisms to represent two adjacent unit soils, and the inward/outward arrows represented the amount of water entering/draining the unit soils in all directions. The center of each unit soil (i.e., the solid dot) was the node participating in the model simulation, and the nodal area was the surrounding area. The edge of the upper unit soil was drawn with red solid lines, and the adjacent lower layer was blue. The interior was filled with the same color as the edge. The color shade displayed the degree of soil moisture at the corresponding time.

In the simulation, the study area was divided into thousands of elements [not triangular prism, see more detailed in Figure 3-b and A-1 of Vanderkwaak’s Ph.D. dissertation Thesis (1999)]. If we compared the storage change of soil water content of all paired nodes between the lower layer and its corresponding upper layer in a set interval (i.e., the occurrence of NSR), the subsurface flow paths could be inferred from the plot scale to catchment scale approximately by the velocity difference of the flow (the control distance of the paired nodes is equal). Considering the not-so-fast velocity of subsurface flow and the monitoring frequency of the soil moisture sensor, the mentioned time interval above was selected as 5 min.

Hence, we used the following equations to implement it. S C j = θ t i + Δ t , j − θ t i , j , (1) D S C = S C j + 1 − S C j , (2) where θ was volumetric water content, m³/m³; the subscripts t i and j were the corresponding time i and the jth node, respectively; the storage change ( S C j ) was the change of soil water content of the jth node in Δ t , m³/m³; the difference of storage change ( D S C ) was the difference between the storage change of node j+1 in Δ t and that of its corresponding upper layer node j, m³/m³. In a certain period ( Δ t ), if D S C > 0, the SC in the lower layer was more massive than that in the upper layer, more water flowed through the lower layers by subsurface lateral path, the NSR occurred.

Many researchers have utilized topographic algorithms, like HAND (Height Above the Nearest Drainage), to identify hydrological similarity which revealed strong correlation between soil water conditions and topography (Renno et al., 2008; Gao et al., 2019; Loritz et al., 2019). However, there would be no gravitationally driven water movement between two hydrologically connected points with the same height in HAND. Considering that the distribution of soil moisture would directly affect the distribution of NSR, we finally chose topographic wetness index [TWI = ln (a/tan b)] as the topographic index. Although it was sensitive to DEM resolution (Sorensen, and Seibert, 2007; Gao et al., 2019), it was more in line with our needs. The TWI was utilized to describe topographic characteristics, which was a combination of the upslope area per unit contour length a and the local slope tan b (Beven and Kirkby, 1979), computing by the multiple flow direction algorithm. In this study, we chose the method of Yang et al. (2000), which calculated the distributions of topographic index from DEM using a single flow direction algorithm offered by ARC/INFO software. Besides, the cardinal direction of the contour length ignored the diagonal of the gird.

The value of a, which indicated the area per unit contour length, was computed for every grid in the catchment as follows: a = ( n u g + 1 ) ∗ A C , (3) where n _ug is the number of upslope grids, A is the grid area, and C is the contour length, given by the following: C = l , (4) where l is the grid size at the cardinal direction.

The magnitude of tan b, which is a measure of the potential drainage from a place, in the steepest downslope direction was calculated as follows: tan ⁡ b = d H d L , (5) where dH is the change in elevation between neighboring grid cells and dL is the horizontal distance between centers of neighboring grid cells.

Rinderer et al., 2014, Rinderer et al., 2019, and Jencso et al., 2009 showed that sites with a TWI>6, a local slope <30%, and an upslope contributing area >600 m² were defined as footslope sites. Upslope sites were defined as sites with a TWI<4, a local sloe >50%, and an upslope contributing area <200 m², whereas midslope sites had a TWI between 4 and 6, a local slope between 30 and 50%, and an upslope contributing area between 200 and 600 m². The slope of the above classification conditions was not considered in the mountainous area. The study area was divided into the channel (CH), hillslope (HI), and upland (UP) zones by TWI and upslope contributing area, with the corresponding average slope of 25.8, 33.2, and 35.8°, respectively (Figure 5).

The two measures of model performance used to test InHM were modeling efficiency (EF) (Nash and Sutcliffe, 1970) and coefficient of determination (R ² ). The mathematical expressions for EF and R2 are as follows: E F = [ ∑ i = 1 n ( O i − O ¯ ) 2 − ∑ i = 1 n ( P i − O i ) 2 ] / ∑ i = 1 n ( O i − O ¯ ) 2 , (6) R 2 = 1 − ∑ i = 1 n ( O i − P i ) 2 ∑ i = 1 n ( O i − O ¯ ) 2 , (7) where O i are the observed values, P i are the predicted values, n is the number of samples, and O ¯ is the mean of observes data. When the observed and simulated values are identical the EF and R ² statistics are both 1.0.

The simulation of hydrological response was continuous in the study. Figure 6 was the rainfall-runoff events used to drive the InHM simulation of hydrograph response for study area. The minutely rainfall record in Figure 6 was a month from June 24 to July 27, 2018, containing 13 observed runoff events. The characteristics of driven rainfalls and model performance are summarized in Table 2. The observed data were mostly concentrated on the recession limb because the flood often occurred at night in the mountainous area. Besides, the image velocimetry performed better at low flows, while the simulated values matched better with the observed discharge. Imaging that a rainfall event occurred at night and far away to the flow gauging station, the measured data and the simulated value (assuming that the rainfall has always covered the whole study area) might have a difference in time and space, resulting in low or out of reality Nash coefficient and R ² . The simulated values were underestimated compared with the observed values. For example, the vegetation on the road was washed away in the image/video of the next morning in the first rainfall event (start from June 24, 2018), which indicated that the flood should brim over the culvert at night. Although the maximum water stage has been used in the calculation, it was still lower than the real value. It could be seen from Table 2 that mean rainfall intensities were less than the soil saturated hydraulic conductivity, indicating that the main mechanism of runoff generation was saturation excess (Dunne) flow in the study area. The time since last rainfall event ranged from 13 to 115 h. The performance of the model before July 5, 2018, was worse than the latter, which might be due to 1) the instability of measured discharge at high flows and 2) the discontinuity of the spatial and temporal distributions of rainfall in mountainous areas.

Figure 7 showed the observed DSC of GND and HS at different depths, from June 24 to July 27, 2018. The fluctuation of GND was more intense than that of HS, indicating the higher occurrence frequency of the NSR. The occurrence of the NSR varied between >0 and 4%. Event 1 had the widest effect depth in longitudinal, and the NSR was monitored even at 100 cm depth below the surface. The NSR of the GND and HS fluctuated greatly in events 1, 4, 7, 10, and 11. Meanwhile, the fluctuation of the 0–50 cm depth was larger than that of 50–100 cm depth, especially at HS. The duration of events above ranged from 4 to 74 h, and the amount between 37 and 233.5 mm. The mean intensity was relatively heavy, and the time since last rainfall event was about 20–30 h.

The amplitude of the simulated DSCs (Figure 8) was smaller than that of the observed. The corresponding nodes of every layer at HS in simulation were hard to match the measured points because of coarse-mesh accuracy for the outlet. But it was worth recalling that the simulation results showed a similar trend with the observation in the events mentioned above, implying that the simulation could reflect the real situation to a certain extent.

From June 24 to June 25, 2018, we recorded a heavy rainfall event. The videos showed that the plants on both sides of the road (Intersecting with the outlet of the study area) were still there at 20 p.m. on June 24, but they have washed away the next morning. It was speculated that the night flood exceeded the road elevation at the weir. The heavy rain seriously destroyed the main road of Longxi River, and we chose the event to study the spatial and temporal distributions of the NSR on a catchment scale. Figure 9 demonstrated the double-peak precipitation simulation started from 15:56 on June 24, 2018. The "+" represented the times corresponding to interest points of the hydrograph along the runoff process (Times 1–11). The "o" denoted the times corresponding to the peaks, the troughs, and about a day after the event along the runoff process (times A to E). In the two obvious runoff peaks, the peak time of the former (∼7.1 h) was longer than that of the latter (∼2.0 h). Table 3 summarized the characteristics of double-peak precipitation. The peak intensity of the two periods was very close. The first period had the larger rainfall depth and the shorter duration compared with the second one, and there was an 8-h interval between them. The mean intensity of the first period was twice bigger that of the second.

The map of spatial DSC, as shown in Figure 11, was pictured at all layers for times 1 to 11 in simulation. The proportion in Table 4 was the ratio of the NSRs area to the whole area. The position was determined according to the location of the NSR corresponding to the topographical zone in Figure 5. The table gave more comprehensive information for analyzing the movement of the NSR. At the beginning of the event, the runoff increased slightly and immediately, and there were almost no NSRs. The depth of the NSR spread from the middle layer to the upper layer in the rising limb of the first period of the event, while the position of the NSR gradually covered the channel zone from the hillslope zone. The proportion of the NSR decreased the channel zone of the upper layer contracted in the falling limb. The development of the NSR in the second peak of the event was similar to the previous one. Under the influence of the pre-event, the proportion of the NSR in the deep layer increased, the distribution expanded from the middle layer to the deep layer (times 6, 8, 9, and 10), and the rainfall condition corresponding to the NSR changed (times 4 and 8). During a period of time after the precipitation stopped, although there was a large area occurrence in the NSRs in the catchment, the position was the channel zone, which was different from the hillslope zone in precipitation.

The duration and intensity of the rainfall strongly affect the occurrence of the NSR, regardless of the terrain (Figure 7). In the condition of similar mean rainfall intensity, the greater the rainfall amount, the more fluctuated the soil water content in the upper layer (Table 2). This finding is in agreement with other studies (Hardie et al., 2011; Hardie et al., 2013; Koestel and Jorda, 2014; Wiekenkamp et al., 2016). For instance, Zhu et al. (2014) demonstrated that the maximum change of soil water content in the event was mainly controlled by the precipitation amount and intensity in Taihu Lake Basin, China. The two monitoring sites had a higher occurrence of the NSR and a higher maximum NSR occurrence (the max change of the DSC) in the upper layer (0–50 cm depth), where soil properties then changed at 50–60 cm depth, forming a relatively impermeable layer. Demand et al. (2019) found that the soil with a higher clay content enhanced the occurrence of the NSR. The NSRs concentrated on the 0–50 cm depth and went on to the deeper layer as the rainfall amount increased. Albertson and Kiely (2001) also found that the root-zone soil moisture variation was influenced by the precipitation intensity in a humid study area in Virginia, United States. In the events on June 24 and July 21, 2018, the rainfall amount and duration were 233 mm, 42.7 h, and 37 mm, 4.4 h, respectively, and the mean intensity was almost an order of magnitude higher than others. The NSR appeared at 100 cm depth during the two events, showing that the occurrence depth of the NSR increased with precipitation intensity when the time since last event was more than 24 h and the total amount of this event exceeded 37 mm in this paper. Similarly, Buttle and Turcotte (1999) argued that they did not find a relationship of preferential flow (nonsequential response flow) with initial soil water content, but with throughfall intensity. The NSR always appeared together with the sequential response (Figures 7, 8). These findings denoted that the rainfall amount and intensity affected the NSR always played a role together and how to separate the individual influence needed the further study.

The hillslope zone was never saturated at all depths during the event (Figure 10); thus, the discharge at the outlet was generated as flow over the zoot-zone layer and lateral flow through the subsurface layer rather than by Dunne or Horton overland flow. The upper soil layer responded shortly with a rise in saturation after the rainfall start, and the runoff ascended rapidly (time 1 in Figure 9). However, the flow over and through the litter layer was a minor component of the runoff and never comprised more than 0.5% of rainfall (Buttle and Turcotte, 1999; Whipkey, 1965). As shown in Figure 9, runoff increased only slightly at time 1. The soil-bedrock interface accumulated water vertically from bottom to top, forming a condition to initiate the subsurface lateral flow (times A to E at 250 cm depth in Figure 10). Furthermore, the lateral flow at the soil-bedrock interface moved from the hillslope of the deep layer to the channel at surface. As reported by Peters et al. (1995), runoff at the trough about 0.1 m above the bedrock surface was the upper portion of flow at the soil-bedrock interface that was recorded when the saturated layer, initiated at the bedrock surface, rose above the elevation of the trough.

As the rainfall went on, more water entered into the soil, some of which infiltrating into the upper layer moved above the relatively impermeable layer (60–50 cm depth), and some seeping into the middle layer tended to meet the water accumulated at the bedrock interface. Meanwhile, the NSR started to occur in the upper and middle layer (times A and B in Figure 11). The water of the upper layer drained out as soon as rainfall stopped due to the small volume. But at the same time, the deep layer was still in the state of receiving the lateral soil water, resulting in the occurrence of the NSR at the hillslope. In the middle layer, the NSR mainly appeared at the steep area because it was conducive to the convergence of the subsurface flow.

During the double-peak precipitation, the initial soil saturation became wet at the beginning of the second peak stage due to the pre-event. More water assembled at the bedrock surface compared to the before. The occurrence region of NSR extended upward in depths (Figure 11, compare times A and C), and soil moisture converged from the hillslope to the depression channel. When the second period rainfall began, the water at the soil-bedrock interface rose toward the relatively impermeable layer, causing a higher occurrence frequency of the NSR at 60–50 cm depth, and a relatively lower in other zones (times 4 and 6 in Table 4 and Figure 11). According to Hardie et al. (2013), high frequency soil moisture monitoring of the ∼90 cm depth between June 23 and July 28, 2008 in Australia demonstrated that occurrence of preferential flow (NSR) became less apparent as each rainfall event increased the antecedent soil moisture content. The higher occurrence frequency of NSR during the first rainfall period may be related with dry soils that develop water repellent (Bouma, 1991).

When precipitation occurred, the soil water content in the upper layer responded quickly. The water through the relatively impermeable layer penetrated down to the soil-bedrock interface vertically, the riparian zone expanded, and the lateral flow moved downslope in the form of a near-saturated wedge (Mosley, 1979; Mosley, 1982). The NSR appeared at deep layer in the channel (Figure 11, compare times B and D), and the runoff rose. After the precipitation stopped, the water on the slope continued to move laterally above the relative impervious layer and the soil-bedrock interface, and the NSR mainly appeared on the hillslope (Figure 11, compare times C and A). Besides, the water in the deep layer was hindered at the relatively impermeable layer when moved to the surface. After passing through this layer from bottom to top, the flow exfiltrated out the surface together with the water originally accumulated here (time E in Figure 10). The NSR appeared in the hillslope zone, and the runoff fell. During the rainfall and a period of time after the rain stopped, the water still accumulated at the soil bedrock interface, and this accumulated water might relate to the subsurface storm flow (Markus Weiler et al., 2006; Sidle et al., 2000; Mirus and Loague, 2013).In case of rainfall event with two peaks, pre-event increased the antecedent soil moisture content in deep layer (250 cm depth in Figure 10), the region observed with nonsequential response expanded (Figure 11, compare times A and C). The soil layer at the interface of bedrock could be saturated quickly, and became saturated upwards. This kind of nonsequential response can be observed on the hillslope (Table 4, compare times 2 and 5) at the beginning of rainfall events, and then found beneath stream channels afterwards (Table 4, compare times four and 6). Furthermore, nonsequential response could also happen after rainfall events the constant accumulated water at the soil-bedrock interface.