Logging intensity drives variability in carbon stocks in lowland forests in Vietnam

Forest degradation in the tropics is generating large carbon (C) emissions. In tropical Asia, logging is the main driver of forest degradation. For effective implementation of REDD+ projects in logged forests in Southeast Asia, the impacts of logging on forest C stocks need to be assessed. Here, we assess C stocks in logged lowland forests in central Vietnam and explore correlations between logging intensity, soil, topography and living aboveground carbon (AGC) stocks. We present an approach to estimate historical logging intensities for the prevalent situation when complete records on logging history are unavailable. Landsat analysis and participatory mapping were used to quantify the density of historical disturbances, used as a proxy of logging intensities in the area. Carbon in AGC, dead wood, belowground carbon (BGC) and soil (SOC) was measured in twenty-four 0.25 ha plots that vary in logging intensity, and data on recent logging, soil properties, elevation and slope were also collected. Heavily logged forests stored only half the amount of AGC of stems ≥10 cm dbh as lightly logged forests, mainly due to a reduction in the number of large (≥60 cm dbh) trees. Carbon in AGC of small trees (5–10 cm dbh), dead wood and BGC comprised only small fractions of total C stocks, while SOC in the topsoil of 0–30 cm depth stored ~50% of total C stocks. Combining logging intensities with soil and topographic data showed that logging intensity was the main factor explaining the variability in AGC. Our research shows large reductions in AGC in medium and heavily logged forests. It highlights the critical importance of conserving big trees to maintain high forest C stocks and accounting for SOC in total C stock estimates.


