Modeling the Effect of Environmental and Topographic Variables Affecting the Height Increment of Norway Spruce Stands in Mountainous Conditions with the Use of LiDAR Data

Differing levels of humidity, sunlight exposure or temperature in different areas of mountain ranges are fundamental to the existence of particular vegetation types. A better understanding of even local variability of trees may bring significant benefits, not only economic, but most of all, nature-related. The main focus of this study was the analysis of relationships between increment in stand height, age and the natural topography in the examined area. Among others, the following were examined with regard to their influence on the growing process: age, altitude above sea level (m a.s.l.), aspect and slope, topographic wetness index (TWI), and topographic position index (TPI) generated from an airborne laser scanning (ALS)-derived elevation model. To precisely calculate forest growth dynamics in mountain conditions for different spruce stands, repeated airborne lidar measurements from 2007 and 2012 were used (with resolution respectively 4 and 6 pts./m2). Detailed information on every stand including species composition, share of individual species, as well as their age, were acquired from the State Forests IT System (SILP). It was proven in this study, that environmental and topographic variables may have an impact on forest growth dynamics on even closely located areas. Apart from the age, the greatest influence on tree growth has an altitude above sea level, aspect and slope. The highest height increment of spruce was observed in the stands of up to 30 years old, those that had grown at an altitude under 850 m a.s.l., on the slopes up to 15 degrees or on those which were on the northeastern exposure. The results obtained show that the physiology of species, even those that are well known, largely depends on local topographic conditions. The proven impact of different topography factors on the growth of spruce may be used while planning economic activities in precision forestry. Additional research with using multiple laser scanning in the context of other regions or other species may bring us better recognition of local growth conditions and in consequence, significantly better planning and higher revenues obtained from the sale of trees.


