Assessing impacts of selective logging on water, energy, and carbon budgets and ecosystem dynamics in Amazon forests using the Functionally Assembled Terrestrial Ecosystem Simulator

Tropical forest degradation from logging, fire, and fragmentation not only alters carbon stocks and carbon fluxes, but also impacts physical land surface properties such as albedo and roughness length. Such impacts are poorly quantified to date due to difficulties in accessing and maintaining observational infrastructures, as well as the lack of proper modeling tools for capturing the interactions among biophysical properties, ecosystem demography, canopy structure, and biogeochemical cycling in tropical forests. As a first step to address these limitations, we implemented a selective logging module into the Functionally Assembled Terrestrial Ecosystem Simulator (FATES) by mimicking the ecological, biophysical, and biogeochemical processes following a logging event. The model can specify the timing and aerial extent of logging events, splitting the logged forest patch into disturbed and intact patches; determine the survivorship of cohorts in the disturbed patch; and modifying the biomass and necromass (total mass of coarse woody debris and litter) pools following logging. We parameterized the logging module to reproduce a selective logging experiment at the Tapajós National Forest in Brazil and benchmarked model outputs against available field measurements. Our results suggest that the model permits the coexistence of early and late successional functional types and realistically characterizes the seasonality of water and carbon fluxes and stocks, the forest structure and composition, and the ecosystem succession following disturbance. However, the current version of FATES overestimates water stress in the dry season and therefore fails to capture seasonal variation in latent and sensible heat fluxes. Moreover, we observed a bias towards low stem density and leaf area when compared to observations, suggesting that improvements are needed in both carbon allocation and establishment of trees. The effects of logging were assessed by different logging scenarios to represent reduced impact and conventional logging practices, both with high and low logging intensities. The model simulations suggest that in comparison to old-growth forests the logged forests rapidly recover water and energy fluxes in 1 to 3 years. In contrast, the recovery times for carbon stocks, forest structure, and composition are more than 30 years depending on logging practices and intensity. This study lays the foundation to simulate land use change and forest degradation in FATES, which will be an effective tool to directly represent forest management practices and regeneration in the context of Earth system models.

1 Introduction

Land cover and land use in tropical forest regions are highly dynamic, and nearly all tropical forests are subject to significant human influence (Martínez-Ramos et al., 2016; Dirzo et al., 2014). While old-growth tropical forests have been reported to be carbon sinks that remove carbon dioxide from the atmosphere through photosynthesis, these forests could easily become carbon sources once disturbed (Luyssaert et al., 2008). Using data from forest inventory and long-term ecosystem carbon studies from 1990 to 2007, Pan et al. (2011) suggested a net tropical forest can be a net source of carbon source of 1.3±0.7 Pg C yr⁻¹ from land use change, consisting of a gross tropical deforestation loss of 2.9±0.5 Pg C yr⁻¹ that is partially offset by a carbon uptake by tropical secondary forest regrowth of 1.6±0.5 Pg C yr⁻¹. These estimates, however, do not account for a tropical forest that has been degraded through the combined effects of selective logging (cutting and removal of merchantable timber), fuelwood harvest, understory fires, and fragmentation (Nepstad et al., 1999; Bradshaw et al., 2009). To date, the effects of forest degradation remain poorly quantified. Recent studies suggested that degradation may contribute to carbon loss 40 % as large as clear-cut deforestation (Berenguer et al., 2014), and the emission from selective logging alone could be equivalent to ∼10 % to 50 % of that from deforestation in the tropical countries (Pearson et al., 2014; Huang and Asner, 2010; Asner et al., 2009). Selective logging of tropical forests is an important contributor to many local and national economies and corresponds to approximately one-eighth of global timber (Blaser et al., 2011). The integrated impact of timber production and other forest uses has been posited as the cause of up to ∼30 % of the difference between potential and actual biomass stocks globally, comparable in magnitude to the effects of deforestation (Erb et al., 2017). Selective logging includes cutting large trees and additional degradation through widespread damage to remaining trees, subcanopy vegetation, and soils (Asner et al., 2004; Asner et al., 2005). Selective logging accelerates gap-phase regeneration within the degraded forests (Huang et al., 2008).

Over half of all tropical forests have been cleared or logged, and almost half of standing old-growth tropical forests are designated by national forest services for timber production (Sist et al., 2015). Disturbances that result from logging are known to cause forest degradation at the same magnitude as deforestation each year in terms of both geographic extent and intensity, with widespread collateral damage to remaining trees, vegetation, and soils, leading to disturbance to water, energy, and carbon cycling, as well as ecosystem integrity (Keller et al., 2004b; Asner et al., 2004; Huang and Asner, 2010).

In most Earth system models (ESMs) that couple terrestrial and atmospheric processes to investigate global change (e.g., the Community Earth System Model or the Energy Exascale Earth System Model), selective logging is typically represented as simple fractions of affected area or an amount of carbon to be removed on a coarse grid (e.g., 0.5^∘). One exception is the representation of wood harvest in the LM3V land model that explicitly accounts for postdisturbance land age distribution, as part of the Geophysical Fluid Dynamics Laboratory (GFDL) Earth system model (Shevliakova et al., 2009). In the ESMs, grid cell fractional areas are typically based on timber production rates estimated from sawmill, sales, and export statistics (Hurtt et al., 2011; Lawrence et al., 2012). This approach, while practical, does not effectively differentiate selective logging that retains forest cover from deforestation.

The realistic representation of wood harvest was absent in most ESMs because the models generally did not represent the demographic structure of forests (tree size and stem number distributions; Bonan, 2008). But progress over the past two decades in ecological theory and observations (Bustamante et al., 2016; Strigul et al., 2008; Hurtt et al., 1998; Moorcroft et al., 2001) has made it feasible to include vegetation demography more directly into Earth system models through individual to cohort-based vegetation in land models (Sato et al., 2007; Watanabe et al., 2011; Smith et al., 2001; Smith et al., 2014; Weng et al., 2015; Baidya Roy et al., 2003; Hurtt et al., 1998; Fisher et al., 2015). These vegetation demography modules are relatively new in land models, so efforts are still underway to improve their parameterizations of resource competition for light, water, nutrients, recruitment, mortality, and disturbance including both natural and anthropogenic components (Fisher et al., 2017).

In this study, we aim to (1) describe the development of a selective logging module implemented into The Functionally Assembled Terrestrial Ecosystem Simulator (FATES), for simulating anthropogenic disturbances of various intensities to forest ecosystems and their short-term and long-term effects on water, energy, carbon cycling, and ecosystem dynamics; (2) assess the capability of FATES in simulating site-level water, energy, and carbon budgets, as well as forest structure and composition; (3) benchmark the simulated variables against available observations at the Tapajós National Forest in the Amazon, thus identifying potential directions for model improvement; and (4) assess the simulated recovery trajectory of tropical forest following disturbance under various logging scenarios. In Sect. 2, we provide a brief summary of FATES, introduce the new selective logging module, and describe numerical experiments performed at two sites with data from field surveys and flux towers. In Sect. 3, FATES-simulated water, energy, and carbon fluxes and stocks in intact and disturbed forests are compared to available observations, and the effects of logging practice and intensity on simulated forest recovery trajectory in terms of carbon budget, size structure, and composition in plant functional types are assessed. Conclusions and future work are discussed in Sect. 4.

