Assessing Canopy Responses to Thinnings for Sweet Chestnut Coppice with Time-Series Vegetation Indices Derived from Landsat-8 and Sentinel-2 Imagery

Abstract

Forest management treatments often translate into changes in forest structure. Understanding and assessing how forests react to these changes is key for forest managers to develop and follow sustainable practices. A strategy to remotely monitor the development of the canopy after thinning using satellite imagery time-series data is presented. The aim was to identify optimal remote sensing Vegetation Indices (VIs) to use as time-sensitive indicators of the early response of vegetation after the thinning of sweet chestnut (Castanea Sativa Mill.) coppice. For this, the changes produced at the canopy level by different thinning treatments and their evolution over time (2014–2019) were extracted from VI values corresponding to two trials involving 33 circular plots (r = 10 m). Plots were subjected to one of the following forest management treatments: Control with no intervention (2800–3300 stems ha⁻¹), Treatment 1, one thinning leaving a living stock density of 900–600 stems ha⁻¹ and Treatment 2, a more intensive thinning, leaving 400 stems ha⁻¹. Time series data from Landsat-8 and Sentinel-2 were collected to calculate values for different VIs. Canopy development was computed by comparing the area under curves (AUCs) of different VI time-series annually throughout the study period. Soil-Line VIs were compared to the Normalized Vegetation Index (NDVI) revealing that the Second Modified Chlorophyll Absorption Ratio Index (MCARI2) more clearly demonstrated canopy evolution tendencies over time than the NDVI. MCARI2 data from both L8 and S2 reflected how the influence of treatment on the canopy cover decreases over the years, providing significant differences in the thinning year and the year after. Metrics derived from the MCARI2 time-series also demonstrated the capacity of the canopy to recovery to pretreatment coverage levels. The AUC method generates a specific V-shaped time-signature, the vertex of which coincides with the thinning event and, as such, provides forest managers with another tool to assist decision making in the development of sustainable forest management strategies.

1. Introduction

Sustainable forest management plays a central role in achieving the 2030 climate targets that were endorsed by the European Council in 2014. It also secures the survival of forest ecosystems for the benefit of present and future generations, with managed forests providing functions that may be social, economic or ecological in nature. Many studies have addressed these different aspects in order to find strategies to preserve the forest for future generations, taking care of them and the resources they provide [1,2]. As a fundamental attribute of vegetation, canopy cover is an essential biophysical variable that is widely used to monitor forest ecosystems in terms of ecology, hydrology, carbon- and nutrient cycling and global change [3]. One standard quantitative tool used to assess foliage levels that helps to understand the connection between canopy structure and ecosystem functions is the Leaf Area Index (LAI) [4], and, indeed, it is defined as an Essential Climate Variable by the Global Climate Observing System [5].

LAI represents or quantifies the total one-sided green leaf area per unit of ground surface as a dimensionless quantity [6]. Its importance is due to the fact that it is related to driving ecosystem functions and the processes involved, such as primary production, water fluxes, energy exchange and carbon cycle [7,8,9]. LAI varies around the world depending on the different biome involved. However, at a regional scale, LAI variation can be impacted by environmental factors such as natural disturbances, site fertility, phenology, management practices, and interactions among all of these factors [10]. LAI field measurement methods have been investigated in depth through both direct and indirect methods (e.g., [11,12,13,14,15,16]). Direct methods involve estimating LAI from the destructive sampling of leaves by harvesting trees or through the use of litter traps [14,15,16]. The indirect methods can be divided into two different categories: allometric equations and optical methods. The former estimates LAI based on empirical regression with other forest variables, such as diameter at breast height (DBH) while the latter uses a logarithmic relationship based on gap fraction measurements or canopy transmittance [17]. Optical methods (e.g., Canopy Light Analyzer LAI-2000/2200/2200C or hemispherical cameras) have become the most commonly implemented in field measurements and a wide range of reviews have covered the devices available for LAI field measurements (e.g., [11,16]). Although direct estimates provide values that are closer to the true LAI, indirect measurements are less time consuming and labor intensive than direct methods and thus allow sampling at larger spatial scales [16,18]. One important limitation of these indirect techniques, however, is the light/sky conditions under which LAI needs to be measured. The ideal scenario is an overcast sky, otherwise errors in LAI calculation may occur [11,19]. Also, while indirect methods can be used to cover a larger area than direct methods, the reality is that the results can only be applied to the scale at which they were measured. This therefore limits the utility of these methods in terms of spatial coverage, and the results cannot be applied at the landscape level [20].

In this sense, free access to satellite images has solved the shortcomings of LAI field measurements allowing for the upscaling of results from the plot level to the landscape level [18], meaning that LAI estimates have been derived at different temporal and spatial scales as a result of the development of multispectral satellite sensors [21]. Many satellites create images from solar radiation reflected by the Earth at different wavelengths. Each satellite has its own set of spectral bands corresponding to different wavelength domains (i.e., blue, green, red) which is used to measure the spectral response of vegetation and record the level of “greenness” or absence of vegetation, which are then used as the basis for vegetation indices (VIs) calculations. VIs are quantitative measurements that maximize the responses of vegetation parameters and minimize the impact of atmospheric effects and soil background reflectance [22,23]. Furthermore, a wide variety of VIs can be computed with the information from the different spectral bands used by the satellites to create images, and VIs are currently the most common approach for analysing vegetation changes over different time scales [24].

