^{3}
ha
^{-1}
yr
^{-1}
, which corresponds to the average growth of the
^{2}
(85%) of the wood volume model for 2011 is similar to that obtained in other studies. However, as expected, the one obtained for the wood volume model for 2015 (80%) is significantly lower.

^{2}ha

^{-1})

^{3}ha

^{-1})

Over recent years, there has been increasing interest in the evaluation of the forestry variables, such as biomass and timber volume, especially within the context of the Kioto protocol (

The link between the parameters obtained in forestry inventories and LiDAR (

Now that LiDAR technology has demonstrated its efficiency and accuracy in estimating most forestry variables, today the challenge is to try to take advantage of this tool in order to predict the annual growth of the plant structure at a large scale. The ease with which data can be updated allows us to readjust the estimations and quickly assess changes or disturbances in forest stands by estimating their development over time, or establish comparisons between different areas with a simple update of the forestry cartography (

This study assesses the application of data from the last LiDAR flight made over the Autonomous Community of the Basque Country in 2012 to evaluate the growth of volume over bark (Vob) during 5 years on stands of

The area based approach (ABA) for forest inventory applications using LiDAR has been utilised in this study, following the methodology applied by other researchers (

The study area lies in the Province of Bizkaia in the Autonomous Community of the Basque Country (ACBC) in the north of Spain. Its center is the municipality of Muxika and data from the surrounding municipalities between the coordinates 2º 30’ – 2º 54’ west and 43º 06’ – 43º 24’ north are also used. The average height is 230 m. The mean annual rainfall is 1277 mm and the mean annual temperature is 13.8 ºC. Frosts are infrequent and no physiological drought is apparent. With regard to the silviculture applied in the area, the rotation cycle is, on average, 35 yrs, varying between 30 and 40 according to the site quality and the price of wood at the time when harvest takes place. Clear cutting are applied followed by repopulation with initial densities of 1300-1600 trees ha
^{-1}
. Between the 7th and 9th year, non-commercial thinning takes place, eliminating 500 trees ha
^{-1}
and pruning up to 2.5 m on the trees that have not been eliminated. Between the 12th and 15th year, the first commercial thinning is carried out, leaving between 500 and 700 trees ha
^{-1}
. Between the 17th and 20th year, the second commercial thinning is applied, leaving between 400 and 500 trees ha
^{-1}
. The third thinning is carried out between the 22nd and 25th year, leaving between 250 and 300 trees for the clear cutting. The total production is on average 550 m
^{3}
ha
^{-1}
in one rotation, and therefore the average growth is 14.85 m
^{3}
ha
^{-1}
yr
^{-1}
. Volumes and growths of 3153 ha occupied by

Two sets of data were used according to the year in which these were gathered. In order to fit the wood volume model from 2011, data was used from the Fourth National Forestry Inventory (IFN4) of the Province of Bizkaia (
^{2}
ha
^{-1}
were also eliminated as these were not considered possible for this species in this area. Following the application of these filters, only 8 additional plots were discarded. Therefore, 90 plots were used finally in this study. Their location can be seen in

IFN4 plots are of variable radius, so that a tree will or will not be measured in accordance with its diameter and the distance to the center of the plot. For each selected plot, the number of trees per hectare (N, trees ha
^{-1}
), the dominant height (Ho, m) defined as the mean height of the 100 thickest trees per hectare, the mean height (Hm, m), the mean diameter (Dm, cm), the basal area (G, m
^{2}
ha
^{-1}
) and the wood volume over bark (Vob, m
^{3}
ha
^{-1}
) were calculated. For the last of these calculations, the Excel add-in CUBIFOR (

In order to be able to assess the possibility of using LiDAR data from a single flight when estimating the forestry resources in subsequent years, 55 field plots were measured during the month of December, 2015. So as to prevent any lack of accuracy in locating the center of the plot from affecting the estimation of the final results, those areas with a sufficient level of uniformity were selected so that a displacement of up to 10 m in the center of the plot would have only a slight effect on the final result.

With the aim of locating uniform areas, two influential variables were selected, the first being Canopy Cover Fraction (CCF), as a percentage of the surface covered by the canopy of trees, and the second being the height of the vegetation. The data gathering area was divided with a 10 m-sided grid and in each 100 m
^{2}
square and the CCF was estimated as a percentage of first returns, higher than 3 m over the total of first returns of the LiDAR points contained in it. For each grid, the typical deviation and the average of its CCF value and those of the 8 surrounding squares were calculated and those squares in which the coefficient of variation (CV) was not greater than 30% were selected if the CCF was lower than 70% and those in which the CV was greater than 10% if the CCF was greater than 70%. To estimate the height of the vegetation, the data gathering area was divided with a 10 m-sided grid and in each grid cell, the height of the 95th percentile of the heights of the LiDAR points contained in it was calculated. For each grid cell, the mean and typical deviation of its value and that of the 8 surrounding cells was calculated. Those grid cells in which the CV was not greater than 20% were selected if the 95th percentile was less than 8 m; those grid cells in which the CV was not greater than 10% were selected if the 95th percentile was between 8 and 12 m, and those in which the CV was not greater than 5% were selected if the 95th percentile was greater than 12 m. Only those grid cells that met both conditions were candidates for containing a plot center for data gathering. In accordance with these conditions, 55 points were found and located in the field with a GPS navigator. From each point, the dominant height was measured as the average of the thickest three trees within a radius of 9.77 m from the point that corresponds to a plot with a surface area of 300 m
^{2}
. The basal area was also measured by relascope inventory from the center of the plot. The measurement of these two variables allows the field work to be done cheaply and quickly. Later, the Vob was calculated for each plot applying the Vob tariff in accordance with Ho and G calculated with the plots from the IFN4. The values and statistics of the variables used in the measurements of the plots from the IFN4 (2011) and the new plots (2015) are shown in

The LiDAR data used came from the flight over the Autonomous Community of the Basque Country in July, 2012, with a mean density of 1 pulse m
^{-2}
and up to 4 returns per pulse.

For the entire data-gathering area, those

In order to make this clip, the PolyClipData tool of the FUSION program (

The same procedure was followed using the LiDAR data for 2011 and the plots measured in the field in 2015. Fifty-five cuts of LiDAR points were obtained with their statistics and 55 Vob data calculated for each plot measured in 2015. Forty plots were selected to fit the model and the validation group was formed by the 15 plots. It should be pointed out that none of the plots measured in 2015 was located in areas in which clear cutting had been performed during the 5-yr period 2011-2015. This meant that the fitted model did not consider the possibility of volume 0. Some plots were indeed located in areas that had been subjected to thinning during the 5-yr period.

Clear cutting and thinning data in the municipality of Muxika between 2011 and 2015 was gathered in order to be able to use this to estimate the growth in wood volume during the period 2011-2015. Data on harvest provided by the Provincial Council of Bizkaia (the administrative office that manages the forest stands of the Province of Bizkaia), can be seen in

In order to adjust the Vob tariff per hectare in accordance with the stand variables G and Ho, 90 plots of the IFN4 were used. Of these, 60 were used to fit the model and 30 to validate this. Plots were assigned to the model group or the validation group according to the random number assigned to them. In order to build the tariff, the STATGRAPHICS (
^{2}
, G, G
^{2}
and Ho*G were used as independent variables (

The best model was chosen with the following indicators: adjusted determination coefficient (
^{2}
), root mean square error (RMSE) and mean absolute error (MAE). The heterocedasticity was examined visually, by plotting residuals as a function of predicted values and using the Breusch-Pagan test (

Linear models were used to establish empirical relationships between field measurements and LiDAR variables. In order to fit the Vob for 2011, the Vob calculated in each IFN4 plot based on the model fitted with the independent variables G and Ho was used as the dependent variable. In the STATGRAPHICS regression model selection procedure, all statistics relating to the height of the points obtained with the FUSION CloudMetrics order were included as independent variables. In order to fit the model, 60 plots from the IFN4 were used. The other 30 plots were used to validate the model. The highest model efficiency in accordance with the LiDAR metrics was obtained by ordering the proposed combinations according to the best Mallows Cp. The same goodness-of-fit test utilized for the models in the previous section was used. Likewise, for the selected models, a residuals analysis was made for the selected models, with the aim of identifying evidence of normality deviations, lack of adjustment and/or heterocedasticity.

In order to fit the Vob for 2015, the Vob was used as dependent variable, calculated in each one of the 55 plots measured in the field in 2015 with the previously adjusted tariff of G and Ho. In order to maintain the uniformity of the results obtained, the model was fitted with the same independent variable as the Vob_2011 model, in other words, the 70 percentile of the height of the points cloud was used. In order to fit the model, 40 plots were used. The other 15 plots were used to validate the model.

In order to apply the wood volume models obtained, the area occupied by
^{2}
. The cells that intersected the polygons corresponding to the surface areas of

The lowest Mallows Cp value was 0.51. This corresponded to the model (3) that had a single independent variable: Ho*G (^{2}
(0.99), although it presents heteroscedasticity (

The lowest Cp value was 3.28 and this corresponded to the model that had the value of the height of the percentile 70 (h
_{70}
) of the points cloud (Eq. [4], ^{2}
(0.85) but was heterocedastic with a

The stem volume model for 2015 based on LiDAR data is shown in

The models obtained allowed us to estimate the total volume over bark with a different degree of accuracy, although in all cases it is quite high. The best regression model corresponds to the Vob model based on stand variables (H
_{0}
, G). This model presents the highest value of the adjusted coefficient of determination
^{2}
_{adj}
= 99), followed by the volume estimation model in accordance with LiDAR metrics (Vob_2011,
^{2}
_{adj}
= 85) and lastly, the model based on the field data obtained in 2015 (Vob_2015,
^{2}
_{adj}
= 80) (^{3}
ha
^{-1}
) and MAE (46 m
^{3}
ha
^{-1}
). Residual plots for Vob prediction indicated that the proposed model was appropriate, as can be seem from the predicted versus the observed values for Vob in

Each of the adjusted equations showed good fit statistics. All the parameters were significant at a confidence level of 95% (

The mean values of volume with bark of the three models are shown in

By applying the volume per hectare of wood with bark models to the entire surface area of Muxika for 2011 and 2015, the results shown in

The results of the distribution of the Vob for 2011 and 2015 and the difference between both can be seen in

The difference between the Vob estimated for 2015 and that estimated for 2011 is 81,603 m
^{3}
for the entire surface area of P. radiata in the municipality of Muxika. Therefore, the annual mean growth of the volume over bark (Vob) obtained during the 5-yr period was 5.14 m
^{3}
ha
^{-1}
yr
^{-1}
(

The forestry harvests data shows that the wood volume extracted during the period 2011-2015 was 135,624 m
^{3}
(^{3}
ha
^{-1}
were extracted. The estimated mean growth for the period 2011 - 2015 is 13.74 m
^{3}
ha
^{-1}
yr
^{-1}
, this being the value obtained by adding the estimated mean growth in the aforementioned period (Vob = 5.14 m
^{3}
ha
^{-1}
yr
^{-1}
) and the mean volume extracted 8.6 m
^{3}
ha
^{-1}
.

As can be seen in ^{3}
ha
^{-1}
. The increases diminish gradually for Vob values close to 300 m
^{3}
ha
^{-1}
and become negative above that value.

Both the
^{2}
and the RMSE of the model fitted for the Vob of 2011 are similar to those published by other authors under similar conditions (
^{2}
and RMSE of the Vob model for 2015 are more unfavorable than those obtained under similar working conditions. This is due to the fact that in the preparation of the Vob_2011 model, the field data and the LiDAR data were almost contemporary,

With regard to adapting the results obtained, the estimated mean growth being 13.74 m
^{3}
ha
^{-1}
yr
^{-1}
, the volume will have increased 68.7 m
^{3}
ha
^{-1}
on average during the period 2011-2015. Therefore, in those places where harvest was more than 68.7 m
^{3}
ha
^{-1}
, the growth has failed to offset the volume extracted and as a result will have less volume in 2015 than in 2011. These are the grid cells with a negative growth value. On the other hand, the mean volume extracted in a thinning operation according to ^{3}
ha
^{-1}
. It is for this reason that the areas that were subjected to thinning between 2011 and 2015 must correspond to cells with a volume difference of less than 29.5 m
^{3}
ha
^{-1}
as the 68.7 m
^{3}
of growth in the 5 yrs, minus the 39.2 m
^{3}
extracted in thinning operations will have allowed this amount to accumulate on average in these areas. Those grid cells with a volume difference of under that amount add up to 1,552 ha. This is the surface area in which, according to the fitted models, some kind of harvest took place during 2011 and 2015 and coincides with those contributed as real registered harvest: 1068.1 ha. As shown in ^{3}
ha
^{-1}
and, therefore, no plot in which clear cutting had taken place was measured. Consequently, the 2015 model cannot detect the clear cutting that took place during this period and assimilates these to intensive thinning, and therefore the surface area of harvest estimated by the models includes both clear cutting and thinning. This circumstance must be borne in mind for subsequent applications of the method and those plots subjected to clear cutting during the period of use of LiDAR data. It can be seen, moreover, that the areas in which the models predict greater harvest intensity correspond to those that had a higher volume in 2011. These areas are the ones in which third thinnings and clear cuttings (being the harvest of greater weight) are normally carried out. Both the assimilation of clear cutting to intensive thinning and the fact that the areas that were actually thinned or clearcut after July 2012 are captured in the raster models, which are based on a Lidar dataset from July 2012, can be explained when analyzing the diagram resulting from a comparison of Vob_2011 and the growth between 2011 and 2015 (

One of the aspects that must be clarified in subsequent studies is the development in growth in annual periods between both ends of the 5-yr interval. In this study it has not been possible to deal with this question due to a lack of field data in this intermediate period. Moreover, due to a lack of funding, locating the IFN plots was not possible for the second gathering of field data and for this reason the samples for 2011 and 2015 were different. Even though this circumstance might reduce the validity of the results obtained, it was decided to begin the study to verify whether this methodology could be used in subsequent studies with higher funding levels. Once sufficient funds are available, field data must be taken on an annual basis and in the same plot sample.

With regard to the multi-annual use and estimation of growth using a single LiDAR flight, no published references have been found concerning similar methods. This is the main contribution of this paper: the use of LiDAR data in successive years through the updating of the mass conditions with fast, reliable and cost-effective annual field data gathering. In this way, the investment involved in LiDAR data gathering is justifiable as this can be amortized in several years.

Fitted models that relate the descriptive variables of the points cloud of a LiDAR flight with the wood volume field data obtained in field plots with a 5-yr interval, provides a good approximation to the growth values produced in masses of
^{3}
ha
^{-1}
yr
^{-1}
), as well as an estimation of the surface area subjected to harvest in this period. This allows the multi-annual use of LiDAR data gathered in a campaign with the consequent decrease in costs. To ensure the highest levels of accuracy in the estimation of the growth and volume of wood during the years following the gathering of LiDAR data, it is necessary to include in the sample of plots to be measured those in which both thinning and clear cutting operations have been carried out.

Our sincere thanks to Hazi Fundazioa for providing essential data for this study.