Introduction
The world's forests play a crucial role in carbon (C) storage and regulating the climate (Bonan, 2008;Malhi and Grace, 2000). By deforesting and degrading forests, a substantial amount of carbon dioxide is released to the atmosphere, thereby contributing to global warming (Le Quéré et al., 2018). With increasing awareness of the importance of forests in mitigating climate change, a mechanism was developed in which financial incentives are provided for Reducing Emissions from Deforestation and forest Degradation, and promoting conservation, sustainable forest management and enhancement of forest C stocks (REDD+) in developing countries (Angelsen et al., 2009).
Forest degradation and disturbances account for 69% of total C losses from tropical forests (Baccini et al., 2017). Degradation can be broadly defined as a human-induced decrease in the capacity of an ecosystem to provide services (Putz and Redford, 2010;Sasaki and Putz, 2009; Thompson et al., 2013). Cross-continental analyses of 46 tropical and sub-tropical countries showed that degradation mainly results from timber extraction and logging, forming over 80% of all sources of forest degradation in Asia (Hosonuma et al., 2012).
Logging can have considerable negative impacts on the environment. Selective logging consists of harvesting the largest, and often oldest, trees of certain species. The felling and extraction of those large trunks generates substantial collateral damage to the surrounding stand (Sist et al., 1999;Sist and Nguyen-The, 2002), thereby affecting soil (Pinard et al., 2000), tree regeneration (Bagchi et al., 2011;Fredericksen and Mostacedo, 2000) and overall C stocks (Rozak et al., 2018). Yet, selectively logged forests can still contain substantial C and timber stocks, biodiversity and associated ecosystem functions (Putz et al., 2012;Edwards et al., 2014).
In several meta-analyses, the logging intensity, i.e. the number of harvested trees or volume of extracted wood per hectare, appeared to be the most important parameter in determining post-logging biomass stocks (Martin et al., 2015;Zhou et al., 2013) and biomass recovery (Piponiot et al., 2016;Rutishauser et al., 2015). Aboveground living carbon (AGC) typically declines with logging intensity (Gerwing, 2002;Martin et al., 2015;Rozak et al., 2018;Zhou et al., 2013). The amount of necromass and trees damaged by logging usually temporarily increases with logging intensity and can represent significant amounts of total C stocks in logged forests (Gerwing, 2002;Rozak et al., 2018;Sist et al., 1998). Logging intensity was also found to be the main factor explaining the variability in belowground living carbon (BGC) stocks (i.e. roots) (Rozak et al., 2018).
While most studies have focused on AGC, living C (above-and belowground) encompasses about half of total C stored in tropical forests (Pan et al., 2011). Other important C pools, such as necromass, litter, BGC and soil organic carbon (SOC) (IPCC, 2006;Pan et al., 2011), are often not included in forest C stock assessments. Only a few studies in Southeast Asia have investigated the effects of logging on C pools beyond AGC (Nam et al., 2018;Pfeifer et al., 2015;Rozak et al., 2018;Saner et al., 2012).
Vietnam experienced intense state logging after the end of the Vietnam War (referred to as American War in Vietnam; McElwee, 2004). Between 1943 and 1993 much of the country's forests were cleared, with forest cover declining from 43% to 28% (FCPF, 2018). Since then, forest cover has increased due to the expansion of plantations and natural forest regeneration. At the same time, remaining natural forests have become increasingly degraded through ongoing logging (Cochard et al., 2017;de Jong et al., 2006;Meyfroidt and Lambin, 2008). This shift from net deforestation to net reforestation in Vietnam in the 1990s is the so-called "forest transition" (Meyfroidt and Lambin, 2008).
Vietnam has pledged strong commitment to REDD+ (FCPF, 2018;Pham et al., 2012), requiring a sound understanding of how logging impacts forest C stocks. Several forest biomass and allometric model studies have been conducted in Vietnam (e.g. Con et al., 2013;Huy et al., 2016aHuy et al., , 2016bHuy et al., , 2016cKralicek et al., 2017), but only a few were carried out in logged forests (Hai et al., 2015;Luong et al., 2015;Nam et al., 2018Nam et al., , 2016. Successive logging bans implemented by the Vietnamese government means that since 2014 all logging in natural forests in Vietnam is illegal (FCPF, 2018;pers. comm. Trai Trong Le, 2019). Limited knowledge of the extent of logging and the impacts on forest C stocks in Vietnam hinders proper implementation of REDD+ activities, notably a reduction in forest degradation through better forest management.
In Vietnam and most tropical forest countries, detailed information on logging history is often lacking. Satellite data can be used to identify historical forest disturbances resulting from selective logging. Although this is challenging, recent work has shown it is possible with high to medium spatial resolution (5-30 m) satellite images (Miettinen et al., 2014). Images from the Landsat satellites have been widely used to identify forest degradation associated with selective logging (Asner et al., 2004;Hirschmugl et al., 2014;Langner et al., 2018;Souza et al., 2013). Although it is not possible to discriminate individually harvested trees using the 30 m resolution of Landsat images, the image resolution is sufficient to identify disturbances associated with selective logging such as logging roads, access tracks, skid trails and log decks (Asner et al., 2004).
Here, we assess C stocks in four main pools, namely AGC, aboveground necromass, BGC and SOC, in logged lowland forests in central Vietnam. We present a novel approach, which combines Landsat analysis and participatory mapping, to estimate historical logging intensities for the prevalent situation when records on logging history are unavailable. Finally, we explore possible correlations between past logging intensities, topography, soil and current estimates of AGC stocks.

Site description
Khe Nuoc Trong forest (KNT) is located in Le Thuy District in Quang Binh Province in north-central Vietnam. KNT covers approximately 20,000 ha of evergreen tropical forests and is part of a larger natural forest network of approximately 500,000 ha (Department for agriculture and rural development, 2010). Elevation in KNT ranges from 120 m to 1220 m, with the majority of the area (90%) < 700 m altitude on hilly terrain (Department for agriculture and rural development, 2010). Climate conditions are similar to the neighbouring Bac Huong Hoa Nature Reserve, where the area has a tropical monsoon climate with hot summers and relatively cold winters and storms from June to September (Mahood and Hung, 2008;pers. comm. Trai Trong Le, 2019). The hot season is from March to June, when the area receives a hot, dry wind originating from Laos ("Foehn"), with average temperatures of 29°C in June and July and frequent temperatures of 39°C when the Foehn is blowing (Mahood and Hung, 2008). In December and January, temperatures can drop to 15°C (> 500 m elevation). Total annual rainfall varies between 2400 and 2800 mm with highest rainfall from August to November. The drier period lasts only for a few months per year, i.e. between December and March (Mahood and Hung, 2008).
KNT is a key biodiversity area and one of the last remaining extensive lowland forests in the Annamese mountain range in Vietnam. A recent study showed that the largest area of deforestation and forest degradation between 2000 and 2010 was found in the north-central region of Vietnam (Khuc et al., 2018). KNT was selectively logged by the state between 1982 and 2007. Since 2007, the forest has been officially protected as a watershed protection forest and logging is forbidden, but widespread illegal logging still occurs (illegal logging also happened during the state logging period). The forest contains areas with various levels of forest degradation, from relatively undisturbed forests to heavily degraded forests. KNT has not been affected by the defoliants during the Vietnam War.

Plot establishment
Twenty-four permanent plots of 0.25 ha each (50 m × 50 m on the ground), 6 ha in total, were established in lowland forests (< 700 m elevation) between April 2016 and June 2017 (Fig. 1). As the poor, medium and rich lowland forests cover similar percentages of area in KNT, sample plots were equally and spatially stratified by the three forest types (n = 8 per class). Within each stratum, sampling was directed towards locations that cover the gradient of living aboveground biomass (AGB) values estimated by Baccini's pantropical biomass map (Baccini et al., 2012; access data via http://whrc.org/publicationsdata/datasets/pantropical-national-level-carbon-stock/) and constrained to accessibility. Within the selected pixels on Baccini's map, plots were randomized in the field.
A previous study in logged forests in Vietnam showed that AGC and BGC of woody stems 1 cm ≤ dbh < 10 cm stored 1.3% of total C stocks, while shrubs (including seedlings with stem < 1 cm dbh) stored only 0.2% of total C stocks (Nam et al., 2018). In logged forests in Vietnam and Indonesia, litter comprised only 0.2-1.6% of total C stocks (Nam et al., 2018;Rozak et al., 2018) and did not differ with logging intensity (Rozak et al., 2018). For these reasons, we focus on the dominant C stocks: AGC in woody stems ≥5 cm dbh, C in dead wood, BGC and SOC.

Living stems
Measurements in the plots were conducted according to the RAINFOR and GEM protocols (Marthews et al., 2014;Phillips et al., 2016; http://www.rainfor.org/en/manuals). In each 0.25 ha plot, living woody stems (trees, shrubs and stranglers, but excluding lianas) ≥10 cm diameter at breast height (dbh; i.e. diameter at 130 cm height or above buttresses or stem deformities, but below major branching points) were tagged and their dbh measured. A botanist identified each stem with Vietnamese and scientific species names in the field and, when needed, collected botanical vouchers for further identification in the herbarium of FIPI in Hanoi, Vietnam. Total living tree heights were measured with a Haga altimeter. Heights of deviating stems, i.e. broken, fallen, leaning, bended, rotten, hollow, resprouted, forked ≤1.3 m height, fluted stems and stranglers, were measured for each individual separately. For trees without these deviations, i.e. "intact" trees, heights were measured from a subset of trees to construct a local height-diameter (HD) model. For each plot, we randomly selected~10 intact trees in each of the following diameter classes (if this amount of trees was present): 5-10 cm, 10-20 cm, 20-30 cm, 30-40 cm, ≥40 cm dbh and the biggest stems. Living saplings, i.e. stems 5-10 cm dbh, were measured in a belt transect of 4 m × 50 m, running from north to south through the centre of each plot. The dbh, height and species were measured and identified as described above.

Dead wood
All fallen and standing aboveground dead wood, i.e. woody stems and branches ≥10 cm dbh (including logged stumps), were measured in each 0.25 ha plot. For each fallen piece (often referred to as coarse woody debris; CWD), the diameter at both extremities, length, shape, a visual estimation of % wood (subtracting hollow spaces) and degree of wood decay was recorded (Baker and Chao, 2011). For fallen pieces that were too heavy to lift to measure the diameter with a diametertape around the stem or branch, the length of two perpendicular sides was recorded. The degree of wood decay was classified in the field into one of five decomposition classes, based on simple wood characteristics (Table A.1; Baker and Chao, 2011). Measurements were only made to the point where the piece was at least 10 cm diameter and length measurements were made to the plot borders. For standing dead wood, the same measurements were done.

Soil collection
In each plot, soil samples were collected for 1) chemical analyses and 2) bulk density, both in two soil layers: 0-10 cm depth and S. M. Stas, et al. Forest Ecology and Management 460 (2020) 117863 10-30 cm depth. Each sample used for chemical analyses was collected with a soil auger from five sampling points within the plot and equal amounts were mixed to obtain a representative plot sample of at least 600 g. In a few cases, the soil was too stony to collect the samples with a soil auger and samples were collected by digging a hole and taking equal amounts of soil over the total depth. Bulk density samples were collected close to the centre of each plot by placing a core of 100 cc volume vertically in the 0-10 cm and 10-30 cm depth layers. Within a few days after collection, samples for the chemical analyses were airdried. The soil samples were analysed for the following properties: pH, SOC, total N, total P, available P, total K 2 O, exchangeable bases (K + , Na + , Ca 2+ , Mg 2+ ), cation exchange capacity (CEC), exchangeable H + and Al 3+ , texture and bulk density. Methods of the soil laboratory analyses are described in the supplementary material.