VIs and LAI are inter-related in that LAI is also based on the reflectance characteristics of vegetation; in other words, linked to spectral measurements [25,26]. There is, consequently, increased international interest in mapping LAI using satellite sensors to determine trends or monitor vegetation behavior, over time, with minimum effort [19,27,28,29,30]. Using indirect methods and remote sensing data, some studies have estimated LAI using data from a particular moment in time, such as the summer months [9,26,30,31] while others consider the entire vegetation period [28]. A number of different VIs have been tested, among the most well-known being: Normalized Vegetation Index (NDVI), Simple Ratio (SR), Soil-Adjusted Vegetation Index (SAVI), Transformed Soil-Adjusted Vegetation Index (TSAVI), and Global Environment Monitoring Index (GEMI). Not only have VIs been used to evaluate a particular moment in time, they have also been used with time series data to improve understanding of how processes change or behave over time [32,33,34]. To understand the processes and drivers of disturbances in ecosystem functions, the first step is to detect changes within a time series. The ecosystem changes that are more commonly observed with remote sensing approaches can be classified into three groups [35]: (1) seasonal or cyclic changes, (2) gradual trend changes and (3) abrupt trend changes. The latter can be the result of natural causes such as a drought event or a wind throw, as well as by human disturbance, including forest management practices. Clearly it is important for forest managers to know how forests react to this latter type of change in order to develop sustainable forest management practices. One such forest treatment is thinning (i.e., the removal of a certain quantity of trees), used to improve spacing and concentrate the resources of a site on a selected set of trees in a forest stand in order to increase their growth rates. This obviously has a considerable impact on forest structure and one important application of LAI is its use as an indicator to evaluate such thinnings. That said, a remote sensing proxy for LAI that is more efficient and cost effective, such as VIs, would be a valuable tool in monitoring the evolution of managed forest stands. However, to our knowledge, there are no studies in the literature that focus on the sensitivity of VIs to the temporal evolution of the canopy following different forest management strategies, and thinning, in particular.

In various European regions, the broadleaf tree sweet chestnut (Castanea Sativa Mill.) has been subject to human management in terms of both fruit and timber production, specifically and traditionally coppicing, which in terms of timber production is the production of poles, vineyard-stakes, sticks, pickets, fuel wood, etc., [36]. However, this traditional management does not include silvicultural treatments such as thinning, and therefore it fails to exploit the full potential of the sweet chestnut, which requires this treatment in order to obtain quality sweet chestnut timber [37]. Consequently, to move from coppice management for ‘traditional’ to ‘quality’ wood, more studies that consider these aspects are needed. Another important aspect of this species is its high resprouting capacity, which, from a remote sensing point of view, can make it very challenging to model the canopy structure accurately in certain species where new sprouts have the same leaf typology.

In the northwest of Spain, although sweet chestnut is common, most traditional coppice stands have been abandoned or their rotation time has been increased considerably. Nevertheless, local authorities and various other stakeholders are becoming increasingly interested in the active and sustainable long-term management of historic chestnut coppice. As such, detailed information is needed at various levels of forest management. In fact, recently, dynamic tools for this species that are able to predict future rates of change in various stand characteristics of sweet chestnut, have been developed [38]. That said, there remains the need to support forest managers through the development of other new tools that help them to understand how forests evolve following different management treatments, in order to inform sustainable forest management strategies, and to facilitate the production of quality wood. Consequently, the general aim of this study was to remotely monitor the evolution of sweet chestnut coppice in the northwest of Spain in the years after thinning treatments. Specifically, this paper sought to: (1) study the effect produced by different thinning treatments and the evolution of forest stands over the years at the canopy level with Landsat-8 and Sentinel-2, (2) identify the optimal remote sensing Vegetation Indices (VIs) that reflect temporal sensitivity to the early response of vegetation after thinning treatments and (3) find a VI time-signature that allows past events in forest sweet chestnut stands to be identified.

2. Materials and Methods

2.1. Forest and Field Data

The research site is located in Redes Natural Park, in the central eastern part of the Principality of Asturias in northern Spain (Figure 1). The average annual temperature of the area is 10–11 °C, with an annual rainfall ranging from 1000 to 1350 mm. The two pilot trial sites are located between 600 and 700 metres above sea level and have a northern orientation and a high degree of slope (45–57%). Both trials were carried out in small stands of similar ages and dominant heights, but of varying densities due to the intrinsic heterogeneity that the sweet chestnut coppice generates when it sprouts. A total of three different treatments were analyzed: (1) Control (2800–3300 stems ha⁻¹), (2) Treatment 1 which consisted of a single thinning that left a stock living density of between 900 and 600 stems ha^−1, and (3) Treatment 2, which involved a more intensive thinning that left a living stock density of around 400 stems ha⁻¹. In Trial 1, three areas were established and each was allocated to either Control, Treatment 1 or Treatment 2 conditions. In Trial 2, the area was divided into two and one area was left as the control while the other was subject to Treatment 1. Forest inventory data (

Table 1. Summary of forest inventory data for Trial 1 and Trial 2.

	Trial 1								Trial 2
	Control		Treatment 1			Treatment 2			Control		Treatment 1
	BT/AT	AT	BT	AT	AT	BT	AT	AT	BT/AT	AT	BT	AT	AT
Year	2015	2019	2015	2015	2019	2015	2015	2019	2015	2019	2015	2015	2019
t	16	19	16	16	19	16	16	19	13	16	13	13	16
N	2843	2660	4484	933	905	3664	439	439	3338	3169	3756	622	545
G	21.5	22.9	21.2	10.6	11.5	29.8	5.8	7.04	13.9	15.4	21.1	4.7	5.4
Dm	9.5	10.2	9.0	11.9	12.6	9.8	12.7	14.0	7.4	7.9	8.3	10.5	11.8
V	92.5	116.3	-	63.9	71.0	-	36.8	45.4	52.5	66.5	-	19.8	26.5

BT: Before Thinning, AT: After Thinning, t: Stand age (years), N: Stems per hectare (stems ha⁻¹), G: Basal area (m² ha⁻¹), D_m: mean diameter at breast height per plot (cm), V: volume (m³ ha⁻¹).

) for the study sites were collected on three occasions: before and after the thinning treatments (after leaf fall) in winter 2015 (installation year) and in summer 2019. Figure 2 shows example images from areas subjected to either treatment or control conditions.

