Monitoring clearcutting and subsequent rapid recovery in Mediterranean coppice forests with Landsat time series

Abstract

Key message

This work analyses the rate of recovery of the spectral signal from clearcut areas of coppice Mediterranean forests using Landsat Time Series (LTS). The analysis revealed a more rapid rate of spectral signal recovery than what was found in previous investigations in boreal and temperate forests.

Context

The rate of post-disturbance vegetation recovery is an important component of forest dynamics.

Aims

In this study, we analyze the recovery of the spectral signal from forest clearcut areas in Mediterranean conditions when the coppice system of forest management is applied.

Methods

We used LTS surface reflectance data (1999–2015).We generated an annual reference database of clearcuts using visual interpretation and local forest inventory data, and then derived the Normalized Difference Vegetation Index (NDVI) and Normalized Burn Ratio (NBR) spectral trajectories for these clearcuts. From these spectral trajectories, we calculated the Years to Recovery or Y2R, the number of years it takes for a pixel to return to within a specified threshold (i.e., 70%, 80%, 90%, 100%) of its pre-disturbance value. Spectral recovery rates were then corroborated using measures of canopy height derived from airborne laser scanning (ALS) data.

Results

The coppice system is associated with rapid recovery rates when compared to rates of recovery from seeds or seedlings in temperate and boreal forest conditions. We found that the Y2R derived from the spectral trajectories of post-clearcut NBR and NDVI provided similar characterizations of rapid recovery for the coppice system of forest management applied in our study area. The ALS measures of canopy height indicated that the Y2R metric accurately captured the rapid regeneration of coppice systems.

Conclusion

The rapid rate of spectral recovery associated with the coppice system is 2–4 years, which  contrasts with values reported in boreal and temperate forest environments, where spectral recovery was attained in approximately 10 years. NBR is an effective index for assessing rapid recovery in this forest system.

1 Introduction

Forests cover 30.6% of the world terrestrial surface (FAO 2015). In Europe, forest area continues to increase with forests now accounting for 33% of the total land area (FOREST EUROPE 2015). The importance of forest ecosystems is related to their capacity to provide a vast array of ecosystem goods and services, which may be altered by climate and land use changes. Within Europe, Mediterranean forests are more potentially at risk of suffering the impacts of climate change because of drought and augmented forest fires (Schröter et al. 2005).

Since the direct economic value of European forests is primarily related to wood production, logging is both one of the main methods for producing ecosystem services, as well as one of the main disturbances to these ecosystems (Seidl et al. 2011), further underpinning the need to manage forest ecosystems sustainably. The main indicator for measuring the sustainability of forest management has conventionally been the ratio between increments and fellings. Nowadays, just as in the last century, growth increments are estimated by National Forest Inventories (NFI), usually comparing repeated field measures in permanent plots on the ground (Tomter et al. 2016; Chirici et al. 2020a).

Traditionally, NFI programs monitor the forested land based on field data collected over a permanent network of sample plots. Although these plots are visited repeatedly through time, there are large temporal gaps between remeasurements (e.g., 5–10 years) such that many forest canopy disturbances go undetected (Schroeder et al. 2014). For this reason, information on felling is instead estimated from different sources such as reports provided by the local authorities responsible for authorizing such activities. The shortcomings of these systems, which do not fully account for all wood removals at the national level, have been identified at least for Germany (Jochem et al. 2015), Canada (White et al. 2017), and Italy (Chirici et al. 2011). Forest fellings, especially those resulting from clearcutting, are clearly visible from satellite images and readily detected with automated algorithms, particularly when time series or multitemporal data are used, and when the spatial resolution of the data is compatible with the size of the harvested areas (Schroeder et al. 2007).

Dating back to 1972, Landsat data constitute the longest record of medium-resolution satellite imagery that have been demonstrated suitable for detecting forest harvesting and post-harvest vegetation recovery (Kennedy et al. 2010; Pflugmacher et al. 2012; Shimizu et al. 2016). The analysis of Landsat Time Series (LTS) increased dramatically with the change to a free and open access data policy (Woodcock et al. 2008), and applications of LTS to characterize forest disturbance and recovery have been regional (Griffiths et al. 2014), national (White et al. 2017), and global (Hansen and Loveland 2012). Moreover, the full potential of LTS for forest monitoring has been extensively reviewed (e.g., Frolking et al. 2009; Banskota et al. 2014).