Introduction
The Western Sudetes is an area located on the border of the Czech Republic, Poland and Germany.As the westernmost range of the Sudetes, they are the first orographic barrier for air masses from the west and southwest.Such a location makes the climatic conditions here more unfavorable in comparison to surrounding areas.In the 1980s, together with air pollution and insects, climatic conditions were one of the triggers of a large-scale deforestation which have stricken spruce stands.Even then, under the extremely hard conditions for forest growth in the whole area of the Sudetes, differences in the stress reaction of forest communities were observed and in consequence, differences in their growth.Such a variability may also be a result of human activity in XVIII and XIX Century, when broadleaf species were used in glass manufactures and replaced by spruces seedlings from different provenances.We know that tree height increment may be subjected to climatic and topographical factors but so far there has been no precise technique to help us to answer the question: what is the precise influence relating to each factor?
Airborne laser scanning has become the standard for acquiring information about the environment for the last two decades.Nearly the total area of Europe has been inventoried in such a manner, and now in some areas, Airborne Laser Scanning (ALS) measurements have been repeated.It should be expected that the forthcoming years will bring a significant decrease in the costs of using this technology, and it will become (in some part of the globe already is) an elementary source of information also in other areas of the world.In this context, an emerging question is whether and how much information regarding the height increment of trees may be delivered through repeated laser scanning.The question is even more significant for areas where land topography is more diverse, with different aspects, slopes, altitudes above sea level, etc., because land relief may have huge impacts on spruces growth [1,2].
Dynamics of forest communities expressed by the height increment are subject to environmental gradients.The tree height increment is potentially subject to temperature and wetness gradients and this relationship may give us a measurable reference for understanding patterns of resource use, especially the physiological mechanism of how trees may react to climate stress.These patterns are the factors that predispose these ecosystems to climate-affected disturbances, the knowledge of which is crucial for sustainable forest management [3].
The dependence between topography and various environmental factors has long been observed [4][5][6].The most frequently used topographic indices are aspect, slope, slope position, elevation, longitude and latitude [1,[6][7][8][9].On-site evaluation of latitude and land relief are used as indirect measures of local and regional climate [1,[10][11][12].Elliott et al. [13] proved in their study that elevation, terrain shape and soil organic matter, were the variables most intensely related to vegetation distribution.Fekedulegn et al. [14] found that all species showed significant differences in growth between northeast and southwest aspects, one exception was the northern red oak.
Influence of topography were confirmed also in studies on spruce forests.The strong correlation between land relief and the growth of spruces in the mountains of Central Europe was shown by Albert and Schmidt [15] and Socha [1].Topography, together with winds and fog, may also have a negative impact on the growth of trees, which was confirmed by Farrelly et al. [16], Watt et al. [17], Błaś [18] and Godek et al. [19].Respectively, Night et al. [20] proved, based on their study of spruces at various levels of humidity in Canada, that there is no difference in height growth between biogeoclimatic subzones (different regions of Canada) of boreal forests up to 75 years in British Columbia.At the same time, growth below breast height after age 75 in drier and cooler climates is faster comparing to growth in wetter and warmer boreal climates.It should be noted that the topographic wetness index (TWI), may be used as an indirect measure of the water status of a given site [21].Although it's known that numerous inter-correlated ecological factors (e.g., precipitation, winds, nutrients and sunshine duration) may have an influence on the height increment of trees.Junttila [22], together with Messaoud and Chen [23], concluded that the most important impact on spruce growth is temperature, or more precisely, temperature has an impact on young spruce shoots and as a result, on their growth.
Due to the difficulties associated with the measurement of the height increment in standing trees, the analysis of the height increment and its response to environmental gradients have so far received less attention and have been mainly limited to relatively small samples collected from cut trees [24][25][26][27].
During the last years, plenty of remote sensing techniques appeared, which are allowed to increase the accuracy of measurements and decrease the costs in the same [28].Light detection and ranging (LiDAR) systems are commonly and successfully used for forest inventory and management purposes [29,30].Repeated ALS observations enable measurements of tree height increment over time [31][32][33].ALS measurements may be biased as a result of several factors like species architecture, flying height, instrument specifications and the measurement method used.These sources of errors appear as measurement inaccuracy that can range from a few centimeters up to even a few meters [31,32,34].At the same time, Hyyppä et al. [35] proved that repeated ALS measurements may be used to measure forest growth and standard error (based on estimation of height increment at the stand level) was less than 5 cm.Airborne LiDAR data, acquired in 1998 and 2003, were also used to detect tree growth in Québec, Canada by St-Onge and Vepakomma [36].In 2008, Hopkinson et al. [37] used ALS data acquired four times on the pine plantation to evaluate the estimation of plot-level mean tree height increment.They proved that the 100th height percentile provides the most solid estimation of field-measured forest growth.In another study, Yu et al. [31,38] provided an estimation of Scots pine (Pinus sylvestris L.) and Norway spruce (Picea abies (L.) Karst.)canopy growth based on two ALS acquisitions (the second one, acquired 21 months later) and proved that precision of plot-level growth was of 10 to 15 cm.The best method to estimate the mean tree height increment was the maximum tree height differencing method.Determining the specific increment specification from ALS data is critically dependent upon this low error because marked change in stand height can only be measured if the height increase is greater than any biases in the LiDAR measurement [29].
The objective of this study on Norway spruce, representing all the elevation and aspect gradients in the Western Sudetes, was to: 1.
detect the site-dependent variability of mean annual height increments, 2.
explore the dominant environmental factors limiting height growth using repeated ALS data obtained for the whole research area.
We assumed that the limiting impact of temperature and water stress on height increment is affected by height above sea level and topography, also, we assumed that ALS data from two period coverage of vast areas should give unbiased height increment trends and increment estimates.Therefore, we revealed that elevation, aspect and topographic indices, through their influence on thermal and moisture conditions, may affect the height increment of Norway spruce.The determined factors limiting height increment may be used as indicators of susceptibility of a given site to disturbances.
In our study we independently applied two different approaches to find the most important topographic factors that influence the growth of the spruce.In the first approach, we conducted classical statistical method multiple linear regression (MLR).In the second one, machine learning technique boosted regression trees (BRTs) [39] was used.

Study Area
The research area included the Western Sudetes (the Jizera Mountains and a part of the Karkonosze Mountains), located in the south-west of Poland, next to the border with Germany and the Czech Republic (Figure 1).The area covers approximately 197.5 km 2 , most of which are Polish State Forests located in the Regional Department of Wrocław, the Świeradów Forest District and the Szklarska Poręba Forest District.Administratively, the area is located in the Dolnośl ąskie Province and the districts of: luba ński, lwówecki, and jeleniogórski.
The Western Sudetes are the westernmost mountain range in Poland, and they constitute one of the orographic barriers for air masses flowing from the southwest and west.Such topography definitely impacts the climate conditions, which in many aspects may be considered as extremely unfavorable.The climate is cold and humid.Along with the increase of altitude above sea level, precipitation increases and the air temperature decreases significantly, it ranges from 600 mm in mountain basins and in the foothills to 1200 mm and above on mountain peaks.The average annual temperature is 6-8 °C in lower parts, and at upper elevations it decreases below 3 °C.Due to such conditions, the growing season is fairly long i.e., the maximum of 210 days in a year it gets shorter as the altitude increases.The vegetation in the Sudetes is characterized by storied arrangement, conditioned by climate differences, where the following vegetation zones may be specified: foothills: <500 m a.s.l., lower subalpine forest: 500-1000 m a.s.l., upper subalpine forest: 1000-1250 m a.s.l., dwarf pine zone: 1250-1450 m a.s.l., alpine zone: >1450 m a.s.l.(in the Śnie żka massif-excluded from the study).

Forest Digital Map Data and Information System
For the work, data from the State Forests IT System (SILP) was to provide detailed information on every stand including species composition, share of individual species, as well as their age.SILP includes precise information about species' structure and the age of each group of trees thanks to yearly updated information by the foresters, verified additionally during the work on the ten year forest management plans.Detailed forest stand borders were specified using data from the Forest Digital Map, which is updated annually, the same as SILP Database by the State Forest's Service employees.
Information on the species composition obtained from the SILP DB, was used to select these forest stands, in which the share of spruce in the upper story ranged from 80% and upward.

LiDAR Data
ALS data were collected in 2007 and 2012 in vegetation periods using an Airborne Laser Terrain Mapper 3100 Optech Laser Scanner and RIEGL Airborne Scanner Lite Mapper 6800i, IGI system.Despite the differences in the type of laser scanner and average point density, similar average errors of altitude and grid coordinates were obtained (Table 1).

Processing the Digital Terrain Model (DTM) and Canopy Height Model (CHM)
The developed digital terrain model (DTM) which is based on data from the ALS was of a spatial resolution equal to 0.5 m.Stere ńczak et al. [40] tested the accuracy of the DTM and proved that the root mean square error (RMSE) ranges from 0.19 m to 0.23 m depending of filtration and interpolation methods.In this study, we used DTM generated using an algorithm based on an active Triangulated Irregular Network model (TerraScan Software) and an Inverse Distance Weighted interpolation method which is implemented in ArcGIS 10.3 software (RMSE equal 0.19 m).The DTM was an entry layer (basic information layer for further analysis) used for generating slope raster.For statistical analysis, the aspect values were transformed using the following equation [41]: K-is the transformed value of aspect, varying from 0 to 2, A-is azimuth (degrees), A max -is the assigned direction (azimuth) of maximum importance equal 45 • .
To generate raster information, DTM (based on ALS data), which had been resampled to a 100 m resolution, was used to verify the dependence of local and regional variability on the total mountain chains' exposure.The use of an original DTM (based on ALS data of high resolution) would only provide the necessary information on the local layer, preventing the analysis of general trends for the whole separations, whole slopes, or their greater fragments.Based on the DTM, the Topographic Wetness Index was also calculated in the SAGA GIS 2.3.2 software.A modified TWI algorithm was implemented in SAGA GIS 2.3.2 and it is based on a modified catchment area calculation.As a result, it predicts a more realistic, higher potential soil moisture for cells situated in valley floors with a small vertical distance to a channel compared to the standard TWI calculation [42].The topographic positioning index (TPI) was used to describe the landform.In this study, the TPI was calculated using a circular neighborhood with a 200 m radius.TPI was calculated using the following equation [43]: where: Z 0 -elevation of the model point under evaluation, Z n -elevation of the grid, n-the total number of surrounding points employed in the evaluation.
The application of ALS data from two different periods (2007 and 2012) in this study enabled a very accurate height increment estimation of the spruce stands in the Western Sudetes.Based on ALS data, digital surface models were produced.Then, using the raster difference (DSM-DTM) in map algebra, Canopy Height Models (CHM) were established.It is essential to note that to create the CHM, a single DTM from 2012 was used for both years.ALS data, applied in this work, as it was also proven in separate research [44], are characterized by high accuracy compared to ground measurements with a hypsometer or measurements based on stereophotography.The application in this study of ALS data from two different periods (2007 and 2012) enabled a highly accurate specification of the height increment of spruce stands in the Western Sudetes Jizera Mountains.In the study period, the error in determining altitude in this manner (±15 cm) was much lower than the increment of stand height, even stands older than 102 years.
The process of adjusting the spatial extent of the analysis area was carried out by performing the calculation of the top height, which was defined as the mean height of the 100 highest trees per hectare.Therefore, for each plot the top height was calculated as the mean height from 25 subplots (10 × 10 m).In the first step, raster data (i.e., CHM) was cut to the geometric range of the forest area.The borders of forest sub-compartments from the Forest Digital Map are characterized by a few-meter shift comparative to the actual course.Therefore, to minimize the effect of geometry on the results of the studies, the authors decided to limit the analyzed area by generating the sub-compartments from the borders.In addition, buffers were generated for the public and forest roads.The size of each buffer was equal to 10 m.Within the borders of the accordingly specified research area, the basic plot grids were generated by the following dimensions: 10 × 10 m and 50 × 50 m.The coordinates of the upper left corner of the basic plot grid for both dimensions were the same.As a result, in every single 50 × 50 m basic plot, there were 25 subplots [45].In forest stands mixed at least with 80% spruce, 1939 grids of 50 × 50 m were located and used for further analysis.The total study area was 484.75 ha.
The top height in individual years was calculated applying the following methodology: every 10 × 10 m basic plot the highest value was assigned based on the CHM raster (2007 and 2012).Furthermore, based on a 10 × 10 m basic plot, the average height of a forest stand was estimated for every 50 × 50 m plot.The increment of tree height in 2012 for every 50 × 50 m basic field plot (Figure 2) was calculated based on the difference in CHM values between the years 2007 and 2012.

Statistical Analyses
Different parametric and nonparametric approaches were used in certain forest applications with remotely sensed data.Parametric modelling techniques such as classical linear regression have widely been used due to their familiarity and practicality [46], however, they require specific conditions which may be violated.Nonparametric approaches have been studied extensively lately, such as nearest neighbor methods [47,48], regression trees [49], and neural networks [50] and boosted regression trees [51].Among different machine learning techniques, boosted regression trees (BRT) were the most effective to find the most important topographic factors that influence the growth of the spruce, so this method was used in our study.We also conducted multiple linear regression (MLR) to compare the results of both methods [46].
The advantage of BRTs is that there is no need to drop prior predictors.They can fit complex nonlinear relationships, they automatically account for interaction effects between predictors, and they can handle different types of response data [39,46].The general idea of the BRT statistical technique is to compute a sequence of simple trees, where each successive tree is built for the prediction residuals of the preceding tree.BRT models are provided on the basis of how good each model explains the observed data (training data correlation) and how well excluded data can be predicted (CV correlation).Relative importance values are generated for all predictors (totaling up to 100% for all included variables) to indicate each variable's relative contribution to the model [39,52].It is based on how often the variable is selected and its ability to improve the model.The fitted function visualization in a BRT model is reached using partial dependence functions that demonstrate the effect of a variable on the after accounting for the average effect of all other variables in the model [39].
The increment of tree height between the years 2007 and 2012 was used as a dependent variable (height increment).The set of predictors consisted of topographic factors-aspect, slope, altitude, TPI, TWI (MLR1 and BRT1).Additionally, models considering all the previous mentioned variables and stand age as an independent factor were analyzed (MLR2 and BRT2).
We used BRT models which were implemented in R version 3.0.3(R Development Core Team, 2015) with the gbm package [53,54] with a Gaussian response distribution and 10-fold cross-validation procedures.Each model was defined with a tree complexity (number of nodes) of 3, a learning rate of 0.001 and 0.75 bag fraction (proportion of data selected at each step).

Results of Multiple Linear Regression
Increments of the tree height during a five-year period were equal to 1.94 ± 0.63 m, the minimum value was approximately 0.64 m and the maximum value was almost 4.5 m.It means that the average growth per year was equal to 0.39 m.For half of the plots, the values of height increment were from 1.47 to 2.34 m.Table 2 indicates the relationships between response and explanatory variables for MLR.The method turned out to be not effective in studying the relationship between topographic variables and height increment.The multiple correlation coefficient between response and topography for MLR1 was very low (R MLR1 = 0.15).In the MLR2 (with additional age factor) correlation between response and explanatory variables was high (R MLR2 = 0.80) but was determined almost by stand age (R age = 0.78) and the role of the other factors was only minor.The reason is that the linear regression method indicates only linear relations and cannot handle interaction effects between predictors.Correlation analysis did not show clear relations between topography and response (Table 3).Topographic position index (TPI), topographic wetness index (TWI).

Boosted Regression Tree (BRT) Results
The results of BRT analysis confirmed that the response depended mostly on stand age (with the highest relative importance percentage approximately 76.06% of the BRT2 model) and analysis of height increment only for the group of topographic variables are not sufficient (Table 4).BRT1 results are not stable, there is a disproportion between training and tested data correlations.Altitude was the most important topographic factor but did not show a clear trend (Figure 3).
Based on the analysis of partial residuals of the BRT model, it could be stated that the height growth of stands growing on steeper slope areas is larger than in stands growing in the plains.With an increase in slope from 0 to 15, the height increment increases.After reaching the value of 15, further increases in slope do not result in an increase in height increment.Aspect was not a key predictor of height increment, but it is worth noticing that response values observed in stands growing on N, NE, E aspects were higher than on S, SW, W aspects (Figure 3).The role of TPI and TWI was relatively small.In the case of TPI, the BRT analysis showed a negative correlation with height increment.A positive association between the topographic wetness index and the growth of spruces up to 5 was observed, after which the lower values of height increment was observed.The pairwise interaction strength of explanatory features for a BRT analysis of height increment was also checked.The most important one was interaction age and altitude (Figure 4).In our dataset we have plots over 100 years with altitude values between 500 and 800 m a.s.l.As a result, lower height increments were noticed for these altitudes.

Discussion
The results of BRT analysis confirmed that, apart from the age, the greatest influence on tree growth has an altitude above sea level.However, as shown by this research, the increase is not linear, and it is much more intensive up to the altitude of 850 m a.s.l.than above this level.It is 50 m higher [1] than in case of the Beskids (800 m a.s.l.), where the mountains are characterized by a slightly milder climate.The result is astonishing taking into account the fact that in the West Beskids forest strata are lower than in the Beskids.At the same time, Zawiła-Niedźwiecki [55] has shown that during the ecological disaster of the 1980 s, the best-preserved stands were the ones below 600 m a.s.l., therefore there is probably no significant relationship between the intensive growth of spruces up to the level of 850 m a.s.l. and their resistance to atmospheric pollution in this case.
The analyses performed confirm the slightly weaker growth of spruces on the southwest exposures in relation to spruces on northeast exposures, which matches the results of Socha [1] for the Beskids.Southwestern slopes are drier and sunnier than northeastern slopes because of the more intense solar radiation.
The faster growth of spruces was also observed in the areas characterized by greater slopes; however, it should be noted here, that flat areas mainly occur in the upper parts of mountains from 800 m a.s.l. up to 1100 m a.s.l.In order to determine errors resulting from the lower accuracy of ALS on slopes [56,57] and limit their impact on the created models, most terrain plots were located there.
TWI used in the research indicate the potential of a given site for saturation.A high value TWI indicate on high and low on lower potential for saturation.We found that in the case of the analyzed area, the role of TPI and TWI on Norway spruce height growth was relatively small.In the case of TPI, the BRT analysis showed a negative correlation with height increment.A positive association between the topographic wetness index and the growth of spruces up to 5 was observed, after which, the lower values of height increment were observed.
In the Sudetes as well as in the Beskids, west and southwest winds are predominant and have been identified as responsible for the negative impact on the tree growth process [18].This is because wind direction is responsible for the greater deposition of pollutants and the occurrence of snow rime.Based on research findings of research using tree ring data, Godek et al. [19] indicated meaningful differences in the growth of 80-year-old spruce trees depending on the exposure and the altitude above sea level.Areas more exposed to pollution were characterized by small annual increments, which, however, is not confirmed in the current research based on data from the years 2007-2012.Due to the data from the pollution monitoring station in Jelenia Góra we know that nitro oxides as well as sulfur oxides previously emitted by the heavy industry located in the Black Triangle (border region shared by Germany, Poland and the Czech Republic, characterized in the past by extremely high levels of pollution) no longer playing any significant role.At present we can say, there is no apparent dependence on the direction of pollution and the height increment of spruce trees.The current observed tendencies match those previously described by Gieruszy ński [58], who observed the highest site indices on the northern and the lowest site indices on southern aspects.It should be noted additionally that the influence of the origin of spruces on their health condition is also unknown.It is not certain whether the trees currently growing in the Sudetes originate from this region or not and what kind of consequences.
Although the relationship between the increment and aspect of trees has been known for a long time [59][60][61][62][63], various species respond to aspects to a different extent and, sometimes, in a different manner [14].It cannot be excluded that, apart from aspect and therefore the resulting solar radiation and humidity, other factors also have a significant impact on differential tree increments [6].As shown by Lee and Sypolt [64] in their study conducted in West Virginia, in certain years, the soil humidity on opposite aspects did not differ significantly.
Certainly, some impact on the results have comes from measurements technique and topography itself.Although ALS data enabled the accurate estimation of the top height of tree stands, it is necessary to remember, that for trees and tree stands that were 10 years old in the year 2012, the applied data could have been insufficiently precise.Perin et al. [65], based on research in Belgium, claims that a modeling site index for spruces younger than 20 years old is invalid due to very high variability, resulting from weather conditions, damage by animals, etc. Undoubtedly, the use of ALS data contributes to the quality of height increment and site index measurements [66,67], particularly if data from two periods are used.
Although, for flat areas or sites that are homogeneous with regards to the geological composition in site growth and site quality modeling, the topography may be omitted [68], for areas of more varied composition it has a more significant role [61,69].At the same time, the research conducted confirmed the great dependence between the growth of spruces and elevation.This factor has already been identified as significant in spruce productivity in rare studies in the USA [70], Canada [71], Sweden [7], Poland [72,73] and France [9].This does not mean that all species demonstrate an equally strong relationship between growth and elevation.For example, for the Engelmann spruce in Canada, the correlation between the growth and the elevation was only −0.29 [74].
The results obtained show that the physiology of species, even those that are well known, largely depends on local topographic conditions.The severe impact of topography on the growth of spruce may be used while planning economic activities in precision forestry.The use of detailed information regarding the growth of millions of trees provides an entirely new dimension to forest research.
Boosted regression trees proved to be an effective method for modelling the influence of environmental and topographic variables on the growth of Norway spruce.BRT generated better results than classical linear regression approaches [39], implicitly accounts for interactions among predictors, and allows for non-monotonous correlations between the response and the explanatory variables [39].BRT has been previously used in ecological studies to predict species distributions [75], disease risks [76] and risks of bark beetle-caused tree mortality [51,52] but, to the best of our knowledge, this is the first study to use it for modelling the growth of trees, especially Norway spruce.
Additional research is needed to use multiple laser scanning not only in the context of this single species, but also with all other species, as well as with regard to their spatial variability.A better recognition of local growth conditions of individual tree species will enable significantly better planning of efforts to increase growth, and, as a consequence, it will enable higher revenues to be obtained from the sale of trees.

Conclusions
As it was proven in this study that even closely located areas may be characterized by different growth dynamics depending on their local topography.Apart from the age, the greatest influence on

Figure 1 .
Figure 1.Research area including two forest districts of the Western Sudetes against the map of Poland.

Figure 2 .
Figure 2. Assumptions of basic fields localization, which were used for further analysis of the top height of forests stands based on the digital terrain model (DTM).

Figure 3 .
Figure 3. Partial dependency plots for the BRT model exhibiting the relationship between age (a) and topographic factors (b-f) and the height increment.Percentages represent the relative contribution of the predictors to the BRT model.The fitted function (y-axis) shows the effect of the predictors on the height increment after all other predictors are averaged (higher fitted functions depict greater height increment).

Figure 4 .
Figure 4. Interaction effect of age and altitude for a BRT analysis of height increment.

Table 1 .
Airborne Laser Scanning data acquisition specifications and characteristics.

Table 2 .
Multiple linear regression (MLR) between height increment and topographic factors without age (MLR1) or with age (MLR2).

Table 3 .
Pearson correlation coefficients between topographic variables and height increment.

Table 4 .
Boosted regression tree predictor importance values, related to the height increment.