Application of the Reformulated Gash Analytical Model for Rainfall Interception Loss to Unmanaged High-Density Coniferous Plantations Laden with Dead Branches

Abstract

Interception loss (IL) by the forest canopy removes a substantial quantity of rainwater within forested ecosystems. The large-scale unmanaged Japanese coniferous plantations with high stand density (SD) in Japan raise concerns about an additional increasing IL as a result of a new influential factor of dead branches under canopies. Thus, evaluating the usage of IL estimation models is vital to regulating the water and environment in such coniferous plantations. This study aimed to examine the applicability of the reformulated Gash analytical model (RGAM) to unmanaged coniferous plantations with high SD laden with dead branches. We established two plots (P1 and P2) laden with dead branches under the same SD of 2250 stems ha⁻¹ but with different numbers of dead branches (56 vs. 47 branches per tree) in an unmanaged Japanese coniferous plantation. Results demonstrated that a large difference was found in canopy storage capacity (S) in P1 and P2 (3.94 vs. 3.25 mm), which was influenced by the different number of dead branches; therefore, the IL ratio to gross rainfall differed considerably (32.7% in P1 and 26.7% in P2) regardless of the SD being the same. The difference in S enables the RGAM to reflect the influence of dead branch structures on IL, leading to an acceptable RGAM performance for both P1 and P2 (“fair” IL relative errors: −20.2% vs. −16.1%) in the present study of unmanaged coniferous plantations with high SD laden with dead branches.

1. Introduction

The forest canopy is the first contact surface for rainwater. When gross rainfall (GR) encounters forest canopies, it is partitioned into net rainfall (throughfall (TF): the portion falling to the forest floor through the canopy; stemflow (SF): the portion draining through the canopy along the branches and stem surfaces) and rainfall loss (interception loss (IL): the portion evaporating back to the atmosphere). The IL from trees is a primary component of evapotranspiration [1,2,3] and exerts a strong influence on the forest water resources and the environment in forest ecosystems [4,5,6,7,8,9]. The IL as a proportion of GR (IL/GR, %) is as much as 30% in mature coniferous plantations and is strongly influenced by forest stand structure parameters, such as stand density (SD), and canopy structure factors (plant area index, PAI; canopy cover, c) [2,4,6,8,10,11]. Therefore, understanding the effect of forest stand structure on IL is essential to regulating the water and environment in forest ecosystems.

Numerous studies of IL have been conducted in diverse forest types worldwide (e.g., [5,9,12,13]). Previous studies of IL conducted in different forest types determined the influence of forest stand structure factors, such as SD [4,11], c [6], PAI [14], and water storage capacity on the tree surface (S, [15]). These findings demonstrate that factors influential on IL depend on tree traits and climate conditions. Thus, studies of IL for the same species with different forest stand structures under the same climate condition help to give a better understanding of the changes in IL due to changes in forest stand structures, which would contribute to the assessment of local water and environmental management in forest ecosystems [2,4,11,16,17].

In Japan, forested areas account for 67% of the total land area, with plantations constituting 41% of the forested area [18]. The dominant plantation species are Japanese cedar (Cryptomeria japonica D.Don) and Japanese cypress (Chamaecyparis obtusa Endl.) with a typical planting density of ≥2500 stems ha⁻¹ [18]. Many plantations, however, have not been managed adequately because of the impacts of economic recession on forestry since the 1980s [19]. Unmanaged plantations, compared with managed plantations, generally feature high SD, high c, and sparse understory vegetation (Figure 1). In addition, unmanaged plantations with high SD have a complex stand structure, such as dead branches under canopies (see Figure 1b, [17,20,21]).

Nanko et al. [22] examined TF drop velocity in an unmanaged Japanese cypress plantation. The authors reported that some drips of TF had a lower velocity than expected probably due to the influences of dead branches around the stem as well as leaves on live branches on TF partitioning. Shinohara et al. [23] pruned some of the many dead branches under canopies in Japanese cypress plantations and observed that dead branches influence raindrop erosivity induced by TF. Jeong et al. [17] detected strong correlations between TF and dead branch structures (a negative relationship between TF and number of dead branches and a positive relationship between TF and vertical distance of dead branches; see Figure 6 of [17]). These findings imply that IL, which is the difference between GR and the sum of TF and SF, may differ between upper-canopy structure and under-canopy structure with dead branches. Thus, evaluation of the accuracy of IL estimation models is vital for their potential application to unmanaged coniferous plantations laden with dead branches.

Various IL estimation models have been developed [13]: empirical models (e.g., [24,25]) and physically based models, such as numerical models (e.g., [26,27]), analytical models (e.g., [28,29]), and stochastic models (e.g., [30]).

Komatsu et al. [11] targeted the dominant plantation species in Japan growing in similar climates for estimation of IL in response to forest management practices, such as thinning. The authors developed an empirical regression IL estimation model based on the relationship between IL/GR and SD using 46 data sets for two dominant plantation species in regions without heavy snowfall. The authors found that none of the other parameters (canopy height, species, and annual precipitation) improved the predictability of IL. Nagano et al. [31] included three additional data sets for sparse Japanese coniferous plantations (i.e., n = 49, SD range: 356 to 3864 stems ha⁻¹ with a mean of 1506 stems ha⁻¹). The equation is formulated as follows:

