Factors controlling throughfall in a Pinus tabulaeformis forest in North China

The factors that control throughfall in Pinus tabulaeformis plantations were investigated using linear and curve analyses based on direct measurements of rainfall, throughfall and stemflow from 36 rainfall events. The results showed the following: (1) there was significant spatial heterogeneity in throughfall rates in P. tabulaeformis plots; (2) the throughfall rate increased with increasing rainfall; and (3) the rate of increase gradually decreased. When rainfall reached approximately 25 mm, the throughfall rate stabilized. The coefficient of variation of the throughfall rate decreased with increasing rainfall, with a peak at approximately 10 mm of rainfall. The coefficient of variation of throughfall stabilized at 20%, and the coefficient of variation of the throughfall rate stabilized at 17%. A linear regression equation (R2 = 0.76) was derived by fitting the P. tabulaeformis average diameter at breast height (DBH), average tree height, average branch height, stand density, canopy thickness, canopy density, and the rainfall and throughfall rate. A highly positive correlation was found between the throughfall rate, canopy density, rainfall class and tree height (P < 0.01). By establishing a quadratic response surface model of the stand structure indicators and the throughfall rate, R2 was increased to 0.85 (P < 0.01). The quadratic regression analysis demonstrated a highly positive correlation between throughfall rate, canopy density and rainfall class.

(1) Monitor rainfall, TF, stemflow and interception loss in study area. We expect in-depth understanding of TF of P. tabulaeformis plantations and the relationship between the four in a variety of circumstances.
(2) Identify, classify and quantify the factors relevant to TF. In order to distinguish the single factor heterogeneity, we refine each factor and especially discriminate the heterogeneity of throughfall in detail under different rainfall class.
(3) Compare all relevant factors, filter and establish model with mathematical method. We expect the most accurate model to determine the factors that have the most important impact on TF, by a variety of mathematical methods to calculate and verify.

