Spatio-temporal relevance and controls of preferential flow at the landscape scale

The spatial and temporal controls of preferential flow (PF) during infiltration are still not fully understood. As soil moisture sensor networks allow us to capture infiltration responses in high temporal and spatial resolution our study is based on a large- 10 scale sensor network with 135 soil moisture profiles distributed across a complex catchment. The experimental design covers three major geological regions (slate, marl, sandstone) and two land covers (forest, grassland) in Luxembourg. We analyzed the responses of up to 353 rainfall events for each of the 135 soil moisture profiles. Non-sequential responses ( NSR ) within the soil moisture depth-profiles were taken as one indication of bypass flow. For sequential responses maximum pore water velocities ( v max ) were determined from the observations and compared with velocity estimates of capillary flow. A measured 15 v max higher than the capillary prediction was taken as a further indication for PF. While PF was identified as a common process during infiltration it was also temporally and spatially highly variable. We found a strong dependence of PF on the initial soil water content and the maximum rainfall intensity. Whereas a high rainfall intensity increased PF ( NSR , v max) as expected, most geologies and land covers showed highest PF under dry initial conditions. Hence, we identified a strong seasonality of both NSR and higher v mat dependent on land cover, revealing a lower occurrence of PF during spring and increased occurrence 20 during summer and early autumn, probably due to water repellency. We observed the highest fraction of non-sequential response ( NSR ) in forests on clay-rich soils (slate, marl). Maximum pore water velocities ranged from 6 cm day -1 to 80640 cm day -1 with a median of 120 cm day -1 across all events and soil moisture profiles. The soils in the marl geology had the highest


Introduction
Preferential flow (PF) in soils describes different flow processes with higher flow velocities than soil matrix flow and heterogeneous flow patterns (Hendrickx and Flury, 2001). Many studies have shown that PF is ubiquitous (Jarvis, 2007) and that "PF is the norm and not the exception" (Weiler 2017). PF can affect water distribution in soil (Ritsema et al., 1996), groundwater recharge (Ireson and Butler, 2011), root water uptake (Schwärzel et al., 2009) and solute transport (Larsbo et al., 5 2014). Since the early work of Beven and Germann (1982), the importance of PF pathways such as macropores (created by roots, earthworms), fissures or cracks is widely recognized. Most studies focusing on different PF processes, such as fingered flow (Selker et al., 1992), macropore flow (Weiler and Naef, 2003) or funnel flow (Kung, 1990), were carried out at the point or plot scale (spatial scale smaller than a few meters). Since PF increases the range of flow velocities in the vadose zone by orders of magnitudes (Nimmo, 2007), it is essential to include this process when modeling water and solute transport in soil. 10 (Blume et al., 2009;Eguchi and Hasegawa, 2008;Germann and Hensel, 2006;Hardie et al., 2013;Kim et al., 2007) or for analyzing the sequence of their response with depth (Graham and Lin, 2011;Lin and Zhou, 2008;Liu and Lin, 2015;Wiekenkamp et al., 2016). Using these methods most studies found a relationship with precipitation characteristics (Liu and Lin, 2015;Wiekenkamp et al., 2016) or initial soil moisture (Blume et al., 2009;Hardie et al., 2013;Liu and Lin, 2015;Wiekenkamp et al., 2016). 5 Even though some of the studies described above show differences in PF occurrence between soils or landscape properties, most of them do not rigorously compare contrasting landscape units at the larger scale. Zhao et al. (2012) tested out-ofsequence responses of the soil moisture sensors as an indication of PF for two contrasting land covers and found much higher occurrence of PF in the forest sites compared to a cropland. However, since both sites also had different soils it could not clearly be attributed to land cover. Most field experiments studying the effect of soil texture and land cover on soil water flow 10 measured infiltration characteristics or hydraulic conductivities of soil cores (Bormann and Klaassen, 2008;Gonzalez-Sosa et al., 2010;Jarvis et al., 2013;Zimmermann et al., 2006). In general, higher infiltration rates and hydraulic conductivities were observed at sites with natural vegetation or forests. These higher infiltration rates were often attributed to the presence of macropores, but not connected to the dynamics of PF occurrence under natural field conditions. Studies linking the spatial and temporal PF occurrence in high resolution and comparing contrasting landscapes under natural initial and boundary conditions 15 are still scarce.
A correct estimation of PF occurrence is important for hydrological predictions (e.g. modeling) and can improve water resource management. Therefore, the main aim of this study is to identify and compare the temporal dynamic of PF occurrence by using profiles of soil moisture sensors in different large-scale spatial units that could potentially be used as representative units for catchment modelling. Since it can be expected that rainfall intensity and soil moisture have a strong influence on the 20 initialization of PF (Beven and Germann, 1982) we will mainly focus on the temporal controls of initial soil moisture and rainfall. More specifically, we attempt to answer the following question: Does PF occurrence increase with rainfall intensity since higher intensity leads more frequently to an exceedance of matrix infiltration capacity? Does PF occur more often under wet conditions since the infiltration capacity is lower? How is the temporal PF dynamic influenced by spatial factors like geology/soil type and land cover? 25