Topography
For each plot, the mean elevation was calculated from the four plot corners and the centre, recorded with a GPS device. The 50 m × 50 m plot was divided into subplots of 10 m × 10 m and the slope of each 10 m side of a subplot was measured. Subsequently, the maximum absolute slope of the north-south and east-west measurement was calculated for each subplot. The mean slope of the plot was calculated by averaging the maximum absolute slope values of the subplots.

Logging assessment
State logging activities started in 1982 and observed AGC stocks are the result of past disturbances and subsequent forest recovery. Turnover rates in tropical forests are high and mean decomposition rates of coarse woody debris in tropical forests vary between 0.11 and 0.46 year −1 (Baker et al., 2007), suggesting that evidences of historical logging could rapidly vanish in the field. To understand better historical logging activities at each site, we interviewed (ex-)employees working for the state logging company and contacted responsible government offices. Unfortunately, no records of historical logging or maps were available anymore for our sites, making it impossible to gather spatial information about logging intensities in KNT at plot-level. Therefore, we assessed logging in three ways: 1) measuring logged stumps in the plots to assess recent logging levels; 2) analysis of satellite images to evaluate the spatial pattern of historical forest disturbances associated with selective logging; and 3) participatory mapping to verify historical logging activities.

Logged stumps
Logged stumps ≥10 cm diameter were measured in each 0.25 ha plot. The diameter of the stumps was measured at the highest point (when buttresses were present, the inner round part of the stem was measured to estimate the diameter).