IL/GR = 0.290 {1 − exp (−0.00110 ∗ SD)} × 100

(1)

Komatsu et al. [11] reported that the estimation error of IL using Equation (1), i.e., the absolute value of IL_modeled/GR − IL_measured/GR, was as small as 4%. Komatsu et al. [11] did not explicitly state that unmanaged plantations laden with dead branches were included in the 49 datasets. Moreover, no previous studies have examined the applicability of Equation (1) to unmanaged coniferous plantations laden with dead branches.

The Gash analytical model was originally developed as a physically based IL estimation model in 1979 [28]. It was modified due to its failure to estimate IL in sparse forests [29], which is generally called the reformulated Gash analytical model (hereafter RGAM). Among physically based IL estimation models, RGAM is the most commonly used IL estimation model for various forest types worldwide (111 previous studies reviewed by Muzylo et al. [13]) because it requires relatively few parameters and variables that can be derived from observation data of forest stand structure, rainfall partitioning (RP), and meteorology [13]. The RGAM has been applied to unmanaged Japanese coniferous plantations before thinning [1] and after thinning [6], which showed that the RGAM is sufficiently robust and reliable for the estimation of IL for unmanaged Japanese coniferous plantations regardless of thinning practice. Although the presence of dead branches is known as the typical under-canopy characteristic in unmanaged Japanese coniferous plantations [20], the characteristics of dead branch structures (e.g., the number of dead branches [17,21]) has been poorly studied. Moreover, to the authors’ knowledge, no study has attempted to apply the RGAM to unmanaged Japanese coniferous plantations with high SD laden with dead branches almost fully distributed on the stems, although some previous studies may have implicitly tested the applicability of the RGAM in high-density coniferous plantations probably with dead branches.

This study aimed to examine the applicability of the RGAM to unmanaged coniferous plantations with high SD laden with dead branches. To do this, (1) we applied the RGAM to two Japanese coniferous plantation plots with the same high SD laden with dead branches but with different numbers of dead branches; (2) we performed a sensitivity analysis to detect essential parameters for RGAM IL estimation; and (3) we assessed the applicability of the RGAM together with Equation (1) to such coniferous plantations.

2. Materials and Methods

2.1. Study Site

The study site is located at the Takada experimental site of Kasuya Research Forest, Kyushu University, Fukuoka, Japan (33°37′58″ N, 130°31′48″ E, 100 m a.s.l.; Figure 2a). The study site is covered with Japanese cypress (Chamaecyparis obtusa Endl.) trees planted in 1985. The plantation has not been managed. The mean annual temperature and mean annual precipitation are 17.2 ± 0.5 (S.D.) °C and 1649.0 ± 337.6 (S.D.) mm, respectively, based on the Fukuoka meteorological observatory’s data for 1986−2017 measured at the nearest long-term meteorological station. The rainy and snowy seasons are from June to September and December to February, respectively [23].

Two 20 m × 10 m plots (hereafter P1 and P2, respectively) were established at the study site. P1 and P2 each contained 50 trees (i.e., 2500 stems ha⁻¹), among which five standing trees were dead in each plot (i.e., in the process of self-thinning). Jeong and Otsuki [21] reported that the stand-scale SF/GR can be overestimated when standing-dead trees are considered alive and included for the calculation of the stand-scale SF/GR. The authors suggested that the stand-scale SF/GR would be reliably estimated with only standing-live trees, and the standing-dead trees should not be included in standing trees for stand-scale SF/GR estimation. Thus, the standing-dead trees were excluded from the SD calculation (i.e., 45 trees in each plot, equal to 2250 stems ha⁻¹) to avoid errors in IL estimation through SF [21]. The plot slopes were relatively steep (26°). The forest structure was categorized into (1) stand structure, (2) upper-canopy structure, and (3) under-canopy structure (i.e., dead branch structure) (Figure 2b and

Table 1. Forest stand structure categorized into general stand structure, stand structure, upper-canopy structure, and under-canopy structure of the unmanaged Japanese cypress plantations laden with dead branches under the same and high stand density. P1 and P2 denote each study plot. The data in this table were retrieved from [17].

Forest Stand Structure	P1	P2
General stand structure
a The number of trees (Ntree)	50 (45 + 5)	50 (45 + 5)
Tree age (years)	33	33
b Stand density (SD, stems ha−1)	2250 (2500 − 250)	2250 (2500 − 250)
Stand structure
c Diameter at breast height (DBH, cm)	17.7 ± 4.1	20.4 ± 4.1
d Tree height (h, m)	12.7 ± 2.0	15.4 ± 2.0
c Basal area (BA, m2 ha−1)	64.9	85.2
e Plant area index (PAI, m2 m−2)	5.6 ± 0.2	6.4 ± 0.2
Upper-canopy structure
f The number of live branches (Nlb)	15.7 (21.8%)	14.7 (23.6%)
d Canopy length (CL, m)	2.8 ± 1.4	3.0 ± 1.1
g Canopy projection area (CPA, m2)	4.8 ± 1.8	4.7 ± 1.9
e Canopy cover (c, %)	99.3 ± 0.1	99.6 ± 0.1
Under-canopy structure
f The number of dead branches (Ndb)	56.3 (78.2%)	47.4 (76.4%)
h Vertical distance of dead branches (Vd, cm)	14.1 ± 1.8	22.9 ± 3.8

Data are given as average ± standard deviation. ^a The numbers in parenthesis denote the number of standing-live trees + the number of standing-dead trees. ^b SD was calculated with only standing-live trees. ^c DBH was measured for all trees in P1 and P2 using diameter calipers; BA was calculated from the DBH. ^d Measured for all trees in P1 and P2 using a hypsometer (Haglȍf Vertex IV-360, Haglȍf Inc., Madison, MS, USA) by dividing the branch distribution (height of treetop, lowest live branch and lowest dead branch). ^e Measured at 30 points in P1 and P2, respectively, using a pair of LAI-2000 plant canopy analyzers (Li-Cor Inc., Lincoln, NE, USA). ^f Counted using a telescope. ^g Measured for all trees in P1 and P2 at the horizontal distance from the tree to the projected edge of the crown along six compass directions (N, NE, SE, S, SW, NW) using a tape measure. ^h calculated as 1/(the height of dead branch length (m)/N_db) × 100. Note that PAI is measured for both leaves and branches at the same time. Accordingly, it was included in the stand structure to divide into the canopy and under-canopy structures.

). Stand structures and under-canopy structures were statistically different, whereas no statistical difference in upper-canopy structures was observed (Table 1). Specifically, stand structure parameters (diameter at breast height (DBH), h, basal area (BA), and PAI) were larger in P2 than those in P1 (t-test, p < 0.001). For under-canopy structures (number of dead branches (N_db) and vertical distance of dead branches (V_d)), N_db in P1 was larger than that in P2, and V_d in P1 was shorter than that in P2 (t-test, p ≤ 0.001). Upper-canopy structure parameters (number of live branches (N_lb), crown length (CL), and crown projection area (CPA)) were almost identical between P1 and P2 (t-test, p > 0.05).

