**Applying LiDAR to Quantify the Plant Area Index Along a Successional Gradient in a Tropical Forest of Thailand**

#### **Siriruk Pimmasarn 1,\*, Nitin Kumar Tripathi 1, Sarawut Ninsawat <sup>1</sup> and Nophea Sasaki <sup>2</sup>**


Received: 31 March 2020; Accepted: 28 April 2020; Published: 6 May 2020

**Abstract:** Long-term monitoring of vegetation is critical for understanding the dynamics of forest ecosystems, especially in Southeast Asia's tropical forests, which play a significant role in the global carbon cycle and have continually been converted into various stages of secondary forests. In Thailand, long-term monitoring of forest dynamics during the successional process is limited to plot scales assuming from the distinct structure of successional stages. Our study highlights the potential of coupling airborne light detection and ranging (LiDAR) technology and stand age data derived from Landsat time-series to track back forest succession, and infer patterns in the plant area index (PAI) recovery. Here, using LIDAR data, we estimated the PAI of the 510 sample plots of a seasonal evergreen forest dispersed over the study area in Khao Yai National Park, Thailand, capturing a successional gradient of tropical secondary forests. The sample plots age was derived from the available Landsat time-series dataset (1972–2017). We developed a PAI recovery model during the first 42 years of the succession process. We investigated the relationship between the model residuals and PAI values with topographic factors, such as elevation, slope, and topographic wetness index. The results show that the PAI increased non-linearly (pseudo-R2 of 0.56) during the first 42 years of forest succession, and all three topographic factors have less influence on PAI variability. These results provide valuable information of the spatio-temporal PAI patterns during the successional process and help understand the dynamics of tropical secondary forests in Khao Yai National Park, Thailand. Such information is essential for forest management and local, regional, and global PAI synthesis. Moreover, our results provide significant information for ground-based spatial sampling strategies to enable more accurate PAI measurements.

**Keywords:** forest succession; LiDAR; leaf area index; plant area index

#### **1. Introduction**

Among the world's tropical forests, specifically Southeast Asia's tropical forests are considered to contribute significantly to carbon storage and climate mitigation, but they also belong to the major deforestation areas [1,2]. After deforestation, these forests feature patches of deforested areas and various stages of recovered forest through succession [3,4], where each stage has different forest structural attributes, species diversity, species composition, and capability in the forest ecosystem [3,5].

Currently, collecting spatiotemporal data, data on patterns, and data on the rate of forest recovery still represent a challenge, especially in tropical forests, while it is required for a better understanding

of the functionality and dynamics of forest ecosystems. For decades, many studies have attempted to explore the recovery of forests after disturbance [3,6,7]. Among the forest attributes, the leaf area index (LAI), defined as the one-sided leaf area per unit ground surface area [8,9], is one of the most significant biophysical variables [10] that is required for quantitative analysis of physical and biological processes related to vegetation dynamics and the modelling of ecosystem processes at local, regional, and global scales [11,12]. Therefore, long-term LAI monitoring following stand regrowth would further enhance our understanding of the dynamic changes in tropical forests and the climate impacts on tropical forest ecosystems.

Data on directly measured LAI in the field, especially in tropical forests, are very rare [13] due to the labor intensive work that is limited in time and space [14,15]. With the advancement of technology, indirect methods are now widely used, mainly through optical sensors such as the LI-COR LAI-2000 Plant Canopy Analyzer (which produced by LI-COR Inc., Lincoln, NE, USA) and digital hemispherical photography (DHP) [15,16]. However, the use of these methods is still limited in inaccessible sites, due to their spatial scale, and they often involve large uncertainties, especially related to the dense, tall, and heterogeneous canopies in tropical forests [14,15]. At regional and global scales, remotely sensed data meet the requirement for spatial and temporal estimates of LAI over large landscape areas. Most of the estimated LAIs via passive remote sensing data rely on the empirical relationships between field-based measured LAIs and spectral information [9,17,18]. Although passive satellite-based approaches are highly suitable for a wide range of observations, they suffer from several factors, such as weather conditions, the signal saturation of spectral reflectance in dense tropical forests, and an inability to capture the vertical LAI [10,19].

The light detection and ranging (LiDAR) technology collects the three-dimensional (3D) forest structure and helps to overcome many shortcomings of passive optical sensors. LiDAR has, hence, attracted much attention because of the advantage of monitoring forest structure and secondary succession of many different types [20–22]. Currently, LiDAR data is being widely and successfully used to derive LAIs [23–25]. In previous efforts, some studies estimated LAI from LiDAR data based on the empirical relationship between the ground-based measured LAI and LiDAR-derived metrics, such as mean height, maximum height, percentile height, and canopy density metrics [23–25]. Including the leaf, stems, twigs, and fine branches, the plant area index (PAI) is usually measured through the use of optical sensors and the value of PAI is comparable to the LAI value because leaf area is generally much larger than branch area, and the majority of branches are shaded by leaves [26]. With the advent of such voxel-based approaches, the term PAI instead of LAI was used. Due to a lack of recent literature on PAI, LAI will also be used for comparison where possible throughout the paper. Recently, voxel-based approaches have been developed to improve the estimation of PAI in tropical forests. With this approach, LiDAR data are reconstructed in the form of a 3D voxelized space of the forest structure and converted into PAIs by applying the Beer-Lambert law [27]. These approaches now become a promising method to obtain the PAI with high accuracy because they are not influenced by the orientation of the sun, or the spatial distribution, size, or shape of canopy components. Such methods can, thus, overcome the problem of the leaf clumping effect and signal saturation.

Although LiDAR data can be used effectively to retrieve the vertical PAI with high accuracy in tropical forests, they are limited in areas intended for long-term monitoring, as this technology is relatively costly. To investigate long-term PAI dynamics in successional gradients, data with a high temporal and spatial resolution is required. There are many studies that investigated the combination of high temporal resolution satellite images (e.g., Landsat) and high spatial resolution LiDAR data for estimating the dynamics of aboveground forest [28,29]. Recently, this approach has also been successfully applied to the assessment of forest aboveground biomass and its resilience in Thailand [30]. Until recently, however, there have not been any studies examining the long-term vertical PAI in successional gradients of tropical forests. Therefore, combining Landsat time-series datasets with high spatial resolution LiDAR data may be a promising way to achieve highly accurate PAI values for the past years along secondary forest successions.

The main aim of this study was to estimate the PAI and model the PAI recovery using LiDAR data combined with a Landsat time-series dataset in the secondary forest site in Khao Yai National Park, Thailand. Here, the landscape contains a mosaic of secondary forest with different ages surrounded by old-growth forest. The effects of topographic factors, including elevation, slope, and topographic wetness index (TWI) on the model residuals were also investigated. These factors are well known to influence forest growth and LAI, especially TWIs, which are widely used as topography-based indicator of soil-moisture. This is the first attempt to study the PAI and its recovery modeling in Thailand and Southeast Asian using LiDAR data in successional gradients.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The study area is located in Khao Yai National Park in central Thailand (Figure 1). It covers an area of approximately 64 km2, with an elevation ranging from 700 to 800 m above mean sea level [7,31]. At this altitude, the forest is seasonal evergreen with a dry season from November to April. The average annual precipitation is ca. 2100 mm, and the annual temperature ranges between 19 to 28 ◦C [7]. This area has been strongly affected by anthropogenic activities since the end of the 19th century to the establishment of the park in 1962, mostly through low-intensity agriculture. Since 1962, some areas have been maintained through the practice of open fires (for space opening to accelerate the growth of natural grasses) by the park managers. Therefore, the study area is a landscape mosaic consisting of forest patches of different ages, even if old-growth or relatively mature forest dominates the landscape [32].

**Figure 1.** Location of the study site in Thailand (**a**) and in the Khao Yai area (**b**). The map (**c**) represents the canopy height in the study site with three stand initiation stage plots (SIS; red dots), three stem exclusion stage plots (SES; blue square), Mo Singto or permanent old growth stage plot (OGS; red rectangle), and 444 Landsat time-series plots (LTS; black squares).

#### *2.2. Field Datasets*

