Optimizing Stand Spatial Structure by Using Neighborhood-Based Quantitative Indicators: A Case of Boreal Forests

: 12 Background: Over the past fifty years, societies have placed increasing demands on forests, and their use has shifted 13 gradually from wood production to the diversified benefits and functions of ecosystem services. The effects of 14 neighborhood-based structural characteristics on regulating growth and promoting sustainability have therefore drawn 15 much attention. However, direction for managing natural mixed forests using neighborhood-based indicators are still 16 not clear. 17 Methods: In this study, a tree-level harvest planning tool that considers four neighborhood-based structural parameters 18 (species mingling, diametric differentiation, horizontal spatial pattern and crowdedness of trees) while concurrently 19 recognizing other operational constraints, was developed using simulated annealing algorithm. The approach was 20 applied to four 1-ha mapped stands in northeast China, namely a natural larch forest (NLF), a natural birch forest (NBF), 21 a natural secondary forest (SEF), and a Korean pine broad-leaved forest (KBF). 22 Results: The tree-level harvest optimization improved the objective function values by approximately 78.33% of NLF, 23 and 134.96% of NBF, and 156.70% of SEF and 252.95%, respectively. The optimal harvest intensities for partial cutting 24 activities varied from 22.16% (SEF) to 26.07% (NBF) of the standing volume. In evaluating the four 25 neighborhood-based structural parameters, both species mingling and crowdedness have the highest priority to be used 26 in structure-based forest management. 27 Conclusions: Our results demonstrated that that the commonly used neighborhood-based structural parameters could be 28 used to control the spatial layout of potential harvest trees, in turn may be conducive to regulate the growth and stability 29 of forests. 30

