Spatial Variability and Optimal Number of Rain Gauges for Sampling Throughfall under Single Oak Trees during the Leafless Period

This study examined the spatial variability of throughfall (Tf) and its implications for sampling throughfall during the leafless period of oak trees. To do this, we measured Tf under five single Brant’s oak trees (Quercus brantii var. Persica), in the Zagros region of Iran, spanning a six-month-long study period. Overall, the Tf amounted to 85.7% of gross rainfall. The spatial coefficient of variation (CV) for rainstorm total Tf volumes was 25%, on average, and it decreased as the magnitude of rainfall increased. During the leafless period, Tf was spatially autocorrelated over distances of 1 to 3.5 m, indicating the benefits of sampling with relatively elongated troughs. Our findings highlight the great variability of Tf under the canopies of Brant’s oaks during their leafless period. We may also conclude that the 29 Tf collectors used in the present study were sufficient to robustly estimate tree-scale Tf values within a 10% error of the mean at the 95% confidence level. Given that a ±10% uncertainty in Tf is associated with a ±100% uncertainty in interception loss, this underscores the challenges in its measurement at the individual tree level in the leafless season. These results are valuable for determining the number and placement of Tf collectors, and their expected level of confidence, when measuring tree-level Tf of scattered oak trees and those in forest stands.


