Modelling the Spatial Structure of White Spruce Plantations and Their Changes after Various Thinning Treatments

: Research Highlights: The spatial distribution of trees results from several ecological processes that can be difﬁcult to measure. We applied a point process modelling approach that uses the diameter and species of neighbouring trees to represent inter-tree interactions through repulsive and attractive processes. Thinning treatments slightly inﬂuence the tree spatial distribution of trees in white spruce plantations. Integrating this “spatialiser” into growth models could help improve stand simulations following various thinning treatments over larger areas and longer periods. It could also allow for the use of spatially explicit models when tree position is not available. Background and Objectives: Tree spatial patterns result from several ecological processes and have important implications in forest ecology and management. The use of spatial information can signiﬁcantly improve our understanding of forest structures. However, this implies intensive ﬁeld work that is rarely integrated into forest inventories. The aims of this study were to develop a spatial distribution simulator of trees in white spruce plantations and to evaluate the inﬂuence of thinning treatments. Materials and Methods: A point process modelling approach was used to represent inter-tree interactions through repulsive and attractive process in white spruce ( Picea glauca (Moench) Voss) plantations in eastern Quebec, Canada, that had been commercially thinned ﬁve years ago. Balsam ﬁr ( Abies balsamea (L.) Mill.) and hardwoods together can represent 30–40% of the basal area of these plantations. Results: The diameter and species of each tree’s two closest neighbours were found to be the most important predictors in explaining the observed distances between trees. Despite the short period since thinning treatments, results showed that the treatment had slight signiﬁcant effects on tree interactions. However, their impact on the global spatial distribution of stands is quite limited. Conclusions: Using only a few readily-available variables (species and diameter of trees), this “spatialiser” will make it possible to assign spatial coordinates to trees and generate realistic stand spatial structures even after various silvicultural treatments.


Introduction
The distance between a tree and its closest neighbours has important implications in forest ecology and management [1,2]. As forest dynamics (growth, mortality, and recruitment) are influenced by competitive processes between neighbouring trees for access to resources [3], the use of spatial information can significantly improve the understanding of forest structure [4]. However, since obtaining tree spatial coordinates requires intensive field work in addition to the usual field measurements (e.g., tree diameter, height, and species), spatial variables are rarely integrated into forest inventories. Assigning spatial managed and unmanaged forests [28,29,36,37]. In addition to their wood production objective, treatments are meant to restore or maintain some of the declining species such as white spruce [34]. Schütz [38] showed that stand conversion can be achieved through a distinct succession of steps that include (1) modifying stand dynamics by altering competition, (2) promoting long-lived species regeneration by creating local opening, (3) managing structural development, and (4) repeating silvicultural interventions over time to create and maintain an irregular/uneven-aged structure. However, the effects of these steps need to be studied in the short and long term. In the Bas-Saint-Laurent region, the management of a portion of the current even-aged stands will shift toward an irregular/uneven-aged path. For this purpose, commercial thinning approaches of various intensities were tested.
The aim of this study was to develop a spatial distribution simulator of trees in white spruce plantations that would work according to a point process approach in order to generate a precise stand spatial structure and to integrate parameterised models for all species. Our first objective was to analyse and to model local intra-and interspecific spatial interactions. Our second objective was to evaluate the influence of thinning treatments on spatial structure and then to use the "spatialiser" to simulate stands treated with various thinning approaches and to determine how these treatments immediately affect stand structure.