Materials and Methods
Site description. The study site is located at the Mulan Forestry Management Bureau in Hebei Province (41°49′35.8″N, 117°35′31″E; elevation 1270 m), which lies in the cross connecting area of the Yinshan Mountains and Swallow Mountain(s) in the upper reaches of the Luanhe River, China. The study site is a typical rocky mountainous area of northern China. Meteorological records indicate that the long-term mean annual air temperature is 3.2 ± 0.3 °C, and February and August are the coldest and warmest months, respectively. During the study period, the average temperature was obtained from meteorological records. The main tree species were Larix principis-rupprechtii, P. tabulaeformis, Betula platyphylla and Populus davidiana; the understory vegetation consisted of Rhododendron micranthum, Spiraea fritschiana, Carex lancifolia, Phlomis umbrosa and Thalictrum minus.
Instrumental setup and data collection. Six standard P. tabulaeformis sample plots were established in the study area, and the hydrological characteristics of these permanent sample plots were investigated from 2010 to 2011. The P. tabulaeformis stand density in the six plots constituted a gradient (2438 plants/hm 2 , 1815 plants/ hm 2 , 1500 plants/hm 2 , 1080 plants/hm 2 , 928 plants/hm 2 and 650 plants/hm 2 ). The investigation included tree, shrub and grass surveys. Table 1 shows the plot characteristics, including the plot area, mean diameter at 1.3 m above the ground (DBH), mean tree height, canopy density, stand density, and mean branch height (Table 1).
A HOBO U30 weather station was established in the wilderness area outside the plots. Air temperature, air relative humidity, wind velocity, solar radiation and rainfall were measured and recorded every 10 minutes. In this experiment, throughfall was collected in water reservoirs constructed of iron sheets with a collection area of 0.156 m 2 (0.2 × 0.78 m), which was equivalent to five standard rain gauges. An opening was made at the bottom and connected to 10-L plastic barrels using a plastic tube.
The canopy density was measured using a digital plant canopy imager CI-110 (USA). Fish-eye photographs were captured above the centre of each water reservoir at 4-9 points, which were evenly distributed in each plot. The photographs were taken 0.8 m above the throughfall collectors on cloudy days to avoid the effects of direct sunlight. Sky image pixels (a) and total image pixels (A) were used to remove the trunk from the image, and Photoshop CS software was used to process the binary image as follows: Canopy density (f) was calculated following (1), and the mean value for each measurement point was calculated. In addition, two water reservoirs placed outside the forest on open ground, at a distance of approximately 50-100 m, served as controls.
The water reservoir devices were located based mainly on the canopy density in the sample plots (Equation 1). To ensure the measurement data were representative of the entire plot, we chose areas with dense canopy, sparse canopy and forest gaps to establish the water reservoirs because some plots contained relatively higher canopy densities. However, for stands with a relatively lower canopy density, we chose standard trees with complete canopies, as well as gap centres, to position the water reservoirs. In this experiment, there were three water reservoirs in sample plot No 1, five in No 2, three in No 3, three in No 4, three in No 5, and three in No 6 ( Fig. 1).
At the same time, 10 × 10 m quadrats were established around each water reservoir in the plots. The height of the water reservoirs was at least 30 cm with an approximately 1° angle relative to the earth's surface. The average tree height, average DBH, average branch height, average canopy thickness, canopy density and stand density were carefully investigated within each 10 × 10 m quadrat.
When considering rainfall rate and stand structure (average tree height, average DBH, average branch height, average canopy thickness, canopy density and stand density), there is a need to determine which factor is the most important. To this end, stepwise regression analysis and a quadratic polynomial equation were adopted. Then, the accuracy of each equation was compared, and a decision was reached regarding which factor was most suitable. Reanalysis of datasets. The SPSS stepwise regression analysis method was used, and we established the optimal linear regression equation between throughfall and stand structure variables (average tree height, average DBH, average branch height, average canopy thickness, canopy density and stand density) to determine the fit of the parameters and throughfall. To improve simulation accuracy, quadratic polynomial equations were used to determine a quadratic response surface model for throughfall and stand structure variables. Regarding the difference in the analysis of rainfall levels, SAS software was used to perform significant difference testing and sorting of the rainfall levels associated with the throughfall rates. First, the Shapiro-Wilk normality test and Levene's homogeneity of variance test were carried out. When the results complied with the normality and homogeneity of variance of the data, we then carried out analysis of variance and Duncan's multiple comparisons. If any term did not match, the Kruskal-Wallis nonparametric test was applied using NPAR WAY. Based on the results of these tests, an assessment was made to identify significant differences in the rainfall level associated with the throughfall rate. If there were no significant differences in the majority of quadrats, it was assumed that the throughfall rate showed little change, and the rainfall level could be marked at the same level in the regression analysis. Based on general rainfall level division standards and data testing and rearrangement, a new level of rainfall was obtained. We then established new rainfall levels ranging from light to heavy to perform the regression analysis.

Results
Rainfall characteristics during the experimental period. To reduce error, we also observed rainfall at six locations. Rainfall at these locations was measured in an open area approximately 10 m from the plots. During the experimental period, effective rainfall was detected 36 times from July 2010 to September 2011. During the study period, there were six times when rainfall exceeded 20 mm and five times when rainfall intensity exceeded 10 mm · d −1 (Fig. 2). Figure 3 shows the relationship between throughfall and rainfall for P. tabulaeformis in the study area. The figure illustrates that rainfall and throughfall have a highly positive correlation (R 2 = 0.98) (Fig. 3a). The throughfall rate increased with rainfall, but the rate of increase gradually decreased and stabilized at approximately 75% (Fig. 3c). The interception rate stabilized at approximately 25% in this region. In addition, the throughfall rate and its coefficient of variation decreased as rainfall increased. The speed was initially fast and then became slow, with a turning point at approximately 10 mm (rainfall). The coefficient of variation of throughfall steadied at approximately 20% (Fig. 3b), and the coefficient of variation of the throughfall rate steadied at approximately 17% (Fig. 3d).