Introduction
Tree canopies significantly affect the terrestrial hydrological cycle by modifying the distribution of water received during rainfall events. The partitioning of gross rainfall (P g ) into throughfall (T f ), stemflow (S f ), and rainfall interception (I) by the forest canopy layer has a major role in structuring the water budget of forest ecosystems [1][2][3][4]. The T f can be further divided into free throughfall, the water reaching the forest floor without coming into contact with vegetation, and release throughfall, which splashes or drips downward from the canopy [5]. Not only does the presence of trees affect the amount of water reaching the forest floor, but it also determines the spatial distribution of T f [6,7].
The spatial variability of T f potentially exercises considerable control of forest hydrology and biogeochemistry dynamics because it can lead to pronounced heterogeneity in edaphic conditions [8,9]. To optimize forest management in terms of the availability of water and nutrition in the soil, a better understanding of the spatial distribution of T f within forests and its controls is imperative [10].
To date, the spatial redistribution of water under single trees (i.e., 'open-grown trees', those developing sufficiently apart from other to form full canopies and maximum branching) remains understudied. As noted above, T f is spatially variable due to heterogeneity imposed by variation in canopy structure (e.g., foliage and branch cover) and the proximity of tree boles onto the free throughfall and/or canopy drip patterns [10]. Nonetheless, literature reviews in hydrology have documented positive [11,12], negative [13,14], and inconsistent relationships [15,16] between where T f is sampled and nearness to the tree bole; therefore, their spatial dynamics can vary substantially across sites and tree types and we cannot generalize these patterns. Many such studies were done in architecturally simple conditions, for example in even-aged stands, yet T f can be highly variable at small spatial scales, namely beneath single trees [16,17]. Although researchers have reported on T f 's variability in an oak forest [16,18], no study has yet focused on its spatial variability in a leafless period in single oak trees.
Beyond the implications of T f spatial heterogeneity for soil processes, this heterogeneity also challenges constraining measuring T f as well as estimating I [19][20][21][22]. Overall, two sampling strategies for measuring T f exist: stationary and roving methods. In the stationary (fixed setup) method, T f collectors are kept at fixed positions during the entire study period [23][24][25][26][27][28]. The roving collector method is a sampling strategy whereby the T f collectors are randomly relocated at regular or variable times [22,[29][30][31][32][33][34]. Fewer gauges are required for estimating T f with high statistical confidence and a low margin of error is desirable for minimizing both cost and effort [7,19,21,35], especially at the tree-based scale characterized by high spatial heterogeneity. Using fixed T f gauges are advantageous because there is no need to reposition them periodically; also, fixed gauges permit a detailed examination of the I process over the course of individual rainstorm events [36,37].
In many temperate regions, especially those with a Mediterranean climate, most of the rain falls when trees are leafless [38]. However, to our best knowledge, few reports studied the spatial variability of T f or estimated its required numbers of collectors to reliably estimate means in a deciduous stand during this leafless period. Hence, the two objectives of this field research, conducted with single Brantii oak trees (Quercus brantii var. Persica) were to (1) evaluate the T f sampling method (stationary method) for finding the optimal number of T f collectors, and (2) evaluate the spatial variability of T f under these sampled trees.

Study Area
The experimental work took place from October 2011 to March 2012 in an oak forest in western Iran (46 • 24 E, 33 • 37 N, 1383 m a.s.l). According to the leaves' phenology, one year was divided into two canopy development stages: the leafed period (from April to mid-October) and the leafless period (from mid-October to the end of March). At the research site, vegetation consists of sparse and scattered Brant's oak trees (Q. brantii var. Persica), whose understory is used for agroforestry. Long-term (1986Long-term ( -2017  Five single, even-aged, healthy, and mature Q. brantii var Persica trees with similar morphologies of tree height, diameter at breast height (DBH), crown projected area (CPA), and height under branch were randomly selected among the trees of similar size, all located in a 2.5-ha tract ( Table 1). The crowns of the five trees did not overlap with any neighboring trees.

Field Measurements
The P g was measured using six cylindrical plastic rain collectors, each 9 cm in diameter, placed in an adjacent clearing about 20 m from the sample trees (with no crown interaction). Rainfall totals were collected on a daily (24 h) basis. The experimental network consisted of a total of 16 T f collectors (identical to the P g collectors) at each of the five trees. These T f collectors were randomly placed in a radial layout centered on each tree stem; along with eight azimuths, two T f collectors were installed at two randomly determined distances (i.e., 'near' and 'far' from the tree's stem; Figure 1). This arrangement was kept throughout the entire experiment (i.e., fixed gauges for the T f sampling). The quantity of water collected from the P g and T f collectors was manually measured using a graduated cylinder, to an accuracy of 1 mL. The averaged value of the rainfall collectors installed in the open area and beneath the oak stand was used to calculate P g and T f , respectively. All the P g and T f collectors were emptied and dried after each rainstorm event. Continuous daily records of wind speed and wind direction were obtained from the Ilam Meteorological Station.  The spatial variability of T f was assessed by a variogram analysis, this being a quantitative descriptive statistic that characterizes the spatial continuity of data sets. Variograms are a useful tool for quantifying how similarity in a response variable relates to proximity to better understand the structure of spatial patterns [15,39]. It expresses continuity as the average of the squared difference between quantities measured at different locations.
The calculation of meaningful directional variograms when having a low number of field observations is reportedly impossible [40,41]. Consequently, this has precluded checks for anisotropy in the spatial variability's structure of the respective variables of interest. Nonetheless, here we sought to know whether an anisotropic spatial structure could be discernable in the T f process using observations from just 16 rain collectors per tree, and these were artificially increased based on the method of Sterk and Stein [40]. More details can be found in Fathizadeh et al. [7,16].

Minimum Number of Throughfall Collectors
We calculated the number of T f samplers of all individual trees required to determine leafless T f averages within acceptable error margins (5%, 10%, 15%, and 20% of the mean cumulative T f ) and three confidence levels (α = 0.01, 0.05, and 0.1), assuming the mean T f of all 16 samplers for each individual tree represented the true value. The minimum number of T f collectors (N min ) needed to estimate T f within a preset percentage of the mean (E) at the 90%, 95%, and 99% confidence interval can be estimated from a coefficient of variation of T f (CV t ), as follows [20]: where z c is the critical value of the 90%, 95%, and 99% confidence levels.

Statistical Analysis
The Kolmogorov-Smirnov test was used to assess the normality of the data set. The quality of each variogram was tested using a cross validation. For this, all the samples were excluded piecemeal (one each time) from the data set and iteratively re-estimated by kriging using the remaining samples. Then the observed data and predicted values were compared to evaluate the kriging results [42]. From the cross validation we estimated the mean absolute errors (MAE) and the relative MAE (RMAE = MAE/T f ). MAE or RMAE measure how close forecasts or predictions are to the eventual outcomes. Cross validation was performed for trees with relatively suitable T f autocorrelation structure. Accordingly, total T f patterns were created for those trees having robust fitted models. Meanwhile, residual sums of squares (RSS) were estimated as a measure of how well the model fits the variogram data. More information about the methods is available in Fathizadeh et al. [16].

Throughfall Characteristics
Over the six-month monitoring period, total gross rainfall (P g ) was 302.4 mm and mean T f , expressed on a crown-projected area basis was 259 mm (SD = 6.8 mm) for 24 rainfall events. Therefore, the average interception loss equaled 14.3% of P g . Snowfall did not occur. The mean P g per event was 12.6 mm (SD = 13.4 mm; SE = 2.7 mm), with high variation ranging from 0.7 to 47.3 mm (CV = 1.06). Mean T f depth per event was 10.8 mm (SD = 0.31 mm) or 75.8% as a proportion of P g (SD = 4.4%).

Throughfall Spatial Patterns
Fitted models in multiple directions using the geostatistical analysis of T f amounts resulted in undetectable anisotropy for trees A, B, and E, but zonal anisotropy was dis-cernable in the data for trees C and D, which the surface variogram maps show as well (Figures 2 and 3). Therefore, we assumed an isotropic correlation structure for the former trees. A spherical model fitted for trees A and B yielded a stable sill emerging at about 2.3 and 3.5 m, respectively. The range over which spatial dependence is apparent was estimated to be~1 m for an exponential model fitted for tree E. The other two trees (i.e., C and D) were fitted with linear models and these had an unclear stable sill ( Figure 3, Table 2). A, B, and E) and anisotropic (for trees C and D) variograms for throughfall (mm) and fitted models in multiple directions for trees C and D. Details for fitting models are available in Fathizadeh et al. [16]. Table 2. Key characteristics of the variogram models for individual oak trees. RSS: residual sums of squares that provide a measure of how well the model fitted the variogram data; A 0 : range parameter; A 1 : the range parameter for the major axis of variation; A 2 : the range parameter for the minor axis; C 0 : nugget effect, and the ratio between structural variance (C) and sill (C 0 + C) parameters. * isotropic model, ** anisotropic model.  Cross validation was performed for five single trees and all 24 rainfall events during the leafless period, producing MAE and RMAE (Table 3). According to the results of previous analysis, described in Fathizadeh et al. [16], only those rainfall events with an RMAE <30% were used to generate the spatial patterns of T f via kriging ( Figure 4). Evidently, there were patchy spatial distributions of T f under trees A, B, and E (i.e., trees with an autocorrelation structure), its values being generally lower near the tree trunk. When quantitatively estimated for all trees, the average cumulative T f was 235.4 mm (SD = 24.3 mm) or 77.8% of P g for collectors near the tree trunk, 17.6% less than average cumulative T f (288.7 mm, or 95.5% of P g ; SD = 19.9 mm) for collectors positioned close to the crown edge (i.e., far from the tree trunk; Table 4). Except for tree B (t-test, F = 1.271, t(14) = −2.436, p < 0.05), the average cumulative T f of collectors near the crown edge and tree bole were not significantly different (t-test; α = 0.05). Considering that all five trees experienced the same meteorological conditions, it seems that tree E has a less physically heterogynous crown than do the others, in that most of its canopy has redistributed T f values for 50-80% of P g (Figure 4). The T f areas exceeding 100% were observed under all single trees, but presentable only for trees A, B, and E via kriged maps, mostly at the crown margins ( Figure 4). The spatial variability of T f can be expressed using the coefficient of variation (CV), i.e., the standard deviation as a proportion of the mean. The mean CV of cumulative T f amount was 25%, ranging from 17.9% for tree E to 28.8% for tree B (Table 3).  for trees A, B, and E, calculated as the ratio between the sum of individual T f events with a relative mean absolute errors (RMAE) <30% and total rainfall during those events, multiplied by 100. At the rainfall-event level, the spatial CV of T f had a tendency of decreasing with an increasing rainfall amount, with a median value of 33.7% ( Figure 5). During the study period, the CV of T f declined from a median value of 72.1% for a 2.1-mm rainfall event to stabilize at a median value of 32.2% for those events with at least 10 mm rainfall for the single studied trees.

Figure 2. Experimental isotropic (for trees
Forests 2021, 12, x FOR PEER REVIEW 1 At the rainfall-event level, the spatial CV of Tf had a tendency of decreasing w increasing rainfall amount, with a median value of 33.7% ( Figure 5). During the period, the CV of Tf declined from a median value of 72.1% for a 2.1-mm rainfall ev stabilize at a median value of 32.2% for those events with at least 10 mm rainfall single studied trees.

Minimum Number of Throughfall Collectors
We present the results of Equation (1) used to calculate the minimum Tf col needed to keep the relative error of the mean below 20%, 10%, or 5%, for the 24 o data sets. The average number of Tf collectors required to ensure the relative error is 20% varies between 5 (α = 0.10) and 14 (α = 0.01). The corresponding numbers for a r error ≤ 15% were 9 (α = 0.10) and 25 (α = 0.01); likewise, the number of collectors n to estimate the mean Tf ± 10% were 20 (α = 0.10) and 55 (α = 0. 01), and for a mean T they were 79 (α = 0.10) and 220 (α = 0. 10). According to Eq. (1), as an example, for and an acceptable error of 10% of the mean Tf, the required number of collectors w mated to be 24 (α = 0.05) (Table 5, Figure 6). These results suggested that 13 collect sufficient to estimate mean cumulative Tf with an error of 15% and a confidence in of 95% (Table 5). The improvement in the mean and CV of Tf volume (as explained error of mean Tf) arising from additional collectors is shown in Figure 6. The mo matic improvement occurred during the addition of the first 10 collectors (Figure 6 erally, once about 20 collectors are used, the relation between the number of collecto the mean Tf error plateau, achieving a near-zero slope ( Figure 6). Hence, above this ber, little improvement would be achieved by adding more collectors to the samplin work beneath trees. Moreover, and perhaps not surprisingly, we found that more

Minimum Number of Throughfall Collectors
We present the results of Equation (1) used to calculate the minimum T f collectors needed to keep the relative error of the mean below 20%, 10%, or 5%, for the 24 original data sets. The average number of T f collectors required to ensure the relative error is below 20% varies between 5 (α = 0.10) and 14 (α = 0.01). The corresponding numbers for a relative error ≤ 15% were 9 (α = 0.10) and 25 (α = 0.01); likewise, the number of collectors needed to estimate the mean T f ± 10% were 20 (α = 0.10) and 55 (α = 0. 01), and for a mean T f ± 5% they were 79 (α = 0.10) and 220 (α = 0. 10). According to Equation (1), as an example, for tree A and an acceptable error of 10% of the mean T f , the required number of collectors was estimated to be 24 (α = 0.05) (Table 5, Figure 6). These results suggested that 13 collectors are sufficient to estimate mean cumulative T f with an error of 15% and a confidence interval of 95% (Table 5). The improvement in the mean and CV of T f volume (as explained via the error of mean T f ) arising from additional collectors is shown in Figure 6. The most dramatic improvement occurred during the addition of the first 10 collectors (Figure 6). Generally, once about 20 collectors are used, the relation between the number of collectors and the mean T f error plateau, achieving a near-zero slope ( Figure 6). Hence, above this number, little improvement would be achieved by adding more collectors to the sampling network beneath trees. Moreover, and perhaps not surprisingly, we found that more T f collectors are needed for trees with a higher CV (tree B, C, and D) ( Figure 6, Table 5).   99  90  95  99  90  95  99  90  95  99  90  95  99   5  64  95  181 102 150 284  96  142 271  92  135  256  40  58  110  10  16  24  46  26  38  71  24  36  68  23  34  64  10  15  28  15  8  11  21  12  17  32  11  16  31  11  15  29  5  7

Throughfall Spatial Pattern
The partitioning of rainfall measured during the six months under single oak trees in the leafless period (T f = 85.7% of P g ) differed considerably from our previous research on the same trees in their leafed period (T f = 68.9% of P g ) [16]. The spatial patterns of T f , as quantified by its spatial continuity measured with variograms clearly varied among the five trees. This could be attributed to the differences in canopy structure factors (e.g., Wood Area Index, branching angle, canopy width) among the trees. The cumulative T f under single oak trees was spatially autocorrelated up to a range of about 1-3.5 m in the leafless period, overlapping to some extent with the 3-5 m that has been observed for the leafed period [16]. Similar variogram studies of small-scale spatial T f autocorrelation are scarce; for example, 1-3 m for the isolated olive trees in Spain [41] and 3-4 m for an individual deciduous beech tree in Belgium [43], and likewise at the stand scale; for example, 6-7 m for a 120-year-old beech forest in Luxembourg [44]; 2-6 m for a 42-year-old mixed-hardwood forest [45]; 2.6, 5.3, and 3.9 m in a teak plantation, young secondary forest, and an old secondary forest, respectively [46]; and 4-6 m in a rubber plantation [47]. However, T f in an old-growth tropical wet forest in Costa Rica indicated an autocorrelation range of 43 m because of the high structural variability of the forests caused by large areas dominated by either individual tree canopies or by treefall gaps [48]. A spatial correlation distance of 12 m and 8 m was detected for spring and summer, respectively, in a Douglas-fir stand in the Netherlands [49], whereas no spatial autocorrelation was observed for T f in a holm-oak forest [50]. As we mentioned in our previous study [16], differences in the spatial correlation of T f can be ascribed to the differences among canopy species and in canopy structure, density, spatial homogeneity, and meteorological phenomena [15]. Staelens et al. [43] related the small spatial range of T f to no spatial dependence existing in these studies because the minimum sampling distance exceeded the possibly present spatial range of T f variability. In other words, using a higher minimum intermediate distance of T f collectors that fails to sample at short lags can result in a lack of information on spatial autocorrelation at scales below 0.5 m. Close to the tree trunks, more branches and leaves intercepted the T f , which resulted in its lower value at this position. Conversely, those sub-canopy areas where T f surpassed incident rainfall were observed in spatial patterns at trees A, B, and E (Figure 4), mostly lying close to the crown edge and main branches, where intercepted water is apt to get channeled downward. This result could also be partly attributed to extreme fog events, wherein fog entrapment by the canopy leads to an enhanced T f volume [51,52]. The fog during the leafless period is very common due to the Mediterranean climate, therefore, the hydrological importance of fog should be considered in future hydrological studies.
For T f volumes, their spatial CV was 25% for oak single trees in the leafless period, similar to that previously observed for their leafed period (CV = 26.7%) [16]. Compared with our study, beneath the crown of one dominant beech (Fagus sylvatica L.) tree in Belgium, Staelens et al. [43] found less spatial CV (8%) for T f during its leafless period. Spatial variations in T f may be explained by canopy cover structure and wind combined [8,15,18,29,43].
Indeed, the high spatial variability of T f in the leafless period in our study could be explained by possible effects from non-vertical rainfall and a non-homogeneous distribution of raindrops driven by turbulent air flow, both above and within the canopy, because of local windy conditions during winter [53]. Due to the presence of wind, raindrops fall inhomogeneous with a tree canopy. Moreover, T f variability has a tendency of decreasing with an increasing P g ( Figure 5) through a constant median value of around 32% for P g ≥ 10 mm. Our findings agree with previous studies reporting that the CV diminishes with a greater magnitude of rainfall events [21,[53][54][55]. A plausible explanation is that the canopy becomes saturated when rainfall is more than 10 mm, such that any extra rainfall is mostly shunted to T f with no interception by the tree canopy [29,33,56,57].

Minimum Number of Throughfall Collectors
It is hard to estimate T f accurately using a minimum number of T f collectors because of difficulties in achieving a high spatial representation of T f measurements [58,59]. To obtain a proper estimate of the average T f , many collectors would be required [15,60], which entails considerable time and effort. Therefore, it is necessary to try to estimate what would be a suitable number of T f collectors for reliably estimating T f at the scale of a single tree. As indicated in Figure 6 and Table 5, if the goal is for 95% of resampled means to lie within a 10% error of the overall mean T f , 29 collectors should be used. Some scholars suggest that a lower number of collectors could be used to estimate the volume of T f , albeit accepting a higher margin of error [61,62]. Our resampling results indicate a progressive narrowing of the variability of cumulative mean T f as the number of collectors used increases. This result is similar to the findings of Puckett [61], who arrived at 11 for a minimum number of collectors to produce mean T f values limited to a 10% error and a 95% confidence interval. In our work with 29 collectors, adding one more collector would reduce considerably the variability of the T f means. Therefore, an effort to implement more collectors in this interval represents a considerable improvement in the determination of the mean cumulative T f . There are many reasons why different minimum numbers of T f collectors have been reported across research sites. Apart from the influence of tree species' traits, a difference in the T f collectors' sizes is likely also a major factor influencing the inferred number of collectors [63].
The optimal design of T f collectors for estimating the mean T f depends on the temporal scale of sampling, as well as forest complexity [46]. The measurement of T f is a critical phase in the evaluation of I because a significant error in the T f term results in corresponding over-or under-estimates of I. Based on our results for Brant's oak during its leafless period, we now also know how an individual crown contributes to lateral water translocation within the tree crown at fine spatial scales. This study can contribute to making an informed determination of the number of how many T f collectors would be needed for the robust estimates of mean T f and their associated error. One of the most important implications of our results concerns the determination of canopy ecohydrological parameters. In most studies that model I, parameters such as canopy saturation point, canopy storage capacity, free throughfall, and evaporation rate during rainfall are derived by plotting T f against P g (reviewed by Friesen et al. [64]). Hence, the accuracy of these derived canopy ecohydrological parameters depends on the accuracy of the P g and T f measurements from which they were derived.

Conclusions
Although within the range reported by similar studies, the spatial variability of throughfall water (T f ) investigate here using the CVs under single Brant's oak trees did not match up well, implying significant seasonal differences between the leafed and the leafless periods. The present study confirms that tree canopies modify the spatial distribution of T f reaching the forest floor during their leafless period. The T f quantity for single oak trees during their leafless period differed significantly from that during their leafed period, probably because of non-vertical rainfall profiles and the non-homogeneous distribution of raindrops above and within the canopy during winter. However, the spatial CV for T f volumes during the leafless period was 25%, similar to the leafed period. Mostly, those T f areas exceeding incident rainfall were found close to the crown edge, where the downfacing branches at the crown edge probably channeled the intercepted water to the forest floor. Conversely, near the tree trunks is where more intercepted water was transported to stem and converted to stemflow rather than T f , leading to lower T f values there.
We conclude that using 29 T f collectors is adequate for estimating the T f value with an accepted error of 10% at the 95% confidence level. The very important tree characteristic for considering T f spatial distribution is the characteristic of the canopy above the observed point, so we recommend measuring this parameter in future research. To gain further insight, long-term studies employing different T f collectors and sizes are now needed to draw comprehensive conclusions.

Data Availability Statement:
The data presented in this study are available on request from the first author.