How Do Ground Litter and Canopy Regulate Surface Runoff?—A Paired-Plot Investigation after 80 Years of Broadleaf Forest Regeneration

: Relatively minimal attention has been given to the hydrology of natural broadleaf forests compared to conifer plantations in Japan. We investigated the impacts of ground litter removal and forest clearing on surface runoff using the paired runoff plot approach. Plot A (7.4 m 2 ) was maintained as a control while plot B (8.1 m 2 ) was manipulated. Surface runoff was measured by a tipping-bucket recorder, and rainfall by a tipping-bucket rain gauge. From May 2016 to July 2019, 20, 54, and 42 runoff events were recorded in the no-treatment (NT), litter removed before clearcutting (LRBC), and after clearcutting (AC) phases, respectively. Surface runoff increased 4 × when moving from the NT to LRBC phase, and 4.4 × when moving from the LRBC to AC phase. Antecedent precipitation index (API 11 ) had a signiﬁcant inﬂuence on surface runoff in the LRBC phase but not in the NT and AC phases. Surface runoff in the AC phase was high regardless of API 11 . The rainfall required for initiating surface runoff is 38% and 56% less when moving from the NT to LRBC, and LRBC to AC phases, respectively. Ground litter and canopy function to reduce surface runoff in regenerated broadleaf forests.


