Statistical Evidence for Managing Forest Density in Consideration of Natural Volatile Organic Compounds

: Rapid deforestation, coupled with the growing population seeking forest therapy, urges the necessity for research on how to maximize forests’ therapeutic functions when cultivating damaged or unmanaged forests. This study was formulated to provide a basis for forest stand density management to maximize the therapeutic effects of forests with a focus on natural volatile organic compounds (NVOCs), a representative component of forest therapy through analysis of variance and regression analyses. The results of this study revealed all studied stand densities yield the highest total NVOC (TNVOC) emissions in summer, especially in the study site which has a forest density of 700/ha. In addition, treeless areas (0/ha) were found to have the most signiﬁcant difference in average NVOC emissions when cultivated at a density of 700/ha. When managing forests with a density of 900/ha to 1000/ha, it has been shown that it is most desirable, in terms of therapeutic function efﬁciency, to control a density of 500/ha to 700/ha. Finally, regression equations for the ﬁve experimental sites with signiﬁcant explanatory power were derived. Based on the results of the conducted analyses, 700/ha of forest density is recommended to maximize the therapeutic effects of forests, compared to other ranges of forest density. and 700/ha, are located within 100 m of the center of the 600/ha site, while the 900/ha and 1000/ha study sites are located within 100 m of each other. The distance between the 500/ha site and the 1000/ha site is about 2 km. The average age of P. koraiensis trees in the site were 31 to 40 years, with an average diameter at breast height of 27 cm and an average height of 21 m. In addition, dominant wind directions in the study sites during the study period were south-southeast (700/ha, 900/ha, and 1000/ha), south (0/ha and 600/ha), and south-southwest (500/ha).


