Development of a Model to Evaluate Water Conservation Function for Various Tree Species

: The ecosystem services of forests, such as the water conservation function, are the combined results of diverse processes, and the modification of one part of a forest affects each ecosystem service separately via complex processes. It is necessary to develop an ecosystem service assessment model for various tree species to ensure proper forest management. In this study, a model to evaluate three ecosystem services, namely, the water supply, wood supply, and carbon sink, for various tree species in Japan is developed using many observation data from the previous literature. The integrated evaluation model consists of the forest model, hydrological model, and carbon stock assessment model. The forest model consists of the forest growth model and LAI estimation model, based on allometry. The results of the simulations for the major tree species yield the following findings: (1) Water supply varies with tree species but decreases until about 40 years of age, after which it is near constant. (2) Although beech has a larger LAI than needleleaf forests, water supply is not significantly different. (3) Broadleaf forests are more affected by thinning than needleleaf forests and tend to receive increased water supply as a result of processes such as thinning. This study enabled the evaluation of water conservation function in watersheds containing various tree species.


Introduction
Forests provide a variety of ecosystem services such as disaster prevention, water supply, water purification, global warming mitigation, timber production, and biodiversity.In order to obtain sustainable ecosystem services from forests, it is necessary to apply effective forest management.Ecosystem services are the combined results of diverse processes, and the modification of one part of a forest affects each ecosystem service separately via complex processes.Thus, there are trade-offs whereby some practices increase the effectiveness of timber supply but reduce water supply, or synergies whereby both services are improved [1].In Japan, there is a tax on forests called the Forest Environment Tax [2], and it is necessary to explain how the ecosystem services obtained will change if forests are improved.Therefore, a model is needed to evaluate how the ecosystem services provided vary with tree species and forest condition.Forest management in Japan has so far focused on timber production.However, a survey conducted by the Cabinet Office in Japan showed that disaster prevention, water resources, and global warming mitigation were the top three functions expected of forests by the public, indicating that the public considers the water Water 2024, 16, 588 3 of 23

Forest Model
The forest model is divided into two parts: the forest growth model, which calculates forest annual changes using growth curves and inverse yield-density effects, and the forest structure model, which calculates forest factors such as the diameter at breast height and LAI (leaf area index) using other forest structural factors such as stand density and tree height.Details of the model parameters are given in Appendix A.
where t is forest age (y), and M, L, and K are parameters determined for each tree species and each area.Variation in growth at the same forest age is adjusted for the growth status.
The growth status is a parameter determined for each location based on sunlight, soil, moisture condition, and other factors.In this study, the tree height determined through Equations ( 1)-( 3) was taken as status 3, and the index of statuses 1-5 given to each area was used to represent the variation in growth.The distribution range of statuses was determined so that 95% of the forest sample said to determine the parameters of the growth curves fell between statuses 1 and 5. Figure 1 shows the plots and the growth curves of Japanese cedar (Cryptomeria japonica) and type L-1 broadleaf forests, such as beech (Fagus crenata), Mizunara oak (Quercus crispula), Keyaki (Zelkova serrata), walnut (Juglans), horse chestnut (Aesculus turbinata), and so on, for samples of statuses 1-5.

Forest Model
The forest model is divided into two parts: the forest growth model, which calculates forest annual changes using growth curves and inverse yield-density effects, and the forest structure model, which calculates forest factors such as the diameter at breast height and LAI (leaf area index) using other forest structural factors such as stand density and tree height.Details of the model parameters are given in Appendices A and B.
where  is forest age (y), and , , and  are parameters determined for each tree species and each area.Variation in growth at the same forest age is adjusted for the growth status.The growth status is a parameter determined for each location based on sunlight soil, moisture condition, and other factors.In this study, the tree height determined through Equations ( 1)-( 3) was taken as status 3, and the index of statuses 1-5 given to each area was used to represent the variation in growth.The distribution range of statuses was determined so that 95% of the forest sample said to determine the parameters of the growth curves fell between statuses 1 and 5. Figure 1 shows the plots and the growth curves of Japanese cedar (Cryptomeria japonica) and type L-1 broadleaf forests, such as beech (Fagus crenata), Mizunara oak (Quercus crispula), Keyaki (Zelkova serrata), walnut (Juglans), horse chestnut (Aesculus turbinata), and so on, for samples of statuses 1-5.
(a) (b)   Stem volume excluding branches and leaves is the so-called "merchantable volume".The average "merchantable volume" (hereinafter "stem volume") per hectare V (m 3 /ha) is obtained from the following inverse formula for the yield-density effect.
where N is the stand density (trees/ha), and b 1 to b 4 are parameters determined for each tree species and each area.The average diameter at breast height D (cm) is obtained from the following equations: (5) where H F is the stand shape ratio, D g is the average section diameter (cm), and α 1 to α 3 and β 1 to β 3 are parameters determined for each tree species and each area.Stand density varies with forest management, but if the stand is not artificially managed, it decreases due to natural mortality.A smaller stand density for the next year based on the natural mortality curve, defined as per Equation ( 8), and maximum density N R f (trees/ha) based on the maximum density curve, defined as per Equation ( 9), are adopted.
where N(t + 1) is the stand density in year t + 1, N 0 is the planted stand density (trees/ha), and v(t) is the average stem volume for a stand (m 3 ) in year t.K 1 to K 4 are determined with the parameters b 1 to b 4 of the inverse formula for the yield-density effect and marginal competition ratio number R f as follows: Figure 2 shows the observed sample forest plots and estimated stand density of Japanese cedar and type L-1 broadleaf forests.The plots with densities less than the natural mortality curves and the maximum density curves indicate locations where densities have decreased due to thinning, etc.