Study Sites
To test our research question we analyzed a dataset of 405 soil moisture sensors at 45 sites distributed across a complex landscape (varying geology and land cover) but under similar climatic conditions. The monitoring network is located in the 30 Attert catchment in the Grand Duchy of Luxembourg. The climate is temperate semi-oceanic with a mean annual rainfall of 845 mm (Pfister et al., 2006) and mean monthly temperatures between 0°C (January) and 17°C (July) and only very few days per year with snow coverage (Wrede et al., 2015). Elevation ranges between 265 and 480 m a.s.l. and the catchment covers three major geologies (Colbach and Maquil, 2003). The northwestern part of the catchment is located at the southern edge of the Ardennes and the geology here is dominated by Devonian Slate bedrock covered by periglacial slope deposits mixed with eolian loess (Juilleret et al., 2011;Moragues-Quiroga et al., 2017). The southern part of the catchment is dominated by 5 sedimentary rocks of the Paris Basin (Wrede et al., 2015) with Jurassic Luxembourg Sandstone at the southern catchment border and Triassic Sandy Marls in the central part of the catchment (Fig.1). The slate region has agriculturally used plateaus between steep forested slopes (~15-25°). The sandstone hillslopes are mostly forested with grasslands only present on the footslopes (Juilleret et al., 2012;Martínez-Carreras et al., 2012). The land cover in the Luxembourgian part of the marl region is mainly characterized by agricultural sites (30 %) and grasslands (41 %, mainly pasture) with gentle slopes (~3°). 10 43 sensors were replaced with SMT100 (TRUEBNER GmbH, Neustadt, Germany) and 9 sensors with GS3 sensors (Decagon Devices/METER Environment, USA) in 2016. Sensors were installed horizontally with minimum disturbance from a 30 cm diameter hole drilled with a power auger. Each sensor was installed slightly shifted in the horizontal direction to the one above, to be unaffected by potential flow path changes by the sensor above. Furthermore, sensor cables were laid downwards in the hole first and led up on the opposite wall to prevent artificial PF along the cables leading to the sensors. In each of the three 5 main geologies, the sensor sites were situated in two different land cover classes, forest and grassland. The selected forest sites were dominated by European beech (Fagus sylvatica) with occurrence of oak (Quercus robur, Quercus petraea) and common hornbeam (Carpinus betulus). Furthermore, rainfall was measured with one tipping bucket (Davis Instruments, USA, 0.2 mm resolution, collection area 214 cm²) at each grassland site and five randomly placed tipping buckets at each forest site to account, to at least some degree, for the spatial variability of throughfall. We defined six different landscape units 10 distinguishing the three main geological formations and the two land covers (forest, grassland) to test our research questions.
The number of soil moisture profiles for the different land cover and geological classes are summarized in Table 1. Additional information and specific site properties are shown in Appendix A.

Rainfall events
A full workflow of the data analysis is depicted in Fig. 2 showing the number of excluded events due to different quality criteria. Rainfall (P) events were defined using the rainfall data with 5-minute resolution individually for each site. For the forest sites the mean of all five tipping buckets for every 5-minute time step was calculated to obtain average throughfall for each site. Forest tipping buckets that measured no rainfall over one hour were excluded (assuming they were clogged), when 10 at least three other buckets observed rainfall during the same timeframe. If the rainfall data contained more than one missing value in a 2-hour period it was excluded from further analysis. Following the approach of Graham and Lin (2011) and Wiekenkamp et al. (2016), a rainfall event was defined as rainfall with a minimum amount of 1 mm. The end was defined as the last monitored response of a rain gauge followed by a specific time period without rain (te). The procedure of determining this time period is described below. 15 Dividing soil water dynamics into single events based on P input is always a trade-off: On the one hand, short rainfall events do not allow for a clear separation of the infiltration signals from different input pulses. On the other hand, long rainfall that is grouped into one event can result in too much information from several consecutive rain input pulses that are merged into one rainfall event. Hence, different rainfall regimes require different threshold values, i.e. hours without rainfall (te) for the identification of event endings. The sensitivity of te on the number of rainfall events and their characteristics in our case was 20 investigated by testing different values of te: 3, 6, 12 and 24 consecutive hours without rain.
For each P event total rainfall amount (Psum), the maximum P intensity in a 5-min time step (Pmax) and the event average rainfall intensity of the entire event (Pint) was determined. Events that were not plausible were excluded by using a threshold method for event P amount (Psum > 100 mm), average event intensity (Pint > 15 mm h -1 ) and maximum P intensity in a 5minute time step (Pmax > 80 mm h -1 ). These implausible events were observed to happen during the reconnecting of the loggers following a logger error (no power etc.) or clogging and release of the clogged water. To exclude snowfall or frozen soil conditions, events with a mean air temperature below 0°C during the event were not included in the analysis. By applying the quality criteria for rainfall events using te = 12 h, 1392 of 32025 rain events (sum of profile rainfall events) were excluded 5 because of the threshold criteria and 426 because the mean temperature was below 0°C during the event.

8
The rainfall event separation method is sensitive to the required number of consecutive hours without rain (te) between the events. Table 2 shows te values with the resulting number of events, mean event duration, rainfall amount (Psum) and event average rainfall intensity (Pint). Shorter te results in more events and decreasing mean event duration. Mean Pint is gradually decreasing with longer te due to longer event durations while mean Psum is increasing. We considered te = 12 h to be sufficient to ensure event separation yielding an appropriate event length and to avoid possible superimposition of soil water flow signals 5 from different input pulses. Therefore, the following analyses are performed with the event definition based on te = 12 h. This results in total rainfall event numbers between 144 and 353 per profile. 54.2 % of all analyzed rainfall events had sums lower than 5 mm and 77.7 % lower than 10 mm. The distribution of rainfall intensities (Pint) shows that 69.2 % of all events had a Pint < 0.4 mm h -1 . The density distributions show slightly higher Pmax for grassland sites but no difference among the geologies (Fig. 3). 10