Introduction
As of 2017, approximately 25 million ha (67%) of land area in Japan is forested. Plantation forests make up approximately 40% of the forested area while the remaining 60% is natural forests [1,2]. From 1966 to 2017, plantation forests have increased by 29% to 10.2 million ha while natural forests shrunk by 13% to 13.5 million ha [1]. While forested areas remain high (twice the world's average of 29%), the proportion of forest type (broadleaf, conifer, mixed, natural, and plantation forests) fluctuates according to market forces, situational demands, and governmental policies [1,3]. Conifer plantation forests are expected to expand, but instead, some conifer forests are being reverted to broadleaf forests [4]. At the same time, some broadleaf and mixed-forests are also being cleared to make way for development [5].
Since hydrology is highly influenced by forests, there is a need to understand how these forests (and their modification) affect runoff processes and water resources. In the past, there have been extensive studies carried out on hydrological processes in cypress plantations, including those investigating the hydrological impacts of forest conversion [6], forest litter [7,8], tree physiology [9,10], and various management practices [11,12]. The was carried out, resulting in the secondary mixed-forest cover (predominantly broadleaf) today [24,25]. However, due to their proximity to urban areas, more and more of such areas are being cleared to make way for development (i.e., residential, commercial, and industrial buildings) in recent years [4,5]. For this reason, it is also important to understand the impacts of clearcutting secondary broadleaf forests on surface runoff. In this study, we extend previous works [13,26,27] by observing at the plot scale, changes in surface runoff after litter removal and clearcutting. We hypothesise that surface runoff will increase with litter removal and clearcutting in a regenerated temperate broadleaf forest. At a larger scale, this may have implication on floods during storm events.

Study Area
Experimental plots were set up in the 13.9-ha Ananomiya Experimental Forest of the Ecohydrology Research Institute, the University of Tokyo Forests, in Seto city, Aichi prefecture ( Figure 1). The general area, traditionally known as the Eastern Owari Hills, is known for a long history of ceramic and pottery production. As a result of intensive timber harvesting to be used as kiln fuel, the area was characterised by denuded hills with bare, exposed subsoil in the past (Figures 2 and 3) [24,25]. Another round of intensive timber exploitation took place during the Meiji era , when industrialisation in Japan was at its peak [25]. High sediment yields and landslides had been a common phenomenon. With auspices from the government, the Ananomiya Experimental Forest (AEF) area was established in 1922 by the Tokyo Imperial University (University of Tokyo, today for the purpose of restoring forests and preventing erosion-related hillslope disasters. Construction of a weir dam was completed in 1924, and hydrological and meteorological observations commenced in 1925 [24]. Although this site was severely denuded in the past, forest cover has now recovered due to restoration-planting efforts (known as "sabo planting"), where Japanese red pine (Pinus densiflora) and Japanese black pine (Pinus thunbergii) were planted as the pioneer species [24,25].     The area has a warm temperate climate characterised by hot summers, moderate winters, and high annual rainfall that occurs mainly in two distinct wet periods-the Baiu season (May-June) with frequent and prolonged rain events; and the typhoon season (September-October) with high-intensity rainfall [27]. Annual rainfall is 1594 mm year −1 (averaged 1 May 2016-30 April 2019). The annual mean temperature is 12-15 °C, but the extremes may reach −12 °C (dawn, in winter) and 39 °C (midday, in summer). Daily averages of relative humidity range between 55 and 99.7% (data from on-site meteorological station). The dominant geology is a deep cretaceous layer overlain by weathered granite. Through soil core sampling around the plots followed by the jar test, percentage sand (15.03-20.61%), silt (26.22-39.05%), and clay (43.21-53.18%), show that the soil texture is of the clay type [29]. Bulk density range is 808.3-1069.2 kg m −3 (mean = 928.7 kg m −3 ).
At present, the area has a temperate mixed-forest composition with dominant canopy species comprising Quercus serrata and Ilex pedunculosa. The sub-canopy layer is dominated by Eurya japonica, Lyonia ovalifolia and Gamblea innovans; while the ground cover, by Sasa nipponica and Dicranopteris dichotoma. Although this is a secondary mixed-forest, the number of pioneer conifer trees (black pine and red pine) is small owing to species succession. Trees in and around the experimental plots are broadleaved; therefore, effectively making this study a comparison in a broadleaf forest. Stand density is 1925 stems/ha, which comprises canopy trees (>9 m in height; 725 stems/ha) and sub-canopy trees (3-9 m in height; 1200 stems/ha).   The area has a warm temperate climate characterised by hot summers, moderate winters, and high annual rainfall that occurs mainly in two distinct wet periods-the Baiu season (May-June) with frequent and prolonged rain events; and the typhoon season (September-October) with high-intensity rainfall [27]. Annual rainfall is 1594 mm year −1 (averaged 1 May 2016-30 April 2019). The annual mean temperature is 12-15 °C, but the extremes may reach −12 °C (dawn, in winter) and 39 °C (midday, in summer). Daily averages of relative humidity range between 55 and 99.7% (data from on-site meteorological station). The dominant geology is a deep cretaceous layer overlain by weathered granite. Through soil core sampling around the plots followed by the jar test, percentage sand (15.03-20.61%), silt (26.22-39.05%), and clay (43.21-53.18%), show that the soil texture is of the clay type [29]. Bulk density range is 808.3-1069.2 kg m −3 (mean = 928.7 kg m −3 ).
At present, the area has a temperate mixed-forest composition with dominant canopy species comprising Quercus serrata and Ilex pedunculosa. The sub-canopy layer is dominated by Eurya japonica, Lyonia ovalifolia and Gamblea innovans; while the ground cover, by Sasa nipponica and Dicranopteris dichotoma. Although this is a secondary mixed-forest, the number of pioneer conifer trees (black pine and red pine) is small owing to species succession. Trees in and around the experimental plots are broadleaved; therefore, effectively making this study a comparison in a broadleaf forest. Stand density is 1925 stems/ha, which comprises canopy trees (>9 m in height; 725 stems/ha) and sub-canopy trees (3-9 m in height; 1200 stems/ha). The area has a warm temperate climate characterised by hot summers, moderate winters, and high annual rainfall that occurs mainly in two distinct wet periods-the Baiu season (May-June) with frequent and prolonged rain events; and the typhoon season (September-October) with high-intensity rainfall [27]. Annual rainfall is 1594 mm year −1 (averaged 1 May 2016-30 April 2019). The annual mean temperature is 12-15 • C, but the extremes may reach −12 • C (dawn, in winter) and 39 • C (midday, in summer). Daily averages of relative humidity range between 55 and 99.7% (data from on-site meteorological station). The dominant geology is a deep cretaceous layer overlain by weathered granite. Through soil core sampling around the plots followed by the jar test, percentage sand (15.03-20.61%), silt (26.22-39.05%), and clay (43.21-53.18%), show that the soil texture is of the clay type [29]. Bulk density range is 808.3-1069.2 kg m −3 (mean = 928.7 kg m −3 ).
At present, the area has a temperate mixed-forest composition with dominant canopy species comprising Quercus serrata and Ilex pedunculosa. The sub-canopy layer is dominated by Eurya japonica, Lyonia ovalifolia and Gamblea innovans; while the ground cover, by Sasa nipponica and Dicranopteris dichotoma. Although this is a secondary mixed-forest, the number of pioneer conifer trees (black pine and red pine) is small owing to species succession. Trees in and around the experimental plots are broadleaved; therefore, effectively making this study a comparison in a broadleaf forest. Stand density is 1925 stems/ha, which comprises canopy trees (>9 m in height; 725 stems/ha) and sub-canopy trees (3-9 m in height; 1200 stems/ha).

Treatment and Data Collection
Two sloping experimental plots, A (7.4 m 2 ) and B (8.1 m 2 ), were established with a 4.2-m buffer in-between ( Figure 4). Boundaries of each plot were created by inserting hard plastic (roofing) sheets into the ground. The boundary was made contiguous to prevent leakage of surface runoff. The mean slope for plots A and B were 35 • and 28 • , respectively. At the bottom boundary of each plot, a gutter was installed to collect overland flow which was then channelled via a PVC pipe to an 'Uizin UIZ-TB200' tipping-bucket gauge. Throughout the study period, no treatment was applied to plot A (control plot) whereas litter removal followed by clearcutting were carried out in plot B.

Treatment and Data Collection
Two sloping experimental plots, A (7.4 m 2 ) and B (8.1 m 2 ), were established with a 4.2-m buffer in-between ( Figure 4). Boundaries of each plot were created by inserting hard plastic (roofing) sheets into the ground. The boundary was made contiguous to prevent leakage of surface runoff. The mean slope for plots A and B were 35° and 28°, respectively. At the bottom boundary of each plot, a gutter was installed to collect overland flow which was then channelled via a PVC pipe to an 'Uizin UIZ-TB200' tipping-bucket gauge. Throughout the study period, no treatment was applied to plot A (control plot) whereas litter removal followed by clearcutting were carried out in plot B. To ensure minimal disturbance to the soil surface, litter removal was carried out from outside the plot boundary and precautionary measures were taken to not step into the plot area. The dry weight of litter removed (including leaf, twigs, fruits, and seeds) was 9.48 kg (oven dried at 60 °C for 72 h). In the clearcutting phase, a total of 13 trees (species: Deutzia crenata, Eurya japonica, Gamblea innovans, Ilex pedunculosa, and Lyonia ovalifolia) with DBH ranging 0.8-3.5 cm as well as 10 small shrubs were removed from plot B. In plot A (control plot), six trees (species: E. japonica and L. ovalifolia) and small shrubs were left as they were. Although having fewer trees than in plot B, one of the trees (E. japonica) in plot A was multi-stemmed (nine stems), resulting in a canopy coverage that was more extensive than that of a normal tree. To minimise external environmental influence, areas adjacent to the plots were made to be as similar as possible to that inside the plots. Areas surrounding plot A (3 metres away from the top, left, and bottom boundary) were left forested. For plot B, a 1-metre buffer from the perimeter (plot boundary) was given the same treatment (litter removal and clearcutting) as inside the plot. Photographs of plot B in the NT, LRBC, and AC phases are shown in Figure 5.
At a nearby weather station (approximately 470 m westwards of the runoff plots), rainfall data was recorded at 5-minute intervals by an 'Ota Keiki OW-34-BP' tippingbucket rain gauge (0.5 mm/tip) connected to a 'Campbell Scientific CR10X' datalogger. Data collection, site observation, and maintenance were carried out in three phases: no treatment phase (NT) from May to October 2016, litter removed before clearcutting phase (LRBC) from March 2017 to June 2018, and after clearcutting phase (AC) from September 2018 to July 2019. To ensure minimal disturbance to the soil surface, litter removal was carried out from outside the plot boundary and precautionary measures were taken to not step into the plot area. The dry weight of litter removed (including leaf, twigs, fruits, and seeds) was 9.48 kg (oven dried at 60 • C for 72 h). In the clearcutting phase, a total of 13 trees (species: Deutzia crenata, Eurya japonica, Gamblea innovans, Ilex pedunculosa, and Lyonia ovalifolia) with DBH ranging 0.8-3.5 cm as well as 10 small shrubs were removed from plot B. In plot A (control plot), six trees (species: E. japonica and L. ovalifolia) and small shrubs were left as they were. Although having fewer trees than in plot B, one of the trees (E. japonica) in plot A was multi-stemmed (nine stems), resulting in a canopy coverage that was more extensive than that of a normal tree. To minimise external environmental influence, areas adjacent to the plots were made to be as similar as possible to that inside the plots. Areas surrounding plot A (3 metres away from the top, left, and bottom boundary) were left forested. For plot B, a 1-metre buffer from the perimeter (plot boundary) was given the same treatment (litter removal and clearcutting) as inside the plot. Photographs of plot B in the NT, LRBC, and AC phases are shown in Figure 5.
At a nearby weather station (approximately 470 m westwards of the runoff plots), rainfall data was recorded at 5-minute intervals by an 'Ota Keiki OW-34-BP' tipping-bucket rain gauge (0.5 mm/tip) connected to a 'Campbell Scientific CR10X' datalogger. Data collection, site observation, and maintenance were carried out in three phases: no treatment phase (NT) from May to October 2016, litter removed before clearcutting phase (LRBC) from March 2017 to June 2018, and after clearcutting phase (AC) from September 2018 to July 2019.

Data Processing and Analysis
The main mode of analysis is comparing storm surface runoff between the NT, LRBC, and AC phases. For this, we considered storm events to be separate if there were at least 6 hours of rain-free period in between events-a common rule used in the region. In all analyses, runoff in plot B was normalised by runoff in plot A. Erroneous data and periods of equipment malfunction were excluded from analysis. In total, 116 runoff events were available for analysis-20 in the NT, 54 in the LRBC, and 42 in the AC phases, respectively.
We explored the data in four different aspects: (i) overall runoff amount and duration, (ii) storm event surface runoff, (iii) the effects of antecedent precipitation on surface runoff, and (iv) rainfall threshold required for the initiation of surface runoff (amount of rain from the start of a rain event to the first recorded surface runoff). For comparison between the treatment phases, the variable of interest (surface runoff, rainfall threshold) and their various influencing factors (treatment, average rainfall intensity, antecedent precipitation index) were fitted via the generalised linear model (GLM). The best-fit GLM was chosen via the Akaike Information Criterion (AIC) and Bayes Factor in the software R v.4.0.3 (required packages: glm2, devtools, flexplot). Differences between phases were assessed via the Mann-Whitney U test.

Data Processing and Analysis
The main mode of analysis is comparing storm surface runoff between the NT, LRBC, and AC phases. For this, we considered storm events to be separate if there were at least 6 hours of rain-free period in between events-a common rule used in the region. In all analyses, runoff in plot B was normalised by runoff in plot A. Erroneous data and periods Before analysis, we ascertained the relationships between surface runoff (Q) and rainfall (P) as well as between surface runoff in plot B (Q B ) and surface runoff in plot A (Q B ). Strong positive relationships were found between Q B and Q A , especially in the NT phase (r 2 = 0.98), which signify suitability for a paired-plot study. After establishing a baseline relationship and ensuring that there were no anomalies, the latter regression was carried forward to be used in the paired-plot comparison.
Following this, mean rainfall intensity (P i ) was also tested for influence on runoff generation. However, it was not statistically significant (p > 0.05, r 2 ≤ 0.1); thus, not included for further analysis.
Because the data were non-normally distributed (Shapiro-Wilk, p < 0.05) and have unequal population variance (Levene, p < 0.05), Q B vs. Q A was provisionally test-fitted using the 'glm2' command (package: glm2) in R using several 'family' (gaussian, poisson, gamma) and 'link' (identity, log, inverse) functions. Both original and log-transformed data were tried. The 'gamma-identity' fit on original data was found to produce the best fit in the NT, LRBC, and AC phases, and was thus used.
Differences in runoff between NT, LRBC, and AC were assessed via the Mann-Whitney U test using calculated treatment effects (Te). To calculate Te, the first step was to establish a calibration equation (using the earlier gamma-identity GLM) between plot B and plot A in the NT phase: where Q NT obs B and Q NT obs A are event surface runoff observed in the NT phase in plot B and A, respectively; and a and b are regression coefficients. After treatment, runoff in plot B in the LRBC and AC phases were estimated using Equations (2) and (3), respectively: where Q LRBC est B and Q AC est B are 'estimated runoff' in plot B in the LRBC and AC phases, respectively, and Q LRBC obs A and Q AC obs A are 'observed runoff' in plot A in the LRBC and AC phases, respectively. Treatment effects (Te) in LRBC (Te LRBC/NT ) and AC (Te AC/NT ) with reference to the NT phase were calculated as follows: In the same way, with reference to the LRBC phase, Te in AC (Te AC/LRBC ) was calculated by first establishing a calibration equation between plot B and plot A in the LRBC phase (Equation (6)); then, using the equation coefficients (c and d) to estimate surface runoff in AC (Equation (7)), and subtracting the estimated values from observed values (Equation (8)):

Effects of Antecedent Precipitation on Surface Runoff
The first step was to assess the relationship between surface runoff and antecedent precipitation index (API), as well as to determine antecedent days with the highest influence. Several combinations of GLM family and link function were tried as with the earlier section. All equations were then inspected visually and compared based on AIC and Bayes Factor. Relationships were weak and monotonic instead of linear. API 11 were found to produce the highest correlation and best fit (AIC: 450.05 via 'gamma-log').
After confirming that API 11 has a negative relationship with Q, the GLM, Q B /Q A against API 11 , were fitted for each treatment phase (NT, LRBC, and AC): where e and f are equation coefficients.

Rainfall Threshold
The rainfall threshold required for the initiation of surface runoff (P T )-defined as the total rainfall from the start of an event to the first tip of the tipping-bucket dataloggerwere also compared between the different treatments. Before comparison, the influence of API 11 and P i on P T were assessed. Relationships were tested using data from individual plots and both plots collectively, in individual phases and all phases collectively, as well as every combination of phases and plots. However, no clear relationships were found. Therefore, the comparison of P T was conducted via the Mann-Whitney U Test.

Overall Surface Runoff and Precipitation
At the monthly timescale, surface runoff in plot B was higher than in plot A after treatment (LRBC and AC phase) ( Figure 6). This differed from the pre-treatment period where surface runoff in plot A was slightly higher. The normalised runoff ratio (C N )defined as the runoff ratio in plot B divided by runoff ratio in plot A-increased from 0.558 in NT to 2.647 and 11.875 in LRBC and AC, respectively ( Table 1). The standard deviation (SD) also increased, signifying increased variability (Table 1).
Bayes Factor. Relationships were weak and monotonic instead of linear. API11 were found to produce the highest correlation and best fit (AIC: 450.05 via 'gamma-log').
After confirming that API11 has a negative relationship with Q, the GLM, QB/QA against API11, were fitted for each treatment phase (NT, LRBC, and AC): where e and f are equation coefficients.

Rainfall Threshold
The rainfall threshold required for the initiation of surface runoff (PT)-defined as the total rainfall from the start of an event to the first tip of the tipping-bucket datalogger-were also compared between the different treatments. Before comparison, the influence of API11 and Pi on PT were assessed. Relationships were tested using data from individual plots and both plots collectively, in individual phases and all phases collectively, as well as every combination of phases and plots. However, no clear relationships were found. Therefore, the comparison of PT was conducted via the Mann-Whitney U Test.

Overall Surface Runoff and Precipitation
At the monthly timescale, surface runoff in plot B was higher than in plot A after treatment (LRBC and AC phase) ( Figure 6). This differed from the pre-treatment period where surface runoff in plot A was slightly higher. The normalised runoff ratio (CN)defined as the runoff ratio in plot B divided by runoff ratio in plot A-increased from 0.558 in NT to 2.647 and 11.875 in LRBC and AC, respectively ( Table 1). The standard deviation (SD) also increased, signifying increased variability (Table 1).   Besides runoff, the runoff frequency represented by the ratio of runoff frequency in plot B to runoff frequency in plot A (T B/A ) also increased with disturbance: 0.977, 2.292, and 7.250 in the NT, LRBC, and AC phases, respectively (Table 1). In the NT phase, surface runoff frequency was slightly higher in plot A as visualised in the flow-duration curve (FDC) in Figure 7. In the LRBC and AC phases, despite a decrease in surface runoff frequency in plot A, surface runoff frequency in plot B have increased. The FDC of plot B also became more convex, which is a sign of increased frequency of high-magnitude surface runoff. Besides runoff, the runoff frequency represented by the ratio of runoff frequency in plot B to runoff frequency in plot A (TB/A) also increased with disturbance: 0.977, 2.292, and 7.250 in the NT, LRBC, and AC phases, respectively (Table 1). In the NT phase, surface runoff frequency was slightly higher in plot A as visualised in the flow-duration curve (FDC) in Figure 7. In the LRBC and AC phases, despite a decrease in surface runoff frequency in plot A, surface runoff frequency in plot B have increased. The FDC of plot B also became more convex, which is a sign of increased frequency of high-magnitude surface runoff.

Event-Scale Surface Runoff in Different Phases
In all phases and in both plots, surface runoff (Q) increased linearly with event rainfall (P) (Figure 8). In the NT phase, strong linear relationships were found (plot A, R 2 = 0.9019; plot B, R 2 = 0.8488).

Event-Scale Surface Runoff in Different Phases
In all phases and in both plots, surface runoff (Q) increased linearly with event rainfall (P) (Figure 8). In the NT phase, strong linear relationships were found (plot A, R 2 = 0.9019; plot B, R 2 = 0.8488). Besides runoff, the runoff frequency represented by the ratio of runoff frequency in plot B to runoff frequency in plot A (TB/A) also increased with disturbance: 0.977, 2.292, and 7.250 in the NT, LRBC, and AC phases, respectively (Table 1). In the NT phase, surface runoff frequency was slightly higher in plot A as visualised in the flow-duration curve (FDC) in Figure 7. In the LRBC and AC phases, despite a decrease in surface runoff frequency in plot A, surface runoff frequency in plot B have increased. The FDC of plot B also became more convex, which is a sign of increased frequency of high-magnitude surface runoff.

Event-Scale Surface Runoff in Different Phases
In all phases and in both plots, surface runoff (Q) increased linearly with event rainfall (P) (Figure 8). In the NT phase, strong linear relationships were found (plot A, R 2 = 0.9019; plot B, R 2 = 0.8488).  Q B = 9.3967·Q A + 0.5503 (AC phase; R 2 = 0.4388, p < 0.01) (12) From the LRBC to AC phase, surface runoff increased 4.3 times. These increases are in the range of that in the monthly timescale. When expressing in terms of treatment effects (Te) the LRBC and AC phases have Te values of 1.341 and 2.589, respectively. From LRBC to AC Te is 1.843 (Table 1). Te ( Figure 10) were statistically significantly different between al phases (Mann-Whitney, p < 0.01). In addition, the y-intercept increased by two orders when moving from NT to LRBC and AC, signifying earlier commencement of surface runoff.  The strong correlation between plot A and B in the NT phase (R 2 = 0.9223) indicates suitability for a paired-plot comparison. Based on regression coefficients, surface runoff in the LRBC and AC phase is 4.4 and 18.7 times higher, respectively, than in the NT phase. From the LRBC to AC phase, surface runoff increased 4.3 times. These increases are in the range of that in the monthly timescale. When expressing in terms of treatment effects (Te), the LRBC and AC phases have Te values of 1.341 and 2.589, respectively. From LRBC to AC, Te is 1.843 (Table 1). Te ( Figure 10) were statistically significantly different between all phases (Mann-Whitney, p < 0.01). In addition, the y-intercept increased by two orders when moving from NT to LRBC and AC, signifying earlier commencement of surface runoff. QB = 2.1971·QA + 0.3336 (LRBC phase; R 2 = 0.3630, p < 0.01) QB = 9.3967·QA + 0.5503 (AC phase; R 2 = 0.4388, p < 0.01) The strong correlation between plot A and B in the NT phase (R 2 = 0.9223) ind suitability for a paired-plot comparison. Based on regression coefficients, surface run the LRBC and AC phase is 4.4 and 18.7 times higher, respectively, than in the NT p From the LRBC to AC phase, surface runoff increased 4.3 times. These increases are range of that in the monthly timescale. When expressing in terms of treatment effect the LRBC and AC phases have Te values of 1.341 and 2.589, respectively. From LRBC t Te is 1.843 (Table 1). Te ( Figure 10) were statistically significantly different betwe phases (Mann-Whitney, p < 0.01). In addition, the y-intercept increased by two orders moving from NT to LRBC and AC, signifying earlier commencement of surface runo  Figure 11 shows surface runoff against API 11 in plots A and B in different treatment phases. In plot A (control), surface runoff varies little with API 11 throughout all phases whereas in plot B (treatment), the effects of API 11 (negative relationship) increased from NT to LRBC to AC.  Figure 11 shows surface runoff against API11 in plots A and B in different treatment phases. In plot A (control), surface runoff varies little with API11 throughout all phases whereas in plot B (treatment), the effects of API11 (negative relationship) increased from NT to LRBC to AC. Naturally, surface runoff is directly proportional to rainfall and inversely proportional to API. In our study, the relationship QB/QA vs. API11 is only significant in the LRBC phase (QB/QA = −0.014861·API11 + 3.863; p>0.01).

Rainfall Threshold for Surface Runoff Generation
The amount of rainfall required for surface runoff generation (rainfall threshold, denoted by "θPSR") was also evaluated ( Figure 12). With decreasing interception (NT to LRBC to AC phase), the rainfall threshold decreased accordingly (Figure 12). Pi and API did not have significant influence on rainfall threshold; thus, rainfall threshold between the different phases was assessed via the Mann-Whitney U test. Rainfall threshold in plot B normalised by that in plot A (θPSR_B/θPSR_A) is 38% lower when moving from the NT to LRBC phase. When moving from NT to AC, and LRBC to AC, θPSR_B/θPSR_A is 73% and 56% lower, respectively. Differences in rainfall threshold were statistically significantly different between the NT, LRBC, and AC phases (Mann-Whitney, p < 0.01).  Naturally, surface runoff is directly proportional to rainfall and inversely proportional to API. In our study, the relationship Q B /Q A vs. API 11 is only significant in the LRBC phase (Q B /Q A = −0.014861·API 11 + 3.863; p > 0.01).

Rainfall Threshold for Surface Runoff Generation
The amount of rainfall required for surface runoff generation (rainfall threshold, denoted by "θP SR ") was also evaluated ( Figure 12). With decreasing interception (NT to LRBC to AC phase), the rainfall threshold decreased accordingly (Figure 12). Pi and API did not have significant influence on rainfall threshold; thus, rainfall threshold between the different phases was assessed via the Mann-Whitney U test. Rainfall threshold in plot B normalised by that in plot A (θP SR_B /θP SR_A ) is 38% lower when moving from the NT to LRBC phase. When moving from NT to AC, and LRBC to AC, θP SR_B /θP SR_A is 73% and 56% lower, respectively. Differences in rainfall threshold were statistically significantly different between the NT, LRBC, and AC phases (Mann-Whitney, p < 0.01). Figure 11 shows surface runoff against API11 in plots A and B in different tr phases. In plot A (control), surface runoff varies little with API11 throughout al whereas in plot B (treatment), the effects of API11 (negative relationship) increas NT to LRBC to AC. Naturally, surface runoff is directly proportional to rainfall and inversely tional to API. In our study, the relationship QB/QA vs. API11 is only significant in th phase (QB/QA = −0.014861·API11 + 3.863; p>0.01).

Rainfall Threshold for Surface Runoff Generation
The amount of rainfall required for surface runoff generation (rainfall thresh noted by "θPSR") was also evaluated ( Figure 12). With decreasing interception LRBC to AC phase), the rainfall threshold decreased accordingly ( Figure 12). Pi did not have significant influence on rainfall threshold; thus, rainfall threshold the different phases was assessed via the Mann-Whitney U test. Rainfall threshol B normalised by that in plot A (θPSR_B/θPSR_A) is 38% lower when moving from th LRBC phase. When moving from NT to AC, and LRBC to AC, θPSR_B/θPSR_A is 7 56% lower, respectively. Differences in rainfall threshold were statistically sign different between the NT, LRBC, and AC phases (Mann-Whitney, p < 0.01).

No-Treatment Phase
In the NT phase, surface runoff increased linearly with precipitation. This follows findings of existing studies that found similar relationships at the catchment scale in a secondary mixed-forest [13]. Antecedent moisture did not significantly affect surface runoff in the NT period, which may have been due to more important natural factors (rainfall, canopy, ground litter) that govern surface runoff in natural forests. Rainfall threshold for runoff generation was the highest in the NT phase, which reflects good interception properties (by ground litter and canopy) compared to in the treatment phases. These characteristics serve to establish a baseline for post-treatment comparisons.

Litter Removed, before Clearcutting Phase
In the LRBC phase, surface runoff was higher than in the NT phase (Figure 9), which demonstrated the importance of ground litter in regulating surface runoff [30]. Although the loss in ground litter interception is perceived to be the immediate cause, other mechanisms may have also contributed. In particular, the loss of protective cover may have caused soil compaction and reduced permeability when rainfall and concentrated throughfall directly impacts the soil surface; additionally, the increased floor evaporation may have resulted in drier and hydrophobic conditions. Both of these reduce infiltration and promote rapid direct runoff [20,31,32]. Compared to the straightforward increase in effective rainfall, soil compaction could have played a more important role because canopy in broadleaf forests are known to concentrate rainwater into larger droplets with higher momentum; thus, compacting the soil surface upon impact [19,31]. This increase in surface runoff is consistent with an experiment at the catchment scale [13]. This present study, however, differed in terms of spatial and temporal scales. Therefore, the results are not directly comparable. A more similar plot study was conducted using simulated rainfall in northwest Beijing, China, but with litter from different tree species (Quercus variabilis and Pinus tabulaeformis) [31]. Our results agree with theirs, but the degree of runoff increase is different. Surface runoff in our study increased over four times as opposed to only 43% in Li et al. [31]. This may be attributed to different species, litter morphology, soil type, and rainfall characteristics [22]. Miyata et al. [19] observed that overland flow increases with silt and clay content, and have hypothesised that this is due to silt and clay being key components in the formation of soil crust. This could be a possible mechanism in the present study considering the high clay content in our plots. As for species and litter morphology, broadleaf litter are known to have higher water storage capacity than conifer litter [33]. The increased flow frequency may have resulted from the detection of throughfall of small events that were previously effectively intercepted by ground litter (higher effective rainfall).
Besides an increase in the amount and duration of surface runoff, sensitivity towards API 11 also increased ( Figure 11). Although dry soils enhance hydrophobicity, resulting in increased surface runoff in forested environments, the effects of this mechanism may have been masked by the presence of ground litter in the NT phase [11,[34][35][36]. In addition, the absence of ground litter in the LRBC phase may have accelerated soil drying due to enhanced floor evaporation; thus, creating hydrophobic conditions that promote direct runoff [33]. Therefore, for a given value of API 11 , actual soil moisture could have been drier in LRBC compared to in the NT phase.
The rainfall threshold required for surface runoff generation was 37.55% lower in LRBC (Figure 12), which reflects the buffering properties of ground litter interception. This is also a sign of increased water repellency usually associated with dry and compacted soils. This indicated that the suggested runoff-on-litter mechanism [20,34], if ever occurred, did not play a significant hydrological role in the NT phase. Li et al. [31] also reported similar results, whereby shorter time was required to initiate surface runoff in bare plots as opposed to litter-covered plots. They also observed that broadleaf litter was marginally more effective in delaying the generation of surface runoff.

After Clearcutting Phase
Surface runoff was the highest in the AC phase. The increase in surface runoff (4× from LRBC to AC) demonstrated the role of forest canopy in intercepting and evaporating rainwater in secondary broadleaf forests ( Figure 9) [14,37,38]. Compared to other hydrological losses, canopy interception (12-24% throughout Japan) may have a distinctively large influence as it directly controls effective rainfall [14,21,39]. Studies from other regions also recorded similar values: 21.3% of rainfall in a tropical broadleaf forest in Kalimantan [40]; 10% and 20% in a primary and regenerated tropical rainforest, respectively, in Sabah, Malaysia [41]; and 22.4% in a conifer forest in California [42]. Our findings differed from what Nanko et al. [43] have suggested, whereby throughfall in forests have larger droplets and higher momentum that compacts the soil and promotessurface runoff (corresponding to LRBC in this study). The increase in flow frequency compared to in the AC phase ( Figure 7) could have been caused by similar mechanisms as when transitioning from the NT to LRBC phase, which are increased amount and duration of effective rainfall. Besides a shift in the position of the FDC, its shape has become more convexed; therefore, signifying an increase in the frequency of large-magnitude runoff events (despite plot B having a lower slope angle). In contrast to the LRBC phase, where soil compaction by throughfall was assumed to be the main cause of increased runoff, increases in the AC phase is attributed to increased effective rainfall due to the removal of canopy.
Surface runoff in AC was generally higher than that in the NT and LRBC phases regardless of API 11 ( Figure 11). This may be caused by high water repellency as a result of increased solar radiation, hence the drying-up of soils and formation of crust. Additionally, the effects of removing the canopy and increasing effective rainfall may have been disproportionately dominant to the extent of exempting influence from soil moisture. Although some studies have attributed this to the latter [42,44,45], at present, we suggest causality to the former based on findings in existing studies [46,47]. The findings of Daikoku et al. [48]vapour pressure deficit and the below-canopy available energy strongly controls forest floor evaporation in Japan-partly support this. Although saturated soils may also override the influence of API, such conditions never occurred in our well-drained sloping plots.
Compared to the NT and LRBC phases, the decrease in rainfall threshold required for surface runoff generation demonstrated the dominant influence of canopy interception ( Figure 12). Although stemflow is known to be a major mechanism that channels canopy-intercepted water to the ground [49][50][51][52], our data has shown that its influence was overridden by the effects of loss in canopy interception. Instead of ground factors, we attribute changes in the AC phase to the loss in canopy interception because of two reasons: (i) following common logging practices, clearcutting has left the lower portion of trunks and roots untouched, hence ground roughness remained unchanged at least in the medium term before root decay takes place; (ii) prior to clearcutting, litter removal in the LRBC phase have discounted the possibility of changes in ground litter being a factor. Therefore, we are able to deduce that it is the increased effective rainfall from the removal of canopy that has caused enhanced surface runoff (amount, duration, and initiation threshold) in the AC phase.

Additional Consideration and Suggestions
The heterogenous rainfall regime in the region resulted in different rainfall input and moisture conditions throughout the year, which may in turn affect surface properties (such as hydrophobicity) that govern surface runoff. Although surface runoff data of different treatments were collected in different phenological periods, differences in rainfall and moisture conditions were accounted for by the control plot and should not affect the results [13,53,54]. Although effort was made to establish the plots to be as similar as possible, there were small differences in tree species, slope, and plot size. Being situated in a natural forest, it was impossible to establish both plots with the exact tree species and size. The uneven microtopography and undulations on the ground resulted in a 0.7 m 2 difference in plot size. This should not affect the results and conclusion because comparisons were performed between different treatment (phases) instead of between different plots. Plot A merely served as a control to account for environmental conditions. Rainfall intensity is known to significantly affect surface runoff generation in various land cover and regions [55][56][57][58]. The reason for the lack of influence from rainfall intensity in our plot study requires further investigation, but this may be provisionally attributed to surface runoff being naturally high in the area as well as the lack of complex hydrological pathways and sinks found in catchment scale studies [8,59].
This study has quantified the role of ground litter and canopy on surface runoff generation at the plot scale. Results from this study are not directly extrapolatable to the hillslope, catchment, and basin scale due to various hydrological and landscape factors. At the hillslope scale, species composition in a patch area may govern canopy and litter interception. At the catchment scale, various factors such as subsurface flow, groundwater, sinks, and hydrological connectivity need to be accounted for. At the basin scale, land-use and regional weather may be the governing factors instead [22]. In agreement with Liu et al. [22], future studies should cover more species, regions, and spatial scales-larger scales to understand environmental impacts; smaller scales to understand ecohydrological mechanisms.
Although we have investigated changes in surface runoff from litter removal and clearcutting, other hydrological components (infiltration, preferential flow, soil moisture) are still less understood. Investigating these (via moisture sensors, dyes, etc.) will give insights on the below-ground hydrology as well as verify the possible mechanisms that have been discussed in this article.

Conclusions
Due to the scarcity of information on the role of ground litter in regenerated broadleaf forests and the effects of clearcutting on surface runoff, we conducted a paired-plot experiment and discussed the possible hydrological mechanisms.
The absence of ground litter increased surface runoff by up to four times. Antecedent moisture had a significant influence on surface runoff generation after ground litter was removed. Without ground litter, 38% less rainfall was required to initiate surface runoff.
Clearcutting increased surface runoff by another four times when compared to the litter-removed period. Without both ground litter and the canopy, antecedent moisture no longer affect surface runoff. Canopy loss resulted in 56% less rainfall required to initiate surface runoff when compared to the ground litter-removed phase.