Remote sensing analysis
Previous studies have identified forest disturbances using different vegetation indices, including the Normalized Difference Vegetation Index (NDVI) and the Normalized Burn Ratio (NBR). NBR was first used to identify and map forest fires, but has recently been used for the identification of tropical forest canopy disturbances caused by selective logging (Grogan et al., 2015;Langner et al., 2018;Lima et al., 2019;Shimizu et al., 2017). NBR values vary between −1 and +1, with larger positive values indicating closed canopy cover and smaller positive or negative values for open canopy (Langner et al., 2018). NBR is particularly suitable for detecting forest disturbance, due to its sensitivity to bare soil and non-photosynthetic vegetation created by forest disturbance events, and the lower contamination from atmospheric haze compared to NDVI (Grogan et al., 2015;Langner et al., 2018).
We used the Earth Explorer and the USGS ESPA ordering system to retrieve all available Landsat images for Path/Row 126/48 and 125/48 with less than 70% land cloud cover between 1988 and 2016 inclusive.
A total of 439 Landsat 4/5, 287 Landsat 7 and 65 Landsat 8 surface reflectance images were obtained. Low quality pixels or those covered in cloud were removed from the analysis via the supplied 'pixel_qa' layer. We calculated NBR as: where NIR is the Near Infrared band (Band 4 in Landsat 4, 5 and 7 and Band 5 in Landsat 8) and SWIR is the Landsat Shortwave Infrared band (Band 7). We identified disturbance as pixels exhibiting a reduction in NBR, with two consecutive images showing NBR < 0.5. In this way, we were able to distinguish between recently disturbed canopy cover and naturally open canopy (Langner et al., 2018). Pixels without tree cover throughout the study period (i.e. rivers) and pixels with plantation forestry (identified through large-scale disturbance at 3-7 year intervals) were removed from the analysis. The Ho Chi Minh highway in the western part was not counted as disturbance (the road was not developed as a logging road). Previous work found that the area of skid trails (Broadbent et al., 2006) or the area damaged by logging (i.e. skid trails and tree fall gaps; Panfil and Gullison, 1998) followed a quadratic increasing function of harvest intensity. Therefore, we assumed a positive relationship between the density of disturbances and the logging intensity, and used the disturbance intensity as a proxy of logging intensity. Previous studies have suggested that the majority of timber extraction occurs within 700-1000 m of a logging road (Bryan et al., 2013;Gaveau et al., 2014). The logging intensity per plot was estimated by summing the number of disturbed pixels within 1000 m of each plot centre, with disturbed pixels within 500 m weighted by a factor 2. If after sufficient reduction in the NBR of a pixel (i.e. a disturbance), the NBR recovered and reduced significantly again in a different year, this was counted as another disturbance. Plots 1-18 were established in 2016 and disturbances were summed from 1988 to 2015. Plots 19-24 were established in 2017 and disturbances were summed from 1988 to 2016.

Participatory mapping
To further assess the intensity of logging across different regions of KNT, and confirm that the identified pixels with disturbance on satellite images represent mainly logging events, interviews were conducted in Kim Thuy commune (Le Thuy District, Quang Binh Province). Interviews were conducted with commune and village leaders as well as with households (both men and women) within two villages in March 2018. During the interviews, participants developed participatory maps identifying the level of logging prior to 2010 within different regions of KNT (logging activities after 2010 were likely still evident in the field by the presence of logged stumps).
First, maps of KNT showing key landmarks (rivers, roads, villages) were discussed and the participants familiarised themselves with the map. After discussion and agreement amongst the participants, each area was allocated a logging classification (light, medium, heavy) and a period when the majority of the logging occurred. After the interviews, the individual maps were combined to create a single map of KNT. Consistent responses from households within both villages and from the commune and village leaders provide confidence in the assigned logging classifications within the broad regions. The forest plots were then allocated a logging classification based on the map.

Estimation of forest C stocks
The BIOMASS package in R (Réjou-Méchain et al., 2017) was used to develop a local HD model, extract wood density (wood specific gravity; wsg) values and calculate plot-level living AGB. For the subset of intact trees from which heights were measured (n = 1065), several HD-models were fitted and the model with lowest residual standard error (a Weibull function) was chosen. As the linear relationship S.M. Stas, et al. Forest Ecology and Management 460 (2020) 117863 between logarithmic height and logarithmic diameter did not differ significantly with logging class, the intact trees from all plots were grouped to develop one HD-model. This model was applied to estimate heights of intact trees without height measurements. Heights of deviating trees were assigned from the individual height measurements in the field. Wood densities were assigned using the global wood density database Zanne et al., 2009). If the species was present in the database, the global species-level average was assigned. Wood densities of missing species were attributed by, in subsequent order, the global genus-or plot-level average. Tree AGB was computed using the generic allometric equation 4 of Chave et al. (2014), including dbh, height and wsg as input parameters (Réjou-Méchain et al., 2017), and summed per plot to retrieve plot-level AGB, expressed in Mg ha −1 . The volume of standing dead trees and fallen pieces of dead wood was calculated using the formula of a conical frustum. Aboveground necromass was calculated by multiplying these dead wood volumes with the estimated wood density for each decomposition class and summed per plot. For each plot, the wood density per decomposition class was estimated by multiplying the mean plot-level wood density of living trees ≥10 cm dbh with a proportional reduction in wood density that was estimated by Chao et al. (2017) for tropical lowland forests in Taiwan (see Tables A.1

and A.2).
Living below ground biomass (BGB), i.e. root biomass, for stems ≥10 cm dbh was estimated using the allometric equation of Nam et al. (2016), developed for trees > 10 cm dbh in logged forests in central Vietnam and using dbh and wsg as input parameters. SOC (in Mg C ha −1 ) was calculated by multiplying the fractional SOC with the bulk density and soil volume of the depth layer. Living AGB, necromass and living BGB were converted into C stocks using a default value of 47% C content (IPCC, 2006). Finally, living AGC ≥5 cm dbh, C in dead wood, living BGC and SOC were summed by sample plot to retrieve total C stocks.