Soil moisture & infiltration events
Signal spikes (strong increase of soil moisture within a 5 minute time step and a decrease to the initial value) in the measured soil moisture time series were removed and data was visually checked for plausibility and long-term consistency. In addition, sensor readings were validated against those of the other sensors in the same depth for each site. No site specific calibration of the soil moisture sensors was conducted and soil moisture values were obtained by the sensor internal θ-permittivity 5 relationship following Topp et al. (1980). For the 5TE sensors the manufacturer gives an absolute sensor accuracy of volumetric water content of ±3 Vol% (DecagonDevices, 2015). For a relative change of 1 Vol% a maximum sensor-to-sensor difference of ± 0.25 Vol% can be found in the very dry range (θ ~ 10 Vol%) (Rosenbaum et al. 2010). Since Rosenbaum et al. (2011Rosenbaum et al. ( , 2012 showed that temperature effects on the sensors and on soil dielectric properties can cancel each other out, permittivity was not corrected for soil temperature. Furthermore, electrical conductivity effects of soil water on permittivity 10 were neglected as bulk electrical conductivity was low (< 0.1 dS m -1 ) for most profiles. Although some marl profiles show higher bulk electrical conductivities, results of soil water content change should not be affected since these profiles do not reveal fast bulk electrical conductivity fluctuations on the event scale.
For each defined rainfall event the soil moisture time series of all sensors in a profile was checked for their response. Infiltration events were defined as a θ increase of ≥1 Vol% of at least one sensor in the soil profile. This threshold was chosen to avoid 15 diurnal fluctuation, caused by e.g. soil temperature, being classified as infiltration events (Graham and Lin, 2011;Wiekenkamp et al., 2016). If a soil moisture event was identified, the timing of first response of every sensor was determined. The first response is defined as the point in time when the θ change is higher than the instrument noise (Lin and Zhou, 2008) that was found to be 0.4 Vol% for the 5TE sensors (Rosenbaum et al., 2010;Wiekenkamp et al., 2016). Linear interpolation was used to calculate the time between two 5 min readings to increase the temporal resolution. The soil moisture response was tracked 20 for up to 48 hours after the end of a rainfall event or until the time a new rainfall event starts.
The chosen rainfall event separation based on te = 12 h already avoids superimposition of consecutive rainfall input signals on the soil water content. However, to have clearly separated soil water flow events that are uninfluenced by a new rainfall event for at least 24 hours, both consecutive infiltration events were excluded if the second rainfall event occurred within 24 hours after the first rainfall event end. In the case of a response later than 24 hours we assumed that the following infiltration event 25 is likely to be triggered by the new rainfall event (Hardie et al., 2013). Only if more than 99 % of the data points for all profile sensors during an infiltration event were usable, they were considered for further analysis (termed completeness criterion).
Furthermore, infiltration events that showed an increase in soil moisture, but were caused by an oscillating signal (not more than four different θ values during one event) were excluded (termed consistency criterion).
From the total of 30207 rainfall events, 15645 could be used for the analysis of the soil moisture, since they allowed for a clear 30 separation of soil water flow by more than 24h without a new rainfall input. 7395 of these events did not meet the quality criteria of completeness and consistency of the soil moisture time series, hence 8250 infiltration events (sum of soil moisture event observations at all 135 profiles) could be used for the analysis. Changing the completeness criterion from 99% usable soil moisture data points during an event to e.g. 95% is only slightly affecting the number of infiltration events (e.g. 8353 events usable in the analysis). This is due to the fact that most exclusions result from long term failure of one sensor of a profile that leads to a complete exclusion of the entire profile. A diagram showing the portion of active (all quality criteria met) 5 profiles on a daily basis can be found in the supplement ( Figure S1).
Various soil moisture and rainfall characteristics were determined for each event. Initial volumetric water content (θini) was defined as the water content before the rainfall event starts. Furthermore, change of θini to the peak water content (Δθmax) of every event and sensor response was calculated. We grouped soil moisture into dry and wet initial conditions using θ quartiles of each profile. Additionally, rainfall amounts and intensities were calculated for the time before the first soil moisture sensor 10 response (∆θ = 0.4 Vol%) of any profile (rPsum, rPint, rPmax). This was done since our infiltration event classification described in the next section is partly based on the first sensor response and later rainfall input is not further influencing the classification.

Soil moisture sensor response by infiltration events
For all soil moisture profiles and rainfall events which met the described quality criteria, the sequence of the first sensor 15 response was classified similar to Liu and Lin (2015) into: (i) not classifiable (NC): none of the sensors in the profile showed a response (≥1 Vol%) or only a 10 cm sensor response was observed (ii) non-sequential response (NSR): events where the first response did not progress in a sequence starting from the surface (e.g., the 30 cm sensor showed a response before the 10 cm sensor) 20 (iii) sequential response (SR): the sensors in the profile showed a response in sequence from the uppermost sensor downwards (e.g., 10 cm to 30 cm to 50 cm or 10 cm to 30 cm) The potential for using these different infiltration responses (SR, NSR) and related parameters as a proxy for PF are described in the following sections. All statistical analysis were performed using Dunn's rank sum test (Dinno, 2017).
Additionally, we estimated how often PF should have be observed based on the classical assumption that rainfall intensity 25 exceeded matrix infiltration capacity (Beven and Germann, 1982). We used matrix saturated hydraulic conductivity (Kmat) as the minimum infiltration capacity and tested how often maximum 5-min rainfall intensity exceeded this threshold (Pmax > Kmat; for the measurements of Kmat see section 2.2.2.2). Furthermore, comparison of maximum water content change during an event (Δθmax) between the infiltration response types can give information on PF processes by showing differing water content depth distributions and can help to estimate the relevance of the different flow processes in terms of transported water quantity.

Non-sequential response (NSR)
The NSR classification indicates non-uniform flow that can be a result of various PF processes (e.g. bypass flow), hence it is taken as a proxy for PF. NSR could also be a result of subsurface lateral flow or groundwater rise before the vertically downward progressing wetting front reaches that depth (Lin and Zhou, 2008). But even in these cases, such responses describe 5 water flow that shows either non-uniform flow or surroundings that infiltrate water faster than the profile. Both can be seen as an indication of PF. None of the profiles showed a permanent water table smaller than 50 cm below ground level, nevertheless some profiles are influenced by groundwater fluctuations and are temporary waterlogged in 50 cm especially during winter.
The length of the time series is adequate for detecting patterns of NSR as Liu and Lin (2015) showed in their analysis that overall sensor response patterns show stable results using >3 years of soil moisture data. The occurrence frequency of NSR 10 was analyzed with respect to initial soil moisture and rainfall characteristics for the landscape units. All NSR analyses were done with pre-response rainfall characteristics (rPsum, rPint, rPmax). Calculated portions of NSR for the landscape units, geologies or land covers for different rPmax, θini or month are always calculated as the sum of NSR responses of the indicated class divided by total number of infiltration events in the same class.

Sequential response (SR)
A sequential response of the sensors in the profile does not necessarily mean that no PF occurred. To get an estimate for the frequency of SR events showing PF, one method is comparing soil matrix (capillary) flow velocities to measured in-situ flow velocities (Germann and Hensel, 2006;Wiekenkamp et al., 2016). A measured flow velocity that is faster than the soil matrix flow velocity can be expected to be influenced by PF. Matrix flow velocity can either be obtained by modeling or with 20 measurements. To determine the in-situ flow velocities we used the approach of Germann and Hensel (2006) where the maximum pore water velocity (vmax) is determined from the first responses of two sensors (often called wetting front velocity).
The upper sensor allows for the definition of a clear starting time of the water flow. Hence, vertical maximum pore water velocities were calculated from the SR for two distinct flow depths: 10 to 30 cm and 30 to 50 cm. It is important to note that vmax represents only the fastest flow components in the sphere of influence around the soil moisture sensor (Hardie et al., 2011). 25 To model matrix flow velocity (vmat), the 1D steady state flow equation according to Darcy's law for unsaturated conditions was used (Hillel, 1998): With q being the vertical volume flux [cm day -1 ], K the hydraulic conductivity [cm day -1 ], ψm the matric potential [cm], H the hydraulic potential [-] and z the depth [cm]. For the vertical 1D case, matrix flow velocity (or piston flow velocity) can be calculated by dividing the volume flux by the volumetric water content θ [-] (Gerke, 2006):

= /
The hydraulic gradient was calculated between two sensors using the matric and gravitation potential (H = ψm + ψg). The 5 maximum gradient between the θ peak of the upper sensor and θini of the lower sensor is calculated to obtain maximum vmat. This is a conservative approach since steady state assumptions are used to calculate flow velocity. For obtaining the matric potential the van Genuchten retention curves (van Genuchten, 1980) were parameterized using the parameter sets of Sprenger Table S2). The van Genuchten parameters of Sprenger et al. (2016) do not need further corrections for matching θ with absolute values of e.g., soil core data since these parameters were calibrated for a shorter period of the 10 same dataset. For those ten sites where no parameters were determined by Sprenger et al. (2016), we simply used the mean for the respective geology. Although these retention parameters were inversely fitted and should therefore account for fast flow components, they more closely represent matrix flow due to the single domain Richards equation and the unimodal nature of the van Genuchten retention function that was used (Durner et al. 1994). In addition, the fit on a daily basis does not allow for fast processes other than matrix flow. A geometric mean hydraulic conductivity was calculated between two sensors located 15 in different depths (Zhu 2008) to obtain the effective unsaturated hydraulic conductivity of the vertical layered soil profile. To again provide a conservative estimation of PF and rather overestimate vmat the moisture content used to calculate this unsaturated hydraulic conductivity was the maximum event water content, determined for both sensor depths individually. The mean of these two maximum event water contents was also used to calculate the matrix flow velocity (vmat) from the volume flux (q). Events that showed an upward hydraulic gradient based on this calculation were excluded from further comparisons. 20 For directly measuring matrix flow velocity we assumed that saturated matrix hydraulic conductivity at the surface is an appropriate threshold for dividing flow into matrix flow and PF (Wiekenkamp et al., 2016). Tension infiltrometer measurements were used to obtain saturated matrix hydraulic conductivity in the field. The tension infiltrometer used in this study is a special type called "hood infiltrometer". The advantage of the hood infiltrometer is that it can be placed directly on the soil surface without need of any contact material (Schwärzel and Punzel, 2007). The derivation of matrix saturated 25 hydraulic conductivity (Kmat) from measured infiltration rates accounts for the three-dimensional nature of flow using the solution of Wooding (1968) (steady state infiltration from a circular source). Measurements were carried out either in the direct vicinity of our sensor sites or within the same geology and land cover class (Appendix A). All values of matrix surface hydraulic conductivity consists of at least three measurement locations (median), except for two sites where the infiltration rate was too high and the hood could not be filled. Hood infiltrometer measurements were not available for grassland sites in 30 the sandstone and hence observed flow velocities of this landscape unit were not compared with measured matrix flow (2) velocities. In total measurements from 66 locations were used to determine Kmat for the different landscape units. For every measurement location infiltration rates with at least three tensions between 0.4 -5.9 hPa were recorded to be able to fit an exponential function to calculate surface hydraulic conductivity at a tension of 6 hPa (Gardner, 1958). At this tension, pores with a diameter ≥ 0.5 mm are excluded from flow and measured hydraulic conductivities represent matrix infiltration capacities (Jarvis, 2007;Schwärzel and Punzel, 2007). Due to the high macroporosity at many forest locations pressure in the hood was 5 difficult to adjust and measurements could only be conducted for maximum tensions of 1-3 hPa. Hence, for some sites matrix saturated hydraulic conductivity is just an extrapolation of the Gardner fit to a tension of 6 hPa.

Infiltration events 10
The number and proportions of classified infiltration event responses (NC, SR, NSR) of the six defined landscape units are shown in Table 3. The absolute number of events in a certain landscape unit and response class, which were included in the different analysis, can be found in the supplement (Table S3) showed a change in soil moisture deeper than 10 cm. Most classifiable infiltration events were of type SR. Under sandstone forest sites they accounted for 24.6 %, whereas under marl grassland sites they accounted for only 13.6 % of all events. Within the group of SR, 47.4 % were observed at a depth of 30 cm, whereas sequential flow to sensors at 50 cm depth was found for 52.6 % of the SR. NSR events occurred in 5.3 % to 16.1 % of all events depending on the landscape unit. The slate and marl forest regions showed the highest proportion (13.3 % and 16.1 %, respectively). In total 48.7 % of the NSR events showed a 20 response in 30 cm first and 23.9 % in 50 cm. 27.4 % of the NSR events reacted in 10 cm first and then in 50 cm without a 30 cm reaction in-between. The NSR variability between the single profiles within a landscape unit was found to be high (Table   3). The site-intern variability of NSR (profiles within the same sites) measured as the median standard deviation was highest in marl (forest: 7.5 %, grassland 6.4 %) followed by slate (forest: 4.2 %, grassland 6.1 %) and sandstone (forest: 1.9 %, grassland 3.0 %). 25 To estimate how often PF should have be observed based on the classical assumption that rainfall intensity exceeded matrix infiltration capacity in the different landscape units we calculated the portion of rainfall events with a Pmax exceeding Kmat.
With the exception of marl grassland (13.8 % Pmax > Kmat), all other landscape units only showed an exceedance rate lower than 2 % (Table 3). To test how much P characteristics and θini influence the different response behaviors, we calculated the median of each parameter for all infiltration events of a certain response type and their corresponding depth (Table 4). We included pre-5 response P characteristics (rP) to show their differences between NSR and SR events. High Psum mainly affect the depth of the soil moisture front during SR. In addition, the Pmax is also increasing with depth of response, which could partly be due to a correlation of Pmax and Psum (Spearman R = 0.54). SR events show similar median θini values for both infiltration depths, which suggests no effect of θini on the flow depth. The rPsum is similar for SR and NSR 30 and 50 cm events, while rPmax is higher for NSR events. NSR10-50, with a response in 10 cm first followed by a 50 cm reaction, shows a different pattern than the other 10 NSR reactions with the lowest rP intensities, but the highest θini.and rPsum. In contrast to SR the median θini of the NSR events is lower and also decreases with increasing depth of first response (30, 50 cm), which indicates that these infiltration response type is sensitive to dry soil moisture conditions.

Water content change
To estimate the relevance of the different response types in terms of the transported water quantity through the soil, the maximum change in water content for every event (Δθmax) has been taken as a proxy which can further indicate differences in response properties. The patterns of Δθmax in each geology were compared among response type and depth. Figure 4 shows 10 violin plots with Δθmax in the two individual depth during SR. For SR the plots include all events that show a response in the respective depth, independent of the maximum response depth. For NSR 30 and 50 cm events only Δθmax of the first response depth was considered in the respective depth. For NSR10-50 only the water content change in 50 cm (first out-of-sequence reaction) was taken into account. Observed median Δθmax values range between 1.8 and 4.3 Vol%. For the SR events, a significant decrease of Δθmax with depth was observed for slate and sandstone sites. Marl sites did not show this damping of 15 the water content signal with depth and exhibited a significant increase of Δθmax at 50 cm depth (SR). For the NSR events no damping of Δθmax with depth was observed. In contrary, NSR in sandstone and marls both had higher Δθmax at 50 cm depth compared to 30 cm. Furthermore, for all geologies Δθmax at NSR 50 cm was similar or even stronger than for NC/SR 10 cm or SR 30 cm responses.

Non-sequential response (NSR)
The fraction of NSR events in dependence of θini and P characteristics was analyzed to reveal the spatial and temporal patterns and possible controls of PF. Pmax and θini of each profile are only weakly correlated (median profile Spearman R: -0.19). An 10 increase of NSR with increasing rPmax was observed (Fig. 5). Especially forested sites in the slate and marl region showed a strong increase of NSR above a threshold of rPmax = 10 mm h -1 . This pattern was only weakly pronounced for the grassland sites. More NSR with higher rPmax in the forests was also found when using maximum rainfall intensity for the whole event (P) instead of the pre-response characteristics (rP).  Figure 6 shows the portion of NSR response for the six different landscape depending on individual θini quartiles for every profile to account for the differences in absolute θini values among landscape units. We observed that the drier the forested 5 sites were, the higher the measured NSR occurrence was. Especially slate and marl sites showed a strong increase in NSR occurrence (up to ~25 % of events) for the driest θini quartile. At slate grassland sites observed NSR occurrence was not responding to drier conditions in the same way as for the forested sites. The fraction of NSR events at the marl grassland sites did not change with initial conditions and at sandstone grassland sites NSR occurrence increased only under wetter conditions.  (Table S4).

Number of events observed in the different classes can be found in the supplement
To test for a seasonal effect on the NSR occurrence we also analyzed the frequency of NSR on a monthly basis. Since land cover seems to play an important role for NSR occurrence (Fig. 5 and 6) the NSR portion for all infiltration events of the two land covers was calculated separately. Forests show a distinct seasonal dynamics (Fig. 7): From March to June NSR showed a constant value slightly higher than 5 % which increases to 13-20 % from July until October and decreases again towards 5 winter. In the same time period θini dropped to its lowest annual values and also rPmax had its maximum in the summer month.
For grasslands this dynamic was less pronounced with the highest value in September.  (Table S5).

Estimating PF by comparison with modeled and measured matrix flow 15
To identify PF from SR we further compared measured maximum pore water velocities (vmax) against measured (hood infiltrometer, Kmat) and modeled matrix flow velocities (vmat). Table 5 indicates the percentage of observed vmax that exceed either the measured infiltrometer or modeled values. Both comparisons indicate that observed water flow is in most of the cases faster than water that is flowing in the soil matrix only. Between 72.9 and 89.0 % of the observed SR responses are faster than the modeled matrix flow velocities. The median difference in flow velocity for the events with vmax > vmat is 114 cm day -1 . 20 The model matches with the exceedance obtained by the hood infiltrometer, except for marl and sandstone forest sites with an exceedance rate of the infiltrometer being only 48.7 and 44.0 %, respectively. This is due to the high surface Kmat values that were measured with the hood infiltrometer for these two landscape units. The high conductive parameters of these two landscape units were not distinct higher in the set of hydraulic parameters used for modeling.

Observed maximum pore water velocities
Since the vmax observed from soil moisture responses (SR) exceeded the modeled or measured matrix values most of the time we examined vmax in more detail. The measured vmax ranged from 6 to 80640 cm day -1 with a median of 120 cm day -1 . Only a 10 weak correlation was found between vmax of the shallow versus the deeper depths (10-30 cm to 30-50 cm; Spearman-R: 0.36).
Median observed vmax values per group ranged between 72 cm day -1 for forested sandstone sites (for the shallow depth 10-30 cm) and 274 cm day -1 for forested marl sites (for the depth 30-50 cm) (Fig. 8). Comparing vmax for all landscape units the marl soils showed more variable flow velocities and higher median values, especially between 30 and 50 cm soil depth. Slate soils do not show a significant difference between the two depths or the land covers. Sandstone exhibited highest flow velocities 15 under grassland sites. Forested sandstone soils had a significant lower SR flow velocity than all other soils. To further evaluate the variability of vmax in respect to θini and Pmax for all observed events, 2D kernel density estimations 5 (KDE) (Venables and Ripley, 2002) are shown in Figure 9 with higher KDE values indicating more events. There is no clear relationship of vmax with θini or Pmax and high maximum pore water velocities can be found over the full range of θini and Pmax. Analyzing the median response of vmax to θini and Pmax for the different landscape units we can see an increase of median vmax for high Pmax for most landscape units (Fig. 10). Furthermore, the median vmax is increasing under dry conditions for marl independent of land cover and for slate grassland (Fig. 11). The other landscape units do not show a clear pattern between vmax and θini.  (Table S6).

10
Although, the relationship of vmax with Pmax and θini is not as clear as with NSR, the seasonal dynamics of median vmax shows an increase during the summer months with the highest flow velocities during times with low θini and high Pmax. In contrast to 22 NSR grasslands showed a stronger increase than forests with a maximum between June and August and a median vmax between 225-325 cm day -1 . For forests a weaker increase in the time between July and August and a stable median vmax of around 200 cm day -1 was seen. The number of observed events furthermore indicates that most SR events are not observed during the times of high vmax but rather during the wet winter month.   (Pmax, θini). The importance of PF during infiltration was highlighted by the fact that observed SR flow velocity (vmax) was most of the time faster than pure soil matrix flow and depended on the landscape unit NSR accounted for 18-44 % of the responses deeper than 10 cm. The variability of response types within the landscape units and even within some sites 15 was high, which highlights soil heterogeneity on such larger scales and shows the importance of small scale soil properties and soil water flow. However, we found a strong variability in PF occurrence that was dependent on spatial and temporal factors which are discussed in section 4.4 and 4.5.
PF is not only important in terms of its occurrence frequency, but also relevant for the quantity of transported water as indicated by the observed water content changes (∆θmax). Especially during NSR the ∆θmax is higher than ∆θmax for SR at the same depth, 5 which implies fast flow of large amounts of water into deeper zones. Furthermore, the marl sites with their high velocities in 50 cm depth also showed the strongest ∆θmax increase in this depth, unlike the other geologies. Similar observations were made by Hardie et al. (2013), who found higher ∆θmax in greater depth during NSR or events with high vmax, and Eguchi & Hasegawa (2008) calculated that high amounts (16 to 27 %) of the total annual drainage was produced by PF.

Observed non-uniform flow (NSR)
In our study, occurrence of NSR for single soil moisture profiles (0-46.2 %) was similar to other studies. Liu and Lin (2015) found profile NSR occurrence varying between <1 and 72.4 % for single years, Graham and Lin (2011) found 18 to 54 % for a three-year period and Wiekenkamp et al. (2016) found 7-51 % also using a three-year time series. However, we found a lower average NSR occurrence (mean of the profiles within one landscape unit) of 5.9-14.6 % for the landscape units in our 15 study (data not shown) compared to 26% in the Shale Hills catchment of Graham and Lin (2011). Until now, most studies on NSR events from soil moisture time series focused on a relatively similar substrate (shale), land cover (forest) and a temperate climate (Graham and Lin, 2011;Lin and Zhou, 2008;Liu and Lin, 2015;Wiekenkamp et al., 2016). The slate forest of our study is the landscape unit most comparable to the studies cited above. It shows a comparable range of NSR occurrence (0-46.2 % for a single profile). As our experimental design targeted not one but 6 different landscape units, we were able to 20 compare responses observed in the shale forest to other environments. Sandstone grassland showed a maximum NSR at a single profile of only 15.6 % of the events. Soil profiles under forest on clayey soils (slate & marl) had a higher occurrence of NSR (based on the landscape units) and a higher maximum NSR occurrence for single profiles within these landscape units compared to sandstone or grassland sites. Zhao et al. (2012) also found difference in land cover (forest vs. cropland) and soil characteristics to affect NSR occurrence. They found lower values with 5.8 -32.4 % NSR in the croplands compared to the 25 nearby Shale Hills forest, but as geology differs between the sites the lower NSR cannot be unequivocally attributed to land cover. Nimmo, 2007). In addition, studies that measured vmax in single sprinkling experiments in the slate forest region of the Attert catchment observed a vmax of 864-19000 cm day -1 using GPR and TDR during a hillslope irrigation experiment with an intensity of 30.8 mm h -1 (Angermann et al. 2017). Jackisch et al. (2017 observed vertical transport velocities of bromide in the range of 2732 cm day -1 with sprinkling intensities of 30 and 50 mm h -1 . The highest maximum pore water velocity of the slate forest landscape unit measured in our study was with 14662 cm day -1 in a similar range. 5

Observed maximum pore water velocities
Most of the studies mentioned above are sprinkling experiments which apply high P intensities (>10 mm h -1 ) and high Psum and thus do not provide information on the response to low intensity events that make up a large portion of the annual rainfall events (see Fig. 3). In his review, Jarvis (2007) found that solute transport studies were either carried out at (near-) saturated conditions or with high irrigation rates (> 10 mm h -1 ). Langhans et al. (2011) found an increase of infiltration capacity with higher rainfall intensity, probably due to the initiation of more macropore flow. This could be an explanation for the higher 10 velocities found by high intensity sprinkling experiments. Therefore, a reason of partly lower vmax observed in this study might be that we are also accounting for low P intensity events due to our focus on natural rainfall events. This assumption is supported by the fact that Hardie et al. (2013) measured a vmax of 24 -960 cm day -1 under natural rainfall conditions which is more in the range of most velocities observed in our study. In summary, it is remarkable that no clear differences in flow velocities between different soil types could be identified (neither in our study nor across all previous studies). Instead, all soil 15 types showed a similarly large range of velocities (10 0 -10 5 cm day -1. ). Furthermore, one can see orders of magnitude difference in vmax between different events but not among the landscape units. A clear reduction of maximum pore water velocity with decreasing θ (dry soils) as predicted by conventional unsaturated hydraulic conductivity relationships (e.g. van Genuchten, 1980) was not observed under field conditions. In contrast, higher flow velocities during driest conditions were observed for most profiles in our study. 20

Temporal controls of PF
We found that both, a low initial soil moisture (θini) and a high maximum rainfall intensity (Pmax) affect the occurrence of PF.
This results in a higher occurrence of PF during summer time. Increased PF (NSR, vmax) during low θini is in contrast to the classical assumption of PF, which should be initiated more often under wet initial conditions with a lower infiltrability. 25 Furthermore, the mismatch of measured PF occurrence (NSR, fast vmax) compared to the prediction based on Pmax exceeding Kmat indicates that initiation processes such as hydrophobicity/water repellency, local microtopographic depressions or channeling of water by vegetation could be the reason of the frequent occurrence of PF (Blume et al., 2008;Doerr et al., 2000;Schwärzel et al., 2012;Weiler and Naef, 2003). Locally, these processes can lead to higher water contents and thereby pressures at the soil surface close to atmospheric pressure which in turn trigger PF. The higher probability of NSR under dryer (2015). Also Hardie et al. (2011) found faster flow velocities under dry conditions, which they concluded, was due to hydrophobicity and resulting finger flow and Blume et al. (2009) found response time of soil moisture and thereby flow velocity to be much faster during summer time. However, Buttle and Turcotte (1999) did not find a relationship of PF and initial soil water content, but on throughfall intensity.
Due to the strong seasonal variation with a maximum in summer and early autumn ( Fig. 7 and 12), the most probable 5 explanation is the influence of water repellency that has frequently been observed on natural surfaces in summer (Doerr et al., 2006;Täumer et al., 2006). Water repellency hinders infiltration and ensures a pressure buildup at the soil surface until pressure reaches a positive water entry potential (Bauters et al., 2000). Gimbel et al. (2016) observed that their clayey and loamy plots developed strong water repellency during a simulated drought field experiment with a 40-year return period and that infiltration patterns changed from homogeneous to preferential flow. Also sandy soils were found to be strongly affected by water 10 repellency (e.g. Ritsema et al., 1997). Wessolek et al. (2008) found from a one year TDR measurements on a pine stand that PF is minor from February to April since the soil was not water repellent. They found a maximum of PF from May to September which matches in general with our observations, just that our observed maximum starts and ends approximately one month later. Furthermore, Täumer et al. (2006) observed a similar seasonal pattern over a 3-year period with a maximum of PF in summer and early autumn and also Rye and Smettem (2015) observed a similar seasonality in Australia. That during these dry 15 and water repellent conditions the P intensity is highest further supports the initialization of PF. In general higher P intensities can lead to water pressures at the soil surface close to the water entry potential (Gjettermann et al., 1997;Jarvis, 2007;Weiler and Naef, 2003).

Clay content
Examining the temporal effects of θini and Pmax between the landscape units in detail, PF dynamics were not the same throughout all landscape units in our study. Especially clayey soils seems to be strongly influenced by low θin and clay content enhance NSR occurrence and vmax. Many studies showed that the clay content increases macroporosity under dry conditions through shrinkage and the subsequent cracking of the soil (e.g. Li and Zhang, 2011;Novák, 1999;Stewart et al., 2016a). Das 25 Gupta et al. (2006) measured high infiltration capacity for the macropore domain of clay soils using a tension infiltrometer.
The higher macroporosity of the clay soil can then further enhance the occurrence of PF, initialized by higher Pmax and hydrophobic condition in summer as observed by (dye) tracers, infiltration and soil moisture measurements (Dekker and Ritsema, 1996;Hardie et al., 2011;Jarvis et al., 2008). Liu and Lin (2015) found clay content to be an important predictor of NSR in the Shale Hills catchment and we also measured higher NSR in clayey landscape units (slate, marl). Furthermore, we content) is probably more attributed to the high abundance of biopores observed in the topsoil of this region. The high flow velocities in the subsoil are in accordance to other studies that showed fastest velocities due to structure development in unsaturated clay soils (Baram et al., 2012;Hardie et al., 2011;Tiktak et al., 2012). Probably ponding of water on top of the clay layer and subsurface initiation of macropore flow could be a reason of higher flow velocities in the subsoil (Weiler and Naef, 2003). Such a process was observed in the field by Hardie et al. (2011). This demonstrates that in the unsaturated zone 5 close to the surface, clay should not be treated as a low conductivity but rather as a high conductivity material.