2 Model description and study site

2.1 The Functionally Assembled Terrestrial Ecosystem Simulator

The Functionally Assembled Terrestrial Ecosystem Simulator (FATES) has been developed as a numerical terrestrial ecosystem model based on the ecosystem demography representation in the community land model (CLM), formerly known as CLM(ED) (Fisher et al., 2015). FATES is an implementation of the cohort-based ecosystem demography (ED) concept (Hurtt et al., 1998; Moorcroft et al., 2001) that can be called as a library from an ESM land surface scheme, currently including the CLM (Oleson et al., 2013) or Energy Exascale Earth System Model (E3SM) land model (ELM) (https://climatemodeling.science.energy.gov/projects/energy-exascale-earth-system-model, last access: 28 June 2020). In FATES, the landscape is discretized into spatially implicit patches, each of which represents land areas with a similar age since last disturbance. The discretization of ecosystems along a disturbance–recovery axis allows the deterministic simulation of successional dynamics within a typical forest ecosystem. Within each patch, individuals are grouped into cohorts by plant functional types (PFTs) and size classes (SCs), so that cohorts can compete for light based on their heights and canopy positions. Following disturbance, a patch fission process splits the original patch into undisturbed and disturbed new patches. A patch fusion mechanism is implemented to merge patches with similar structures, which helps prevent the number of patches from growing too big. In addition to the ED concept, FATES also adopted a modified version of the perfect plasticity approximation (PPA) (Strigul et al., 2008) concept by splitting growing cohorts between canopy and understory layers as a continuous function of height designed for increasing the probability of coexistence (Fisher et al., 2010). An earlier version of FATES, CLM(ED), has been applied regionally to explore the sensitivity of biome boundaries to plant trait representation (Fisher et al., 2015).

In this study, we specified two plant functional types (PFTs) in FATES corresponding to early successional and late successional plants, representative of the primary axis of variability in tropical forests (Reich, 2014). The early successional PFT is light demanding and grows rapidly under high-light conditions common prior to canopy closure. This PFT has low-density woody tissues, shorter leaf and root lifetimes, and a higher background mortality compared to the late successional PFT that has dense woody tissues, longer leaf and root lifetimes, and lower background mortality (Brokaw, 1985; Whitmore, 1998) and thus can survive under deep shade and grow slowly under closed canopy.

The key parameters that differentiate the two PFTs in FATES are listed in Table 1, including specific leaf area at the canopy top (SLA₀), the maximum rate of carboxylation at 25 ^∘C (V_cmax25), specific wood density, background mortality, leaf and fine root longevity, and leaf C:N ratio. The parameter ranges were selected based on literature for tropical forests. Specifically, it has been reported that SLA values ranges from 0.007 to 0.039 m² g C⁻¹ (Wright et al., 2004) and V_cmax25 ranges between 10.1 and 105.7 µmol m⁻² s⁻¹ (Domingues et al., 2005). The specific wood densities were set to be 0.5 and 0.9 g cm³ , and the background mortality rates were set to 0.035 and 0.014 yr⁻¹ for early and late succession PFTs respectively, consistent with those used in the Ecosystem Demography Model version 2 for Amazon forests (Longo et al., 2019) . For simplicity, leaf longevity and root longevity were set to be the same for each PFT (i.e., 0.9 and 2.6 years for early and late successional PFTs) following the range in Trumbore and Barbosa De Camargo (2009).

Table 1 FATES parameters that define early and late successional PFTs.

Given that both SLA₀ and V_cmax25 span wide ranges and have been identified as the most sensitive parameters in FATES in a previous study (Massoud et al., 2019), we performed one-at-a-time sensitivity tests by perturbing them within the reported ranges. Based on these tests, it is evident that these parameters not only affect water, energy, and carbon budget simulations, but also the coexistence of the two PFTs. In the version of FATES used in this study (Interested readers are referred to the “Code and data availability” section for details), coexistence of PFTs is not assured for all parameter combinations, even if they are both within reasonable ranges, on account of competitive exclusion feedback processes that prevent coexistence in the presence of large discrepancies in plant growth and reproduction rates (Fisher et al., 2010; Bohn et al., 2011). In order to demonstrate FATES' capability in simulating water, energy, carbon budgets, and forest structure and composition in a holistic way, we chose to report results based on a set of parameter values that produces reasonable, stable fractions of two PFTs, as reported in Table 1. Nevertheless, we have included a summary of all sensitivity tests performed in the Supplement for completeness. The sensitivity tests demonstrated that, by tuning SLA₀ and V_cmax25 for the different PFTs, FATES is not only capable of capturing coexistence of PFTs, but also capable of reproducing observed water, energy, and carbon cycle fluxes in the tropics.

Figure 1 (a) Landscape components of selective logging; (b) location of the Tapajós National Forest in the Amazon; and (c) a typical logging block showing tree-fall location, skid trail, road, and log deck coverages. Panels (b) and (c) are from Asner et al. (2008).

2.2 The selective logging module

The new selective logging module in FATES mimics the ecological, biophysical, and biogeochemical processes following a logging event. The module (1) specifies the timing and areal extent of a logging event; (2) calculates the fractions of trees that are damaged by direct felling, collateral damage, and infrastructure damage and adds these size-specific plant mortality types to FATES; (3) splits the logged patch into disturbed and intact new patches; (4) applies the calculated survivorship to cohorts in the disturbed patch; and (5) transports harvested logs off-site by reducing site carbon pools and adds remaining necromass to coarse woody debris and litter pools.

The logging module structure and parameterization are based on detailed field and remote sensing studies (Putz et al., 2008; Asner et al., 2004; Pereira et al., 2002; Asner et al., 2005; Feldpausch et al., 2005). Logging infrastructure including roads, skids, trails, and log decks are conceptually represented (Fig. 1). The construction of log decks used to store logs prior to road transport leads to large canopy openings, but their contribution to landscape-level gap dynamics is small. In contrast, the canopy gaps caused by tree felling are small, but their coverage is spatially extensive at the landscape scale. Variations in logging practices significantly affect the level of disturbance to tropical forests following logging (Pereira et al., 2002; Macpherson et al., 2012; Dykstra, 2002; Putz et al., 2008). Logging operations in the tropics are often carried out with little planning and typically use heavy machinery to access the forests accompanied by construction of excessive roads and skid trails, leading to unnecessary tree fall and compaction of the soil. We refer to these typical operations as conventional logging (CL). In contrast, reduced impact logging (RIL) is a practice with extensive preharvest planning, where trees are inventoried and mapped out for the most efficient and cost-effective harvest and seed trees are deliberately left on-site to facilitate faster recovery. Through planning, the construction of skid trails and roads, soil compaction, and disturbance can be minimized. Vines connecting trees are cut, and tree-fall directions are controlled to reduce damages to surrounding trees. Reduced impact logging results in consistently less disturbance to forests than conventional logging (Pereira et al., 2002; Putz et al., 2008).

The FATES logging module was designed to represent a range of logging practices in field operations at a landscape level. Both CL and RIL can be represented in FATES by specifying mortality-rate-associated direct felling, collateral damages, and mechanical damages as follows: once logging events are activated, we define three types of mortality associated with logging practices – direct-felling mortality (lmort_direct), collateral mortality (lmort_collateral), and mechanical mortality (lmort_mechanical). The direct-felling mortality represents the fraction of trees selected for harvesting that are greater than or equal to a diameter threshold (this threshold is defined by the diameter at breast height (DBH) = 1.3 m denoted as DBH_min); collateral mortality denotes the fraction of adjacent trees that are killed by felling of the harvested trees; and the mechanical mortality represents the fraction of trees killed by construction of log decks, skid trails, and roads for accessing the harvested trees, as well as storing and transporting logs off-site (Fig. 1a). In a logging operation, the loggers typically avoid large trees when they build log decks, skids, and trails by knocking down relatively small trees as it is not economical to knock down large trees. Therefore, we implemented another DBH threshold, DBH${}_{{\max }_{\mathrm{infra}}}$, so that only a fraction of trees ≤ DBH${}_{{\max }_{\mathrm{infra}}}$ (called mechanical damage fraction) are removed for building infrastructure (Feldpausch et al., 2005).

To capture the disturbance mechanisms and degree of damage associated with logging practices at the landscape level, we apply the mortality types following a workflow designed to correspond to field operations. In FATES, as illustrated in Fig. 2, individual trees of all plant functional types (PFTs) in one patch are grouped into cohorts of similar-sized trees, whose size and population sizes evolve in time through processes of recruitment, growth, and mortality. For the purpose of reporting and visualizing the model state, these cohorts are binned into a set of 13 fixed size classes in terms of the diameter at the breast height (DBH) (i.e., 0–5, 5–10, 10–15, 15–20, 20–30, 30–40, 40–50, 50–60, 60–70, 70–80, 80–90, 90–100, and ≥100 cm). Cohorts are further organized into canopy and understory layers, which are subject to different light conditions (Fig. 2a). When logging activities occur, the canopy trees and a portion of big understory trees lose their crown coverage through direct felling for harvesting logs, or as a result of collateral and mechanical damages (Fig. 2b). The fractions of the canopy trees affected by the three mortality mechanisms are then summed up to specify the areal percentages of an old (undisturbed) and a new (disturbed) patch caused by logging in the patch fission process as discussed in Sect. 2.1 (Fig. 2c). After patch fission, the canopy layer over the disturbed patch is removed, while that over the undisturbed patch stays untouched (Fig. 2d). In the undisturbed patch, the survivorship of understory trees is calculated using an understory death fraction consistent with the default value corresponding to that used for natural disturbance (i.e., 0.5598). To differentiate logging from natural disturbance, a slightly elevated, logging-specific understory death fraction is applied in the disturbed patch instead at the time of the logging event. Based on data from field surveys over logged forest plots in the southern Amazon (Feldpausch et al., 2005), the understory death fraction corresponding to logging is now set to be 0.65 as the default but can be modified via the FATES parameter file (Fig. 2e). Therefore, the logging operations will change the forest from the undisturbed state shown in Fig. 2a to a disturbed state in Fig. 2f in the logging module. It is worth mentioning that the newly generated patches are tracked according to the age since disturbance and will be merged with other patches of similar canopy structure following the patch fusion processes in FATES in later time steps of a simulation, pending the inclusion of separate land use fractions for managed and unmanaged forest.

Figure 2 The mortality types (direct felling, mechanical, and collateral) and patch-generating process in the FATES logging module. The white fraction in panels (c), (d), and (f) indicates mortality associated with other disturbances in FATES. (a) Canopy and understory layers in each cohort in FATES; (b) mortality applied at the time of a logging event; (c) the patch fission process following a given logging event; (d) canopy removal in the disturbed patch following the logging event; (e) calculation of the understory survivorship based on the understory death fraction in each patch; (d) the final states of the intact and disturbed patches.

Figure 3 The flow of necromass following logging.

Logging operations affect forest structure and composition as well as also carbon cycling (Palace et al., 2008) by modifying the live biomass pools and flow of necromass (Fig. 3). Following a logging event, the logged trunk products from the harvested trees are transported off-site (as an added carbon pool for resource management in the model), while their branches enter the coarse woody debris (CWD) pool, and their leaves and fine roots enter the litter pool. Similarly, trunks and branches of the dead trees caused by collateral and mechanical damages also become CWD, while their leaves and fine roots become litter. Specifically, the densities of dead trees as a result of direct felling, collateral, and mechanical damages in a cohort are calculated as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{l}{D}_{\mathrm{direct}}={\mathrm{lmort}}_{\mathrm{direct}}\times \frac{n}{A}\\ D_{\mathrm{collateral}}={\mathrm{lmort}}_{\mathrm{collateral}}\times \frac{n}{A}\\ D_{\mathrm{mechanical}}={\mathrm{lmort}}_{\mathrm{mechanical}}\times \frac{n}{A},\end{array}\end{array}$$

where A stands for the area of the patch being logged, and n is the number of individuals in the cohort where the mortality types apply (i.e., as specified by the size thresholds, DBH_min and DBH${}_{{\max }_{\mathrm{infra}}}$). For each cohort, we denote $D_{\mathrm{indirect}}=D_{\mathrm{collateral}}+D_{\mathrm{mechanical}}$ and $D_{\mathrm{total}}=D_{\mathrm{direct}}+D_{\mathrm{indirect}}$.

Leaf litter (Litter_leaf, kg C) and root litter (Litter_root, kg C) at the cohort level are then calculated as

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{\mathrm{Litter}}_{\mathrm{leaf}}=D_{\mathrm{total}}\times B_{\mathrm{leaf}}\times A,\text{(3)}& {\displaystyle }& {\displaystyle }{\mathrm{Litter}}_{\mathrm{root}}=D_{\mathrm{total}}\times \left(B_{\mathrm{root}}+B_{\mathrm{store}}\right)\times A,\end{array}$$

where B_leaf and B_root are live biomass in leaves and fine roots, respectively, and B_store is stored biomass in the labile carbon reserve in all individual trees in the cohort of interest.

Following the existing CWD structure in FATES (Fisher et al., 2015), CWD in the logging module is first separated into two categories: aboveground CWD and belowground CWD. Within each category, four size classes are tracked based on their source, following Thonicke et al. (2010): trunks, large branches, small branches, and twigs. Aboveground CWD from trunks (CWD_{trunk_agb}, kg C) and branches and twigs (CWD_{branch_agb}, kg C) are calculated as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{\mathrm{CWD}}_{\mathrm{trunk}\mathrm{\_}\mathrm{agb}}=D_{\mathrm{indirect}}\times B_{\mathrm{stem}\mathrm{\_}\mathrm{agb}}\times f_{\mathrm{trunk}}\times A,\text{(5)}& {\displaystyle }& {\displaystyle }{\mathrm{CWD}}_{\mathrm{branch}\mathrm{\_}\mathrm{agb}}=D_{\mathrm{total}}\times B_{\mathrm{stem}\mathrm{\_}\mathrm{agb}}\times f_{\mathrm{branch}}\times A,\end{array}$$

where B_{stem_agb} is the amount of aboveground stem biomass in the cohort, and f_trunk and f_branch represent the fraction of trunks and branches (e.g., large branches, small branches, and twigs). Similarly, the belowground CWD from trunks (CWD_{trunk_bg}, kg C) and branches and twigs (CWD_{branch_bg}, kg C) are calculated as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathrm{CWD}}_{\mathrm{trunk}\mathrm{\_}\mathrm{bg}}=D_{\mathrm{total}}\times B_{\mathrm{root}\mathrm{\_}\mathrm{bg}}\times f_{\mathrm{trunk}}\times A,\text{(7)}& {\displaystyle }& {\displaystyle }{\mathrm{CWD}}_{\mathrm{branch}\mathrm{\_}\mathrm{bg}}=D_{\mathrm{total}}\times B_{\mathrm{root}\mathrm{\_}\mathrm{bg}}\times f_{\mathrm{branch}}\times A,\end{array}$$