2.2. Measurement of Throughfall and Stemflow

TF was monitored using 30 storage-type TF rain collectors (210-mm-diameter funnel) in P1 and 28 counterparts in P2. The TF rain collectors were randomly arranged and maintained throughout the study period. Note that the relative standard errors of the mean stand-scale TF between the two plots, assessed using the following Student’s t-value analysis at the significance level of 0.05, were almost identical (5.0% vs. 5.4%, [17,32]).

$$\epsilon =\frac{\left(t_{0.05,n-1}\times CV\right)}{\sqrt{n}}$$

(2)

where ε is the error of the mean, t_0.05,n−1 is the Student’s t-value, CV is the coefficient of variation, and n is the number of TF rain collectors.

Nine standing live trees were selected to measure SF according to the DBH classes in each plot. SF was collected using a hose, plastic seat, and silicon sealant at the height of 1.5 m, which drained into two 90 L connected bins to prevent overflow. The stand-scale SF (mm) was estimated with the following equation [16,17,21,33]:

$$SF=\frac{{{\displaystyle \sum }}^{\text{}}SF{V}_i}{n_{tree}}\times \frac{N_{all}}{PA}$$

(3)

where SFV_i is the sum of the SF volume (L) of the observed trees, n_tree and N_all are the number of observed trees (n = 9) and all standing-live trees (n = 45) in each plot, respectively, and PA is the plot area (200 m²).

TF and SF were collected on a weekly basis but on a day with no rain during the growing season from April to October 2017. The range of the periods was 5 to 9 days, with an average of 7 days, and thus we call the period “weekly” hereafter. We obtained 25 weekly RP data and used them for the analysis. IL was calculated as the difference between GR and the sum of TF and SF as follows:

$$IL=GR-\left(TF+SF\right)$$

(4)

More detailed descriptions of the study site and measurement of RP are reported in our previous study [17].

The growing season from April to October is coincident with the rainy season in this region. The precipitation at the Fukuoka meteorological observatory from April to October during 1986−2017 was 621.5−1879.0 mm with a mean of 1244.0 ± 300.1 (S.D.) mm, accounting for 75.4% of the mean annual precipitation. During the study period from April to October 2017, total GR was 1068.4 mm, which accounted for 85.9% of the mean precipitation for the same period, which implies that the measurement period was under the normal meteorological condition of this region. Two previous studies reported no relation between annual GR and IL/GR for Japanese coniferous plantations in Japan [11,16] and also support this assumption that our measurements for IL estimation could be representative of general meteorological conditions in this region.

2.3. Measurement of Meteorological Data

Meteorological conditions were measured on a 2 m tower installed in an opening 50 m from P1 and P2. Before measuring meteorological conditions in the opening, we removed all obstructions such as trees and shrubs within 45° around the 2 m tower for similar meteorological conditions to the above-canopy conditions. We continuously managed the opening area during the study period. GR was measured using a tipping bucket with a resolution of 0.2 mm (HOBO RG3-M rain gauge, Onset Computer, Bourne, MA, USA) with a data logger (HOBO Event, Onset Computer, Bourne, MA, USA). There were missing data for 1 week from the 0.2 mm rain gauge (GR_0.2), which were substituted with data from a tipping bucket with a resolution of 0.5 mm (GR_0.5) measured in the opening (Ota Keiki Seisakusho Co. Ltd., Tokyo, Japan) (GR_0.5 = 1.00 × GR_0.2; R² > 0.97). Air temperature (T, °C) and relative humidity (RH, %) were measured with a temperature/humidity probe (S-THB-M002, Onset Computer, Bourne, MA, USA). The data were recorded at 10 min intervals with a data logger (HOBO U30 Weather Station; Onset Computer, Bourne, MA, USA). Solar radiation (R_s) was measured in an opening site approximately 1 km away from the study site using a solar radiometer (ML-020VM, EKO Instruments, Tokyo, Japan). Wind speed (WS) was measured at Fukuoka airport, located approximately 7 km southwest of the study site, by the Japan Meteorological Agency.

