Throughfall spatial patterns translate into spatial patterns of soil moisture dynamics – empirical evidence

Throughfall heterogeneity induced by the redistribution of precipitation in vegetation canopies has repeatedly been hypothesized to affect the variation in the soil water content and runoff behavior, especially in forests. However, observational studies relating the spatial variation in the soil water content directly to net precipitation are rare, and few confirm modeling hypotheses. Here, we investigate whether throughfall patterns affect the spatial heterogeneity in the soil water response in the main rooting zone. We assessed rainfall, throughfall and soil water content (at two depths, 7.5 and 27.5 cm) on a 1 ha temperate mixed-beech forest plot in Germany during the 2015 and 2016 growing seasons using independent, high-resolution, stratified, random designs. Because the throughfall and soil water content cannot be measured at the same location, we used kriging to derive the throughfall values at the locations where the soil water content was measured. We first explored the spatial variation and temporal stability of throughfall and soil water patterns and subsequently evaluated the effects of input (throughfall), soil properties (field capacity and macroporosity), and vegetation parameters (canopy cover and distance to the next tree) on the soil water content and dynamics.

Throughfall spatial patterns were related to canopy density. Although spatial autocorrelation decreased with increasing event sizes, temporally stable throughfall patterns emerged, leading to reoccurring higher- and lower-input locations across precipitation events. Linear mixed-effects model analysis showed that soil water content patterns were poorly related to spatial patterns of throughfall and that they were more influenced by unidentified, but time constant, factors.

Instead of the soil water content itself, the patterns of its increase after rainfall corresponded more closely to throughfall patterns: more water was stored in the soil in areas where throughfall was elevated. Furthermore, soil moisture patterns themselves affected the local soil water response, probably by mediating fast drainage and runoff. Locations with a low topsoil water content tended to store less of the input water, indicating preferential flow. In contrast, locations with a high water content stored less water in the subsoil. Moreover, the distance to the next tree and macroporosity modified how much water was retained in soil storage. Overall, throughfall patterns imprinted less on the soil water content and more on the soil water dynamics shortly after rainfall events; therefore, percolation rather than the soil water content may depend on the small-scale spatial heterogeneity in canopy input patterns.

1 Introduction

Over the past decades, in order to understand and predict hydrological processes related to vegetation structure and land use change, there has been an increased interest in how water input at the soil surface is affected by vegetation canopies (Western et al., 2004; Savenije, 2004; Murray, 2014; Guswa et al., 2020; Oda et al., 2021). Due to interception losses, less water arrives below the canopy compared with above (Horton, 1919, and references therein; Carlyle-Moses and Gash, 2011), which has implications for the soil water balance (Durocher, 1990; Bouten et al., 1992; Schume et al., 2003; Klos et al., 2014; Metzger et al., 2017) and the overall water budget at the catchment scale (Brown et al., 2005; Oda et al., 2021).

Besides interception loss, the interaction of precipitation with the vegetation canopy causes the spatial redistribution of incoming water. This leads to characteristic spatial heterogeneity in the dripping (throughfall) and flowing (stemflow) below-canopy precipitation, causing enhanced water input to the soil surface at the local scale. For example, dripping point hot spots (enhanced water flow from peculiarities in the canopy; Falkengren-Grerup, 1989; Keim et al., 2005; Staelens et al., 2006; Voss et al., 2016) and stemflow hot spots (Levia and Germer, 2015; Carlyle-Moses et al., 2018) have been well documented. The available research suggests that both throughfall patterns and stemflow spatial distributions are reoccurring (Keim et al., 2005; Staelens et al., 2006; Zimmermann et al., 2008; Wullaert et al., 2009; Guswa and Spence, 2012; Metzger et al., 2017; Van Stan et al., 2020).

The observed persistence of spatial patterns of below-canopy precipitation has created a strong expectation that the aforementioned patterns affect soil water content patterns (Schume et al., 2003; Wullaert et al., 2009; Rosenbaum et al., 2012; Zehe et al., 2010) and hot spots of percolation or preferential flow (Bouten et al., 1992; Schume et al., 2003; Blume et al., 2009; Bachmair et al., 2012) in forests soils. However, this has only been partly confirmed with observations: for stemflow-affected locations, soil moisture microsites have repeatedly been demonstrated (Pressland, 1976; Durocher, 1990; Liang et al., 2007; Germer, 2013; Metzger et al., 2021). Stemflow can create a substantial funneling of water to the forest floor, and water availability on the forest floor can be locally enhanced by a factor of 10–100 (Levia and Germer, 2015; Carlyle-Moses et al., 2018; Metzger et al., 2021).

Much less research is available regarding how the less pronounced but still spatially persistent pattern of throughfall shapes soil water dynamics. Modeling has suggested that throughfall patterns influence the root zone soil moisture pattern (Coenders-Gerrits et al., 2013; Guswa, 2012). However, soil moisture patterns are also influenced by several other factors that create substantial heterogeneity, such as the heterogeneity in soil properties, local micro-topography, litter thickness or root water uptake (Bouten et al., 1992; Schume et al., 2003; Schwärzel et al., 2009; Gerrits and Savenije, 2011; Rosenbaum et al., 2012; Liang et al., 2017; Molina et al., 2019), and these factors are typically not fully captured in virtual experiments. In comparison to modeling work, observational studies have found that throughfall and root zone soil moisture are not related (Shachnovich et al., 2008; Rodrigues et al., 2022), only occasionally related (Metzger et al., 2017) or weakly related (Molina et al., 2019). In contrast, Klos et al. (2014) found a relation below the rooting zone by strategically sampling at high and low throughfall positions, and several authors (Zehe et al., 2010; Rosenbaum et al., 2012; Metzger et al., 2017) have found indirect evidence by interpreting the change in the spatial variation in the soil water content after precipitation events.

In light of the substantial heterogeneity in other influencing factors, one of the reasons for the limited direct observational evidence of the effect of throughfall on the soil water content may be the limited number of studies investigating the relation between below-canopy precipitation and soil water patterns in a dedicated and coordinated fashion. The characterization of spatial patterns, such as those of throughfall, requires a large number of throughfall samplers (Kimmins, 1973; Lloyd and Marques, 1988; Zimmermann et al., 2010; Van Stan et al., 2020), and the same is true for belowground observations. Furthermore, a fundamental challenge is that the soil water input and soil water content cannot be assessed at the same location, as throughfall measurements disturb infiltration into the soil.

Therefore, the objective of this study is to compare the patterns of soil water content, soil properties and throughfall using a dedicated high-spatial-resolution sampling design to reveal whether the input, distance to the next tree or soil properties affect the spatial variation in the soil water content and soil water response. We used independent designs for above- and belowground observations and applied kriging to derive the throughfall values at the locations where the soil water content was measured. The aims of the study were to (a) explore the spatial heterogeneity and temporal stability in throughfall and the soil water content and (b) evaluate the influence of soil properties (field capacity and macroporosity), vegetation parameters (canopy cover and distance to the next tree) and input variation (throughfall) on the variation in the soil water content and soil water content increase after precipitation.

