Remote assessment of extracted volumes and greenhouse gases from tropical timber harvest

Timber harvest from tropical regions generates seven billion dollars annually in exports and is estimated to occur across 20% of the area of remaining tropical forests. This timber harvesting is estimated to account for more than one in eight of all greenhouse gas emissions from tropical forests. Yet there is currently no means to independently estimate extracted volumes and associated greenhouse gas emissions. In this study, we built upon an earlier paper that used an automated algorithm applied to LiDAR to accurately identify area of timber harvest impact in the categories of roads/decks, skid trails and gaps. This algorithm was applied to 2014 harvest areas in four concessions in Kalimantan, Indonesia. In two of these concessions, total harvested timber volumes and greenhouse gas emissions were measured and calculated in the field using data from 188 harvested and extracted trees. In order to relate remote sensing data with the estimated extracted volumes, we calculated factors that linked extracted timber volumes with greenhouse gas emissions, and applied three different regression equations. The parameters of the most accurate equation were the areas of roads, skid trails and gaps, explaining 87% of the variation in the data. For situations where rivers are used in place of roads for extracting timber and for instances of non-mechanized, often illegal logging, a second equation was created in which only skid trail and gap attribute data were used, and in this equation 86% of the variation was accounted. The final equation, intended for use in scenarios where LiDAR data are not available but moderate resolution imagery could be used, associated length of roads only with extracted volumes. In this case, 78% of the variation was explained. Application of the first equation permitted estimation of extracted volumes and associated greenhouse gas emissions from two additional logging concessions. We discuss the application of these equations to areas that have been identified as illegal logging concessions, and propose that these may be applied to larger regions across the country. These equations offer a way to estimate volumes of timber extraction when no ground data is available, and to calculate greenhouse gas emissions associated with extracted volumes, providing a simple methodology useful across forested tropical countries.


Introduction
Timber from tropical forests is an important resource, land use, and income generator, as well as a significant source of greenhouse gas emissions. The International Tropical Timber Organization (2015) estimated annual harvests of 270 million cubic meters from tropical forests, which generated seven billion US dollars in exports. Asner et al (2009a) estimated that during the 2000s, 2.98 million square kilometers of humid tropical forests were subject to timber removal, which amounts to just over 20% of the total forest area. Pearson et al (2017) identified timber harvest as responsible for 53% of forest degradation greenhouse gas emissions across the tropics, or 13% of the combined emissions from deforestation and forest degradation.
Yet independent, verifiable information on timber harvest is sparse. Tropical forest regions are extensive and often difficult to access, making data troublesome to obtain, and even more challenging to verify. National statistics as used by Pearson et al (2017) rely on completeness and accuracy, both from governments and from concessionaires reporting to governments. Such statistics necessarily omit illegal logging volumes. Illegal logging has been estimated to be about 40% of all tropical logging (Contreras-Hermosilla et al 2007), and as high as 61% in Indonesia, 65% in Ghana or 72% in the Brazilian Amazon (Lawson and MacFaul 2010).
Remote sensing based approaches for timber harvest monitoring faces challenges with regards to temporal frequency, spatial resolution, and sensor penetration. The timing of repeat observations is critical because treefall gaps and infrastructure associated with timber harvest are rapidly overgrown by surrounding canopy (Stone and Lefebvre 1998) and thus timber harvesting can easily be under reported. Spatial resolution is also an important criterion because sensors with coarser than 1 m resolution are limited to capturing only primary logging roads, rather than the logging gaps and skids trails where logs are being extracted (Asner et al 2004, Souza andBarreto 2000). Even high-resolution data are subject to limitations, specifically in canopy penetration; much of the damage from logging occurs sub-canopy, especially skid trails and thus optical imagery is not able to detect this impact. LiDAR imagery has the potential to capture such data, but extensive surveys using airborne LiDAR with the associated required temporal frequency becomes costprohibitive (Ellis et al 2016). Melendy et al (2018) developed an automated method to measure the extent of selective logging damage using airborne LiDAR data. The approach leads to cost savings through avoided manual classification of imagery, allowing efficient and automated assessment of large areas. Yet two gaps remain: (1) How can assessment occur beyond the path of LiDAR data collection? and (2) How can timber volumes and greenhouse gas emissions be associated with the extent of selective logging damage? Pearson et al (2014) developed a methodology that calculates emission factors for the extracted logs (ELEextracted log emissions), the damage created by logging gap formation (LDF-logging damage factor), and for the infrastructure (LIF-logging infrastructure factor), each correlated with the volume of the extracted log. There is a need for similar factors with comparably high confidence that could be associated with remotelysensed parameters such as area of canopy gaps, skid trails, and logging roads. This paper builds on the work of Melendy et al 2018 pairing LiDAR measurements with ground data. The approach allows for the development of new emissions factors from tropical timber harvest that subsequently can be extrapolated across entire LiDAR data collection areas. Our paper directly addresses: (1) the development of cost-effective methods that can assess greenhouse gas emissions from tropical timber harvest across broad areas, applicable to both mechanized and non-mechanized logging, and with the capability to estimate emissions from moderate resolution imagery; and (2) calculations of extracted volumes and greenhouse gas emissions due to logging across the 13 336 hectares of LiDAR imagery captured and analyzed by Melendy et al (2018) in Indonesia.

