The effect of spatial throughfall patterns on soil moisture patterns at the hillslope scale

Improving the understanding of the controls on subsurface stormflow generation has been the goal of numerous experimental and modeling studies. However, the effect of the spatial variability of throughfall on soil moisture patterns and subsurface stormflow (SSF) generation has not yet been studied in detail. The objectives of this study are threefold: (1) to investigate the influence of a spatially variable throughfall pattern on soil moisture; (2) to investigate if soil moisture patterns reflect a balance between a throughfall and bedrock topography patterns; and (3) to investigate how this balance changes when soil depth, storm size and slope angle are varied. Virtual experiments are used to address these questions. A virtual experiment is a numerical experiment driven by collective field intelligence. It provides a learning tool to investigate the effect of individual processes in a complex system. In our virtual experiment we combined spatial throughfall data from the Huewelerbach catchment in Luxembourg with the topography of a well-studied hillslope within the Panola Mountain Research Watershed, Georgia, USA. We used HYDRUS-3D as a modeling platform. The virtual experiment shows that throughfall patterns influence soil moisture patterns, but only during and shortly after a storm. With a semi-variogram analysis we showed how the effective range of the soil moisture pattern (i.e., the main descriptor of a spatial pattern in case of a small nugget to sill ratio), is similar to the effective range of the throughfall pattern during the storm and gradually returns to the effective range of the bedrock topography after throughfall has ceased. The same analysis was carried out to investigate how this balance changes due to changes in storm size, soil depth, and slope. The analysis showed that the throughfall pattern is more important during large storms on gentle slopes. For steeper slopes the bedrock topography becomes more important.


Introduction
The question of how rainfall makes its way through the canopy to the forest floor has been studied for over a century (DeWalle, 2011).Following the description of interception characteristics of several tree species by Horton (1919), many papers documented interception characteristics for different sites and tree species around the world (Navar et al., 1999;Bruijnzeel and Wiersum, 1987;Calder, 1990;Toba and Ohta, 2005).Whilst interception studies were most commonly done in midlatitude and tropical regions (Bruijnzeel, 2005;Cuartas et al., 2007;Ziegler et al., 2009) interception is also largely determined by climatic parameters.Some studies extended this work into semi-arid regions and snow dominated regions (Storck et al., 2002).
Although the process of forest interception, throughfall initiation, and their influence on point scale infiltration are well documented in the forest hydrology literature, the effect of spatial throughfall patterns on soil moisture patterns is poorly understood (Zehe et al., 2010;Ivanov et al., 2010).This is because the measurement of soil moisture patterns within a forest is much more difficult than the measurement of throughfall.Measuring such patterns jointly is a difficult experimental challenge.Those studies that attempted to explore the effect of throughfall patterns on soil moisture dynamics have been highly equivocal in their findings.For instance, Bouten et al. (1992) found that soil moisture patterns are primarily dependent on soil physical properties, and less on throughfall; on the other hand, Liang et al. (2007) found that topography and stemflow influence soil moisture patterns; Jost et al. (2004) found that mainly the antecedent soil moisture content is important; Raat et al. (2002) concluded that throughfall patterns and forest floor thickness determine soil moisture patterns; Sansoulet et al. (2008) related soil moisture patterns to stemflow and throughfall; and Shachnovich et al. (2008) found that throughfall is not correlated with soil moisture at the point scale.
Paired catchment studies (Bosch and Hewlett, 1982) can be used to understand throughfall influences on soil moisture patterns, but the outcome of paired catchment studies requires enormous data collection and the results would always remain site and case specific.Virtual experiments, however, could overcome these issues.Weiler and McDonnell (2004) defined virtual experiment as "numerical experiment with a model driven by collective field intelligence".A virtual experiment enables us to investigate the response of spatially variable input on hillslope behavior.
Recently, Hopp and McDonnell (2011) investigated the effect of different throughfall patterns on subsurface stormflow (SSF) generation.Their work built upon previous examination of the interactions between slope angle, bedrock permeability, soil depth, and storm size on subsurface stormflow production (Hopp and McDonnell, 2009).Their work used multiple realizations of a fine-scale throughfall pattern and concluded that SSF was controlled by bedrock topography and that throughfall had a limited influence on the amount and timing of SSF generation.The latter might be caused by the fact that throughfall often only influences the near-surface soil moisture patterns (as shown by Liang et al., 2011) and therefore has a limited influence on SSF generation.On the other hand, Liang et al. (2011) also found that high stemflow concentrations influence soil moisture patterns at deeper depths.Since high stemflow concentrations might behave similar as high concentrations of throughfall we do not explore a fine-scale throughfall pattern as Hopp and McDonnell (2011) did, but we explore an observed throughfall pattern with a distinct "hotspot" (as shown in many process studies of throughfall in the field - Germer et al., 2006;Ziegler et al., 2009;Gerrits et al., 2010).Using a larger-scale throughfall pattern can also influence the drainage behavior as hypothesized in Hopp and McDonnell (2011).Furthermore, we only focus on one single throughfall pattern (as selected in the Supplement) and not on differences between throughfall patterns (which was the focus of Hopp and McDonnell, 2011), because Gerrits et al. (2010) showed that the spatial characteristics of throughfall did not differ much throughout the seasons.The tree structure (stem and configuration of the branches) is more dominant than the phenology.We only focus on soil moisture patterns and less on SSF.
We use a geostatistical and visualization approach to quantify relative correlation length scales.
The questions addressed in this paper are: -What is the influence of a certain spatially variable throughfall pattern with clear hotspots on soil moisture?
-Do soil moisture patterns reflect a balance between a certain throughfall and bedrock pattern?
-How does this balance change when soil depth, storm size and slope angle are varied?