2.4. Application of RGAM

In the RGAM, GR is deemed to occur as a series of discrete events and is separately analyzed using the threshold gross rainfall (GR′) to determine whether the canopy is saturated or not. Based on the rainfall and canopy structure, GR′ was determined as follows:

$$G{R}^{\prime }=-\frac{\overline{R}}{\overline{E_{\mathrm{c}}}}{s}_{\mathrm{c}}ln\left(1-\frac{\overline{E_{\mathrm{c}}}}{\overline{R}}\right)$$

(5)

$$S_{\mathrm{c}}=\frac{S}{c}$$

(6)

where $\overline{R}$ is mean GR rate during the saturated canopy conditions (mm h⁻¹), $\overline{E}$ is mean wet canopy evaporation rate under saturated canopy conditions per unit ground area (mm h⁻¹), $\overline{E_{\mathrm{c}}}$ is mean wet canopy evaporation rate under saturated canopy conditions from the canopy covered area (mm h⁻¹), S is canopy storage capacity (mm), and S_c is canopy storage capacity per canopy cover (mm). Based on the assumption of Gash [28], $\overline{R}$ and $\overline{E_c}$ were calculated for all the hourly periods with GR greater than 0.5 mm h⁻¹, which could represent the average rainfall and evaporation rates under saturated canopy conditions, respectively.

For a given rainfall event, the components of IL were calculated as follows:

$$IL=\{ \begin{array}{c}I{L}_{insuf},\hspace{1em}\hspace{1em}GR<G{R}^{\prime }\\ I{L}_{suf},\hspace{1em}\hspace{1em}GR\ge G{R}^{\prime }\end{array}$$

(7)

where IL_insuf is the interception (mm) during a small rainfall event that was insufficient to saturate the canopy and IL_suf is the interception (mm) during a large rainfall event that was sufficient to saturate the canopy.

IL_insuf was calculated as follows:

$$I{L}_{insuf}=c{\displaystyle \sum }_{j=1}^{m}G{R}_j$$

(8)

where m is the number of storms that are insufficient to saturate the canopy.

IL_suf consists of four IL components and was calculated as follows:

$$I{L}_{suf}=I{L}_w+I{L}_s+I{L}_d+I{L}_t$$

(9)

$$I{L}_w=nc\left(G{R}^{\prime }-S_c\right)$$

(10)

$$I{L}_s=c\frac{\overline{E_c}}{\overline{R}}{\displaystyle \sum }_{j=1}^n(G{R}_j-G{R}^{\prime })$$

(11)

$$I{L}_d=nc{S}_c$$

(12)

$$I{L}_t=q{S}_t+p_t{\displaystyle \sum }_{j=1}^{n-q}G{R}_j\hspace{1em}(G{R}_j>P_t{}^{\prime })$$

(13)

$$P_t{}^{\prime }=S_t/p_t$$

(14)

where IL_w is IL for wetting the canopy (mm), n is the number of storms sufficient to saturate the canopy, IL_s is IL for wet canopy evaporation during storms (mm), IL_d is IL for wet canopy evaporation after storms cease (mm), IL_t is IL for evaporation from trunks (mm) for q storms that saturate the trunks, q is the number of storms sufficient to saturate the trunks, S_t is the trunk storage capacity (mm), p_t is the proportion of rainfall diverted to SF, n−q is the number of storms insufficient to saturate the trunks, and P_t′ is the threshold gross rainfall sufficient to saturate trunks (mm).

For calculation of the wet canopy evaporation rate under saturation canopy condition E_c (mm per unit of covered area, and assumed to be equal for the whole plot area), the Penman–Monteith equation setting canopy resistance (r_s) to zero was used as follows:

$$E_c=\left(\Delta R_n+\frac{\rho c_{p}D}{r_a}\right)/\lambda \left(\mathsf{\Delta }+\gamma \right)$$

(15)

where Δ is the rate of change of saturated vapor pressure with temperature (Pa K⁻¹), R_n is the net radiation on the canopy (W m⁻²), ρ is the density of dry air (kg m⁻³), c_p is the specific heat of air at constant temperature (J kg⁻¹ K⁻¹), D is the vapor pressure deficit (Pa), r_a is the aerodynamic resistance (s m⁻¹), λ is the latent heat of vaporization of water (J kg⁻¹), and γ is the psychometric constant (Pa K⁻¹). Given that we did not measure net radiation, R_n was assumed to be 0.8 times R_s [4,8], where R_s is the incident solar radiation (W m⁻²).

The RGAM can be run both on an event basis (e.g., [6,8]) and on a daily basis assuming one storm per rain day [29]. In this study, the RGAM was run on an event basis. The drying time for canopy wetting in Japanese coniferous plantations was assumed to be 6−8 h in the humid temperate climate of Japan (e.g., [6,8,16]). Consequently, a rainfall event was classified as a separate event when the period without rainfall after its cessation exceeded 8 h. Among 25 weekly data sets, 6 data sets included a single rainfall event in each event, and 19 data sets included 2–5 rainfall events. The event-based IL estimates derived from the RGAM were added up for the period of weekly RP measurement to obtain the estimated weekly IL. We then compare it with the measured weekly IL.

2.5. Parameter Estimation for RGAM

Weekly RGAM parameters (S, S_t, and p_t) were derived from our weekly RP measurement as in [8]. S is a crucial parameter for determining IL [14,34,35,36,37,38,39]. We calculated the values of S as the negative intercept of the linear regression of the sum of TF and SF against GR [38], instead of using the minimum method [17,40], because the latter may provide unrealistically low values [38]. We calculated the values of S_t and p_t as the negative intercept and slope of the linear regression of SF against GR [28].

Hourly RGAM parameter r_a was estimated from the hourly wind speed and stand structure, employing the momentum method as follows:

$$r_a={\left(k^{2}u\right)}^{-1}\ln {\left[\frac{z-d}{z_o}\right]}^2$$

(16)

$$z=h+2$$

(17)

$$z_o=\left\{0.264 \left(1.0-\frac{d}{h}\right)\right\}\times h$$

(18)

$$d=\left(1.0-\frac{1.0-\exp \left(-\alpha SD\right)}{\alpha SD} \frac{1.0-\exp \left(-\beta PAI\right)}{\beta PAI}\right)\times h$$

(19)

where k is the von Kármán constant (0.41), u is the wind speed (m s⁻¹), z is the reference height (m), h is tree height (m). z_o is the roughness length (m), d is the zero-plane displacement (m), and α and β are empirical coefficients (α = 0.000714 and β = 0.273), respectively. The value of z was taken as 2 m higher than tree height. The values of z_o and d were calculated using SD, PAI, and h in accordance with [41].

The above-mentioned parameters derived in this study were compared with those from previous studies of Japanese coniferous plantations with different forest stand structures.

2.6. Model Evaluation

To evaluate the RGAM, the root mean square error (RMSE) was employed as follows:

$$\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}^{\text{}}{\left(P_{\mathrm{j}}-O_{\mathrm{j}}\right)}^2}{w}}$$