include the uniform angle (W), mingling index (M), dominance index (D), and crowdedness index (C). The 48 measurements of spatial relationships reflect the distribution, diversity, size variation, juxtaposition of trees, 49 respectively, within a unit of area (Hui and Gadow, 2003). The selection of neighbouring trees (e.g., how many 50 neighbors) is a critical issue in the measurement of stand spatial structure (Wang et al., 2016). Methods have been 51 devised to arrive at the reasonable number of neighbouring trees, such as nearest neighbour, fixed radii and 52 Voronoi tessellations (Pommerening and Stoyan, 2006). However, Hui and Gadow (2003) and Wang et al. (2016) 53 both indicated that four neighbours were enough for these measurements, when one considers sampling accuracy 54 and costs simultaneously. Expanding from the tree-level to the stand-level, average values of these measures can 55 help explain the structural differences between stands composed of different forest types, positioned at different 56 successional stages, and influenced by different management activities (e.g., Wan et al., 2019). 57 58 Figure 1 Schematic diagram of a structural unit, in which N 1~N4 are the four nearest neighbours of a reference 59 tree, N 5~N6 are also close to the reference tree but are not the nearest neighbours, α 1~α4 are the angles between 60 6 (123°20'-124°21'E, 52°16'-52°47'N) in the cold-temperate zone, where the mean annual air temperature is -2.8°C 105 and the mean annual precipitation level is about 450 mm. Precipitation events occurs primarily in summer in this 106 part of the Heilongjiang Province. Due to the higher latitude of this location, as compared to the other case studies, 107 species abundance is relatively lower. Here, larch and birch species usually coexist in NLF forests, while birch is 108 usually most frequent tree species in NBF forests. Other tree species that may be found in these forest types 109 include spruce (Picea asperata), Scots pine (Pinus sylvestris), and poplar (Populus davidiana). 110 The third forest is composed of Korean pine (Pinus koraiensis) and other broadleaved tree species (KBF). 111 This type of forest is reflective of the historical climax community of the Xiaoxing'an Mountains 112 (129°11′-129°26′E, 46°32′-46°39′N). The forest is situated in an area that has a temperate continental monsoon 113 climate, where the mean annual air temperature is 2°C, and the mean annual precipitation level is 600 mm. The 114 research plot that was measured in the third forest has been subjected to different intensities of commercial 115 thinning activities. These were conducted in the 1990s by local communities for the purpose of generating income, 116 however this forest has been developing naturally since. The vegetation in the research plot consists of more than 117 20 tree species that are native Xiaoxing'an Mountain flora. The main species are Korean pine (P. koraiensis), oak 118 (Quercus mongolica), birch (B. platyphylla), other hardwoods (e.g., Fraxinus mandshurica, Phellodendron 119 amurense, Juglans mandshurica, and Ulmus pumila) and softwood species (e.g., P. davidiana). 120 The fourth forest is representative of a naturally regenerated second growth forest (SEF), and is situated in a 121 low-altitude area of the Zhangguangcai Mountains (127°18′-127°41′E, 45°2′-45°18′N). The climate here is 122 temperate, with a mean annual air temperature of 3.1°C and a mean annual precipitation level of 700 mm. 123 Following a clearfelling activity in 1920s, the research site was naturally colonised by trees; few management 124 activities have occurred here since regeneration began. However, to promote regeneration in forest gaps, Korean 125 pine seedlings were planted in 2004 at a density of approximately 500 trees per hectare. Species abundance in this 7 study area is variable, and can also be as large as 20. F. mandshurica, P. davidiana, U. pumila, P. amurense, P.  To conduct the optimization simulation, a 1-ha (100 m×100 m) permanent plot was established within each 131 of the case study forests in the summer of 2016-2017. All trees with a diameter at breast height (DBH) ≥5 cm 132 were marked, and their locations, species, DBH, total tree height (HT), crown width (CW) in four directions (i.e., 133 east, west, south and north), living branch height, and status (e.g., survival, diseases or insects) were recorded. The 134 four case study areas were composed of forests that had tree densitities ranging from about 800 to 1400 trees per 135 hectare, average diameters ranging from 10.5 to 15.8 cm, and tree heights ranging from about 10 to 13 m ( Spatial structure parameters 139 The spatial structural measurements for each structural unit (each reference tree, or averaged to the entire stand) in which are usually defined as the proportion of n nearest neighbours that are a different species from each 144 reference tree. The dominance index characterizes the size (e.g., diameter, height) differentiations between each 145 reference tree and its n nearest neighbours, which can reflect the competition status between the reference tree and 146 its n nearest neighbours. The uniform angle index is used to measure the regularity degree of the n nearest 147 neighbours around each reference tree. The crowdedness index reflects the overlapping between the canopy of 148 each reference tree and its n nearest neighbours. The formulations of the four indices were stated in Dong et al. 149 (2020) as: 150 where M i , D i , W i and C i are the values of mingling degree, diameter dominance, uniform angle, and crowdedness 157 for reference tree i, respectively; n is the number of neighbours of reference tree i, where n=4 in this analysis; and 158 v ij , k ij , z ij and y ij are discrete variables. Further, sp i and sp j represent the tree species of reference tree i and 159 neighbour tree j, respectively; dbh i and dbh j represent the DBH of reference tree i and its neighbour tree j, 160 respectively; and L ij represents the Euclidean distance between the boles of reference tree i and its neighbour tree j. 161 The following logic was employed in equations 1-4: 162 163  If neighbour tree j is not the same tree species as reference tree i, v ij =1, otherwise v ij =0 164  If the DBH of neighbour tree j is smaller than the DBH of reference tree i, k ij =1, otherwise k ij =0 165  If the angle α of two neighbour trees is smaller than the expected standard angle α 0 , z ij =1, otherwise z ij =0 166  If the sum of the crown width of neighbour tree j (CW j ) and reference tree i (CW i ) is less than the Euclidean 167 distance between the boles of neighbour tree j and reference tree i (L ij ), y ij =1, otherwise y ij =0 diversity, a disadvantaged growth status, clumped tree patterns, and a dense growth space. To avoid the edge 180 effects (Pommerening and Stoyan, 2006), trees located within the core area of each study site were considered as 181 either reference trees or neighbour trees, while trees located around the outer edges (with 5 m of the boundary) of 182 each study site were only considered to be neighbour trees. 183

Optimization formulations 185
As was suggested in section 2.2, the ecological complexity and character of a forest usually increases with 186 increases of species mingling (M), and decreases with increases of diameter dominance (D) and crowdedness (C). 187 However, the contribution of uniform angle index (W) has a reversed-U form, namely the closer the W value is to 188 0.5, the more reasonable the stand structure would be. Thus, the tree-level operational plan for a single stand is 189 designed as a non-linear optimization problem that maximizes the average structural characteristics of the four

Simulated annealing algorithm 231
Since the objective function and constraints are all computed and evaluated in a post-harvest manner, the 232 neighbors of each reference tree may always be in a state of change (Figure 1), therefore it would be 233 computationally burdonsome to know every potential structural unit a priori. Thus, the optimization procedures 234 13 should dynamically obtain the correct combination of neighboring trees to construct the potential structural unit 235 prior to solving the model. Thus, a s-metaheuristic algorithm (simulated annealing) was employed to solve the 236 non-linear integer planning problem, rather than an exact method (mixed integer programming) where these 237 potential relationships need to be known and listed in the problem formulation. 238 Simulated annealing (SA) was initially described by Kirkpatrick et al. (1983) although its foundations can be 239 traced back several decades earlier (Metropolis et al., 1953). As a heuristic methods for solving combinatoiral 240 forest planning problems, SA has been shown to be able to produce results very close to optimal (e.g., mixed value) according to a given cooling rate, which acts to decrease the probability of accepting non-improving moves 250 as the search progresses. Near the end of the heuristic search process, T is so low that the probability of 251 acceptance for non-improving moves is nearly 0, and the process acts as a hill-climbing heuristic. 252 In our implementation of SA, the search process is initially guided by a user-defined harvesting intensity (i.e., 253 30% of the trees) from which a feasible solution (i.e., a set of assigned harvest trees) would be generated randomly. 254 From this point, three alternative moves could be made to diversify the potential harvesting intensities. The first 255 strategy randomly schedules a previously unscheduled tree to the current harvest plan, which will increase the 256 14 harvest intensity. These second strategy randomly unschedules a previously scheduled tree from the current 257 harvest plan, which will decrease the harvest intensity. For the third strategy, a tree previously scheduled for 258 harvest would be randomly replaced by a previously unscheduled selected tree, which keeps the harvest intensity 259 consistent. The selection of these three strategies during a run of the SA algorithm was also random. 260 The tree-level harvest planning formulations and simulated annealing algorithm were programmed in R 261 environment (R Core Team, 2021). The procedures was run on a 2.6 GHz Core i5 processor with 4GB RAM and 262 Windows 7 operating system. Parameters that were used to control the search processes of SA include the starting 263 temperature, stopping temperature, cooling rate and the number of iterations at each temperature. Based on a set 264 of trial-and-error tests, the values were set as 100, 1, 0.999 and 10, respectively. These parameters imply that each 265 independent run will result in approximately 46 031 potential iterations (assessments of potential moves). Ten 266 solutions were generated for each forest and only the best solution (as identified as having the highest objective 267 function value) was employed to perform an analysis of differences amongst the case study forests. 268

269
The SA algorithm produced results that indicate the objective function values of each type of forest increased

283
The neighborhood-based structural characteristics for each of the tested forests before-and after-thinning 284 indicate that parameter M seemed to be the most important contributors to the objective function values for NBF 285 forests (Table 3) between the values of W after thinning and an expect value (namely 0.5) were decreased slightly for NBF (0.4720 295 vs 0.4818, in a term of before-vs after-thinning) and KBF (0.5070 vs 0.4976), while were increased for NLF 296 16 (0.4904 vs 0.4874) and SEF (0.5107 vs 0.4856). However, the values of W for the four forests implied that the 297 horizontal distribution of the remaining trees all belonged to the categories of random, especially for the NBF that 298 were adjusted from a regular to a random distribution. 299 300

302
Since the Q-values were only employed as constraints rather than objectives, the variations of diameter 303 distribution were all not well-marked as expected (Table 3). However, slight decreases in the Q-values were 304 17 observed for all of the four forests, indicating the increases on the proportion of smaller trees, yet all forests still 305 had Q-values in reasonable ranges for natural uneven-aged forests. The assigned harvests mainly focused on the 306 trees with DBH less than 15 cm, namely 76.38% of NBF, and 61.31% of KBF, and 88.49% of NLF, and 74.35% 307 of SEF, respectively (Figure 4).

314
Spatial distribution of the assigned harvest trees of optimal solutions for the four forests illustrates that the 315 assigned harvest trees usually had the same tree species and dense crowdedness degrees within the structural units 316 ( Figure 5). The statistical results indicated that the assigned harvest trees mainly focused on B. platyphylla (97%) 317 of NBF, while L. gmelinii were dominated the harvest trees of NLF (80%). However, the distribution on tree 318 18 species of the assigned harvest trees for SEF and KBF were much more complex. The top five species scheduled 319 for harvest in the SEF were F. mandshurica (26.44%), and P. koraiensis (24.90%), and P. davidiana (14.94%), and 320 U. pumila (8.81%), and P. amurense (6.13%). The top five species scheduled for harvest in the KBF were P. 321 koraiensis (40.83%), and B. platyphylla (12.43%), and A. fabri (8.88%), and other softwoods (6.51%), and Q.

330
One purpose of forest management is to generate or maintain a complex stand structure and cultivate a 331 NLF NBF SEF KBF