where B_croot (kg C) is the amount of coarse root biomass in the cohort. Site-level total litter and CWD inputs can then be obtained by integrating the corresponding pools over all the cohorts in the site. To ensure mass conservation, the total loss of live biomass due to logging, ΔB (i.e., carbon in leaf, fine roots, storage, and structural pools), needs to be balanced with increases in litter and CWD pools and the carbon stored in harvested logs shipped off-site as follows:

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }B=\mathrm{\Delta }\mathrm{Litter}+\mathrm{\Delta }\mathrm{CWD}+\mathrm{trunk}\mathrm{\_}\mathrm{product},\end{array}$$

where Δlitter and ΔCWD are the increments in litter and CWD pools, and trunk_product represents harvested logs shipped off-site.

Following the logging event, the forest structure and composition in terms of cohort distributions, as well as the live biomass and necromass pools, are updated. Following this logging event update to forest structure, the native processes simulating physiology, growth, and competition for resources in and between cohorts resume. Since the canopy layer is removed in the disturbed patch, the existing understory trees are promoted to the canopy layer, but, in general, the canopy is incompletely filled in by these newly promoted trees, and thus the canopy does not fully close. Therefore, more light can penetrate and reach the understory layer in the disturbed patch, leading to increases in light-demanding species in the early stage of regeneration, followed by a succession process in which shade-tolerant species dominate gradually.