Rainfall level evaluation.
In this research, the variability of observational data at the 5-mm level of rainfall penetration was tested during the study period. Rainfall was lower than 50 mm·d −1 in the observed area. According to the rainfall level division standards of the China Meteorological Administration, daily rainfall in the 0~10 mm range represents light rain, 10~25 mm represents moderate rain, and 25~50 mm represents heavy rain. To determine whether there were differences in the throughfall rates of light rain, the light rain stage division was set to 5 mm, followed by a difference analysis. The same procedure was performed for moderate rain. With respect to throughfall associated with heavy rain, the rate mostly stabilized and was no longer divided ( Table 2).
The analysis results showed the following ( Table 2): throughfall rates in 17 of 20 small quadrats of P. tabulaeformis showed no significant differences when the rainfall accumulation was 0~10 mm. Therefore, we believe that throughfall rate in the P. tabulaeformis forest at the 0~10 mm·d −1 rainfall level was not significantly altered. Thus, this rainfall level (0~10 mm·d −1 ) was assigned a value of 1 to unify logistics data for regression analysis. We also found that 12 of 20 small quadrats exhibited no significant differences when the rainfall accumulation was 10~25 mm · d −1 ,  representing the vast majority of events. Thus, rainfall accumulation in the 10-25 mm · d −1 class was assigned a value of 2. Similarly, rainfall accumulation in the 25-50 mm · d −1 class was assigned a value of 3 (Table 3).
Relationship between throughfall and related factors. The throughfall rate of artificial P. tabulaeformis plantations was designated the dependent variable. Stand structure factors and rainfall factors closely related to the throughfall rate were chosen as independent variables ( Table 3).
The throughfall associated with each rainfall level and the corresponding stand structure indicator data were analysed using stepwise regression (Table 4), and a linear regression model for the throughfall rate was fitted. The throughfall rate model was significant (P < 0.0001), and the F value was high ( Table 4). The main factors affecting the throughfall rate were rainfall level (X 7 ), canopy density (X 6 ) and average tree height (X 2 ). When the stepwise regression model included only the rainfall level (X 7 ), the multiple correlation coefficient R 2 was 0.55, and when canopy density (X 6 ) was introduced, R 2 reached 0.74. The results indicate that the rainfall level (X 7 ) and canopy density (X 6 ) contributed significantly to the multiple correlation coefficient; they were the main factors that affected the throughfall rate. Finally, when average tree height (X 2 ) was introduced, R 2 only increased to 0.76, indicating that the average tree height was not a significant factor affecting the throughfall rate (P = 0.11); it was the second most important factor. Thus, the optimal linear regression model for the throughfall rate (TF r ) of the P. tabulaeformis forest is as follows (Equation 2): 3 83 X 106 07 X 11 57 X 82 33 (2) r 2 6 7 The linear stepwise regression equation resulted in a multiple correlation coefficient (R 2 ) of 0.76 (Fig. 4a), but the relative error mean was 14.89%. The standard deviation was 11.65%, and the relative error of the larger observation sample size accounted for 30% (Fig. 4c). Clearly, it was not the best model.
To generate a more accurate relationship, we fitted those factors that most significantly affected the throughfall rate -canopy density (X 6 ) and rainfall class (X 7 ) -using a quadratic response surface model ( Table 5). The model multiple correlation coefficient increased to 85% (Fig. 5a), and the average relative error decreased to 11.74%. The standard deviation was 9.01% (Fig. 5c) The distribution of a random variable always involves a measure of randomness. When the residual plot is spread around e = 0, there is no bias. Cases where the standard deviation of the residual is not uniform represent a violation of homoscedasticity, and parametric regression cannot be used. This finding indicates that the model was not appropriate. The approximated residuals of the throughfall rate are scattered points and present a Λ-shaped distribution (Fig. 5b), indicating that the model was not suitable. By contrast, the residual plots of the quadratic response surface exhibit a random distribution. Thus, this model was relatively ideal (Fig. 5).

