SPATIOTEMPORAL DYNAMICS OF FOREST PLANTATION CLEARCUTTING AT LANDSCAPE LEVEL

This study is about clearcutting sizes, their spatial distribution, and how this affects the structure of the whole landscape when they are assessed over time. We introduced the concept of "extended clearcutting patches" (ECPs), which are formed aggregating adjacent patches harvested in consecutive years. The main purpose of the study was to evaluate the spatiotemporal pattern of clearcutting and ECPs at landscape level. The chosen landscape was the coastal commune of Constitución, Maule Region, in Chile, dominated by Pinus radiata D. Don plantations with a high land cover dynamic due harvesting activities. We used Landsat 5, 7 and 8 images (Surface Reflectance Tier 1) to produce yearly land cover maps from 1999 to 2016 and the approach proposed by Zhao et al. (2016) to classify them. We used 32,516 control points in the classification step obtaining a mean global accuracy of 0.94. The spatial distribution of the harvested patches is aggregated in single o accumulate years showing the long-term spatiotemporal pattern if this forestry practice. In average, ECPs have double the size than yearly harvested parches but they can be much bigger in same cases. This intensive forestry practices creates some spatiotemporal harvested complex, or extended clearcutting patches, that can have a bigger impact in functional aspects of the landscape and the wildlife present in this type of landscapes.


Clearcutting effects on forests
Forest commercial plantations are mainly used for the production of wood and biomass but, during the past two decades, there has been a growing interest in society and among timber companies in improving the role of these industrial forests in the conservation of biodiversity (Lindenmayer, Hobbs, 2004). Most of them are managed using intensive silvicultural schemes that involve a harvesting-replanting cycle. This process is usually carried out by clearcutting, removing all the trees in a given area, and then replanting it the following season (Gerding, 2009;Niklitschek, 2015). Clearcutting affects both the structure and composition of biodiversity, increasing the abundance of open space species until the second year after disturbance. Harvesting also causes the loss of structural similarity with the native forest stands, changes in microclimatic variables, increase of the edge effect (Pawson et al., 2006) and the decrease of permeability to the movement of forest species (Allen et al., 1995). Many other negative effects of clearcutting have been reported, including impacts on water yield (Sun, Li, 2005;Brown et al., 2005), nutrient cycles (Prescott, 1997), erosion (Iroumé et al., 2006) and promoting invasive species (Selmants, Knight, 2003).
Probably, the most controversial topic regarding the design of clearcutting operations is the size of individual cutblocks (Pawson et al., 2006). If a timber company requires to cut "X" ha in a given year, it can do so by harvesting just one cutblock of "X" ha or several smaller patches whose areas sum up to "X". Larger clearcutting sizes are more economical because more timber is removed at one time and fewer roads need to be built. Thus, normally the Net Present Value of plantations increases with average clearcutting size (Boston, Bettinger, 2001). On the other hand, larger clearcutting sizes produce undesirable effects such an increase in erosion and runoff levels (Iroumé et al., 2006). For that reason, in many countries limits are established for the maximum size of individual clearcuttings. For example, in the United States a single harvest area should not exceed 40.46 hectares (FSC-US, 2012). In terms of time dynamics, when using satellite images (i.e. Landsat-8 or Sentinel-2) to map this areas, a harvested stand can be sensed as a "harvested area", i.e. "none vegetation", for many years after the clearcutting due the time vegetation need to be sufficiently recovered to influence the spectral signature in each pixel. When two or more adjacent stands are harvested in consecutive years, the condition of "harvested area" can be applied to the whole area, including all adjacent harvested stands. In the scope of this study we call these areas "extended clearcutting patches" or ECPs. At landscape level, this patter of plantation clearcutting practice creates a mosaic of highly contrasting patches (dense forests vs recently harvested areas) that are habitat to different biological communities. Most wildlife living in commercial forests either die or are forced to relocate because of clearcutting operations (Escobar et al., 2015).
Since the appearance of satellite images, there have been numerous studies that search to obtain maps of land cover and quantify the changes over time. These include different scales, time periods, and frequencies (e.g Altamirano et al., 2013;Zhao et al., 2016;Gong et al., 2013;Hansen et al., 2000;Hansen et al., 2013). Nevertheless, few studies had included harvested areas as a different category (they are even considered as bare soil) or the time scale does not adjust to the spatiotemporal dynamics of this type of land coverage.