Methods
The design of the study included LiDAR imagery collection and analysis paired with field data collection. The field data allowed estimation of extracted volumes and greenhouse gas emissions that could be correlated with the LiDAR data, to allow extrapolation to the additional concessions, and ultimately across Kalimantan.

Study site
The study focused on four timber concessions in the 2014 cutting permit areas in Kalimantan, Indonesia, on the island of Borneo (figure 1, table 1).
Topography, accessibility, history and reduced impact logging (RIL) status varied between these concessions (table 1). Two of the timber concessions, Roda Mas and Timberdana, were visited for field data collection.

Ground data collection
Field data were collected between December 2014 and February 2015 in the Timberdana and Roda Mas concessions. Data collection followed the methods of Pearson et al (2014), to estimate per-gap extracted volume, and emissions from the extracted log, emissions incurred during treefall and emissions resulting from logging infrastructure (table 2).
A total of 188 felled timber trees were measured, plus 963 trees that were damaged incidentally, and additional measurements of infrastructure to be correlated with remote sensing measurement of logging infrastructure (table 3).
Reports on actual harvest statistics were collected directly from both Timberdana and Roda Mas (table  4).

LiDAR data collection
Details of the LiDAR data collection and analysis are given in Melendy et al (2018). In summary, airborne LiDAR data and high-resolution aerial photos were collected between October and December of 2014. The LiDAR data achieved complete coverage of a subset of the 2014 logging blocks of the concessions. Melendy et al 2018 developed and applied a method for automated classification across the imagery to determine separately the logging gaps, skid trails, and combined   Emissions directly resulting from logs extracted from the forest. Including milling waste and emissions from products subsequent to retirement.
Volume and density of extracted log. Assumed milling waste percentage and percentage of product remaining sequestered in excess of 100 years. Logging damage factor-LDF t CO 2 m −3 Dead biomass generated during tree felling. Includes both the top and stump of the timber tree and trees incidentally snapped or uprooted adjacent to the timber tree during the process of tree felling.
Top and stump calculated as the difference between the biomass of the extracted log and the total tree biomass (determined through allometry). Biomass of incidentally damaged trees calculated through allometry of observed damaged trees. Logging infrastructure factor-LIF t CO 2 m −3 Emissions resulting from infrastructure needed to extract logs from the forest. Includes roads, skid trails and decks where logs are stored prior to road hauling.
For areas with known timber extraction volumes, measurement of area and/or length of infrastructure with per area (or per length) emissions. Features classified by the algorithm as roads/decks and skid trail corresponded well with those same features identified in the field, with agreement measures ranging between 69% and 99% when adjusting for GPS location error. The automated algorithm attained the lowest agreement measures when detecting logging gaps, due in part to challenges associated with collecting high quality field GPS boundaries of the irregularly shaped gaps. Particularly as compared with the linear and systematic features created by mechanical means as usually seen with roads/decks and skid trails. Details of these comparisons are in Melendy et al (2018).

