Spatial Patterns and Temporal Stability of Throughfall in a Mature Douglas-fir Forest

Forest plays a key role in spatial distribution of rainfall and nutrients at fine spatial scales. Areas of localized rainfall and nutrient input at the soil surface may have a large effect on several hydrological and biogeochemical processes. In this paper, a Douglas-fir stand was revisited to evaluate the changes in the throughfall spatial distribution and its temporal stability due to forest growth and thinning. We used 32 funnel-type collectors distributed in a random stratified array within a 0.2 ha plot to measure throughfall amounts from February to November 2015. The throughfall variability was much lower as compared to the values reported ~25 years ago in the same site. We further assessed the spatial patterns of throughfall in spring and summer. We detected a spatial correlation length of 12 m and 8 m for spring and summer, respectively, which are higher than the values reported for other mature Douglas-fir forests in similar climatic conditions. Temporal stability plots confirmed that detected spatial patterns were stable in time.


Introduction
Forest ecosystems play a significant role in the redistribution of water and nutrients that reaches the terrestrial surface by means of throughfall (Tf ) and stemflow (Sf ) [1][2][3][4].Throughfall and stemflow patterns below forest canopies have been shown to be highly variable in time and space [2, 5,6].The spatial variability of Tf is a critical issue at different temporal and spatial scales, because it has a direct impact on several hydrological and biogeochemical processes such as soil water dynamics [7], water chemistry [1,8,9], nutrient cycling [3,10], root water uptake [7], and root growth [11].Moreover, the spatial variability of throughfall is known to affect the accuracy of estimates on stand-scale interception losses [12,13].Thus, evaluations of spatial variability of Tf are of practical importance for effective measurements in the estimation of stand-scale Tf.
Spatial variability of throughfall has been traditionally expressed using coefficient of variation (CV), e.g., [1,[13][14][15][16].Several studies demonstrated that CV decreases asymptotically with the increase in size of the rainfall events, e.g., [13,15,[17][18][19].Similarly, CV decreases with the increase of the total rainfall depth in studies that aggregate Tf from several events at larger temporal resolution (i.e., weekly, monthly) [6].Less common in the literature is the description of spatial patterns of throughfall using geostatistical methods, which can be used to produce Tf distribution maps as long as structured variograms are found.However, in cases in which the spatial variability has been studied in terms of correlation length, the results have been highly variable.For instance, Keim et al. [20] quantified patterns of Tf using variograms and identified the spatial correlation lengths of throughfall of ~5 to 10 m in coniferous stands.Also Staelens, et al. [14] found spatial correlation lengths of 3 to 4 m in the leafed period of a beech stand.In contrast, some researchers did not find any spatial autocorrelation in Tf data [17,21].Such large differences in spatial patterns among forest ecosystems can

Rainfall
Gross rainfall (PG, mm) was measured in a nearby, well-exposed clearing (ca.250 m from the center of the plot) using two tipping bucket rain gauge (Rain Collector II, Davis Instruments, Hayward, CA, USA) with a resolution of 0.2 mm per tip.The orifice of the rain gauge was positioned at 1.5 m above the ground to avoid ground-splash effects.The automatically recorded data were stored by a HOBO event logger at 1-min interval (Onset Computer Corporation, Bourne, MA, USA).Gross rainfall was also collected at the top of the 47 m scaffolding tower operated by University of Twente (ITC-UT) (at ca.200 m distance from the clearing) using two tipping bucket rain gauges (Onset HOBO-RG3, resolution 0.2 mm).The data at the top of the tower were only used to fill a few gaps (from 23 July 2015 to 12 August 2015; 24 May 2016 to 9 June 2016) in the data at the clearing using a linear regression equation that linked 10-min rainfall totals at the two locations (R 2 = 0.93, n = 1000).
Because the instrumentation for PG was installed on 01 June 2015, and throughfall was started earlier than that (see below), cumulative rainfall data for the first six periods of measurement from 17 February 2015 to 12 June 2015 were obtained by inverse distance weighted (IDW) interpolation using hourly cumulative rainfall from the three nearest KNMI weather stations in de Bilt, Lelystad, and Deelen located at 39, 25, and 25 km, respectively, from the study area.

Throughfall
A stratified random sampling approach was used for the throughfall measurement (Tf, mm) to ensure an even spread of sampling locations.Plot size was decided based on (i) previous studies carried out in the same plot cf.[1,7] to allow a direct comparison and (ii) previous studies with similar

Rainfall
Gross rainfall (P G , mm) was measured in a nearby, well-exposed clearing (ca.250 m from the center of the plot) using two tipping bucket rain gauge (Rain Collector II, Davis Instruments, Hayward, CA, USA) with a resolution of 0.2 mm per tip.The orifice of the rain gauge was positioned at 1.5 m above the ground to avoid ground-splash effects.The automatically recorded data were stored by a HOBO event logger at 1-min interval (Onset Computer Corporation, Bourne, MA, USA).Gross rainfall was also collected at the top of the 47 m scaffolding tower operated by University of Twente (ITC-UT) (at ca.200 m distance from the clearing) using two tipping bucket rain gauges (Onset HOBO-RG3, resolution 0.2 mm).The data at the top of the tower were only used to fill a few gaps (from 23 July 2015 to 12 August 2015; 24 May 2016 to 9 June 2016) in the data at the clearing using a linear regression equation that linked 10-min rainfall totals at the two locations (R 2 = 0.93, n = 1000).
Because the instrumentation for P G was installed on 01 June 2015, and throughfall was started earlier than that (see below), cumulative rainfall data for the first six periods of measurement from 17 February 2015 to 12 June 2015 were obtained by inverse distance weighted (IDW) interpolation using hourly cumulative rainfall from the three nearest KNMI weather stations in de Bilt, Lelystad, and Deelen located at 39, 25, and 25 km, respectively, from the study area.

Throughfall
A stratified random sampling approach was used for the throughfall measurement (Tf, mm) to ensure an even spread of sampling locations.Plot size was decided based on (i) previous studies carried out in the same plot cf.[1,7] to allow a direct comparison and (ii) previous studies with similar approaches in similar forest type elsewhere cf.[19,20].Based on the above-mentioned criteria and logistic limitations, we defined a plot size of 32 × 64 m that was divided into 32 square sub-plots of 8 m × 8 m each.Considering the study of Keim, et al. [20], in which the plot size was selected to be at least three times of the crown size, the width and length of our plot were sufficient (3 to 6 times larger than the crown size, respectively) to perform spatial analysis.Collectors (32 in total) were placed at some distance from each marked point by generating random values for an azimuth angle and distance from the grid point.The azimuth angles ranged from 0 to 360 degrees and the distances from 0 to 4 m (Figure 2).In case the randomly selected position coincided with the position of a stem, the azimuth angle was maintained while the distance was adjusted until the collector was located next to the tree base and the adjusted distance was recorded.The funnel-type collectors consisted of a 2 L collector and a funnel (165 cm 2 orifice area).The orifices of the gauges were positioned 50 cm above the forest floor to avoid splash-in from the ground.The funnel-type collectors were read ~bi-weekly.Measured Tf volumes were converted to equivalent depth (in mm) by dividing the volume of water in each gauge by the orifice area.From 17 February 2015 to 2 November 2015, Tf was measured 15 times (herein after referred to as Periods (Prd); Table 1).For the Prd-1 to Prd-10, we applied the roving method by randomly generating new positions and placing the funnel-type collectors after each Tf measurement.For the Prd-11 to Prd-15, we used the non-roving technique (i.e., the gauges were not relocated after measurement).The purpose for applying the non-roving technique during the last stage of the study was to examine the temporal stability of the Tf in the plot.
The 15 periods of measurements were grouped by season as winter (Prd-1 and Prd-2), spring (Prd-3 to Prd-6), summer (Prd-7 to Prd-11), and fall (Prd-12 to Prd-15).The CV of each period was calculated in order to compare with previous estimations in the same site by Raat, et al. [1].Periods of Tf measurements were also characterized by the number of rainfall events that occur within the period.An event was defined as the period of time with at least 0.4 mm of rain, separated by a dry period of at least 3 h cf.[27].
Water 2018, 10, x FOR PEER REVIEW 4 of 15 approaches in similar forest type elsewhere cf.[19,20].Based on the above-mentioned criteria and logistic limitations, we defined a plot size of 32 × 64 m that was divided into 32 square sub-plots of 8 m × 8 m each.Considering the study of Keim, et al. [20], in which the plot size was selected to be at least three times of the crown size, the width and length of our plot were sufficient (3 to 6 times larger than the crown size, respectively) to perform spatial analysis.Collectors (32 in total) were placed at some distance from each marked point by generating random values for an azimuth angle and distance from the grid point.The azimuth angles ranged from 0 to 360 degrees and the distances from 0 to 4 m (Figure 2).In case the randomly selected position coincided with the position of a stem, the azimuth angle was maintained while the distance was adjusted until the collector was located next to the tree base and the adjusted distance was recorded.The funnel-type collectors consisted of a 2 L collector and a funnel (165 cm 2 orifice area).The orifices of the gauges were positioned 50 cm above the forest floor to avoid splash-in from the ground.The funnel-type collectors were read ~bi-weekly.Measured Tf volumes were converted to equivalent depth (in mm) by dividing the volume of water in each gauge by the orifice area.From 17 February 2015 to 2 November 2015, Tf was measured 15 times (herein after referred to as Periods (Prd); Table 1).For the Prd-1 to Prd-10, we applied the roving method by randomly generating new positions and placing the funnel-type collectors after each Tf measurement.For the Prd-11 to Prd-15, we used the non-roving technique (i.e., the gauges were not relocated after measurement).The purpose for applying the non-roving technique during the last stage of the study was to examine the temporal stability of the Tf in the plot.The 15 periods of measurements were grouped by season as winter (Prd-1 and Prd-2), spring (Prd-3 to Prd-6), summer (Prd-7 to Prd-11), and fall (Prd-12 to Prd-15).The CV of each period was calculated in order to compare with previous estimations in the same site by Raat, et al. [1].Periods of Tf measurements were also characterized by the number of rainfall events that occur within the period.An event was defined as the period of time with at least 0.4 mm of rain, separated by a dry period of at least 3 h cf.[27].Note: a Prd-11 to Prd-15 had stationary arrangement.b Rainfall amounts for Prd-1 to Prd-6 were interpolated (IWD) from the three nearest KNMI weather stations (de Bilt, Deelen, and Lelystaad).c An event was defined as the period of time with at least 0.4 mm of rain, separated by a dry period of at least 3 h.d Mean of rainfall intensity (event based).For Prd-1 to Prd-6 was derived from data of KNMI weather station (Deelen).e Mean of wind speed (event based).For Prd-1 to Prd-6 was derived from data of KNMI weather station (Deelen).f Longest duration of an event is the maximum duration in hours of a rainfall event during the period of Tf collection.

Spatial Variability of Throughfall
The regionalized variable theory [28] was applied to model spatial correlation structures by means of variograms.The variogram model was estimated using the method-of-moments (MoM) that consists of two steps: calculating the experimental variogram and fitting a variogram model.The experimental variograms are usually calculated according to Matheron [29] equation which was applied as follows: in which γ(h) is the semi variance, T f i and T f i+h are the throughfall observations located at i and i+h, respectively, and N(h) is the number of pairs separated by the lag distance h.The variogram is characterized by three parameters: sill, range, and the nugget.The sill is the limit of the variogram with h tending to infinite.The range is the lag distance at which the variogram reaches the sill.The range describes the length scale over which the observations are correlated; for this reason it is also known as spatial correlation length.The nugget is the discontinuity of the variogram at the origin.Spatial analysis requires a minimum number of observations to obtain a proper variogram (at least 100) [30].To overcome the lack of enough observations, Sterk and Stein [31] proposed a method to extend the analysis into the space-time domain and join data from multiple storms to map wind blow mass transport.The same methodology has been successfully applied to study spatial variability of Tf at fine-scale in Spain and Iran [32,33].Following Sterk and Stein [31], the first step of pooling the spatial variation of Tf is to standardize the observations to the overall mean with the following equation: in which T f i,j is the standardized Tf observation at the location i and for the period j, and T f j and µ are the period's mean Tf and the overall mean of the pooled observations, respectively.The variograms were then calculated using the standardized data of several periods: in which T f i,j and T f i+h,j denote the pair of standardized observations, separated by the distance h, for the period j.Using this formulation, each pair of observations within the same period contributes to the estimation of the standardized variogram, irrespective of the time that the measurements were made.
We pooled Tf observations from 4 periods (Prd-3, Prd-4, Prd-5, and Prd-6) in spring, and from 3 periods (Prd-7, Prd-9, and Prd-11) in the summer season.Prd-8 and Prd-10 were not included in summer pool, because the first one presented a strongly skewed distribution (out of the [−1 1] range suggested by [34]) and the second presented a relatively low P G that was not comparable with the other pooled periods.
Equation (3) was applied to generate two experimental standardized variograms, one for the spring and another for the summer season.We used a maximum lag distance of 36 m, that is, one half of the largest distance within the plot [21].Experimental standardized variograms were fitted using 4 standard variograms models widely used to represent the spatial distribution of throughfall (namely exponential, Gaussian, spherical, and pure nugget), i.e., [21,35,36].Models were evaluated at different lag sizes ranging from 1 to 4 m at intervals of 0.1 m.To select the best fitted model, the criteria of the smallest sum of square errors was used [34].The model with the minimum sum of square errors Water 2018, 10, 317 7 of 15 and its respective lag size were considered for further calculations.For the data analysis we used R software, version 3.2.3(R Development Core Team, 2015).Additionally, gstat package [37] was used for geostatistical analysis.
Each of the two fitted models was converted to 4 and 3 period-variogram models for spring and summer, respectively.To achieve this, we considered that the range of spatial dependence of Tf was independent of the P G , which means that same range was inferred for all the periods of the same season.The sill parameter (C) and the nugget effect (C 0 ) were estimated using the overall mean µ and the Tf mean per period, T f j , with the following equations: in which C s and C s 0 are the sill and nugget effect parameters of the standardized variograms, respectively.
Maps of Tf were created using ordinary kriging with the corresponding period-variogram.We assessed the period models using a leave-one-out cross-validation.The cross-validation consists of predicting values of throughfall at each location i by means of ordinary kriging but leaving out the observed value on the location i.From the cross validation we estimated the mean absolute errors (MAE) and the relative MAE (RMAE = MAE/T f j ).To depict the spatial variation of Tf in the plot, we created maps of Tf as percentage of rainfall (Tf-p) by dividing the kriged Tf map to the respective P G of the period.

Temporal Persistence of Throughfall
Time stability plots (TSP) were used to investigate the temporal persistence of throughfall cf.[1,20].
To make a comparison among different periods, we calculated a "normalized" throughfall T f i,j cf.[20] as: in which T f i,j is Tf collected at location i during period j, and T f j and S j are the respective mean and the standard deviation calculated for every period.We computed these values for i = 32 locations, and j = 5 stationary periods of the data collection (Prd-11 to Prd-15).We created TSP by ranking the values of normalized throughfall from the smallest to the largest, showing the relative deviation of Tf from the mean.TSPs describe two types of temporal persistence: extreme persistence and general persistence.Extreme persistence is identified in steep tails in the plot, and its presence indicates either sites persistently very wet or very dry [1,20].General persistence is described by the slope of the line in the middle quantiles and indicates whether the collectors are persistently wetter or drier than the mean.

Throughfall
Between 17 February 2015 and 2 November 2015, Tf was measured for 15 periods; we excluded Prd-1 and Prd-2 from the analysis, because some records of snowfall were detected in the nearest weather station in Deelen.During these periods (Prd-3 to Prd-15), a total of 572 mm of rain was recorded.Similarly, measured Tf totals for the same period was 350 mm, representing 61.2% of gross rainfall.Tf totals for the spring, summer, and autumn seasons were 33.9 mm (42% of P G ), 210.2 mm (67% of P G ), and 106.3 mm (60% of P G ), respectively.
The median CV for the 13 periods of throughfall collection was 16%.Variations in CV appeared to be mainly associated with the variations in P G .For example, Prd-9 had the highest Tf depth of 74.1 mm and had the lowest CV of 10.1%, while the Prd-15 had the lowest Tf depth of 2.6 mm and had the highest CV of 29.3% (Table 1).

Spatial Patterns of Throughfall
The similarity between the mean and a median values of Tf (Table 1) and estimated skewness coefficients within the limits [−1, 1] indicates that Tf observations exhibited a distribution close to normal.Only Prd-8 presented a skewness coefficient larger than 1; then, it was excluded from subsequent analysis.Distance to the nearest tree was not a good predictor for any period of measurements (R 2 < 0.2).
One hundred twenty five (spring) and 94 (summer) Tf (standardized) observations were used to estimate combined variograms.The best-fitting models (based on the sum of square errors) for the combined variograms were exponential for both seasons.The lag sizes that minimize the sum of square error in the fitting were 3 m and 2.7 m for spring and summer, respectively.The effective range parameter was about 12 m and 8 m for spring and summer variograms, respectively.Partial sill effect was detected only for spring season (1.5), and the sill parameter was 3.8 m 2 and 59.6 m 2 for spring and summer, respectively (Table 2, Figure 3).Water 2018, 10, x FOR PEER REVIEW 8 of 15

Throughfall
Between 17 February 2015 and 2 November 2015, Tf was measured for 15 periods; we excluded Prd-1 and Prd-2 from the analysis, because some records of snowfall were detected in the nearest weather station in Deelen.During these periods (Prd-3 to Prd-15), a total of 572 mm of rain was recorded.Similarly, measured Tf totals for the same period was 350 mm, representing 61.2% of gross rainfall.Tf totals for the spring, summer, and autumn seasons were 33.9 mm (42% of PG), 210.2 mm (67% of PG), and 106.3 mm (60% of PG), respectively.
The median CV for the 13 periods of throughfall collection was 16%.Variations in CV appeared to be mainly associated with the variations in PG.For example, Prd-9 had the highest Tf depth of 74.1 mm and had the lowest CV of 10.1%, while the Prd-15 had the lowest Tf depth of 2.6 mm and had the highest CV of 29.3% (Table 1).

Spatial Patterns of Throughfall
The similarity between the mean and a median values of Tf (Table 1) and estimated skewness coefficients within the limits [−1, 1] indicates that Tf observations exhibited a distribution close to normal.Only Prd-8 presented a skewness coefficient larger than 1; then, it was excluded from subsequent analysis.Distance to the nearest tree was not a good predictor for any period of measurements (R 2 < 0.2).
One hundred twenty five (spring) and 94 (summer) Tf (standardized) observations were used to estimate combined variograms.The best-fitting models (based on the sum of square errors) for the combined variograms were exponential for both seasons.The lag sizes that minimize the sum of square error in the fitting were 3 m and 2.7 m for spring and summer, respectively.The effective range parameter was about 12 m and 8 m for spring and summer variograms, respectively.Partial sill effect was detected only for spring season (1.5), and the sill parameter was 3.8 m 2 and 59.6 m 2 for spring and summer, respectively (Table 2, Figure 3).The quality of the period-variograms was better for the summer periods than for the spring periods.We found MAE ≤ 1.5 mm for the spring periods and MAE < 7 mm for summer periods.The RMAE were larger for the spring periods (~20%) than for the summer periods (~10%) The quality of the period-variograms was better for the summer periods than for the spring periods.We found MAE ≤ 1.5 mm for the spring periods and MAE < 7 mm for summer periods.The RMAE were larger for the spring periods (~20%) than for the summer periods (~10%) Frequent occurrence of rainfall events is expected to enhance interception loss (and reduce Tf ) due to frequent canopy wetting and drying.This effect can be observed in Figure 4, although total P G was larger in Prd-4; Tf-p was larger for Prd-6.This could be explained largely by the number of rainfall events during those periods (4 in Prd-6 vs. 7 in Prd-4; Table 1).However, this effect is not consistent in other periods, and considerable scatter was evident (Table 1).This could be attributed to the differences in other factors (e.g., rainfall intensity, wind speed) among the periods.Frequent occurrence of rainfall events is expected to enhance interception loss (and reduce Tf) due to frequent canopy wetting and drying.This effect can be observed in Figure 4, although total PG was larger in Prd-4; Tf-p was larger for Prd-6.This could be explained largely by the number of rainfall events during those periods (4 in Prd-6 vs. 7 in Prd-4; Table 1).However, this effect is not consistent in other periods, and considerable scatter was evident (Table 1).This could be attributed to the differences in other factors (e.g., rainfall intensity, wind speed) among the periods.

Temporal Persistence of Throughfall
Temporal persistence of Tf was analyzed for Prd-11 to Prd-15 (stationary position).The time stability plot (Figure 5) showed a tail, indicating persistence of extremely dry collectors.Funnel collectors #19 and #22 were positioned at extremely very dry sites, while none of the collectors was positioned at any extremely wet site.
The general persistence is related to the slope of the line in the middle quartiles.A gentle slope indicates that Tf observations are not different from the mean.The general persistence was moderate considering that 15 collectors were not significantly different from the mean (t-test, α = 0.05), 3

Temporal Persistence of Throughfall
Temporal persistence of Tf was analyzed for Prd-11 to Prd-15 (stationary position).The time stability plot (Figure 5) showed a tail, indicating persistence of extremely dry collectors.Funnel collectors #19 and #22 were positioned at extremely very dry sites, while none of the collectors was positioned at any extremely wet site.
The general persistence is related to the slope of the line in the middle quartiles.A gentle slope indicates that Tf observations are not different from the mean.The general persistence was moderate considering that 15 collectors were not significantly different from the mean (t-test, α = 0.05), 3 observations were excluded from the analysis due to disturbances, and 9 collectors were significantly wetter and 5 significantly drier than the mean (t-test, α = 0.05).
From field observations, we found that some of the collectors located close to the tree stem have low Tf.However, the linear correlation between the distance to the nearest tree (DNT) and the averaged normalized Tf was weak (Figure 6).The observations indicated that over a distance of 0-1 m to the nearest tree, two of the three extremely dry collectors were found, while in the same range of distance one wet collector was found.The driest collector was located at a distance of 1.3 m from the nearest tree.
From field observations, we found that some of the collectors located close to the tree stem have low Tf.However, the linear correlation between the distance to the nearest tree (DNT) and the averaged normalized Tf was weak (Figure 6).The observations indicated that over a distance of 0-1 m to the nearest tree, two of the three extremely dry collectors were found, while in the same range of distance one wet collector was found.The driest collector was located at a distance of 1.3 m from the nearest tree.From field observations, we found that some of the collectors located close to the tree stem have low Tf.However, the linear correlation between the distance to the nearest tree (DNT) and the averaged normalized Tf was weak (Figure 6).The observations indicated that over a distance of 0-1 m to the nearest tree, two of the three extremely dry collectors were found, while in the same range of distance one wet collector was found.The driest collector was located at a distance of 1.3 m from the nearest tree.

Throughfall Variability
The average variability in Tf (~16%) for the fifteen periods was about 15% less than the values reported by Bouten, et al. [7] for the same plot when the forest was denser (885 trees ha −1 ) and shorter (mean height ~20 m) than at the time of the present study.Likewise, the average CV found in the present study was ~5% lower than the values reported by Raat, et al. [1] for the same plot.Effects of forest thinning or seasonal changes in stand structures on spatial variability of throughfall (CV) are documented by few other studies [14,15].For example, Sun, et al. [15] reported a significant decrease in CV of Tf in Japanese cypress plantation because of strip thinning.Similarly, Staelens, et al. [14] found a significantly higher mean CV of Tf in a leafy period than that observed in a leafless period.One possible explanation for the reduction in throughfall CV in the studied plot is that thinning practice occurred in winter 1995-1996.After forest thinning or canopy pruning, openness increases, which, in turn, leads to decrease in canopy storage capacity (i.e., lower rainfall is needed to stabilize the Tf rate) and increase in throughfall rate cf.[15].
Another possible explanation for the lower CV of Tf found in the present study as compared to the other two studies can be related to the amount of throughfall collected across the periods.In many studies it was found that aggregation of small events can cause higher variation in Tf than an aggregation of large events [1,6,7,10].The latter is closely related to the temporal resolution of the measurements.Staelens, et al. [14] reported that the CV decreased with increasing length of the sampling interval (i.e., increase in the amount of throughfall received).A similar effect was observed in the present study: in periods with a higher amount of rainfall, CV ranged from 10% to 17%, whereas in periods with a lower amount of rainfall, CV values ranged from 28.8% to 29.3% (Table 1).The weekly period of Tf observation used by Bouten, et al. [7] (as opposed to ~15 day intervals in the present study) with lower accumulated values of Tf may thus result in higher CV.The sampling interval in the present study was similar to that used by Raat, et al. [1], but their study was concentrated in the months with low rain intensities.The CV found by Raat, et al. [1] was indeed lower than the one estimated by Bouten, et al. [7], and higher than the CV in the present study.

Spatial Patterns of Throughfall
Based on the skewness coefficient, we found that most of our periods of measurements presented Gaussian distribution (Table 1).Many of the cases found in the literature are event-based studies; they showed skewed frequency distributions of throughfall.However, there are also some studies, mainly in temperate forest, that show a Gaussian behavior, e.g., [17].In a maritime pine stand Loustau, et al. [17] reported normal distribution of throughfall (52 collectors) for most of 32 individual storms.Likewise, in two (young and old) Douglas-fir stands Keim, et al. [20] found non skewed distributions of Tf (94 collectors) 8 storms (P G > 19 mm).
Similarly, skewness evaluated in our study was within the normality limits for all but in one period (Prd-8) (Table 1).Such normality in the observations distribution seems to be an effect of the aggregation period used for the measurement of Tf.A similar finding was reported by Zimmermann and Zimmermann [6] in their study in a 12-year old teak plantation in Panamá (i.e., octile skew range [−0.2, 0.2] for weekly measurements).
Previous descriptions of throughfall patterns through variograms have shown contrasting results (see comparisons in [21,35]).The study of Bouten, et al. [7] at Speulderbos forest with 36 collectors and weekly measurements reported that spatial patterns of throughfall around trees can be determined effectively through variograms, but they did not report the obtained parameter values.This is an indicator that, at least in a younger state of the forest and before thinning, the spatial correlation of throughfall can be determined.In a storm based study over two Douglas-fir stands in the Pacific Northwest, Keim, et al. [20] found no strong but evident spatial correlations of throughfall equal to 5 m for both, a mid-age (60 year, crown size of 5 m), and an old (600 years, crown size of 10 m) stand, respectively.However, they did not fit a variogram model, and the reported range was inferred from the experimental variogram plots.
Voss, et al. [35] reported that for tropical forest in Panama, the accuracy of the variogram parameters is largely influenced by the number gauges used.In addition, the very low spatial variability of Tf (Table 1) found at our study site indicates that the forest type (e.g., mono species or species-rich) and presence or absence of understory shrubs play an important role in the redistribution rainfall in the forest floor.Therefore, the number of throughfall gauges to capture the Tf spatial variability properly as recommended by Voss, et al. [35] may not necessarily be applicable in all the forest types.Moreover, Voss, et al. [35] further recommended that the extent of the plot should exceed several times the range, and our study fulfilled that requirement.
Two other studies at fine scale have effectively used the methodology proposed by Sterk and Stein [31] to study spatial distribution of Tf.Gómez,et al. [33] detected spatial correlation structure in Tf distribution in mature olive trees.Also, Fathizadeh, et al. [32] reported a spatial correlation length of Tf measurements of approximately 5-6 m beneath Persian oak trees.However, Fathizadeh, et al. [32] found high values RMAE (~40% to 100%) and attributed these to the lack of observational points at short distance.They recommended that random relocation of the gauges between collection periods would reduce the uncertainty.
The application of Sterk and Stein [31] methodology allowed us to generate information about spatial distribution of Tf.For instance, spring periods showed a homogeneous distribution of Tf for the Prd-4 (Figure 4a) while it was more variable during the Prd-6 (Figure 4b).Similar analysis can be done for summer periods in case similar spatial distribution of Tf-p for both periods (Prd-9 and Prd-11) can be observed (Figure 4c,d).During summer, it seems that the rainfall duration is large enough to saturate the canopy and produce large Tf-p of ~70%.Also, in the map of Prd-11, there are several blue spots indicating larger Tf-p (~80%).Those can also be largely due to the dynamic interplay between event frequency, duration, and intensity vis á vis canopy storage capacity within the period of measurement.
Our study supports the findings of Keim, et al. [20], who detected spatial correlation structure of Tf beneath Douglas-fir canopy in a temperate climate.The lack of large number of points can be overcome by pooling Tf observations from periods of similar characteristics.Our study had the advantage of random relocation of gauges between collection periods that facilitated the observational points at short distance.The relatively low RMAE (ranging from 8% to 22%) found in our study could also be due to the longer temporal scale of sampling (~biweekly).

Temporal Persistence of Throughfall
The majority of studies on persistence of throughfall have reported a systematic temporal trend (i.e., persistently very wet or dry) over time [1,14,15].The persistence of wet and dry area may well affect the spatial patterns of soil moisture in the forest floor, thereby affecting, for example, root water uptake [7], nitrification-denitrification rates [1,38] and soil micro-flora and fauna distributions [39].
The results of the temporal stability showed throughfall in the present study was (slightly) more stable than it was ~25 years ago [1].Although the differences can partly be related to the sampling designs, the presence of extremely dry collectors in this study can largely be attributed to the change in canopy structure due to the enlargement of the living crown.Moreover, we investigated whether there is any relationship between distance to the nearest tree and the amount of Tf.We did not find a clear correlation between the distance and the normalized Tf (Figure 6), although extreme dry collectors (#9 and #12) were located at less than 1.5 m to the nearest tree.The literature is not conclusive about this relationship.Some authors found evidence of this relationship [8,40] and some others reported from weak to no evidence [17,18,20,41].In their analysis of throughfall before and after forest thinning, Sun, et al. [15] hypothesized that lower throughfall recorded by some of the gauges close to the tree after thinning was due to enhanced interception loss through the enhanced ventilation after thinning.However, evidence supporting this hypothesis is quite scattered and contradictory.For example, Keim, et al. [20] reported few persistently wet gauges close to tree stems in a 60-year old conifer, whereas the result was opposite in an uneven age old conifer stand.Likewise, Sato, et al. [42] reported large spatial variability of throughfall close to trees in the eucalyptus plantation.Such results suggest that not only the distance to the stems but also the foliage aggregation above the positioned gauges could play a larger role in the Tf spatial variability.

Conclusions
The study examined spatial variability of throughfall in Douglas-fir stand near the settlement Garderen, the Netherlands.Over the 10 months of the study period, measured throughfall was 431 mm (61% of P G ).The spatial variability of Tf was lower than that reported by similar studies.This can largely be attributed to the homogeneity of the stand, the absence of understory, and, to some extent, the time scale of aggregation used in the study.The spatial correlation lengths detected in spring and summer season were about 12 m and 8 m, respectively.By using ordinary kriging, relative Tf (throughfall normalized by incoming rainfall) maps were created for 4 periods.Such maps depicted the homogeneous distribution of Tf within the studied plot, indicating comparatively low spatial variability of Tf in the study area.Throughfall measurements were further analyzed using time stability plots that showed, in general, a persistent pattern with few noticeable very dry spots.The latter was probably due to relatively high canopy density above the collectors.

Figure 1 .
Figure 1.Study area in Speulderbos: (a) location map of the study site, (b) top view of the Douglasfir canopy, (c) funnel-type collectors used for throughfall measurements.

Figure 1 .
Figure 1.Study area in Speulderbos: (a) location map of the study site, (b) top view of the Douglas-fir canopy, (c) funnel-type collectors used for throughfall measurements.

Figure 2 .
Figure 2. Throughfall sampling scheme.The black crosses (+) represent the marked grid spaced at 8 m, on the x and y direction.Black triangles represent trees with symbol size scaled to DBH.The black circles represent the position of the funnel-type collectors in one period (Prd).Each collector was located by randomly selecting an azimuth (Az) between 0 and 360°, and a radius distance (r), between 0 and 4 m (see upper right corner).

Figure 2 .
Figure 2. Throughfall sampling scheme.The black crosses (+) represent the marked grid spaced at 8 m, on the x and y direction.Black triangles represent trees with symbol size scaled to DBH.The black circles represent the position of the funnel-type collectors in one period (Prd).Each collector was located by randomly selecting an azimuth (Az) between 0 and 360 • , and a radius distance (r), between 0 and 4 m (see upper right corner).

Note: a
Model applied to fit the variogram: Exponential (Exp); b Effective ranges (a) were calculated for exponential model (Effective range = range × 3).

Figure 3 .
Figure 3.Estimated standardized variograms of Tf for (a) spring season with fitted exponential model; (b) summer season with fitted exponential model.Each point was labeled with the respective number of pairs per lag distance.

Figure 3 .
Figure 3.Estimated standardized variograms of Tf for (a) spring season with fitted exponential model; (b) summer season with fitted exponential model.Each point was labeled with the respective number of pairs per lag distance.

Figure 4 .
Figure 4. Tf-p maps for (a,b) spring periods: Prd-4 and Prd-6; (c,d) summer periods: Prd-9 and Prd-11.Filled black triangles represent tree position; the size is proportional to the diameter at breast height, and black circles represent the funnel-type collectors.

Figure 4 .
Figure 4. Tf -p maps for (a,b) spring periods: Prd-4 and Prd-6; (c,d) summer periods: Prd-9 and Prd-11.Filled black triangles represent tree position; the size is proportional to the diameter at breast height, and black circles represent the funnel-type collectors.

Figure 5 .
Figure 5.Time stability plot of normalized Tf for 32 stationary funnel-type collectors.Black dots are values of normalized Tf for each funnel-type collector of the stationary periods (Periods 11-15).Black asterisks are the averaged values of normalized Tf.The numbers on the x-axis are field labels for the collectors sorted by lowest averaged normalized Tf value.

Figure 6 .
Figure 6.Normalized Tf versus distance to nearest tree.Black dots represent the normalized throughfall value for each funnel-type collector.Each symbol is labeled with the collector number.

Figure 5 .
Figure 5.Time stability plot of normalized Tf for 32 stationary funnel-type collectors.Black dots are values of normalized Tf for each funnel-type collector of the stationary periods (Periods 11-15).Black asterisks are the averaged values of normalized Tf.The numbers on the x-axis are field labels for the collectors sorted by lowest averaged normalized Tf value.

Figure 5 .
Figure 5.Time stability plot of normalized Tf for 32 stationary funnel-type collectors.Black dots are values of normalized Tf for each funnel-type collector of the stationary periods (Periods 11-15).Black asterisks are the averaged values of normalized Tf.The numbers on the x-axis are field labels for the collectors sorted by lowest averaged normalized Tf value.

Figure 6 .
Figure 6.Normalized Tf versus distance to nearest tree.Black dots represent the normalized throughfall value for each funnel-type collector.Each symbol is labeled with the collector number.

Figure 6 .
Figure 6.Normalized Tf versus distance to nearest tree.Black dots represent the normalized throughfall value for each funnel-type collector.Each symbol is labeled with the collector number.

Table 2 .
Parameters of the fitted standardized-variograms and converted period-variogram models.
Note: a Model applied to fit the variogram: Exponential (Exp); b Effective ranges (a) were calculated for exponential model (Effective range = range × 3).