(20)

where P_j is the estimated IL (mm), O_j is the measured IL (mm), and w is the number of testing data (dimensionless). The RMSE value of zero indicates a perfect model.

In addition, we classified RGAM accuracy using the cumulative IL relative error consistent with the criteria of [13]:

• extremely good (error ≤ 1%)
• very good (1% < error ≤ 5%)
• good (5% < error ≤ 10%)
• fair (10% < error ≤ 30%)
• bad (error > 30%).

Acceptable relative errors for IL estimation can be as high as 20% [13]. This classification for IL estimation model evaluation has been widely employed (e.g., [6,35,36]). In this study, model parameter calibration and validation were performed in each stand (P1 and P2) using all the data derived from each stand to independently test the applicability of the RGAM in the two stands with different under-canopy structures under the same SD.

2.7. Sensitivity Analysis

To detect the essential parameters in the RGAM for IL estimation, we performed sensitivity analysis by specifically categorizing parameters into:

(1)

forest structure parameters

stand-related: z_o, canopy-related: c, trunk-related: p_t

canopy storage-related: S, trunk storage-related: S_t

evaporation-related: $\overline{E_c}$

(2)

climate parameter: $\overline{R}$.

We increased or decreased these parameter values by as much as 40% of their initial value except for c, which we decreased by as much as 40% and did not increase owing to its almost complete canopy coverage of ≈1 (0.993 for P1, 0.996 for P2; Table 1).

3. Results

3.1. Interception Loss

During the study period (total GR: 1068.4 mm), RP of P1 and P2 differed considerably under the same SD of 2250 stems ha⁻¹, resulting in a difference in IL/GR of 6.0% (Figure 3). For P1, TF and SF were 493.0 mm (46.1%) and 226.0 mm (21.2%), respectively, resulting in IL of 349.4 mm (32.7%). For P2, TF and SF were 570.4 mm (53.4%) and 212.5 mm (19.9%), respectively, resulting in IL of 285.5 mm (26.7%).

Weekly IL/GR of P1 ranged from 21.2% to 94.6% with a mean of 43.9% and P2 ranged from 17.2% to 94.7% with a mean of 38.9% (Figure 4). The value logarithmically decreased with the increase in weekly GR amount (mm) for P1 (R² = 0.75) and P2 (R² = 0.82) (Figure 4a). The value also logarithmically decreased with weekly mean GR intensity (mm h⁻¹) for P1 (R² = 0.63) and P2 (R² = 0.60) (Figure 4b). However, a weak relationship was observed between weekly IL/GR and weekly GR duration (h) for P1 and P2 (R² < 0.05) (Figure 4c). The difference in weekly IL/GR between P1 and P2 was greater with increase in weekly GR amount and weekly mean GR intensity (Figure 4a,b).

3.2. Parameters for RGAM