Study Site
The plantations used in this study are located in the eastern balsam fir-yellow birch bioclimatic subdomain of the boreal mixed wood zone [39] in eastern Quebec, Canada. In this area, forests are characterised by mixed stands of yellow birch (Betula alleghaniensis Britton), balsam fir, white spruce, and eastern white cedar [40]. Some other hardwoods that can be locally observed include trembling aspen (Populus tremuloides Michx.), paper birch (Betula papyrifera Marshall), red maple (Acer rubrum L.), and sugar maple (Acer saccharum Marshall). A total of 66 sample plots were used to model tree spatial distribution. These plots are part of 2 commercial thinning trials that compare various treatments (TrT): a thinning from below (TrT 1/3 ), in which the smallest trees were cut while ensuring equal spacing between the remaining trees; a crop tree release thinning (TrT CT ), in which the competition was removed 3 m around either 50 or 100 dominant trees per hectare [41]; and a thinning with priority selection of balsam fir (TrT BF ), in which all balsam fir trees were harvested. Both experiments also included an untreated control (TrT 0 ).

Data Acquisition
Considering the large area of the plots, the quantity of trees to be positioned, and the difficulty of manually positioning them, the use of terrestrial lidar was preferred to obtain the most precise spatial coordinates of all trees (including saplings). Terrestrial lidar is a powerful technology for capturing the three-dimensional structure of forests with a high level of precision [42,43]. All plots were scanned immediately after treatment with a Faro ® Focus 3D terrestrial laser scanner (TLS). Thirteen scans were done per plot to obtain full plot coverage and to minimise occlusion. Spherical targets (at least 3 seen from each scan position) were placed throughout the plot and used for scan registration. These multiple scans were then aligned to create a single 3-dimensional point cloud with the Faro ® Scene 5.0 software (FARO Technologies, Rugby, Warwickshire, UK).
We used the point cloud to map the coordinates of each tree in each plot (Figure 1). To do so, we imported the point clouds in the R programming environment [44]. We produced a digital terrain model from the points with the lowest altitude and reconstructed the ground surface with a Delaunay triangulation using the lidR package [45]. Then, we extracted a 20 cm slice of points centred at 1.30 m above the ground. To eliminate returns due to the branches, the 20 cm slice was rasterised on the XY plan into 1 × 1 cm rasters. We took advantage of the fact that rasters with points from a stem have a higher density than those with branches to derive a threshold for cleaning the data from raster that did not contain stems. The remaining points were clustered using density-based algorithms implemented in the DBSCAN package [46]. Finally, we automatically adjusted a circle to the points in the cluster to measure the diameter at breast height (DBH; measured at 1.3 m) and to obtain the XY coordinates of all the trees with a DBH greater than 1 cm. We corrected some DBH mis-adjustments (e.g., to add missing trees or delete clusters not belonging to trees) following a manual visual validation. Since no efficient algorithms to recognise species via TLS were available, species were conventionally tallied and then merged into 3 groups (Gsp): white spruce (WS), balsam fir (BF), and merchantable hardwoods (HW). The HW group included paper birch, yellow birch, sugar maple, red maple, ash (Fraxinus sp.), and trembling aspen. Some rare species (e.g., eastern white cedar, Thuja occidentalis L.) and non-commercial hardwoods (e.g., elderberry (Sambucus nigra ssp. canadensis (L.) R. Bolli), rowan (Sorbus aucuparia L.), and speckled alder (Alnus incana ssp. rugosa (Du Roi) R.T. Clausen)) were not considered in the analysis. A total of 10,696 trees were kept for the analysis. In order to deal with border effects, we used all border trees (i.e., those located within 1.5 m of the plot edges) as neighbours to calculate the distances between neighbours. These distances were also used to calculate the Clark and Evans index [6,47] (Equation (1)), which describes the aggregation of trees at the stand level. This index was calculated for each Gsp (R WS , R BF , and R HW ) and all trees (including saplings and merchantable trees) (Table 1). A value below 1 indicated a tendency to cluster, a value close to 1 as associated with a random distribution, and a value above 1 indicated a regular distribution.
where N represents the number of trees for a given Gsp and r i is the distance between tree i and its closest neighbour of the same Gsp. Stand characteristics are summarized in Table 1. The variables used in this study and their abbreviations are listed in Table 2.