Discussion and Conclusion
Throughfall and rainfall. In this study, the throughfall rate of P. tabulaeformis increased with increasing rainfall, and the growth rate decreased gradually with increasing rainfall. When rainfall reached 25 mm, the throughfall rate was essentially stable, in agreement with the results of previous studies. CV TF was strongly affected by rainfall at levels lower than 20 mm and was fairly stable when rainfall exceeded 20 mm (Fig. 4) Table 3. Throughfall rate for each rainfall class in each quadrat and corresponding stand structures in the study area. Note: Number of 10 × 10 m quadrats in the plots.
of decrease was initially fast and then became slow, and the turning point occurred at approximately 10 mm. This phenomenon was observed because at the beginning of a rainfall event, the majority of water saturates the canopy. The throughfall was initially low and then gradually increased before stabilizing after the canopy was completely saturated (Fig. 4). The average CV Tf values, i.e., 11-70%, were within the range of CV Tf values (21-53%) measured in other forests and deciduous conifers. In this study, CV Tf was higher than that reported for broadleaf forests in Japan, i.e., 12% for a mixed white oak forest 36 . When rainfall was heavier, the average CV Tf values were higher; otherwise, the values were lower. These differences in the spatial variation of throughfall among stands can be expected due to differences in canopy species, canopy structure, density, spatial homogeneity and meteorological phenomena. Differences in experimental design, such as those related to the size of plots, number, size and spatial density of collectors, and sampling time scale, make it difficult to directly compare our results with those reported in previous studies. Nevertheless, it appears reasonable to expect a relationship between stand characteristics and the spatial variation of throughfall under the same climatic conditions.
Factors that impact throughfall. The canopy features affecting throughfall include the following: the shape and size of trees, the roughness and thickness of the canopy, average tree height, average DBH, average branch height, stand density, branching pattern, leaf angle and LAI. Some studies have shown that the understory variation of throughfall and rainfall exhibit a negative correlation [37][38][39][40] . Studies have shown no significant correlation between throughfall and LAI 41,42 . In this study, the throughfall rate had a highly positive correlation with canopy density and rainfall, in agreement with previous studies 43,44 . The multiple correlation coefficient (R 2 value) of the stepwise linear regression equation was 0.76. This finding indicates that canopy density, rainfall and height had important effects on throughfall. Interception losses almost always depend on rainfall type and other meteorological conditions 20 . Throughfall increases as canopy density decreases and increases with decreasing rainfall and height. Interception values increase with increasing canopy density and decrease with decreasing canopy density.
This study compared the relationship between throughfall and stand factors under different rainfall levels. A highly positive correlation between throughfall and rainfall level was found (Equations 2 and 3). At lower rainfall levels, rainfall interception increased rapidly with increasing rainfall. However, at higher rainfall levels, rainfall  Table 4.
Stepwise linear regression analysis for the throughfall rate of the P. tabulaeformis forest in the study area. interception increased slowly with increasing rainfall, and the coefficient of variation of throughfall (CV Tf ) decreased drastically, reflecting the limitations of canopy interception. A denser canopy resulted in larger interception and CV Tf, with higher rainfall levels, but the impact weakened gradually. When estimating the spatial variability of throughfall, the placement of water reservoirs was used as an index. The largest throughfall amounts were measured around the canopy edge, near the trunk, or midway between the trunk and canopy edge 45 . In this study, we chose a dense canopy, sparse canopy and forest gap to establish the water reservoirs. Doing so could decrease the variability of throughfall due to the position of reservoirs. The stepwise linear regression analysis equation indicated that there was a relationship between throughfall rate and tree height. Few studies exist on the relationship between throughfall rate and tree height. In this study, we found a limited relationship between the throughfall rate and tree height.  Table 5. Results of the response surface quadratic (nonlinear) model for the throughfall rate in the P. tabulaeformis forest in the study area.