Study sites
We used field data from two sites: the Panola hillslope in Georgia, USA, and the Huewelerbach catchment in Luxembourg.The hillslope geometry and soil hydraulic properties were derived from the well-studied Panola hillslope.Storm characteristics were also obtained from Panola; however, the spatial pattern of throughfall to the soil surface was derived from the experimental interception plot in the Huewelerbach catchment in Luxembourg (since this information was not available for the Panola hillslope at the desired temporal scale).

Panola hillslope
The Panola hillslope is located in the Panola Mountain Research Watershed (PMRW) in Georgia Piedmont, 25 km southeast of Atlanta.The site has been described in detail by Freer et al. (2002) and others and here only a brief description is given.The climate is humid and subtropical with an average temperature of 16.3 • C and average rainfall of 1240 mm yr −1 .The hillslope faces southeast and has a slope of 13 • .The surface topography is relatively planar, but the permeable saprolite bedrock (soft disintegrated granite) is highly irregular (Fig. 1a).This results in highly variable soil depths ranging from 0 to 1.86 m, and an average soil depth of 0.63 m.The soil consists of loamy sand with a 0.15 m deep organic-rich horizon.At the lower hillslope boundary a 20 m wide trench has been excavated, where subsurface flow is measured by ten 2 m wide sections.Further details on the Panola hillslope were described by Tromp-van Meerveld andMcDonnell (2006a), andTromp-van Meerveld andMcDonnell (2006b).Tromp-van Meerveld and McDonnell (2006a) presented an analysis of 147 storms at the Panola hillslope and showed that subsurface stormflow only occurred for storm events larger than 55 mm.For this study we selected the best studied rainstorm of the Panola data set of 6-7 March 1996 (Burns et al., 2001;Freer et al., 2002).This storm had a return period of 3 yr and a total storm depth of 87 mm in 31 h divided over two events (see Fig. 1b).It can be seen in Fig. 1b      first rainfall event almost entirely went into storage, while the second event generated the runoff peak.This clearly shows the threshold behavior of the Panola hillslope.

Interception plot Huewelerbach
The interception plot is located in the Huewelerbach catchment in Luxembourg, 20 km northwest of Luxembourg city.The site has been described in detail previously by Gerrits et al. (2010) and only the key details are repeated here.The experimental plot has a total area of 596 m 2 and consists of beech trees (Fagus sylvatica L.) with an average stand density of 168 trees ha −1 .The climate is modified oceanic with mild winters and temperate summers.The average rainfall is 845 mm a −1 and the average temperature is 8 • C (Pfister et al., 2005).In the plot throughfall is measured with 81 manual rainfall collectors installed in a grid with an average distance of about 3 m (Fig. 2a).The collectors are read out at a bi-weekly interval.In an open valley close to the plot gross rainfall is measured by a tipping-bucket rain gauge.Further details on the interception plot can be found in Gerrits et al. (2007) and Gerrits et al. (2010).
From the throughfall data set we selected a random period with full canopy development, because we learned from a time stability analysis that the spatial pattern does not vary much in time (Gerrits et al., 2010).We selected the period starting from the 10 May 2007 until the 25 May 2007 (see Fig. 2b).Total rainfall in this period was 33 mm.
An important hotspot (i.e., point of throughfall funneling) can be seen in Fig. 2b at coordinates (15 m, 15 m) where throughfall exceeds precipitation.This location is around a tree that creates hotspots of high throughfall with lower throughfall values in a ring around the tree.This phenomenon is also observed in field experiments of for example Germer et al. (2006) and Ziegler et al. (2009).
For the analysis we use five classes of throughfall and defined them based on the percentage of maximum throughfall (i.e., 75 mm; Table 1).We used only five classes, because the selected modeling platform HYDRUS-3D (see Sect. 3.1) only has four "variable fluxes" plus one "atmospheric flux" to model spatial input patterns for throughfall.The main difference with respect to the classes as defined by Hopp and McDonnell (2011) is the variance in the throughfall data.The standard deviation of the percentage throughfall of precipitation in Hopp and McDonnell (2011) is 12 % (c v = 0.17), while in this case the standard deviation is 23 % (c v = 0.31).The higher standard deviation in the large-scale pattern is caused by the hotspot.The coefficient of variation, c v of the large-scale throughfall pattern is now in the same order as the c v of the bedrock topography anomaly (i.e., deviation of bedrock topography from a planar slope) of Panola (c v = 0.34).Subsequently, the temporal dynamics of the throughfall are taken from the Panola storm of 6-7 March 1996 (P panola (t)).This storm is scaled for each class i to represent the Huewelerbach spatial pattern: 3 Methods and materials