Modelling Spatial Stand Structure
We used a spatial distribution simulator to describe the spatial stand structure in WS plantations. This simulator includes models at the tree level describing intra-and interspecific interactions between neighbouring trees. These models are presented in the following paragraphs.
For each model developed in this study, we tested a list of potential stand-level explanatory variables including silvicultural treatment, aggregation index, and stand density, as well as potential tree-level variables including DBH, Gsp, and interactions between them. Variables that were highly correlated were excluded (variance inflation factor (VIF) > 10 [48]).
We used the glm stats library of the R statistical programming environment [44] to fit the models, and we selected variables using a backward elimination process based on the corrected form of Akaike's information criterion (AICc) [49].

Clark and Evans index
At the stand level, spatial structure can be summarised by the Clark and Evans index, which provides a single value by Gsp following a linear model (Equation (2)).
where X is the design matrix of covariates, ω is the fixed-effect parameter vector to be estimated, and ε is the error term.

Spatial Interactions between Individual Trees
At the individual tree level, spatial structure was characterised on the basis of interactions between neighbouring trees. We calculated the distance between a tree and its neighbours by a nearest neighbour analysis using the RANN package [50], from which we obtained the distance to the 2 closest neighbours (MinDist T1 and MinDist T2 ). Then, we modelled these 2 distances as a function of the available inventory stand variables and of those from the targeted tree and its 2 nearest neighbours. As these distances are continuous and positive variables, we used a gamma regression model [51] (Equation (3)).
where Y is the distance to a neighbour (MinDistT1 or MinDistT2), Z is the design matrices of covariates, γ is the fixed-effect parameter vector to be estimated, and ε is the error term.

Spatial Interactions within Species Groups
For Gsp with a R Gsp ≥ 1, no intraspecific interaction was modelled. However, for those Gsp that tended to form clusters (R Gsp < 1), we first determined the number of clusters observed in a stand using the method of affinity propagation clustering [52] (Figure 2). To describe these aggregations, we characterised 2 variables: the number of clusters per Gsp per hectare (NClusHa Gsp ) and the distance between trees of the same Gsp within a cluster. As NClusHa Gsp is a count variable (positive and discrete), we chose a Poisson regression model [51]. The general model equation form was similar to Equation (3), where Y is the number of clusters per hectare and per Gsp (NClusHa Gsp ), Z is the covariate design matrix, and γ the fixed effect parameter vector to be estimated.
Within a cluster, we computed the minimal distance between a tree and its closest neighbour of the same Gsp (Dist Gsp ). We calculated this distance by a nearest neighbour analysis using the RANN package [50]. Since Dist Gsp is a continuous and positive variable, A gamma regression was choose to model Dist Gsp [51].

Model Validation
We used a repeated 5-fold cross validation to calculate the predictive accuracy of the models [53]. The stands were randomly split into 5 subsets of equal size: the first four were used to calibrate the model, and the fifth was used to validate it. The root mean square error (RMSE) and the determination coefficient (R 2 ) of the predictions were then calculated using the validation subset. The process was repeated until all subsets were used to validate the model. The random segregation of the plots into the subsets was repeated 50 times. The validation result was the average of all the repetitions. After building these models, we developed a spatial stand structure generator (hereafter referred to as a "spatialiser") to generate coordinates for a non-spatial inventory. This "spatialiser" generates random coordinates that are evaluated to determine if positioning is possible based on the inter-and intraspecies models presented above. To avoid edge effects (i.e., trees close to the plot border having no competitors on the exterior side), the simulated plot has to be embedded in a torus in order to keep the same plot area and have no border effects. We followed these steps to generate a list of XY coordinates for each tree: 1.
Trees were ordered by DBH.

2.
For the Gsp where R Gsp < 1, the number of clusters (NClusHa Gsp ) was calculated.