Land cover
The question arises why NSR is much more often observed in forests during summer compared to grassland and why vmax is higher in grassland. In general, forests tend to have highly connected macroporosity caused by roots (Alaoui et al., 2011;Gonzalez-Sosa et al., 2010;Lange et al., 2009). Furthermore, higher soil organic carbon content in forest can enhance aggregate 10 stability and hence interaggregate porosity in clayey soil (Lado et al., 2004;Six et al., 2002). However, the sole presence of a higher macroporosity in forests does not explain the higher NSR occurrence. That higher macroporosity results in more NSR could also be caused by more laterally directed pathways in forests created by roots as observed by Bachmair et al. (2009).
Funneling of rainfall by stemflow (not measured in this study) may support this mechanism (Schwärzel et al., 2012). In contrast, the stronger increase of vmax in grasslands during summer could be an indication of a seasonally changing 15 macroporosity due to high temporal variation of biopores created by the soil fauna (e.g. earthworms), as observed in our study region (Schneider et al., 2016). Biopores such as earthworm burrows were frequently found to enhance vertical PF (Reck et al., 2018;Weiler and Flühler, 2004;Zehe and Flühler, 2001).

Conclusions 20
Our results demonstrate that infiltration is strongly controlled by PF phenomena. As expected a higher maximum rainfall intensity increases the occurrence of PF, but different from common theory a higher soil moisture decreases the PF occurrence.
However, the here studied landscape units show a high spatial heterogeneity and high temporal variation with different PF processes involved, such as more fast PF in grasslands and more non-uniform flow (NSR) in forest. Clay-rich soils showed to increase both, non-uniform PF (NSR) and fast PF (high vmax). By systematically comparing the dynamics of different landscape 25 units we were able to identify that beside the amount of connected macropores such as cracks (influenced by a high clay content and low soil moisture) or biotic macropores (roots channels, earthworm borrows), PF strongly depends on initiation processes (water repellency, rain intensity). This leads to a strong seasonal dynamics with more non-uniform flow and highest flow velocities in summer and early autumn due to dry soils, high rainfall intensities and hydrophobic soil surfaces.
Furthermore, the amounts of transported water are higher during non-uniform flow. This can have a potential impact solute 30 transport during summer months and should be considered in water management.
We could show that soil texture is not the main driver of water flow velocity during infiltration in the vadose zone as we typically assume. We suggest to include dynamic flow, dynamic initialization processes and varying macroporosity into physically based hydrological models rather than static hydraulic conductivities derived from soil cores or soil maps. Therefore it needs easily transferable relationships or pedotransfer functions, which can help to find structure-related PF parameters similar to retention parameters. More effort is necessary to find or adapt already existing approaches of measuring and 5 monitoring PF in diverse landscapes. We further suggest implementing large-scale sensor networks under different climatic settings, substrates, topographies, and land covers worldwide and to create standardized approaches for analyzing soil moisture datasets. Our approach can be expanded by combining it with groundwater response time series and stable isotope methods to identify and understand flow patterns in the vadose zone at the landscape scale.

Author contribution
DD prepared the data, developed and performed the analysis strategy and planned and conducted the fieldwork. MW and TB designed the sensor cluster setup, were involved in their installation and contributed to the data analysis strategy. DD prepared the manuscript with contributions from all co-authors.