Development of extrapolation factors
Extrapolation factors were developed to allow the estimation of: -Extracted volumes from LiDAR measurements of gap, skid and road areas -Greenhouse gas emissions from reported timber volumes -Extracted volumes and greenhouse gas emissions from measurements of road length.
The relationship between extracted timber volumes and the areas identified as gap, skid trail, and road/deck damage in the Timberdana and Roda Mas concessions were explored with a multivariate regression analysis with a confidence level of 95%. The differences between actual concession volumes and predicted volumes using the regression were determined with the student t-test. The equation resulting from the significant linear regression was then used to estimate the extracted timber volume of the Suka Jaya Makmur and Sari Bumi Kusuma sites, the two concessions we had road/deck, skid trail, and gap area data (table 5), but no data on timber volumes extracted. We also developed an equation with a multivariate regression analysis for the relationships between extracted timber volumes and the gaps and skid trail damage (95% confidence level). The calculations were performed using the data from Timberdana and Roda Mas, with the intention of applying it to non-mechanized (often illegal) logging areas where no roads and decks are usually created.
The greenhouse gas emissions associated with the concessions' timber volumes were estimated using the ELE, LDF, and LIF factors (table 2). These factors were developed with the logging data from Roda Mas and Timberdana, and applied to the estimated timber volumes of Suka Jaya Makmur and Sari Bumi Kusuma.
The relationship between road lengths and total timber volume extracted in Timberdana and Roda Mas was determined with a univariate regression with a confidence level of 95%. To develop this regression, the concession blocks sharing a road were grouped, and their corresponding extracted timber volumes and road lengths were summed, resulting in three composite data points for Timberdana and three for Roda Mas.

Field data
From the 183 measured timber trees, the mean stump diameter was 100.1 cm in Roda Mas and 86.7 cm in Timberdana, while the mean log lengths were 19.1 m in Roda Mas and 17.3 m in Timberdana (table 6). Extracted logs, however, had basal diameters as large as 210 cm and lengths as long as 39.3 m. Extracted volumes were higher in Roda Mas than Timberdana, leading to higher extracted biomass and greater total damage. Skid trail and road widths did not vary greatly between the sites.

Emission factors
We calculated greenhouse gas emission factors for Timberdana and Roda Mas (table 7), following Pearson et al (2014). The LDF was the highest source of emissions in both concessions. The factors were not significantly different between sites (p-value > 0.05, F(2,5) = 18.513), so we combined factors to generalize and allow broader extrapolation.