3.
For all trees, we repeated the following steps, depending on Gsp.
• For trees with R Gsp ≥ 1, we generated a random position, identified the 2 closest competitors and, knowing the characteristics (Gsp and DBH) of the neighbours, calculated the theoretical minimum distance with these 2 trees (i.e., MinDist T1 and MinDist T2 ). If the 2 measured distances were greater than the 2 theoretical values, the point was kept as a potential valid position.

•
For trees with R Gsp < 1, if the number of trees of the same Gsp already positioned was less than the predicted NClusHa Gsp , the tree was positioned randomly. If the number of trees of the same Gsp was greater than NClusHa Gsp , we randomly selected a tree of the same Gsp that was already positioned and used it as an anchor for a cluster. From this anchor position, we generated a random point around the anchor point at a distance equal to the modelled Dist Gsp . Finally, this position was then evaluated in the same way as in the previous point, i.e., by comparing the distances between the target tree and the 2 closest neighbours.

4.
The generator could start from an empty stand. However, a plantation scheme describing the spacing between planted trees and the presence of planting rows could be used to place planted trees species in these positions. A thinning path could also be added.

Performance Tests of the "Spatialiser"
In order to assess the effectiveness of the "spatialiser" to generate a realistic spatial stand structure, we selected the 15 control plots from the dataset, simulated the tree coordinates with the "spatialiser," and then calculated the Clark and Evans index for each Gsp.
We also simulated the same plots by assigning random coordinates to the trees (i.e., by not considering the interactions between trees). Hence, we simulated the 15 control plots in both methods ("spatialiser" vs. random) and calculated the Clark and Evans index for each Gsp and both simulations. These simulations were repeated 100 times.
Afterward, we performed an ANOVA to evaluate the difference between the index of the 2 simulations (random vs. "spatialiser"). As a final step, we calculated the bias between observed index and simulated index values.

Thinning Treatment Simulations
As TrT CT , TrT 1/3 , and TrT BF were used as independent variables in the models, the "spatialiser" could directly simulate a recently thinned stand with one of these 3 thinning treatments. However, the main use of the "spatialiser" (once it was integrated into the growth model) was to simulate a realistic spatial tree distribution in WS plantations where a thinning treatment could be applied. The applied treatment could be like those used in this study or of another type.
To illustrate this application, we chose 2 interventions. The first was the TrT CT thinning treatment already used in this study. The second was the creation of gaps. Gaps are carried out operationally in certain cutblocks in eastern Quebec to create openings for wildlife. The 2 selected treatments create openings in the stand that can promote regeneration by seed trees, as well as recreate other stand attributes such structural diversity. As skidding trails are an essential component of commercial thinning operations, they were added to both simulations (about 33 m apart and 4 m wide).
The first treatment (TrT CT ) aims to free the selected dominant trees from competition. This treatment is the first step of a series of interventions aimed at changing the silviculture path from even-aged to uneven-aged stands. It is characterized by two parameters: the number of selected crop trees (50-100/ha) and the thinning radius around each tree (1-5 m). We randomly selected the thinned trees among those of merchantable size within the stand, with the only requirement being that 2 selected trees had to be at least 9 m apart.
The second treatment (Gap creation) aims to create larger openings in the stand canopy in order to promote regeneration in the opening [54]. We performed this gap simulation to evaluate the behaviour of the spatial model for a treatment in which the choices of the number of gaps, their position in the stand, and their size are very important. In these simulations, the gap positions were chosen randomly. However, in an operational context, their choice should depend on the proximity of selected seed trees, the autecology of the species to be regenerated, the autecology of competing or harmful species, and the position of the skidding trails. The number of gaps varies from 1 to 4, and total gap area varies from 5% and 40% of the total stand area (i.e., 500-4000 m 2 ).
To simulate the effect of these 2 treatments, we selected 15 un-thinned stands (TrT 0 ) in which we upscaled plot size to 1 hectare and positioned all trees with the "spatialiser." We then applied both treatments and calculated the Clark and Evans index for each Gsp in these simulated plots. For each treatment, the simulation was repeated 100 times. We performed an ANOVA to compare the means of the different groups obtained after the simulations and to assess the impact of the treatment parameters on the spatial distribution (i.e., the number of selected crop trees and thinning radius around each tree for TrT CT and the number of gaps and total percentage of surface cut for gap creation).