Introduction
Since 2010, there has been a steady loss of forests globally at a rate of ten million hectares per year [1]. When considering the amount of newly cultivated forest areas and combining it with the total loss, then the rate is only brought down to 4.7 million hectares per year [1]. With this continuing rate of loss each year, countries are now focusing their efforts on forest cultivation to grow healthy and superior forests. Researchers are conducting diverse studies on various ways to afforest forests suitable for their needs and purposes, such as forests that maximize wood production or therapeutic effects [2][3][4][5]. The cultivating forestry plantations are projects to cultivate and grow forests so that artificial forests and natural forests can grow to be healthier and more superior than untouched forests. Depending on the age and condition of the forest, several things can be performed, such as weeding, planting young trees, thinning, natural forest nurturing, and removing winders [6]. Forest cultivation has several positive effects, such as enhancing water supply functions, increasing carbon dioxide absorption, improving the ecological health of forests, raising the economic value of forests, and preventing wildfires [6][7][8][9][10][11][12][13][14][15]. Additionally, forests grown through the cultivating forestry plantation projects have a positive function for human health. Previous studies have shown that tended forests have more of a positive physical and mental impact on visitors when compared to wild forests. Furthermore, visitors to the tended forests experience increased metabolic activity and increased physical vitality [16,17]. Prior studies of forest thinning, among several forest management methods, have shown that it is a viable method for expanding the benefits that a forest provides to both human health and the natural ecosystem. The forest thinning projects are well known for its positive effects, such as improving forest habitat value, increasing forest carbon absorption and storage rates, and facilitating forest water supply by controlling the leaf area index; thus, allowing forests to respond more flexibly to climate change and preventing forest fires from spreading rapidly [18][19][20][21][22]. As such, forest stand density control is beneficial for forest ecosystems and its economic value protection. Furthermore, density control through forest thinning was found to have positive effects on human health which include physiological relaxation, stress reduction, anxiety and anger relief, depression improvement, vitalization, emotional stability, and mental health recovery [23][24][25][26][27][28][29][30][31][32]. Therefore, it is expected that the adjustment of forest density through forest cultivation will maximize the forest healing effects suitable for the characteristics and needs of participants who partake in forest therapy programs.
Forest therapy includes a plethora of factors, and phytoncide is a forest component that, according to current research, appears to be relevant in the health effects induced by forest exposure [33,34]. Phytoncide is a complex mixture of chemicals secreted by trees to protect themselves from external harmful factors and is also a representative therapeutic factor in nature [35][36][37]. Phytoncide can be classed as a natural volatile organic compounds (NVOC), and terpenes, a class of chemicals found in phytoncide mixtures, are found within NVOCs. Terpenes are notable for improving people's health, which occurs when individuals are exposed to the atmosphere where terpenes are scattered. Phytoncide has been shown to benefit both physiological and psychological health [38][39][40][41][42][43][44]. Phytoncide increases activities of natural killer (NK) cells [38], has antioxidant and antibacterial effects [39], and lowers blood pressure, pulse rate, and cortisol levels [40][41][42]. As for the psychological effects, it reduces stress and improves one's mood [40,[42][43][44]. When considering the effects of the factors that make up phytoncide, alpha-pinene accounts for most of the monoterpene produced in forests, and several prior studies have demonstrated anti-inflammatory, antioxidant, and antianxiety effects [45][46][47]. In addition, several monoterpenes, including limonene and beta-pinene, have been proven to have a positive effect when it comes to diseases by increasing immunity, alleviating cardiovascular disease, and improving depression [48][49][50][51][52][53][54]. As such, phytoncide is attracting attention as an important factor in forest therapy due to its various health benefits, with studies also being conducted on NVOCs emitted by forests.
NVOCs, which have these various healing effects, have been identified through several prior studies. Studies measuring and analyzing forest emissions of NVOCs have shown that monoterpenes are significantly correlated with temperature and solar radiation, with highest concentrations in the summer [33,55]. In addition, phytoncide concentrations are generally the highest in the early morning and noon, with diurnal cycle [33,34,[55][56][57]. On the other hand, the emission characteristics of NVOCs vary depending on the density of forests; there were several previous studies analyzing the relationship between forest density control and NVOC emissions conducted for similar purposes as the aim of this study. Their findings show that when forest density in dense forests was adjusted, monoterpene emissions increased significantly, and the highest concentration of phytoncide has been found to be emitted from forest environments that are neither too dense nor too open [58][59][60]. As mentioned above, since NVOCs have different emissions and compositions depending on the forest stand density, several studies have conducted experiments attempting to link specific forest density and natural environments with their corresponding forests' therapeutic effects. However, most of the relevant studies presented were only one-off studies, and phytoncide was measured at a low frequency of about once a season or once a day for one range of forest density. Furthermore, studies conducted in conjunction with changes in forest density through forest cultivation and phytoncide emissions are very rare, and studies with a focus on microclimate factors and phytoncide's link to forest density are also rarely conducted. As mentioned earlier, forest density plays a vital role in creating forest therapy spaces, and the size of a forest's therapeutic effect varies greatly depending on the level of forest density. However, exactly what mechanism or causative material increases a forest's therapeutic effects after the forest density has been controlled has not yet been identified. Therefore, this study aims to demonstrate the relationship between forest density control and forest therapy effects by conducting long-term follow-up studies on several ranges of forest density, and from those findings, select the most suitable forest density for cultivation of forestry plantation, with a focus on phytoncide. As a result of this study, a forest density of 700/ha is the most suitable forest density to maximize the therapeutic effects of forests when focusing on NVOC emissions.

Study Site
The study was conducted in Pinus koraiensis Siebold et Zucc. forests of Gwangneung Experimental Forest, managed by National Forest Research Institute located in Pocheon, Gyeonggi-do, South Korea. The geographical locations of the study sites are 37 • 47 12.786 N, 127 • 11 21.0876 E, and they are located at 150 m to 170 m above sea level. Gwangneung Experimental Forest has a total area of 1072 ha, of which 515 ha of the coniferous forest, including P. koraiensis forests, accounts for 48% of the total area. Gwangneung Experimental Forest was designated as a UNESCO (United Nations Educational, Scientific and Cultural Organization) Biosphere Reserve in 2010 and is recognized worldwide for its importance [61]. In addition, Gwangneung Experimental Forest is home to more than a thousand species of plants and over 3000 species of animals [62]. In this study, the research was conducted on five experimental areas with 500, 600, 700, 900, and 1000 tree counts per hectare of P. koraiensis forests, and control areas with zero tree counts per hectare outside the research site ( Figure 1). The study sites, including the control site, 500/ha, 600/ha, and 700/ha, are located within 100 m of the center of the 600/ha site, while the 900/ha and 1000/ha study sites are located within 100 m of each other. The distance between the 500/ha site and the 1000/ha site is about 2 km. The average age of P. koraiensis trees in the site were 31 to 40 years, with an average diameter at breast height of 27 cm and an average height of 21 m. In addition, dominant wind directions in the study sites during the study period were south-southeast (700/ha, 900/ha, and 1000/ha), south (0/ha and 600/ha), and south-southwest (500/ha).
phytoncide's link to forest density are also rarely conducted. As mentioned earlier, forest density plays a vital role in creating forest therapy spaces, and the size of a forest's therapeutic effect varies greatly depending on the level of forest density. However, exactly what mechanism or causative material increases a forest's therapeutic effects after the forest density has been controlled has not yet been identified. Therefore, this study aims to demonstrate the relationship between forest density control and forest therapy effects by conducting long-term follow-up studies on several ranges of forest density, and from those findings, select the most suitable forest density for cultivation of forestry plantation, with a focus on phytoncide. As a result of this study, a forest density of 700/ha is the most suitable forest density to maximize the therapeutic effects of forests when focusing on NVOC emissions.

Study Site
The study was conducted in Pinus koraiensis Siebold et Zucc. forests of Gwangneung Experimental Forest, managed by National Forest Research Institute located in Pocheon, Gyeonggi-do, South Korea. The geographical locations of the study sites are 37° 47' 12.786'' N, 127° 11' 21.0876' E, and they are located at 150 m to 170 m above sea level. Gwangneung Experimental Forest has a total area of 1072 ha, of which 515 ha of the coniferous forest, including P. koraiensis forests, accounts for 48% of the total area. Gwangneung Experimental Forest was designated as a UNESCO (United Nations Educational, Scientific and Cultural Organization) Biosphere Reserve in 2010 and is recognized worldwide for its importance [61]. In addition, Gwangneung Experimental Forest is home to more than a thousand species of plants and over 3000 species of animals [62]. In this study, the research was conducted on five experimental areas with 500, 600, 700, 900, and 1000 tree counts per hectare of P. koraiensis forests, and control areas with zero tree counts per hectare outside the research site ( Figure 1). The study sites, including the control site, 500/ha, 600/ha, and 700/ha, are located within 100 m of the center of the 600/ha site, while the 900/ha and 1000/ha study sites are located within 100 m of each other. The distance between the 500/ha site and the 1000/ha site is about 2 km. The average age of P. koraiensis trees in the site were 31 to 40 years, with an average diameter at breast height of 27 cm and an average height of 21 m. In addition, dominant wind directions in the study sites during the study period were south-southeast (700/ha, 900/ha, and 1000/ha), south (0/ha and 600/ha), and south-southwest (500/ha).

Measurement Methods
In this study, NVOC emissions and microclimate environment factors data were measured in six sites (five experimental and one control) to determine the impact of forest stand density control on phytoncide and microclimate factors. In order to reduce measurement errors for each survey site, experiments were conducted at three different plots with a radius of 10 m of each survey site. This study was conducted from September 2015 to December 2017, and measurements were taken twice a month. NVOCs were measured from March to December, and microclimate environment factors were measured from March to September due to the winter entry policies of the studied forest. Detailed indicators measured are shown in Table 1, and a total of 17 NVOC compounds were selected out of 19 total detected NVOC compounds. Two compounds, α-Cedrene and Geranyl acetate, were excluded because the concentration was too small; therefore, most collected samples could not detect it. NVOCs were collected every two hours from morning (9:00) to sunset (17:00) in consideration of the peak visitation times for the forests. According to geographical characteristics, NVOC measurements were carried out at three different circular plots located within five meters of the central tree at each of the six study sites. Five pumps were installed per each plot, considering the vegetation characteristics. The adsorption tube method was used to collect the samples. Adsorption was carried out in tubes containing 150 mg of Tenax TA from Markes (Sacramento, CA, USA). NVOC was measured using a µg/m 3 unit, where m 3 represents the volume of the measurement sites' surrounding environment. The total air volume of 9 L collected over an hour was converted to m 3 , and a detailed description follows. The sample capture system was a mini pump (MP-30KN, Sibata, Japan), and the calibration was preceded by an adsorption error measurement before using the flow meter. A total volume of 9 L of NVOC was collected at a flow rate of 150 mL/min. A previous study focusing on the sampling amount, conducted to make NVOC measurements in the forest more accurate and efficient, discovered that 9 L produced the most efficient results when compared to 1, 3, 6, 9, 12, 24, and 48 L; thus, this study was also conducted with an air volume of 9 L [63]. The sampling equipment was mounted on a tripod 1.5 m above the ground, and the average value was calculated by repeating the process at each site. To avoid artificial errors when in contact with the tube, disposable polyethylene gloves and antibacterial masks were used during the experiment. The sample tubes were kept at a temperature below 4 • C for 48 h after sampling before being analyzed ( Table 2). To reduce the possibility of error, the value of the tubes collected without Tenax TA inside were also reflected.
The samples were analyzed qualitatively and quantitatively using a gas chromatographymass spectrometer (7890N-5975, Agilent, Santa Clara, CA, USA) equipped with a thermal desorption device (GC-MSD, Gerstel TDS, Gerstel, Germany). The adsorption tube substances are concentrated in a low-temperature cryofocusing system, which takes highpurity helium gas from a thermal desorption device at a rate of 1 mL/min. The system desorbed the gas for 3 min at 210 • C while maintaining a temperature of −30 • C. After, the compounds were then heated for 3 min at 220 • C before being injected into a GC spectrometer and detected with an MSD. Temperature, humidity, wind speed, solar radiation, and photosynthetically active radiation (PAR) were measured at each of the six study sites in order to collect microclimate environment factors. Additionally, the site's direction and slope were calculated in terms of locational environment. Using a portable multifunction meter, the physical characteristics of the site environment were recorded at 5 min intervals (HOBO-U23 V2, Onset, Bourne, MA, USA). Solar radiation sensors (S-LIB-M003, Onset, Bourne, MA, USA) and photosynthetically active radiation sensors (S-LIA-M003, Onset, Bourne, MA, USA) were placed in the same location and tracked throughout the experiment.
In order to obtain wind velocity data at each designated site, a wind monitoring sensor (Wind Monitor O5103-45, R.N.Y., Logan, UT, USA) was also mounted in consideration of the geological features. The meter was balanced at a height of 1.5 m on a tripod approximately 5 m away from a mini pump, and the digitalized measurement results were saved and converted for the study. The findings were evaluated using the HOBO-ware Pro software (Onset, Bourne, MA, USA). To reduce the possibility of measurement errors, data saved 5 min before and after each measurement were excluded from the study.

Calibration Curve
A few measures can be used to validate both the analysis device and the procedures. The calibration curve was created using 20 different species of standard chemicals, including α-pinene and β-pinene. Using the calibration curve to calculate each element's mass number and the square of its rate of diluting standard materials, the majority of the materials have a linearity greater than 0.997. Examples include α-pinene (R 2 = 0.997), β-pinene (R 2 = 0.998), and d-limonene (R 2 = 0.999). Experiments with these materials are highly reproducible in terms of the linear correlation coefficient, making them ideal for research.

Analysis Methods
The analysis was conducted by using R 4.0.5 and RStudio. "DescTools", "lmtest", "gvlma", "lm.beta", "pheatmap", "PMCMRplus" and other R packages were used for the analysis. In this study, data were collected for a total of 17 NVOC compounds and 5 microclimate environmental factors, with samples of 2351 and 47,810, respectively ( Figure 2). The data observed at three different measurement plots for each forest density were aggregated into one, and the data was used for the analysis of the results through preprocessing processes such as missing value and outlier elimination. For NVOCs, we calculated total NVOC (TNVOC) values that could represent 17 different NVOC compounds and carried out descriptive statistics including the average and standard deviation of NVOCs and microclimate environmental factors. aggregated into one, and the data was used for the analysis of the results through prepro cessing processes such as missing value and outlier elimination. For NVOCs, we calcu lated total NVOC (TNVOC) values that could represent 17 different NVOC compound and carried out descriptive statistics including the average and standard deviation o NVOCs and microclimate environmental factors. In order to understand the NVOC emission characteristics of P. koraiensis forests, an nual TNVOC emissions at the site were examined, and differences of annual TNVOC emissions by forest density were also investigated. In addition, microclimate environmen analysis of P. koraiensis forests was carried out according to forest density. Subsequently correlation analyses between TNVOC and microclimate factors by forest density was con ducted to explore the relationship between NVOCs and microclimate factors by fores density of P. koraiensis forests.
Next, to find out the effect of forest density on TNVOC emissions, a one-way analysis of variance (one-way ANOVA) was conducted. Since there were more than 30 samples the normality test was omitted by the central limit theorem, and the Bartlett's test of ho mogeneity of variance was performed, but the result did not satisfy the homoscedasticity Therefore, Welch's one-way ANOVA was executed, and results confirmed significant dif ferences between forest densities. Accordingly, Dunnett's T3 test for post hoc analysis was conducted to determine the statistical significance between the control group with zero tree counts per hectare and experimental groups with 500, 600, 700, 900, and 1000 tree counts per hectare. In order to understand the NVOC emission characteristics of P. koraiensis forests, annual TNVOC emissions at the site were examined, and differences of annual TNVOC emissions by forest density were also investigated. In addition, microclimate environment analysis of P. koraiensis forests was carried out according to forest density. Subsequently, correlation analyses between TNVOC and microclimate factors by forest density was conducted to explore the relationship between NVOCs and microclimate factors by forest density of P. koraiensis forests.
Next, to find out the effect of forest density on TNVOC emissions, a one-way analysis of variance (one-way ANOVA) was conducted. Since there were more than 30 samples, the normality test was omitted by the central limit theorem, and the Bartlett's test of homogeneity of variance was performed, but the result did not satisfy the homoscedasticity. Therefore, Welch's one-way ANOVA was executed, and results confirmed significant differences between forest densities. Accordingly, Dunnett's T3 test for post hoc analysis was conducted to determine the statistical significance between the control group with zero tree counts per hectare and experimental groups with 500, 600, 700, 900, and 1000 tree counts per hectare.
Finally, multiple regression analysis using the bidirectional procedure method was conducted to investigate the effects of microclimate factors on phytoncide emissions by forest density, and a phytoncide emission prediction equation was formulated using micro-climate factors by forest density. It is also expected that the general public will be easily able to calculate the amount of phytoncide emissions of the selected study site with forest density that can maximize the forest healing effects. Multiple regression analyses were performed for each forest density, totaling six times. Model 1, a regression analysis containing all five microclimate factors, and Model 2, a regression analysis conducted by selecting only significant microclimate factors through bidirectional procedures, were established, and F-tests were conducted to determine whether microclimate factors removed from Model 2 had statistically significant effects. Multicollinearity tests were also conducted to determine the correlation between independent variables in Model 2, and finally, Durbin-Watson statistics were performed to verify autocorrelation. Based on the regression results validated in this way, we establish a final regression equation by forest density.

Characteristics of NVOCs at P. Koraiensis Forests by Forest Density
The average monthly total NVOC (TNVOC) emissions from P. koraiensis forests are shown in Figure 3. The highest concentration of TNVOC was found to be emitted in summer from all the surveyed sites, and July recorded higher phytoncide emissions than other months. In addition, as a result of adding up monthly TNVOC emissions from each survey sites, 0, 500, 600, 700, 900, and 1000 tree count per hectare were 1.84 µg/m 3 , 2.92 µg/m 3 , 2.79 µg/m 3 , 3.72 µg/m 3 , 1.87 µg/m 3 , and 1.41 µg/m 3 , respectively. The standard deviation of the above calculated results of each survey sites were 0.19 µg/m 3 , 0.25 µg/m 3 , 0.26 µg/m 3 , 0.23 µg/m 3 , 0.14 µg/m 3 , and 0.11 µg/m 3 , respectively. Annual TNVOC emissions in 700/ha were the highest, followed by 500/ha, 600/ha, 900/ha, 0/ha, and 1000/ha. In general, it can be seen that the yield of NVOC emission is high at forest densities between 500 and 700 tree count per hectare.
performed for each forest density, totaling six times. Model 1, a regress taining all five microclimate factors, and Model 2, a regression analysis lecting only significant microclimate factors through bidirectional proc tablished, and F-tests were conducted to determine whether microcl moved from Model 2 had statistically significant effects. Multicollineari conducted to determine the correlation between independent variables finally, Durbin-Watson statistics were performed to verify autocorrelati regression results validated in this way, we establish a final regression e density.

Characteristics of NVOCs at P. Koraiensis Forests by Forest Density
The average monthly total NVOC (TNVOC) emissions from P. kor shown in Figure 3. The highest concentration of TNVOC was found to be mer from all the surveyed sites, and July recorded higher phytoncide emi months. In addition, as a result of adding up monthly TNVOC emission vey sites, 0, 500, 600, 700, 900, and 1000 tree count per hectare were 1.84 μ 2.79 μg/m 3 , 3.72 μg/m 3 , 1.87 μg/m 3 , and 1.41 μg/m 3 , respectively. The sta of the above calculated results of each survey sites were 0.19 μg/m 3 , μg/m 3 , 0.23 μg/m 3 , 0.14 μg/m 3 , and 0.11 μg/m 3 , respectively. Annual TNV 700/ha were the highest, followed by 500/ha, 600/ha, 900/ha, 0/ha, and 10 it can be seen that the yield of NVOC emission is high at forest densities 700 tree count per hectare.  Figure 4, which explains the monthly deviation in TN at the surveyed sites by forest density, shows the largest deviation in mo emissions at the 700/ha site. In addition, the regions of 0, 900, and 100 hectare were found to have relatively small deviations compared to othe emissions at the 700/ha site. In addition, the regions of 0, 900, and 1000 tree count per hectare were found to have relatively small deviations compared to other study sites.
Atmosphere 2021, 12, x FOR PEER REVIEW 8 of Regarding the types of NVOCs detected for each study site, generally α-Pinene, Ca phene, and β-Pinene were detected at high concentrations ratio as shown in Figure 5. Sa inene and Myrcene were detected at higher rates in the 0/ha region than in the other gions, and Farnesene was detected at higher rates in the 700/ha region than in other gions. In addition, while p-Cymene was detected in other sites, no p-Cymene was detect in 900/ha areas, and no 3-Carene and d-Limonene were detected in 1000/ha sites. At t 500/ha, 600/ha, and 700/ha sites, all 17 types of NVOCs were detected, although the pr portions of phytoncide components were different. Regarding the types of NVOCs detected for each study site, generally α-Pinene, Camphene, and β-Pinene were detected at high concentrations ratio as shown in Figure 5. Sabinene and Myrcene were detected at higher rates in the 0/ha region than in the other regions, and Farnesene was detected at higher rates in the 700/ha region than in other regions. In addition, while p-Cymene was detected in other sites, no p-Cymene was detected in 900/ha areas, and no 3-Carene and d-Limonene were detected in 1000/ha sites. At the 500/ha, 600/ha, and 700/ha sites, all 17 types of NVOCs were detected, although the proportions of phytoncide components were different.

Characteristics of Microclimate Environments at P. Koraiensis Forests by Forest Density
The results of examining fluxes in forest density and monthly microclimate factors in P. koraiensis forests are described in Figure 6. In the case of temperature, there was no significant difference depending on forest density, and in case of humidity, the 700/ha surveyed site was significantly higher than other research sites with different forest density. In the case of solar radiation, it tended to appear at high levels in areas with low forest density, 0/ha and 500/ha, while the 700/ha area was notably high in May and June. Regarding levels of photosynthetically active radiation, 600/ha and 700/ha sites generally showed lower values than other study sites. In addition, wind speeds were generally higher in areas with low forest density, except for April and June.
inene and Myrcene were detected at higher rates in the 0/ha region than in the other regions, and Farnesene was detected at higher rates in the 700/ha region than in other regions. In addition, while p-Cymene was detected in other sites, no p-Cymene was detected in 900/ha areas, and no 3-Carene and d-Limonene were detected in 1000/ha sites. At the 500/ha, 600/ha, and 700/ha sites, all 17 types of NVOCs were detected, although the proportions of phytoncide components were different.

Characteristics of Microclimate Environments at P. Koraiensis Forests by Forest Density
The results of examining fluxes in forest density and monthly microclimate factors in P. koraiensis forests are described in Figure 6. In the case of temperature, there was no significant difference depending on forest density, and in case of humidity, the 700/ha surveyed site was significantly higher than other research sites with different forest density. In the case of solar radiation, it tended to appear at high levels in areas with low forest density, 0/ha and 500/ha, while the 700/ha area was notably high in May and June. Regarding levels of photosynthetically active radiation, 600/ha and 700/ha sites generally showed lower values than other study sites. In addition, wind speeds were generally higher in areas with low forest density, except for April and June.

Correlation Analysis of NVOCs and Microclimate Environments
The Pearson correlation coefficient analysis between TNVOC emissions and microclimate factors by forest density is shown in Figure 7. The TNVOC emissions data and the microclimate environmental data were individually measured at each site surveyed (0/ha, 500/ha, 600/ha, 700/ha, 900/ha, and 1000/ha). After processing the collected data in units

Correlation Analysis of NVOCs and Microclimate Environments
The Pearson correlation coefficient analysis between TNVOC emissions and microclimate factors by forest density is shown in Figure 7. The TNVOC emissions data and the microclimate environmental data were individually measured at each site surveyed (0/ha, 500/ha, 600/ha, 700/ha, 900/ha, and 1000/ha). After processing the collected data in units of one hour, the analysis was performed separately for each study site. Temperature was significantly positively correlated in 500/ha, 900/ha, and 1000/ha study sites. The control site with zero tree count per hectare showed no strong correlation with any microclimate factors, and PAR and wind speed did not show much correlation in 900/ha and 1000/ha experimental sites. In addition, temperature and humidity generally showed positive tendencies in all the study sites except the control site. Solar radiation, PAR, and wind speed were noted to manifest negative tendencies in experimental sites.
tmosphere 2021, 12, x FOR PEER REVIEW experimental sites. In addition, temperature and humidity generally tendencies in all the study sites except the control site. Solar radiation speed were noted to manifest negative tendencies in experimental sites.

One-Way Analysis of Variance of NVOCs and Forest Density
One-way analysis of variance (one-way ANOVA) was performed to NVOC emissions from P. koraiensis forests in one region were the same re density control as shown in Table 3. The study site with no tree count control group, and the forest density-controlled regions were set as exp (500/ha, 600/ha, 700/ha, 900/ha, and 1000/ha). To more accurately analy sions, the analysis was conducted on 2351 TNVOC raw data points col survey period, instead of monthly phytoncide emissions by survey site. of samples is more than thirty, normality tests were omitted by central li homogeneity test was performed prior to the one-way ANOVA. Bartle geneity of variances rejected the null hypothesis with p < 0.001 and was moskedasticity. As a result, Welch's one-way ANOVA was carried out cant level of p < 0.001, the null hypothesis was rejected, indicating tha

One-Way Analysis of Variance of NVOCs and Forest Density
One-way analysis of variance (one-way ANOVA) was performed to answer whether NVOC emissions from P. koraiensis forests in one region were the same regardless of forest density control as shown in Table 3. The study site with no tree count, 0/ha, was set as control group, and the forest density-controlled regions were set as experimental groups (500/ha, 600/ha, 700/ha, 900/ha, and 1000/ha). To more accurately analyze TNVOC emissions, the analysis was conducted on 2351 TNVOC raw data points collected during the survey period, instead of monthly phytoncide emissions by survey site. Since the number of samples is more than thirty, normality tests were omitted by central limit theorem, and homogeneity test was performed prior to the one-way ANOVA. Bartlett's test of homogeneity of variances rejected the null hypothesis with p < 0.001 and was not satisfied homoskedasticity. As a result, Welch's one-way ANOVA was carried out, and at a significant level of p < 0.001, the null hypothesis was rejected, indicating that the population means among independent variables were not all equal. A post hoc analysis was conducted to determine the differences between specific groups, and Dunnett's T3 test was performed. By analyzing the differences between control groups and experimental groups, only the 700/ha experimental site was found to have a significant difference in NVOC emissions between itself and the control group (0/ha). Furthermore, significant differences in NVOC emissions were found when 900/ha and 1000/ha of dense forests were adjusted to a density of 500/ha to 700/ha.

Regression Analysis of NVOCs and Microclimate Environments by Forest Density
In order to discover the impact of microclimate factors on NVOC emissions and to figure out which microclimate factors significantly affect phytoncide emissions for each forest density, an NVOCs prediction equation was established. Through this prediction equation, it is expected that the general public can easily calculate the phytoncide emissions of the selected study site with an optimized forest density that would maximize the forest healing effects.

Control Site (0 Tree/ha)
Prior to the analysis of the results, Durbin-Watson statistics were calculated to verify its D-W value, and the resulting value was 1.43, close to 2.00, and it was determined to be suitable for a multiple regression analysis model. However, it was not possible to derive a significant regression equation with a significance probability of 0.645 for Model 1, and the explanatory power of which the dependent variable is explained by the independent variable was also very low at 9.7%. For Model 2, where Model 1 was processed according to the bidirectional procedure method, the significant regression equation could not be derived with a significance probability of 0.162, and the explanatory power was very low at 5.9%. Furthermore, since there was only one remaining independent variable in Model 2, the process of verifying multicollinearity among independent variables was omitted, and it was concluded that significant regression analysis could not be derived by the multiple regression model.

Experimental Site 1 (500 Trees/ha)
Regression results for the first experimental group, the 500/ha region, are shown in Table 4. Prior to the analysis of the results, Durbin-Watson statistics were calculated to verify its D-W value, and the resulting value was 1.57, close to 2.00, and it was determined to be suitable for a multiple regression analysis model. The explanatory power of Model 2 was 58.2%, and the significance probability of Model 2 was 0.036, confirming that at least one independent variable had a significant effect on the dependent variable. Furthermore, a multicollinearity test on Model 2 confirmed that all independent variables have a variance inflation factor (VIF) of less than 10, indicating there is no multicollinearity. Among the independent variables included in the regression model, temperature (p < 0.05) had a significant effect on the dependent variable, TNVOC emissions, in the 500/ha region. In addition, regression analysis shows that higher temperature (B = 0.047) result in higher TNVOC emissions at the surveyed site, and higher solar radiation (B = −0.009) results in lower TNVOC emissions. Among the independent variables of Model 2, temperature and solar radiation, temperature (β = 0.568) was shown to have a greater effect on TNVOC emissions in the 500/ha site. The regression equation of Model 2 for prediction of TNVOC emissions in the 500/ha study site is shown below. (1)

Experimental Site 2 (600 Trees/ha)
Regression results for the second experimental group, the 600/ha region, are shown in Table 5. Prior to the analysis of the results, Durbin-Watson statistics were calculated to verify its D-W value, and the resulting value was 1.52, close to 2.00, and it was determined to be suitable for a multiple regression analysis model. The explanatory power of Model 2 was 56.6%, and the significance probability of Model 2 was 0.043, confirming that at least one independent variable had a significant effect on the dependent variable. Furthermore, a multicollinearity test on Model 2 confirmed that all independent variables have a variance inflation factor (VIF) of less than 10, indicating there is no multicollinearity. Among the independent variables included in the regression model, temperature (p < 0.05) had a significant effect on the dependent variable, TNVOC emissions, in the 600/ha region.
In addition, regression analysis shows that higher temperature (B = 0.012) and humidity

Experimental Site 3 (700 Trees/ha)
Regression results for the third experimental group, the 700/ha region, are shown in Table 6. Prior to the analysis of the results, Durbin-Watson statistics were calculated to verify its D-W value, and the resulting value was 1.77, close to 2.00, and it was determined to be suitable for a multiple regression analysis model. The explanatory power of Model 2 was 76.4%, and the significance probability of Model 2 was 0.017, confirming that at least one independent variable had a significant effect on the dependent variable. Furthermore, a multicollinearity test on Model 2 confirmed that all independent variables have a variance inflation factor (VIF) of less than 10, indicating there is no multicollinearity. Among the independent variables included in the regression model, temperature and humidity (p < 0.05) had a significant effect on the dependent variable, TNVOC emissions, in the 700/ha region.
In addition, regression analysis shows that higher temperature (B = 0.041) and humidity (B = 0.029) result in higher TNVOC emissions at the surveyed site. Among the independent variables of Model 2, temperature and humidity, humidity (β = 0.472) was shown to have a greater effect on TNVOC emissions in the 700/ha site. The regression equation of Model 2 for prediction of TNVOC emissions in the 700/ha study site is shown below. TNVOC = −0.410 + 0.041 × (Temperature) + 0.029 × (Humidity)

Experimental Site 4 (900 Trees/ha)
Regression results for the fourth experimental group, the 900/ha region, are shown in Table 7. Prior to the analysis of the results, Durbin-Watson statistics were calculated to verify its D-W value, and the resulting value was 2.07, close to 2.00, and it was determined to be suitable for a multiple regression analysis model. The explanatory power of Model 2 was 70.3%, and the significance probability of Model 2 was 0.003, confirming that at least one independent variable had a significant effect on the dependent variable. Furthermore, a multicollinearity test on Model 2 confirmed that all independent variables have a variance inflation factor (VIF) of less than 10, indicating there is no multicollinearity. Among the independent variables included in the regression model, temperature (p < 0.01) and wind speed (p < 0.05) had a significant effect on the dependent variable, TNVOC emissions, in the 900/ha region. In addition, regression analysis shows that higher temperature (B = 0.046), solar radiation (B = 0.020), and PAR (B = 0.007) result in higher TNVOC emissions at the surveyed site, and higher wind speed (B = −0.664) results in lower TNVOC emissions. Among the independent variables of Model 2, temperature, wind speed, solar radiation, and PAR, temperature (β = 0.731) was shown to have a greater effect on TNVOC emissions in the 900/ha site. The regression equation of Model 2 for prediction of TNVOC emissions in the 900/ha study site is shown below. TNVOC = −1.192 + 0.046 × (Temperature) − 0.664 × (Wind Speed) + 0.020 × (Solar Radiation) + 0.007 × (PAR) (4)  Regression results for the last experimental group, the 1000/ha region, are shown in Table 8. Prior to the analysis of the results, Durbin-Watson statistics were performed to verify its D-W value, and the resulting value was 1.96, close to 2.00, and it was determined to be suitable for a multiple regression analysis model. The explanatory power of Model 2 was 94.6%, and the significance probability of Model 2 was 0.000, confirming that at least one independent variable had a significant effect on the dependent variable. Furthermore, a multicollinearity test on Model 2 confirmed that all independent variables have a variance inflation factor (VIF) of less than 10, indicating there is no multicollinearity. Among the independent variables included in the regression model, temperature (p < 0.001), wind speed (p < 0.001), solar radiation (p < 0.01), and PAR (p < 0.001) had a significant effect on the dependent variable, TNVOC emissions, in the 1000/ha region. In addition, regression analysis shows that higher temperature (

Discussion
This study was conducted with the aim of finding suitable forest density for maximizing the beneficial effects of forests on human health in the new afforestation of forests or cultivated forests devastated by rapid urbanization and industrialization. As revealed in previous papers, the cultivated forests are said to be of greater physiological and psychological help to humans, compared to wild forests [16,17], and they also show different therapeutic effects depending on the extent of their cultivation [23][24][25][26][27][28][29][30][31][32]. At a time when the number of people visiting forests to experience their therapeutic function is increasing, due to the COVID-19 pandemic, it is essential to optimize the therapeutic function through forest therapy programs in properly cultivated forests [64][65][66][67]. However, despite this situation, forest density control methods have rarely been studied with a focus on forests' therapeutic functions, and only few studies have investigated healing factors for various forest densities, researched relationships with forest microclimate factors, and suggested suitable forest densities. Therefore, to compensate for the research gaps and limitations of these prior studies, we collected and analyzed NVOCs and microclimate factor data for forestry environments with six different ranges of forest density and present the most optimal forest density for forest therapy.
This study was conducted based on a total of 17 NVOCs and five microclimate factor data collected from September 2015 to December 2017. The TNVOC analysis, which combined 17 types of NVOC emissions, showed that it was the highest in the summer, especially in July. Additionally, α-Pinene, Camphene, and β-pinene were detected at higher rates in every study site. These three compounds were also emitted at higher rates in other prior studies which researched volatile organic compounds under the forest canopy [33,57,68,69]. Additionally, the annual TNVOC emissions were the highest in 700/ha regions, which leads one to infer that the therapeutic effect through phytoncide would be the highest in 700/ha regions. However, since the deviation of phytoncide emissions in the above survey site was the largest, further investigation into phytoncide emissions in the 700/ha area is needed. Observing the changes in the annual microclimate factor, solar radiation was high in the forests with low forest density. This may be due to the low density of the forest, resulting in a greater amount of solar radiation reaching the ground. Wind speeds were also high in areas with low forest density, which is expected due to the low density allowing winds to pass quickly through the forests without lingering. According to the results of the Pearson correlation coefficient analysis, there were correlations between microclimate factors and TNVOC emissions by forest density, and several significant results could be obtained. Temperatures were shown to have significant correlations between 500/ha, 900/ha, and 1000/ha regions. The control site without trees (0/ha) had no significant correlation with any microclimate factor. This result is also consistent with prior studies which demonstrated relationships between temperature and terpenes emissions under forest canopy [33,56,57,[70][71][72]. Meanwhile, in the Pearson correlation coefficient analysis, wind speed showed slight negative tendencies with phytoncide emissions in all study sites, and in July, wind speed in all study sites decreased dramatically. Thus, as shown in Figure 4, it can be inferred that the increase in phytoncide emissions in July may have been affected by wind speed. This study establishes that there is a relationship between forest density control and the extent of therapeutic effects, based on one-way ANOVA and post hoc analyses conducted to determine changes in TNVOC emissions due to forest density control. Analyses found that when conducting a forest cultivation project on bare land without any trees (0/ha), it was most desirable to cultivate up to a density of 700/ha. Furthermore, considering the therapeutic effects of forests, cultivating forests down to levels of 500/ha to 700/ha is desirable when it comes to denser forests (900/ha, 1000/ha). In some previous studies, changes of physiological and psychological health conditions of visitors to forests of different ranges of forest density were measured to reveal the therapeutic effects of forests according to forest stand densities. However, in most cases, surveyed forests with different forest stand density in these studies did not comprise the same species and trees of the same age, and there was a limitation in that the measurement was only one-off [73][74][75]. Therefore, this result is expected to be used as an important basis for presenting the appropriate directions of forest density control in the planning stage of future forest cultivation projects that aim to maximize the forest's therapeutic functions. Finally, this study analyzed microclimate factors that significantly affect TNVOC emissions by forest density and performed regression analyses to predict TNVOC emissions by forest density. As a result, the study failed to derive significant regression equations for the control site, the 0/ha region. However, for the five experimental sites of 500/ha, 600/ha, 700/ha, 900/ha, 900/ha, and 1000/ha, regression equations with explanatory power of 58.2%, 56.6%, 76.4%, 70.3%, and 94.6%, respectively, were derived. In all five surveyed sites, temperature had a significant positive impact on TNVOC emission forecasts, and in all four sites except 700/ha, temperature had the largest impact compared to other microclimate factors. The 600/ha and 700/ha regions had significant positive effects on humidity, which is consistent with prior research that focused on different forest densities of 500/ha, 600/ha, and 700/ha [60]. One interesting result was that wind speeds in the 900/ha and 1000/ha regions had a significant negative effect on TNVOC emission forecasts, the second largest after temperature. This can be found in the report regarding guidelines for the management and creation of fine dust-reducing urban forests distributed by the Korea Forest Service [76]. According to the report, forests with densities of 800/ha to 1000/ha are said to maximize the length of time that wind circulates in the forest by properly inducing the airflow in forests. In addition, there are several prior studies that proved that dense forests significantly reduce wind speed and increase turbulence in the air [77,78]. Therefore, the 900/ha and 1000/ha surveyed sites in this study are expected to be affected by the prolonged wind retention and turbulence, making it difficult to capture phytoncide and dispersing the released phytoncide. In conclusion, these derived predictive equations are expected to enable forest visitors to easily calculate phytoncide emissions by measuring several microclimate factors without directly measuring phytoncide in complex and challenging ways. In particular, the equation for predicting the amount of phytoncide in the 700/ha survey site, which was found to have excellent forest healing effects, is expected to be of great assistance for follow-up studies.
It is important to note that this study has several limitations. First, this regression is derived from data that measured P. koraiensis forests in Pocheon, South Korea, and the threshold remains that predictive accuracy can be reduced in areas with different climates and different physical or geographical elements. Second, as many other factors besides microclimate elements can affect NVOC emissions of P. koraiensis forests, outside factors also need to be considered, including NVOCs transportation through wind flow within forests, soil nutrient availability, herbivore damage to plants, and microbial interactions. Third, it is limiting that the NVOCs and microclimate data of P. koraiensis forests were not measured continuously, but rather from March to December and March to September, respectively. In addition, since only few studies have recognized changes in NVOC emissions and their relationship with microclimate factors by forest density, there was a limit to comparing the results of this study with prior research. In addition, 900/ha and 1000/ha regions are relatively far from other study sites; therefore, there can be a limitation that the microclimate environment in these areas may have been notably different from other study sites. Lastly, there is also a limitation that as NVOCs and microclimate factors data were not measured directly by controlling the forest density of each site surveyed, but by setting up areas with different forest densities in the same forest as experimental groups, the results may differ from the data measured immediately after direct forest cultivation; thus, in subsequent studies, data collection on different forest types in varied regions after direct forest cultivation, with continuous observations throughout every month, will allow for a more expansive study. Additionally, measuring the difference in NVOC emissions by further segmenting the forest density will enable to provide deeper insights. Furthermore, since this study identified different microclimate factors which significantly affect NVOC emissions at each forest density range, it is expected that further studies will be able to consider these specific microclimate factors at each forest density in their research (instead of all microclimate factors). Considering forest density control methodology for maximizing forest healing functions, it is suggested that 700/ha is the optimal density for forest therapy compared to other ranges of forest density.
The modern world, which has serious problems with massive forests being devastated and disappearing every year, is focusing its efforts on forest restoration and forest cultivation. In addition, as the recent COVID-19 pandemic has led to a significant increase in the number of forest visitors to experience forests' therapeutic functions, research on how to restore and cultivate damaged and disappeared forests to maximize forest healing functions is urgent; yet very few studies have been conducted to link these forest therapeutic functions with forest cultivation methods. Therefore, based on NVOCs and microclimate data for the different ranges of forest densities, this study, which derived prediction equations of phytoncide emissions for each forest density and suggested suitable forest density based on statistical analyses, is significant in that it provides directions and guidelines for maximizing forest therapy functions in future forest cultivation trials. It is also expected to be used as a basis for designing adequate locations for future forest therapy programs.

Conclusions
This study aimed to analyze NVOC emissions and microclimate factors according to forest density and determine the forest density that is most suitable for forest therapy when considering phytoncide emissions among various other forest therapy components. This study was conducted in the P. koraiensis forest in Pocheon, South Korea, where the area with 0 tree count per hectare was set up as a control group, and five study sites of 500, 600, 700, 900, and 1000 tree counts per hectare were set as experimental groups. Between September 2015 and December 2017, 2351 NVOCs data and 47,810 microclimate data were collected, based on a total of 17 NVOCs data and five microclimate indicators. Based on the collected data, the analysis of TNVOCs resulted in the release of the highest concentration of phytoncide in the summer from all the target sites, and the highest annual phytoncide emissions was observed at 700/ha out of the six study sites (including both control and experimental groups). The Pearson correlation analysis indicated that temperatures were positively correlated in study sites of 500/ha, 900/ha, and 1000/ha. The results of a one-way ANOVA of the collected data showed the most significant difference in average NVOC emissions, which was found when a treeless area (0/ha) was cultivated at a density of 700/ha, and when managing dense forests (900/ha, 1000 ha), it has been shown that it is most desirable in terms of therapeutic function to control densities of 500/ha to 700/ha. As a result of regression analysis for building regression equations that can predict TNVOC emissions and finding out microclimate indicators that significantly affect TNVOC emissions by forest density, for the five experimental sites of 500/ha, 600/ha, 700/ha, 900/ha, and 1000/ha, regression equations with explanatory power of 58.2%, 56.6%, 76.4%, 70.3%, and 94.6%, respectively, were derived. Based on these analysis results, considering the therapeutic function of forests focused on phytoncide, 700/ha of forest density is expected to maximize the therapeutic effects of forests, compared to other forest densities. In addition, this study is expected to be used as an important basis and guideline for controlling forest density when restoring devastated forests or creating forest therapy spaces in the future.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available since they are collected for internal research purposes of the National Institute of Forest Science, Republic of Korea.