Table 2 Distributions of stem density (N ha⁻¹), basal area (m² ha⁻¹), and aboveground biomass (Kg C m⁻²) before and after logging at km83, separated by diameter of breast height (normal text) and aggregated across all sizes (bold text).

Data in this table obtained based on inventory during the LBA period (Menton et al., 2011; de Sousa et al., 2011).

2.3 Study site and data

In this study, we used data from two evergreen tropical forest sites located in the Tapajós National Forest, Brazil (Fig. 1b). These sites were established during the Large-Scale Biosphere-Atmosphere Experiment in Amazonia (LBA) and are selected because of data availability including those from forest plot surveys and two flux towers established during the LBA period (Keller et al., 2004a). These sites were named after distances along the BR-163 highway from Santarém: km67 (2^∘51^′ S, 54^∘58^′ W) and km83 (3^∘3^′ S, 54^∘56^′ W). They are situated on a flat plateau and were established as a control-treatment pair for a selective logging experiment. Tree felling operations were initiated at km83 in September 2001 for a period of about 2 months. Both sites are similar with a mean annual precipitation of ∼2000 mm and mean annual temperature of 25 ^∘C, on nutrient-poor clay Oxisols with low organic content (Silver et al., 2000).

Prior to logging, both sites were old-growth forests with limited previous human disturbances caused by hunting, gathering Brazil nuts, and similar activities. A comprehensive set of meteorological variables, as well as land–atmosphere exchanges of water, energy, and carbon fluxes, have been measured by an eddy covariance tower at an hourly time step over the period of 2002 to 2011, including precipitation, air temperature, surface pressure, relative humidity, incoming shortwave and longwave radiation, latent and sensible heat fluxes, and net ecosystem exchange (NEE) (Hayek et al., 2018). Another flux tower was established at km83, the logged site, with hourly meteorological and eddy covariance measurements in the period of 2000–2003 (Miller et al., 2004; Goulden et al., 2004; Saleska et al., 2003). The towers are listed as BR-Sa1 and BR-Sa3 in the AmeriFlux network (https://ameriflux.lbl.gov, last access: 28 June 2020).

These tower- and biometric-based observations were summarized to quantify logging-induced perturbations on old-growth Amazonian forests in Miller et al. (2011) and are used in this study to benchmark the model-simulated carbon budget. Over the period of 1999 to 2001, all trees ≥35 cm in DBH in 20 ha of forest in four 1 km long transects within the km67 footprint were inventoried, as well as trees ≥10 cm in DBH on subplots with an area of ∼4 ha. At km83, inventory surveys on trees ≥55 cm in DBH were conducted in 1984 and 2000, and another survey on trees >10 cm in DBH was conducted in 2000 (Miller et al., 2004). Estimates of aboveground biomass (AGB) were then derived using allometric equations for Amazon forests (Rice et al., 2004; Chambers et al., 2004; Keller et al., 2001). Necromass (≥2 cm diameter) production was also measured approximately every 6 months in a 4.5-year period from November 2001 through February 2006 in a logged and undisturbed forest at km83 (Palace et al., 2008). Field measurements of ground disturbance in terms of number of felled trees, as well as areas disturbed by collateral and mechanical damages were also conducted at a similar site in Pará state along multitemporal sequences of postharvest regrowth of 0.5–3.5 years (Asner et al., 2004; Pereira et al., 2002).

Table 3 Cohort-level fractional damage fractions in different logging scenarios.

¹ DBH_min=50 cm. ² DBH${}_{{\max }_{\mathrm{infra}}}=\mathrm{30}$ cm. ³ Applied to the new patch generated by direct felling and collateral damage.

Table 2 provides a summary of stem density and basal area distribution across size classes at km83 based on the biomass survey data (Menton et al., 2011; de Sousa et al., 2011). To facilitate comparisons with simulations from FATES, we divided the inventory into early and late succession PFTs using a threshold of 0.7 g cm⁻³ for specific wood density, consistent with the definition of these PFTs in Table 1. As shown in Table 2, prior to the logging event in year 2000, this forest was composed of 399, 30, and 30 trees per hectare in size classes of 10–30, 30–50, and ≥50 cm respectively; following logging, the numbers were reduced to 396, 29, and 18 trees per hectare, losing ∼1.3 % of trees ≥10 cm in size. The changes in stem density (SD) were caused by different mechanisms for different size classes. The reduction in stem density of 2 ha⁻¹ in the ≥50 cm size class was caused by timber harvest directly, while the reductions of 3 and 1 ha⁻¹ in the 10–30 and 30–50 cm size classes were caused by collateral and mechanical damages. Corresponding to the loss of trees in logging operations, the basal area (BA) decreased from 3.9, 4.0, and 12.9 m² ha⁻¹ to 3.8, 3.9, and 10.8 m² ha⁻¹, and aboveground biomass (AGB) decreased from 3.8, 2.3, and 10.4 kg C m⁻² to 3.8, 2.2, and 8.7 kg C m⁻² in the 10–30, 30–50, and ≥50 cm size class, respectively.

2.4 Numerical experiments

In this study, the gap-filled meteorological forcing data for the Tapajós National Forest processed by Longo et al. (2019) are used to drive the CLM(FATES) model. Characteristics of the sites, including soil texture, vegetation cover fraction, and canopy height, were obtained from the LBA-Data Model Intercomparison Project (de Gonçalves et al., 2013). Specifically, soil at km67 contains 90 % clay and 2 % sand, while soil at km83 contains 80 % clay and 18 % sand. Both sites are covered by tropical evergreen forest at ∼98 % within their footprints, with the remaining 2 % assumed to be covered by bare soil. As discussed in Longo et al. (2018), who deployed the Ecosystem Demography Model version 2 at this site, soil texture and hence soil hydraulic parameters are highly variable even with the footprint of the same eddy covariance tower, and they could have significant impacts on not only water and energy simulations, but also simulated forest composition and carbon stocks and fluxes. Further, generic pedotransfer functions designed to capture temperate soils typically perform poorly in clay-rich Amazonian soils (Fisher et al., 2008; Tomasella and Hodnett, 1998). Because we focus on introducing the FATES logging, we leave for forthcoming studies the exploration of the sensitivity of the simulations to soil texture and other critical environmental factors.

Table 4 Comparison of energy fluxes (mean ± standard deviation) between eddy covariance tower measurements and FATES simulations.

CLM(FATES) was initialized using soil texture at km83 (i.e., 80 % clay and 18 % sand) from bare ground and spun up for 800 years until the carbon pools and forest structure (i.e., size distribution) and composition of PFTs reached equilibrium, by recycling the meteorological forcing at km67 (2001–2011) as the sites are close enough. The final states from spin-up were saved as the initial condition for follow-up simulations. An intact experiment was conducted by running the model over a period of 2001 to 2100 without logging by recycling the 2001–2011 forcing using the parameter set in Table 1. The atmospheric CO₂ concentration was assumed to be a constant of 367 ppm over the entire simulation period, consistent with the CO₂ levels during the logging treatment (Dlugokencky et al., 2018).

We specified an experimental logging event in FATES on 1 September 2001 (Table 3). It was reported by Figueira et al. (2008) that, following the reduced-impact-logging event in September 2001, 9 % of the trees greater than or equal to DBH_min=50 cm were harvested, with an associated collateral damage fraction of 0.009 for trees ≥ DBH_min. DBH${}_{{\max }_{\mathrm{infra}}}$ is set to be 30 cm, so that only a fraction of trees ≤30 cm are removed for building infrastructure (Feldpausch et al., 2005). This experiment is denoted as the RIL_low experiment in Table 2 and is the one that matches the actual logging practice at km83.

Figure 4 Simulated energy budget terms and leaf area indices in intact and logged forests compared to observations from km67 (left) and km83 (right) (Miller et al., 2011). The dashed vertical line indicates the timing of the logging event. The shaded areas in panels (a)–(f) are uncertainty estimates based on u∗-filter cutoff analyses in Miller et al. (2011).

We recognize that the harvest intensity in September 2001 at km83 was extremely low. Therefore, in order to study the impacts of different logging practices and harvest intensities, three additional logging experiments were conducted as listed in Table 3: conventional logging with high intensity (CL_high), conventional logging with low intensity (CL_low), and reduced impact logging with high intensity (RIL_high). The high-intensity logging doubled the direct-felling fraction in RIL_low and CL_low, as shown in the RIL_high and CL_high experiments. Compared to the RIL experiments, the CL experiments feature elevated collateral and mechanical damages as one would observe in such operations. All logging experiments were initialized from the spun-up state using site characteristics at km83 previously discussed and were conducted over the period of 2001–2100 by recycling meteorological forcing from 2001 to 2011.

3 Results and discussions

3.1 Simulated energy and water fluxes

Simulated monthly mean energy and water fluxes at the two sites are shown and compared to available observations in Fig. 4. The performances of the simulations closest to site conditions were compared to observations and summarized in Table 4 (i.e., intact for km67 and RIL_low for km83). The observed fluxes as well as their uncertainty ranges noted as Obs67 and Obs83 from the towers were obtained from Saleska et al. (2013), consistent with those in Miller et al. (2011). As shown in Table 4, the simulated mean (± standard deviation) latent heat (LH), sensible heat (SH), and net radiation (Rn) fluxes at km83 in RIL_low over the period of 2001–2003 are 90.2±10.1, 39.6±21.2, and 112.9±12.4 W m⁻², compared to tower-based observations of 101.6±8.0, 25.6±5.2, and 129.3±18.5 W m⁻². Therefore, the simulated and observed Bowen ratios are 0.35 and 0.20 at km83, respectively. This result suggests that, at an annual time step, the observed partitioning between LH and SH is reasonable, while the net radiation simulated by the model can be improved. At seasonal scales, even though net radiation is captured by CLM(FATES), the model does not adequately partition sensible and latent heat fluxes. This is particularly true for sensible heat fluxes as the model simulates large seasonal variabilities in SH when compared to observations at the site (i.e., standard deviations of monthly mean simulated SH are ∼21.2 W m⁻², while observations are ∼5.2 W m⁻²). As illustrated in Fig. 4c and d, the model significantly overestimates SH in the dry season (June–December), while it slightly underestimates SH in the wet season. It is worth mentioning that incomplete closure of the energy budget is common at eddy covariance towers (Wilson et al., 2002; Foken, 2008) and has been reported to be ∼87 % at the two sites (Saleska et al., 2003).

Figure 5 Simulated (a) aboveground biomass and (b) coarse woody debris in intact and logged forests in a 1-year period before or after the logging event in the four logging scenarios listed in Table 3. The observations (Obs_intact and Obs_logged) were derived from inventory (Menton et al., 2011; de Sousa et al., 2011).

Figure 4j shows the comparison between simulated and observed (Goulden et al., 2010) volumetric soil moisture content (m³ m⁻³) at the top 10 cm. This comparison reveals another model structural deficiency; that is, even though the model simulates higher soil moisture contents compared to observations (a feature generally attributable to the soil moisture retention curve), the transpiration beta factor, the down-regulating factor of transpiration from plants, fluctuates significantly over a wide range and can be as low as 0.3 in the dry season. In reality flux towers in the Amazon generally do not show severe moisture limitations in the dry season (Fisher et al., 2007). The lack of limitation is typically attributed to the plant's ability to extract soil moisture from deep soil layers, a phenomenon that is difficult to simulate using a classical beta function (Baker et al., 2008), and potentially is reconcilable using hydrodynamic representation of plant water uptake (Powell et al., 2013; Christoffersen et al., 2016) as in the final stages of incorporation into the FATES model. Consequently, the model simulates consistently low ET during dry seasons (Fig. 4e and f), while observations indicate that canopies are highly productive owing to adequate water supply to support transpiration and photosynthesis, which could further stimulate coordinated leaf growth with senescence during the dry season (Wu et al., 2016, 2017).

Figure 6 Simulated carbon fluxes in intact and logged forests compared to observed fluxes from km67 (left) and km83 (right). The dashed black vertical line indicates the timing of the logging event, while the dashed red horizontal line indicates estimated fluxes derived based on eddy covariance measurements and inventory (Miller et al., 2011). The shaded areas in panels (a)–(f) are uncertainty estimates based on u∗-filter cutoff analyses in Miller et al. (2011). Panels (g)–(i) show comparisons between annual fluxes as only annual estimates of these fluxes are available from Miller et al. (2011).

Figure 7 Trajectories of carbon pools in intact (left) and logged (right) forests. The dashed black vertical line indicates the timing of the logging event. The dashed red horizontal line indicates observed prelogging (left) and postlogging (right) inventories respectively (Menton et al., 2011; de Sousa et al., 2011).

3.2 Carbon budget as well as forest structure and composition in the intact forest

Figures 5–7 show simulated carbon pools and fluxes, which are tabulated in Table 5 as well. As shown in Fig. 5, prior to logging, the simulated aboveground biomass and necromass (CWD + litter) are 174 and 50 Mg C ha⁻¹, compared to 165 and 58.4 Mg C ha⁻¹ based on permanent plot measurements. The simulated carbon pools are generally lower than observations reported in Miller et al. (2011) but are within reasonable ranges, as errors associated with these estimates could be as high as 50 % due to issues related to sampling and allometric equations, as discussed in Keller et al. (2001). The lower biomass estimates are consistent with the finding of excessive soil moisture stress during the dry season and low leaf area index (LAI) in the model.

Table 5 Comparison of carbon budget terms between observation-based estimates^* and simulations at km83.

^* Source of observation-based estimates: Miller et al. (2011). Uncertainty in carbon fluxes (GPP, ER, NEE) are based on u∗-filter cutoff analyses described in the same paper.

Combining forest inventory and eddy covariance measurements, Miller et al. (2011) also provides estimates for net ecosystem exchange (NEE), gross primary production (GPP), net primary production (NPP), ecosystem respiration (ER), heterotrophic respiration (HR), and autotrophic respiration (AR). As shown in Table 5, the model simulates reasonable values in GPP (30.4 Mg C ha⁻² yr⁻¹) and ER (29.7 Mg C ha⁻² yr⁻¹), when compared to values estimated from the observations (32.6 Mg C ha⁻² yr⁻¹ for GPP and 31.9 Mg C ha⁻² yr⁻¹ for ER) in the intact forest. However, the model appears to overestimate NPP (13.5 Mg C ha⁻² yr⁻¹ as compared to the observation-based estimate of 9.5 Mg C ha⁻² yr⁻¹) and HR (12.8 Mg C ha⁻² yr⁻¹ as compared to the estimated value of 8.9 Mg C ha⁻² yr⁻¹), while it underestimates AR (16.8 Mg C ha⁻² yr⁻¹ as compared to the observation-based estimate of 23.1 Mg C ha⁻² yr⁻¹). Nevertheless, it is worth mentioning that we selected the specific parameter set to illustrate the capability of the model in capturing species composition and size structure, while the performance in capturing carbon balance is slightly compromised given the limited number of sensitivity tests performed.

Consistent with the carbon budget terms, Table 5 lists the simulated and observed values of stem density (ha⁻¹) in different size classes in terms of DBH. The model simulates 471 trees per hectare with DBHs greater than or equal to 10 cm in the intact forest, compared to 459 trees per hectare from the observed inventory. In terms of distribution across the DBH classes of 10–30, 30–50, and ≥50 cm, 339, 73, and 59 N ha⁻¹ of trees were simulated, while 399, 30, and 30 N ha⁻¹ were observed in the intact forest. In general, this version of FATES is able to reproduce the size structure and tree density in the tropics reasonably well. In addition to size distribution, by parametrizing early and late successional PFTs (Table 1), FATES is capable of simulating the coexistence of the two PFTs and therefore the PFT-specific trajectories of stem density, basal area, canopy, and understory mortality rates. We will discuss these in Sect. 3.4.

3.3 Effects of logging on water, energy, and carbon budgets

The response of energy and water budgets to different levels of logging disturbances is illustrated in Table 4 and Fig. 4. Following the logging event, the LAI is reduced proportionally to the logging intensities (−9 %, −17 %, −14 %, and −24 % for RL_low, RL_high, CL_low, and CL_high respectively in September 2001; see Fig. 4h). The leaf area index recovers within 3 years to its prelogging level or even to slightly higher levels as a result of the improved light environment following logging leading to changes in forest structure and composition (to be discussed in Sect. 3.4). In response to the changes in stem density and LAI, discernible differences are found in all energy budget terms. For example, less leaf area leads to reductions in LH (−0.4 %, −0.7 %, −0.6 %, −1.0 %) and increases in SH (0.6 %, 1.0 %, 0.8 %, and 2.0 %) proportional to the damage levels (i.e., RL_low, RL_high, CL_low, and CL_high) in the first 3 years following the logging event when compared to the control simulation. Energy budget responses scale with the level of damage, so that the biggest differences are detected in the CL_high scenario, followed by RIL_high, CL_low, and RIL_low. The difference in simulated water and energy fluxes between the RIL_low (i.e., the scenario that is the closest to the experimental logging event) and intact cases is the smallest, as the level of damage is the lowest among all scenarios.

Figure 8 Changes in total stem densities and the fractions of the early successional PFT in different size classes following a single logging event on 1 September 2001 at km83. The dashed black vertical line indicates the timing of the logging event, while the solid red line and the dashed cyan horizontal line indicate observed prelogging and postlogging inventories respectively (Menton et al., 2011; de Sousa et al., 2011).

As with LAI, the water and energy fluxes recover rapidly in 3–4 years following logging. Miller et al. (2011) compared observed sensible and latent heat fluxes between the control (km67) and logged sites (km83). They found that, in the first 3 years following logging, the between-site difference (i.e., logged–control) in LH reduced from 19.7±2.4 to 15.7±1.0 W m² and that in SH increased from 3.6±1.1 to 5.4±0.4 W m². When normalized by observed fluxes during the same periods at km83, these changes correspond to a −4 % reduction in LH and a 7 % increase in SH, compared to the −0.5 % and 4 % differences in LH and SH between RL_low and the control simulations. In general, both observations and our modeling results suggest that the impacts of reduced impact logging on energy fluxes are modest and that the energy and water fluxes can quickly recover to their prelogging conditions at the site.

Figure 9 Changes in basal area of the two PFTs in different size classes following a single logging event on 1 September 2001 at km83. The dashed black vertical line indicates the timing of the logging event, while the solid red line and the dashed cyan horizontal line indicate observed prelogging and postlogging inventories respectively (Menton et al., 2011; de Sousa et al., 2011). Note that for the size class 0–10 cm, observations are not available from the inventory.

Figures 6 and 7 show the impact of logging on carbon fluxes and pools at a monthly time step, and the corresponding annual fluxes and changes in carbon pools are summarized in Table 5. The logging disturbance leads to reductions in GPP, NPP, AR, and AGB, as well as increases in ER, NEE, HR, and CWD. The impacts of logging on the carbon budgets are also proportional to logging damage levels. Specifically, logging reduces the simulated AGB from 174 Mg C ha⁻¹ (intact) to 156 Mg C ha⁻¹ (RIL_low), 137 Mg C ha⁻¹ (RIL_high), 154 Mg C ha⁻¹ (CL_low), and 134 Mg C ha⁻¹ (CL_high), while it increases the simulated necromass pool (CWD + litter) from 50.0 Mg C ha⁻¹ in the intact case to 73 Mg C ha⁻¹ (RIL_low), 97 Mg C ha⁻¹ (RIL_high), 76 Mg C ha⁻¹ (CL_low), and 101 Mg C ha⁻¹ (CL_high). For the case closest to the experimental logging event (RIL_low), the changes in AGB and necromass from the intact case are −18 Mg C ha⁻¹ (10 %) and 23.0 Mg C ha⁻¹ (46 %), in comparison to observed changes of −22 Mg C ha⁻¹ in AGB (12 %) and 16 Mg C ha⁻¹ (27 %) in necromass from Miller et al. (2011), respectively. The magnitudes and directions of these changes are reasonable when compared to observations (i.e., decreases in GPP, ER, and AR following logging). On the other hand, the simulations indicate that the forest could be turned from a carbon sink (−0.69 Mg C ha⁻¹ yr⁻¹) to a larger carbon source in 1–5 years following logging, consistent with observations from the tower suggesting that the forest was a carbon sink or a modest carbon source ($-\mathrm{0.6}\pm \mathrm{0.8}$ Mg C ha⁻¹ yr⁻¹) prior to logging.

The recovery trajectories following logging are also shown in Figs. 6, 7, and Table 5. It takes more than 70 years for AGB to return to its prelogging levels, but the recovery of carbon fluxes such as GPP, NPP, and AR is much faster (i.e., within 5 years following logging). The initial recovery rates of AGB following logging are faster for high-intensity logging because of increased light reaching the forest floor, as indicated by the steeper slopes corresponding to the CL_high and RIL_high scenarios compared to those of CL_low and RIL_low (Fig. 9h). This finding is consistent with previous observational and modeling studies (Mazzei et al., 2010; Huang and Asner, 2010) in that the damage level determines the number of years required to recover the original AGB, and the AGB accumulation rates in recently logged forests are higher than those in intact forests. For example, by synthesizing data from 79 permanent plots at 10 sites across the Amazon basin, Ruttishauser et al. (2016) and Piponiot et al. (2018) show that it requires 12, 43, and 75 years for the forest to recover with initial losses of 10 %, 25 %, or 50 % in AGB. Corresponding to the changes in AGB, logging introduces a large amount of necromass to the forest floor, with the highest increases in the CL_high and RIL_high scenarios. As shown in Fig. 7d and Table 5, necromass and CWD pools return to the prelogging level in ∼15 years. Meanwhile, HR in RIL_low stays elevated in the 5 years following logging but converges to that from the intact simulation in ∼10 years, which is consistent with observations (Miller et al., 2011; Table 5).

Figure 10 Changes in mortality (5-year running average) of the (a) early and (b) late successional trees in different size classes following a single logging event on 1 September 2001. The dashed black vertical line indicates the timing of the logging event.

Figure 11 Changes in mortality (5-year running average) of the (a) canopy and (b) understory trees in different size classes following a single logging event on 1 September 2001. The dashed black vertical line indicates the timing of the logging event.

Table 6 The simulated stem density (N ha⁻¹) distribution at km83.

Table 7 The simulated basal area (m² ha⁻¹) distribution at km83.

3.4 Effects of logging on forest structure and composition

The capability of the CLM(FATES) model to simulate vegetation demographics, forest structure, and composition while simulating the water, energy, and carbon budgets simultaneously (Fisher et al., 2017) allows for the interrogation of the modeled impacts of alternative logging practices on the forest size structure. Table 6 shows the forest structure in terms of stem density distribution across size classes from the simulations compared to observations from the site, while Figs. 8 and 9 further break it down into early and late succession PFTs and size classes in terms of stem density and basal areas. As discussed in Sect. 2.2 and summarized in Table 3, the logging practices, reduced impact logging and conventional logging, differ in terms of preharvest planning and actual field operation to minimize collateral and mechanical damages, while the logging intensities (i.e., high and low) indicate the target direct-felling fractions. The corresponding outcomes of changes in forest structure in comparison to the intact forest, as simulated by FATES, are summarized in Tables 6 and 7. The conventional-logging scenarios (i.e., CL_high and CL_low) feature more losses in small trees less than 30 cm in DBH, when compared to the smaller reduction in stem density in size classes less than 30 cm in DBH in the reduced-impact-logging scenarios (i.e., RIL_high and RIL_low). Scenarios with different logging intensities (i.e., high and low) result in different direct-felling intensity. That is, the numbers of surviving large trees (DBH ≥ 30 cm) in RIL_low and CL_low are 117 ha⁻¹ and 115 ha⁻¹, but those in RIL_high and CL_high are 106 and 103 ha⁻¹.

In response to the improved light environment after the removal of large trees, early successional trees quickly establish and populate the tree-fall gaps following logging in 2–3 years as shown in Fig. 8a. Stem density in the <10 cm size classes is proportional to the damage levels (i.e., ranked as CL_high>RIL_high>CL_low>RIL_low), followed by a transition to late successional trees in later years when the canopy is closed again (Fig. 8b). Such a successional process is also evident in Fig. 9a and b in terms of basal areas. The number of early successional trees in the <10 cm size classes then slowly declines afterwards but is sustained throughout the simulation as a result of natural disturbances. Such a shift in the plant community towards light-demanding species following disturbances is consistent with observations reported in the literature (Baraloto et al., 2012; Both et al., 2018). Following regeneration in logging gaps, a fraction of trees win the competition within the 0–10 cm size classes and are promoted to the 10–30 cm size classes in about 10 years following the disturbances (Figs. 8d and 9d). Then a fraction of those trees subsequently enter the 30–50 cm size classes in 20–40 years following the disturbance (Figs. 8f and 9f) and so on through larger size classes afterwards (Figs. 8h and 9h). We note that, despite the goal of achieving a deterministic and smooth averaging across discrete stochastic disturbance events using the ecosystem demography approach (Moorcroft et al., 2001) in FATES, the successional process described above, as well as the total numbers of stems in each size bin, shows evidence of episodic and discrete waves of population change. These arise due to the required discretization of the continuous time-since-disturbance heterogeneity into patches, combined with the current maximum cap on the number of patches in FATES (10 per site).

As discussed in Sect. 2.4, the early successional trees have a high mortality (Fig. 10a, c, e, and g) compared to the mortality (Fig. 10b, d, f, and h) of late successional trees as expected given their higher background mortality rate. Their mortality also fluctuates at an equilibrium level because of the periodic gap dynamics due to natural disturbances, while the mortality of late successional trees remains stable. The mortality rates of canopy trees (Fig. 11a, c, e, and g) remain low and stable over the years for all size classes, indicating that canopy trees are not light-limited or water-stressed. In comparison, the mortality rate of small understory trees (Fig. 11b) shows a declining trend following logging, consistent with the decline in mortality of the small early successional tree (Fig. 10a). As the understory trees are promoted to larger size classes (Fig. 11d, f), their mortality rates stay high. It is evident that it is hard for the understory trees to be promoted to the largest size class (Fig. 11h); therefore the mortality cannot be calculated due to the lack of population.

4 Conclusion and discussions

In this study, we developed a selective logging module in FATES and parameterized the model to simulate different logging practices (conventional and reduced impact) with various intensities. This newly developed selective logging module is capable of mimicking the ecological, biophysical, and biogeochemical processes at a landscape level following a logging event in a lumped way by (1) specifying the timing and areal extent of a logging event; (2) calculating the fractions of trees that are damaged by direct felling, collateral damage, and infrastructure damage and adding these size-specific plant mortality types to FATES; (3) splitting the logged patch into disturbed and intact new patches; (4) applying the calculated survivorship to cohorts in the disturbed patch; and (5) transporting harvested logs off-site and adding the remaining necromass from damaged trees into coarse woody debris and litter pools.

We then applied FATES coupled to CLM to the Tapajós National Forest by conducting numerical experiments driven by observed meteorological forcing, and we benchmarked the simulations against long-term ecological and eddy covariance measurements. We demonstrated that the model is capable of simulating site-level water, energy, and carbon budgets, as well as forest structure and composition holistically, with responses consistent with those documented in the existing literature as follows:

1. The model captures perturbations on energy and water budget terms in response to different levels of logging disturbances. Our modeling results suggest that logging leads to reductions in canopy interception, canopy evaporation and transpiration, and elevated soil temperature and soil heat fluxes in magnitudes proportional to the damage levels.

2. The logging disturbance leads to reductions in GPP, NPP, AR, and AGB, as well as increases in ER, NEE, HR, and CWD. The initial impacts of logging on the carbon budget are also proportional to damage levels as results of different logging practices.

3. Following the logging event, simulated carbon fluxes such as GPP, NPP, and AR recover within 5 years, but it takes decades for AGB to return to its prelogging levels. Consistent with existing observation-based literature, initial recovery of AGB is faster when the logging intensity is higher in response to the improved light environment in the forest, but the time to full AGB recovery in higher intensity logging is longer.

4. Consistent with observations at Tapajós, the prescribed logging event introduces a large amount of necromass to the forest floor proportional to the damage level of the logging event, which returns to prelogging level in ∼15 years. Simulated HR in the low-damage reduced-impact-logging scenario stays elevated in the 5 years following logging and declines to be the same as the intact forest in ∼10 years.

5. The impacts of alternative logging practices on forest structure and composition were assessed by parameterizing cohort-specific mortality corresponding to direct felling, collateral damage, mechanical damage in the logging module to represent different logging practices (i.e., conventional logging and reduced impact logging), and intensity (i.e., high and low). In all scenarios, the improved light environment after the removal of large trees facilitates the establishment and growth of early successional trees in the 0–10 cm DBH size class proportional to the damage levels in the first 2–3 years. Thereafter there is a transition to late successional trees in later years when the canopy is closed. The number of early successional trees then slowly declines but is sustained throughout the simulation as a result of natural disturbances.

Given that the representation of gas exchange processes is related to, but also somewhat independent of, the representation of ecosystem demography, FATES shows great potential in its capability to capture ecosystem successional processes in terms of gap-phase regeneration; competition among light-demanding and shade-tolerant species following disturbance; and responses of energy, water, and carbon budget components to disturbances. The model projections suggest that while most degraded forests rapidly recover energy fluxes, the recovery times for carbon stocks, forest size structure, and forest composition are much longer. The recovery trajectories are highly dependent on logging intensity and practices, the difference between which can be directly simulated by the model. Consistent with field studies, we find through numerical experiments that reduced impact logging leads to more rapid recovery of the water, energy, and carbon cycles, allowing forest structure and composition to recover to their prelogging levels in a shorter time frame.

5 Future work

Currently, the selective logging module can only simulate single logging events. We also assumed that for a site such as km83, once logging is activated, trees will be harvested from all patches. For regional-scale applications, it will be crucial to represent forest degradation as a result of logging, fire, and fragmentation and their combinations that could repeat over a period. Therefore, structural changes in FATES have been made by adding prognostic variables to track disturbance histories associated with fire, logging, and transitions among land use types. The model also needs to include the dead tree pool (snags and standing dead wood) as harvest operations (especially thinning) can lead to live tree death from machine damage and windthrow. This will be more important for using FATES in temperate, coniferous systems, and the varied biogeochemical legacy of standing versus downed wood is important (Edburg et al., 2011, 2012). To better understand how nutrient limitation or enhancement (e.g., via deposition or fertilization) can affect the ecosystem dynamics, a nutrient-enabled version of FATES is also under testing and will shed more lights on how biogeochemical cycling could impact vegetation dynamics once available. Nevertheless, this study lays the foundation to simulate land use change and forest degradation in FATES, leading the way to direct representation of forest management practices and regeneration in Earth system models.

We also acknowledge that as a model development study, we applied the model to a site using a single set of parameter values, and therefore we ignored the uncertainty associated with model parameters. Nevertheless, the sensitivity study in the Supplement shows that the model parameters can be calibrated with a good benchmarking dataset with various aspects of ecosystem observations. For example, Koven et al. (2020) demonstrated a joint team effort of modelers and field observers toward building field-based benchmarks from Barro Colorado Island, Panama, and a parameter sensitivity test platform for physiological and ecosystem dynamics using FATES. We expect to see more of such efforts to better constrain the model in future studies.