Clark and Evans Index
The spatial distribution of WS tended to be regular (R WS = 1.34; Figure 3 and Table 3), which was expected, given that the stands were plantations. Conversely, BF and HW tended to form clusters, with HW following a stronger trend (R BF = 0.92 and R HW = 0.79; Figure 3). The final model form for predicting the aggregation index was found to be: where a 0 -a 6 are the parameter estimates for the fixed effects found in Table 3.  Table 2 for definitions of Gsp and treatment abbreviations. Table 3. Parameters estimations (standard error) of the model selected to predict aggregation index at the stand level (Equation (4)) (* (p < 0.05); ** (p < 0.001); *** (p < 0.0001)). Fit statistics show that the proportion of total variation explained by the model (R 2 ) was 0.84 and that the RMSE was 0.12. The RMSE and R 2 were 0.08 and 0.65 for R WS , 0.15 and 0.58 for R BF , and 0.14 and 0.64 for R HW , respectively.

Coefficient
For each species, the Clark and Evans index value was inversely proportional to the total stand density (a 2 = −0.00030; Table 3). At the species level, R WS was proportional to WS tree density (a 1 = 0.00022), whereas R BF and R HW were inversely proportional to BF and HW density (a 5 = −0.00030 for BF and −0.00031 for HW). On average and compared to the control plots, commercial thinning slightly decreased the aggregation index value, making a regular distribution more random and a random distribution more aggregated. The strongest effect was for TrT BF (a 4 = −0.91933), followed by TrT CT (a 4 = −0.47603), and TrT 1/3 (a 4 = −0.44839. However, for TrT BF , the aggregation index significantly increased with total density (a 6 = 0.00034).

Coefficient
Variable We observed that the minimum distance between the target trees and the first competitor increased with the aggregation index (b 1 = 0.60337; Table 4) and with the DBH of the target tree and its closest neighbour (T0: b 2 = 0.00018 and T1: b 4 = 0.00015). The distance was smaller for HW or BF, and the effect was greater for HW (T0: b 3 = −0.22557 and T1: b 5 = −0.27708). The greater the total density and the density of trees of the same Gsp as the target tree, the smaller the distance (NHa Gsp : b 6 = −0.00008 and NHa tot : b 7 = −0.00011).

Spatial Interactions within Species Groups
As R WS was superior to 1, we only modelled intraspecific interaction for BF and HW. The number of clusters per hectare and per species was predicted by Equation (7) (R 2 = 0.41, RMSE = 1.7): where c 0 , c 1 , c 8 , and c 11 are the parameter estimates for the fixed effects found in Table 5. Table 5. Parameters estimations (standard error) of the models selected to predict intra-species spatial distribution (Equations (7)- (9)). (* p < 0.05; ** p < 0.001; *** p < 0.0001). The number of clusters per hectare was greater for HW (c 1 = 0.42423; Table 5), with an average value of 30 clusters for BF and 57 for HW for all treatments. TrT 1/3 was the only thinning treatment that significantly reduced the number of clusters (c 8 = −0.31400) compared to the control treatment, but this was only for BF, since the interaction between HW and TrT 1/3 (c 11 = 0.32693) negated this effect for HW. The interaction HW with TrT BF (c 11 = −0.32710), for its part, indicated that there were more HW clusters in this treatment than in the control.

Coefficient
The shape parameters of the gamma distribution for minimal distance between trees of the same species within a cluster were 1.204 for BF and 0.959 for HW, while the scale parameters were 0.357 for BF and 0.403 for HW. For BF, the final model for predicting the minimal distance is presented in Equation (8) (R 2 = 0.19, RMSE = 3.22): where c 0 , c 2 , c 4 , c 5 , c 6 , and c 8 are the parameter estimates for the fixed effects found in Table 5. The final model is slightly different for HW (Equation (9), R 2 = 0.25, RMSE = 2.39): log(Dist HW ) = c 0 + c 3 × DBH T0 + c 4 × DBH T1 + c 5 × R Gsp + c 6 × NHa Gsp + c 7 × NHa tot + c 9 × DBH T0 × DBH T1 + c 10 ×NHa Gsp NHa tot + ε (9) where c 0 , c 3 , c 4 , c 5 , c 6 , c 7 , c 9 , and c 10 are the parameter estimates for the fixed effects found in Table 5. Inside a cluster, the distance between two trees of the same species increased as the DBH of the closest neighbour increased (c 4 = 0.00109 for BF and c 4 = 0.00215 for HW; Table 5). The DBH of the target tree also had a positive effect for HW (c 3 = 0.00287). For BF, the larger the difference in DBH between two neighbour trees, the larger the distance between the two trees (c 2 = 0.00113), whereas for HW, the higher the DBH T0 :DBH T1 ratio, the shorter the distance between the two trees (c 9 = −0.00002). The distance between two trees of the same Gsp was smaller when the density of the Gsp was high (c 6 = −0.00107 for BF and c 6 = −0.00279 for HW). For HW, this distance was also reduced in stands with a high density (c 7 = −0.00022).
For BF, the TrT 1/3 and TrT BF treatments increased the minimal distance between two trees compared to the control (c 8 = 0.12117 for TtT 1/3 and 0.16893 for TrT BF ; Table 5).

Performance Tests of the "Spatialiser"
Based on the Clark and Evans index calculated after the simulations, we observed that the "spatialiser" worked better than randomly assigning coordinates to trees to generate a realistic spatial stand distribution (Figure 4). We measured the correlation and bias between observed and simulated data. On average, for the "spatialiser," R 2 was 0.63 and RMSE was 0.19, whereas for the random simulations, R 2 was 0.01 and RMSE was 0.30. However, we observed disparities between Gsp (Figure 4). For WS, R 2 was 0.29 and RMSE was 0.18 with the "spatialiser," (R 2 = 0.01 and RMSE = 0.36 for the random simulations); for BF, R 2 was 0.28 and RMSE was 0.18 with the "spatialiser" (R 2 = 0.02 and RMSE = 0.31 for the random simulations). For HW, R 2 was 0.69 and RMSE was 0.21 with the "spatialiser" (R 2 = 0.10 and RMSE = 0.43 for the random simulations). Figure 5 presents the results of one simulation over a 1-hectare area before treatment and after applying either the TrT CT or gap creation treatment.

Thinning Treatment Simulations
With the TrT CT treatment, an increase of the thinning radius around each tree and the number of trees freed from competition significantly affected the aggregation index value. The link between the aggregation index and these two parameters was negative, except in WS, for which there were slight positive effects on the number of crop trees (R WS = 0.0001 for the number of crop trees and −0.0158 for thinning radius, R BF = −0.0004 for the number of crop trees and −0.0026 for thinning radius, and R HW = −0.0002 for the number of crop trees and −0.0136 for thinning radius). These negative parameters indicated that greater treatment intensity decreased the aggregation index value. However, the effect on spatial structure at the stand level was very limited ( Figure 6A).  Table 2 for definitions) and the indexes predicted from the "spatialiser" simulations and from random simulations.   Table 2 for definitions of Gsp and treatment abbreviations.
With gap creation treatment, increases in both the number of gaps and the total percentage of the area that was cut had significant effects on the aggregation index value. The link was positive with the number of gaps but negative with the percentage of the surface cut (R WS = 0.0001 for the number of crop trees and −0.6415 for the percentage of surface cut, R BF = 0.0044 for the number of crop trees and −0.3965 for the percentage of surface cut, and R HW = 0.0040 for the number of crop trees and −0.3049 for the percentage of surface cut). Despite its significance, the effect of the number of gaps was found to be very small at the stand scale. By comparison, increasing the cut surface had a much greater effect on stand spatial structure, especially for WS ( Figure 6B).