Studies for the characterizations of regenerating forests with LTS analysis have been conducted primarily in temperate (Kennedy et al. 2007, Kennedy et al. 2012) and boreal forest environments (Olsson 2009), and only a few studies have focused on post-harvest recovery (Schroeder et al. 2007; White et al. 2017; White et al. 2018; White et al. 2019). In Schroeder et al. (2007), trajectories of Landsat spectral signals after clearcut in Douglas fir forests in Oregon (USA) were analyzed. The authors classified as “fast recovery” those classes where the spectral behavior was completely recovered after 12–14 years from the time of clearcut, also highlighting that the rates of forest regrowth after disturbance in western Oregon can be highly variable. Cohen et al. (2010) confirmed for the same region and forest type an average recovery time (time for recovering the same spectral behavior as before harvesting) of 10.6 years. White et al. (2017) analyzed the entire forested ecosystems of Canada for the period 1985–2010 and found that on average, it took 6.6 years (with σ of 3.9 years) for a forest harvested pixel to reach 80% of its pre-disturbance spectral behavior. In White et al. (2018), a similar analysis was carried out in boreal forests of southern Finland, with airborne laser scanning (ALS) data used to define recovery, determined as benchmark thresholds for canopy cover (> 10%) and tree height (> 5 m). White et al. (2018) found that < 5% of harvested forests recovered 100% of their pre-disturbance spectral properties in less than 10 years and achieved the ALS benchmark thresholds for cover and height. By comparison, 30% of pixels recovered in 10–17 years, and 28% recovered in > 17 years and achieved the ALS benchmarks. However, 27% of pixels did not recover 100% of their pre-disturbance spectral properties by the end of the time series, although these pixels did attain the ALS benchmark values for cover and height. When the authors applied an 80% spectral threshold instead of 100% threshold to determine spectral recovery, estimates of time to recovery were more realistic: 51% of pixels spectrally recovered in < 5 years and met the ALS benchmarks, with 36% of pixels recovering in 10–17 years, and 2% in > 17 years. Under such an 80% threshold scenario, no pixels were considered as non-recovered by the end of the time series.

This study originated from the idea that post-harvesting recovery time in Mediterranean forests can be consistently different from previous studies carried out in Oregon (Schroeder et al. 2007; Cohen et al. 2010), Canada (White et al. 2017) and southern Finland (White et al. 2018, 2019) as a result of differences in species composition, cultivation methods, and climatic and edaphic conditions. Especially in regions like Italy where the silvicultural system is coppice and vegetation recovery are dominated by fast vegetative resprouting from stumps. Typically, coppice forests in Italy are managed by retaining standard trees (Mariotti et al. 2017), that is, even-aged stands with 40–150 trees ha⁻¹ (the standards) of two to three times the rotation age, which are released at coppicing. Coppice felling is carried out by clearcutting at the end of the rotation period (usually 15–35 years), on areas ranging from a few square hundred square meters up to 10–20 ha, but clearcut size is most usually in the range of 1–5 ha (Bottalico et al. 2014).

The short rotation and the rapid asexual regeneration of trees that is typical of these coppice systems make them more challenging for monitoring with LTS data. Therefore, the objectives of this study were twofold. The first objective was to develop a reference database of clearcuts in Mediterranean coppice forests that can provide insights on the temporal intervals necessary for detecting and monitoring clearcutting in these forest types with automated approaches applied to LTS data. The second objective was to assess the utility of spectral recovery measures in rapidly regenerating coppice systems, in order to determine the capacity to monitor coppice regeneration with LTS data and supply information that is useful, to determine the ratio between increments and fellings and thus the sustainability of forest management practices. This is crucial to monitor the impact of climate on forest regeneration capacity.

2 Materials and methods

2.1 Study areas

Three study areas were selected in the Tuscany territory, each with an area of 225 km² within the Worldwide Referencing System (WRS-2) Path 192 Row 030, as shown in Fig. 1. The study areas were selected following three criteria: (i) areas without Scan Line Corrector (SLC) failure in Landsat 7 ETM+ since the 31st of May 2003 (USGS 2017), (ii) areas in different latitudinal zones, (iii) area in different elevation zones. So, the path of Landsat scene without SLC problem was divided in three strata on the basis of latitudinal and elevation zones to represent changes in forest types, and one squared study area per stratum was randomly positioned within the strata.

Fig. 1