In the comparison of the model parameters with those of previous studies for sensitivity analysis of RGAM IL estimation [6,8], z_o (stand-related parameter), p_t (trunk-related parameter) and S (canopy storage-related parameter) in the present study differed from those of the previous studies, whereas the other parameters were compatible (

Table 2. Interception parameters (forest structure parameters: stand-related, canopy-related, trunk-related, canopy storage-related, trunk storage-related, and evaporation-related parameters; and climate parameter) in the reformulated Gash analytical model (RGAM) for the present two study plots (P1 and P2), including the previous studies of Japanese coniferous plantations. SD and PAI denote stand density and plant area index, respectively.

				Forest Structure Parameters						Climate Parameter
Study	Species	SD	PAI	Stand-Related	Canopy-Related	Trunk-Related	Canopy Storage-Related	Trunk Storage-Related	Evaporation-Related
				zoa	cb	ptc	Sd	Ste	$\overline{E_c}$ f	$\overline{R}$ g
		Stems ha−1	m2 m−2	m			mm	mm	mm h−1	mm h−1
Present study
This study (P1)	Cypress	2250	5.6	0.85	0.99	0.24	3.94 1	1.19	0.23	2.84
This study (P2)	Cypress	2250	6.4	0.95	1.00	0.21	3.25 1	0.58	0.24	2.84
Previous studies
[6] *	Cypress	2198	−	1.60	0.97	0.12	2.42 2	0.39	0.38	2.37
[6] **	Cypress	1099	−	1.60	0.76	0.09	1.54 2	0.34	0.32	3.46
[8] *	Cedar	1300	4.8	4.41	0.98	0.16	1.34 1	0.61	0.28	2.75
[8] **	Cedar	600	2.4	1.63	0.86	0.07	0.83 1	0.42	0.06	1.44

Asterisks * and ** denote pre-thinning and post-thinning, respectively. Forest structure parameters: ^a Stand-related parameter: z_o, roughness length (m). ^b Canopy-related parameter: c, canopy cover. ^c Trunk-related parameter: p_t, the proportion of rain diverted into stemflow. ^d Canopy storage-related parameter: S, canopy storage capacity (mm). ^e Trunk storage-related parameter: S_t, trunk storage capacity (mm). ^f Evaporation-related parameter: $\overline{E_c}$ (mm h⁻¹): mean evaporation rate during rainfall under the saturated canopy conditions per unit of canopy-covered area. Climate parameter: ^g $\overline{R}$: mean gross rainfall rate (mm h⁻¹) during the saturated canopy conditions. Superscripts 1 and 2 denote the mean method [39] and the individual storm method [37], respectively.

). The value of z_o was smaller because tree height was lower, p_t was larger because SF/GR was larger, and S was larger because SD was larger compared with those of the previous studies.

In the comparison of the parameters between the two study plots, S and S_t (trunk storage-related parameter) distinctly differed, whereas the other parameters were almost identical. The parameters of P1 were larger than those of P2 (Table 2).

3.3. Performance of RGAM

The predominant component of IL estimated using the RGAM for both P1 and P2 was wet canopy evaporation after storms ceased (I_d) (48.1% for P1 and 47.4% for P2), followed by wet canopy evaporation during the storms (I_s) (26.2% for P1 and 32.7% for P2) (

Table 3. Components of estimated interception loss (IL) using the reformulated Gash analytical model (RGAM) during the growing season in the two study plots of the unmanaged Japanese cypress plantations laden with dead branches. P1 and P2 denote the study plot 1 and plot 2.

Components	P1			P2
	Modeled IL (mm)	Proportion of IL (%)	Estimated Error (%)	Modeled IL (mm)	Proportion of IL (%)	Estimated Error (%)
For storms GR < GR′ (ILinsuf)
Evaporation from canopy	26.2	9.4		22.3	9.3
For storms GR ≥ GR′ (ILsuf)
Wetting up the canopy (ILw)	5.8	2.1		5.1	2.1
Wet canopy evaporation during the storms (ILs)	73.1	26.2		78.3	32.7
Evaporation after storms ceased (ILd)	134.0	48.1		113.6	47.4
Evaporation from trunks (ILt)	39.7	14.2		20.3	8.5
Total estimated IL	278.7	100.0	−20.2	239.6	100.0	−16.1
Total measured IL	349.4	−		285.5	−

). These components accounted for 75.1% and 80.7% of the estimated IL for P1 and P2, respectively.

The RGAM underestimated cumulative IL for both P1 and P2 (Figure 5), in which the accuracy of the RGAM for P1 was lower than that for P2. For P1, the cumulative measured and estimated IL was 349.4 mm (IL/GR = 32.7%) and 278.7 mm (IL/GR = 26.1%), respectively, which resulted in a 6.6% underestimate with a cumulative IL relative error of 20.2%. For P2, the corresponding values were 285.5 mm (IL/GR = 26.7%) and 239.6 mm (IL/GR = 22.4%), respectively, which resulted in a 4.3% underestimate with a cumulative IL relative error of 16.1%. Although the estimation error of P1 was larger than that of P2, the performance of the RGAMs for P1 and P2 based on their cumulative IL relative errors were both classified as “fair” (10% < error ≤ 30%).

Lower reproducibility in P1 than P2 was observed on a weekly basis in RGAM IL estimation (Figure 6). The coefficient of determination (R²) derived from the relationship between the measured and estimated IL was 0.737 for P1 and 0.833 for P2, and the RMSE was 5.72 mm for P1 and 3.73 mm for P2.

3.4. Sensitivity Analysis