2.2. LAI Field Measurements

Field measurements of LAI were conducted in July 2019 in 33 circular plots (r = 10 m) (Figure 1). The total number of LAI plots per trial was dependent on the size of the trial site, with each plot centre being located using a differential GNSS (Global Navigation Satellite System) with submetric accuracy (a similar system and data collection procedure to that described in [39]).

The canopy measurements for LAI estimation were carried out in the field with a LAI-2200C plant canopy analyzer (

Table 2. Mean LAI_e field measurements for the different treatments in Trials 1 and 2.

	Trial 1			Trial 2
	Control	Treatment 1	Treatment 2	Control	Treatment 1
LAIe	2.78	1.57	1.56	2.75	2.11
n LAI	2	6	8	4	23

n LAI: number of LAI plots, LAIe: Effective leaf area index.

). LAI indirect measurement methods are widely used in similar studies, and there are several extensive review papers covering the use of this device for LAI field measurements (e.g., [4,11]). It is a portable instrument that directly provides the LAI value, which is referred to in the literature as Effective LAI (LAI_e) as it describes an indirect measurement of LAI that may include other parts of the vegetation besides leaves [6]. Hereafter LAI_e is used to refer to LAI measurements made with the LAI-2200C.

This instrument is composed of five concentric rings with central zenith angles of 7, 23, 38, 53 and 68 degrees in a “fish-eye” optical sensor, which measures light in the blue band (320–490 nm) under the canopy. Also, the light level that reaches the top of trees is simulated for calibration by measuring it in clearings without trees under the same circumstances (and close to where the other measurements are taken), so as to determine the ratio between the two transmittances (above and below the canopy). Then, the gap fraction is estimated from each sensor ring. Finally, LAI is computed through the inversion of the Poisson model comparing the transmittances from above and below the canopy (LI-COR Inc., 2015).

In this study, data acquisition was performed using an optical sensor equipped with a 90° view restrictor in accordance with the guidance of the LAI-2200C manual, and a control unit that was used to record readings above and below the canopy. All measurements were performed under uniform overcast skies to reduce the effect of scattered light in the canopy.

In the forest, LAI_e field measurements were made using the plot layout employed in similar previous studies within circular plots [40]. Briefly, in each plot, one measurement is collected at the plot centre, and four more are taken at each cardinal point, three metres from the centre. Due to differences in the canopy, to guarantee low standard error of LAI_e, supplementary measurements were collected at random points within the LAI plot to reach a total of 18 measurements per plot. Additionally, hemispherical photographs were taken within each plot using a digital camera (Nikon Coolpix 5400) equipped with a fish-eye lens (FC-E9, Nikon), which was held in place horizontally. All field measurements with both devices were made 1.5 m above ground in order to avoid ground vegetation influencing LAI values.

The sampling points with the camera were the same as those used for the LAI-2200C in order to rule out possible LAI_e points where there was a wide range of sweet chestnut sprouts, which might overestimate LAI_e. As a consequence, LAI_e was recomputed after field measurements by discarding the points that, in the hemispherical photographs, showed a high degree of influence of the sweet chestnut sprouts and those points were considered outliers.

2.3. Satellite Data