Clearcutting in Chilean Plantations
Chile has close to 3 million hectares of forest plantations with exotic species, where Pinus radiata D. Don covers the largest extent. These are concentrated in the coastal communes between the Maule and the Araucanía Region and are all harvested through clearcutting. Contrary to other countries, in Chile, the law only establishes that in the mentioned regions, harvested areas that exceed 500 ha/year should enter the environmental impact assessment system, where it is decided if it can be done or not without a formal limit value (FSC-Chile, 2010). Further references can be found in Niklitscheck (2015).
The purpose of this study was to evaluate the spatiotemporal pattern of forest plantations clearcutting and extended clearcutting patches at landscape level. In particular, we aimed to assess the sizes and spatial distribution of harvested and extended harvested areas.

Study Area
This study was carried out in the coastal community of Constitución, Maule Region, Chile ( Figure 1). This commune of 134,000 hectares, our focal landscape, is dominated by forest plantations, mostly monocultures of Pinus radiata covering approximately 70% of its surface (CONAF-CONAMA-BIRF, 1999).

Methods
For the study area, yearly land cover maps were obtained between 1999 and 2016, covering 18 years. We applied supervised classification using Google Earth Engine platform (Woodcock et al., 2008). Surface Reflectance Tier 1 of Landsat 5, 7 and 8 images were used for classification. Collection Tier 1 are atmospheric corrected and radiometrically calibrated, being suitable for time series analysis (USGS, 2018).
All image collections were filtered by date (January 1 to December 31 of the corresponding year), by location (polygon of the study area) and by percentage of cloudiness of the scene (<5%). 328 images were finally selected. Subsequently, to reduce the relief effect by applying a topographic correction using the approach proposed by Soenen et al., (2005). The SRTM Digital Elevation Model was used for this purpose.
Our approach follows the methods proposed by Zhao et al. (2016), developed for Chilean land cover classes which apply a multi-temporal approach combined with learning machine algorithms (i.e. random forest). For all classifications, we used the following variables as predictors:  Annual average reflectance of blue, green, red, NIR, SWIR1 and SWIR2 bands from Landsat images.  NDVI Maximum-value composites (Zhao et al. 2016).  Tasseled Cap brightness, greenness and wetness bands (Crist, Cicone, 1984;Huang et al., 2002;Baig et al., 2014).  Phenological components Phase and Amplitud from adjusted harmonic functions of NDVI annual time series (Shumwa, Stoffer, 2011).  Elevation, slope and aspect derivated from SRTM Digital Elevation Mode.
The following land cover classes were included: PL: Plantation, HV: harvested area, NF: native forest, SC: scrublands, WT: water, UB: urban, SA: sands and AG: agriculture lands. To assign control points to each class, every year, we used aerial photointerpretation on Google Earth. These yearly sets of sampling points were necessary because of the highly dynamic rate of land cover change of this type of forest landscapes dominated by commercial pine plantations.
The accuracy of classifications were assessed through a confusion matrix (Story, Congalton, 1986) which was built using ground control points, considering a minimum of 60 per class (Congalton, Green, 2007). Subsequently, using GIS environment, the harvested stands areas were obtained by calculating the change raster for each period, selecting only the pixels that changed from plantation to harvest. Meanwhile, extended clearcutting patches were obtained directly from the land cover maps.
Spatial analysis of the harvested stands areas and extended clearcutting patches and for each year was carried out. Metrics of number of patches, average patch area, Aggregation Index (AI) and Nearest Neighbor Index (NNI) between patches for the harvest class were obtain. Furthermore, Moran Index was calculated using periods from 2 up to 6 years. This procedure was done with Fragstat 4.2 (McGarigal et al., 2002).