Modelling Spatial Stand Structure
A stand's spatial structure results from interactions between trees; it is one of the most important ecological characteristics influencing its development. In plantations, human actions such as (initial spacing and thinning) can be added and modify these processes. The model we developed describes local intra-and interspecific interactions and links point patterns with ecological/silvicultural processes. Despite the absence of precise biological information such as crown characteristics, which are known to strongly influence the behaviour between two closely located trees, it was possible to develop a "spatialiser" that can reproduce attractive behaviours among trees of the same species and repulsive behaviours among trees that are too close. Choi et al. [55] showed that adding more detailed competition variables provided only slightly better predictions of spatial patterns for hardwoods. In addition, in most situations where data are limited to standard inventories (e.g., species, DBH, and height), the "spatialiser" would need to rely on allometric relationships to predict the needed crown attributes [56]. This would generate a source of uncertainty and error propagation. Hence, coupling a "spatialiser" to standard inventory data could be an interesting alternative to measuring the location of trees in the field and would allow for the use of models working with spatially explicit plantation growth simulators. The spatial information could also be used to calculate structural habitat characteristics and shed light on the stand community structure and relative abundance of animal species [57].
Some limitations and simplifications of this spatial simulator should be considered. The models were calibrated on plantations that have less complex structures than natural stands, even if plantations in eastern Quebec have a high proportion of regeneration by other species [58]. In addition, we chose to group all the different hardwood species together and to remove some species of coniferous very rarely present as well as noncommercial species from the analysis. These species, however, could have a significant effect on competition or be of great interest to understand the diversity of species. Additionally, saplings and merchant trees were studied together. However, saplings are generally present in large numbers, whereas merchant trees are the result of competition over a longer term. These differences could be reflected in the spatial distribution of these two groups (Table 1). However, despite these limitations, the combination of tree-scale models considering DBH and species of neighbouring trees, as well as a plot-scale model describing the type of aggregation of a species, seems flexible enough to adapt to different types of stands.