Testing differences in forest structure between logging classes
Plots were ordered by the number of disturbances identified on the Landsat images and equally divided over three logging classes (n = 8 plots per logging class): light, medium and heavy logging intensity. The logging classes derived from the Landsat analysis and participatory mapping were compared. In the analyses we used the logging intensities calculated from the Landsat analysis, where possible the continuous data or otherwise the derived three logging classes. Differences in forest structures were tested amongst logging classes using ANOVA or, when distributions within groups were not normally distributed, the Kruskal-Wallis rank sum test. For parameters that showed differences among logging classes, a Tukey test (for parametric tests) or Dunn test (for nonparametric tests) was performed to identify which logging classes differ.

Explaining variances in AGC
Linear mixed effects models (lmm; using package lme4 in R; Bates et al., 2015) and multiple linear regression models (mlm) were fitted, in which the AGC ≥10 cm dbh was predicted by various fixed effects (for lmm) and independent parameters (for mlm) and their interactions, i.e. logging intensity, elevation, slope, basal area (B.A.) of logged stumps and soil properties. In the lmm's, the site, i.e. plots located in the same area, was treated as a random effect to reduce pseudo-replication. Soil parameters that explain most of the variances in the soil dataset were selected through a principal component analysis (PCA), in which both depth layers were combined. The selected soil parameters were averaged over the 0-10 cm and 10-30 cm depth layer and added to the mlm's and lmm's. Models were generated and the best-fitted models were selected using model averaging based on Bayesian Information Criterion (BIC). The importance of each parameter in the initial model was calculated, showing the proportional frequency that each parameter was selected in fitted models. Package relaimpo in R (Grömping, 2006) was used to calculate the relative importance of the selected parameters in explaining the variance in AGC in the final model. All analyses were carried out using R (version 3.5.3; R Core Team, 2019) and results expressed per hectare.  Table B.1 for logging classifications per plot).