Location of the three study areas within WRS-2 Path 192/ Row 030 Landsat scene in central Italy. In the quadrats are reported the detailed information of the three study areas: (i) clearcut maps obtained by photointerpretation and (ii) ALS coverage and (iii) the forest area in light green

Study area 1 is located in southern Tuscany with hilly terrain ranging between 0 and 650 m a.s.l. The area is dominated by Mediterranean vegetation with Pinus pinea, Pinus pinaster, and Quercus ilex on the coast and Castanea sativa and Quercus cerris in the inner part. Study area 2 is in central Tuscany, ranging between 100 and 1000 m a.s.l., and it is dominated by Quercus cerris, Quercus ilex, Ostrya carpinifolia, and Castanea sativa. Study area 3 is in the northern pre-Appennine area; the elevation ranges between 35 and 930 m a.s.l., and it is dominated by Quercus pubescens, Quercus cerris, and Castanea sativa. Six out of the 14 European Forest Types (Giannetti et al. 2018) are represented in the study area. More information about the orography and growing stock volume amount of Tuscany forests can be found in Chirici et al. (2020a). In all the three areas, the silvicultural system for broadleaf tree species is based on coppice with standards, managed with clearcut with reserves. In the study areas, trees are cut with a rotation period of approximately 18–20 years. Following the regional forest authority regulation, clearcutting can be carried out only during the winter season (1st of November–30th of April). The coppice system is not applied to coniferous forests, and clearcutting in areas dominated by coniferous species is followed by natural or artificial regeneration. Coniferous dominated stands in the study areas represent less than 6% of the total forest area, so the coppice system is the dominant forest management regime in these areas.

2.2 Forest mask

A fine-resolution forest mask of the three study areas representing circa 2013 conditions was generated from reclassifying a local land use/land cover map, available in vector format at a reference scale of 1:10,000 for the whole Tuscany Region. The forest mask has a minimum overall accuracy at 2013 with a minimum overall accuracy of the 85% (Regione Toscana 2013). In the three study areas, the forest mask was resampled with a pixel size of 30 m to exactly match the geometric resolution of Landsat scene. The forest mask was used as a reference to estimate the total forest area in the three study areas.

2.3 Forest harvesting database

To characterize the post-harvest spectral recovery of vegetation, we constructed a series of annual maps of clearcuts. These reference data were generated using visual interpretation of the LTS in an approach similar to that implemented in TimeSync for the same aim (Cohen et al. 2010). Photo interpreters delineated the spatial extent of each clearcut and recorded the year of harvest, adopting a minimum mapping unit of 0.1 ha. In addition to the LTS, the interpreters used high-resolution digital aerial photos available on-line via a Web Map Server service. First acquired in 1999, repeated acquisitions of high-resolution aerial photos of the study areas were available approximately every 3 years (Regione Toscana 2018). In the case of doubts during the photo interpretation phase, the clearcuts were visited in the field together with the local forest administration. As a result of the quality control implemented, we considered the resulting forest harvesting geodatabase to accurately represent the location and timing of clearcutting within our study sites (1999–2015). The resulting clearcut maps for each of the three study areas are shown in Fig. 1 (Chirici et al. 2020b).

The clearcut maps were intersected with the Regional Forest Inventory (RFI) data to determine dominant species (by area). The RFI field observations are available in squared units of 400 m × 400 m (Regione Toscana 2014).

2.4 Landsat time series data

For this study, we constructed a LTS of mostly cloud free images (cloud cover < 5%) from Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI images acquired in the period 1999–2015 (Table 1) and downloaded from the USGS web service https://earthexplorer.usgs.gov/. Manual interpretation was used to mask out areas with clouds and cloud shadows, which accounted for less than 1% of the total study areas. All images were acquired during the late spring/summer period and, to standardize as far as possible the spectral response of the vegetation, the data were transformed to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Masek et al. 2006).

Table 1 Time series of Landsat imagery used

Multiple spectral indices and metrics have been applied in the literature for monitoring post-disturbance vegetation recovery (Gitas et al. 2012; Chu and Guo 2014). Some studies have used multidimensional transformations of original bands such as Tasseled Cap Indices, principal component analysis, or Fourier transformation (Franklin 2001). On the basis of Pickell et al. (2016) and the applications presented in White et al. (2017, 2018), we selected the NDVI and NBR indices:

$$ \mathrm{NDVI}=\frac{\mathrm{NIR}-\mathrm{RED}}{\mathrm{NIR}+\mathrm{RED}} $$

$$ \mathrm{NBR}=\frac{\mathrm{NIR}-\mathrm{MIR}}{\mathrm{NIR}+\mathrm{MIR}} $$

2.5 ALS data

To characterize the post-harvest development of forest structure (e.g., height) and to evaluate the derived spectral recovery trajectory, we used open-access airborne laser scanning (ALS) data acquired by the Tuscany Regional Administration and by Italian Ministry of Environment (http://www.regione.toscana.it/geoscopio). The ALS data covering the three study areas were acquired in different years between 2005 and 2012 (Fig. 3). The ALS data were acquired under leaf-on conditions in the late spring or in the summer season (between May and August) using an ALTM OPTECH 3033 laser scanner. Data were acquired with a maximum scan angle of 20° and a flying height of approximately 2000 m. The point density of the different acquisitions ranged between 0.5 points/m² and 3 points/m². ALS-derived products were downloaded, including a 1-m digital terrain model (DTM) and a 1-m digital surface model (DSM). On the basis of DTM and DSM, we estimated the canopy height model (CHM) to characterize canopy heights above ground, by subtracting the DTM from the DSM. The CHM was used to provide estimates of stand-level heights, as the top of canopy height is known to be well estimated with ALS (R² of 85–90% between measured and estimated tree heights; Popescu et al. 2002).

2.6 Extraction of remote sensing data of forest harvested areas

To analyze the post-harvest recovery time of the vegetation, we normalized the year of disturbance estimated by the photointerpreter to time since disturbance (Borrelli et al. 2014). The estimated year of clearcut harvest was reclassified according to the time since disturbance, with negative numbers for years before the disturbance and positive numbers for the years after the disturbance.

On the basis of this information, for each clearcut polygon, we calculated the distance (in years) since the disturbance and calculated the stand-level average and standard deviation of the NDVI and NBR indices for all the images in the LTS stack, for each year.

For each clearcut polygon, we also extracted the median and maximum height of forest vegetation based on the ALS-derived CHM. We also calculated time since disturbance for the ALS data, by determining the number of years between the year of the disturbance and the year of the ALS acquisition.

2.7 Spectral trajectories analysis

For each polygon, we analyzed the spectral trajectories of NDVI and NBR to monitor changes in temporal spectral values related to clearcut and recovery of forest vegetation. We followed the approach of Kennedy et al. (2007) and used fitted temporal trajectories for our analyses. As reported by Cohen et al. (2018), the signal associated with spectral indices can be affected by temporal variations in spectral response as a function of different atmospheric conditions, vegetation phenology, sun angle variation, and sensor degradation. Therefore, trend-fitted values are preferred for analyzing spectral trajectories (Schroeder et al. 2007).

Kennedy et al. (Kennedy et al. 2007; Fig. 3) reported that there are four possible hypothesized temporal trajectories for forests:

• 1. Simple disturbance

• 2. Disturbance followed by exponential revegetation

• 3. Revegetation from disturbance prior to the observation record

• 4. Revegetation from prior disturbance, reaching a stable point during the observation period

For each of four hypothesized temporal trajectories, initial estimates of trajectory shape parameters are made, and then sent to a fitting function to adjust these initial parameters to find the best fitted trajectory (Kennedy et al. 2007). The discrepancy between the best fitted trajectory and the observed raw trajectory is summarized in terms of a standard F-statistic, and the probability of that F-statistic (p value) is calculated. The model with the lowest p value is selected. The model parameters describe key aspects of the disturbance regime, including year of disturbance, change in spectral value at disturbance (a proxy for intensity of disturbance), and rate of recovery of spectral values (a proxy for recovery rate). This process was repeated for each clearcut polygon. The fitted trend trajectories model represents appropriate functions to describe hypothesized trajectories, and the parameters describing those functions themselves capture the key characteristics of disturbance and regrowth. We refer to Kennedy et al. (2007) for more details on the trend-fitted trajectory models.

In our case, before fitting the models, we adopted a despiking approach similar to that of Kennedy et al. (2010), Bolton et al. (2015), and White et al. (2017) to avoid noisy observations derived by inter-annual phenology or atmospheric conditions, where noisy observations are detected by examining them in relation to their previous and subsequent spectral values using a smooth approach implemented in the R-cran packages oce (Kelley 2018). Then, the fitted trend models were obtained by searching for three breakpoints in the time series regression (i.e., where the coefficient estimates of linear regression models shifted from one stable regression relationship to a different one) by minimizing the residual sum of squares (Bai 1994) (Fig. 3). Hence, breakpoints may be defined as positions in time by which a significant shift in the time series values is detected. We adopted three breakpoints because for our purposes we needed to identify three different times in the LTS spectral trajectory fitted trend: (i) the year before clearcut (see Fig. 3), (ii) the year of the clearcut, and (iii) the point where the trajectory again reaches a stable condition after the clearcut.

Fitted trends maintained the capacity of expressing the original vegetation temporal dynamics and were well correlated with the original spectral reflectance values (0.92 ≤ R² ≤ 0.98 for both NBR and NDVI).

On the basis of fitted trend trajectories, for each of the polygons where the fitted models were associated to the hypothesis “disturbance followed by exponential revegetation,” we calculated the Years to Recovery (Y2R) metric following the approach of White et al. (2017, 2018). The Y2R indicates the number of years required for a pixel to attain 70%, 80%, 90%, and 100% (i.e. Y2R70%, Y2R80%, Y2R90%, Y2R100%) of its pre-disturbance fitted trend NDVI or NBR value. In our study, the pre-disturbance value used to define the Y2R was calculated as the average of the fitted NDVI or NBR values for the 2 years prior to disturbance, consistent with the approach applied in White et al. (2017).

For each recovery scenario (i.e., Y2R70%, Y2R80%, Y2R90%, Y2R100%), as proposed by White et al. (2018), we determined the CHM median height of the polygon when the NBR and NDVI trajectories reached the specified spectral recovery thresholds, in order to determine if the indication of spectral recovery corresponded to the minimum height value proposed by Bartels et al. (2016) (i.e., height > 5 m) for assessing forest recovery, and that also corresponded to the minimum value required to satisfy the FAO’s definition of a forest (FAO 2012).

3 Results

3.1 Trends in clearcut area

Overall, the three study areas cover 67,500 ha, of which almost 71% (47,745 ha) are forested. In the study period (1999–2015), a total of 9450 ha were clearcut, representing 14% of the total area and almost 20% of the total forest area, or approximately 1.2% of the forest area annually. As shown in Fig. 2, the most commonly logged forest types are turkey oak (Quercus cerris), holm oak (Quercus ilex), hop-hornbeam (Ostrya carpinifolia), and chestnut (Castanea sativa). Turkey oak was the most commonly harvested species, with more than 30% of the forest area dominated by this species during the study period. Turkey oak clearcuts represented 41% of the total logged area in study area 2, and 18% and 19% in study areas 1 and 3, respectively.

Fig. 2

In the main part of the figure the total forest area and the felled area collected as field truth within the three study areas, for the dominant tree species in the study period 1999–2015. In the smaller barplot, the felled area for each considered year for the four main dominant tree species

The temporal trend of clearcut area for the main species is reported in Fig. 2. No increasing or decreasing trends can be denoted, with the exception of a fluctuation for Quercus cerris with a period of 3–4 years.

3.2 NBR and NDVI spectral trajectories

Fitted and raw trajectories of the spectral signals for both NDVI and NBR (Fig. 3) clearly denote the estimated change due to the vegetation removal associated with the clearcut with a decrease, relative to average values for undisturbed forest areas, in NDVI of − 20% and a decrease in NBR of − 27%. All the mapped clearcut polygons, on the basis of the best fitted trend trajectories models, were associated with “disturbance followed by exponential revegetation”. In fact, the immediate vegetative regrowth from stumps in clearcut areas, which is typical of the coppice system, is evident in the year following clearcutting with a change of + 16% relative to the year of disturbance for NDVI and of + 20% for NBR. This can be seen in Fig. 4 where the first derivative of trends are plotted against time from disturbance. Of note, the temporal trajectories of NDVI and NBR have consistent trends with an R² = 0.922.

Fig. 3

Example of fitted trend models of spectral trajectories. The full line represents the temporal trend average value of the two vegetation indices NDVI and NBR directly measured on Landsat time series images. The dotted line represents the trend-fitted models obtained searching for three breakpoints (i.e., black dots) in the time series regression (i.e., where the coefficients shifted from one stable regression relationship to a different one) that identify the spike of clearcut disturbance

Fig. 4

Normalized differential percent annual changes of NDVI and NBR of clearcut areas depending on the time since logging. In the bottom right figure, the regression between the annual NBR values is displayed

The spectral trajectories of NDVI and NBR were similar, although NBR had a greater change magnitude (Figs. 3 and 4). While NBR has demonstrated capacity to characterize both short- and long-term recovery (Kennedy et al. 2012, White et al. 2017), NDVI is known to saturate rapidly making it better suited to shorter-term recovery assessments (Chu and Guo 2014; Pickell et al. 2016). In the context of coppice forest management presented herein, we observed no marked difference in the temporal trends or assessment of spectral recovery achieved using either the NDVI or NBR (Fig. 3). For rapid recovery, both indices were effective for assessing Y2R.

Differential NDVI and NBR trends were estimated for the different dominant species in the study area without revealing specific differences in spectral recovery by species. The average Y2R for each scenario was 0.51, 0.53, 0.57, and 0.59, for NDVI and 0.50, 0.53, 0.56, and 0.59 for NBR for the Y2R70%, Y2R80%, Y2R90%, and Y2R100%, respectively. With the Y2R80% threshold, 60% of clearcut areas had spectrally recovered after 1 year, 80% were recovered after 2 years, and 90% had recovered after 3 years. With the Y2R70% threshold, we found that 90% of the clearcut areas had likewise spectrally recovered within 3 years, while for theY2R90% and Y2R100% scenarios, we found that 100% of clearcuts had spectrally recovered after 4 years for both vegetation indices (Fig. 5).

Fig. 5

Percentage distribution of disturbed areas on the basis of their recovery spectral behavior. Trends calculated through the indicator Years To Recovery (Y2R) with different thresholds at 70%, 80%, 90%, and 100% of pre-disturbed conditions. The red line on the basis of NBR, in blue for NDVI. In the x-axis, the recovery time, while in the y-axis, the percentage of total area of clearcut that reached the indicated recovery threshold

Analyzing the trends of the Y2R indexes for the four different main tree species indicated that all species have similar trends in rapid spectral recovery. The Y2R80% indicated that almost 100% of the clearcut areas had spectrally recovered within 5 years, regardless of dominant species. Only Ostrya carpinifolia had a slower regrowth process in the first 2 years after the clearcut, compensated by a faster rate of recovery in the subsequent 3 years (Fig. 6).

Fig. 6

Same information as for Fig. 5 but here presented for the four main tree species in the study areas. Red line for NBR while in blue for NDVI. The x-axis gives the time needed for a clearcut to recover, while in the y-axis gives the percentage of total area of clearcut that reached to the indicated recovery threshold

3.3 ALS trajectories

In total, 1323 clearcut polygons were covered by the ALS data, representing 5076 ha or 51.2% of the total clearcut area. Normalized by time since disturbance, temporal trends in forest vegetation height in undisturbed forests were estimate from the CHM ALS data (Fig. 7). These trends relative to the year of clearcutting indicate that the estimated average maximum height was 19.3 m, ranging from a minimum of 16.4 m to a maximum of 22.4 m, and remaining relatively stable over the study period. In the year of clearcutting, there was a 90% decrease in the estimated median stand height (from 16.8 to 2.4 m), followed by regrowth at a rate of 0.94 m per year during the first 10 years post-clearcutting (Fig. 7). Temporal trends in the estimated average median height indicate rapid regrowth post-clearcut, we found that the estimated median CHM stand height after 2 years, 5 years, and 10 years was about 4.5 m (σ = 1.3 m), 8.3 m (σ = 1.7 m), and 12.1 m (σ = 1.2 m) respectively (Fig. 7), while we found that the estimated maximum heights in the clearcut areas remained relatively stable over time, likely as a function of the retention trees that remained within the clearcuts post-harvest (Fig. 7).

Fig. 7

Maximum height (red) and median height (blue) for all the mapped clearcuts covered by ALS data by years from clearcut. Whiskers represent the 1st and 99th percentiles of the CHM considered height

Moreover, we found that only 13% of clearcut areas that reached the recovery thresholds (both NDVI and NBR) in the scenario of Y2R70% also reached at the same time the threshold height of 5 m, while 51% of the clearcut area deemed recovered in the Y2R80% scenario had reached the threshold height of 5 m. By comparison, 87% and 93% of the clearcut area that was considered spectrally recovered according to the Y2R90% and Y2R100% scenarios respectively had also attained the 5 m height benchmark (Fig. 8).

Fig. 8

Box plots show the median CHM height of the pixel that reached the recovery threshold for each scenario: Y2R70%, Y2R80%, Y2R90%, Y2R100%. Whiskers represent the 1st and 99th percentiles of the CHM median height. The blue dashed line represents the recovery height threshold proposed by Bartels et al. (2016), and the thresholds correspond to the minimum values required to satisfy the FAO’s definition of forest (FAO 2012)

Of note, NDVI and NBR post-clearcut spectral recovery trajectories were strongly related with height trajectories derived from the ALS data (for NDVI: R² = 0.87 and for NBR R² = 0.72). Overall, we found that the use of a more conservative threshold (i.e., Y2R90% or Y2R100%) ensured that the majority of clearcut areas that were deemed to be spectrally recovered had exceeded the > 5-m height benchmark, with the median height in these areas of 8.3 m. Conversely, lower spectral recovery thresholds (i.e., Y2R70%) were attained more rapidly, but were not as strongly associated with a concomitant achievement of the specified height benchmark (i.e.  > 5 m).

4 Discussion

A large effort has been devoted in the last decades to develop methods for automatically detecting disturbances to forest ecosystems using remotely sensed data. In turn, the capacity to characterize the recovery of vegetation at these disturbance locations has been enabled by increased computational capacity, novel processing approaches, and free and open access to the Landsat archive (Woodcock et al. 2008). Recent research has concentrated on temperate and boreal forests demonstrating that in such conditions: (i) stand-replacing forest loggings can be accurately mapped by optical remote sensing because of their unique spectral and spatial characteristics relative to other disturbance types (Schroeder et al. 2011; Hermosilla et al. 2015b, 2016) and (ii) the number of years required to recover to pre-harvest spectral conditions ranges from 7 to 20 years (Cohen et al. 2010; Schroeder et al. 2007; White et al. 2017; White et al. 2018).

The results from three different study areas in Mediterranean conditions reported herein confirmed that forest harvesting is an important disturbance in these areas, with 21% of the forest area logged in the period 1999 to 2015. Both the NBR and NDVI were useful for estimating the spectral changes associated with clearcutting activities in these forests, with the NBR registering a larger magnitude change associated with the harvesting. This is consistent with previous studies that have found that the SWIR bands used in the NBR are better suited to relating changes in forest structure (Kennedy et al. 2010; Pickell et al. 2016). For the rapid recovery found within the coppice forest management system, we found no marked difference between the NDVI- and NBR-derived assessments of spectral recovery (Fig. 4). This result demonstrates that the NBR, which has been widely used for post-disturbance recovery assessments in other forest environments, is as effective for assessing rapid recovery in coppice systems as it is for assessing longer-term recovery in boreal and temperate forests.

With this experiment, we demonstrated that in Mediterranean conditions, when forest management is based on the coppice system, the post-harvest spectral behavior of these areas is different from that reported in previous studies in boreal and temperate forest conditions. In our study, the Y2R80% metric indicated that almost 80% of the clearcut areas were spectrally recovered within 2 years, compared to the ~ 11 years reported by White et al. (2017) in Canada and by White et al. (2018) in Finland for boreal and temperate forest conditions.

Such rapid spectral recovery in Mediterranean conditions may be attributed to two factors. First, clearcut areas in Tuscany must have reserve trees (the standards), whereby a minimum of 100 standards per ha is retained. The reserve trees in coppice forest continue the photosynthetic activities contributing to spectral trajectory regrowth, while in boreal countries, after a clearcut, recovery is either from seed already present in the soil or from nursery seedling that will be replanted within 2 years of the cut. Second, vegetative regeneration from stumps in the coppice systems is extremely fast, as confirmed by ground-based studies. For instance, in a chestnut coppice with standards, Manetti and Amorini (2012) found that the mean height of vegetative resprout was about 3 m, 8 m, and 10 m after 2 years, 6 years, and 10 years from clearcut, respectively. These data are confirmed by our results based on ALS metrics. In fact, we found that the median CHM height of clearcut areas after 2 years, 5 years, and 10 years was, respectively, approximately 4.5 m (σ = 1.3 m), 8.3 m (σ = 1.7 m), and 12.1 m (σ = 1.2 m) (Fig. 7). By comparison, growth rates in boreal forests are much slower, and that is reflected in longer recovery times (White et al. 2019).

These results are in keeping with expected growth rates associated with coppice forest systems, whereby vegetative regeneration from stumps accounts for the majority of regeneration in the stand (Müllerová et al. 2016). Unlike regeneration from seed, regeneration from stumps competes less with herbs and other vegetation that can quickly establish at a site. We observed that 50% of clearcut areas attained 80% of their pre-disturbance spectral value in 1 year, with a median height of 2.9 m (δ = 12.8 m).

These results demonstrate the potential of using LTS data for the automated detection of clearcut areas in coppice forests, and moreover, to monitor the rate and efficacy of post-clearcut regeneration. Most of the algorithms used for automated disturbance detection have been developed in boreal or temperate forest environments, where the persistence of the change event is sufficient to enable detection of disturbances (e.g., Kennedy et al. 2010; Huang et al. 2010a, b; Hermosilla et al. 2015a, b). But such approaches can be confounded by disturbance followed by rapid regeneration, particularly if temporal noise filtering approaches are aggressive, resulting in the misidentification of these events as noise. In fact, the characteristic 1-year spike of NBR and NDVI spectral trajectories that we found (Fig. 4) is similar to the spectral behavior that other studies may identify as noise spikes due to residual clouds, snow, smoke, or shadow in the analyzed images (Hermosilla et al. 2016; Kennedy et al. 2010). For example, Kennedy et al. (2010) affirmed that a real land cover change should be more persistent than a 1-year spike, because in cover changes, the post-spike spectral value cannot return to the pre-spike value so quickly. Our results provide instead a different interpretation of such rapid changes, which in our cases are due to the very fast post-harvesting recover of vegetation. Our findings suggest that the calibration of temporal noise filtering procedure optimized for temperate and boreal forest conditions cannot be used in Mediterranean coppice forest if the objective is to establish an automated system for disturbance monitoring (Masek et al. 2013), and to report the related trends in biomass (Pickell et al. 2014).

However, time series from Landsat eventually integrated by more recent missions such as Sentinel-2 offer a unique opportunity to provide a more complete, consistent overview of forest disturbances from natural or anthropogenic causes (White et al. 2017), especially because several authors have reported already that traditional forest logging statistics constructed from administrative records on felling authorization may not be consistent or spatially explicit, (e.g., Canada:Gillis and Leckie (1996); Italy:Chirici et al. (2011)). The integration of these data streams in a virtual constellation concept (Wulder et al. 2015) increases the density of the time series and, when coupled with increased computational capacity (Hansen and Loveland 2012), enables automated approaches to disturbance detection and attribution (e.g. Hermosilla et al. 2015a, b), even in complex European land use contexts (White et al. 2018; Senf et al. 2017). Our results demonstrated that clearcut areas in Mediterranean coppice systems can be very accurately predicted from LTS; thus, operational systems could be implemented to produce on an annual basis clearcut area statistics, which at the moment are not produced in Italy. LTS can also be used to monitor post-harvested vegetation recover over large areas in order to provide reliable information on the impact of climate change on the regeneration capacity of coppice forests.

5 Conclusions

Based on the analysis of a LTS acquired between 1999 and 2015 for three study areas located in Tuscany (Italy), the following conclusions can be drawn relative to the defined research objectives:

• Clearcutting in coppice systems can be readily detected and mapped via visual interpretation of LTS data in conjunction with high-resolution aerial imagery, thus confirming the results from Chirici et al. (2011) in a different Italian region. This work contributes to the development of a reference database of clearcuts in Mediterranean coppice forests that can provide insights on the temporal intervals necessary for detecting and monitoring clearcutting in these forest types with automated approaches applied to LTS data.

• The NBR index and associated recovery metric (the Years to Recovery or Y2R) is appropriate for studying spectral change and rapid recovery trajectories of forests in Mediterranean conditions.

• Clearcut areas in coppice forests have a very rapid post-harvest recovery with 80% of the disturbed areas spectrally recovering within an average of 2 years post-disturbance. These results are independent of the tree species composition and are confirmed by temporal analysis of ALS data.

• ALS data confirmed prior in situ observations relating to the vertical regrowth of coppiced trees, which grow at a rate of approximately 1 m per year in the first 10 years.

• Methods and approaches for monitoring post-disturbance recovery using automated LTS analysis systems developed for boreal and temperate forests have to be adapted to the very fast spectral recovery rates of Mediterranean coppice forests.