Sensitivity analysis confirmed that $\overline{R}$ (a climate parameter) showed nonlinear negative sensitivity for RGAM IL estimation, as reported previously (Table 2 [1]). In contrast, all forest structure parameters produced almost linear positive sensitivities, and the changes for both P1 and P2 were similar (Table 2 and Figure 7). Among the model parameters, z_o (stand-related parameter), S (canopy storage-related parameter), and $\overline{E_c}$ (evaporation-related parameter) were sensitive to RGAM IL estimation.

Among the above-mentioned three dominant parameters for IL estimation, S showed the greatest difference in sensitivity between the two study plots: when the value of S was reduced by 20%, the difference in the sensitivity of IL is 4.0% (Figure 7). In contrast, the differences in sensitivities of the parameters z_o and $\overline{E_c}$ were negligible between the two study plots: when the values of z_o and $\overline{E_c}$ were reduced by 20%, the differences in sensitivity were 1.3% and 1.3%, respectively.

4. Discussion

4.1. Applicability of RGAM IL Estimation

This study examined the applicability of the RGAM, the most commonly used IL estimation model across various forest types worldwide (as revealed in a literature review of 111 studies [13]), to unmanaged Japanese coniferous plantations with high SD laden with dead branches almost fully distributed on the stems. In the current study, “fair” performances of the RGAM were confirmed with IL relative errors of −20.2% for P1 and −16.1% for P2.

Previously, several studies applied the RGAM to Japanese coniferous plantations. Shinohara et al. [8] reported that the RGAM performed well for IL estimation in Japanese cedar plantations before thinning but may not be applicable after high-intensity selective thinning (reduction of SD by 54% and BA by 50%). However, Sun et al. [6] successfully estimated IL using the RGAM for unmanaged Japanese cypress plantations before and after high-intensity strip thinning (reduction of SD by 50% and BA by 48%) and reported “good” performances with IL relative errors of −5.7% (before thinning) and 6.8% (after thinning). Although neither the presence nor features of dead branches are described in [6], the performance could be partly underestimated by the dead branches influential on IL, which is the typical under-canopy characteristic of unmanaged Japanese cypress plantations with high SD [17,20,21].

To compare the applicability of the RGAM with the empirical model, Equation (1) was applied to P1 and P2. Komatsu et al. [11] reported that SD was the predominant parameter in IL estimation for Japanese coniferous plantations and introduced Equation (1) with a typical IL estimation error of 4.0%. Based on Equation (1), we estimated the IL/GR in P1 and P2 under the same SD of 2250 stems ha⁻¹ in both plots to be 26.6%. The IL/GR estimate for P1 was 5.3% underestimated with a cumulative IL relative error of 18.8%, which was classified as “fair” performance. By contrast, the IL/GR estimate for P2 was only approximately 0.2% underestimated with a cumulative IL relative error of 0.6% and thus was classified as “extremely good” performance.