Model description of base case scenario
To simulate SSF on the Panola hillslope, we used the finite element model HYDRUS-3D, version 1.10 (Simunek et al., 2006).HYDRUS-3D solves the Richards equation for water flow in variably saturated porous media.We used the model as described in detail by Hopp and McDonnell (2009).Here we only briefly describe the crucial information.
The model domain covers an area of 28 m by 48 m.The surface and bedrock topography have been derived from a survey with a spatial resolution of 2 m.From this survey a digital elevation model (DEM) has been generated with a spacing of one meter by linear interpolation; subsequently a mesh of triangular prisms was generated based on this DEM.The mesh consists of 10 layers, each with 1715 nodes.Layers 1 to 5 represent the bedrock sublayer, layers 6 to 10 the soil sublayer.
The model has only been parameterized to the trench outflow on an event basis, and performed in a way consistent with field observations of spatially distributed soil matrix potential (Fig. 1 in Hopp and McDonnell, 2009).The model is parameterized with a uniform rainfall pattern.A spatially distributed pattern would possibly have resulted in a different parameter set, as shown by e.g., Arnaud et al. (2002), Zehe et al. (2005), andTromp-van Meerveld andWeiler (2008).However, the model is not meant to represent the exact behavior of the Panola hillslope, the model is used as a benchmark model for comparison.Important complex processes of Panola like preferential flow (Tromp-van Meerveld and Weiler, 2008) are also not included.We chose not to include macropores in the model description.Soil pipes were observed in the trench face of the Panola hillslope, however, there is no further information on distribution, length and size of soil pipes for the entire hillslope.The study of Hopp and McDonnell (2009) also showed that the simulated distributed response of state variables (pressure head time series) agreed very well with measured pressure head dynamics, indicating that the model could reproduce hydrologic behavior satisfactorily.The model has also not been tested for long-term hydrological modeling.
Soil hydraulic properties are described by the van Genuchten-Mualem model (Van Genuchten, 1980).Because the study of James et al. (2010) showed that including field observed soil core data (Tromp-van Meerveld et al., 2008) did not lead to reasonable model performance, α and n are calibration parameters.The residual water content (θ r ), saturated water content (θ s ), and saturated hydraulic conductivity (K s ) are determined based on long-term field observations (McIntosh et al., 1999;Tromp-van Meerveld and McDonnell, 2006c;Tromp-van Meerveld et al., 2007).The values for the residual and saturated water content are determined by the minimum and maximum observed soil moisture.In Table 2 the soil hydraulic properties are given.We chose to keep the soil hydraulic properties in the horizontal plane homogeneous, however in the vertical direction they differ.For a more detailed explanation on how soil parameters were specified see Hopp and McDonnell (2009).
The boundary conditions of the model domain for the upper and side boundaries are defined as "no flux".The downslope boundary is different for the bedrock and for the soil layers.The downslope bedrock boundary is defined as "no flux", thereby assuming negligible lateral flow in the bedrock.We realize that this is a limitation, because tracer experiments by Tromp-van Meerveld et al. (2007) showed that water can move quickly through the bedrock.However, we chose not to include preferential flow in lateral direction.The downslope boundary of the soil layers is treated as a "seepage flux" allowing water to leave the domain through saturated parts of the boundary.The bottom boundary is assigned as "free drainage", meaning a unit total vertical gradient.For the surface boundary we used generated throughfall patterns as will be described in Sect.3.2 (Fig. 3).As the initial conditions for the entire domain a pressure head of −0.7 m is used, followed by a 7 days drainage period where no rainwater enters the domain.This corresponds to the actual weather conditions that preceded the storm event of 6-7 March 1996.
In this paper we consider the outflow over the entire hillslope width (28 m) and not only the outflow from the excavated trench (20 m) as described by Freer et al. (2002).However we refer to "trench" if we mean the entire downslope boundary for simplicity in writing.The trench is divided atial throughfall pattern ('Upper Right-2.1') on Panola hillslope.Class 1: magreen; class 3: yellow; class 4: dark blue; class 5: light blue.  in 13 segments, where segment 1 equals the outflow from 1-3 m, segment 2 outflow from 3-5 m, segment 3: 5-7, etc. (Fig. 1).