LAI Estimation Model
LAI (leaf area index) is the leaf area per unit area (m 2 /m 2 ).LAI is obtained with the equation for the relationship of dry leaf weight and leaf area and the equation for the relationship of the breast height diameter and dry leaf weight.Those equations are defined for each tree species.Total leaf weight of a tree is more strongly related to the crosssectional area under the living branch than the diameter at breast height, including dead conduits, as shown by Shinozaki et al. [21].Shinozaki's pipe model is a well-known model in the field of ecology, but very few forest survey data are available, and in many cases, equations have been proposed using the breast height diameter as an alternative to the cross-sectional area under the living branch, for example, in [22,23].Therefore, in this study, LAI estimation models using the cross-sectional area under the living branch were developed for Japanese cedar and cypress (Chamaecyparis obtusa), for which survey data are abundant, and models using the breast height diameter for other species as follows: The cross-sectional area under the living branch is related to the height under the branch  (m).Although the height under the living branch changes due to forest management practices such as pruning, in a case where there is no artificial care, it is estimated using the following equation proposed by Kanazawa et al. [24]: where  is the relative spacing index, defined as follows: where  and  are parameters determined for each tree species.Equation ( 14) is modified as follows: Figure 3 shows the results of plotting data extracted from several previous studies [22,[25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].In this study, the parameters  and  of Equation ( 16) are defined as  = 0.02039 and  = 1.384 for cedar and  = 0.01694 and  = 1.439 for cypress from Figure 3.The correlations between the observed and the estimated  , using Equation ( 16) and the obtained parameters, are  = 0.967 and 0.943 for cedar and cypress, respectively, indicating that the parameters were estimated with good accuracy.