Logging intensities in KNT
No disturbed pixels were detected within either of the 24 forest plots (Fig. 2). We calculated the disturbance intensity for each plot as the sum of disturbed pixels within 1000 m of the plot, with disturbances within 500 m of the plot weighted by a factor 2 (see Material and methods). At least three disturbed pixels were identified within 1000 m of all plots, with summed weighted disturbance varying by a factor 65 between the most heavily disturbed and most lightly disturbed plots. For most plots, disturbances within 1000 m occurred on multiple different years (mean = 8 different years, maximum 14 different years), suggesting ongoing disturbance rather than isolated harvesting years. Plots were classified into light, medium and heavy logging intensity based on the calculated disturbance intensity. Overall, there was good agreement between the logging classifications from the remote sensing analysis and participatory mapping (Spearman's ρ = 0.66), with a broadly similar distribution of plots identified as light, medium and heavy logging intensity (Table B.1). This agreement provides additional confidence that our remote sensing analysis is able to quantitatively assess historical logging activities in this region. Table 1 compares the forest plots according to the logging classifications derived from the Landsat analysis (Table B.2 shows details for each plot). The B.A. of logged stumps still identifiable in the plots, indicating the amount of logging in recent years, was significantly higher in medium logged compared to lightly logged forests (Table 1). Maximum stem diameter, B.A. of living stems and AGC of stems ≥10 cm significantly differed among logging classes, with bigger trees, higher B.A. and higher AGC stocks in lightly logged plots compared to medium and heavily logged plots. AGC of stems ≥10 cm dbh in lightly logged plots was double that in heavily logged plots (Fig. 3a), with consistent results using logging classifications derived from participatory mapping (Fig. 3b). For the further analyses, the logging intensity according to the Landsat analysis was used. The AGC of regenerating stems of 5-10 cm dbh, C in dead wood and SOC (both 0-10 cm and 10-30 cm depth) did not differ between the logging classes (Table 1). BGC and total C were significantly higher in lightly logged areas compared to forests logged with medium and heavy logging intensity.

Forest C stocks for different logging classes
Average total C stocks equalled 229, 173 and 163 Mg C ha −1 in lightly, medium and heavily logged forests, respectively (Table 1). The main difference in C stocks among logging classes was found in AGC stocks ≥10 cm dbh, i.e. heavily logged forests stored half the amount of AGC as lightly logged forests. This difference is largely attributed to the amount of AGC stored in stems ≥60 cm dbh: this C pool represents 17% of total C stocks in lightly logged forests, while these large stems were almost absent in medium and heavily logged forests (4 and 6% of total C stocks, respectively; Fig. 4a and b. See Figs. C.1 and C.2 for absolute and proportional C stocks per plot). Small amounts of C were stored in AGC of regenerating trees (5-10 cm dbh), dead wood and BGC. SOC in the top layer (0-30 cm depth) stored a large fraction of total C, representing respectively 43%, 54% and 55% in lightly, medium and heavily logged forests.
While AGC of stems ≥10 cm dbh and BGC differed with logging intensity, AGC of stems 5-10 cm dbh, C in dead wood and SOC did not vary significantly with logging intensity (Table 1). Hereafter, we explore the impacts of logging intensity, topography and soil on AGC of stems ≥10 cm dbh and assess the importance of logging intensity in predicting AGC stocks.

Explaining variances in AGC
We used a PCA analysis to select the soil parameters that explain most of the variances in the soil dataset. The six most dominant soil parameters were included in the models, i.e. % clay, exchangeable H + , exchangeable K + , exchangeable Ca 2+ , % sand and available P, however these and the next highest ranked soil parameters were almost all equally dominant (Fig. D).
To explore which parameters explain the variability in AGC of stems ≥10 cm dbh, linear mixed effects models (lmm) and multiple linear regression models (mlm) were fitted according to the following initial model, with the random effect "site" only in the lmm's: The disturbance parameters and topography values per plot are listed in Table B.1, the soil values in Table E. The multiplication in Eq.
(1) indicates that the interaction between parameters is also included.
The best lmm predicted AGC by logging intensity and elevation (fixed effects) and site (random effect) (R 2 = 0.81, BIC = 211). The best mlm predicted AGC by logging intensity and elevation (R 2 = 0.80, p < 0.001, BIC = 207). The mlm model had the lowest BIC and was therefore chosen to predict AGC. The model was fitted as follows: From the 24 mlm's that were fitted, logging intensity was selected in each model (1.0) and elevation in 19 models (0.83) (Fig. 5). The fitted AGC (using Eq. (2)) versus observed AGC shows a good relationship (Fig. 6). The best mlm (Eq. (2)) explained 80% of the variance in AGC of stems ≥10 cm dbh, of which 77% was explained by logging intensity and 23% by elevation. AGC ≥10 cm dbh declined significantly with logging intensity (Fig. 7).

Discussion
While most Vietnamese forests experienced logging activities in the past, detailed information on logging extent and intensity is often lacking. Here, we propose a novel approach to trace back historical logging activities and estimate logging intensities for sites where these data are absent, using a combination of Landsat analysis and participatory mapping. In this study, we assessed the effects of logging intensity on main forest C pools in a lowland forest in Vietnam. We found lower AGC of stems ≥10 cm dbh in medium and heavily logged compared to lightly logged forests, due to the reduction of stems ≥60 cm dbh, which were likely harvested. SOC in the 0-30 cm top layer stored up to 55% of total C stocks and is therefore an important C pool at our site. Fitting various models in which AGC of stems ≥10 cm dbh was related to logging intensity, topography and soil showed that logging intensity was the main factor explaining the variability in AGC between plots.

Logging intensities
Logged stumps were observed in 21 of 24 plots, indicating that most areas have been recently logged. No relationship was observed between the B.A. of these recently logged stumps and the AGC of stems ≥10 cm dbh, suggesting that recent logging has not caused the difference in AGC stocks and that logging in earlier periods may have been important. Disturbed locations identified with the Landsat analysis are mostly located along rivers and near the border of KNT, or are organised into linear patterns, suggesting a human origin. This, in combination with the high similarity in logging classifications based on the participatory mapping, are likely evidences that most disturbances identified on Landsat images are due to logging events.
Most plots in the west of KNT have been lightly logged, while the plots in the east have been logged with medium or heavy logging intensity. The lightly logged plots were located at slightly higher elevation (mean elevation of 370 m) compared to the plots logged with medium and heavy intensity (mean elevation of 233 and 254 m, respectively). Erythrophleum fordii, Sindora tonkinensis and Sindora siamensis, the most valuable timber species, usually grow in areas with lower elevation (pers. comm. Trai Trong Le and Anh Tuan Pham, 2018). More intense logging in the east can be explained by the fact that the state logging company was located just outside the eastern border of KNT and the area is more easily accessible due to the flatter terrain, lower elevations and the presence of big rivers.
Testing various models showed that the logging intensity and elevation are the most important predictors for plot-level AGC of stems ≥10 cm dbh, with highest AGC values in lightly logged forests and forests at higher elevation. The relationship between logging and AGC is consistent with the results of meta-analyses, which showed that logging intensity was the best predictor for post-logging changes in AGC Table 1 Parameters describing the forest structure in lightly (n = 8), medium (n = 8) and heavily (n = 8) logged forests. Logging intensity was classified according to the results of the satellite analysis. Total C is a summation of the C stocks in AGC of stems ≥5 cm dbh, C in dead wood, BGC and SOC 0-30 cm depth. Values represent the mean ± standard deviation. Significant p-values are annotated with *, different letters within a given row indicate significant differences between logging classes. KW = Kruskal-Wallis rank sum test.    Fig. 6. Fitted AGC (using Eq. (2)) versus measured AGC of stems ≥10 cm dbh. and changes in AGC were negatively correlated with logging intensity (Martin et al., 2015;Zhou et al., 2013). A study in Indonesia, similar to our study, tested the effects of logging intensity, topography and soil variables on several main forest C stocks and also found that logging intensity was the main factor in explaining the variability in AGC of trees > 20 cm dbh (Rozak et al., 2018). Previous work in unlogged forests showed that AGC decreased or had a weak relationship with elevation (Aiba and Kitayama, 1999;Álvarez-Dávila et al., 2017;Kitayama and Aiba, 2002;Spracklen and Righelato, 2014). In contrast, in our study we observed a positive relationship between AGC and elevation. Although the interaction between logging intensity and elevation was not selected in our final model, i.e. any possible interaction between these two parameters may only explain a small fraction of the variation in AGC stocks, it is possible that greater AGC at higher elevation is a result of less logging in these plots rather than a direct result of elevation.

C stocks in logged forests
Differences in total C stocks between logging classes were largely found in the AGC stock of stems ≥10 cm dbh. Big trees (≥60 dbh) were largely reduced in medium and heavily logged forests, most likely removed by logging. Big trees store large amounts of AGC and their density and survival after logging largely determines C dynamics (Rozak et al., 2018;Sist et al., 2014;Slik et al., 2013).
The mean AGC of trees ≥10 cm dbh in our study varied from 108 Mg C ha −1 in lightly logged, 56 Mg C ha −1 in medium logged to 52 Mg C ha −1 in heavily logged forests. Our AGC stocks are similar to other logged forests in Vietnam and Malaysia (Hai et al., 2015;Luong et al., 2015;Nam et al., 2018;Saner et al., 2012) and old growth forests in Vietnam (Con et al., 2013), but considerably lower than AGC stocks in logged Dipterocarp forests in Indonesia (Rozak et al., 2018). Trees from the Dipterocarpaceae family consist of large canopy and emergent trees (Ashton, 1982;Whitmore, 1984) and are important C stocks, but Dipterocarps did not occur in our plots. In our study, the AGC of stems ≥10 cm dbh in medium and heavily logged forests was, respectively, 48% and 52% lower than the AGC in lightly logged forests. Several studies found a considerable decline in AGC with logging intensity, but the exact reduction varies between studies, which is likely caused by variations in logging intensity, logging method and time since logging. For instance, forests logged with medium and high logging intensity in Indonesia stored, respectively, 9% and 33% less AGC of stems > 20 cm dbh than lightly logged sites (Rozak et al., 2018). In Brazil, the AGC of trees ≥10 cm dbh of an intact forest decreased by 23% after logging with moderate intensity and by 43% after logging with heavy intensity (Gerwing, 2002). AGC stocks of stems > 10 cm dbh in a forest 22 years post-logging were 28% lower compared to an unlogged forest in Malaysian Borneo (Saner et al., 2012). Overall, a meta-analysis found that selective logging reduced AGC stocks on average by 43% (Zhou et al., 2013). Small trees of 5-10 cm dbh stored only small amounts of AGC in our sites, i.e. 3-4 Mg C ha −1 , and did not differ with logging intensity. Similar values were found in logged forests in Vietnam and Malaysia (Nam et al., 2018;Saner et al., 2012; Table 2). Like our study, no correlation was observed between the logging intensity and AGC of small stems (5-20 cm dbh) in Indonesia (Rozak et al., 2018). However, in forests in Brazil, the AGC of small trees (> 2 m tall and < 10 cm dbh) decreased in heavily logged forests (Gerwing, 2002).
Logging can cause a temporal increase in necromass, which reduces over time when decomposition proceeds. Mean C stocks in dead wood at our sites ranged from 3 to 6 Mg C ha −1 and did not differ with logging intensity. A similar amount of C in necromass was found elsewhere in Vietnam (Hai et al., 2015;Nam et al., 2018; Table 2). Considerably higher values of C in necromass were observed in logged forests in Malaysia  and Indonesia (Rozak et al., 2018) (Table 2). Further, in forests in Indonesia and Brazil, the amount of necromass increased with logging intensity (Gerwing, 2002;Rozak et al., 2018). Previous work showed that necromass is sensitive to wood density values of living trees, i.e. necromass is positively related to plotlevel average living wood density (Chao et al., 2009). In our study, mean wood density of living stems ≥10 cm dbh was low in each plot, ranging from 0.505 to 0.604 g cm −3 , and the majority of the logging occurred more than nine years before field measurements. Therefore, it is likely that much of the necromass caused by logging has already decayed, which contributed to low necromass values and no significant differences in necromass with logging intensity in our forests.
BGC stocks of our forests ranged from 10 to 15 Mg C ha −1 and are similar to values in other Southeast Asian forests (Nam et al., 2018;Saner et al., 2012; Table 2). Higher BGC stocks were found in Dipterocarp forests in Indonesia (Rozak et al. 2018; Table 2), in which BGC was estimated using an allometric model based on dbh. Mature Dipterocarp trees often have large diameters, which likely explains the high BGC values found in these forests. In our study, BGC significantly decreased with logging intensity. We calculated BGC using an allometric equation with dbh and wsg as input parameters and the reduction of large trees in medium and heavily logged forests likely explains the reduced BGC values in those forests. Consistent to our results, the BGC Fig. 7. The relationship between AGC ≥10 cm dbh and the logarithmic logging intensity. S.M. Stas, et al. Forest Ecology and Management 460 (2020) 117863 of trees > 20 cm dbh decreased with logging intensity in Indonesian forests (Rozak et al., 2018). Our plots are located in one geographical area, which resulted in similar soil properties across plots and none of the selected soil parameters contributed to explain sufficient variation in AGC. SOC ranged from 93 to 99 Mg C ha −1 in the top 0-30 cm layer, did not differ significantly between the logging classes and comprised more than half of total C stocks in our medium and heavily logged forests. Several studies estimated SOC values in Southeast Asia, but the depth layer varies per study (Table 2). SOC values in our plots are similar to other studies in Vietnam (Hai et al., 2015;Nam et al., 2018;Sang et al., 2013) and larger than most soils in tropical Southeast Asia (Rozak et al., 2018;Saner et al., 2012) (Table 2). This suggests that Vietnamese forests in general have higher SOC storage, possibly because the forests are in the sub-tropics rather than tropics, leading to slower decomposition rates. In Indonesian forests, no relationship was observed either between logging intensity and SOC (Rozak et al., 2018).

Implications for forest management
Our results suggest that historical logging has reduced the C storage of forests in a protected natural forest in Vietnam. Large areas of natural forests in Vietnam have been degraded by logging and are likely to store less C than unlogged forests. If no further logging occurs in these forests, they will recover and sequester C, providing a substantial potential for Vietnamese forests to mitigate climate change. However, ongoing illegal logging threatens the recovery of AGC in our forests and is likely to be an issue more widely across Vietnam. Efforts to reduce illegal logging whilst supporting livelihoods of local communities are now required. Vietnam is starting to establish a REDD+ programme and is one of the countries selected to participate in the Forest Carbon Partnership Facility. By certifying the project activities under these programmes or through projectbased schemes such as the Climate, Community and Biodiversity Alliance (CCBA; CCBA, 2013), it is possible that further forest degradation will be prevented and natural recovery favoured.

Conclusions
We presented an approach to estimate historical logging intensities by combining Landsat analysis and participatory mapping, which could be more widely applicable for sites that lack historical logging data. Heavily logged forests stored only half the amount of AGC of stems ≥10 cm dbh as lightly logged forests, mainly due to a reduction in large (≥60 cm dbh) trees. Conserving big trees thus remains important to maintain high forest C stocks. SOC in the 0-30 cm top layer stored around half of total C stocks in these forests, showing the importance of accounting for this C stock. While topography and soil properties can influence AGC stocks, our study found that logging intensity was the most important factor explaining the variability in AGC of stems ≥10 cm dbh in lowland forests in Vietnam. Timber demand sharply increased after the Vietnam War, which led to intense state logging in forests throughout the country. Our study contributes to the understanding of the effects of long-term logging activities on forest C stocks in central Vietnam. These data are crucial to develop REDD+ approaches to prevent further degradation and allow eventual recovery of Vietnamese logged forests.

Author contributions
SMS, DVS, OLP, MvK and ER designed the research; SMS, TCL, HDT, TTL and DVS organised and/or conducted the fieldwork; DTN and AVL led the participatory mapping; TTAT and TTHH analysed the soil samples; AvO and SMS wrote the soil laboratory methods; BDS and DVS conducted the satellite analysis; SMS and ER analysed the data; SMS wrote the manuscript, with contributions from all authors.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. * Poor, medium and rich forests refer to the classification used by the Vietnamese government and is based on timber reserves of standing trees. The forests in the study of Hai et al. (2015) varied from 10 to 30 years post-logging and/or shifting cultivation to low disturbed old growth forests. Pfeifer et al. (2015) studied lightly or illegally logged forests and twice-logged forests (~10 years after the last logging cycle).

Acknowledgments
This study was financially supported by World Land Trust (UK) and the Newton Fund Institutional Links programme of the British Council (UK; grant number 216372155). DVS acknowledges the Natural Environment Research Council (NERC, UK; NE/M003574/1) and a Philip Leverhulme Prize of The Leverhulme Trust (UK). We are very grateful for all the support that we received from Roger Wilson, Anh Tuan Pham and other staff at World Land Trust and Viet Nature Conservation Centre. We thank Stephan Mantel, Imam Basuki and Yves Laumonier for their help in developing the methods for soil collection and Wolfgang Buermann and Tim Kasoar for creating maps. We are grateful to Tue Van Ha and Tuan Manh Le for tree species identifications, Heinjo During for statistical advice and the reviewers and editor for their constructive comments.

Appendix A
The decay classes used in our study (Baker and Chao, 2011) and the reduction in wsg for each decay class relative to the wsg in living trees, as estimated by Chao et al. (2017) in tropical lowland forests in Taiwan (Table A.1). For each plot, the wsg per decay class was estimated (Table A.2) by multiplying the mean plot-level wsg of living trees ≥10 cm dbh with the proportional reduction in wsg according to Chao et al. (2017).

Decay class Description
Relative reduction wsg (%) 1 Solid wood, recently fallen, with intact bark and fine branches (< 10 cm) still attached 32 1.5 Solid wood, but with no fine branches, and bark starting to fall off 42 2 Non-solid wood, in poorer condition, but where it is still difficult to push a nail into the wood by hand 51 2.5 Soft, rotten wood, where a nail can be pushed into the wood easily 54 3 Soft, rotten wood, which collapses easily when stepped on 61 Plot Mean plot-level wsg living trees ≥ 10 cm dbh (g cm −3 ) wsg class 1 (g cm −3 ) wsg class 1.5 (g cm −3 ) wsg class 2 (g cm −3 ) wsg class 2.5 (g cm −3 ) wsg class 3 (g cm −3 )

Appendix B
Parameters describing disturbances, site information and topography (Table B.1) and forest structure (Table B.2) per plot, ordered by the logging intensity derived from the Landsat analysis. Total C includes AGC ≥5 cm dbh, C in dead wood, BGC and SOC 0-30 cm depth (    Stas, et al. Forest Ecology and Management 460 (2020) 117863