Approach
To investigate the effect of spatially variable throughfall on soil moisture, we selected the throughfall pattern from the Huewelerbach catchment and used this as input to a numerical model of the hillslope.Hopp and McDonnell (2009) already developed a finite element model of the Panola hillslope.We used the same model domain and identical parameters and combined it with the large-scale Huewelerbach throughfall pattern for our virtual experiment.Since the model domain of Panola is larger than the spatial throughfall pattern we needed to expand the throughfall pattern in a way that the spatial characterization remained the same.Since we did not want to enlarge the pattern, we mapped the pattern in eight different ways onto the Panola hillslope.The eight patterns were derived by rotating and mirroring the original pattern of the Huewelerbach.All eight patterns had a total storm size of 63 mm.These patterns were compared with a uniform pattern with also a total storm size of 63 mm.Finally, the pattern with the largest influence on SSF was used for this study.Although this suggests that this pattern may also have a maximum effect on soil moisture, the differences between the patterns on SSF are small.In Fig. 3 the selected pattern ("Upper Right-2.1") is shown.We verified that this pattern has similar geostatistical properties as the initial pattern.Hence we conclude that using the simple method of mirroring (instead of using a geostatistical method) to enlarge the Huewelerbach pattern is in this case justified.See Supplement for details on the creation and selection of the eight throughfall patterns.
Mean soil moisture patterns (depth average) are analyzed with semi-variograms (Cressie, 1993;Chilès and Delfiner, 1999) to investigate if the soil moisture pattern reflects a balance between throughfall and bedrock patterns, since spatial scale is an important controlling factor (Nicótina et al., 2008;Mandapaka et al., 2009).A semi-variogram represents the variance of two points separated by a certain distance in a spatial field and explains to which distance observations are still correlated (Keim et al., 2005): Where h is the lag distance, n(h) is the number of measurement pairs in the data set that are a distance h apart, and N x,y is the anomaly of the spatial data at measuring point (x, y): with σ being the standard deviation.We used normalized data to be able to compare spatial patterns.Important features of a semi-variogram are the nugget, sill, and range.The nugget is a measure for the randomness of observations at one and the same location.The sill is the limit of the semi-variogram, where no autocorrelation exists anymore.The range is a measure of the distance over which there is significant spatial correlation.In case of a small nugget to sill ratio, the range is a good descriptor of a spatial pattern.
We calculated the semi-variograms of the anomaly of the bedrock topography (i.e., deviation of DEM from a plane with slope A), the throughfall, and the soil moisture pattern (average of nodal values of mesh layers 5-10, simulated by HYDRUS).We fitted an exponential model (Chilès and Delfiner, 1999) to find the effective range (r) of the semivariogram, which is the correlation length.
The effective range (r) is defined as the lag (h), where the variance (γ ) is 95 % of the sill (c), and is a measure for the correlation between the points.High spatial correlation between the throughfall collectors causes the effective range to be large and vice versa.In other words, it is the length over which the data points are still spatially correlated with each other.We did not consider anisotropy, hence the effective range is equal to the modulus of the lag.Because we found that the nugget (n) was small and to save computation time (limiting the degree of freedom), we assumed a nugget of zero.Furthermore, we investigate if the balance between throughfall and bedrock patterns changes if storm size (R), slope angle (A) and soil depth (S) are changed according to Table 3.We selected these three characteristics, because Hopp and McDonnell (2009) found that these are important controls for subsurface stormflow generation.We realize that changing the storm size also alters the spatial characteristics of the throughfall pattern due to the fact that interception is a threshold process (Gerrits et al., 2009;Savenije, 2004).Knowledge on how spatial throughfall patterns change during a storm event requires high temporal resolution data on spatial throughfall.Throughfall studies generally do not provide this (e.g., Bouten et al., 1992;Keim et al., 2005;Staelens et al., 2006;Zimmermann et al., 2008;Gerrits et al., 2010), hence we chose to neglect this effect in our study.By changing the slope angle we also have to adjust the throughfall pattern because trees grow vertically and not perpendicular  to the hillslope.For steeper slopes the same number of trees would result in a denser canopy coverage and thus a shorter effective range.Therefore, we adjusted the effective range (r ) of the throughfall pattern by multiplying the effective range of the base case (r (base case) ) with the cosine of the slope angle, and took the effective range of the base case scenario as a reference: 4 Results

Throughfall effects on SSF
In Fig. 4a subsurface stormflow along the downslope trench of the base case scenario (R = 63 mm, A = 13 • , S = 0.62 m) with uniform input is shown.The shape of modeled total subsurface flow is similar to the observed hydrograph at Panola (see Fig. 1).As can be seen in Fig. 4a, subsurface flow varies strongly along the trench (variance Q s /Q t = 10.4 × 10 −2 ); segment 6 drains the major part of the hillslope.This segment is on the transition of the very shallow soil to the thicker soil and has a relatively large contributing upslope surface area.In Fig. 4b the modeled subsurface flow response of the distributed input ("Upper Right-2.1") is shown.The hydrograph of this pattern is significantly different from the uniform pattern mainly caused by a different response from segments 6 and 7.While the uniform pattern has a rather smooth recession curve, pattern "Upper Right-2.1"leads to a double peak for segments 6 and 7. Modeled SSF for all other segment hydrographs do not differ much from the uniform case.

Throughfall effects on soil moisture
In Fig. 5 the soil moisture pattern of distributed input is shown for the five soil mesh layers (mesh layers 10-6) and the bedrock interface layer (mesh layer 5).The mesh layers are equally distributed between the surface layer and the bedrock.Thus the soil moisture patterns do not represent the soil water content at the same depth, which makes it more complicated to compare the results with field observations or to interpret the patterns.The soil moisture patterns are plotted for t = 5 h (first storm peak), t = 23 h (just before second storm peak), t = 25 h (second storm peak), t = 56 h (highest deviation in SSF between uniform and distributed case), and t = 190 h (end of drainage).
The patterns show that at the depressions (i.e., large soil depth) in the bedrock (see Fig. 1a for location of bedrock depressions) the soil moisture content is low whereas it is high in the shallow soil (i.e., small soil depth) in the near-surface layers.The shape of the bedrock and the main drainage channels are clearly visible at t = 56 h at the bedrock interface layer (mesh layer 5).As can be seen at t = 5 h (during the first peak of the storm), only the top near-surface layers are influenced by the throughfall pattern.The four hotspots are clearly visible.During the second peak (t = 25 h) the hotspots are again visible in the top layers, but now also at the deeper soil layers.

Soil moisture comparison of uniform and distributed input patterns
The mean soil moisture content over mesh layers 5-10 does not differ much between uniform input and distributed input (Fig. 6).Only at t = 56 h there is a significant difference between the two input patterns, which lasts for about a day, similar to the modeled SSF (Fig. 4).Also the maximum mean soil moisture content is higher just after the second storm peak.The spatial differences per mesh layer are shown in Fig. 7.At the beginning of the storm the distributed input pattern only influences on the top layers and slowly moves down in time.During the second peak of the storm the effect of the distributed pattern is largest, because the effect of the previous storm peak has not yet vanished.Also the size of the hotspots in the soil moisture pattern are larger than at t = 5 h.After the storm, the difference between uniform and distributed become negligible except for the shallow thin soils in the lower left corner.At t = 190 h also this difference disappears.