Two forest inventory datasets with different successional stages were used in this study to estimate PAI change patterns with forest age. The first dataset included six inventory of 60 m × 80 m (0.48 ha) size, established from March to May 2013 [32]. These plots were located across the study area in different successional stages, including three plots in the stand initiation stage (SIS) and three plots in the stand exclusion stage (SES). The second dataset comprised a 30-ha Mo Singto plot, which is defined as old growth stage (OGS). The Mo Singto plot is a permanent plot established in 1996 in the Mo Singto area located in a large area of old-growth forest in the central part of Khao Yai National Park. This permanent plot is part of the global network of large forest plots of the Center for Tropical Forest Science (CTFS) of the Smithsonian Tropical Research Institute (STRI). The plot age was estimated by Chanthorn et al. [32], based on Landsat imagery and interviewing the old rangers, as 15–20 years for the SIS stage, 35–40 years for the SES stage, and more than ~200 years for the OGS stage. To homogenize its scale, the Mo Singto plot was subdivided into subplots of 0.5 ha, resulting in 60 plots with 50 m × 100 m.

#### *2.3. Forest Age Datasets*

The age of secondary forest pixels in the study area was obtained from Jha et al. [30], whose study was based on LiDAR and a Landsat time-series dataset (LTS). Jha et al. [30] employed a random forest algorithm to classify the Landsat time-series (1972 to 2017) dataset using training pixels derived from the mean height of the Canopy height Model (CHM) from LiDAR (2017) at a 60-m resolution. In total, 34 images from Landsat 1-3 MSS (1972–1983), Landsat 4-5 TM (1984–2011), and Landsat 8 OLI and TIRS (2013–2017) (http://glovis.usgs.gov) were classified into forest and non-forest areas with over 90% accuracy. Using 34 classified Landsat time-series images and applying the quality filters, Jha et al. [30] further selected 550 secondary forest pixels with various recovery ages, in which pixels experienced a shift from non-forest to forest areas. Methodological details on the classification and selection of secondary pixels with respect to age are given in Jha et al. [30].

Based on Jha et al. [30], and after excluding pixels that were not located in our footprint of 1 m LiDAR-derived CHM, we obtained 444 secondary forest pixels (LTS pixels) for use in this study.

#### *2.4. LiDAR Data Acquisition*

Airborne discrete return LiDAR (ALS) data was collected by Asian Aerospace Services Limited of Bangkok on 10 April 2017 over an area of 64 km<sup>2</sup> (Figure 1c). LiDAR data was acquired using a RIEGL LMS Q680i full-waveform laser scanner (RIEGL Laser Measurement Systems GmbH, Horn, Austria), installed into a Diamond Aircraft "Airborne Sensors" DA-42 fixed-wing plane (Asian Aerospace Services, Bangkok, Thailand). The flight altitude was approximately 500–600 m above the ground level with a 60◦ field of view. The pulse repetition frequency was 400 kHz. The average point density was 22 pts m2. Post processing of LiDAR data, including point cloud classification into ground, non-ground, and noise removal were done by AAS engineering (Aerospace Services Limited) using TerraScan, Terrasolid Version 14 (Terrasolid Ltd., Helsinki, Finland).

#### *2.5. Overall Methodology*

Firstly, all the sample plot coordinates were used to extract the point clouds for estimating the PAI. Secondly, all sample plot point clouds were used to estimate the PAI using the Amapvox software [27] and R programming software (which created by Ihaka, R. and Gentleman, R., University of Auckland, New Zealand). Then, the PAI was used to generate a recovery model by plotting against stand age. Simultaneously, the point clouds were also used to generate the digital surface model (DSM), digital terrain model (DTM), and canopy height model (CHM). The topographic factors were extracted from the DTM using the RSAGA package in R (which developed by Alexander Brenning, Jena, Germany). Finally, the topographic factors were used to assess their impact on the variability of the PAI recovery

model, while PAI data were used to assess the correlation with the CHM. The overall methodology is illustrated in Figure 2; details of the methodology are explained in the next sections.

**Figure 2.** Overall methodology.

#### 2.5.1. LiDAR Data and Topographic Factors Extraction

The point clouds classified as ground were interpolated to generate a DTM at1mresolution. The point clouds classified as non-ground points were used to generated a DSM with the same spatial resolution. A1mresolution CHM was then computed by subtracting the DTM from the DSM. Finally, we used the CHM to extract the canopy height at the plot level for assessing the correlation between forest canopy height and PAI data.

For creating topographic factors, the elevation was derived directly from the DTM, while both slope and TWI were computed from the DTM as the primary attribute using the RSAGA package in R. The slope was defined at each point in the DTM as a function of a gradient in the X and Y direction [33]:

$$\text{Slope} = \arctan\sqrt{\left(\text{fx}\right)^2 + \left(\text{fy}\right)^2} \tag{1}$$

TWI was calculated by combining the slope angle and specific catchment area (SCA), where SCA is considered as the factor that describes the tendency of the area to receive water [34]. The equation for TWI is computed as follows:

$$\text{TWI} = \ln \left( \frac{\text{SCA}}{\tan \phi} \right) \tag{2}$$

where *ln* is the natural algorithm, SCA is the specific catchment area, and φ is the slope angle.

#### 2.5.2. Estimation of LiDAR-Based Plant Area Index

We estimated the PAI of 510 sample plots, including 444 60 × 60 m pixels for which a forest starting recovery date was available, six 80 × 60 m successional sample plots, and 60 old-growth plots with an area of 50 × 100 m. For each plot, we first extracted all ALS points whose projected coordinates fell within the pixels using R Second, each plot containing point clouds was contructed in a form of a 3D voxelized space (m3). Then, a local vegetation transmittance of each voxel was computed as the ratio of the sum of existing energy and the sum of entering energy normalized by the mean optical path length in each voxel, which was implemented in AMAPvox (version 1.3.5) [27], available at http://amap-dev.cirad.fr/projects/amapvox/files?sort=filename%2Csize. The equation for estimating transmittance P of an individual voxel is:

$$P = \left[\frac{\sum\_{i}^{n} PFOut\_{i} \times S\_{i} \times l\_{i}}{\sum\_{i}^{n} PFEut\_{i} \times S\_{i} \times l\_{i}}\right]^{\frac{1}{\frac{1}{T}\sum\_{i}^{n} li}}\tag{3}$$

where P is the transmittance, *PFEnti* is the incoming fraction of pulse *I*, *PFOuti* is the exiting fraction of pulse *i*, *li* is the length of pulse *i* optical path, and *Si* is the cross section of a pulse at the voxel center.

To control the uncertainty in plant area density (PAD) estimation, which is related to low sampling intensity, the transmittance of each voxel was then refined using a hierarchical linear mixed model in R. Individual voxels, which were nested in a 5 × 5 × 5 m neighborhood, were treated as "random" effects, while the neighborhood was treated as a "fixed" effect. Subsequently, the PAD was computed for each voxel from the transmittance values by applying Beer-Lambert's turbid medium approximation, assuming isotropic transmittance as shown in Equation (4) [27]:

$$P(\theta) = \varepsilon^{-G(\theta).\text{PAD.I}} \tag{4}$$

where *P*(θ) is the gap probability for inclination θ (view angle/shooting angle), *G*(θ) is the ratio of foliage area projected in direction θ to actual, and *l* is the optical path length.

Finally, the PAI values were obtained from the vertical integration of the PAD profiles (Figure 3).

**Figure 3.** Example of PAD voxelization of LiDAR data point clouds applying AMAPVOX software.

#### 2.5.3. Statistical Analyses

Nonlinear regression analyses were used to explore the change in PAI along the forest chronosequence. The PAI recovery model was, thus, of the form:

$$\text{PAI} = a + b(T^{\mathfrak{c}}) + \mathfrak{e} \tag{5}$$

where PAI is the Plant Area Index, *a*, *b*, *c* are model parameters to be inferred, *T* is the forest age, and ε represents the model residuals.

For the statistical test, pseudo-R2 was used to assess the fitted non-linear model. To evaluate the effects of topographic factors (elevation, slope, and TWI) on the variation of the PAI (in terms of PAI residuals and PAI values) over 42 years and the correlation between PAI and CHM, correlation through stepwise multiple regression analysis were performed. PAI residual values of the PAI recovery model were individually regressed through ordinary least square regressions. All of the statistical analyses were performed using the statistical program R version 3.6.2.

#### **3. Results**

#### *3.1. Forest Successional Structure and Spatial PAI Variability*

An example of the three successional stand characteristics with the LiDAR-derived CHM and the frequency of canopy height per plot is shown in Figure 4. The LiDAR-derived CHM characteristics are clearly differed in each of the successional sample plots. The crown canopy size and tree height increased from the initiation to the old-growth stage following regeneration, while the canopy was dense in the exclusion stage only but sparse in the initiation and old-growth stages due to disturbed areas and gap formation. In addition, the canopy indicates the characteristics of forest succession by the vertical stratification along with its height. Both initiation and old-growth stages had a heterogeneous upper canopy surface due to canopy gaps and variations in height, while the exclusion stage showed a smooth upper canopy surface, with a lack of understory.

**Figure 4.** Example of the characteristics of the three forest successional stages. Each row represents the Canopy height Model (CHM), three-dimensional (3D) point clouds, and canopy height distribution of each stage, respectively. (**a**–**c**) represent the characteristics of SIS. (**d**–**f**) represent the characteristics of SES. (**g**–**i**) represent the characteristics of OGS.

Table 1 shows the characteristics of the successional sample plots, including mean PAI, mean canopy height, and stand age. We found that the PAI values varied with forest successional stage, canopy height, and age. The PAI values were lowest in the initiation stage plots and higher in the stem exclusion and old-growth stage plots, respectively. Moreover, the tree canopy height indicated that the stand height increased following the stand development stage. This result indicates that PAI values increase when the tree height increases. For all sample plots (secondary successional plots, Mo Singto sample plots, and LTS plots), the values of the secondary successional plots and those of OGS plots

seem to fit into the trend of the LTS sample plots (R<sup>2</sup> 0.75, *p* < 0.001; Figure 5). The mean PAI was 7.07 m<sup>2</sup> m−<sup>2</sup> and ranged from 2.69 to 11.02 m2 m−<sup>2</sup> (s.d. =1.85 m<sup>2</sup> m<sup>−</sup>2), while the mean canopy height was 16.57 m and ranged from 5.34 to 28.41 m (s.d. = 5.56 m).


**Table 1.** Means and standard deviations of the PAI, the canopy height, and the age.

**Figure 5.** Relationship between PAI and canopy height model for all datasets. The regression model is illustrated by the black line. Red dots represent the PAI in SIS from the successional sample plots, blue dots represent the PAI in SES from successional sample plots, green dots represent PAI in OGS from Mo Singto sample plots, and black dots represent the PAI from LTS sample plots.

#### *3.2. PAI Recovery Analysis*

The LiDAR-derived PAI recovery along the successional gradient is illustrated in Figure 6. The LiDAR-derived PAI of the plots ranged from 1.36 to 11.02 m2 m−2, with a mean of 6.81 m<sup>2</sup> m<sup>−</sup>2. We found that the PAI accumulation increased non-linearly through time during the 42 years. The relationship between PAI and age was best modeled with a non-linear power model and an exponent higher than 1. However, as shown in Figure 6, the exponent value indicated a close to linear relationship with a pseudo-R2 of 0.56, indicating an increase in the PAI with recovery time during the 42 first years of the succession. Note that the model predicted a non-null PAI in year zero because we defined forests as areas with a mean top canopy height greater than 5 m, following the Food and Agriculture Organization of the United Nations (FAO) [35] definition of forests. After 20 years, the PAI was predicted to be 6.30 m<sup>2</sup> m<sup>−</sup>2, and 7.58 m2 m−<sup>2</sup> after 40 years. For the neighboring old-growth forest in which the forest age was more than 200 years, the median of the estimated PAI was 8.87 m2 m−<sup>2</sup> and ranged from 8.05 to 10.63 m<sup>2</sup> m−<sup>2</sup> (Figure 6).

**Figure 6.** Relationship between the PAI estimated from the LiDAR and past time after disturbance (grey dots). The red line represents the fitted power model. Red dots represent the PAI of SIS in successional sample plots. Blue dots represent PAI of SES in successional sample plots. The box plots represent the PAI values of the OGS in Mo Singto plots with a forest age of ~ 200 years. The box edges indicate the 25th and 75th percentiles, the solid line indicates the median.

Topographic factors showed mixed effects on both the PAI values and residual values. As shown in Figure 7 and Table 2, a weak trend was observed in the relationship between both the PAI values and residual values, and topographic factors such as TWI, slope, and elevation. Results from the stepwise multiple regression analysis are shown in Table 2 and an interpretation is described below Table 2.

**Figure 7.** The left colmn showed the relationship between PAI residuals and (**a**) topographic wetness index (TWI), (**c**) elevation, and (**e**) slope. The right colmn showed the relationship between the PAI values and (**b**) TWI, (**d**) elevation, and (**f**) slope. The color grediants represent the the past time after disturbance.


**Table 2.** Effects of topographic factors on PAI values and residual values using stepwise multiple regression analysis.

Note: \* Significant coefficients at the 95% confidence level; \*\* Significant coefficients at the 99% confidence level. S.E. = Standard error.

By the backward stepwise elimination method, the variable TWI was dropped from the PAI model, suggesting that TWI did not significantly influence PAI value. On the contrary, the slope had a significant and positive linear effect on the PAI value. When the slope increases by one unit, PAI is predicted to increase by 0.1016 unit on average, while controlling elevation.

Elevation had a significant quadratic effect on PAI (Table 2), where the effect of elevation on PAI varies by the level of elevation. Using differential calculus, the effect of elevation on PAI was expressed by the following equation:

$$\text{Marginal effect of elevation on PAI} = 2 \times 0.0001977 \times \text{Elevation} - 0.3090521 \tag{6}$$

This implies that at elevation = 700, the effect of elevation was −0.0323, indicating that when elevation increases from 700 to 701, PAI is expected to decrease by 0.0323 on average, while controlling for slope. Likewise, we obtained the effect of elevation as −0.0125 at elevation = 750, 0.0073 at elevation = 800, 0.0270 at elevation = 850, and 0.0468 at elevation = 900. As elevation increases from 700, the effect on PAI is initially negative, but increases and reaches zero when elevation is between 780 and 790. Thereafter, the effect turns positive and continues increasing. In other words, the PAI value took its minimum when elevation was between 780 and 790.

Likewise, the backward stepwise elimination method was used to analyze the PAI residuals after the time variable was accounted for. In this analysis, the slope variable was dropped, suggesting that slope did not significantly influence PAI residual. This also implies that time and slope were correlated.

TWI had a significant and negative linear effect on PAI residual. When TWI increases by one unit, PAI residual is predicted to decrease by 0.1104 unit on average, while controlling elevation. Elevation was not statistically significant at 5% significance level. This can mean that elevation had no significant effect on PAI residual. If we consider a 10% significance level, then quadratic effects were identified, where the effect of elevation on PAI residual varies by the level of evelation. Using differential calculus, the effect of elavation on PAI residual was expressed by the following equation:

$$\text{Marginal effect of elevation on PAI} = 2 \times (-0.0001096) \times \text{Elevation} + 0.1583402 \tag{7}$$

This implies that at elevation = 700, the effect of elevation was 0.0049, indicating that when elevation increases from 700 to 701, PAI residual is expected to decrease by 0.0049 on average, while controlling for TWI. Likewise, we obtained the effect of elevation as -0.0061 at elevation = 750, −0.0170 at elevation = 800, −0.0280 at elevation = 850, and −0.0389 at elevation = 900. As elevation increases from 700, the effect on PAI residual is initially positive, but decreases and reaches zero when elevation is between 720 and 730. Thereafter, the effect turns negative and continues increasing. In other words, PAI residual took its maximum when elevation was between 720 and 730.

#### **4. Discussion**

#### *4.1. Spatial PAI Variability*

The spatial variability in the PAI was derived from the LiDAR-derived PAI for all sample plots. For successional sample plots, we found that the PAI values differed greatly among the three successional stages, which related to canopy height and stand age. The PAI values were higher in old-growth plot than in stem exclusion and stand initiation stage plots (Table 1). The old-growth stage included a larger number of tall trees with random gaps, while the stem exclusion stage included trees with a smaller crown size and lower height, and the stand initiation stage included the smallest and lowest number of trees (Figure 4). This result indicates that differences in forest structure have a strong impact on the PAI value. For the LTS plots, the results also showed that the older areas had higher PAI values (Figure 8a,b). For all 510 plots, the PAI varied from 2.69 to 11.02 m<sup>2</sup> m−<sup>2</sup> (mean = 7.07). The coefficient of variation was 26.12%, indicating a high spatial variability of the PAI at the investigated sites. Compared to other studies, our PAI values were in the range (2.66–12.94) reported by Clark et al. [36], who conducted the first direct measurements in the heterogeneous landscape of the tropical forest, and Vincent et al. [27] and Tang et al. [22], who both obtained PAI via LiDAR technology. With respect to the tree height, the relationships between the PAI and CHM were significant (R2 = 0.75, *p*-value = 0.001) (Figure 5), indicating that the stand height is a key variable in understanding the spatial distribution and variability of the PAI over large areas.

**Figure 8.** (**a**) Map showing the past time after disturbance of LTS sample plots. (**b**) Map showing the PAI values of LTS sample plots.

#### *4.2. Long-Term PAI Accumulation Through Succession*

The LAI is considered to be one of the most valuable determinants of forest growth and productivity. At the same time, long-term monitoring during the succession of the LAI in tropical forests is very scarce. We demonstrated that the combination of LiDAR and data of the classified forest age obtained from time-series Landsat data can quantify the long-term PAI accumulation over 42 years of succession in a tropical moist forest in Khao Yai National Park, Thailand. As a result, our analyses showed that a non-linear power model with an exponent greater than 1 was the best fit to our data, indicating an increase in the PAI during the first stage over the 42 years of succession (Figure 6). In agreement with the model, the estimated PAI of the successional field plots, which were used for validating our model, also showed that the PAI increment following the succession stages and their values were close to our predicted PAI (Figure 6). In addition, our study showed that the PAI increased gradually in the first 20 years and continued to increase twofold after 40 years. Compared to the PAI estimated for adjacent old-growth forests (mean = 8.92 m<sup>2</sup> m<sup>−</sup>2, more than 200 years), the average PAI (mean = 8.07 m2 m−2) after 42 years of succession decreased by 0.85 m2 m−2. These differences were very small when compared to the deviation in forest age. Accordingly, our results suggest that the PAI would increase rapidly during the first 42 years of forest succession, while it would only slightly increase thereafter and remain constant afterwards. In tropical forests, LAI predictions along successional gradients, assuming the distinct structure of successional stages or ages, are limited and were only reported on the plot scale (mostly ~0.5–1 ha). Chanthorn et al. [7] reported that the estimated LAI

differed considerably between successional stages, i.e., 3.83, 5.47, and 3.81–5.74 m2 m−<sup>2</sup> for very young forests, smooth young forests, and old-growth forest, respectively. Although we conducted our study in the same study area, the LAI is very different due to the different methodology and sensor used. Tang et al. [22] reported that LAI values were estimated to be 1.74, 5.20, and 5.62 m<sup>2</sup> m−<sup>2</sup> for open pasture, secondary forest, and old-growth forests, respectively. These studies showed a similar pattern when compared to our study, i.e., the LAI increased rapidly during successional stand development, and, then, the increase weakened during the transformation of secondary to mature forests. However, the estimated LAI values of these studies differs greatly from that of our study (Table 1), and the rate of recovery remains unclear because the LAI values were estimated with a different methodology, successional stage definition, and in a different study site.

#### *4.3. Topographic Factors Influence the Variation in the PAI Increase*

Topographic factors were found to be critical factors in determining the variation in the PAI across the landscape in several previous studies [37]. In those studies, topographic factors were proposed to be the most important factors contributing to tree growth and LAI variation in tropical forests. Moser et al. [38] found that the LAI decreased by 40%–60% in the elevation range of 1000–300 m. Unger et al. [39] reported that the LAI significantly decreased in the elevation range of 500–2000 m by about 1.1 per 1000 m. Liu et al. [40] reported a positive correlation between the LAI and soil moisture. However, this positive correlation was only found in clay-rich soil, which indicates that soil texture is one of the relevant factors determining the soil moisture-LAI relationship [41]. At our site, all three topographic factors, TWI, slope, and elevation showed week relationship with the PAI residual values and PAI values, indicating that the variability in the PAI in this study area was less influenced by these factors (Figure 7, Table 2). It has been hypothesized that the variability in the PAI is due to the variation in soil nutrients and properties during succession. In addition, our study site is a special case in that it is a natural ecosystem, affected by anthropogenic activities, such as agriculture, before the establishment of the park, and maintained as pasture land by fire. The study area comprises a mosaic of secondary forest patches, with different age and surrounded with old growth forest. Therefore, some of the sample plots are in the transition zone and/or large trees were partly left in the secondary forests, which probably affected the PAI variability in the recovery model.

#### **5. Conclusions**

Study on PAI patterns in the successional forests after disturbance is rare in the tropics. This study demonstrated the combined use of LiDAR data and stand age data derived from Landsat time-series to determine PAI patterns in a forest succession in Khao Yai National Park, Thailand. The PAI values non-linearly increased following forest successional stages with a rapid increase during the first 20 years, twofold increase during the first 40 years or, so but a constant slow increase after 200 years onward. Effects of the topographic factors on PAI values and residuals are less significant or weak in our study.

These findings provide an information of the long-term PAI recovery patterns during successional processes and the spatial variation in the PAI in heterogeneous tropical moist forests following disturbance. This information is very important for understanding the vegetation regrowth and the rate of change in a disturbed area as the large areas of tropical forests have experienced overexploitation, forest fires, and forest degradation. Our study findings can also provide the information on PAI values in tropical forests in Southeast Asia, where such values are still limited. The information on the spatial variation of the PAI in in the successional forests as a consequence of forest disturbance is critically needed for the development of ecological surveys and ground-based sampling strategies for ecological research and conservation. Further study on the effects of soil conditions and seasonal variations of vegetation on the PAI values could improve our understanding of the PAI patterns under various environmental and topographic factors.

**Author Contributions:** Conceptualization, S.P. and N.K.T.; methodology, S.P.; validation, S.P., N.K.T., S.N. and N.S.; formal analysis, S.P. and N.S.; investigation, S.P., N.K.T., S.N. and N.S.; resources, S.P. and N.K.T.; data curation, S.P. and N.K.T.; writing—original draft preparation, S.P.; writing—review and editing, S.P., N.K.T., S.N. and N.S.; visualization, S.P.; supervision, N.K.T.; project administration, S.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors gratefully thank the support of the National Science and Technology Development Agency (Thailand), Kasetsart University, the Department of National Parks, Wildlife and Plant Conservation (DNP), and the Institute of Research for Development (IRD). Takuji W. Tsusaka of the Natural Resource Management program at the Asian Institute of Technology is acknowledged for helping with the statistical analysis.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Shrub Biomass Estimates in Former Burnt Areas Using Sentinel 2 Images Processing and Classification**

**José Aranha 1,2,\*, Teresa Enes 1,2, Ana Calvão <sup>3</sup> and Hélder Viana 1,4**


Received: 1 April 2020; Accepted: 12 May 2020; Published: 14 May 2020

**Abstract:** Shrubs growing in former burnt areas play two diametrically opposed roles. On the one hand, they protect the soil against erosion, promote rainwater infiltration, carbon sequestration and support animal life. On the other hand, after the shrubs' density reaches a particular size for the canopy to touch and the shrubs' biomass accumulates more than 10 Mg ha<sup>−</sup>1, they create the necessary conditions for severe wild fires to occur and spread. The creation of a methodology suitable to identify former burnt areas and to track shrubs' regrowth within these areas in a regular and a multi temporal basis would be beneficial. The combined use of geographical information systems (GIS) and remote sensing (RS) supported by dedicated land survey and field work for data collection has been identified as a suitable method to manage these tasks. The free access to Sentinel images constitutes a valuable tool for updating the GIS project and for the monitoring of regular shrubs' accumulated biomass. Sentinel 2 VIS-NIR images are suitable to classify rural areas (overall accuracy = 79.6% and Cohen's K = 0.754) and to create normalized difference vegetation index (NDVI) images to be used in association to allometric equations for the shrubs' biomass estimation (R2 = 0.8984, *p*-value < 0.05 and RMSE = 4.46 Mg ha<sup>−</sup>1). Five to six years after a forest fire occurrence, almost all the former burnt area is covered by shrubs. Up to 10 years after a fire, the accumulated shrubs' biomass surpasses 14 Mg ha<sup>−</sup>1. The results described in this paper demonstrate that Northwest Portugal presents larger shrubland areas and greater shrub biomass accumulation (average 18.3 Mg ha<sup>−</sup>1) than the Northeast (average 7.7 Mg ha<sup>−</sup>1) of the country.

**Keywords:** sentinel 2; landsat; remote sensing; GIS; shrubs biomass; bioenergy; vegetation indices

#### **1. Introduction**

Portugal is an European country with a constituent land mass and 4 separate archipelagos. The former is located in the east of the Iberian Peninsula with an area of approximately 90,000 km2. From the mainland area (52%) there are: forest stands (39%), dense shrubland (12%), and sparse shrubland (1%) [1,2]. Between the mid-1980s and 2020, due to increasing human rural abandonment and edaphoclimatic conditions, a large number of forest fires occurred in mainland Portugal during the summer. The intensity of these fires increased dramatically each decade [3–7]. The same edaphoclimatic conditions that make the territory prone to wildfire occurrences, however, also create suitable ecological conditions for shrub regrowth after the fires. Previously published results [8–12] demonstrate that,

five years after wildfire occurrence, the fire scars are no long visible because they have been covered by shrubs as well as with the growth of scattered trees from self seedling processes.

Shrublands assume several diametrically opposed roles. On the one hand, they constitute the vegetable fuel that will eventually burn, and a social-economic problem. Portuguese law [13], specifies that it is mandatory to cut shrubs 10 m alongside the road network and 50–100 m beside other man-made structures on a regular basis. This has led to the use of fire as a means to eliminate the remaining cut area. On the other hand, these shrubs are also a source of carbon [12,14,15], they promote water and nutrients circulation in forested areas [8,16,17], protect the soil from erosion processes [16–19], support animal life, and promote biodiversity [8,19]. In essence, the ecological benefits to these shrublands could be summarised in two words: ecosystem services [20–22]. This is a central measure within European Community [22] and used as a way to assess forested areas' values. Thus such shrublands could take on an unanticipated new economic role potentially generating biofuel for power plants [23,24].

Analysis of shrubland location and its biomass accumulation is therefore important as it could influence the working processes for several stakeholders: forest management [25,26], wildfire hazard reduction, ecosystem surveying and biomass harvesting for energy production.

The calculation of forest biomass can be achieved using destructive processes, such as cutting and weighing vegetation in sampling plots. Subsequently, the results obtained can be analysed using appropriate geostatistics processes [23] that generate indicative biomass maps.

The results obtained through these destructive processes can later be used to adjust allometric equations enabling the estimation of biomass weight based on the volume collected. This is a non-destructive process. The results can also be used in conjunction with relevant satellite images. This can also support the calculation used in the allometric models for accumulated biomass estimation. This is achieved by comparing the relevant bands from the satellite images as well as analyzing the vegetation indices calculated by means of those same bands. This is a good example of a non-destructive method of estimating biomass.

It is important to note, however, that the first phase of a process to calculate and estimate accumulated biomass always begins with a destructive method.

The assessment of land cover dynamics in former burnt areas of forest as well as shrubs' regrowth must be carried out over vast areas of territory, and on an annual basis. It also requires the use of appropriate computer technology. Firstly, remote sensing techniques (RS) can be used for image processing and classification to create updated land cover maps. Secondly, geographic information systems (GIS) can be employed to record, manipulate, and present data. In addition, GIS allows the combination of multiple data sources enabling spatial analysis and can enhance the sampling process too. It is recognized that annual fieldwork for data collection is very expensive and time consuming, thus the use of RS and GIS provides a cheaper and appropriate sampling mechanism. Previous research in this subject area has generated several RS based approaches that used multitemporal satellite image classification and comparison.

A literature review on the use of GIS, RS, and combined RS/GIS for forest biomass and shrubs' biomass [25–29] enabled the consideration of different approaches to biomass estimation and mapping on given dates. This resulted in the use of particular regression models that were based on specific vegetation indices [30–34], and are presented in Table 1.

It should be noted that almost all of these previously presented data were not based on a shrubs' biomass time series sampling system thus do not account for any variation due to elapsed time over the former burnt areas. The use of allometric equations adjusted for a given geographic area requires local validation before it is used elsewhere, i.e., necessitating field work for data collection and the use of mathematical models for data analysis.


**Table 1.** Biomass regression models based on vegetation indices.

After 18 years of carrying out fieldwork to measure the volume and the weight of shrublands, the authors' main aim now is to present a suitable methodology that enables the estimation of the accumulated shrubs' biomass. This process takes into account the elapsed time after wildfires, based on satellite imagery processing and classification. The methodology is non-destructive and does not require fieldwork for data collection, thus allows accurate estimates when used in conjunction with RS techniques. It also enables stakeholders to perform dynamic analysis using satellite images in time series processing. To achieve this main aim, the authors adopted a methodology using Sentinel 2 images processing and classification as a way to identify former burnt areas, shrubland and to adjust an allometric equation that enables to estimate shrubs' biomass through using NDVI images.

This methodology also incorporates the elapsed time after identifying any wildfire occurrence effect on the shrubs' regrowth as an estimate. This then enables the creation of accurate maps related to the shrubs' biomass accumulation. It also established that, if the growth rate in the Northwest area of Portugal is different from that of the Northeast area (comparison was using annual satellite images), then this may have an influence on the adjustment of allometric equations. In addition this methodology has used dynamic models to identify forest fire hazards and also used design logistics models for biomass harvesting for energy purposes. These can indicate the areas where the intervention of forest managers (necessary to comply with the law) has taken place. Another unexpected result of this methodology was to help design post fire ecosystem recovery actions.

#### **2. Materials and Methods**

#### *2.1. Study Area Characteristics and Sampling Plots Location*

The study area is located in North Portugal (Figure 1) and comprises a forested area of 813,846 ha (432,000 ha forest stands and 381,846 ha shrubland).

**Figure 1.** Study area (North Portugal) and Portugal world geographical location.

This is a very fragmented landscape and the small forest areas are often side by side with shrubland and agricultural areas. It should be noted that, for cultural reasons, the rural population often uses fire as an instrument for pruning as well as a method for removing any residues. The result of this approach is that a higher number of rural fires every year occurs than is actually appropriate for this type of landscape.

The study area in Northern Portugal includes many morphological and edaphoclimatic conditions typical of this region. Altitude ranges from sea level (0 m) to 1546 m in the Gerês mountains (Figure 2). Mean annual accumulated precipitation ranges from 1000 to 2400 mm in the Northwest areas and from 600 to 1200 mm in the Northeast areas. Mean annual temperature ranges from 12.5 ◦C to 15.0 ◦C in the Northwest areas and from 7.5 ◦C to 12.5 ◦C in the Northeast areas.

**Figure 2.** Sampling plots location.

#### *2.2. Data Sources*

#### 2.2.1. Using GIS in the Project

The project aims were to analyse former burnt areas, assess potential vegetation regrowth and to estimate the shrubs' biomass taking into account the elapsed time since any wildfire took place. The initial process used the burnt areas data. It is possible to download a set of vector files representing the burnt areas' perimeter by year of occurrence from the Portuguese Forest Services website [35] or the European Forest Fire Information System [36].This information has been used to establish one of the layers within the GIS project since 2000. Every year, the burnt areas' vector file is updated with new burnt areas boundaries. Spatial analysis then enabled new calculations to be made such as identifying fire recurrence areas as well as assessing the time since the last occurrence of a fire in the same areas.

For the period between 2000 and 2016, 10 sampling plots per year after the last fire were selected. This resulted in 170 sampling points, dispersed throughout North Portugal. In 2017 and 2018, after the severe rural fires that occurred in those years, the GIS project was updated and new sampling plots were added, increasing the number of samples to 234.

These sample points (Figure 2) were then used to create a survey GIS project (e.g., ArcPad, Survey 123, QField) that was transferred to a DGPS receiver (Trimble Inc., Sunnyvale, CA, USA). All the sample points were also marked (based on the GIS layout) on the Topographic Plan of Portugal on a 1/25000 scale. These were then printed to support fieldwork in areas with no GNSS signal.

Thus the data generated from 2000 to 2018 were used to identify the sampling points on the ground enabling us to record any shrub regrowth since the last known fire occurrence.

Circular 500 m<sup>2</sup> sampling plots (12.62 m in radius) were used along with the cross transect method. Two 25.24 m fiberglass tape measures were stretched perpendicularly across each sampling plot. Then, the shrubs intersecting each fiberglass tape were measured in 3 dimensions: length, width, and height. This enabled a calculation of the volume, assuming that the shrubs canopy was a sphere, and by using Equation (1).

$$\text{Vsh} = 1/6\,\pi \,\text{L W H} \tag{1}$$

where Vsh = shrubs canopy volume (m3); L = length (m), W = width (m), H = height (m) (Formula demonstration in Appendix A).

After measuring the total shrubs' volume along the 2 transects within the sampling plot, 10 shrub plants, 5 per transect (at the edges of the plot, halfway from the center and in the center) were cut in order to be weighed. They were then placed in plastic bags, brought to the lab, put to dry in the shade and weighed after reaching 30% moisture. The achieved results for the 500 m2 circular sampling plots were then extrapolated to an area of 1 hectare and the amount of shrub biomass per plot was estimated using Equation (2).

$$\text{Biormass} = \text{V } \text{W } 200 \text{\AA} \text{xw} \tag{2}$$

where Biomass in Mg ha<sup>−</sup>1, V = total shrubs volume along the 2 transects (m3), W = average shrubs weight (Mg m−<sup>3</sup> at 30% moisture), Mxw = maximum width measured along both transects (Formula demonstration in Appendix A).

2.2.2. The Processing and Classification of Sentinel 2 Images

Sentinel 2 images were freely downloaded either from the Copernicus Hub website [37] or from the Glovis website [38].

As the study area is not covered by a single image, eight Sentinel 2 images were used, two per year, for the years of 2016, 2017, 2018, and 2019. All of them were recorded in the summer season.

It was not possible to download Sentinel-2S2A images for the years 2017, 2018 and 2019. Only Sentinel-2L1C were available. When processing these images for the various dates, it was noticed, after computing the NDVI images, that the calculated values for the water surfaces (e.g., dams) were not consistent. This situation led to performing Sentinel-2 SEN2COR280 Processor (7.0.0) analysis by SNAP. For this reason, an images atmospheric calibration was performed, based on the spectral signature of the water collected by the research team with a spectra radiometer Ocean Optics VIS NIR (Spectrecology, Inc., St Petersburg, FL, USA) and on the atmospheric scattered model proposed by Chavez (1988) [39]. This procedure was carried out in order to have the same spectral signature for all water surfaces. Each image was also submitted for atmospheric correction because solar elevation and

the state of the atmosphere introduce differences in the radiation detected by the sensor for the same area on different dates [40–43].

To ensure that differences in reflectance are due to changes in land cover, and not caused by radiometric distortions, it was also necessary to apply a radiometric correction. One of the most used models for atmospheric correction is a process called dark object subtraction (DOS) proposed by reference [39]. This process is based on an atmospheric scattering model and reduces the haze effect by calculating the expected minimum for a given band after atmospheric correction. This was carried out in relation to the following criteria: at-satellite radiances were converted to surface reflectance by correcting for both solar and atmospheric effects. Then, at-satellite radiance values were converted into surface reflectance using a DOS approach [39]. This assumes no atmospheric transmittance loss and no downward diffuse radiation. The surface reflectance of the dark object was assumed to be 1%, and thus the path radiance was assumed to be the dark-object radiance minus the radiance contributed by 1% surface reflectance [39–43].

Spectral reference signatures, such as water, bare soil and dense shrubland, were created after dedicated work using an Ocean Optics VIS NIR spectra radiometer.

After the images processing operations, a RGB false colour composition image and an NDVI image was created for all dates. The Sentinel 2 RGB482 composition was used, because it employs the near infrared band in the green channel which allows the highlighting of vegetation thus enabling an identification of the burnt areas, both recent and old.

Based on the analysis of these new images, and with the support of the GIS project and the orthophotomaps made available by Bing Maps, spectral signatures were created to support the images supervised classification. In a second stage, spectral signatures for the main rural and forest land cover classes were created, namely:


These spectral signatures were then used to perform supervised classification techniques using the 10 m Sentinel 2 bands: B2, B3, B4, and B8. The minimum distance, maximum likelihood, and random forest were tested.

After the supervised classification process completed, a new image was created indicating the burnt areas and shrubland. This is a necessary step whereby a raster mask is created then applied to the NDVI image in order to estimate the shrubs' biomass that has regrown on former burnt areas.

This raster mask is required because the vector files represent all of the annual burnt areas created. When placed over the relevant satellite images showing the burnt areas boundary or perimeter, they do not consider the unburnt 'islands', nor the rocky areas [44–46]. Thus, to calculate biomass estimates by means of a vector mask may lead to overestimations. In order to be able to calculate the amount of shrub biomass that regenerated in the former burnt areas, it is necessary to create a raster mask that represents only the areas that were actually burnt.

Finally, the sampling points attribute table in the GIS was updated with the NDVI values using a known technique that enables extracting raster values to point-type vector files. Lastly, this attribute table was exported in dBase format and processed in Excel. This way, it was possible to analyse the relationship that could be established between the measured shrubs' biomass, and that which was calculated using dedicated field work. In addition an NDVI value was also calculated using satellite images processing techniques.

Figure 3 outlines a summary of all the working stages used in this research.

**Figure 3.** Work flow chart.

#### **3. Results**

#### *3.1. Landcover Characterization*

The landcover characterization work carried out was based on false colour RGB482 composition visual analysis (Figure 4) and used NDVI images (Figure 5) for interpretation (all dates). The ensuing results for 2019 are shown in Figures 4 and 5.

**Figure 4.** Sentinel 2 images (summer 2019) for the study area using false colour composition RGB482.

**Figure 5.** Sentinel 2 images (summer 2019) for the study area using an NDVI calculation.

As healthy vegetation has its maximum spectral reflectance in near-infrared wavelength, the Sentinel 2 band 8 (near-infrared) was used and coloured in green. Thus, green tonalities depicted in the RGB482 images indicated vegetation density and consequently the darker the green colour indicated the denser the vegetation.

The vegetation index NDVI was calculated using the normalized difference between the near-infrared and the red images. As the vegetation red reflectance is always lower to the near-infrared reflectance, the positive NDVI achieved values could also indicate vegetation density. Thus, the higher the NDVI values the denser the vegetation. Each of the dots depicted in Figure 2 represent a sampling point within a former burnt area. Analysing the NDVI image, Figure 5, it appears that the old burnt areas are in various states of vegetation recovery. It seems, therefore, that the Northwest area of Portugal has more dense vegetation and less burnt areas scars than the Northeast area.

It appears, however, that neither the green colour intensity shown in the RGB482 images (Figure 4) nor the NDVI values (Figure 5) enable a classification of the vegetation type (e.g., burnt areas, forest land, shrubland). These two images alone only infer the potential density of the vegetation cover. Thus, it was necessary to carry out further analysis using Sentinel 2 images and an assisted classification process to classify the land cover in classes.

Using the results from the Sentinel 2 images supervised classification process as well as the Minimum Distance Classifier, enabled us to state that the main land cover features are suitable to be classified accurately as presented in Table 2. This classification accuracy is particularly high for forest features. For example, burnt areas and shrubland was easy to identify and classify with an accuracy over 80%.

It was possible, therefore, to create a mask image for use with an associated NDVI image in order to isolate burnt areas and shrubland. This mask was later used to 'cut' the NDVI image enabling a shrub biomass calculation to be carried out along with the allometric equation application.

With an overall accuracy of 79.6% and a Cohen's K coefficient of 0.754, it can be stated that the Sentinel 2 images were found to be suitable for use in forestry applications as well as in the dynamic analysis of former burnt areas too. The Sentinel 2 classified images created a raster mask that was used subsequently to isolate the burnt areas and the shrubs areas.


**Table 2.** Achieved results after confusion matrix for classification accuracy from Sentinel-2 images.

N: Number of ground control points; Pa: Producer's accuracy; Ua: User's accuracy; Ce: Commission error; Oe: Omission error.

#### *3.2. Allometric Model for Shrub Biomass Estimation*

For the allometric equation adjustment, 110 pairs (NDVI, shrub biomass) were used: 46 extracted from 2016 image, 33 from 2017 image and 31 from 2018 image. In a first approach to data processing, descriptive statistics were calculated for each date and region, as presented in Table 3.


**Table 3.** Descriptive statistic for the sub-samples.

Subsequently, student-t tests were performed in order to verify that the sample points for the NW area of Portugal are different to those of the sample points for the NE area. As no significant differences were found, all the sampling points for each year were then merged into a single sample file.

The NDVI and shrub biomass values were centered and reduced in order to verify if this composite sample had a normal distribution. Descriptive statistic and accumulated probability values were calculated. Initially, the descriptive statistic for NDVI and shrubs biomass was calculated, shown in Table 4. Then, accumulated probability values were calculated as depicted in Figure 6.


**Table 4.** Descriptive statistic for the total sample.

**Figure 6.** Normal cumulative curves to NDVI (top) and to shrubs biomass (bottom).

The achieved results show that both distributions presented a normal distribution (Figure 6).

In the third stage of the process, a XY graphical representation was used in order to analyse the relationship that could be established between NDVI values and shrub biomass (Figure 7). The resulting graphic shows a narrow points cloud on the left for the minimum values and a scattered cloud on the right for the maximum values. This indicates that there were constraints in the regression analysis for the allometric equation calculation, possibly suggesting that there were too few options available.

During the regression analysis processing, it was noted that the NDVI tends to saturate at 0.7, suggesting that the shrubs' growth process maybe asymptotic to approximately 50 Mg ha<sup>−</sup>1. This maybe due to the nature of the plants and the space they occupy over periods of time. As a consequence of the approach adopted, the ensuing model led to better results than anticipated. It is presented in Equation (3).

$$\text{Shrub biomass} = 70.078 \text{ NIDVI}^{2.8113} \tag{3}$$

R<sup>2</sup> = 0.8984 (*p*-value < 0.05) and RMSE = 4.46 Mg ha−<sup>1</sup>

**Figure 7.** Relationship between NDVI and shrub biomass.

#### *3.3. Shrub Biomass Estimation Using NDVI Image Processing*

When the regression analysis completed, the adjusted allometric equation was used to generate the shrubs' biomass estimation using the NDVI. It was also applied to the 2019 Sentinel 2 images before the final calculation. General NDVI images were submitted to mask extraction in order to create new images indicating the former burnt and shrub areas, as depicted in Figure 8.

**Figure 8.** Accumulated shrubs' biomass in former burnt areas in the North of Portugal, estimated using an NDVI image summer 2019 and an allometric equation.

The results show that it is possible to account for an extent of some 172,022 ha and of 1,323,222 Mg of accumulated shrubs' biomass in the Northeast area. Likewise some 209,824 ha and 3,835,047 Mg can be identified in the Northwest area. On average, the Northeast area has approximately 7.7 Mg ha−<sup>1</sup> and the Northwest has 18.3 Mg ha−<sup>1</sup> of accumulated shrubs' biomass identified in the former burnt areas and shrubland.

#### **4. Discussion**

#### *4.1. Sentinel 2 Images*

The Sentinel 2 images, bands B2, B3, B4 and B8 (VIS–NIR), were deemed to be suitable for use in rural areas characterization and mapping, mainly in forested and shrubland areas, as demonstrated by the calculated overall accuracy (OA = 79.6) and the Cohen's K coefficient (k = 0.754) described earlier.

It was possible to classify, with good accuracy, the forest features in deciduous, coniferous and shrubs areas too. It was not possible, however, to classify mixed forest stands of deciduous and coniferous woodlands as the classification methods only isolate deciduous clusters from coniferous clusters. It was also not possible to classify forest stands by species.

Due to their spectral signatures, burnt areas were identified and classified with an users' accuracy of 89%. It was only possible to classify these former burnt areas, i.e., younger than a year and a half, after the fire because this timeframe indicates when a vegetative period has taken place and the shrubs were beginning to grow in these burnt areas. For example, if two full years have passed since the fire, the former burnt areas classification through satellite images starts to present results that confuse these areas with agricultural land and rocky areas with scattered vegetation.

The B4 and B8 bands enabled us to calculate NDVI images which proved to be adequate for the shrubs' biomass characterization and quantification. This was demonstrated using the statistical values of correlation between NDVI and shrubs biomass (r = 0.853) and the determination coefficient for the allometric equation (R<sup>2</sup> = 0.8984).

#### *4.2. Shrubland Characterization*

Vegetation, mainly shrubs, has a great potential to regrow on former burnt areas. The capability to colonize the space and to produce biomass is closely related to local morphology and edaphoclimatic conditions. As previously presented in Table 3, after calculating descriptive statistics for sub-samples and after adjusting allometric equations to each year, no statistically significant differences were found. In order to guarantee the accuracy of the estimates obtained by the general allometric equation now presented, however, it is necessary to verify if the growth rate of the shrubs is different in the two study areas. To achieve this, an allometric equation per area to the pairs was developed using: age and shrubs' biomass. This is presented as Equations (4) and (5).

$$\text{NW region: Shurb biomass} = -0.0062 \text{ t}^3 + 0.2089 \text{ t}^2 - 0.2738 \text{ t} + 2.279 \tag{4}$$

R<sup>2</sup> = 0.7349 (*p*-value < 0.05) and RMSE = 4.9 Mg ha−<sup>1</sup>

$$\text{NE region: Strub biomass} = -0.0072 \text{ t}^3 + 0.2211 \text{ t}^2 - 0.2738 \text{ t} + 2.043 \tag{5}$$

R<sup>2</sup> = 0.6921 (*p*-value < 0.05) and RMSE = 3.7 Mg ha−<sup>1</sup>

The resulting shrubs' estimates, by means of these two equations, showed that there was no statistically significant differences found between regrowth rate in either of the study areas. This is shown in Figure 9.

The Northwest study area was found to have better morphological and edaphoclimatic conditions for shrubs' regrowth after fire than the Northeast area, possibly because the latter had retained burnt area scars for a longer time. After 30 years of forest fires, many of the mountainous areas had lost almost all vegetation and, as a consequence, top soil too. Figures 4 and 5 show that these areas are now characterized by small forested spaces or shrub zones surrounded by rocky extents with scattered shrubs.

**Figure 9.** Estimations for shrubs regrowth rates in the study area.

The achieved results were classified regarding the average value of accumulate shrubs' biomass according to the elapsed time after the last fire event and the sampling date. This reclassification process enabled us to calculate the accumulated shrubs' biomass amount per class and also to create a histogram for the number of available hectares per class. The results are shown in Figure 10.

**Figure 10.** Accumulated shrubs biomass in summer 2019.

Although the growth rates in the two areas are not significantly different, it appears that the Northwest area is more densely vegetated than the Northeast. This was expected since this former area is facing the Atlantic Ocean, has a milder climate and greater rainfall than the Northeast landlocked territory.

It can be noted that for both of the study areas, within two years, the vegetation was capable to regrow enough to disguise the black landscape caused by fire. After five years, almost all the former burnt area was covered by vegetation and the accumulated shrubs' biomass grew up to 5 Mg ha−<sup>1</sup> (30% moisture). Thus it can be demonstrated that in up to five years after a fire occurrence, the ecosystem will recover as evidenced by the fire hard index ranges from very low to medium. It appears that between five and 10 years after a fire, the accumulated shrubs' biomass can reach 14 to 18 Mg ha−1. In terms of the ecosystem, the situation is favorable, but in terms of fire danger less so. Between 10 and 15 years after a fire, the shrubs' accumulated biomass can reach 26 Mg ha<sup>−</sup>1. It must be noted, however, that this biomass already has a large wood structure that gives it properties suitable for its use as fuel for thermoelectric power plants. This means that this shrubs biomass has a potential economic value that could change very dense shrubs areas from a severe fire hazard issue into a green fuel source.

#### *4.3. Allometric Equation for Shrub Biomass Estimation*

The accumulated shrubs' biomass estimation was made using an allometric equation based on the elapsed time after a fire, as previously presented in Section 4.2 and also shown in Equations (4) and (5) as well as Figure 9. In addition the NDVI values derived from satellite image processing, also played a significant role (described earlier in this paper).

When working in a large area that presents different morphological and edaphoclimatic characteristics and which requires the use of multiple satellite scenes, it is appropriate to verify in advance if there are differences in the shrubs' growth rate and if there are differences in relation to the year under study. It was verified, as previously stated, that no statistically significant differences were found in relation to the shrubs' growth rate of the bush in the two study areas.

To verify the second hypothesis, an allometric equation for each date was adjusted. This is presented in Table 5.


**Table 5.** Allometric equations adjusted for the 3 years in analysis.

Each of the allometric equations used the 110 sampling points and consequently the results generated were very similar with no significant differences found between estimates. These are shown in Figure 11.

**Figure 11.** Relationship between NDVI and shrubs biomass estimates.

As no statistically significant differences were found between the shrubs' biomass estimations using the three previously presented equations, a general equation was developed using the 110 sampling points (Equation (3)).

The adjusted allometric equation described here proved to be suitable for assessing the shrubs' biomass estimation using NDVI values. Employing this method also made it possible to monitor the shrubs' biomass regrowth in the former burnt areas on a regular basis by means of Sentinel 2 image processing and classification.

As neither the shrubs' growth rates are significantly different for either study area or the use of the allometric equations (adjusted for each of the years: 2016, 2017, and 2018) led to statistically significant estimates between them and to the general equation, it can be state that the methodology was suitable to be used in this area for any year. This has been demonstrated by using accurate estimates of the accumulated shrubs' biomass (based on satellite imagery) and on a regular basis. It also enables the monitoring of the shrubs' biomass variation as well as calculating the biomass gains and losses. Biomass gains can be converted into sequestered carbon and used to analyse the ecosystem's state of health as well as its production capability [47,48]. It can also be used to update the fire hazard indices calculation and to identify the places most prone to be burnt by large fires and, therefore, requiring special attention. Biomass losses can be converted into carbon released into the atmosphere by forest fires [9] or used to calculate the intensity of forest fires. The difference between post-fire gains and losses can be used to calculate the fire severity [46,49] in any given spot.

This equation was adjusted for the North Portugal study areas and it is suggested that it could be used in other territories with similar morphological and edaphoclimatic conditions. It is important to remember, however, that such work requires a validation process for any area.

#### **5. Conclusions**

Free Sentinel 2 images were an asset to derive multi temporal and dynamic studies about land cover, for monitoring former burnt areas and to estimate the shrubs' biomass accumulation.

The achieved results enable us to state that the methodology presented in this manuscript proved to be robust and that the NDVI derived from Sentinel 2 images can be used to calculate accurate and dynamic estimates of accumulated shrubs biomass.

The allometric equation presented here also allowed us to estimate the shrubs' biomass using Sentinel 2 images without depending on the vector files provided by EFFIS or by ICNF (Portuguese Institute for Nature and Forest). Comparing the shrubs' biomass estimates achieved through the two allometric equations, using NDVI and using elapsed time after fire, demonstrated that no statistically significant differences were found. Thus, it can be stated that the allometric equation presented in this manuscript incorporated the effect of elapsed time after the fire.

The combined use of GIS and RS techniques, complemented by regression analysis proved to be useful for monitoring the shrubs' regrowth in the former burnt areas. It was also useful to analyse the land cover dynamics and also to quantify the accumulated shrubs' biomass. GIS supported data records, management, and the sampling system development was helpful whilst RS supported multitemporal land cover analysis and biomass estimation using associated satellite image processing. Classification also offered major benefits too.

Regression analysis and allometric equation adjustments were found to be suitable processes to assign the biophysical data that was collected via field work as well as the use of the freely available satellite images. This led to the calculation and estimatation of the shrubs' biomass.

Although the NDVI saturates only measured 0.7, it was still possible to obtain good estimates of the shrubs' biomass before the complete canopy closure which is when the accumulated values reach their maximum. It was discovered that the NDVI values were not, however, specific for all types of vegetation. Consequently it was always necessary to create a raster mask in advance that identified the type of vegetation or area under analysis in order to define the estimates for those particular places of interest.

From a forest management perspective, it was found that, after five years, the accumulated shrubs' biomass starts to be a fire hazard related issue as it creates a horizontal continuous coverture which encourages any fire to spread. If it was 10 years since the last fire occurrence, the amount of accumulated shrubs' biomass was found to be over 14 Mg ha−<sup>1</sup> which led to a severe wild fires spread. This is commonly referred to as a 10 year of fire recurrence cycle in Portugal.

The methodology presented in this paper was found to be suitable for use in forest land management and also served a number of unexpected different purposes. From an ecological persepective it has been demonstrated that, over a two-year period, the vegetation was capable of enough regrowth to minimize erosion actions and to support animal life. In a five-year period, it appears that almost all the former burnt areas are covered by vegetation. From this point of view it may be possible to use the shrubs' biomass for energy purposes but it was found that only after a 10-year period that the amount of accumulated shrubs' biomass became economically valuable to cut

and transport elsewhere. This was determined by the monetary biomass value for any potential power plant location, the man hours cost involved to capture it as well as the necessary transportation costs.

**Author Contributions:** Conceptualization, J.A.; performed the experiment and analysis, J.A., T.E., H.V. and A.C.; performed the statistical analysis, J.A., T.E., A.C. and H.V.; resources, H.V., A.C. and J.A.; writing—original draft preparation, T.E., J.A.; writing—review and editing, J.A., H.V.; project administration, J.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is supported by National Funds by FCT—Portuguese Foundation for Science and Technology, under the project UIDB/04033/2020.

**Acknowledgments:** The authors would like to express their sincere thanks for the critical review of the manuscript carried out by Teresa Connolly.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Additional information:


Formulae demonstration:

$$\text{Vsh} = 1/6\,\pi \,\text{L W H} \tag{A1}$$

where Vsh = shrubs canopy volume (m3), L = length (m), W = width (m), H = height (m), and the volume of a sphere = 4/3 π r3.

In order to use the shrubs dimensions measured along the transect, the equation can be rewritten as:

Volume of shrub canopy= 4/3 π Length/2 Width/2 Height/2

Volume of shrub canopy= 4/24 π Length Width Height

Volume of shrub canopy= 1/6 π Length Width Height

$$\text{Biomass} = \text{V } \text{W } \text{200/Mxw} \tag{A2}$$

where Biomass in Mg/ha V = Total shrubs volume along the 2 transects (m3) W = Average shrubs weight (Mg/m<sup>3</sup> at 30% moisture) Mxw = Maximum width measured along both transects.

In a 500 m<sup>2</sup> sampling plot, the plot radius is 12.62 m. This way, each transect has 25.24 m. The 2 transect account for 50.48 m.

When measuring the shrubs along these transects, the sum of shrubs width plus shrubs length define the area occupied by shrubs within the 2 cross section transects. Using the maximum measured shrubs width and the length of both transects is possible to calculate the maximum area for the 2 transects.

Sampling area = Maximum width 50.48 m

For converting the measured shrubs biomass within the sampling area, it is necessary to convert this area to 1 hectare.

10,000/Maximum width 50.48 ≈ 200/Maximum width

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