2 Methods

2.1 Study area

The study was carried out in the Hainich Critical Zone Exploratory (CZE Hainich; see Küsel et al., 2016), run by the Collaborative Research Centre “AquaDiva”. The site is located in an unmanaged beech-dominated forest in the Hainich National Park, central Germany. The mean annual temperature is around 7.5–9.5 ^∘C, depending on the position on the small mountain range. The total annual precipitation drops from 900 to less than 600 mm from ridge to valley (Küsel et al., 2016). The monitoring site as well as measurements of precipitation and soil moisture have been described in Metzger et al. (2017), and the important information is repeated here for completeness. The site covers an area of 1 ha and is situated at 365 m above sea level. The study area contains of 581 individual trees (diameter at breast height of ≥5 cm), representing a heterogeneous age structure. The soils in this area are dominantly Luvisols (Schrumpf et al., 2014; Kohlhepp et al., 2017). The species assemblage consists of 70 % European beech trees (Fagus sylvatica) and smaller proportions of sycamore maple (Acer pseudoplatanus), European ash (Fraxinus excelsior), European hornbeam (Carpinus betulus), large-leaved linden (Tilia platyphyllos), Norway maple (Acer platanoides) and Scots elm (Ulmus glabra). The weathered bedrock is at a depth of 15–87 cm (median depth 37 cm).

2.2 Precipitation measurements and processing

The precipitation sampling followed the procedure given in Metzger et al. (2017). Gross precipitation (P_g) and throughfall (P_TF) were measured manually on a per-event basis in spring 2014, 2015 and 2016. The analysis of the effects of throughfall on the soil water content covers the period when soil moisture sensors were also active, i.e., from 18 June to 28 July 2015 and from 31 May to 14 July 2016. The installed throughfall collectors consisted of circular funnels (A=0.011 m²), the opening of which was placed about 37 cm above the ground surface. A table tennis ball was placed in the opening of the funnel to minimize evaporation.

Throughfall collectors were arranged in a stratified sampling design (Zimmermann et al., 2016). For this, the 1 ha plot was divided into 100 subplots (Fig. 1), each 10 m×10 m, and equipped with two randomly located throughfall samplers. Of those samplers, we selected 50 randomly and added another sampler in their immediate vicinity (0.1 m distance), thereby creating a “short transect”. Furthermore, we added four more samplers to 25 randomly selected short transects at 0.5, 1, 2 and 3 m from the first sampler of the short transect, thereby forming “long transects”. The direction of all transects was also randomly chosen. In total, we sampled n=350 throughfall positions. Sampling started 2 h after the end of rainfall and was carried out by collecting the volume of all sampling containers using graduated cylinders. Gross precipitation was measured at an adjacent (distance 250 m) open grassland using five funnels of the same type as the throughfall collectors.

Figure 1 Experimental setup showing the 1 ha forest plot subdivided into a 10 m×10 m grid, yielding 100 subplots, and the positions of the throughfall samplers (pink crosses) and 49 soil water content subplots (blue) measured in a stratified random design with transects (see Sect. 2 for more details). The figure was sourced from Metzger et al. (2017).

In this research, we are specifically interested in re-occurring spatial patterns of throughfall and whether they translate to soil water dynamics. To characterize the spatial pattern, we decompose the measured throughfall into the spatial median (${\widehat{P}}_{\text{TF},j}$) of event j and deviations from the median at the specific location i (δP_TF,ij). δP_TF,ij is calculated as follows:

$$\begin{array}{}\text{(1)}& \mathit{\delta }{P}_{\text{TF},ij}={\displaystyle \frac{P_{\text{TF},ij}-{\widehat{P}}_{\text{TF},j}}{{\widehat{P}}_{\text{TF},j}}},\end{array}$$

where δP_TF,ij represents the relative deviations of the spatially distributed throughfall (P_TF,ij) at location i and event j from the spatial median ${\widehat{P}}_{\text{TF},j}$ of that event j. Equation (1) is a slight modification of the widely used concept of temporal stability introduced by Vachaud et al. (1985). As throughfall can contain outliers, we used the median (${\widehat{P}}_{\text{TF}}$) instead of the arithmetic mean for normalization, as already done by Zimmermann et al. (2008) and Wullaert et al. (2009a). Negative or positive values of δP_TF,ij indicate locations with below- or above-average throughfall, respectively, while repeatedly low or high δP_TF,ij values indicate persistent cold or hot spots of canopy throughfall across events, respectively. Equation (1) allows one to disentangle the temporal variation, e.g., the event size given by the event spatial median ${\widehat{P}}_{\text{TF},j}$, from the spatial dispersion, characterized by δP_TF,ij. Both are independent of each other and both are used in parallel as predictors of the soil water content and soil water content increase below. In the following, we omit the indices i and j for simplicity and refer to δP_TF as the spatial pattern of throughfall for an event or event class.

To investigate the temporal persistence of the spatial pattern of throughfall, we derived temporal stability plots (Zimmermann et al., 2008; Wullaert et al., 2009) by looking at δP_TF,ij across events at the specific location i. Additionally, we calculated the Spearman rank correlation coefficients between observations of different events, where high correlations indicate strong persistence (or temporal stability) of the throughfall pattern. We paired all events falling into a given rain event class according to Metzger et al. (2017): small (P_g≤3 mm), medium ($\mathrm{3}\mathrm{mm}<P_{\mathrm{g}}\le \mathrm{10}$ mm) and large (P_g>10 mm).