Spatial correlation of soil moisture content
From the model results it appears that the soil moisture content (e.g., for the soil-bedrock interface: mesh layer 5) is highly correlated to the bedrock topography (Fig. 8a) when it has been dry for a long time (Fig. 8c) and that during a large rainfall event, or shortly after, the soil moisture pattern represents the rainfall pattern (see Fig. 8b).
In Fig. 9a the semi-variogram of normalized throughfall is shown; and as can be seen the sill is not constant.The semivariogram of the average soil moisture pattern for different time steps (Fig. 9b) appears to change between the semivariogram of the throughfall pattern (Fig. 9a) and the bedrock topography (Fig. 9c).To test this hypothesis, we choose to look at the effective range of the semi-variogram, r, as the main characteristic of the spatial pattern.This is valid if the nugget to sill ratio is small.For the throughfall we found an effective range of 5 m and for the bedrock anomaly 20 m.We hypothesize that between rainfall events the soil moisture pattern has similar spatial characteristics as the topography,    but during a rainfall event the pattern becomes more similar to the throughfall pattern.
In Fig. 10 this can indeed be seen for the base case scenario.The effective range of the average soil moisture (Fig. 10a) starts at 13 m and drops to 11 m during the first rainfall event.After rainfall ceases, the effective range returns back in the direction of the effective range of the topography.When the second rainfall starts, the effective range of the soil moisture pattern again drops towards the effective range of the throughfall pattern.And finally, after the rain stops, it again moves back in the direction of the effective range of the bedrock anomaly.Hence the change of the effective soil moisture range (blue line in Fig. 10a) looks similar to a hydrograph, and can be called a "geostatistical hydrograph".The fact that the effective start and end range differ, is a consequence of the (nonrealistic) initial soil moisture pattern.
In Figs.10b and c the same plot is generated for the shallow soil layers and the deep soil layers, respectively.As can be seen, the first event has a larger effect on the shallow layers than on the deeper layers.On the other hand, the second event affects both the shallow and the deep soil layers.The effect of drainage is also visible.The peak in the effective range starts earlier for the shallow soil layers than for the deeper soil layers.
If the geostatistical hydrograph is generated for all combinations of storm size, slope angle, and soil depth, we can investigate if the hillslope attributes change the spatial pattern of the soil moisture.In Fig. 11 we present the interplay between storm size, slope angle, and soil depth on the coefficient of determination (R 2 ; Fig. 11a), the equilibrium soil moisture range (Fig. 11b), the effective range of the second peak (Fig. 11c), and the time to the second peak (Fig. 11d).For visibility reasons we only show slices through the 3dimensional cubes.To compensate for the fact that trees grow  2) for a) throughfall pattern, b) the profile average soil moisture pattern at time steps t =0h (start), t =5h (during first storm peak), t =23h (just before second storm peak), t =25h (second storm peak), t =56h (highest deviation in SSF between uniform and distributed case), and t =190h (end of drainage), and c) the bedrock anomaly pattern.
38 Fig. 9. Semi-variograms (Eq.2) for (a) throughfall pattern, (b) the profile average soil moisture pattern at time steps t = 0 h (start), t = 5 h (during first storm peak), t = 23 h (just before second storm peak), t = 25 h (second storm peak), t = 56 h (highest deviation in SSF between uniform and distributed case), and t = 190 h (end of drainage), and (c) the bedrock anomaly pattern.
vertically and not perpendicular to the hillslope, we corrected the effective range of the throughfall by multiplying it with the cosine of the slope and standardized it to the base case slope.Hence the effective range of the throughfall pattern is 5.3, 5.2, 4.8, and 4.0 m for angles of 6.5, 13, 26, and 40 • , respectively.
First, we looked at the mean performance of the fitted semi-variogram model (Eq.4) for mesh layers 5-10 and t > 0. In Fig. 11a the model performance is expressed by the coefficient of determination (defined as the ratio of the sum of squares of the regression and the total sum of squares).Although the overall performance is good with a mean R 2 of 0.8 ± 0.1, it appears that the R 2 is not as high for the steep slopes.For steep slopes, the soil moisture pattern has mostly slope-parallel, straight flow paths just after the rainfall event, which cannot be described well with an exponential model due to anisotropy, causing the relatively bad performance of the model.
Second, we looked at the equilibrium state.This is the final effective range of the soil moisture at t = 190 h and is shown in Fig. 11b for all cases.The final effective range is scaled between the effective range of the topography (20.4 m) and the throughfall (5.3, 5.2, 4.8, and 4.0 m for angles of 6.5, 13, 26, and 40 • , respectively).The interpolated cube shows that with increasing slope the final effective range becomes larger.Hence with increasing slope the topography becomes more important.Furthermore, there is a slight increase in final effective range with increasing storm size.This was also observed by Zehe et al. (2010), who found that with increasing storm size the precipitation patterns exert a more dominant control on soil moisture.Also Western et al. (2004) observed this and concluded that this is likely inherent to the method, where spatial patterns of rainfall are analyzed with the effective range.
Third, we analyzed the effective range of the second peak of the geostatistical hydrograph in Fig. 11c.Similar to Fig. 11b, we scaled the effective range between the effective range of the topography and the throughfall.As can be seen, the effective peak range is related to storm size and slope angle.The bigger the storm the more the throughfall pattern influences the modeled soil moisture pattern and the steeper the slope, the more the bedrock topography influences the modeled soil moisture pattern.This is because the high gradient drains the rainwater so quickly that the throughfall pattern only remains for a short period.The soil depth appears not to have any influence as long as the soil layer is thick enough.Only for very thin soil layers does the throughfall pattern have a larger influence on the modeled soil moisture pattern.
In Fig. 11d the interplay of the hillslope attributes and the time to peak is shown.The time to peak is defined as the time between the start of the rain and the peak of the effective range in the soil moisture pattern.For deep soils the peak occurs faster with increasing storm size and slope angle; however, for a mean soil depth of 1.22 m it appears that slope angle and storm size do not have an influence.For the very thin soil layers the pattern is similar to the thick soil layers.However, for gentle slopes and small storms the modeled time to peak is short.