Classification accuracy
The total number of control point used for the digital classification task was 32,516. Statistics by land cover class are presented in Table 1. Table 2 shows producer accuracies. The global mean accuracy was 0.94 and the minimum was 0.82 for scrubland class. The large number of control points we used for all classes every year explains these high accuracies. We wanted all classifications to be as good as possible because all further landscape metrics depend on them.

Extended Clearcutting Patches
Considering the results for land cover maps for every year, we spatially integrated two or more adjacent patches, classified as "Harvested" in previous years, into single ones that we so called "Extended Clearcutting Patches" (ECP).  The spatial distribution of the harvested patches is aggregated when considering a particular year but also when a longer period is considered (2-6), showing the long-term spatiotemporal pattern of this forestry practice. In Figure 3, the patch area of both harvested and extended harvested patches are shown. In average, ECPs have double the size than yearly harvested parches but they can be a much bigger in same cases ( Table 3). For both types of area, patches and ECPs, a temporal stability of the total amount of area they covered, can be easily observed I Figure 3. This means that, at landscape level, even with local changes in land cover types occurring every year, a spatial compensation in other parts of the landscape is occurring, obtaining a "structural metastability". The mean AI, which range from 0 (maximally disaggregated) to 100 (maximally disaggregated), of yearly harvest patches was 80.43. The NNI was < 1 for all years. The NNI measures the spatial distribution from 0 (clustered pattern) to 1 (randomly dispersed pattern) to 2.15 (regularly dispersed /uniform pattern). Additionally, Moran´s index was above cero in every year, confirming the clustered pattern. All statistics were statistically significant with p-value < 0.01.

Type of area
From the results, we can see that at landscape level, the total harvested area obtained directly from the land cover maps exhibited a soft variation over time reaching the highest values in 2000-2001 (14,119 and 13,113 ha, respectively) and again in 2010 (10,497 ha) (see Figure 3)

Functional implications at landscape level
Our results show that the aggregated composition of the forest landscapes, dominated by radiata pine, exhibited meta-stability, particularly in terms of the number of patches in land cover classes. This landscape "state" was introduce by Turner et al. (1993) under the concept of "landscape equilibrium". Considering the "shifting mosaic steady-state" concept proposed by Bormann & Likens (1979), the vegetation present in different stands on the landscape changes (i.e. adult plantation to harvested areas), but the proportion of them remains relatively constant at landscape level, this is, it is in equilibrium when considered over a large area or long time period ( Figure 5). Several studies have used meta-populations approach in dynamic habitats (Gu et al., 2002;Wimberly, 2006) and few have explicitly addressed how species respond to different yearly rates of habitat dynamics. Our results, indicate that the shifting mosaic steady-state model should not only consider the annual variation in land cover classes but also the accumulate effect of adjacent harvested in consecutive years. The spatiotemporal effects of these large time-accumulated harvested patches on wildlife remains mostly unknown.

CONCLUSIONS
The method approach generated reliable and highly accurate results. The spatial allocation of the harvest stand in a particular year is not random tending to be spatially adjacent to harvested areas from previous year. The forestry behaviour creates some spatiotemporal harvested complex, or extended clearcutting patches, that can have a bigger impact in functional aspects of the landscape, such as biological connectivity or even fragmentation. Nevertheless, the landscape exhibited temporal stability in terms of the total amount of area that is consider harvested, falling in the shifting-mosaic steady-state concept. Considering the former conclusion, new research question opens in order to analyse the relationship between clearcutting spatiotemporal schemes and the metastability of wildlife abundance in industrial forest landscapes, among other ecosystem functions.