Google Earth Engine (GEE) is an online cloud platform of geospatial datasets at the planetary scale that provides access to historical satellite imagery, and it is continuously updated [41]. Landsat-8 (L8) and Sentinel-2 (S2), a paired satellite, are two medium resolution Earth Observation satellites that have become popular in the study of forest variables, including LAI_e estimates, due to them having better resolution than other satellites such as MODIS. L8 was launched in 2013 and has a spatial resolution of 30 m across nine bands with a revisit time of 16 days. S2 satellite, launched in 2015, has 13 spectral bands of 10–30 m spatial resolution, a five-day revisit frequency. For this work, a temporal imagery dataset from L8 available in GEE were used to analyze the study area from 2014 until the end of 2019. The thinning treatments were carried out at the end of 2015, after the full leaf fall of the sweet chestnut coppice, meaning that 2014 and 2015 are considered years without intervention and 2016 to 2019 are the years where the effect of treatments become visible in terms of canopy evolution. The L8 image collection selected was the atmospherically corrected surface reflectance data from the Landsat 8 OLI/TIRS sensors that feed into the USGS Landsat 8 Surface Reflectance Tier 1 with a precision and terrain correction of level-1TP (for a more detailed explanation of the data, see [42]). These images contain five visible and near-infrared (VNIR) bands, two short-wave infrared (SWIR) bands and two thermal infrared (TIR) bands. In order to have a suitable dataset for the time series, the images were visualized one by one in the GEE in order to select only those not affected by the presence of clouds or snow. As a result, there were considerable lapses in continuity of the data and, in fact, only 82 out of a total of 232 images were suitable for the research needs. However, the information is still valuable because the extensive time frame of the available L8 images allows for the comparison of the situation before and after thinning. S2 data, in contrast, is only available from 2016 onwards (i.e., after the study started). Despite this, S2 images from 2016 and 2019 were used in the study so as to compare them with those from L8 due to their better spatial resolution, as well as to enable valuable tools to be developed for both satellites to assist in sustainable forest management. The S2 image surface reflectance collection corresponds to the MultiSpectral Instrument, Level-1C from Copernicus program (for a more detailed explanation, see [43]). S2 images contain 13 spectral bands representing top of the atmosphere (TOA)reflectance. Again, all the images were visualized to select only those without snow and clouds (124 out of 524) to use for the analysis of the canopy evolution over time. Once the L8 and S2 datasets were debugged, different VIs were calculated. To do this, GEE information related to the spectral bands necessary to calculate each VI was extracted at the plot level (i.e., using the same 10-m radius plots that were established for the LAI_e field measurements) weighted in accordance with the percentage of the pixel involving canopy cover.

The Normalized Difference Vegetation Index (NDVI) is one of the most commonly applied VIs in forest studies and it is often related to LAI [44,45]. Despite its intensive use, however, NDVI has the drawback of being sensitive to the reflectivity of the soil in which the plant is rooted, which limits its potentiality as does the fact that it saturates in situations of dense canopy [45]. Therefore, to account for changes in soil optical properties and minimize the background effect, a group of Soil-Line VIs has been developed, such as the Soil-Adjusted Vegetation Index (SAVI), the Modified Soil-Adjusted Vegetation Index (MSAVI) and the Optimized Soil-Adjusted Vegetation Index (OSAVI), among others. VIs of this type facilitate improvements in canopy structure estimates and minimize the impact of the background and the atmosphere (e.g., [46,47,48]). The authors of [49] designed various new indices: the First Modified Triangular Vegetation Index (MTVI1), the Second Modified Triangular Vegetation Index (MTVI2), the First Modified Chlorophyll Absorption Ratio Index (MCARI1) and the Second Modified Chlorophyll Absorption Ratio Index (MCARI2), all of which are less sensitive to chlorophyll content variations and are linearly related to LAI_e, and also resistant to soil and atmosphere effects. Thus, in this study apart from NDVI, a group of Soil-Line VIs was also evaluated.

Table 3. Vegetation Indices included in this study and their formulations. R refers to reflectance and the sub-indices indicate the value of the central wavelength in nanometres.