Discussion
Our modeling results show that a nonuniform throughfall pattern affects SSF.Although the spatial pattern does not have a significant effect on total SSF, it does influence the distribution along the trench and more importantly the shape of the hydrograph.Hopp and McDonnell (2011) also investigated the effect of throughfall on SSF; however, they did not find a clear impact.This is likely caused by the throughfall pattern; Hopp and McDonnell (2011) used a finescale throughfall pattern, while in this study a pattern with a distinct hotspot is used.Depending on the location of the hotspot the throughfall pattern influences modeled SSF.If the hotspot is located above a "channel" of high flow accumulation, this causes quick drainage.This is likely related to the connectivity of the saturated areas as found by Tromp-van Meerveld and McDonnell (2005) to be the main cause for SSF on many hillslopes.Tromp-van Meerveld and McDonnell (2005) found that SSF is strongly correlated to hillslope average soil moisture; however, also that there is a lack of correlation between the (near surface) soil moisture pattern and subsurface saturation due to soil depth and bedrock topography.They found that the soil moisture pattern at the soil-bedrock interface is most important for subsurface saturation.In our study, we found   by means of the geostatistical hydrograph how the soil moisture pattern is influenced by the bedrock topography and the throughfall pattern during a storm.We showed that before the simulated storm the soil moisture pattern reflects the bedrock topography and that during the storm the througfall hotspots are clearly visible.This is both the case for the near-surface layers and the deeper layers.Since our study neglects the effect of transpiration, macropores, spatial variability in soil organic matter, bulk density, soil texture, and antecedent soil moisture, it is possible that in reality this relation is not always so evident.However, similar to our study, other studies that considered hotspots (Liang et al., 2007;Sansoulet et al., 2008;Jost et al., 2004;Raat et al., 2002) also found that hotspots have a big influence on the soil moisture distribution during and shortly after an event.Liang et al. (2007) and Sansoulet et al. (2008) considered in their study only shallow soil moisture, but Jost et al. (2004) and Raat et al. (2002) also looked at deep (profile average) soil moisture.Shachnovich et al. (2008) did not find this relation.However, they compared the change in soil moisture over a week, while the dynamics of throughfall patterns has a much shorter timescale.
During dry periods, the bedrock topography becomes more important as shown in the geostatistical hydrograph and thus soil physical properties are important, as also found by Bouten et al. (1992).Also Tromp-van Meerveld and McDonnell (2006c) found that during the dry state, soil moisture was correlated (although not much) to bedrock topography at the Panola hillslope.The relatively low correlation between the bedrock and the observed soil moisture patterns compared to our modeled results are likely due to transpiration (the study of Tromp-van Meerveld and McDonnell, 2006c, took place in summer), which is not part of our analysis.In reality, the soil moisture pattern will not represent the bedrock topography after a long time, since processes like transpiration and macropores will change the soil moisture pattern as well.Western et al. (1999) and Grayson et al. (1997) also found that soil moisture patterns were influenced by topography, although they considered a seasonal timescale and not the event scale.They found that during the wet state (winter period) the soil moisture pattern was dominated by lateral flow and thus was organized according to the topography.During the dry state (summer period) vertical water flow was dominant and hence the soil moisture was less organized by the topography.Our study is comparable to a "wet season".
The geostatistical-hydrograph analysis may help to understand and predict soil moisture patterns based on throughfall and bedrock patterns; however, one should be careful with using the effective range as the descriptor of a spatial pattern.Two completely different spatial patterns can have the same effective range.Furthermore, in this study we only looked at the balance between throughfall patterns and bedrock topography.However, in reality, other processes also influence soil moisture patterns.For example, our study neglects transpiration, but also its related feedback mechanisms, which can be important.Bouten et al. (1992), for example, stated that the short influence of throughfall on soil moisture is due to feedback mechanisms of drainage and water uptake.Trees optimize their root system to water availability, causing higher water uptake where throughfall is high ("preferential uptake").At Panola most trees indeed grow in the "drainage channels" where the soil depth is large (see Fig. 12 in Trompvan Meerveld and McDonnell, 2006c).However, we think that where throughfall hotspots exist, throughfall will -on the short timescale, and during wintertime-be dominant over water uptake.This is also shown by Liang et al. (2007) and Sansoulet et al. (2008) who did not study hotspots by througfall, but focused on hotspots caused by stemflow, which is in principle the same.The only difference between stemflow and throughfall hotspots is the influence of roots.As shown by Liang et al. (2007), bypass flow occurs especially close to the stem of the trees.For throughfall hotspots this might be less important, but antecedent soil moisture conditions become more important.When it has been dry for a long time, cracks can develop, which facilitate macropore flow (Jost et al., 2004).Unfortunately, the effect of macropore flow has not been considered in our modeling study, because of the lack of data to properly model this.Macropores could lead to a quicker drainage of the throughfall, resulting in a shorter duration of impact of the throughfall pattern on the soil moisture pattern.As a result the saturation pattern at the bedrock interface will be different (Tromp-van Meerveld and Weiler, 2008).However, even if macropores exist the throughfall pattern will influence the soil moisture pattern.This should be investigated in the future.
The influence of additional drivers like macropore flow and root water uptake can in principle also be investigated with the presented geostatistical-hydrograph analysis.Instead of a balance between two drivers (throughfall and bedrock topography), the analysis can be extended to a balance between multiple drivers (e.g., macropore flow and root water uptake).