LAI Estimation Model
LAI (leaf area index) is the leaf area per unit area (m 2 /m 2 ).LAI is obtained with the equation for the relationship of dry leaf weight and leaf area and the equation for the relationship of the breast height diameter and dry leaf weight.Those equations are defined for each tree species.Total leaf weight of a tree is more strongly related to the cross-sectional area under the living branch than the diameter at breast height, including dead conduits, as shown by Shinozaki et al. [21].Shinozaki's pipe model is a well-known model in the field of ecology, but very few forest survey data are available, and in many cases, equations have been proposed using the breast height diameter as an alternative to the cross-sectional area under the living branch, for example, in [22,23].Therefore, in this study, LAI estimation models using the cross-sectional area under the living branch were developed for Japanese cedar and cypress (Chamaecyparis obtusa), for which survey data are abundant, and models using the breast height diameter for other species as follows: The cross-sectional area under the living branch is related to the height under the branch H B (m).Although the height under the living branch changes due to forest management practices such as pruning, in a case where there is no artificial care, it is estimated using the following equation proposed by Kanazawa et al. [24]: where S r is the relative spacing index, defined as follows: S r = 10000 where a and b are parameters determined for each tree species.Equation ( 14) is modified as follows: Figure 3 shows the results of plotting data extracted from several previous studies [22,[25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].In this study, the parameters a 0 and a 1 of Equation ( 16) are defined as a 0 = 0.02039 and a 1 = 1.384 for cedar and a 0 = 0.01694 and a 1 = 1.439 for cypress from Figure 3.The correlations between the observed and the estimated H B , using Equation ( 16) and the obtained parameters, are R = 0.967 and 0.943 for cedar and cypress, respectively, indicating that the parameters were estimated with good accuracy.Assuming that the trunk shape is trapezoidal, the following equation is proposed between the cross-sectional area under the living branch  (m 2 ) and  [27,41]: where  is the cross-sectional area at breast height (m 2 ), and  is breast height, which is typically 1.3 m [41].In this study, the following equation, modified from Equation (17), was used to adjust for small differences among tree species: where  and  are parameters determined for each tree species.Figure 4 shows the results of plotting data extracted from several previous studies [26,29,31,34,36,[42][43][44][45].In this study, the parameters  and  of Equation ( 18  Many equations have been proposed for estimating dry leaf mass, including one that uses the cross-sectional area and another that uses the square of the diameter.In this study, the following equations, which were relatively accurate for all species analyzed, are proposed: Assuming that the trunk shape is trapezoidal, the following equation is proposed between the cross-sectional area under the living branch A B (m 2 ) and H B [27,41]: where A is the cross-sectional area at breast height (m 2 ), and H 0 is breast height, which is typically 1.3 m [41].In this study, the following equation, modified from Equation ( 17), was used to adjust for small differences among tree species: where c 0 and c 1 are parameters determined for each tree species.Figure 4 shows the results of plotting data extracted from several previous studies [26,29,31,34,36,[42][43][44][45].In this study, the parameters c 0 and c 1 of Equation ( 18 Assuming that the trunk shape is trapezoidal, the following equation is proposed between the cross-sectional area under the living branch  (m 2 ) and  [27,41]: where  is the cross-sectional area at breast height (m 2 ), and  is breast height, which is typically 1.3 m [41].In this study, the following equation, modified from Equation ( 17), was used to adjust for small differences among tree species: where  and  are parameters determined for each tree species.Figure 4 shows the results of plotting data extracted from several previous studies [26,29,31,34,36,[42][43][44][45].In this study, the parameters  and  of Equation ( 18  Many equations have been proposed for estimating dry leaf mass, including one that uses the cross-sectional area and another that uses the square of the diameter.In this study, the following equations, which were relatively accurate for all species analyzed, are proposed: Many equations have been proposed for estimating dry leaf mass, including one that uses the cross-sectional area and another that uses the square of the diameter.In this study, the following equations, which were relatively accurate for all species analyzed, are proposed: where W L is the total dry leaf mass of one tree (kg), and d 0 to d 3 are parameters determined for each tree species.Figure 5 shows the relationships of W L and A B with data extracted from several previous studies [22,[25][26][27][29][30][31][32][33][34]36,[40][41][42][43][44][45][46][47].In this study, the parameters d 0 and d 1 of Equation ( 19) are defined as d 0 = 0.9601 and d 1 = 2.759 for cedar and d 0 = 1.031 and d 1 = 2.897 for cypress from Figure 5.The correlations between A B and W L are R = 0.858 and 0.914 for cedar and cypress, respectively, indicating that the parameters were estimated with good accuracy.
(a) (b) Figure 6 shows the relationships  and  with extracted data from several previous studies for beech, Konara oak, and birch [45,[48][49][50][51][52][53][54][55].Table 1 shows the obtained parameters and their correlations, as determined from Figure 6.Any correlations are more than 0.85, so good parameters were obtained.Many equations have been proposed for estimating the leaf area from the dry weight of the leaf.In this study, the following equations, which were relatively accurate for all species analyzed, are proposed: where A L is the total leaf area of one tree (m 2 ), and a is a parameter determined for each tree species.Figure 7 shows the relationships of A L and W L with data extracted from several previous studies for cedar, cypress, beech, Konara oak, and birch [26,32,33,36,45,49,51,[56][57][58][59][60][61][62][63].
Table 2 shows the parameters obtained and their correlations, as determined from Figure 7.
Any correlations are more than 0.95, so good parameters were obtained.Many equations have been proposed for estimating the leaf area from the dry weight of the leaf.In this study, the following equations, which were relatively accurate for all species analyzed, are proposed: where  is the total leaf area of one tree (m 2 ), and  is a parameter determined for each tree species.Figure 7 shows the relationships of  and  with data extracted from several previous studies for cedar, cypress, beech, Konara oak, and birch [26,32,33,36,45,49,51,[56][57][58][59][60][61][62][63].Table 2 shows the parameters obtained and their correlations, as determined from Figure 7. Any correlations are more than 0.95, so good parameters were obtained.Finally, LAI is obtained using the stand density as follows: where LAI is LAI (m 2 /m 2 ) and N is stand density (trees/ha).As described in this section, the LAI estimation model is proposed based on forest structural factors such as forest age, tree height, and stand density for several tree species, including broadleaf forest species.

Hydrological Model
The time variation of the basin storage volume S is expressed using runoff volume Q, precipitation R, and evapotranspiration E as follows: where t is time.Considering the long-term water balance, dS/dt = 0. Thus, Equation ( 23) If the annual water supply is considered as the water conservation function, the estimation of annual evapotranspiration and observation of annual precipitation are important to evaluate the water conservation function for the specific basin.In Japan, rainfall radar observation networks are well developed, so it is possible to estimate annual basin precipitation amounts.Annual evaporation is estimated using the following procedure, which combines the estimation of evapotranspiration and interception evaporation: Evapotranspiration is obtained using potential evapotranspiration according to the Makkink equation [64], which considers albedo, and correction factors as follows: where E p is daily potential evapotranspiration (mm/d), γ is the psychrometric coefficient, ∆ is the gradient of the saturated vapor pressure curve (hPa/ • C), R s is the total solar radiation (MJm −2 d −1 ), ι is latent heat of evaporation (MJ/kg), α is albedo, and a and b are parameters for each location.The difference between actual and possible evapotranspiration varies depending on the condition of the ground surface, vegetation, and so on.Kondo [65] proposed the following correction factors using LAI: where LAI is LAI (m 2 /m 2 ), and E t is daily actual evapotranspiration (mm/d).Equations ( 25) and ( 26) are used for forest area and grassland, respectively.Interception evaporation is estimated according to the following method proposed by Kondo et al. [66][67][68]: where I is interception evaporation (mm), I POT is the potential interception evaporation (mm), τ is the duration of a rainfall event (h), and S is the water storage capacity of tree surfaces (mm).I POT is defined using the following equation: Water 2024, 16, 588 where R n is the net radiation (W/m 2 ), G is the ground heat flux (W/m 2 ), T is the air temperature (K), rh is relative humidity, and σ is the Stefan-Boltzmann constant.T + , q + , and J are defined using the following equations: where q sat (T) is the saturated specific humidity at air temperature T, u is wind speed, c P is the constant pressure specific heat of air (1006 J/kg/K at 1.0 a.t.m.), ρ is the wet gas density (kg/m 3 ), and C H is the bulk exchange coefficient for sensible heat.The water storage capacity S is defined using the following equation: where BAI is the branch area index (BAI), which is the branch area per unit area (m 2 /m 2 ), and SAI is the stem area index (SAI), which is the stem area per unit area (m 2 /m 2 ).s lea f , s branch , and s stem are the water storage capacity of the leaf, branch, and stem (mm/LAI, mm/BAI, mm/SAI), respectively.BAI was estimated using the cedar and cypress tree shape model proposed by Ebisu et al. [69].For other species, the mean of the parameters for cedar and cypress was used instead.SAI considers the tree shape as a triangular pyramid and calculates the surface area of the triangular pyramid with the cross-sectional area at breast height as the base area and the height of the tree as the height of the triangular pyramid.
Here, the probability of precipitation impacting a tree body Ω * is defined using the following equation: where Ω is the canopy closure rate, and f is the factor of leaf tilt, usually 0.5.In the case of Ω * ×P g > I POT × (τ/24) + S, canopy interception can be estimated as I = I POT × (τ/24) + S; otherwise, I = Ω * ×P g , where P g is gross precipitation (mm).Ω is expressed based on the Beer-Lambert law as follows: where k is the extinction coefficient, affected by the angle of the leaves.In this study, the relationship of 1 − Ω and LAI was plotted in Figure 8 [70,71], and the value of k = 0.4035 (R = 0.892) was obtained.

Carbon Stock Assessment Model
In this study, the carbon stock in living biomass in forest C F (t-C/ha) was calculated using the following equation, based on a greenhouse gas inventory report by the CGER, Japan [72]: where V is the stem volume per hectare (m 3 /ha), D d is the dry wood density (t/m 3 ), BEF is the biomass expansion factor for the conversion of merchantable volume, R root is the rootto-shoot ratio, and CF is the carbon fraction of dry matter (t-C/t).Thus, V•BEF•(1 + R root ) reflects the total volume of trees, including leaves, branches, roots, and stem.Parameters for the major tree species are listed in Table 3.

Carbon Stock Assessment Model
In this study, the carbon stock in living biomass in forest  (t-C/ha) was calculated using the following equation, based on a greenhouse gas inventory report by the CGER, Japan [72]: where  is the stem volume per hectare (m 3 /ha),  is the dry wood density (t/m 3 ),  is the biomass expansion factor for the conversion of merchantable volume,  is the root-to-shoot ratio, and  is the carbon fraction of dry matter (t-C/t).Thus,  •  • 1 +  reflects the total volume of trees, including leaves, branches, roots, and stem.Parameters for the major tree species are listed in Table 3.

Application of Each Species
From the above, a model that evaluates water conservation function, carbon stock, and woody production by calculating forest parameters based on tree species, forest age, stand density, and growth status was developed.Figure 9 shows the simulated changes in LAI, stand density, and tree height with forest age in a case where the initial number of trees planted is 4000 trees/ha and the growth status is 3, using the proposed model in this study.The target area is assumed to be Gujo City, Gifu Prefecture, Japan.The density of cedar trees decreases rapidly due to natural mortality, while the density of other tree species decreases relatively slowly.The LAI of beech has a maximum of about 7, while the

Application of Each Species
From the above, a model that evaluates water conservation function, carbon stock, and woody production by calculating forest parameters based on tree species, forest age, stand density, and growth status was developed.Figure 9 shows the simulated changes in LAI, stand density, and tree height with forest age in a case where the initial number of trees planted is 4000 trees/ha and the growth status is 3, using the proposed model in this study.The target area is assumed to be Gujo City, Gifu Prefecture, Japan.The density of cedar trees decreases rapidly due to natural mortality, while the density of other tree species decreases relatively slowly.The LAI of beech has a maximum of about 7, while the LAI of Japanese cedar and Japanese cypress has a maximum of about 5. On the other hand, the LAI of birch is up to 4. This shows that the LAI varies depending on the species, or even among the same broadleaf trees.
Figure 10 shows the simulated changes in LAI, stand density, and tree height with forest management as 30% thinning at 25, 35, 45, and 55 years, respectively, versus Figure 9. Cedar is less affected by thinning, except for a small increase in LAI when thinned.Thinning of cypress increases LAI.This is the reason why a decrease in the living branch height and an increase in the living branch cross-sectional area are caused by a decrease in tree density.The LAIs of beech and birch decrease with thinning.Although not the same for all forest management and forest conditions, thinning like that for the management of cedar and cypress is excessive for beech and birch.Excessive thinning greatly reduces the number of trees, resulting in a large decrease in LAI.The example analysis presented in this section assumes growth status 3 and the initial planting of 4000 trees/ha and does not simulate a real forest.Nonetheless, the development of this model enables us to evaluate the impact of forest management on a forest with diverse tree species.
LAI of Japanese cedar and Japanese cypress has a maximum of about 5. On the other hand, the LAI of birch is up to 4. This shows that the LAI varies depending on the species, or even among the same broadleaf trees.Figure 10 shows the simulated changes in LAI, stand density, and tree height with forest management as 30% thinning at 25, 35, 45, and 55 years, respectively, versus Figure 9. Cedar is less affected by thinning, except for a small increase in LAI when thinned.Thinning of cypress increases LAI.This is the reason why a decrease in the living branch height and an increase in the living branch cross-sectional area are caused by a decrease in tree density.The LAIs of beech and birch decrease with thinning.Although not the same for all forest management and forest conditions, thinning like that for the management of cedar and cypress is excessive for beech and birch.Excessive thinning greatly reduces the number of trees, resulting in a large decrease in LAI.The example analysis presented in this section assumes growth status 3 and the initial planting of 4000 trees/ha and does not simulate a real forest.Nonetheless, the development of this model enables us to evaluate the impact of forest management on a forest with diverse tree species.Figure 11 shows the simulated changes in the water supply evaluation factors, such as annual water supply, annual evapotranspiration, and annual interception evaporation.The target area is assumed to be Gujo City, Gifu Prefecture, Japan.The water supply for each species is simulated using observation data from the Hachiman gauge station, Japan Meteorological Agency.Meteorological data from 1994, the most severe drought year in Figure 11 shows the simulated changes in the water supply evaluation factors, such as annual water supply, annual evapotranspiration, and annual interception evaporation.The target area is assumed to be Gujo City, Gifu Prefecture, Japan.The water supply for each species is simulated using observation data from the Hachiman gauge station, Japan Meteorological Agency.Meteorological data from 1994, the most severe drought year in the past 30 years, are used.Water supply varies with tree species, but generally decreases until about 40 years old, after which it is near constant.Although beech forests have a larger LAI than needleleaf forests, water supply is not significantly different.The reasons for this are not only differences in LAI but also differences in the water storage capacity per unit area and differences in albedo between the broadleaf forest and needleleaf forest.The modified Makkink equation allows us to estimate that the possible evapotranspiration of the broadleaf forest is less than that of the needleleaf forest due to differences in albedo.When applying this calculation, the albedos of the broadleaf and needleleaf forest are α = 0.20 and 0.15, respectively.Water supply does not change much in cedar, decreases in cypress, and increases in beech and birch as a result of LAI changes due to thinning.Yamamura et al. [5] also reported in their model analysis that the water conservation function decreases with forest growth up to a certain forest age and then tends to be almost constant.Vertessy et al. [73] also reported that the water balance and discharge in the watershed of mountain ash forests after forest burns decreased with forest growth up to a certain age, and it becomes almost constant thereafter.Kojima [74] analyzed daily water discharge over about 20 years in a small forested watershed with a cypress priority, and reported a slight decline trend in discharge in a well-grown forest.Our proposed model reproduced the trend reported in the above studies that "the water conservation function declines up to a certain forest age and then becomes almost constat", which can evaluate water conservation function with high accuracy.Figure 12 shows the simulated changes in the stem volume and carbon stock per unit area.Figure 12 shows that needleleaf forests accumulate a very large amount of carbon stock, while broadleaf forests accumulate less carbon.For other species besides cedar, the accumulated carbon stock is reduced via thinning.When considering forests as carbon sinks, planting cedar trees was found to be the most effective using this model.

LAI Estimation
The model is applied to the Kori River Basin in Omura City, Nagasaki Prefecture.This basin (54.7 km 2 ) has a relatively mild climate with an average annual temperature of about 17 • C and annual rainfall of about 1800 mm, with the Kayase Dam (total storage capacity: 6.8 Mm 3 ) in its upper reaches.Figure 13 shows the distribution of tree species and land cover types in the target basin.The basic forest information, such as species, age, and growth status, was obtained from the Omura City Forest Inventory for FY2020, the National Forest Inventory for FY2021, and Natural Environmental GIS by the Biodiversity Center of Japan [75].The forest conditions from 2006 to 2020 were estimated using the forest model described above, based on the forest management history, such as thinning, clear-cutting, and harvest rate.Figure 14 shows the distribution of LAI in 2006 and 2020, estimated with the forest model.The LAI increases in some areas between 2006 and 2020 due to growth, while the LAI decreases in other areas due to forest management such as logging, clear-cutting, and thinning.
Center of Japan [75].The forest conditions from 2006 to 2020 were estimated using the forest model described above, based on the forest management history, such as thinning, clear-cutting, and harvest rate.Figure 14

Hydrological Application and Validation
Using the hydrological model described above, daily effective rainfall was calculated by subtracting daily interception evaporation and evapotranspiration from total precipitation, and daily discharge at the Kayaba Dam was calculated using a simple runoff Center of Japan [75].The forest conditions from 2006 to 2020 were estimated using the forest model described above, based on the forest management history, such as thinning, clear-cutting, and harvest rate.Figure 14

Hydrological Application and Validation
Using the hydrological model described above, daily effective rainfall was calculated by subtracting daily interception evaporation and evapotranspiration from total precipitation, and daily discharge at the Kayaba Dam was calculated using a simple runoff

Hydrological Application and Validation
Using the hydrological model described above, daily effective rainfall was calculated by subtracting daily interception evaporation and evapotranspiration from total precipitation, and daily discharge at the Kayaba Dam was calculated using a simple runoff model.Figure 15 compares the observed and simulated inflow hydrographs at the Kayaba Dam in 2017 with a drought trend and in 2020 with a heavy rainfall event.Peak flow rates were slightly underestimated, but generally good results were obtained.Table 4 shows validation results based on the Nash-Sutcliffe efficiency (NSE).An NSE of more than 0.7 was obtained for all years, indicating that the flow rate could be calculated with good estimation accuracy.The annual water supply in the Kori River Basin was evaluated based on the LAI distribution for 2006 to 2020.Meteorological conditions were unified for the year 2020.Figure 16 shows the variation in annual water supply.We found that changes in LAI associated with forest management can result in slight variations in water supply with the same meteorological data.However, the change is very small, at less than 1% of the total, suggesting that the forest management practices currently in place do not have a significant impact on water supply.
on the LAI distribution for 2006 to 2020.Meteorological conditions were unified for the year 2020.Figure 16 shows the variation in annual water supply.We found that changes in LAI associated with forest management can result in slight variations in water supply with the same meteorological data.However, the change is very small, at less than 1% of the total, suggesting that the forest management practices currently in place do not have a significant impact on water supply.

Conclusions
In this study, we developed an integrated tool to evaluate how ecosystem services such as the water conservation function, wood production, and carbon absorption are changed by forest management practices such as clear-cutting, thinning, pruning, and tree species shifting for sustainable forest management.Previous studies have classified conifers and broadleaf trees, and evergreen and deciduous trees, but have been insufficient to cover the diversity of tree species in Japan.In Japan, there are many observations and developed models for cedar and cypress for timber production purposes, but models to  on the LAI distribution for 2006 to 2020.Meteorological conditions were unified for the year 2020.Figure 16 shows the variation in annual water supply.We found that changes in LAI associated with forest management can result in slight variations in water supply with the same meteorological data.However, the change is very small, at less than 1% of the total, suggesting that the forest management practices currently in place do not have a significant impact on water supply.

Conclusions
In this study, we developed an integrated tool to evaluate how ecosystem services such as the water conservation function, wood production, and carbon absorption are changed by forest management practices such as clear-cutting, thinning, pruning, and tree species shifting for sustainable forest management.Previous studies have classified conifers and broadleaf trees, and evergreen and deciduous trees, but have been insufficient to cover the diversity of tree species in Japan.In Japan, there are many observations and developed models for cedar and cypress for timber production purposes, but models to

Conclusions
In this study, we developed an integrated tool to evaluate how ecosystem services such as the water conservation function, wood production, and carbon absorption are changed by forest management practices such as clear-cutting, thinning, pruning, and tree species shifting for sustainable forest management.Previous studies have classified conifers and broadleaf trees, and evergreen and deciduous trees, but have been insufficient to cover the diversity of tree species in Japan.In Japan, there are many observations and developed models for cedar and cypress for timber production purposes, but models to evaluate species differences have been lacking for other trees, especially for broadleaf trees.By referring to many previous reports and combining data, we developed a model for multiple tree species, which can evaluate changes that may occur as a result of species conversion and forest management in each species.
The results of our simulations for the major tree species produced the following findings: • Water supply varies with tree species but generally decreases until about 40 years of age, after which it is near constant.

•
Cedar is less affected than other tree species by thinning.The water supply is barely affected by thinning.• The thinning of cypress trees increases the LAI due to reduced density, and water supply is reduced in thinned forests compared to in abandoned forests, influenced by an increased LAI.
• Although beech forests have a larger LAI than needleleaf forests, water supply is not significantly different.

•
Birch has a smaller LAI and a much larger water supply than other species.

•
Broadleaf forests are more affected by thinning than needleleaf forests and tend to receive an increased water supply as a result of thinning.
Our model also enables the evaluation of the water conservation function at the basin scale.Application to the Kori River Basin, Omura City, Nagasaki Prefecture, Japan, showed that large-scale clear-cutting and species conversion may have an impact, but the forest management currently practiced by Omura City does not result in significant changes in forest function.
To summarize, in this study, we developed a model to evaluate forest modification and changes in function for diverse tree species.In the future, we would like to use this model to examine the degree to which forest conditions change and whether that is to the extent that significant alterations in forest function occur, as well as the relationship between the costs of forest management and the ecosystem services that change.1) Equation ( 3) Equation ( 3) Equation ( 3) Equation ( 1) Equation ( 1) Equation ( 2) Equation ( 3) Equation ( 1   The relationship between leaf area and dry leaf mass for each broadleaf tree species is shown in Figure A2.Table A3 shows the model parameters and determination coefficients for each tree species, as obtained from Figure A2.Any correlations are more than 0.85, so good parameters were obtained.However, for red pine and larch, the equations proposed by Kumagai [77] and Fellner et al. [78], respectively, are used.

Figure 1 .
Figure 1.Forest sample plots and growth curves: (a) cedar (Cryptomeria japonica); (b) type L-1 broadleaf forest.Thick solid curves are the growth curves for growth status 3. Other solid curves represent variation in growth by status.

Figure 1 .
Figure 1.Forest sample plots and growth curves: (a) cedar (Cryptomeria japonica); (b) type L-1 broadleaf forest.Thick solid curves are the growth curves for growth status 3. Other solid curves represent variation in growth by status.

Figure 2 .
Figure 2. Observed forest sample plots and estimated .Thick solid curves are the estimated natural mortality curves for growth status 3, and the thin solid curves are the maximum density curves  : (a) Japanese cedar, in the case of  = 4000; (b) type L-1 broadleaf forest, in the case of  = 6000.

Figure 2 .
Figure 2. Observed forest sample plots and estimated N. Thick solid curves are the estimated natural mortality curves for growth status 3, and the thin solid curves are the maximum density curves N R f : (a) Japanese cedar, in the case of N 0 = 4000; (b) type L-1 broadleaf forest, in the case of N 0 = 6000.
) are defined as  = 0.9558 and  = −0.09077for cedar and  = 0.9682 and  = −0.08349for cypress from Figure4.The correlations between the observed and the estimated  are  = 0.980 and 0.983 for cedar and cypress, respectively, indicating that the parameters were estimated with good accuracy.
) are defined as  = 0.9558 and  = −0.09077for cedar and  = 0.9682 and  = −0.08349for cypress from Figure4.The correlations between the observed and the estimated  are  = 0.980 and 0.983 for cedar and cypress, respectively, indicating that the parameters were estimated with good accuracy.

Figure 9 .
Figure 9. Simulated time series variation of LAI, stand density, and tree height in the case of an initial number of planted trees of 4000 trees/ha and growth status of 3 for each tree species: (a) cedar; (b) cypress; (c) beech; (d) birch.

Figure 9 .Figure 10 .
Figure 9. Simulated time series variation of LAI, stand density, and tree height in the case of an initial number of planted trees of 4000 trees/ha and growth status of 3 for each tree species: (a) cedar; (b) cypress; (c) beech; (d) birch.Water 2024, 16, 588 13 of 24

Figure 10 .
Figure 10.Simulated time series variation of LAI, stand density, and tree height with thinning in the case of an initial number of planted trees of 4000 trees/ha and growth status of 3 for each tree species: (a) cedar; (b) cypress; (c) beech; (d) birch.

Figure 11 .
Figure 11.Simulated time series variation of annual water supply  (mm/y), actual evapotranspiration  (mm/y), and interception evaporation  (mm/y) for each tree species in the case of abandoned and thinned forest: (a) cedar; (b) cypress; (c) beech; (d) birch.

Figure 11 .
Figure 11.Simulated time series variation of annual water supply Q (mm/y), actual evapotranspiration E t (mm/y), and interception evaporation I (mm/y) for each tree species in the case of abandoned and thinned forest: (a) cedar; (b) cypress; (c) beech; (d) birch.

Figure 11 .Figure 12 .
Figure 11.Simulated time series variation of annual water supply  (mm/y), actual evapotranspiration  (mm/y), and interception evaporation  (mm/y) for each tree species in the case of abandoned and thinned forest: (a) cedar; (b) cypress; (c) beech; (d) birch.

Figure 12 .
Figure 12.Simulated time series variation of carbon stock C F (tC/ha) and stem volume V (m 3 /ha) for each tree species in the case of abandoned and thinned forest: (a) cedar; (b) cypress; (c) beech; (d) birch.
shows the distribution of LAI in 2006 and 2020, estimated with the forest model.The LAI increases in some areas between 2006 and 2020 due to growth, while the LAI decreases in other areas due to forest management such as logging, clear-cutting, and thinning.

Figure 13 .Figure 14 .
Figure 13.Distribution of tree species and land cover types in the Kori River Basin.

Figure 13 .
Figure 13.Distribution of tree species and land cover types in the Kori River Basin.
shows the distribution of LAI in 2006 and 2020, estimated with the forest model.The LAI increases in some areas between 2006 and 2020 due to growth, while the LAI decreases in other areas due to forest management such as logging, clear-cutting, and thinning.

Figure 13 .Figure 14 .
Figure 13.Distribution of tree species and land cover types in the Kori River Basin.

Figure 15 .
Figure 15.Comparison of observed and simulated inflow hydrographs at the Kayaba Dam: (a) in 2017 with drought trend and (b) in 2020 with a heavy rainfall event.

Figure 16 .
Figure 16.Variation in annual water supply in the Kori River Basin using LAI distribution for 2006 to 2020.Meteorological data fixed in 2020.

Figure 15 .
Figure 15.Comparison of observed and simulated inflow hydrographs at the Kayaba Dam: (a) in 2017 with drought trend and (b) in 2020 with a heavy rainfall event.

Figure 15 .
Figure 15.Comparison of observed and simulated inflow hydrographs at the Kayaba Dam: (a) in 2017 with drought trend and (b) in 2020 with a heavy rainfall event.

Figure 16 .
Figure 16.Variation in annual water supply in the Kori River Basin using LAI distribution for 2006 to 2020.Meteorological data fixed in 2020.

Figure 16 .
Figure 16.Variation in annual water supply in the Kori River Basin using LAI distribution for 2006 to 2020.Meteorological data fixed in 2020.

Table 4 .
Accuracy of simulation results for each year.

Table 4 .
Accuracy of simulation results for each year.

Table 4 .
Accuracy of simulation results for each year.

Table A1 .
List of forest growth model parameters.