Vegetation Index	Equation	Reference
Normalized Difference Vegetation Index	$NDVI=\frac{R_{800}-R_{670}}{R_{800}+R_{670}}$	[50]
Soil adjusted vegetation index	$SAVI=\frac{\left(1+0.5\right)\left(R_{800}-R_{670}\right)}{R_{800}+R_{670}+0.5}$	[46]
Modified Soil-Adjusted Vegetation Index	$MSAVI=\left(1+L\right)\frac{R_{800}-R_{670}}{R_{800}+R_{670}+L}$	[48]
Optimized Soil-Adjusted Vegetation Index	$OSAVI=\left(1+0.16\right)\frac{R_{800}-R_{670}}{R_{800}+R_{670}+0.16}$	[47]
First modified triangular vegetation index	$MTVI1=1.2\left[1.2\left(R_{800}-R_{550}\right)-2.5\left(R_{670}-R_{550}\right)\right]$	[49]
Second modified triangular vegetation index	$MTVI2=\frac{1.5[1.2\left(R_{800}-R_{550}\right)-2.5\left(R_{670}-R_{550}\right)}{\sqrt{{\left(2R_{800}+1\right)}^2-\left(6R_{800}-5\sqrt{R_{670}}\right)-0.5}}$	[49]
First Modified Chlorophyll Absorption Ratio Index	$MCARI1=1.2\left[2.5\left(R_{800}-R_{670}\right)-1.3\left(R_{800}-R_{550}\right)\right]$	[49]
Second Modified Chlorophyll Absorption Ratio Index	$MCARI2=\frac{1.5[2.5\left(R_{800}-R_{670}\right)-1.3\left(R_{800}-R_{550}\right)}{\sqrt{{\left(2R_{800}+1\right)}^2-\left(6R_{800}-5\sqrt{R_{670}}\right)-0.5}}$	[49]

shows the group of VIs that were chosen, because of their theoretical sensitivity to canopy structure.

In order to select the VIs, as proxies of LAI_e, that best describe the temporal sensitivity in the evolution of canopy cover over time after thinning treatments in the sweet chestnut coppice forest, a simple linear regression between LAI_e and the VIs from the most recent available date (Summer 2019) was used. The LAI_e-VI relationships were investigated using the coefficient of determination (R²), which shows the proportion of variation explained by the model, and the root mean square error (RMSE), which is expressed in the same units as the dependent variable and gives an idea of the mean error when using the model.

2.4. Comparing Canopy Evolution before and after the Treatments

To compare canopy development and its evolution over a number of years following thinning treatment, this study used the Area Under the Curve (AUC) (Figure 3) to develop a method which aims to detect historical changes in the evolution of forests.

In the AUC method, the first step is to generate a time sequence for the values of each VI for all the images available in the study period, as illustrated in Figure 3a. All of these VI values are joined to form a polyline that represents the dynamics of the canopy’s evolution across the period studied. Each treatment has its own line, and the data for the Control plots serve to enable the detection of changes and facilitate comparison with the treatment sites. As a deciduous species, the foliation period of sweet chestnut is from the end of April to mid-June and leaf-fall occurs in November. Consequently, all of the evaluated VIs tend to have higher values when there are leaves than when there are none, meaning that curves will have their highest point in June, and reach values close to zero in the winter months.

The second step is to calculate the AUC for the corresponding curve (Figure 3b.) using a time interval of each natural year covered by the study, i.e., 1st January to 31st December. The different AUC values (see Figure 3) for each treatment with respect to the control plots are then compared to detect the differences and changes in the evolution of the canopy over the years. Because the VIs used in this study have different ranges of variation, they were normalized, from 0 to 1, to facilitate comparisons. In fact, the AUC detects how the VIs change over time and the behavior of the different VIs with and without treatment.

Lastly, in the third and final step (Figure 3c), the differences in AUC (denoted by AUCD) between each treatment and the control is calculated. AUCD can mitigate some of the variations in the data, in our case foliage, in different years, in order to render different time series more easily comparable. In Figure 4, Year 1 (pale orange) shows the effect of having images missing in or close to winter, which leads to a potential increase in AUC values, Year 2 (green) reflects when there are missing images near summer, meaning that the AUC decreases, and Year 3 (blue) represents the situation where there is a lack of images in both winter and summer, resulting in compensations in AUC values. This method was applied to all of the LAI plots, treatment plots and control plots, for the period analyzed in this study, that is, from 2014 to 2019 for L8 and from 2016 to 2019 for S2.

In addition, to test the statistical significance of the differences between treatments identified through the graphical analysis, an analysis of variance (ANOVA) comparison between the average AUC data for each trial in each year and between years was carried out, where p-value < 0.05 was considered statistically significant.

3. Results

3.1. Regression Between VIs as Proxies of LAI_e

The performance of the Soil-Line VIs differs to that of NDVI (

Table 4. Relationships between LAI_e and the VIs estimated from L8 and S2 for Trial 1 and Trial 2. Second Modified Chlorophyll Absorption Ratio Index (MCARI2) (the best performance among the Soil-Line VIs) and NDVI are shown in bold.

	Trial1				Trial 2
	L8		S2		L8		S2
VI	R2	RMSE	R2	RMSE	R2	RMSE	R2	RMSE
MCARI2	0.74	0.20	0.72	0.21	0.61	0.29	0.75	0.24
MTVI2	0.74	0.20	0.72	0.22	0.61	0.29	0.75	0.24
OSAVI	0.69	0.22	0.70	0.22	0.59	0.29	0.77	0.23
SAVI	0.68	0.22	0.63	0.24	0.55	0.31	0.74	0.24
MTVI1	0.67	0.22	0.65	0.23	0.57	0.30	0.71	0.26
MCARI1	0.67	0.22	0.65	0.23	0.57	0.30	0.71	0.26
MSAVI	0.67	0.23	0.63	0.24	0.54	0.31	0.74	0.25
NDVI	0.24	0.34	0.75	0.20	0.26	0.40	0.80	0.22

R²: coefficient of determination and RMSE: root mean square error.

). In the case of the Soil-Line VIs for predicting LAI_e, with both L8 and S2 data, the LAI_e-VI linear models performed slightly better in Trial 1 (R² values ranging from 0.63 to 0.74) than in Trial 2 (R² values ranging from 0.54 to 0.75), but in all cases, more than 50% of the variation was explained. Furthermore, the indices MCARI2 and MTVI2 showed the highest coefficients of determination with LAI_e in both Trials. However, NDVI models also performed well with S2 data, explaining more than 70% of variation, while in contrast, with L8 data less than 30% of variation is explained.

3.2. Canopy Evolution before and after the Treatments

For both satellites (L8 and S2), the VIs selected as proxies for LAI_e were the MCARI2 (given that, as explained in Section 3.1, of all the Soil-Line VIs it had the best agreement with LAI_e) compared with traditional formulations such as the NDVI (Figure 5, Figure 6, Figure 7 and Figure 8).

Figure 5 and Figure 6 show the behavior of the VIs before and after each treatment for L8. Figure 7 and Figure 8 show the canopy regrowth after the thinning treatment for S2. In all these figures for both satellites (L8 and S2), image availability is indicated by small black dots along the bottom axis, while lack of images due to cloud or snow is denoted by white spaces. Also, in Figure 5, Figure 6, Figure 7 and Figure 8 the seasonality that characterizes the sweet chestnut is clearly visible. At the end of each year, the VIs have a value close to zero and then the following year the VI values increase until reaching their maximum value in the summer months, which is when the sweet chestnut reaches its peak in terms of LAI_e, indicating the appropriateness of using VIs as a proxy of LAI_e.

The thinning treatments were carried out after leaf fall at the end of 2015, so 2014 and 2015 provide pretreatment data. It can be seen in Figure 5 and Figure 6 that for this period of nonmanagement there is very little differentiation between the lines in the graphs relating to the VI data for the control and treatments, unlike in 2016 (the year immediately after the thinning) where the effects of the thinning treatments are clearly shown. In Figure 5, Figure 6, Figure 7 and Figure 8, plots with thinnings are shown in different colors: gray for Control plots, orange for Treatment 1 and blue for Treatment 2 (only in Trial 1).

For L8 in Trial 1 and 2, the VI data throughout the time series is shown in the top of Figure 5 and Figure 6, respectively. The MCARI2 values demonstrate the progressive decrease in the visible effects of the thinning in terms of canopy cover (in the satellite images) over time, particularly in the third and fourth year after the thinnings (2018 and 2019). However, at the end of the time series, the differentiation can still be seen but the difference is less than the previous years: the lines intermingle, and control and treatments have similar values. In the case of NDVI, however, the thinning effect can only be seen in the year after the thinning (2016), while in subsequent years, the thinning effect seems to have disappeared since the lines for the control (green lines) and the treatments (orange and purple lines) show similar values, particularly Control and Treatment 1, in both trials. In winter months, however, the graphic relating to NDVI shows some anomalies (winter peaks), as this is when the sweet chestnut loses its leaves, so the values of the VIs are the lowest of the year. For S2 the VI data (top of Figure 7 and Figure 8) results show the same tendency as described for L8 just after the thinning, and, once again, the behavior is different for MCARI2 and NDVI. In Trial 1, MCARI2 reflects how the influence of treatment decreases over the subsequent years, whereas for NDVI, the effect of the thinning is only visible in the summer months for the two subsequent years, but from 2017 this tendency changes and the differences are apparent in winter months. In Trial 2, both indices have the same behavior as in Trial 1. For both satellites the configuration of the lines for each treatment that describes the canopy behavior is impacted by image availability. For example, in 2019 for MCARI2 in the case of L8 (Figure 5 and Figure 6), the canopy behaves differently in comparison with the years before: 2019 is characterized by the truncation of the summer peak. However, this effect disappears in 2019 for S2 (Figure 7 and Figure 8).

Looking at the AUC data (central area of figure) in Figure 5, Figure 6, Figure 7 and Figure 8, the effect of the lack of images can be seen. For example, for L8 (Figure 5) in the year 2018 the lack of images over winter led to an increase in the AUC value. Conversely, in 2019 there was a lack of images during summer, and as a consequence, the AUC decreased. However, for S2, this phenomenon hardly ever occurs due to the greater availability of valid images.

The average AUC data reveal a complementary vision of the canopy behavior compared to that of the VI data without further processing. With AUC, the behavior of the VIs throughout the time series is always easily seen. In these AUCs, for L8 and S2 and in both trials, the greatest difference between treatments is in 2016. However, yet again, the various VIs behave differently. For both satellites (L8 and S2), using MCARI2 the AUC values once again showed how the effect of the thinning diminishes over time to the extent that AUC is similar for control and treatment plots at the end of the time series (2018 and 2019). However, the behavior of the treatments in 2019 is clearer in the S2 data (Figure 7 and Figure 8) due the greater availability of images previously mentioned. In fact, in Figure 7 Treatment 1 has a higher value for MCARI2 than for the other two treatments. In the case of NDVI, on the other hand, for L8 in Trial 1 for 2018 (Figure 5) both treatments demonstrate only slight differences in values, while in 2019 and in Trial 2, this difference was not evident at all (Figure 6). For S2, in Figure 7 NDVI again showed greater differences between the control plots and treatment plots than did MCARI2. The AUC shows that for Trial 1 after 2017, Treatment 2 has higher values than either Treatment 1 or the Control, the opposite of the results for MCARI2. For Trial 2 (Figure 8), on the other hand, the difference between the Control and Treatment 1 is maintained until the end of the time series data, the same as the MCARI2 data.

The results of the ANOVA conducted on the AUC data (

Table 5. ANOVA data comparing AUCs for MCARI2 and NDVI for both satellites (L8: Landsat 8; S2: Sentinel 2), for all treatments in both trials. X indicates where the difference between treatments (C: Control; T1: treatment 1; T2: treatment 2) is significant (p < 0.05), and the dash (-)., where it is not.

		Trial 1						Trial 2
		T1-C		T2-C		T1-T2		T1-C
		L8	S2	L8	S2	L8	S2	L8	S2
MCARI2	2016	X	X	X	X	X	X	X	X
	2017	-.	-.	X	X	X	X	X	X
	2018	-.	-.	X	X	-.	-.	-.	X
	2019	-.	-.	-.	-.	-.	-.	-.	X
NDVI	2016	-.	X	X	X	X	-.	X	X
	2017	X	-.	X	-.	-.	-.	X	X
	2018	X	-.	X	-.	-.	-.	-.	X
	2019	-	-	-	-	-	-	-	X

) provide a complementary evaluation of the differences between treatments. They reveal the same tendencies shown in Figure 5, Figure 6, Figure 7 and Figure 8. Differences between AUCs for the year 2016 were significant (p-value < 0.05) for both satellites and all VIs, except for NDVI (L8 and Trial 1). For MCARI2, in Trial 1 significant differences were maintained for longer between Treatment 2 and Control than between Treatment 1 and Control. In terms of satellite differences for Trial 1, both satellites behave the same, but in Trial 2, the significant difference between treatments for S2 is maintained until the end of the period analyzed.

Figure 5 and Figure 6 also show that, in the case of L8, the values of the VIs in the years before the thinning (2014 and 2015) are similar to when LAI_e field measurements were made (2019), indicating the recovery of the canopy to previous levels of cover. The AUCD (bottom of Figure 5, Figure 6, Figure 7 and Figure 8) shows the magnitude of the change in Control compared to Treatment plots through minimizing the impact of the lack of images and showing how the canopy develops over the years. For L8 (Figure 5), MCARI2 detects a greater difference between treatments 1 and 2 with respect to the control than NDVI (12.8% and 6.9%, respectively for Trial 1, and 21.09% and 7.55% respectively for Trial 2). This is a trend that continues until 2019 when the effect of treatments becomes barely discernible. However, this trend occurs faster with NDVI, where from 2018 there is no difference between Treatment 1 and Treatment 2 with respect to the Control. In terms of S2 data, the difference between treatments and control levels out in 2018 for both indices.

4. Discussion

The satellites Landsat-8 [51] and Sentinel-2A/B [52] are continuously acquiring and recording huge amounts of image data, with different revisit times, which is valuable for monitoring the changing earth surface. In this study, L8 and S2 satellite data were used for estimating different VIs for use as proxies of LAI_e in order to detect changes over time in canopy regrowth after thinning treatments.

In the north of Spain, only 35% of the L8 images and 30% of the S2 images available for the time period analyzed could be used, as the remainder were severely affected by cloud or snow. As a result, there are long stretches of time with no images. For instance, with the L8 data (as shown in Figure 5 and Figure 6) in 2019, the expected peak in summer is truncated by the absence of images. This limitation is, to a great extent, related to the revisit time of 16 days for this satellite. However, with S2 data, the problem of not being able to see the peak generated in 2019 did not happen (Figure 7 and Figure 8) because S2 has a shorter revisit period, meaning the number of usable images in the summer of 2019 was 44 for S2, in contrast to 15 for L8. As a result, in the study of the canopy evolution of the sweet chestnut, the time in the season when a lack of images occurs, has differing implications, i.e., it is less serious a problem in winter when the tree has no leaves than it is in summer, at the maximum point of the species’ ‘greenness’. These two satellites have been compared in a number of other studies, which have confirmed the capacity of S2 to provide comparable data with L8 data because of their use of similar spectral bands [26,30].

The best prediction models for Soil-Line VIs as proxies of LAI_e were obtained using S2 (R² 0.63–0.72 in Trial 1 and 0.71–0.75 in Trial 2), albeit that there are only small differences with the L8 data (R² 0.67–0.72 in Trial 1 and 0.54–0.64 in Trial 2). In contrast, NDVI showed good prediction for S2 (R² 0.75 and 0.80 in Trial 1 and 2, respectively), but not so good for L8 (R² 0.24 and 0.26). Among the tested VIs, the study findings here are in accordance with several other studies that have reported that introducing the properties of the underlying soil in a VI minimizes the effects of the soil background on the vegetation signal [45,48,49,53]. Although NDVI is also a strong indicator of forest greenness and biomass [54,55,56,57] and it is often related to LAI [44], it has some problems in estimating the vegetation cover. NDVI at moderate values of LAI has difficulties distinguishing canopy from green understory, which causes noncontrasting background reflectance [58,59], while, contrastingly, it lacks sensitivity to green in situations of low vegetation cover. For this reason, it can over- or underestimate vegetation cover [60]. However, although there is a wide range of studies focused on the behavior of different VIs, there is a lack of VI research data related to forest management and temporal sensitivity in the early response of vegetation after thinning.

The change caused by the thinning is clearly detectable using L8 data, which covers a two-year period before the thinning treatment as well as four years after it, while the available S2 time series only starts a few months before the thinning, so the direct changes produced cannot be analyzed from these data, although it does allow for the complete analysis of the evolution of the canopy after the event. The sensitivity of MCARI2 and NDVI to forest management has, however, shown behavioral differences in monitoring forest canopy changes. MCARI2 data from both L8 and S2 reflect the trend of how the influence of treatment on the canopy cover decreases over the years, to the point that in 2019, it shows similar values to those before the treatment. Soil line indices, such as MCARI2, considerably reduce the influence of soil, especially for homogeneous plant canopies [61] such as in our study site, where all of the new sprouts of the sweet chestnut have the same type of leaves as the big trees left after the thinning. In contrast, NDVI showed some irregularities. In L8 (Figure 5 and Figure 6), in winter when there are no leaves, the NDVI has some peaks for the time data series analyzed (i.e., it behaves as if there were some leaves remaining) and in S2 the differences between treatments can mainly be seen in the winter months rather than the summer, again simulating the presence of leaves in the trees. Sweet chestnut is a late foliation species so in spring or autumn the understory can influence the dynamics perceived by the NDVI. These anomalies support the previously mentioned issue of NDVI not performing as well as some Soil–Line VIs when there are low levels of vegetation canopy.

This work demonstrates that the AUC method for analyzing satellite imaging time series data with respect to annual or seasonal changes is useful for monitoring annual or seasonal changes in the canopy of sweet chestnut coppice in the climatic zone of the study area. Such events generate a time signature that indicates a past event in the forest stand. The change considered in the study (i.e., thinning) is identified in the graphical representation of the time series data by a V-shaped dip (e.g., Figure 5 and Figure 6 in 2016). The VIs tested in this study demonstrate that Soil-Line VIs (such as MCARI2) seem to be more explicative than NDVI, with MCARI2 tending to exhibit the V-shaped dip, while NDVI shows some anomalies (Figure 5, Figure 6, Figure 7 and Figure 8). In addition, MCARI2 clearly shows the tendency of LAI_e recovery in terms of canopy regrowth (comparing the values of the VI before and after the thinning for both trials using L8 data), but there is no clear pattern in the NDVI time signature. The ANOVA results support the differences found between treatments using MCARI2 and NDVI throughout the study, especially in 2016, the year after the thinnings. However, while analyzing statistical differences between treatments is essential, in this case study it is also important to remember that new sprouts influence VI values. Consequently, graphical analysis is also necessary, as in the case of the MCARI2 data, where, despite the ANOVA not being significant, the graph showed clear differences in forest canopy.

The AUCD helps to describe what happened in the evolution of each treatment with respect to the control plots. The same tendency for the MCARI2 to demonstrate apparently logical canopy evolution is once again shown in Figure 5, Figure 6, Figure 7 and Figure 8. Canopy cover seems to recover, and the difference in the AUCD lines between the treatment and the control becomes smaller over the period analyzed, especially in 2018 and 2019. In other words, LAI_e in these years is similar to that before the thinning. Also, in general the AUCD eliminates the part of the variation that is due to the availability or lack of availability of images, thereby stabilizing the value of the control, and the slope of AUCD becomes more homogeneous after the treatment.

LAI_e field measurements were made four years after the thinning. However, as shown in Figure 5 and Figure 6, the value of the VIs in this year (2019) is almost the same as before the thinning treatment, although the forest distribution is different. This finding was visually verified in the field, where it was seen that sweet chestnuts that had been cut down in the thinning operation resprouted and their growth filled the gaps between the trees between the trees that had been left standing after the thinning. Each thinning operation results in an increase in the diameter growth of the remaining trees, which is compensated for after a few years by the trend of growth decreasing as trees age [2]. In the cases in this study, in Trial 1 for Treatment 1 the diameter increase was almost 2.5 cm more than the Control, and for Treatment 2 almost 4 cm, while in Trial 2, the growth associated with Treatment 1 was almost 4 cm. The extent to which the remaining growing stock increases is largely determined by the intensity of thinning [62]. It is important to know which VIs are most effective in monitoring the impact of forest management treatments. Thinning interventions provide better quality roundwood at final harvesting [37], which can be used for long-term products [63], it is a costly operation, although using the most effective monitoring system can offset some of the costs of thinning operations. Also, [64] showed a positive influence on the quantity of carbon stored in long- and medium-term products, which store carbon for a long time and have more and better market possibilities than the short-term products that are usually obtained from unmanaged sweet chestnut sites. As such, this study presents a tool for the sustainable management of forests, which can assess the recovery of the sweet chestnut forest stand following thinning operations.

5. Conclusions

This study presents useful information for forest managers of sweet chestnut coppice by providing an approach for detecting and monitoring changes in this particular forest ecosystem. In this case, the forest change analyzed is directly related to forest management, i.e., thinning treatments. Two different thinning intensities were considered with respect to a control site in two trials using time series data from L8 and S2 satellites. Different VIs were used as indicators of canopy regrowth-rate for the monitoring of the impact of thinning treatments.

The L8 data allowed for the comparison of the performance of the VIs both before and after thinning, and the change brought about by the treatment is clearly observable: the impact is far greater in the years immediately following the treatment, then the effect becomes less, to the extent that, at the end of the time series, the difference between treatment and control plots is relatively small, with VI values being similar to those before the thinning treatment. In addition, the frequent revisit time of S2 helped to better describe the behavior of the canopy in periods of the year where there was a lack of images from L8 due to cloud or snow. The images from both satellites evidenced that the effect of the thinning treatment is clearly detectable the following year, but not so much four to five years later, once yearly images do not allow past thinning events to be detected, therefore stressing the need to use time series to identify and analyze the effects of such events.

Of the Soil-Line VIs considered, MCARI2 was found to be a better predictor and proxy of LAI_e than NDVI. The behavior of MCARI2 is logical given the species involved and treatments carried out, unlike NDVI, which showed some anomalies (for example, in winter it shows a difference in canopy cover between treatments even though this is not logical because the sweet chestnut has no leaves). MCARI2 data from both L8 and S2 reflect the trend for the influence of thinning on canopy cover to decrease over the years, with significant differences between Control and treatments being evident in both 2016 and 2017 (the degree of difference depending on the treatment and trial). This is to the point that in 2019, it shows similar values to those before the treatment. For NDVI, however, the thinning effect can only be seen in the year after the thinning (2016).

The results of this study reveal that the AUC method used for testing time-series of remote sensing data is a useful tool for monitoring changes in sweet chestnut coppice stands. It generates a specific V-shaped time-signature, the vertex of which coincides with the thinning event. MCARI2 detects a greater difference between treatments 1 and 2 with respect to the control. In fact, MCARI2 data tended to exhibit the logical V-shaped tendency, but in the case of NDVI, once again, some anomalies appeared. As such, this AUC method for characterizing past events such as the one analyzed in this study helps to monitor canopy development over the years, thus providing forest managers with another tool to assist decision making in the development of sustainable forest management strategies. In this sense, this tool, by monitoring time-series, allows forest managers to detect when thinning interventions took place in the past, by finding the V-shaped dip in the data curve. As a result, it will help forestry managers to develop strategies for sweet chestnut coppice on the basis of the timing of previous management interventions, information that, in many cases, would otherwise be impossible to ascertain.