Spatial Interactions between Individual Trees
The construction of the "spatialiser" using several models has the advantage of allowing for a direct interpretation of inter-and intraspecific interactions. At the stand level, the Clark and Evans aggregation index (R) [47] showed that in the studied plantations, WS had a regular structure while BF and HW tended to aggregate. For the planted WS seedlings, this was expected. However, for the WS saplings that were naturally established later on, the structure was more aggregated (Table 1). The Clark and Evans index that was chosen here is a simple and easy to measure variable for analysing the spatial behaviour of a species at the plot scale. The simulator models were of the tree scale. The unique value of the index then allowed us to simply assess the results of the simulation at the plot level. In order to apply the simulation approach to other types of stands with more species involved or used as surrogate measures to quantify biodiversity, some other indexes could be used to describe or validate the stand spatial structure [59].
The grouping of hardwood species in our analysis prevented us from examining the behaviour of each individual hardwood species. Some of these, such as paper birch, tended to aggregate. In the case of aspen, vegetative reproduction by suckering could explain its tendency to group [10]. Other species, such as red maple and sugar maple, appear to be less aggregated, although red maple can form clusters in cases where it forms stump sprouts [60]. Both BF and HW tended to aggregate; however, the observed clusters were different. For BF, clusters were fewer, but they tended to include both merchantable-sized trees and saplings. For HW, clusters were smaller, more numerous, and consisted mainly of saplings.
Currently, in the Bas-Saint-Laurent region, natural stands are mainly dominated by BF and HW, to the detriment of WS [58]. In plantations, BF and HW can represent up to 25% of the trees. In order to limit competition and to favour the growth of WS, a species that is becoming rarefied [34], good knowledge regarding this species' spatial behaviour and natural regeneration, as well as the development of competing species, is needed. The tendency to BF and HW to cluster can help managers choose approaches that limit their spread to the detriment of WS when they prescribe clearing, precommercial thinning, or commercial thinning treatments.