To relate the general precipitation and soil moisture conditions during the observation period to the average climate, we compared them with precipitation data from a nearby weather station (Mühlhausen-Windeberg, 20 km to the northeast) of the German Weather Service (DWD Climate Data Center, https://www.dwd.de/cdc, last access: 1 October 2019, ID 5593).

2.3 Soil water content measurements

The soil water measurements were first described in Metzger et al. (2017). The volumetric soil water content was monitored using a wireless sensor network (SoilNet; Bogena et al., 2010) equipped with SMT100 frequency domain sensors (TRUEBNER GmbH, Neustadt, Germany). Overall, 210 soil water content measurement points were distributed in a stratified random design in the blue subplots shown in Fig. 1: within each blue subplot, two sampling points were placed randomly. Additionally, transects with three additional measurement points (at 0.1, 2.0 and 6.0 m from the position) were added to a subset of 24 randomly selected points. Furthermore, 40 locations were added as transects near tree stems. At each soil moisture measurement location, sensors were installed at two depth: in the topsoil at 7.5 cm and in the subsoil at 27.5 cm. For this analysis, we used the data collected during the throughfall measurement campaigns from 18 June to 28 July 2015 and from 31 May to 14 July 2016.

At each location, we used soil moisture measurements for the hour preceding the observed rain event (θ_pre,ij) to characterize the soil initial moisture in the drained state and the maximum soil water content just after the rain event (θ_post,ij) to characterize the post-event state. We also assessed the soil water content increase due to rainfall by calculating the change in the soil water content (Δθ_ij) for each event j and each location i as follows:

$$\begin{array}{}\text{(2)}& \mathrm{\Delta }{\mathit{\theta }}_{ij}={\mathit{\theta }}_{\text{post},ij}-{\mathit{\theta }}_{\text{pre},ij}.\end{array}$$

Positive values of Δθ_ij indicate a soil water content increase. In the following, we refer generally to Δθ (with indices omitted for simplicity) as “soil water content increase” or “soil moisture response” due to rainfall.

Equivalently to throughfall, we decomposed the soil water content into the event spatial median (${\widehat{\mathit{\theta }}}_{\text{pre},j}$, ${\widehat{\mathit{\theta }}}_{\text{post},j}$) and relative deviations from that median (δθ_pre,ij, δθ_post,ij) with indices for event j and location i omitted for simplicity in the following. As for throughfall, using the relative deviations of the soil water content alongside the medians in the statistical models (see below) provides us with two independent measures for one variable: one relating to the spatial pattern (δθ_pre,ij, δθ_post,ij) and the other to the temporal variation (${\widehat{\mathit{\theta }}}_{\text{pre},j}$, ${\widehat{\mathit{\theta }}}_{\text{post},j}$).

2.4 Canopy and soil property measurements

At the time of soil sensor installation, undisturbed soil samples were collected using metal-ring cylinders with a volume of 100 cm³. The distance between the sensor position and the soil sample collection was approximately 0.5 m. Soil properties were treated as if they were measured directly at the soil sensor location i. In order to determine the field capacity (θ_FC,i), the samples were saturated, left to drain in a sandbox with a hanging water column imposing a pressure of −60 hPa for 72 h and then weighed. The soil cores were subsequently dried for 24 h at 105 ^∘C and weighed again to obtain the dry weight (m_dry,i). The volumetric soil water content at field capacity (θ_FC,i) was derived from the weight difference of the sample at −60 hPa and m_dry,i, assuming a density of water of $D_{\mathrm{w}}=\mathrm{1}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$. Bulk density (D_bd,i) was calculated from the soil dry weight and volume. Soil apparent porosity (φ_i) was calculated from the bulk density and assuming a constant density of the soil mineral component ($D_{\mathrm{m}}=\mathrm{2.66}\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$) as follows:

$$\begin{array}{}\text{(3)}& {\mathit{\phi }}_i=\mathrm{1}-{\displaystyle \frac{D_{\text{bd},i}}{D_{\mathrm{m}}}}.\end{array}$$

Macroporosity (θ_MP,i, also called the air capacity or air-filled porosity) was then determined as follows:

$$\begin{array}{}\text{(4)}& {\mathit{\theta }}_{\text{MP},i}={\mathit{\phi }}_i-{\mathit{\theta }}_{\text{FC},i}.\end{array}$$

To characterize the canopy density, we counted the number of branches (canopy cover) above the throughfall samplers in 2014. However, these data were not available for soil water measurement locations.

2.5 Statistical analysis

All statistical analyses were processed with R (version 3.2.3; R Core Team, 2016). For the geostatistical analysis (detailed below), we used the geoR (Ribeiro and Diggle, 2001), georob (Papritz and Schwierz, 2020) and gstat (Pebesma, 2004; Gräler et al., 2016) packages. Linear mixed-effects models were implemented using the lme4 (Bates et al., 2015) and lmerTest (Kuznetsova et al., 2017) packages. The variance explained by fixed and random factors (conditional R²) and by only fixed effects (marginal R²; Nakagawa and Schielzeth, 2013) for the final model were calculated with the MuMIn package (Bartón, 2020).

2.5.1 Geostatistical estimation of throughfall

Throughfall was estimated at the soil water content measurement locations using kriging. The overall procedure for obtaining the variograms closely follows Zimmermann et al. (2016) with some adaptations taken from Voss et al. (2016). Important steps and decisions of the exploratory data and geostatistical analysis are as follows (see also visualization in Fig. S1 in the Supplement):

1. Exploratory analysis. First, we determined the skewness using the octile skew. None of the throughfall events had an octile skew larger than 0.2 or smaller than −0.2; therefore, we did not transform the data. If a spatial trend existed (p≤0.150), we used the residuals of the spatial regression model for the coordinates x and/or y instead of the real data in the following.

2. Variogram estimation using the method of moments (MoM). This step was performed to derive the best possible preliminary variogram with outliers present, which will be used for spatial outlier detection in Step 3. For this, we first calculated four empirical throughfall variograms using four different variogram estimators. We used both non-robust (Matheron, 1962) and robust (Cressie and Hawkins, 1980; Dowd, 1984; Genton, 1998) estimators as well as the sample.variogram function in the georob package in R. We chose lags centered at 0.125, 0.375 and 0.75, followed by a step size of 1 m up to 50 m. Next, we fitted three models (spherical, exponential and pure nugget) to each empirical variogram using the fit.variogram.model function in the georob package and chose the model with the lowest residual sum of squares for each variogram. This yielded four models, one for each of the variogram estimators stated above. We then assessed the fitted models with leave-one-out cross-validation using kriging. Based on this, we calculated the normalized kriging error (Θ_i,) and selected the empirical variogram with the median Θ nearest to the expected value of 0.455 for Step 3 (Lark, 2000), as done in Zimmermann et al. (2009)

3. Identification and spatial outlier removal. Before final variogram estimation using residual maximum likelihood (REML) in Step 4, spatial outliers were removed based on kriging and cross-validation using the provisional variogram obtained in Step 2. To identify a spatial outlier at location i, we used the standardized error of cross-validation, ε_s,i (Bárdossy and Kundzewicz, 1990; Lark, 2002). To classify an outlier, we used the Z statistics. Sampled points with ${\mathit{\epsilon }}_{\mathrm{s},i}<-\mathrm{2.576}$ ($\mathit{\alpha }/\mathrm{2}=\mathrm{0.005}$) were removed (Zimmermann et al., 2016).

4. Variogram estimation by REML. After outlier removal, we applied REML to fit the theoretical variogram model, including the spatial trend if necessary, using the likfit function in the geoR package. We used the initial estimates from the provisional variogram (Step 2) for the sill, nugget and range parameters. The range relates to the distance over which the observations are still spatially correlated. In the following, we will use the term correlation length to refer to the effective range, e.g., the distance at which the variogram function approaches the sill asymptotically to 95 %. A high effective range indicates a high spatial correlation between two throughfall collectors. We checked the reliability of the final model with the Θ statistic (see above).

5. Kriging. Using the final theoretical variogram from Step 4, we applied ordinary kriging to predict throughfall values at the soil water content measurement locations. Locations where the kriging variance exceeded 95 % of the spatial variance were removed from further analysis.

2.5.2 The coefficient of quartile variation (CQV)

For event-scale assessments, we used quantile-based statistical metrics for descriptive statistics and correlation in order to avoid bias related to the event size or general soil moisture state. Both throughfall and soil moisture patterns can be skewed slightly (Famiglietti et al., 1998; Zimmermann and Zimmermann, 2014), even if the octile skew is less than 0.2, depending on the soil moisture state or event size. Moreover, as mentioned above, throughfall typically includes outliers due to dripping points (Falkengren-Grerup, 1989; Keim et al., 2005; Staelens et al., 2006; Voss et al., 2016). For the coefficient of variation, we used the CQV (Bonett, 2006) as an alternative:

$$\text{CQV}={\displaystyle \frac{Q\mathrm{3}-Q\mathrm{1}}{Q\mathrm{3}+Q\mathrm{1}}},$$

where Q1 and Q3 represent the respective first and third quartiles. Like the classical coefficient of variation, the CQV is dimensionless statistical measure that describes the relative degree of scattering of the sample.

2.5.3 Linear mixed-effects models calculation

We applied linear mixed-effects (LME) models with a repeat-measurement structure to evaluate the influence of potential drivers explaining the soil water content or soil water content increase. We present results on the following dependent variables (see Table 1 for an overview): the spatial pattern of pre-event (δθ_pre) and post-event (δθ_post) soil water content and soil water content increase ($\mathrm{\Delta }\mathit{\theta }=\mathit{\delta }{\mathit{\theta }}_{\text{post}}-\mathit{\delta }{\mathit{\theta }}_{\text{pre}}$).

Table 1 Overview of variables and symbols used in the statistical models.

The independent variables (fixed effects) for δθ_pre were as follows: gross precipitation (P_g), distance to the nearest tree (d_tree), macroporosity (θ_MP), field capacity (θ_FC) and throughfall of the preceding event (P_TF,prec). The independent variables (fixed effects) for Δθ and δθ_post were as follows: gross precipitation (P_g), spatial median of soil pre-event water content (${\widehat{\mathit{\theta }}}_{\text{pre}}$), spatial pattern of soil pre-event water content (δθ_pre), distance to the nearest tree (d_tree), macroporosity (θ_MP), field capacity (θ_FC), spatial median of throughfall (${\widehat{P}}_{\text{TF}}$) and spatial pattern of throughfall (δP_TF). Year, day of year and sensor position were implemented as random effects, thereby accounting for repeated measurements. To avoid model overfitting, it is important that there are no strong correlations between the explanatory variables (Graham, 2003). To detect multicollinearity and to avoid potentially spurious models, we calculated the Spearman rank correlation coefficients (ρ) for all pairs of predictors (Table S1 in the Supplement). Before the analysis, we removed one of a pair of highly correlated predictors: gross precipitation (P_g, strong correlation with ${\widehat{P}}_{\text{TF}}$) and field capacity (θ_FC, strong correlation with θ_MP). All variables (predictor and response) were scaled to center around 0 and have a standard deviation of 1 (z transformation). In this way, the fitted slopes of the model indicate how strongly a change in the predictor within its own range affects the response variable within its own range; hence, the slope estimate gives insight into the effect strength. Scaling has no effect on model selection. To obtain the minimal adequate models for the response variables, we started with the maximum model and removed all nonsignificant terms in a stepwise manner based on the Akaike information criterion (AIC). The main effects included in significant interactions were retained in the model.

3 Results

3.1 Precipitation, throughfall and soil water content pattern

The summer rainfall (May–October) for the last 30 years (1986–2016) shows an average of 352 mm (Mühlhausen-Windeberg). During the two summer periods of this study (2015 and 2016), the annual rainfall was below the long-term mean (276 and 303 mm, respectively). However, the summer of 2015 was the third driest of the last 30 years (Metzger et al., 2017). The final winter months of 2014 were the driest of the 30 years period, and the 2014–2015 hydrological year was the second driest. The 2015–2016 hydrological year and the final winter months of 2015 received an average amount of precipitation.

Descriptive statistics of the throughfall and soil water content (topsoil and subsoil) are given in Table 2. We observed 14 rainfall events in 2015 and 10 events in 2016. The gross precipitation ranged between 1.6 and 35.2 mm, with three small events, six medium events and five large events in 2015 and one medium event and nine large events in 2016. For both years, the soil water content increased with soil depth (Table 2). The soil water content increase (the difference between the post-event and pre-event soil water content, Δθ) was always higher in the topsoil compared with the subsoil. For smaller rainfall events, an increase in the soil water content was mainly limited to the topsoil, whereas an increase in the soil water content was noted at both soil depths following larger rainfall events.

Table 2 Overview of observed rainfall event properties: event date, gross precipitation (P_g), spatial statistics of throughfall (P_TF), soil water content before (θ_pre) and after (θ_post) the rain event, and the soil water content increase (Δθ) in the topsoil and subsoil. The following statistics are presented: spatial median (med), coefficient of quartile variation (CQV), interquartile range (IQR) and effective range (Range).

3.2 Spatial pattern of throughfall and soil water content

The relation between the event size and correlation length of throughfall is shown in Fig. 2. More detailed information about the spatial distributions and the variogram parameters is given in Tables S2 and S3 in the Supplement. Throughfall correlation lengths increased with decreasing event size from on average 6.2 m for large events to 7.5 m for medium and 9.5 m for small events. In comparison, the canopy density correlation length was 7.5 m, i.e., similar to medium events. Throughfall and canopy density had a small nugget and a strong spatial dependence (nugget / sill ratio <25 %) for all events (Table S3). For both years, throughfall decreased significantly with increasing canopy density (Table S4 in the Supplement), although most of the spatial variance in throughfall was related to unknown random effects. While canopy density had no spatial trend (Table S2), throughfall had a spatial trend in less than half of the events (Table S2). These trends changed direction with time, were of varying strength, and occurred in small and large events.

Figure 2 Correlation length, given as effective range, derived from the throughfall variogram calculated for small (P_g<3 mm), medium ($\mathrm{3}\mathrm{mm}<P_{\mathrm{g}}<\mathrm{10}$ mm) and large (P_g>10 mm) events.

The spatial variation in throughfall (interquartile range) increased with event throughfall, but the coefficient of quartile variation (CQV), which normalizes by event size, decreased (Table 2). The high Spearman rank correlation coefficient indicates a strong similarity of the spatial distribution of throughfall between individual events of the same size class (Fig. 3). Thus, throughfall produced persistent wet and dry spots, as also confirmed by time stability plots (Fig. S2 in the Supplement).

Figure 3 Temporal stability of the spatial throughfall patterns. Shown are the pairwise (Spearman) correlation coefficients between throughfall (normalized deviation from the plot median, δP_TF) from different precipitation events, grouped by the event size class: small (n=8), medium (n=11) and large (n=21).

The soil water content spatial variation coefficient (CQV) values decreased with increasing soil water content (expressed as the spatial median) and, consequently, with increasing soil depth (Table 2 and Fig. S3 in the Supplement). In the topsoil, the highest variation occurs in the post-event soil water content (Fig. S3) and is substantially higher than the pre-event soil water content, indicating that the event response enhanced the soil water content variation, especially under moderately dry (summer) conditions in topsoil. However, the highest CQV values by far were observed for the increase in the soil water content after rain (Δθ).

The pairwise correlation coefficients indicating the temporal stability of the spatial patterns were high for the pre-event (drained) soil water content (θ_pre) in both topsoil (Fig. 4a) and subsoil (Fig. 4b), with ρ≈0.78. For the post-event soil water content (θ_post), they were significantly lower in the topsoil (ρ=0.70; Fig. 4c) than in the subsoil (ρ=0.77; Fig. 4d) (Mann–Whitney U test: $Z=-\mathrm{3.15}$, p=0.002). In the topsoil, they decreased with increasing event size, revealing that patterns were less similar after large precipitation events (Fig. 4a, c). In contrast, the spatial distribution of the soil water content increase after rain events ($\mathrm{\Delta }\mathit{\theta }={\mathit{\theta }}_{\text{post}}-{\mathit{\theta }}_{\text{pre}}$) changed much more between events (Fig. 4e, f), leading to an overall lower correlation between the patterns. However, the similarity of the spatial distribution of Δθ increased with event size, especially in topsoil (Fig. 4e), confirming reoccurring wetting patterns, especially following larger events.

Figure 4 Temporal autocorrelation of the spatial patterns of the pre- and post-event soil water content and increase in the soil water content after rainfall, calculated as pairwise (Spearman, ρ) correlation coefficients between all of the different precipitation events within the following event size classes: small (n=3), medium (n=7) and large (n=13). The figure panels show the following variables: (a, b) pre-event soil water content (θ_pre), (c, d) post-event soil water content (θ_post), and (e, f) increase in the soil water content ($\mathrm{\Delta }{\mathit{\theta }}_i={\mathit{\theta }}_{\text{post}}-{\mathit{\theta }}_{\text{pre}}$) in the (a, c, e) topsoil (−7.5 cm) and (b, d, f) subsoil (−27 cm). The differences between the event size classes were examined using a Duncan multiple-range test. The letters above the bars indicate a significant difference (p≤0.05) between the groups.

Table 3 Factors affecting the pre-event soil water content patterns (δθ_pre) in the topsoil and subsoil. Results are shown for the best linear mixed-effects model; significant effects are highlighted in bold.

3.3 Factors influencing the soil water spatial distribution

3.3.1 Soil water content

In order to identify the basic drivers of the soil water content patterns in the drained state (δθ_pre), we used mixed-effects model selection. The resulting best models for the top- and subsoil are given in Table 3. The variance explained by fixed effects (marginal R²) was low, whereas the variance explained by fixed and random effects together (conditional R²) was high. The model for the subsoil showed an even higher marginal R² compared with the topsoil as well as a somewhat higher influence of fixed effects. The most important effect identified for the topsoil and subsoil was macroporosity, with a lower soil water content (δθ_pre) related to locations of higher macroporosity (Table 3). In the topsoil, the throughfall of the preceding precipitation event also slightly affected the soil moisture pattern. The results for the soil water content itself in the drained state (θ_pre) are similar to those of δθ_pre, except that fixed effects explain even less variation (not shown).

The results of the best linear mixed-effects model relating the soil water content after a precipitation event to potential drivers are given for all events (Fig. 5a, c) and large events only (events with P_g>6 mm; Fig. 5b, d). The median initial soil water content (soil water content before the rain event, θ_pre) and its spatial pattern (δθ_pre) were the major drivers of either the absolute values of the spatially distributed soil water content after the rain event (θ_post; Fig. 5a, b) or the spatial pattern of the soil water content after the rain event (δθ_post). Other fixed (${\widehat{P}}_{\text{TF}}$, δP_TF, ${\widehat{\mathit{\theta }}}_{\text{pre}}$, θ_MP and d_tree) and random effects contributed only a little, especially when small and medium events were excluded (Fig. 5b).

Figure 5 Factors influencing (a, b) the local post-event soil water content (the peak soil water content after rainfall, θ_post) and (c, d) the local soil water content increase due to rainfall (the difference between the soil water content after and before each event, Δθ). Slope estimates are shown for the best linear mixed-effects model including (a, c) all events and (b, d) large events only (P_g>6 mm). Significant effects are shown using thick lines, and the asterisks indicate the level of significance: ${}^{\ast \ast \ast }$ p<0.001, ${}^{\ast \ast }$ $\mathrm{0.001}\le p<\mathrm{0.01}$ and ^∗ $\mathrm{0.01}\le p<\mathrm{0.05}$. All variables (predictors and responses) are z-scaled such that the slope indicates the effect strength. Pseudo-R² values and a summary of all models, including those of small and medium events, are given in the Supplement (Tables S5 and S6).

3.3.2 Soil water response (Δθ)

The slope estimates of models explaining the soil water content increase (Δθ), i.e., how much water was locally added to the soil after rain, are visualized in Fig. 5c and d for all and large events, respectively. In general, a detectable (>1 %) change in Δθ was limited to large rainfall events (Table 2). The spatial patterns responded to several drivers (fixed effects) in the final model. There, the variance explained by fixed effects (marginal R²) was generally higher for subsoil compared with topsoil; moreover, it typically increased with event size and was highest for the models including all event sizes (Table S6). Therefore, in the following, we focus on the effects emerging from models including all events. We also visualize models for events falling in the large-event class, as it covers more than 80 % of the cumulated net precipitation received. The results for the individual event size classes are given in the Supplement (Tables S5 and S6).

Overall, local soil water content increase (Δθ_i) depended not only on the event median throughfall (${\widehat{P}}_{\text{TF}}$) but also on whether more or less throughfall (δP_TF) was received locally and on local soil moisture conditions (δθ_pre). In the subsoil, macroporosity additionally directly affected the soil moisture response. All main effects are also included in interactions, meaning that a third variable influenced the relationship between an independent and dependent variable. For example, locally elevated throughfall enhanced the soil water increase, although more so with increasing event size (see the interaction of ${\widehat{P}}_{\text{TF}}\times \mathit{\delta }{P}_{\text{TF}}$ in Fig. 5c and d and Fig. 6a and b).

Figure 6 Marginal plot of the interaction between the throughfall pattern (deviation of throughfall from the spatial median) and the median event throughfall ($\mathit{\delta }{P}_{\text{TF}}\times {\widehat{P}}_{\text{TF}}$) affecting the soil water content increase (Δθ) for all events in topsoil (a) and subsoil (b). Shown is the influence of the throughfall pattern (δP_TF) on the local soil water content response (Δθ_i), grouped by the event size (the legend indicates the group central median event throughfall, ${\widehat{P}}_{\text{TF}}$). Note that all values are z-scaled such that the slope indicates the effect strength. For example, the yellow line highlights small events, where the local soil water moisture response depends little on the spatial distribution of throughfall input. For large events, marked in blue, the soil moisture response is stronger in locations of above-average throughfall.

Spatial patterns of the pre-event (or initial or drained) soil water content (δθ_pre) notably affected top- and subsoil differently, making it the only factor yielding inverse effects on the soil water content increase at different soil depths. In topsoil, drier locations stored less water per event than moister spots (positive slope estimates), whereas the opposite was the case in subsoil (negative slope estimates). Notably, the slope of the interaction changes with the overall soil water conditions, consistently at both depths (see the interaction of ${\widehat{\mathit{\theta }}}_{\text{pre}}\times \mathit{\delta }{\mathit{\theta }}_{\text{pre}}$ in Fig. 5c and d and Fig. 7): under very dry soil conditions (summer, topsoil), locally moister soil admitted more water into the soil (positive slope in Fig. 7a), but less water was admitted under overall wet conditions (negative slope in spring topsoil and permanently in subsoil). This influence of the pre-event soil moisture patterns on moistening increased with event size (significant interaction ${\widehat{P}}_{\text{TF}}\times \mathit{\delta }{\mathit{\theta }}_{\text{pre}}$, not visualized).

Figure 7 Marginal plot of the interaction between the initial soil water content pattern (deviation of the pre-event soil water content from the spatial median) and the median soil water content ($\mathit{\delta }{\mathit{\theta }}_{\text{pre}}\times {\widehat{\mathit{\theta }}}_{\text{pre}}$) affecting the soil water content increase (Δθ) for all events in topsoil (a) and subsoil (b). Groups indicate the overall soil moisture conditions (the legend indicates the group central spatial median of the pre-event soil water content, ${\widehat{\mathit{\theta }}}_{\text{pre}}$), and graphs show the relation between the initial soil water pattern and the local soil moisture response. All values are z-scaled such that the slope indicates not only the direction but also the strength of the effect. For example, the blue line shows how, under generally moist conditions (e.g., early spring), the soil moisture response to rain is dampened in moister locations (high values of δθ_pre), more prominently so in the subsoil. Additionally, when the topsoil is dry (e.g., summer, red and yellow lines), dry locations also store less water. In combination, this leads to a change in the slope direction in the topsoil over the year. The subsoil water content is always higher than that in topsoil; negative slopes (dampening in moist locations) are shown throughout.

Additionally, in the topsoil, the distance to the next tree affected the soil water response. Locations near trees reacted more strongly to event precipitation than those further away (${\widehat{P}}_{\text{TF}}\times d_{\text{tree}}$ interactions), although only under generally moister soil conditions (interaction of ${\widehat{\mathit{\theta }}}_{\text{pre}}\times d_{\text{tree}}$ in Fig. 5c and d). In the subsoil, higher macroporosity (θ_MP) dampened the soil water response (negative slope; Fig. 5c and d), more so when or where throughfall was high (interactions of ${\widehat{P}}_{\text{TF}}\times {\mathit{\theta }}_{\text{MP}}$ and θ_MP×δP_TF in Fig. 5c and d) and in drier locations (interaction of θ_MP×δθ_pre in Fig. 5c and d).

4 Discussion

4.1 Strengths and weaknesses of the approach

In this analysis, we used extensive spatial data of canopy cover, throughfall and soil water content in order to assess the role of canopy processes in the belowground soil water response to precipitation. For this, we measured precipitation and the soil water content at different locations in order to avoid the disturbance of soil water dynamics due to precipitation measurement and to provide independent random measurement designs. To be able to relate observations at different locations, we used geostatistical methods to predict throughfall values at locations where the soil water content was measured. Throughfall spatial prediction was based on an extensive dataset of a substantial number (n=350) of samplers of comparatively small size (support is A=0.011 m²) in a stratified random design, spanning an extent more than 10 times the correlation lengths of most events. Throughfall showed strong spatial autocorrelation which was reflected by 90 % of the nugget / still ratios being lower than 10 % (Table S3). With the tight sampling at the shortest lag distance (at least 50 locations covered two samplers located directly next to each other), the nugget effect, or unresolved variance, can be attributed to the dispersion variance across the sampler (support) and to the spatial field shifting slightly within the rain event. The choice of the support scale affects the variance estimate, as demonstrated for throughfall by Zuecco et al. (2014). Here, the scale of the sampler and the soil moisture sensors is roughly the same; thus, the variance is appropriately captured for the intended purpose. The spatial correlation depended on the event size in that the correlation length decreased with increasing event rainfall (see the examples in Fig. S4 in the Supplement). In larger events, this decreased the range within which throughfall could be predicted and increased the number of locations with high kriging variance that were removed from the analysis. As a result, this decreased the sample size for large events compared with small and medium sized events. Regardless, for all sampled events, we could still rely on datasets of 59 points on average.

In order to estimate a reliable variogram, we had to remove spatial outliers contaminating the sample, which demonstrates that the throughfall spatial field is not entirely smooth. Outlier locations comprised on average 2.2 % of the spatial sample, were rarely recurring, and were not only related to throughfall hot spots (although they were mainly related to throughfall hot spots for large events). Throughfall hot spots have been related to interrupted flow lines leading to dripping points (Crockford and Richardson, 2000; Herwitz, 1987; Nanko et al., 2022; Staelens et al., 2006; Zimmermann et al., 2009) that can also move around slightly (Keim et al., 2005). In removing these outliers, our analysis neglects the role of dripping points in the soil moisture response. Additionally, kriging tends to smooth the estimates compared with the actual data (Oliver and Webster, 2014). However, the predicted throughfall showed approximately the same median and spatial variance here as the measured data (without outliers), indicating that this was not a concern and that the real variation was captured after the prediction procedure. Unfortunately, there is no perfect way to relate measurements obtained at different locations to each other. However, a large sample size of throughfall, a stratified design, outlier detection with robust variogram estimators and the use of residual maximum likelihood (REML) are established tools for estimating the variogram reliably (Lark, 2000; Voss et al., 2016). The latter provides an essential basis for kriging (Oliver and Webster, 2014). Cross-validation of the kriging estimates against observations (Lark, 2000; Oliver and Webster, 2014) provides further confidence in the variogram and kriging procedure for interpolating the aboveground data to the belowground locations.

In our analysis, we quantified only throughfall input and omitted the role of stemflow, which may play a role in locations near stems. Extrapolating stemflow input to soil moisture locations entails more prediction steps compared with throughfall. The spatial variation in stemflow depends not only on the species, tree and canopy size, neighborhood and individual morphology of the trees (Bellot and Escarre, 1998; Fan et al., 2015a; Levia et al., 2014; Levia and Germer, 2015; Van Stan et al., 2016; Metzger et al., 2019; Magliano et al., 2019) but also on the precipitation intensity and soil conditions determining the infiltration area (Herwitz, 1986; Carlyle-Moses et al., 2018; Metzger et al., 2021). Such a prediction would not only introduce a great deal of uncertainty but would also deviate from the main purpose of this study, which is to evaluate the role of throughfall heterogeneity. Therefore, in the model analysis, microsites near stems were accounted for by including the distance to the stem as a fixed effect in the model. This takes into account, to some extent, the potential influence of stemflow in the interpretation.

4.2 General (temporal and spatial) patterns of throughfall

In agreement with previous studies, the spatial variation coefficients of throughfall decrease with event size (Aussenac, 1970; Loustau et al., 1992; Llorens et al., 1997; Su et al., 2019; Metzger et al., 2017; Carlyle-Moses, 2004; Staelens et al., 2008; Van Stan et al., 2020). Several other studies have suggested that the spatial variation in throughfall depends on the amount of precipitation as well as on canopy characteristics (Loustau et al., 1992; Carlyle-Moses, 2004; Keim et al., 2005; Park and Cameron, 2008; Hsueh et al., 2016; Zimmermann et al., 2009; Herwitz and Slye, 1995). Similarly, at our site, canopy cover was a significant driver of the throughfall spatial distribution for all event size classes, although a small driver compared with the random effects. The spatial trends showed no clear pattern with event size and may have been related to a combination of slope aspect and wind conditions. The correlation length (effective range) of throughfall decreased with increasing event size and roughly corresponded to that of canopy cover for medium events. The change in the spatial pattern with event size underlines that not only canopy storage per se but also other processes like turbulence, wind shadows, the arrangement of canopy gaps or the formation of canopy dripping points can add persistent spatial organization to below-canopy precipitation (Carlyle-Moses, 2004; Keim et al., 2005; Park and Cameron, 2008; Staelens et al., 2008; Zimmermann et al., 2008; Wullaert et al., 2009; Li et al., 2016; Van Stan et al., 2020). Additionally, canopy features also affect within-canopy redistribution (André et al., 2011; Herwitz, 1987; Levia and Frost, 2006; Nanko et al., 2022), which could lead to reoccurring patterns not reflected by canopy density.

Overall, despite the changes in throughfall correlation lengths for different event size classes, throughfall patterns were temporally stable, indicating the existence of permanent hot and cold spots of throughfall, and these hot and cold spots were consistent across small, medium and large events. This is in line with several previous studies that have reported the temporal stability of throughfall patterns (Keim et al., 2005; Staelens et al., 2006; Wullaert et al., 2009; Zimmermann et al., 2009; Fathizadeh et al., 2014; Fan et al., 2015a; Metzger et al., 2017; Molina et al., 2019; Zhu et al., 2021; Rodrigues et al., 2022), even over several years (Wullaert et al., 2009; Rodrigues et al., 2022), although phenology and canopy development have also been observed to deteriorate spatial stability (Zimmermann et al., 2008; Fathizadeh et al., 2014). Furthermore, although spatial variation coefficients are smaller for large events compared with small events, absolute values vary much more in large events such that they have an arguably higher potential to induce spatial patterns in the soil water content or dynamics.

4.3 General soil water content patterns and potential drivers

The mean soil water content was generally lower in the topsoil compared with the subsoil. At our site, the shallow soil is underlain by undulating weathered, calcareous bedrock (Kohlhepp et al., 2017) of low hydraulic conductivity, and it may be locally broken through by tree roots. While the topsoil is well drained (i.e., saturated to field capacity in winter and much lower in summer), the deeper and finer-textured soil layer (Metzger et al., 2021) is influenced by the much less conductive regolith and a generally moister soil water content which very occasionally exceeds field capacity in winter (Metzger et al., 2017).

Much in agreement with previous studies in humid regions (Brocca et al., 2007; Korres et al., 2015; Rosenbaum et al., 2012; Metzger et al., 2017), the spatial variation in the soil water content increased in both the topsoil and subsoil under drier summer soil conditions. In an earlier study at the same site, a strong but short-lived increase in the spatial variation in the topsoil water content in summer was related to precipitation events (Metzger et al., 2017). Regardless, we found that the main controlling factor of the post-event soil water content was the spatial pattern of the pre-event soil water content, while the average throughfall and spatial pattern of throughfall, tree distance, and macroporosity were additional, although much less important, drivers. In other words, although the soil water content variation increases strongly after events, this variation can only be traced back to input patterns in a very limited fashion. This may, in part, be due to the small inputs of water compared with the overall soil water storage, leading to a strong memory effect of the pre-event soil water conditions on the post-event patterns. Furthermore, preferential flow already taking place during the event itself can blur the throughfall pattern within the soil storage (see below).

The soil water content spatial patterns in a drained state, in turn, were strongly driven by random effects. Those are factors that were not described by the measurements but that are temporally stable. These so-called local soil conditions are potentially related to soil hydraulic properties, root water uptake and microtopography (Famiglietti et al., 1998; Vereecken et al., 2007; Fan et al., 2015b). The mixed-effects models confirm, although with a very weak influence, that locations of higher macroporosity were drier at both depths, confirming the role of water retention in soil water patterns (Metzger et al., 2017) at this site. Moreover, the throughfall patterns of the previous event slightly affected the topsoil pre-event soil water content. Thus, an imprint of the throughfall pattern was carried over to the next pre-event soil conditions, but this is barely detectable and is negligible compared with the other sources of variation in the soil water content in a drained state.

4.4 Drivers of the soil water response (Δθ) to rainfall

In contrast to the absolute soil water content values discussed above, the local soil water response (i.e., increase in the soil water content following rainfall events) was clearly affected by the spatial throughfall pattern in both topsoil and subsoil, although this effect was modified by soil properties. As we tested the effect of the spatial pattern (δP_TF) separately from the absolute values of event throughfall (${\widehat{P}}_{\text{TF}}$), we are able to demonstrate the influence of spatial throughfall fields specifically. Among all of the drivers tested, the influence of spatial throughfall variation was not necessarily the strongest, but it was a consistently reoccurring factor. It appeared at both observed soil depths with a similar influence and was more pronounced for larger events. In other words, the spatial variation in the throughfall was a consistent driver of soil wetting.

Measurements ascertaining that the soil water content response is related to canopy drainage patterns are comparatively rare. Metzger et al. (2017) previously reported, for the same site but using a smaller dataset, that increases in the soil water content were correlated with event spatial throughfall patterns during larger rainfall events. Furthermore, Molina et al. (2019), using high-temporal-resolution soil moisture measurements, found a relationship between the average throughfall pattern and the soil water content response in the topsoil of a Mediterranean oak-dominated forest plot but not in a pine plot. Notably, Klos et al. (2014) showed that locations of high and low soil water content below the main rooting zone in a tropical rainforest corresponded to the end-members of high and low throughfall, whereas the soil water content above and below this depth was more homogenous. They concluded, from additional modeling, that preferential flow may have contributed to bypassing the main rooting zone. On the other hand, several studies, such as Raat et al. (2002), Shachnovich et al. (2008), and, more recently, Zhu et al. (2021) using lower-temporal-resolution soil water content measurements (incidentally all in coniferous forests) found no relationships between the spatial patterns of the soil water content and throughfall. All authors reported that the throughfall patterns were pronounced and stable in time and suggest that the forests floor impeded the transfer to soil water patterns. An additional explanation could be that the effect of spatial net precipitation patterns on the soil water content was so short-lived due to preferential flow (Metzger et al., 2017) that it was not observed by infrequent hand measurements. Overall, the stronger soil water response at sites with above-average throughfall indicates that throughfall hot spots and cold spots (Levia and Frost, 2006; Van Stan et al., 2020; Zimmermann et al., 2009) have an impact on the soil water dynamics, although they go almost unnoticed in the soil water content pattern (see above).

In addition to the throughfall pattern, the soil water response after large rainfall events also depended on the pattern of the pre-event soil water content at both depths. Notably, the slope of the relationship changes direction, making it the only factor that shows inverse effects in the top- and subsoil. This can be attributed to its interdependence on the soil water content and the difference in moisture between the two measurement depths. Especially under dry (summer) conditions, wetter topsoil locations stored more of the incoming precipitation water, while drier sites remained dry. This is a strong indication of preferential flow in dry soils, where, for example, hydrophobic conditions, cracks and low hydraulic conductivity of the matrix can enhance preferential flow (Hillel, 1998; Nimmo, 2021; Beven and Germann, 2013). On the other hand, the dampened water response in the wetter subsoil could be due to enhanced hydraulic conductivity and less free pore space (Vereecken et al., 2007; Hagen et al., 2020). It is noteworthy that effects under dry conditions are much stronger than under wet conditions, suggesting a stronger trigger for preferential flow there. Moreover, only at intermediate soil moisture conditions, the spatial distribution of pre-event soil water content had no influence on the spatial drainage behavior.

The soil water response also depended on the distance to the nearest tree in the topsoil and on soil properties (macroporosity) in the subsoil. The enhanced moistening of soils near stems is likely related to stemflow production (Metzger et al., 2019), which was not accounted for as input. Stemflow production generally increases with event size (Levia and Germer, 2015; Metzger et al., 2019), explaining the interaction in the model. The additional modification due to the soil water conditions can be explained by the systematically lower soil water content near tree trunks at the same site (Metzger et al., 2017, 2021), owing to lower soil water retention and likely enhanced drainage there.

Taken together, our data strongly suggest that, in addition to the spatial distribution of throughfall, the spatial pattern in drainage behavior affects the local soil water response to rainfall. Thus, both dry and wet locations can, water supply permitting, act as percolation hot spots, depending on the overall soil conditions. Bypass flow in forests has been repeatedly observed (e.g., Schume et al., 2003; Schwärzel et al., 2009; Bachmair et al., 2012; Blume et al., 2009; Demand et al., 2019), especially under dry summer conditions (Schume et al., 2003; Bachmair et al., 2012; Demand et al., 2019). The spatial variation in the infiltration water supply and intensity, such as is the case for below-canopy precipitation (Keim and Link, 2018), has been suggested as a potential cause of the initiation of finger flow (Nimmo, 2021), which is promoted by dry soil conditions. Moreover, hydrophobicity has been suggested to contribute to maintaining recurring finger flow paths (Blume et al., 2009). Furthermore, macropore flow along biopores (Lange et al., 2009; Nespoulous et al., 2019) and soil cracks (Schume et al., 2003) can be enhanced under dry forest soil conditions due to soil shrinking (Baram et al., 2012). While both finger flow and macropore flow may have contributed to the observed patterns in soil water response, macropore flow (more than finger flow) could explain enhanced matter export (Lehmann et al., 2021) as well as the fast response following the strong storms observed in the shallow aquifers of the AquaDiva Critical Zone Observatory (Lehmann and Totsche, 2020).

Overall, our results confirm that the effect of throughfall on the soil water content is weak, but this effect is stronger on the soil water response. This contrasts with previous modeling work (Coenders-Gerrits et al., 2013) that did not account for preferential flow. As the effect of the throughfall pattern on the soil water response also depends on local conditions related to hydraulic properties, its fate is much more likely to be found in the drainage fluxes, rather than the soil water storage. The further destiny of the net precipitation pattern arguably depends on the deeper subsurface hydrogeological setting. We deduce, however, that net precipitation hot spots have a strong chance of producing patterns of preferential flow below the main rooting zone, which is in line with previous work (Klos et al., 2014) and backs earlier hypotheses (Bouten et al., 1992; Schume et al., 2003).

5 Conclusion

In this study, we collected an extensive dataset to investigate the effect of throughfall spatial heterogeneity on the soil water response and checked which other factors (e.g., pre-event soil water content, macroporosity, and tree distance) modified the result. We first confirmed that throughfall patterns were stable in time and found that they related to the vegetation canopy density, although additional and partly unknown factors strongly affected throughfall distribution. We found that the post-event soil water content did have a very weak intrinsic relation to throughfall, although the variation in the soil water content clearly increased in the aftermath of rain events. The post-event soil water content pattern was overwhelmingly determined by the strong memory effect of the soil water storage and only slightly affected by soil properties, like macroporosity. In contrast, the soil water response showed a clear relation with the throughfall input pattern. In other words, our setup allowed us to confirm experimentally that throughfall patterns do imprint on soil water content dynamics, at least shortly after rain events. However, we also identified locations where the soil water response was dampened, likely due to enhanced fast drainage. These locations could be either very dry locations, likely promoting preferential flow, especially in the topsoil, or wet locations, promoting faster release of the incoming water. Our results demonstrate that throughfall spatial patterns leave a stronger imprint on soil water dynamics than on the soil water content directly, and they explain why the aboveground influence on soil hydrology has been so difficult to lay open in the past. Our findings are in line with previous research and contribute a more general process understanding of the belowground consequences of precipitation redistribution by forests. Most importantly, our results strongly suggest that throughfall patterns induce fast soil water flow with repeating spatial patterns. These patterns would, therefore, already be triggered within the canopy.