Conceptual framework
As stated in the introduction, the literature is equivocal on the importance of throughfall on soil moisture content.Although their results may be contradictory, they may not be so different if one considers the conditions under which the experiments were carried out.Throughfall conspires with other factors, such as slope angle or storm size; the interactions between factors shape hillslope hydrological responses (Hopp and McDonnell, 2009).
We extended the conceptual model presented by Hopp and McDonnell (2009) by looking at the effect of spatially variable throughfall input on modeled soil moisture.In Fig. 12 the outline of our conceptual model is shown.We divided the model into three parts: (1) fixed hillslope configuration; (2) time variable input; and (3) the hillslope response (or state).In part 2 we added the interception threshold, which causes spatially variable throughfall patterns.These throughfall patterns influence the soil moisture content, depending on the slope and the storm size.On a relatively flat hillslope the impact of throughfall on soil moisture is larger in terms of the magnitude of the effective range (Fig. 11c) and duration before it bounces back to the bedrock pattern (Fig. 11b), because lateral drainage along the soil-bedrock interface is slow due to the small slope angle.Storm size mainly influences the magnitude of the effective range: the bigger the storm, the more the soil moisture pattern reflects the throughfall pattern (Fig. 11c).However, in reality the spatial pattern will be affected by the size of the storm, due to intensity smoothing of the canopy.Hence, not only the storm size, but also the intensity is of importance (Keim et al., 2006).
In case hotspots exist (either for throughfall or stemflow), the spatial throughfall pattern can also influence SSF (Fig. 4).It really matters where the hotspot is located in relation to the soil depth distribution and the bedrock topography.First of all, the shape of the hydrograph can be different, but also the distribution along the trench is influenced by hotspots.
As concluded by Hopp and McDonnell (2011), throughfall is not a first order control for SSF in comparison to slope angle and storm size.Keim et al. (2006) also found that storm size has a limited effect on SSF, however, they also concluded that next to storm size also the storm intensity is important.The effect of intensity (and intensity smoothing) is not considered in this study.Our study also shows that throughfall is indeed not a first order control; however, the importance of throughfall patterns on SSF highly depends on the spatial variability of the throughfall pattern in relation to the topography of the bedrock.
In our virtual experiment, we did not include the influence of water uptake by plants.Transpiration will likely reduce the effect of throughfall patterns on soil moisture patterns, because roots optimize their root system to water availability (Bouten et al., 1992).Also antecedent soil moisture conditions and the influence of macropores have not been taken into account, while they do have an influence as shown by Jost et al. (2004), Liang et al. (2007), and Sarkar and Dutta (2012).The difficulty with including water uptake and macropore flow is the mutual interdependence of these processes and the possible feedback between them.Furthermore, macropores are in general difficult to implement in hydrological models and require detailed information on the macropore distribution.Future work can focus on these feedback mechanisms.

Conclusions
The virtual experiments show that throughfall patterns influence modeled soil moisture patterns, but only during and shortly after a storm event.By means of a geostatistical analysis we investigated if the soil moisture pattern reflects the spatial throughfall pattern or that of the bedrock.As an indicator we used the effective range of the semi-variogram.We found that during a storm the soil moisture pattern has a similar effective range as the throughfall pattern, but gradually returns to the effective range of the bedrock pattern after throughfall has ceased.Furthermore, we looked at the impact of hillslope controls and storm size on the geostatistical analysis.It appeared that the throughfall pattern is most important during large storms on gentle slopes.For steeper slopes the bedrock topography becomes more important.The mean soil depth appears to have no significant impact.These findings have been included in a conceptual model, which can be used to evaluate the effect of spatially varying throughfall on SSF and soil moisture patterns.
Overall, we can conclude that interception has a minor influence on SSF generation and a larger impact on the modeled soil moisture patterns during and shortly after rain events.Geostatistical analysis can help to understand the relationship between soil moisture patterns, throughfall patterns and subsurface characteristics.However, more research is necessary to investigate other hillslope variables, such as antecedent wetness, macroporosity, rainfall intensity (Keim et al., 2006), soil evaporation, and transpiration.A next research step would be to confront our model results with observations in an experimental setup.However, field data of high spatial and temporal resolution are rare, and therefore our virtual experiments remain valuable.
Fig. 2. a) Interception plot Huewelerbach: location of the 81 rain gauges to measure throughfall.b) Spatial throughfall pattern (interpolated with a triangle-based cubic interpolation technique and a 10 cm mesh) of May 2007.The black contour lines represent the five throughfall classes, and the white grid the selected area for the analysis.The hotspot is located near the tree at coordinates 15m, 15m.31 Fig. 2. (a) Interception plot Huewelerbach: location of the 81 rain gauges to measure throughfall.(b) Spatial throughfall pattern (interpolated with a triangle-based cubic interpolation technique and a 10 cm mesh) of May 2007.The black contour lines represent the five throughfall classes, and the white grid the selected area for the analysis.The hotspot is located near the tree at coordinates 15 m, 15 m.