Thinning Treatment Simulations
At the tree level, species and DBH were the most important parameters influencing the distribution of BF and HW. At the stand scale, however, thinning caused only slight changes in the aggregation index, except for TrT BF . This was normal, because the purpose of these treatments was to reduce competition in order to increase mean tree wood production. In addition, TrT CT could be considered as a first step to convert plantation structure from regular/even-aged to irregular/uneven-aged if the objective is to reduce the differences between natural and managed forests [61]. However, after five years, this treatment had no significant effect on diameter distribution or spatial structure. Dupont-Leduc [61] suggested that a five-year interval is probably too short to observe such differences. The lack of structural change may also suggest that the thinning intensity of the TrT CT treatment was not sufficient or too evenly applied [61]. Indeed, for the two other thinning treatments, cuts were only done around 50-100 dominant trees with a clearance limited between 2.5 and 3 m (TrT CT ) or mainly with the aim of promoting dominated trees (TrT 1/3 ). Nevertheless, TrT CT could also allow for the long-term development of seed trees and thus the maintenance of WS rather than its replacement by BF. Though it is unclear whether the small observed differences in stand structure compared to the control will be conserved or will evolve in time, the reduction of the number of trees per hectare and of competition could help trees to grow faster and create more heterogeneous forests. Gauthier et al. [62] showed that with a TrT 1/3 thinning, stand structural heterogeneity increased 10 years after the harvest as a result of the better growth of bigger trees and retention of saplings. However, according to Schütz [29], many entries are necessary for the conversion effect to become apparent from regular/even-aged stands to irregular/uneven-aged stands. Thus, several years of stand monitoring will be required to evaluate the real long-term impact of these anthropic disturbances [63].
The "spatialiser" can be used to generate control stands with localised trees where "silvicultural" treatments can be simulated. This would allow for the evaluation of the chosen approach's effectiveness and the testing of silvicultural treatments [64]. Such an exercise suggests that, for example, TrT CT could be applied with a larger clearance distance around targeted trees, since simulations showed that the current TrT CT treatment had a fairly limited impact and that a rather large clearance around trees was needed to modify the spatial distribution at the stand scale. Hence, to emulate natural disturbances, the treatment should create larger gaps. According to the simulations, the creation of larger gaps (e.g., cut area: 5-40% of the total) could be another effective approach to modify stand spatial structure. However, for both these examples, the thinning effect was observed directly after application. An evaluation of the long-term effectiveness of the chosen approach is also required. To do this, the "spatialiser" could be harnessed for a growth model and help forecast evolution at both the individual tree and stand scales.

Conclusions
This study shows that the spatial structure of WS plantations in eastern Quebec, Canada can be modelled and simulated. The diameter and species of the neighbour trees were the most important factors that explained the observed distances between trees. A point process model proved effective at representing inter-tree interactions through repulsive and attractive processes. Analysing all hardwood species in a single group limits the understanding of each species' spatial behaviour. However, since BF and hardwoods together could represent 30-40% of a stand's basal area, the models allow for the better understanding of their interactions with WS.
Our results showed that thinning treatments can have slightly significant and shortterm effects on spatial structure. If the objective is to convert to another stand structure type, linking this "spatialiser" to growth models would allow for the simulation of treatments over larger areas and longer periods. It would also allow for parameters to be modified in order to evaluate the effectiveness of the chosen approach. Thus, silvicultural treatments could be adjusted to better reach management objectives.