Although the accuracy of RGAM IL estimation in P1 and P2 was lower than that reported by Sun et al. [6], and almost the same in P2 and lower in P1 than those estimated by Equation (1), the IL estimation errors in this study could be considered acceptable because the errors lay within the “fair” performance class (10% < error ≤ 30%) [13]. The empirical model of Equation (1) showed “fair” and “extremely good” performances for IL estimation for Japanese coniferous plantations with high SD laden with dead branches, and implied high practicality with a single parameter of SD for the assessment of IL in response to forest management practices such as thinning [2,11]. This high practicality might not be surprising because of the dataset of Komatsu et al. [11]. The dataset includes several data of unmanaged Japanese cypress plantations laden with dead branches (data number 42–45 of Table 1 in [11]), which is a typical characteristic of unmanaged Japanese cypress plantations [17,20,21]. However, the relatively larger error in P1 than that in P2 using Equation (1) could occur for instances of different dead-branch structures (e.g., the number and vertical distance of dead branches [17]. The physically based RGAM showed similar accuracy as Equation (1) and could be applied for the assessment of IL in response to forest management practices affecting not only SD but also other forest stand structure parameters.

However, relatively larger errors compared with that of Sun et al. [6] were observed for RGAM IL estimation on a weekly basis in P1 and P2 where the dead branches were densely distributed (Figure 6). Measurement error can be derived from the non-above canopy meteorological data [42] and the accuracy of a relative humidity (RH) sensor [43,44]. In this study, our meteorological conditions were measured in an opening, the RH accuracy is within ±5% for RH ≥ 90%, and the wind speed was measured 7 km away from our study site. This could have partly caused our relatively larger errors. Therefore, further research is required (1) to parameterize the dead branch characteristic and (2) to measure RP on an event basis and meteorological conditions above canopies with high accurate instruments for more accurate parameters (e.g., vapor pressure deficit and aerodynamic resistance) sensitive to the Penman-Monteith equation and RGAM [45,46]. This could enable us to develop a model for a more precise estimation of IL. In this regard, we recommend TF partitioning studies (free TF, splash TF, and drip TF) measured using a disdrometer [47,48,49], and precise above-canopy measurements with calibrated instruments (e.g., rain gauges, [50]; humidity sensors, [43,44]).

4.2. Characteristics of IL Processes in the Unmanaged Plantations

The contribution of I_d (IL after rainfall cessation) to the total estimated IL was predominant in both P1 and P2, accounting for approximately half of the total IL estimates (48.1% for P1 and 47.4% for P2; Table 3). This indicates that the trees store a large volume of rainwater and then evaporation of the stored rainwater is favored. This result was inconsistent with a previous study of an unmanaged Japanese cypress plantation that reported the main component of the estimated IL to be I_s (I_s = 62.9% [1]). Note that we assume that rainfall characteristics, such as storm size, duration, and frequency, were no different between this study and previous studies conducted in humid temperate regions.

The large S observed in the present study (Table 2) could be the reason for the predominant I_d because S is the primary controller of I_d in the RGAM [15,34,35]. As expected, S was the most influential sensitive parameter for IL estimation in the current study (Figure 7). Previous studies similarly reported that S is among the most sensitive parameters in the RGAM (e.g., [1,8,29,35,51]).

Compared with previous studies of Japanese coniferous plantations, the values of S observed in the present study were considerably higher (Table 2) and within the upper range of S compared with other coniferous forests (1.2–4.3 mm in Douglas fir, 0.5–1.0 mm in Pinus spp., and 1.2–2.0 mm in Picea sitchensis [37]). The substantially higher S values in the present study were possibly due to the densely distributed dead branches under canopies accounting for a considerable proportion of S [52]. The regression method used in this study has been widely used for an S estimate in diverse forest types worldwide (e.g., [8,35,38,52]). However, it might not be optimal because of the following reasons. First, it often involves a subjective selection of rain events and dry interval durations, leading to different model parameter value sets [53]. Second, the regression method does not separate S from S_t, resulting in higher values of S than expected if SF is not negligible [54]. Thus, further S derivation approaches for a more accurate S estimate such as automatic model parameterization procedures with optimization of S (e.g., [53]) and a new procedure to estimate S, S_t, and p_t separately for any rain event [54] would reduce the estimated error using RGAM and improve our understanding of the impact of under-canopy structure with dead branches on S in unmanaged Japanese coniferous plantations.

Between the two study plots, IL in P1 was larger than that in P2 (Figure 3). The upper-canopy structures in P1 and P2 were statistically identical (see Section 2.1) and thus were unlikely to account for the difference in IL between the two plots. Although stand structure parameters (i.e., DBH, h, BA, and PAI) differed between the two study plots (see Section 2.1), z_o (a stand-related parameter) and $\overline{E_c}$ (an evaporation-related parameter), which are affected by stand structure through aerodynamic controls, did not differ between P1 and P2. Although these parameters were sensitive for RGAM IL estimation, the difference in sensitivity of these parameters between P1 and P2 was negligible. A possible reason for this result is that the difference in stand structure between P1 and P2 was insufficient for differentiation of the aerodynamic controls on IL.

The under-canopy structures of P1 and P2 were statistically different (see Section 2.1). The values of S (a canopy storage-related parameter) affected by both upper-canopy and under-canopy structures were larger in P1 than in P2 (Table 2). Although S was the dominant parameter sensitive for RGAM IL estimation together with z_o and $\overline{E_c}$, S alone showed a difference in sensitivity for RGAM IL estimation between P1 and P2: when the value of S was reduced by 20%, the difference in sensitivity is 4.0% (Figure 7). The higher S in P1 than in P2 was probably due to the difference in number (56 and 47 branches per tree, respectively) and narrower vertical space (14.1 and 22.9 cm, respectively) of dead branches (Table 1 and Figure 2b). S was estimated separately from RGAM in this study. It was obtained from the measured RP data influenced by forest stand structures during rainstorms (i.e., stand structure, upper-canopy structure, and under-canopy structure with dead branches) [17,52]. Since the RGAM has no direct parameter quantifying dead branch structures, another approach such as adding new parameters of dead branch structures to the RGAM formula is needed to improve the RGAM performance for forests laden with dead branches.

Microtopography could result in differences in IL between flat areas and mountainous areas [8]. For example, Saito et al. [33] estimated IL/GR at 26% under SD of 1107 stems ha⁻¹ in a mountainous area in a mature Japanese cedar plantation. Shinohara et al. [8] estimated IL/GR under SD of 1300 stems ha⁻¹ in a plain area in a mature Japanese cedar plantation. The authors reported approximately half of the IL/GR of Saito et al. [33] although the two study sites were 15 km away. These results imply that IL could be different in unmanaged Japanese coniferous plantations laden with dead branches under different microtopography. Future studies examining differences in IL processes and modeling are required for evaluating the potential impacts of microtopography in unmanaged Japanese coniferous plantations laden with dead branches on forest water resources and the environment.

5. Conclusions

This study fills a knowledge gap in quantifying and estimating IL using the RGAM for Japanese coniferous plantations with high SD laden with dead branches. Although SD is a predictable parameter for IL estimation in response to forest management practices, such as thinning, for Japanese coniferous plantations, the total IL/GR between P1 and P2 differed by 6.0% (32.7% for P1 vs. 26.7% for P2) despite the same SD in the two plots. The S values in P1 and P2 were extremely high compared with those of previous studies but differed between the plots (3.94 and 3.25 mm, respectively), which were influenced by the different number of dead branches (56 and 47 branches per tree, respectively). The difference in S enables the RGAM to reflect the influence of dead branch structure on IL, leading to an acceptable RGAM performance for both P1 and P2 (“fair” IL relative errors: −20.2% vs. −16.1%) in the present study of unmanaged coniferous plantations with high SD laden with dead branches.