Extrapolating estimates of extracted timber volume and emissions to additional concessions
Multivariate regressions were developed with the timber volume data (dependent variable) and damaged areas due to roads, skid trails, and gaps (independent Table 5. Areas and road lengths recorded across the concessions using the automated method of Melendy et al (2018 variables) for the Roda Mas and Timberdana concessions. Two equations were produced (table 8); one included roads, skid trails and gaps, and a second included just skid trails and gaps for use where timber is extracted using rivers and/or non-mechanized logging situations without road building. The equations estimate extracted timber volume in a given year associated with the areas of new roads, skid trail and gaps in that year.
Using the equation incorporating all factors (equation 1), we calculated the total volume of both concessions, and they were estimated to be 59 452 m 3 , similar to the actual extracted volume for both (60 827 m 3 ). Comparing predicted and actual volumes of both concessions individually, the estimates using the equation are not significantly different (p-value > 0.05) to the actual extracted volumes per block in the concession.
We used equation 1 to predict the timber volume extracted in the other two concessions for which we did not have timber volume data (Suka Jaya Makmur and Sari Bumi Kusuma), and used the combined Roda Mas-Timberdana emission factors (table 7) to estimate their emissions. The resulting timber volume and total concession emissions are summarized in table 9 and figure 2.
Examining emissions per hectare (figure 3), we find that Timberdana had the highest emissions rate of the four concessions, followed by Roda Mas. These four concessions had, on average, an emissions rate of 106 ± 11 t CO 2 e ha −1 , ranging from 42-151 t CO 2 e ha −1 .

Using road lengths to estimate extracted timber volumes
The length of the roads predicts 77.65% of the variability in extracted timber volume estimates (p-value < 0.001; F(5, 6) = 99.740; R 2 = 0.78; figure 4), with a linear fit following the equation (3): where V is the volume of extracted timber (m 3 ) in the concession blocks, and L is the road length (m) measured in the block in the LiDAR imagery for the 2014 harvest area. We used this equation to predict the volume of timber extracted and the combined Roda Mas-Timberdana emission factors (table 7) to estimate the emissions in the four concessions for the timber harvesting year assessed (table 10).

Discussion
Our aim was to develop and test an approach for associating timber volumes and greenhouse gas emissions resulting from timber extraction. Field measurements were collected from a total of 183 felled trees and 963 incidentally-damaged trees across two concessions. These measurements were paired with LiDAR data, and relationships were assessed and applied to LiDAR data that had been collected across 13 336 ha including two other concessions in Kalimantan, Indonesia. Correlation with impact areas (roads, skid Where: V is the volume of extracted timber (m3), R is the road and deck damage area (ha), S is the skid trail damage area (ha), and G is the gap damage area (ha).      trails and gaps) from LiDAR data produced an equation to estimate volumes that explained 87% of the variation in input data. This tight relationship was paired with remotely collected LiDAR data to estimate extraction volumes without the costs and logistical challenges of collecting ground reference data. A second relationship correlating with just skid trails and gaps would apply to concessions extracting timber using rivers, and to instances of non-mechanized (often illegal) logging. This relationship was almost as significant as the first, with the multivariate regression explaining 86% of the variation in input data. LiDAR data provides accuracy, but introduces logistical complexity and costs that might prevent widescale application. Therefore, we identified a relationship between road length and timber volume extraction, as road lengths can be determined from freely-available moderate resolution satellite imagery. The road length relationship with timber volume predicted 78% of the variability in input data.
Volumes were estimated for the two sites for which we had remote sensing data but no volume reports or ground data. Application of emission factors to the extracted timber volumes produced emissions estimates with narrow confidence intervals, with the half width of the intervals equal to 10.5% of the estimated value. It is important to note that this simplified relationship linking road length to volume extracted will not accurately capture activity from non-mechanized logging, which typically lack road infrastructure and rely on transportation via waterways. Pearson et al (2017) demonstrated the global importance of greenhouse gases from selective timber harvest, which accounted for 53% of forest degradation emissions and 13% of total forest emissions across tropical forest nations. Pearson et al (2014) developed a method that could be applied to production statistics from countries, but could not offer oversight independent of the statistics, and could not operate in situations where the statistics are unknown or are known to be poor. Ellis et al (2016) used LiDAR data across concessions in East Kalimantan to record logging impacts, both on canopy cover and on carbon stocks, but lacked volume estimates or emissions associated with extracted volumes, and thus do not provide a comparison.

Comparison with other studies
Many studies exist that use LiDAR to determine timber volumes (e.g. Naesset 1997, Clementel et al 2012, Kamaruzan Jusoft 2008, Latip et al 2013 especially in coniferous forest where strong correlations exist with tree height. Similarly, studies exist that determine forest biomass using LiDAR (e.g. Takagi et al 2015, Ioki et al 2014, Laurin et al 2016. These studies (particularly in non-coniferous forests) rely on extensive field work to derive new relationships with crown area, which is highly challenging to derive in natural closed forest. And none of these studies, nor any that we have found, examine extracted volumes or greenhouse gases as a result of timber harvest. Asner et al (2005) applied the Carnegie Landsat Analysis System (CLAS) to detect and quantify logging in the Brazilian Amazon. CLAS uses moderate resolution imagery (30 m resolution) with automated image analysis to estimate area of selective logging. Asner et al (2005) combined CLAS data with ground data and extraction statistics to estimate volumes of extracted timber and the corresponding logging emissions. However, while CLAS can collate imagery from multiple sources and thereby allows the simultaneous mapping of deforestation and forest degradation, it cannot differentiate between anthropogenic and natural forest change (Asner et al 2009b), and thus CLAS users still rely on ground data to verify emissions associated with a forest degradation agent. CLAS data relies on 30 m resolution optical data, which may be able to determine areas where logging is likely to have occurred, but does not offer the differentiation of actual logging damage that LiDAR does. CLAS is therefore, not being broadly used by countries or the private sector for assessment of forest degradation and logging impacts, likely due to the low resolution paired with excessive costs and complexities of accompanying field data. Our approach actually captures the three-dimensional measurements of gaps, skids trails and roads at a sub-1 m resolution, which makes the creation and application of extrapolation factors significantly stronger.

Potential uses of the methods
The availability of a cost-effective, simple and easily applied methodology has significant potential benefits to developing countries with selective timber harvest. Governments, private industry, and civil society can benefit from such methodology.
-For governments, the methods allow independent oversight of concessionaires to verify correct implementation of harvesting plans and application of taxation. In addition, data for governments on greenhouse gas emissions from timber harvest will be important for international commitments, including Nationally Determined Contributions and REDD+.
-Private industry can verify harvests undertaken and improve future planning, as well as demonstrate compliance with Governmental requirements for practices or participation in national programs.
-For civil society, an oversight function gives powerful additional information to assess how forests are managed, and to bring poor practice and corruption to light.
Understanding of illegal logging is an additional critical function of remote assessment methods. Where illegal logging is associated with concession overharvesting, the methods developed here are ideal. However, illegal logging may be dominated by timber extraction during illegal land conversion (Lawson and MacFaul 2010) and a different approach would be required. For high-value timber, illegal harvest may also occur on a selective basis (e.g. for mahogany; Lawson and MacFaul 2010). In those cases, it is unlikely that new roads are being built, and thus the remote sensing will only be able to capture logging gaps and (likely) minimally invasive logging trails; as such, a LiDAR analysis similar to that of Melendy et al (2018) would be necessary.

Estimating extracted timber volume and emissions from illegal logging
To test the method of Melendy et al (2018) for nonmechanized logging (illegal logging or legal situations where rivers are used and there are no roads), LiDAR imagery was captured over an area outside of legal concessions near Tarantang in Central Kalimantan (see square in figure 1). Over 36 ha, the Melendy et al (2018) method recorded 1.7 ha (4.7%) as skid trail and 1.1 ha (3.1%) as logging gaps. The method did not identify any roads in the area (figure 5), which lies outside any concession, and so the conclusion can be drawn that this is a site of illegal logging.
We applied equation 2 to the skid trail and gap areas of the illegal logging plot, obtaining an estimate of timber volume extraction of 730 m 3 , or 20 m 3 ha −1 . These timber volumes would be equivalent to 1623 ± 171 t CO 2 e, or 45.1 ± 0.5 t CO 2 e ha −1 , a rate about half of the average emissions rate of the four legal concessions. These lower numbers likely reflect the harvest of smaller trees, and the absence of roads and decks due to the logs being floated out.

Application of method across kalimantan
Global Forest Watch (www.globalforestwatch.org) lists a total of 270 concessions for 2012 in Kalimantan, with a total area of 12.5 million hectares, or about a quarter of the forested area in Kalimantan. The four concessions in this study, therefore, represent just 3.8% of the total concession area. For Sari Buma Kusuma, the 2014 logging area represents 1/35th of the total concession area, representing a 35 year rotation (which fits with Indonesian government requirements for natural forest management (Klassen 2005). Assuming the same across concessions in Kalimantan would give a total 2014 logging area of 358 thousand hectares. The mean road length per area, as determined by Melendy et al (2018), for the four concessions was 20.3 m ha −1 . Applied across the 358 thousand hectares, this is 7300 km. Using the equations developed here, that would represent an extracted volume of 8.6 million cubic meters, and an emission of 40.8 million t CO 2 . This is an overestimation, given the FAO (2009) estimate of just 5.6 million cubic meters of extraction from natural forests across Indonesia in 2006, likely due to the significant proportion of the assigned concession area not being actively harvested.

Conclusions
Data is difficult to obtain on the management of forests in many countries around the world. The fact that tropical forests are repositories of vast stores of carbon, are the sites of globally important biodiversity, a critical source of fresh water and resources for the development of nations, elevates the importance of gleaning information pertaining to these forests. Melendy et al (2018) have demonstrated that an automatic algorithm applied to LiDAR imagery is an accurate, effective and efficient means of assessing areas of timber harvest. This paper took the next step, providing approaches and equations that allow estimation of extraction volumes and greenhouse gas emissions. Resulting from this paper, an important research need will be to apply the methods developed here more broadly both with tests in additional concessions in Kalimantan and application to new regions; these tests should include application of the road-based method with validation with the LiDAR approach and/or field data and application of the LiDAR assessment of illegal logging with ground data validation.