Fig. 4 .
Fig. 4. Subsurface storm flow for the entire width of the hillslope (28 m).The upper graphs show the hydrographs of the 13 segments along the trench (Q s ), the lower left the total outflow (Q t ) and the lower right the variability along the trench.The bar indicates the average flow in time.(a) Subsurface storm flow of the base case scenario with uniform input; (b) subsurface storm flow of the base case scenario with spatially variable input "Upper Right-2.1"(see Fig.3).
Soil moisture patterns[-]  of distributed input for soil mesh layers 10 (top) to bedrock interface) at time steps t =5h (during first storm peak), t =23h (just before seco torm peak), t =25h (second storm peak), t =56h (highest deviation in SSF between u orm and distributed case), and t =190h (end of drainage).The x and y axis are the late nd upslope extend of the hillslope.

Fig. 5 .
Fig. 5. Soil moisture patterns (-) of distributed input for soil mesh layers 10 (top) to 5 (bedrock interface) at time steps t = 5 h (during first storm peak), t = 23 h (just before second storm peak), t = 25 h (second storm peak), t = 56 h (highest deviation in SSF between uniform and distributed case), and t = 190 h (end of drainage).The x-and y-axis are the lateral and upslope extent of the hillslope.

Fig. 6 .Fig. 6 .
Fig. 6.Profile average soil moisture content over mesh layers 5-10 in time for uniform input and distributed input.The dashed lines indicate the minimum and maximum of the profile average soil moisture content over mesh layers 5-10.35 Fig. 6.Profile average soil moisture content over mesh layers 5-10 in time, for uniform input and distributed input.The dashed lines indicate the minimum and maximum of the profile average soil moisture content over mesh layers 5-10.
Difference in profile average soil moisture content between uniform input and d ributed input at time steps t =5h (during first storm peak), t =23h (just before second stor eak), t =25h (second storm peak), t =56h (highest deviation in SSF between uniform a istributed case), and t =190h (end of drainage).The x and y axis are the lateral a pslope extend of the hillslope.

Fig. 7 .
Fig. 7. Difference in profile average soil moisture content between uniform input and distributed input at time steps t = 5 h (during first storm peak), t = 23 h (just before second storm peak), t = 25 h (second storm peak), t = 56 h (highest deviation in SSF between uniform and distributed case), and t = 190 h (end of drainage).The x-and y-axis are the lateral and upslope extent of the hillslope.
Fig. 8. a) Flow accumulation map and location of high intensity throughfall input; b) shortly after a rain storm (t=8 hour) at mesh layer 5 (i.e.bedrock interface layer) c) soil moisture pattern of the soil-bedrock interface after a long dry period (t=190 hour).

Fig. 8 .
Fig. 8. (a) Flow accumulation map and location of high intensity throughfall input; (b) shortly after a rain storm (t = 8 h) at mesh layer 5 (i.e., bedrock interface layer); (c) soil moisture pattern of the soil-bedrock interface after a long dry period (t = 190 h).

Fig. 11 .
Fig. 11.a) Interplay of hillslope attributes on model performance R 2 (average for mesh layers 5-10 and t >0h); b) Interplay of hillslope attributes on equilibrium effective range scaled between the effective range of the topography and the throughfall; c) Interplay of hillslope attributes on peak value of the mean effective range of the soil layers scaled between the effective range of the topography and the throughfall.Blue indicates that the effective range is very similar to the effective range of the throughfall pattern, and red indicates similarity to the effective range of the bedrock pattern; d) Interplay of hillslope attributes on the time to peak (i.e., time between start of rain and peak in effective range). 40

Fig. 11 .
Fig. 11.(a) Interplay of hillslope attributes on model performance R 2 (average for mesh layers 5-10 and t > 0 h); (b) interplay of hillslope attributes on equilibrium effective range scaled between the effective range of the topography and the throughfall; (c) interplay of hillslope attributes on peak value of the mean effective range of the soil layers scaled between the effective range of the topography and the throughfall.Blue indicates that the effective range is very similar to the effective range of the throughfall pattern, and red indicates similarity to the effective range of the bedrock pattern; (d) interplay of hillslope attributes on the time to peak (i.e., time between start of rain and peak in effective range).

Table 1 .
Definition of the five classes and the main characteristics, with T f throughfall sum, P gross precipitation sum (=33 mm), and Tf the mean throughfall sum in a class.The spatial cover percentage is only calculated for the initial spatial throughfall pattern.

Table 2 .
Soil hydraulic properties of the van Genuchten-Mualem model.Asterisks indicate the base-case scenario with mean soil depth S = 0.62 m.

Table 3 .
Variations of input and topography.Asterisks indicate the base case scenario.