Investigating the Global Biogeophysical Impact of Area and Mass Based Wood Harvest in a Vegetation Demography Model

Abstract

Wood harvesting alters land surface properties and energy redistribution, but there is a lack of studies estimating these changes on a global scale. We coupled a vegetation demographic model, the Functionally Assembled Terrestrial Ecosystem Simulator, with the E3SM land model to perform offline model simulation to investigate the land biogeophysical responses, including canopy coverage, leaf area index, albedo, surface roughness length, and energy fluxes, to historical wood harvest on the global scale. In this study, we found 50% less harvested carbon (C) when choosing the area‐based harvest rate as driving data that has not been spatially harmonized, compared to reharmonized mass‐based harvesting. By considering the uncertainty from reconstruction of historical wood harvest time series and the choice of wood harvest approach in the model, continuous wood harvest (1850–2015) results in 5%–10% of canopy coverage loss, contributing 0.5%–1% increase of albedo over disturbed land, which is much stronger than a non‐demographic land surface model. Changes in energy flux from the wood harvest are negligible (<1%), but the responses of land surface properties vary (up to 30%) due to differences in model structure between the single canopy, sun‐shade leaf model and vegetation demographic model.

Introduction

The global demand for woody biomass as raw material for production and traditional bioenergy has driven a continual rate of wood harvest that has modified the regional landscape and ecosystem function, contributing CO₂ emissions into the atmosphere and causing longer term climate impact. Deforestation, wood harvest or logging activities contributed an average global emission of 1.9 GtC yr⁻¹ reported by the most recent global carbon budget (Friedlingstein et al., 2022).

Forest wood harvest, including clear-cut logging, selective logging, thinning, and forest management, is one of the major land use activities which disrupts the global land cover and changes the physical properties of the land surface, including albedo, canopy coverage (the fraction of land area occluded by canopy), surface leaf area and roughness. Biogeophysically, these changes influence the ratio of transpiration versus evaporation and alter the surface energy balance and near surface air temperature, and modulate boundary layer dynamics and therefore cloud initiation and formation at local to regional scales (Bonan, 2008; Dickinson & Kennedy, 1992; Hoffmann & Jackson, 2000; Pongratz et al., 2021). Wood harvesting is a major land management activity that alters carbon sequestration potential and forest carbon (C) stock at regional (Calle et al., 2016; Canadell et al., 2009; Huang & Asner, 2010; Law et al., 2018; Pearson et al., 2014) and global scales (Friedlingstein et al., 2022; Houghton & Nassikas, 2018; Shevliakova et al., 2013; Stocker et al., 2014; Zeng et al., 2013), and is estimated as one of the largest contributors of anthropogenic CO₂ emission (Peng et al., 2023). Plenty of site level studies have investigated the biogeophysical effects from wood harvest in the recent decades, including but not limited to: the impact from changes in plant structure, age and post-harvest dynamics on albedo (Bourque et al., 1995; Halim et al., 2019; Kalliokoski et al., 2020; Kellomäki et al., 2021; Kuusinen et al., 2014), and changes in leaf mass/area or forest area (Fourrier et al., 2013; Seedre & Chen, 2010; Strukelj et al., 2015; Wagner et al., 2011), surface roughness length (Maurer et al., 2013; Nakai et al., 2008; Yuan et al., 2021), surface energy budget (Duveiller et al., 2018; Korkiakoski et al., 2019; Leppä et al., 2020; Mamkin et al., 2016; Williams et al., 2014) and surface air temperature (Liao et al., 2020; Su et al., 2019). The global impact of wood harvest on surface energy balance has received less attention compared to other land use change (LUC) activities (De Hertog et al., 2023; Luyssaert et al., 2014; Pongratz et al., 2021), with site level studies finding divergent conclusions due to differences in site specific environmental conditions and management practices. Wood harvest activity disturbs the land and brings post-harvest ecosystem succession and recovery processes, including the formation of forest gaps and regrowth as secondary forest, and can force the disturbed forest toward shorter-statured and younger stands that can cause extended legacy effects on land surface properties and energy fluxes under high harvest intensity (Liu et al., 2011). It can take 20 years (or longer depending on ecosystem type) for an ecosystem to reach the carbon compensation point, when initial losses of carbon stock are compensated through regrowth (Covington, 1981; Peichl et al., 2023), thus short-term field experiments can miss full recovery processes that can be estimated with land surface models.

While models are useful tools to investigate long-term impacts of continuous logging and global biogeophysical impact, land surface models often lack representation of forest heterogeneity (e.g., size-structure, and species-specific composition and disturbance dynamics) leading to relatively small model estimates of biogeophysical effects due to wood harvest (De Hertog et al., 2023). Studies show that the change of forest structure, including the loss of canopy layer area, the formation of forest gaps, and the reduced leaf area of the regrown secondary forest disturb land surface biogeophysical properties to varying degrees at different geographical locations (Alibakhshi et al., 2020; Otto et al., 2014). The classic big-leaf model approach and sun-shade single-canopy model can simulate the change of the leaf area index (LAI) induced by logging activities but ignore the canopy structure dynamics and formation of gaps, therefore are not sufficient to represent the impact on land biogeophysical properties. In addition, land surface models with no plant successional dynamics to represent secondary forest tend to either ignore or merge the wood harvest from secondary mature and young forest with primary forest (Huang et al., 2020; Lawrence et al., 2012), which also causes the biased estimation of biogeophysical impact.

Accounting for both the correct harvested area and correct C harvest and emission from wood harvest is challenging. The widely accepted approach is to drive land surface models with certain globally gridded-resolved wood harvest rates in units of area harvest per time, where the harvest area is based on using fixed forest C densities so that the integrated harvest amount aligns with reported harvest amounts at national or sub-national inventories. However, the assumed forest carbon stocks are generally not consistent with the biomass predicted by the land surface models themselves (Lawrence et al., 2012; Meiyappan et al., 2015), thus the total carbon harvested does not equal the original inventory data. An alternative is to drive the land models with the harvested carbon mass and then use the model-predicted carbon stocks to calculate the area of harvest required to generate this mass. The area-based and mass-based harvest rates will generally not be consistent with each other. Mass-based harvest rates allow the model to more accurately account for the C emission from wood harvest, though area-based harvest rates are easier to apply, and may have more accurate estimates of harvest area if model predictions of biomass are not accurate. With the spread of harvest rates from different data sources and approaches being applied to models (Lawrence et al., 2012), the calculated biogeophysical effects from models can also diverge.

In this work, we use a vegetation demography model (E3SM-FATES, Fisher et al., 2018) which includes plant demography and dynamic plant competition coupled with a land surface model to examine the global biogeophysical impact from wood harvest. E3SM-FATES comprises the Energy Exascale Earth system model (E3SM) Land Model (ELM) coupled with the Functionally Assembled Terrestrial Ecosystem Simulator (FATES), called ELM-FATES. Our objectives are to evaluate (a) the overall ELM-FATES simulated biogeophysical and plant physiological properties and biogeophysical impact caused by clear-cut wood harvest and (b) the related uncertainty caused by different wood harvest approaches (area-based and mass-based) and different historical reconstruction of wood harvest amounts and (c) the major hypothesis on how the inclusion of demographic processes for simulating canopy coverage change amplifies the magnitude of modeled biogeophysical responses to wood harvest activities. We hypothesize that continuous logging (or wood harvest) can lead to a persistent reduction of the canopy coverage, thus modifying LAI and surface roughness length to affect surface energy redistribution after accounting for the vegetation demographic information instead of the single-canopy parameterization. As a consequence of the more accurate representation of post-harvest vegetation successional dynamics, the vegetation demography model FATES generates forest gap and amplifies the response of land surface albedo and energy fluxes to wood harvesting compared to the single-canopy sun-shade model that lacks gap-phase disturbance and successional processes (Figure 1). This study focuses on the impact on land surface properties and energy fluxes without accounting for the feedback from the atmosphere, that is, the non-local biogeophysical effects (Pongratz et al., 2021; Winckler et al., 2017) are not considered.

[IMAGE OMITTED. SEE PDF]

Methods

Model Description

Functionally Assembled Terrestrial Ecosystem Simulator (FATES)

FATES (FATES Development Team, 2024; Fisher et al., 2015; Holm et al., 2020; Koven et al., 2020) is a demographic and dynamic vegetation model simulating dynamic recruitment, growth, and mortality as functions of plant physiology, phenology and plant competition under climate and environmental forces. Recent works have shown the capability of FATES to reproduce energy fluxes, plant physiology, phenology and demographic and plant functional dynamics for tropical (Holm et al., 2020; Koven et al., 2020), temperate forests (Buotte et al., 2021), cold-regions (Lambert et al., 2023), and global application of FATES (Needham et al., 2025). Full, detailed documentation on approaches and equations that govern FATES processes can be found in the model technical documentation (). FATES uses the cohort concept (Moorcroft et al., 2001) as the basic dynamic element to represent vegetation with the same plant functional type, size class and age class. At the leaf level, FATES keep the same separated sun-shade leaf fractions as ELM, with the addition of the canopy structure being multi-layered for biogeophysical calculations including photosynthesis, canopy energy and water fluxes. In addition, FATES also adopts the perfect plasticity approximation (Purves et al., 2008) to account for light penetration through a multi-layered forest canopy, based on cohort height (Fisher et al., 2015). FATES applies allometric equations to govern the allocation of assimilated carbon and growth of leaf area, coarse and fine root, aboveground biomass and height of plant cohorts and promotes and demotes upgrades/downgrades the canopy positioning of cohorts following the corresponding “target biomass” based on plant's diameter.

Individual plants are represented within a cohort as an average specimen of a given plant type and size (height and diameter). Plant cohorts are recruited, grow and diminish in abundance via mortality as they compete with other cohorts for resources and space in discrete patches of area that are defined by their time since disturbance. In current FATES we account for three different types of disturbances: tree fall (natural), fire, and logging. On the daily time step, disturbed cohorts are transferred to a new patch, and patches sharing similar age and plant size classes after the regrowth from disturbances will be merged. Patch dynamics allow FATES to track the age, area, and land use class of each patch after disturbance, which enables distinct wood harvest in primary, secondary mature, and secondary young forest (Figure 1).

Wood harvest in FATES can be either selective logging that allows the younger and understory plants to survive based on a prescribed criteria (Huang et al., 2020) or clear cut that removes all plants. In our study we only applied the clearcut option for two reasons: (a) clear-cut is the most widely used silvicultural treatment globally (Shorohova et al., 2019) and (b) we lack globally harmonized forest management data to drive the selective logging module. Wood harvest, as the area fraction of a cell per day, is one type of anthropogenic disturbance and is determined as the minimum between model derived forest inventory supply and harvest rate demand from a user-input forcing data set. A new secondary patch with either the remaining biomass (selective logging) or bare ground (clearcut) is created and the donor patch inherits the biomass from the undisturbed fraction of a gridcell. The new patch is independent of the donor patch, with separate records of all physiological, biogeochemical and demographic states.

A unique feature of FATES is the representation of size-structured vegetation that extends the linkages between vertical canopy structure and surface energy calculations via coupling to the host land model. Canopy radiative transfer is scaled to multiple layers to account for the variation of absorbed and reflected radiation in the overstory and understory (Bonan, 2008). Albedo at the top canopy layer is here calculated based on the Norman radiation scheme (Norman, 1979) and integrates the impact from non-vegetated bare soil for an open canopy forest. The displacement height and surface roughness are proportional to the canopy area weighted tree height across all the patches within the grid cell, which is governed by the allometric growth curve and plant population residing within patches. Stomatal conductance (Ball-Berry model used here, Collatz et al., 1991) is calculated at multiple layers and then integrated into a single value as a boundary condition to be passed to the ELM host land model.

FATES supports multiple different model settings to enable reduced model complexities for parameter calibration and/or idealized experiments that can address user-specific objectives. For example, the new reduced complexity settings can fix one or several key parameters or processes with either observations or testable constraints. The complexity of FATES calculation can be reduced by activating one or more of these modes. We use the “fixed biogeography” with “no competition” mode which constrains the PFT's geographical distribution in FATES with a user-defined plant biogeographical map. The no competition mode halts the competition between different PFTs for the land area, but intra-PFT competition can still occur. By turning off the prognostic biogeography, this configuration allows us to reduce dynamic interactions among different PFTs, thus improving the robustness of the simulated biogeophysical states and biomass. Due to the fixed area fraction for each PFT, each patch is assigned a unique PFT type, thus logging disturbance only accounts for patches with tree PFTs and the clear-cut will remove all the plant biomass. Although cross-PFT competition is deactivated, most of key demographic features are still active, including: (a) the multiple-layer canopy for radiative transfer, plant photosynthesis, respiration and stomata control, (b) recruitment and mortality to govern the population and age and (c) plant growth along an allometric curve to prognostically calculate the plant height, canopy coverage and diameter at breast height.

E3SM Land Model v2 (ELMv2)

A host land surface model is required to couple with FATES to simulate hydrologic and soil biogeochemical processes, surface energy redistribution in vegetation, snow, and soil, and pass the supportive boundary conditions to FATES. Here we use ELMv2, which describes energy and water dynamics under different vegetation properties determined through FATES. For example, the calculation of surface energy fluxes in ELM, sensible heat flux (SHF) and latent heat flux (LHF), are governed by aerodynamic resistance and the corresponding temperature and humidity gradient between the atmosphere and land surface, which is a function of roughness length determined by canopy height. The calculation of momentum and evapotranspiration is handled by ELMv2 and does not account for the interlayer exchange within plant canopy.

In addition to the lack of plant vertical structure and demography, ELMv2 is different from ELM-FATES as follows: (a) the calculation of surface radiation fluxes from ELMv2 is based on sun-shade canopy assumptions, with both sunlit and shaded parts covering the whole corresponding PFT fractional area; (b) canopy radiative transfer is based on a two-stream approach instead of the default Norman radiation scheme in FATES (however the latest version of FATES includes the two-stream radiation approach); (c) wood harvest activity represented ELMv2 is one type of litter loss from plant and no secondary forest information is accounted for (Figure 1). Apart from these different processes, no changes were made in ELMv2, except for running the model in carbon-only mode without nitrogen-phosphorus cycling.

Implementation of Mass-Based Harvest Approach

In ELM and other land surface models that participated in CMIP6 (Curasi et al., 2024; Lawrence et al., 2012; Lindeskog et al., 2021; Littleton et al., 2020; Nabel et al., 2020; Yue et al., 2018) the standard wood harvest approach is based on the rate of harvested area ($${H}_{A}$$). The model is driven by annual harvest rate input and removes the biomass C equivalent to the biomass within the harvested area calculated from the model. Here, we added a mass-based wood harvest approach as an alternative method to improve the representation of harvested carbon in the model. In this approach, the removal of biomass is directly determined by the annual demand of harvested C ($${H}_{c}$$), while the disturbance rate is equivalent to the fractional harvested area per year ($${H}_{A}$$) and is calculated as follows: 1 $${H}_{A}=\frac{{H}_{c}}{{C}_{\text{avail}}}$$

Here $${C}_{\text{avail}}$$ (kgC) is the total forest C available for harvest within a gridcell and determined based on a list of criteria to constrain the age and size of plants permitted for logging (Huang et al., 2020). In this study, we do not account for forestry management, thus $${C}_{\text{avail}}$$ equals to the total aboveground forest C. Global Time Series (TS) of $${H}_{c}$$ is required as an input to drive the model $${C}_{\text{avail}}$$.

Due to the mismatch of $${C}_{\text{avail}}$$ between FATES and the data source used for calculations, logging disturbance may exceed 100%/yr in certain cases, indicating the model's available carbon supply to be less than the demand from the data source. To prevent complete depletion of forests under harvest stress—an unlikely scenario in the real world—we opt to leave the land undisturbed in such cases. This approach allows for forest regrowth, but leads to smaller harvested carbon than expected. To mitigate the deficiency of harvested carbon in our simulation, we designed a re-harmonization algorithm to re-distribute the spatial distribution of $${H}_{c}$$ following the method in Section 2.1.4.

Re-Harmonization of Wood Harvest Rate Spatial Pattern With FATES Forest C Inventory

The overall strategy of re-harmonization is to obtain the deficit between wood harvest demand from driving data and the FATES-calculated $${C}_{\text{avail}}$$, which we call harvest debt, and reduce the harvest debt by redistributing $${H}_{c}$$ to nearby gridcells and maintaining the original spatial distribution of wood harvest rate forcing data as much as possible (Figure 2).

[IMAGE OMITTED. SEE PDF]

Based on the land disturbance history, forest is separated into primary (no disturbance), secondary mature (last disturbance since 94 years ago) and secondary young. To re-harmonize primary and secondary harvest rates respectively, $${C}_{\text{avail}}$$ of primary and secondary forest are required and can be calculated from model simulation $$f()$$: 2 $${C}_{\text{prim},\text{avail}}=f\left({H}_{c,\text{prim}}=0,{\text{CO}}_{2}=284.7\,\text{ppm}\right)$$

Here we aim to obtain a low value from spinup simulation under fixed pre-industrial CO₂ concentration to exclude the increment of biomass due to CO₂ fertilization. The reason is to redistribute $${H}_{c,\text{prim}}$$ using a low $${C}_{\text{prim},\text{avail}}$$ map, thus $${H}_{c,\text{prim}}$$ is more likely to be lower than actual $${C}_{\text{avail}}$$, when forced by the historical CO₂ TS.

Wood harvest from primary forest permanently reduces primary forest area without forest regrowth, thus $${C}_{\text{prim},\text{avail}}$$ is a good estimation of the total wood product supply over historical period. We accumulate the harvest rate of the whole historical period (850–2015) from forcing data to obtain total wood product demand, then calculate harvest debt, and perform adjustment once for the whole TS. For each grid cell with harvest debt, we initialize an increasing search radius and transfer the debt to nearby grid cells with harvest surplus (i.e., $${C}_{\text{prim},\text{avail}}-{H}_{c,\text{prim}} > 0$$) and end the search when the debt is fulfilled. The amount transferred is proportional to the harvest surplus. $${H}_{c,\text{prim},\text{corr}}$$ are used as forcing in a transient FATES simulation restart from spinup in 1700 to calculate the secondary forest inventory. We only consider secondary harvest from year 1800, allowing the first 100 year simulation to provide an initial secondary forest area: 3 $${C}_{\text{secd},\text{avail}}[\text{year}]=f\left({H}_{c,\text{prim}}={H}_{c,\text{prim},\text{corr}}[\text{year}],{H}_{c,\text{secd}}=0,{\text{CO}}_{2}=284.7\,\text{ppm}\right)$$

For secondary forest harvest rate, re-harmonization is performed year by year to account for the secondary forest regrowth, which is different from primary forest. $${C}_{\text{secd},\text{avail}}[\text{year}]$$ is split into secondary mature forest inventory, defined as carbon from forest patches older than 94 years, which is the global mean age of secondary land at 1900s (Hurtt et al., 2011), and the rest is treated as secondary young forest inventory. Both the adjustment of secondary mature and young forest harvest rates follow the same algorithm as primary harvest rate does.

Specifically for this study, Land Use Harmonization v2 (LUH2) historical reconstructed wood harvest rates (see Section 2.2.2) for primary, secondary young and secondary mature forest are used as the original harvest rate forcing. We merge the primary non-forest to primary forest and secondary non-forest to secondary young forest based on its definition in the LUH2 data set that small and young trees are treated as non-forest. Harvest debt, especially from secondary young forest, still exists (Figure S1 in Supporting Information S1) using the re-harmonized harvest rates forcing due to the different settings between the FATES simulation without real time secondary forest harvest and the historical FATES simulation, leading to inconsistent $${C}_{\text{secd},\text{avail}}[\text{year}]$$ and age information used in the algorithm. Coupling the re-harmonization code in FATES to account for the real time harvest and spatial redistribution could provide precise $${C}_{\text{secd},\text{avail}}[\text{year}]$$ values to resolve this issue, but is difficult to implement due to the parallelization design of ELM, which is beyond the scope of this study.

Data Sets

Global Soil Wetness Project Phase 3 (GSWP3) Atmospheric Forcing

In this study we run an offline ELM (I1850ELMFATES, Golaz et al., 2022) driven by GSWP v3 (Kim, 2017) climate forcing. GSWP v3 is a climate forcing product from 1901 to 2014, dynamically downscaled from the Twentieth Century Reanalysis (20CR) product (Compo et al., 2011) and bias-corrected using the following station level observations: Global Precipitation Climatology Center Version 6 (for precipitation), Climatic Research Unit TS 3.21 (for air temperature), and Surface Radiation Budget project (SRB, for solar shortwave and atmospheric downward longwave radiation). Seven GSWP v3 climate forcing variables at the 2-m level are used: solar shortwave radiation, atmospheric longwave radiation, precipitation, air temperature, air pressure, specific humidity and wind speed.

Land Use Harmonization v2 (LUH2) Historical Reconstructed Land Cover and Wood Harvest Rates

Three different annual mass-based and area-based wood harvest rates and land cover type data sets are collected to drive four different wood harvest scenarios in the model (baseline, high, low, and re-harmonized scenarios, Hurtt et al., 2020). The land cover type data is transferred into ELM PFTs through the land use translator (Di Vittorio et al., 2020; Lawrence et al., 2012). LUH2 wood harvest rates are generated based on national statistics for 199 countries from the History database of the Global Environment v3.2 (HYDE 3.2) population data before 1961 and based on FAO national wood harvest volume data after 1961. The TS of mass-based harvest rates are traced back to 850 using population changes and the country-level per capita wood harvest rate after accounting the linear transition of energy source from fuel wood to fossil fuel between 1800 and 1920. Spatial allocation of mass-based wood harvest rates from country level is based on the proximity to the gridcell with “significant human presence” by assuming an equivalent proximity to transportation infrastructure or local markets. The available forest biomass is calculated based on the MIAMI-LU ecosystem model (Hurtt et al., 2011) and is used to determine the actual harvested C amount in each gridcell and calculate the area-based harvest rates. The major difference among scenarios is the high harvest rate prior to 1920 in the high-end scenario, resulting in a 10 times larger secondary forest area in 1850.

Other Time Invariant Auxiliary Data Sets and Observations for Model Evaluations

Auxiliary data sets used in this study include: (a) Historical CO₂ TS from 1700, (b) soil texture and soil color, and (c) initial PFT coverage maps. Soil texture is based on International Geosphere-Biosphere Programme soil data set products and soil color is derived from Lawrence and Chase (2007). For the initial PFT coverage map, we set a stagnant crop and pasture distribution at 1850 level based on the corresponding LUH2 land cover data. The reason for choosing the 1850 distribution is to maximize the forest coverage to avoid large deficiency of forest inventory for wood harvest, especially during the simulation period before 1950 when the wood harvesting per year rose drastically.

We use multiple observational data products obtained from the International Land Model Benchmarking system data repository (Collier et al., 2018) for ELM-FATES model evaluation on a global scale. Observational products used here are: (a) biomass products: “Tropical” (JPL tropical Forest Biomass, Saatchi et al., 2011), GeoCarbon (Global forest live biomass carbon version 5.1, Saatchi et al., 2011), ESACCI (ESA CCI above-ground biomass product, S, Santoro & Cartus, 2021), Global Carbon (Saatchi et al., 2011), Thurner (Thurner et al., 2014), XuSaatchi-2021 (Xu et al., 2021); (b) albedo product: GEWEX (NASA GEWEX SRB radiation data set, Collier et al., 2018); (c) LAI products: AVHRR (Advanced Very High Resolution Radiometer observation, Collier et al., 2018), MODIS (monthly climatology MODIS total LAI, Collier et al., 2018), AVH15C1 (Collier et al., 2018); (d) gross primary productivity (GPP) products: Fluxcom (Jung et al., 2017), WECANN (Alemohammad et al., 2017); (e) SHF products: CLASS (Conserving Land-Atmosphere Synthesis Suite v1.1, Hobeichi et al., 2020), Fluxcom (Jung et al., 2017), WECANN (Alemohammad et al., 2017), and (f) ET products: MODIS (Collier et al., 2018), GLEAM (Martens et al., 2017). All data sets are upscaled to global 4 by 5° spatial resolution to match the grid used here.

Experiment Design and Model Output Analysis

We perform historical simulations (Shu, 2024) with the same simulation period as CMIP6 (1850–2015). The time steps of simulations are half-hourly for most physics processes (hydrology, energy balance, turbulent transport, radiation scattering, photosynthesis) and daily for plant dynamics (growth allocation, phenology, mortality, disturbance). In addition, we perform a spin-up from 1700 to 1850 to reconstruct the secondary forest inventory before 1850, which is a trade-off against the full simulation from 850 AD (Ma et al., 2022) due to the computational cost. None of the simulations described here account for fire disturbances and all follow these four steps:

• Spinup phase 1. CO₂ concentration is fixed at 284.7 ppm and no logging activities occur for 200 years. This step aims to build and force the forest biomass to reach a pre-industrial steady state.

• Spinup phase 2: This simulation is initialized with the land-surface state (vegetation demographics, soil carbon) from the end of Phase 1, and is driven by primary forest and non-forest harvesting rates from 1700–1799 to build potential secondary forest. Harvested C during 850 AD–1700 AD was accumulated and applied one time in 1700 to build the correct mature (patch age >94 years) secondary forest area.

• Spinup phase 3: This is initialized with the land-surface state from the end of Phase 2, and is driven by harvest rates from all harvest categories during 1800–1849 to build young secondary forest and diminish the legacy effect from the large logging disturbance in 1700.

• Historical phase (1850–2015): This is initialized with the land-surface state from the end of Phase 3, and is forced by historical atmospheric forcings and harvest rate TS, including a transient (changing) atmospheric CO₂ signal. For simplicity, we apply clear-cut once every year at the end of January.

We design experiments by varying two factors contributing to uncertainties one at a time: (a) two wood harvest approaches (area based and mass-based) and (b) three harvest rate forcings (historical baseline, historical low, and historical high). Under each harvest rate forcing scenario, we perform another reference run with no harvest forcing for the further analysis. To simplify the description of each case, we annotate all experiments using three letters describing the forcing and model followed by one letter describing harvest approach (Table 1). The historical mass-based wood harvest simulation driven by re-harmonized baseline rates (“har_c”) is defined as the base case. In addition, we run one case using ELM without coupling to FATES (ELM-only) to represent the traditional sun-shade canopy model.

Table 1 Annotations and Settings of All Experiments Performed in This Study

Case name	Land surface model	Harvest rate driving data	Wood harvest approach	Purposes
har_c	ELM-FATES	Reharmonized historical baseline	Mass-based harvest	Base case for evaluation of wood harvest influence
his_c	ELM-FATES	LUH2 historical baseline	Mass-based harvest	Cases for evaluation of uncertainty from LUH2 wood harvest driving data and harvest approach (area-based vs. mass-based)
his_n	ELM-FATES	LUH2 historical baseline	No harvest
his_a	ELM-FATES	LUH2 historical baseline	Area-based harvest
hlo_c	ELM-FATES	LUH2 historical low	Mass-based harvest
hlo_n	ELM-FATES	LUH2 historical low	No harvest
hlo_a	ELM-FATES	LUH2 historical low	Area-based harvest
hhi_c	ELM-FATES	LUH2 historical high	Mass-based harvest
hhi_n	ELM-FATES	LUH2 historical high	No harvest
hhi_a	ELM-FATES	LUH2 historical high	Area-based harvest
elm_a	ELM only	LUH2 historical baseline	Area-based harvest	Cases for evaluation of sensitivity caused by model structure
elm_n	ELM only	LUH2 historical baseline	No harvest

We focus on regional outputs because globally integrated measures of wood harvest impact will be less significant due to spatially scattered and temporally non-persistent global wood harvest signals (Noblet-Ducoudré et al., 2012). We explore only the regions that have undergone wood harvest in the historical simulation period. For each simulation case, grid cells with no harvest activities will be excluded from the analysis of the wood harvest impact. We also perform aggregation separately for two zones: tropical (30S–30N) and extratropical (<30S and >30N), which have different biogeophysical responses to LUC (Pongratz et al., 2021).

The overall wood harvest impact from the year wood harvesting started (1700) on a variable ($$\mathit{{\increment}}{\text{var}}_{\text{acc}}$$) at year ($$i$$) can be estimated using the difference between the harvest case and control case: 4 $$\mathit{{\increment}}{\text{var}}_{\text{acc}}(i)={\text{var}}_{{\text{wt}}_{hrv}}(i)-{\text{var}}_{{\text{wot}}_{hrv}}(i)$$

Wood harvest impact contains consequences from both the direct loss of C and the recovery of C from regrowth, which contains legacy effects from the beginning year of simulation, thus we call it accumulated wood harvest impact. In addition, we separate the direct loss ($$\mathit{{\increment}}{\text{can}}_{\text{direct}}$$) and gain of canopy coverage from the regrowth for each year. The direct loss is calculated from the canopy coverage multiplied by harvest rate. The regrowth impact after wood harvest ($$\mathit{{\increment}}{\text{can}}_{\text{regrowth}}$$) is estimated by subtracting direct from the overall ($$\mathit{{\increment}}{\text{can}}_{\text{total}}$$) impact during the year ($$i$$): 5 $$\mathit{{\increment}}{\text{can}}_{\text{direct}}(i)={\text{can}}_{\text{jan}}(i)\times {H}_{A}(i)$$ 6 $$\mathit{{\increment}}{\text{can}}_{\text{total}}(i)={\text{can}}_{\text{jan}}(i+1)-{\text{can}}_{\text{jan}}(i)$$ 7 $$\mathit{{\increment}}{\text{can}}_{\text{regrowth}}(i)=\mathit{{\increment}}{\text{can}}_{\text{total}}(i)-\mathit{{\increment}}{\text{can}}_{\text{direct}}(i)$$

This method does not work well with regions with very tiny wood harvest rate, which have stronger climate variability than word harvest signals. Here we only apply the method to analyze the change of canopy coverage and filter out grid cells with a small harvested C (0.01 kgC m⁻² yr⁻¹) due to the relatively smaller signal than natural variability.

Results

Evaluation of ELM-FATES Simulation

We first examine the capability of ELM-FATES to estimate global scale vegetation carbon stock, plant productivity and energy fluxes by comparing the latitudinal mean of the har_c case to observational products from multiple sources in the 2000s (i.e., mean of 2001–2010). The aim is not to present a rigorous cell-level comparison due to the current model configuration excluding full LUC, which would require a more detailed model benchmarking exercise that is out of scope for this paper, but to show if ELM-FATES can capture a general latitudinal pattern of energy related variables, forest productivity, carbon storage with acceptable bias, thus adding confidence to our estimation of wood harvest impact.

Six variables are evaluated (Figure 3): forest biomass C, surface albedo, summer-time (i.e., June-July-August) LAI, GPP, SHF and evapotranspiration (ET). ELM-FATES estimated global total forest carbon of the 2000s (defined as the mean of 2001–2010) to be 675 PgC, which is at the high end but close to the ESA CCI product. By descending order, 631 PgC from ESA CCI, 454.7 PgC from GeoCarbon, 367.7 PgC from Saatchi et al. (2011), and the updated version of the same product by Xu et al. (2021) has a total live carbon of 381.7 PgC. All products vary in their definition (e.g., live biomass or total biomass) and approaches, the method to estimate belowground biomass, thus can only be provided as references. The comparison through latitudinal zones shows the total biomass at the high end in the southern midlatitudes and tropics and at the low end in the northern midlatitudes. Model simulation shows low biomass over northern high latitudes with up to 7 PgC/2.1 kgC m⁻² difference between FATES and the mean of observational products. However, high latitudes are not intensively harvested, contributing only 7.8% of total harvested C in the 2000s, thus introducing a relatively small bias to our results. Through checking the spatial distribution (Figure S5 in Supporting Information S1), FATES successfully simulates the hotspot areas with high vegetation carbon, including tropical rain forest, south-eastern North America and Asia and central Europe, except sub-arctic boreal forest which has a relatively low forest C.

[IMAGE OMITTED. SEE PDF]

Global GPP (Figure 3d) from ELM-FATES is 125.2 PgC m⁻² yr⁻¹, close to the global value of both observational products, Fluxcom (114.4 PgC m⁻² yr⁻¹) and WECANN (115.1 PgC m⁻² yr⁻¹). The mean GPP of each latitudinal zone also shows a good match of general pattern against both products. The Northern subtropical region has around 0.6 gC m⁻² d⁻¹ underestimation and the Southern subarctic overestimates GPP by about 1 gC m⁻² d⁻¹. Similar to the forest biomass, the bias of ELM-FATES estimation of GPP with respect to the mean of all observational products is lower than the difference between observational products. The spatial pattern of FATES GPP has a better match to forest C (Figure S8 in Supporting Information S1). FATES estimation is similar to WECANN rather than Fluxcom in the sense of spatial pattern.

We also examine the GPP-age relationship (Tang et al., 2014) against the observations from Fluxnet products (Besnard et al., 2018) to evaluate the secondary forest regrowth (Figure 4). We borrow the same NEP-age relationship originally developed to represent the temporal patterns of annual GPP-to-ER ratio in Besnard et al. (2018) and assume the GPP-age to have the same statistical relationship, which is valid with trees younger than 100 years (Tang et al., 2014). The available fluxnet sites used here are globally representative because they are located across a gradient of mean annual temperature from tropical to boreal regions. The scatter plot (Figure 4) shows similar GPP and age ranges between the global simulation and site level observations. The fitted spline has a low R-squared (0.05 for Fluxnet data and 0.24 for ELM-FATES) due to substantial contributions from other environmental factors to the spatial and temporal variability of GPP, but still can reflect the general age dependency of tree productivity, that is, regrowth ability. The fitted spline of ELM-FATES shows a close match to Fluxnet data for younger trees (<20 years) but a steeper and more continuous regrowth ability for older trees (>20 years) compared to Fluxnet data, partly caused by fewer boreal and temperate samples available at this age range.

[IMAGE OMITTED. SEE PDF]

Annual mean LAI (Figure 3c) is evaluated against remote sensing products to check the model representation of carbon use efficiency, allocation and plant phenology. The annual mean LAI from ELM-FATES shows a better match in both midlatitudes, with a bias less than 0.5 m² m⁻² compared to observational products. The peak LAI from the tropical region reaches 2.9 m² m⁻² from ELM-FATES estimation, where other products show the range between 3.6 and 4.1 m² m⁻². The spatial pattern of LAI (Figure S7 in Supporting Information S1) shows an overall lower value compared to other products, although the other three data products show large variability.

Land surface energy fluxes are key variables investigated in the remainder of this study for understanding the response from wood harvesting (Figures 3e and 3f). Due to the impact from snow cover and limited daylight length during winter, we only use albedo during summer (JJA). ELM-FATES estimation of summertime albedo (Figure 3b) ranges from 0.17 in the tropics to 0.68 around 60N and matches the latitudinal mean pattern from GEWEX SRB radiation observation while reflecting certain biases in both high latitudinal zones. Biases are more severe in the northern high latitudes due to the area of land snow coverage and timing of snow melt, which is a similar issue found in other benchmark studies (Collier et al., 2018; He et al., 2014). Both SHF and ET estimation show better matches with observations than albedo, partly due to the high accuracy of the prescribed atmospheric forcing used in the offline simulations. Global mean annual SHF from land is 32.3 W m⁻², which is within the range of three observational products: CLASS (25.2 W m⁻²), Fluxcom (29.2 W m⁻²) and WECANN (34.6 W m⁻²). For ET, the global annual mean of ELM-FATES estimation is 526.8 mm yr⁻¹, higher than the global annual mean of GLEAM (460.6 mm yr⁻¹) and MODIS (390.7 mm yr⁻¹) but closer to the GLEAM product. ET bias is largest in the tropical and southern midlatitude regions. The ET overestimation in southern midlatitudes can be attributed to transpiration from the high forest biomass estimation in ELM-FATES. The spatial pattern of the summertime albedo (Figure S6 in Supporting Information S1), sensible heat (Figure S9 in Supporting Information S1) and evapotranspiration (Figure S10 in Supporting Information S1) can illustrate correctly the general pattern compared to other data products. But overall the FATES land surface energy is slightly lower over the tropics and boreal regions, which may be linked to the relatively high albedo over these two regions.

Estimation of Global Secondary Forest Area and Harvested C

Large variation of estimated accumulated harvested C (1850–2015; Figure 5) and secondary forest areas in 2015 (Figure 6, Figure S4 in Supporting Information S1) from different simulation cases is a consequence of the uncertainty from harvest approaches and harvest rate driving products (Figure S2 in Supporting Information S1). The har_c case estimation of 1850–2015 accumulated harvested C (100 PgC) is the closest to the corresponding LUH2 data (119 PgC). Accumulated harvested C calculated from other ELM-FATES cases show 50% less (60 PgC) compared to LUH2 historical baseline (119 PgC) while ELM-only estimation of forest biomass is higher (157 PgC), indicating that harvested C is very sensitive to harvest approach and land model $${C}_{\text{avail}}$$. The mismatch of $${C}_{\text{avail}}$$ between FATES and MIAMI-LU (used for LUH2) exceeds 10 kgC m⁻² over tropical regions (Figure S18 in Supporting Information S1), thus leading to a large difference in spatial pattern of wood harvest rates before and after re-harmonization (Figure S2 in Supporting Information S1). Besides the har_c case, all area-based harvest cases show lower harvested C compared to the corresponding mass-based harvest cases, showing the area based approach would most likely lead to underestimation of C emission from wood harvest. The variation of harvested C among the three harvest level scenarios is much smaller than the wood harvest approach (area-based vs. mass-based) based on model structure. It is important to note that the mass-based high harvest rate scenario (hhi_c) has lower actual harvested C and secondary land area than the other two scenarios despite the higher prescribed harvest rate. The major cause is that hhi_c harvest rates before 1850 are much higher than available forest carbon, leading to less harvest loss of primary forest area and C due to assuming that no harvest is applied if $${C}_{\text{avail}}$$ is less than the one-time harvest rate.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Based on the LUH2 historical baseline, secondary forest area increases from 0.5 in 1850 to 1.2 × 10⁷ km² in 2015 when accounting for only wood harvest activity. However, the uncertainty in harvest-driven secondary forest area ranges from 0.4 to 1.4 × 10⁷ km² as determined by the difference between hlo_a and hhi_a cases at 1850. The secondary forest area from the his_a case generally matches the LUH2 baseline, illustrating that FATES-ELM simulation captures the correct harvest area described by forcings, while the small bias is caused by the gridcell where forest areas are exhausted. The mass-based harvest method shows a larger increase of secondary forest area compared to the area-based harvest method from the 1850s–2010s, a direct consequence from the lower biomass in FATES than the estimation from LUH2, leading to the simulated higher harvested areas than area-based approach and LUH2 baseline. The har_c case driven by the re-harmonized harvest rates shows the highest increase of secondary forest area among all cases, with a large fraction of tropical forest and boreal forest from Siberia disturbed compared to other cases (Figure S3 in Supporting Information S1). The tradeoff between meeting the correct harvest biomass demand and accurately representing the harvested area depends on biomass prediction capabilities and has implications for both biogeochemical and biogeophysical impacts on the Earth system.

Historical Biogeophysical Response to Wood Harvesting Estimated With ELM-FATES

We examined the changes of canopy coverage (Figure S5 in Supporting Information S1), LAI, albedo, surface roughness length (Figure 7) and surface energy fluxes (SHF and LHF, Figure 8) caused by wood harvest disturbance. Most variables show less than 1% of changes averaged over disturbed forest regions in the 1850s, but the accumulated impacts increase and reach up to more than 5% in the 2000s, and more than 60% of the impacted regions show significant change of canopy coverage compared to the 1850s (Figure S5 in Supporting Information S1). Wood harvest increases surface albedo to 0.5%–1% over disturbed forest land. Surface roughness length has the largest relative decrease among the variables because the basal-weighted canopy height decreases more than canopy coverage and LAI. We find a smaller increase of albedo caused by wood harvest than that estimated from LUC (Luyssaert et al., 2014). The response from extra-tropical and tropical regions share the same direction, but with differences in the magnitude. Temperate regions lose more forest area compared to the tropics, reflected by a larger reduction of canopy coverage and LAI. While tropical regions show weaker loss of canopy coverage, a larger reduction of roughness length compared to temperate forest is caused by the loss from tropical cohorts with larger size and height.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Compared with ELM, ELM-FATES has a wider distribution of energy flux responses to changes in surface roughness length due to harvest in intensely disturbed forest areas (Figure 9), even though global energy fluxes are relatively similar across cases (Figure 8). Loss of forest coverage raises the contribution of SHF instead of LHF due to the decrease of plant transpiration, though the magnitude is less than 0.1 W m⁻². Compared to ELM, the sensitivity of SHF to the change of roughness length is largely increased in ELM-FATES, described as the larger slope of the regressed trendline (Figure 9). ELM has a larger sensitivity of LHF to the change of roughness length, though the gross change of LHF from ELM-FATES is larger (Figure 8) because of the larger decrease of roughness length. Changes in the following variables by incorporating ELM-FATES can change the sensitivity of SHF and LHF to surface roughness length compared to ELM: (a) vegetation and land surface temperature/humidity, (b) displacement height and (c) different PFT fractions. Change of displacement height has the same direction in both models, thus it is likely the increase of vegetation temperature (Figure S19 in Supporting Information S1) and consequently the increase of leaf saturated vapor pressure leads to the increase of slopes, although PFT fraction would change their sensitivity and the further analysis is beyond the scope of this study.

[IMAGE OMITTED. SEE PDF]

The spatial patterns of six variables' responses to wood harvest in the 2000s (Figures S11–S16 in Supporting Information S1) are highly related to the harvest rate pattern. Regions with high harvest rates (Figure S2 in Supporting Information S1) show obvious biogeophysical impact with statistical significance. High loss of canopy coverage (Figure S11 in Supporting Information S1) can be found over multiple regions in temperate zones, which are also the hotspot regions contributing changes of albedo, LAI, roughness length and land surface heat fluxes. Over the tropics, the loss of canopy coverage cannot be clearly seen, but the reduction of LAI (Figure S13 in Supporting Information S1) and roughness length (Figure S14 in Supporting Information S1) at the outskirts of the tropical forest area is more obvious, and are the major cause of further change in land surface heat fluxes (Figures S15 and S16 in Supporting Information S1).

We separate the direct canopy loss caused by wood harvest and the post-harvest canopy recovery in order to isolate their individual contributions from the total effect of wood harvest and overall change in canopy coverage. The impact from accumulated historical wood harvesting shows different patterns of canopy loss and recovery between the tropics and extra-tropics (Figure 10). In the extra-tropics (Figure 10a), harvest rates decrease over time with direct canopy loss rate ranging from 0.11% year⁻¹ to 0.02% year⁻¹. Canopy loss and recovery rates are approximately balanced, as shown by the data points falling close to the 1:1 line. In the 1850s, post-harvest recovery rate of the forest compensated and made up for harvesting losses. However, by the 2000s, recovery rates show little to no compensating increase and even slight decrease in certain cases compared to the 1850s, thus switching to an average net loss of 10% forest coverage. For the tropics (Figure 10b), direct canopy loss and recovery are closer to quasi-balance status than extra-tropics. Overall, the regrowth from both regions slightly exceeds the direct loss in the 1850s, but higher wood harvest loss than regrowth in the 2000s generally leads to a net loss of canopy coverage.

[IMAGE OMITTED. SEE PDF]

Inter-Case Variation of Biogeophysical Response to Wood Harvesting

ELM-FATES simulations show clear differences in biogeophysical responses to harvest due to the harvest approach (Figure 7). All cases driven by the mass-based show larger loss of canopy coverage, LAI and surface roughness length, compared to the corresponding cases using the area-based harvest approach. This generally corresponds with more carbon (Figure 5) and area harvested (Figure 6) under mass-based approach (Figure 5), except for the high harvest scenario. As described above, the mass-based high harvest scenario has much less secondary forest than the area-based one in the 1850s, but the increase of harvest rate is faster, which increases the impact of harvest on structural properties in the 2000s. Higher wood harvest rates in the 2000s cause more significant differences from the 1850s in structural responses for the mass-based cases than for the area-based cases, reflected by more instances of student's t values lower than 0.05 (Tables S1 and S2 in Supporting Information S1). When requiring the harvested C to match the LUH2 historical baseline, that is the har_c case, the responses of most structural variables to harvest in both regions are significant and the response in extra-tropics are double of that in the tropics. The spatial distribution of biogeophysical responses show similar patterns for simulations with the same harvest approach. Generally, cases with mass-based harvest approach show a higher response than area-based harvest.

Although the different wood harvest rate scenarios cause large differences in the spatio-temporal pattern of carbon cycle disturbance rate among different cases, the choice of harvest approach leads to more differences in the biogeophysical responses than the harvest rate scenarios (Figure 7; Tables S3 and S4 in Supporting Information S1). One reason is that the historical harvest rate does not change much between LUH2 historical baseline and historical low cases (Figures 5 and 6), especially when compared to the large variation of harvested C under different harvest approaches. Also, the change of the response with time, from the 1850s–2000s, is smaller than the difference due to the wood harvest approach, although the harvest rate per grid cell is higher in the 2000s than in the 1850s (Figure S6 in Supporting Information S1). We found no clear cross-case discrepancy between the biogeophysical responses in tropical and extratropical regions, but the impact in the extratropical region is slightly stronger among different cases compared to the tropical region (Figure 7).

The magnitudes of canopy structural responses from har_c cases (e.g., −4.8%. −4.5%, −6.4% for canopy coverage, LAI and roughness length in the 2000s extra-tropical region) are more than five times larger than those estimated by the ELM-only model (0%, −0.6%, −1.3% for canopy coverage, LAI and roughness length, Figure 7), a consequence from both a larger disturbed area (Figure S4 in Supporting Information S1) and the new demographic processes. The comparison of ELM-FATES his_a case to ELM-only suggests that even under similar disturbed areas, the inclusion of demographic model processes can substantially affect the canopy structural variables thus creating varied biogeophysical responses. The ELM-only response to wood harvesting is less than 1% in LAI and 1.6% in surface roughness length, a consequence of lacking dynamic canopy geometry and plant height, thus leading to smaller changes in albedo and relative surface energy redistribution. Notably, the calculated changes in energy fluxes (Figure 8) show much less variation between FATES enabled and ELM-only simulations, with all cases giving the response at the same order of magnitude. One reason for this small variation is related to the offline simulation with identical atmospheric forcings. Another reason is that both models share the same energy calculation as performed by ELM, with the energy balance of the FATES cohorts aggregated to the ELM grid cell, which obfuscates the finer-scale differences from FATES.

Discussion

Comparison With Previous Studies

Previous studies examining biogeophysical responses to wood harvesting are either limited to large scale vegetation removal like clear-cuts at the site-level or global studies with reduced complexity land surface models. Results from published global studies (De Hertog et al., 2023; Lawrence et al., 2012) are similar to the result from our ELM-only simulation and show almost negligible responses. ELM-FATES simulations advance our interpretation through applying a more detailed model incorporating plant physiology and demographic processes in place of the sun-shade canopy component in traditional land surface models. Since this work is the first global study of wood harvesting using a demographically structured land model, we focus more on the magnitude and direction of the responses from previous studies instead of a more comprehensive, quantitative comparison that was not feasible due to lack of studies.

The direct wood harvesting impact is the change of forest leaf biomass and canopy coverage. Regrowth of forest after wood harvesting is considered to be one major factor determining the magnitude of impact on canopy coverage. Post-logging recovery is faster compared to other disturbances such as post-fire due to the surviving saplings helping regeneration (Miller et al., 2011; Seedre & Chen, 2010). A study of boreal forest in eastern Canada with long term chronosequences of post-logging samples show around 9% recovery of the pre-harvest biomass after 9 years recovery (Strukelj et al., 2015), partially supporting our finding of persistently less canopy coverage under continuous harvesting stress compared to canopy closure status. Site level studies lack enough long-term and large area extent experiments to study the canopy coverage change under continuous wood harvest, which presents one of the major difficulties to evaluating our results.

Change in albedo depends on the regrowth rate of the canopy. Remote sensing based studies on the correlation between boreal forest age and albedo for different seasons shows decrease in albedo from secondary sapling to mature forest of 0.03 (summer), 0.04 (spring) and 0.3 (winter), respectively (Halim et al., 2019; Kuusinen et al., 2014). This age-dependent decrease in albedo may be related to an increase of leaf and canopy area or the change of canopy structure when plants get mature, which in turn will influence the albedo response of clear-cut. The larger decrease of albedo during winter is related to the land snow coverage, which has a much higher albedo magnitude and variability than vegetation. Albedo changes in tropical forest are affected more by the groundwater level than seasonality and other environmental factors, and the regrowth from post-logging shows very small albedo variation (Ohkubo et al., 2021). Consistent to these site level studies, ELM-FATES calculated albedo increases in all cases with a larger magnitude in extratropical than in tropical zones, but the numbers cannot be evaluated directly since the simulated grid cell covers a very large area (4 × 5°) with smaller harvested fractional area than the site studies.

Changes in roughness length and land surface energy fluxes are important consequences from wood harvesting. Roughness length directly links to canopy coverage and the averaged tree height and is highly simplified in the sun-shade canopy model but has sophisticated parameterizations in FATES. Studies show a large contribution of fractional canopy coverage and roughness length to the surface energy partition (Nakai et al., 2008; Yuan et al., 2021), and surface roughness change has impacted land surface temperature by changing the aerodynamic resistance and surface wind profile (De Hertog et al., 2023). Roughness length under the early successional stage is lower compared to the mature forest and has a close correlation with LAI (Maurer et al., 2013), supporting our results of decreasing LAI with decreasing roughness length. However, the response of energy fluxes to clear-cut is not coherent but site-specific. Substantial increases of LHF and decreases in SHF are found in a temperate broadleaf environment (Williams et al., 2014) after a one-time clearcut. A further study shows that the source of the increase of LHF is from the surface vegetation after clearcut in a boreal peatland forest (Korkiakoski et al., 2019), which receives more direct radiation post-harvesting. Another boreal peatland forest site-level study shows a decrease of both SHF and LHF (Leppä et al., 2020). A global study of transition of vegetation to other ecosystem types through data-driven assessment suggests a decrease of LHF and increase of SHF (Duveiller et al., 2018). Divergent response patterns are also reflected in our work, with most of the simulation cases showing a decrease of LHF compensated by an increase of SHF while some cases show the opposite response. Complicated linkages among plant physiology, demography and environmental envelope can cause divergent responses and further studies on the geographical distribution and mechanism of these responses are in demand.

Uncertainties From Data, Model Settings and Harvest Approaches

This study shows divergent biogeophysical responses to wood harvesting as a result of different harvest approaches and application of driving data that has been reharmonized or not, reflected by large disagreement on harvested carbon amount and area (Figure 5). Currently, ELM-FATES ingests harvest rates as a spatially-resolved data set, that mismatched the model derived forest inventory supply. One solution to this issue is to couple the re-harmonization algorithm with FATES and build a dynamic, spatial downscaling submodule that uses regional or country level harvest data. Compared to the current re-harmonization method (Figure 2), a coupled system could adjust the forest inventory and demographic information, such as age and size, at runtime to avoid mismatches of harvest rates and area (Figure S1 in Supporting Information S1) and to mitigate the discrepancy between the area-based and mass-based harvest approach. However, the re-harmonization algorithm may require lateral sharing of information across grid cells, thus increasing the difficulty for model parallelization.

Besides the model structure, this work shows a significant difference of the calculated harvested C and responses of all biogeophysical variables between different harvest approaches. Previous works (Lawrence et al., 2012) have used the area-based harvest rate to drive land surface models and focus more on studying the land use carbon emission from wood harvest. Our study highlights the possible deviation of estimated harvested carbon can reach up to 50 PgC and the harvested forest area up to 110%, depending on the forest C calculated from the land surface model. To reduce this uncertainty source, a better approach is to only use regional or country level harvested forest carbon and determine the disturbance rate caused by wood harvest through model calculated forest inventory. Meantime, the estimation of forest biomass is crucial for more accurate land use C emission under this approach, since the deviation of forest inventory amplifies the differences in the calculated harvested C.

Without demographic processes, the harvest influence is small, while the simulations that account for demographic processes will amplify the biogeophysical responses. A large portion of the uncertainty can be caused by the bias of FATES C calculation to the true observation values. Compared to evaluation of latitudinal patterns, FATES simulated forest C in temperate and boreal zones are at lower end, thus cases driven by mass-based harvest rate would possibly overestimate harvest responses, especially at early simulation years since the disturbed area per year could be overestimated. Conversely, the responses in the tropical zone may be underestimated due to higher FATES C. Re-harmonization algorithms further modify the harvest rate spatial pattern and tend to exaggerate the area affected by wood harvest when forest C is underestimated, which corresponds to the results of har_c case. The increase of disturbed area leads to a stronger direct response to wood harvest, while the indirect response (e.g., plant-soil interactions, edge effects, etc.) on forest growth depends on whether the increased area is small enough that the site can still recover from the disturbance. Nonetheless, the advantage of mass-based approach is that it corresponds directly with the source data, which is harvested wood C, and the harvested area is estimated from carbon/biomass data, similar to reality. This also suggests that improving the vegetation carbon content in models and observational estimates is a good way to improve estimates of wood harvest impacts on the Earth system.

Biogeophysical responses are also affected by the configuration of FATES of choice that does not consider competition between PFTs. Competition between early and late successional plants with diverse functional traits can lead to different regrowth trajectories, depending the choice of parameters, thus affecting the recovery of canopy coverage, but this additional complexity brings more parameters and is difficult to constrain based on currently available global data sets. In this study, FATES uses fixed biogeography from historical PFT distribution maps to balance the tradeoff between model efficiency (Li et al., 2023) and complexity (Koven et al., 2020) and to isolate the effects of wood harvest from land cover change. Biogeography components in FATES introduce more parameters that need calibration to accurately represent PFT dynamics in diverse ecosystems, which requires additional work on global scale model testing.

Coarse spatial resolution also reduces the capability of the model to detect regional or local biogeophysical responses. In current land models, the coarser resolution holds multiple PFTs in one gridcell but share the same column level soil and energy states, thus when only a fraction of the forest is harvested the responses are obscured by the unaffected fraction and become smaller than if having finer resolution with forest-only cells. As a consequence, one variable reflecting this issue is the displacement height from wood harvesting. A site level experiment (Maurer et al., 2013) shows a reduction of displacement height in response to the clear-cut of forest. In our study, we find the same direction of response but the reduction of displacement height under coarse resolution simulation is tiny and requires further high resolution study to evaluate.

Model Parameterization to Account for Successional Processes and Forestry Management

FATES currently represents plant demography and structural dynamics, and ongoing model improvements may increase the accuracy of wood harvest impacts. For example, enabling nutrient cycling (Knox et al., 2024) in response to logging may affect regrowth. Another issue is that although FATES has multi-layer canopy representation, ELM shares the same calculated stomatal conductance with the version of FATES used here, which may be a major reason for the lack of pronounced differences in surface energy changes after logging between the two models. Within canopies, leaves of varying inclination angle interact with, respond to, and produce different radiation environments that can benefit from 3-D radiative transfer over the current scheme.

Second, land-atmosphere feedbacks are expected to strengthen the wood harvest impacts on energy fluxes. Due to the strong roughness length change, the actual near surface air temperature may change, though the relative magnitude compared to LUC and other land management requires more advanced studies. Offline land simulations, like in this study, have smaller surface temperature and humidity gradients than simulations with coupled land and atmosphere, which can reduce energy fluxes (Laguë et al., 2019). Not including such feedback may also underestimate the response of energy fluxes to wood harvesting.

Finally, more detailed forest management will directly influence the estimated impacts of wood harvest (Lindeskog et al., 2021). Logging varies across the globe due to many factors, and it has changed over time, so estimated biogeophysical impacts due to clearcut may not be representative of divergent management practices across regions and over time (Nolet et al., 2018; Pereira et al., 2002; Rutishauser et al., 2015). Studies suggest that different forest management plans can have a distinct impact on carbon stock (Ameray et al., 2021). Instead of clear-cut, selective logging can partially maintain forest coverage but usually disturbs a larger area than clear-cut and thus tends to increase the overall impact around 3–5 times assuming similar production of round wood (Nolet et al., 2018). Intensive and extensive forest management strategies can yield multiple harvest options (e.g., clear-cut, partial-cuts, selective logging and thinning etc.) and rotations, which can substantially alter forest structure and regrowth rates, both of which influence biogeophysical responses. However, new implementations of more detailed management practices may lead to a stronger impact on the carbon cycle (Lee et al., 2002; Strukelj et al., 2015) than on some of the structural aspects such as canopy coverage, albedo and surface roughness (De Hertog et al., 2023). New silvicultural options for wood harvesting such as variable retention (Xing et al., 2018) may help reduce the biogeophysical impact but require further investigation.

Conclusions

Biogeophysical responses from logging were thought to be a minor factor among all land cover change and land use management practices. However, we show a stronger impact can be found when considering the detailed demographic processes in land surface models. One reason for this is that our demographic model decreases canopy coverage and surface roughness length up to 5% more over harvested area globally when compared to a common sun-shade canopy model.

This study contributes to a global perspective of the importance of wood harvest activities in reshaping the forest structure, albedo, and surface energy balance. The results also suggest that there is a substantial change of land surface properties due to global wood harvest and we recommend ESMs to resolve plant structure, representing secondary forest, and applying accurate mass-based harvest rates in land surface models for better estimating the impacts of land management and LUCs. The addition of plant demographics and canopy structure to land models can further our understanding of human influences on the Earth system as additional land practices are proposed and implemented in response to global change.

Acknowledgments

This research was supported by the Energy Exascale Earth System Modeling (E3SM, ) Project, funded by the Office of Biological and Environmental Research of the US Department of Energy Office of Science. This research also used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231. Lawrence Berkeley National Laboratory (LBNL) is managed by the University of California for the U.S. Department of Energy under contract DE-AC02-05CH11231. We sincerely thank the two anonymous reviewers for their valuable feedback, insightful comments, and time spent reviewing and improving the quality of manuscript.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

Data Availability Statement

All data, code and models that support the findings of this study are openly available. The code version of E3SM and FATES producing results in this manuscript are available from the following github reporsitories: E3SM: , and FATES: . A copy of the E3SM code can also be accessed from Foucar et al. (2025). A copy of the FATES code can also be accessed from FATES development team (2025). The code for re-harmonization of LUH2 historical harvest rate and producing figures can be accessed from the following github repository: . All model outputs can be accessed from S. Shu (2024).

References

Alemohammad, S. H., Fang, B., Konings, A. G., Aires, F., Green, J. K., Kolassa, J., et al. (2017). Water, energy, and carbon with artificial Neural Networks (WECANN): A statistically based estimate of global surface turbulent fluxes and gross primary productivity using solar‐induced fluorescence. Biogeosciences, 14(18), 4101–4124. https://doi.org/10.5194/bg‐14‐4101‐2017

Alibakhshi, S., Naimi, B., Hovi, A., Crowther, T. W., & Rautiainen, M. (2020). Quantitative analysis of the links between forest structure and land surface albedo on a global scale. Remote Sensing of Environment, 246, 111854. https://doi.org/10.1016/j.rse.2020.111854

Ameray, A., Bergeron, Y., Valeria, O., Montoro Girona, M., & Cavard, X. (2021). Forest carbon management: A review of silvicultural practices and management strategies across boreal, temperate and tropical forests. Current Forestry Reports, 7(4), 245–266. https://doi.org/10.1007/s40725‐021‐00151‐w

Besnard, S., Carvalhais, N., Arain, M. A., Black, A., Bruin, S. D., Buchmann, N., et al. (2018). Quantifying the effect of forest age in annual net forest carbon balance. Environmental Research Letters, 13(12), 124018. https://doi.org/10.1088/1748‐9326/aaeaeb

Bonan, G. B. (2008). Forests and climate change: Forcings, feedbacks, and the climate benefits of forests. Science, 320(5882), 1444–1449. https://doi.org/10.1126/science.1155121

Bourque, C. P.‐A., Daugharty, D. A., Dickison, R. B. B., & Arp, P. A. (1995). Changes in Albedo of a Northern Hardwood Forest following clearcutting. Forest Science, 41(2), 268–277. https://doi.org/10.1093/forestscience/41.2.268

Buotte, P. C., Koven, C. D., Xu, C., Shuman, J. K., Goulden, M. L., Levis, S., et al. (2021). Capturing functional strategies and compositional dynamics in vegetation demographic models. Biogeosciences, 18(14), 4473–4490. https://doi.org/10.5194/bg‐18‐4473‐2021

Calle, L., Canadell, J. G., Patra, P., Ciais, P., Ichii, K., Tian, H., et al. (2016). Regional carbon fluxes from land use and land cover change in Asia, 1980–2009. Environmental Research Letters, 11(7), 074011. https://doi.org/10.1088/1748‐9326/11/7/074011

Canadell, J. G., Raupach, M. R., & Houghton, R. A. (2009). Anthropogenic CO2 emissions in Africa. Biogeosciences, 6(3), 463–468. https://doi.org/10.5194/bg‐6‐463‐2009

Collatz, G. J., Ball, J. T., Grivet, C., & Berry, J. A. (1991). Physiological and environmental regulation of stomatal conductance, photosynthesis and transpiration: A model that includes a laminar boundary layer. Agricultural and Forest Meteorology, 54(2), 107–136. https://doi.org/10.1016/0168‐1923(91)90002‐8

Collier, N., Hoffman, F. M., Lawrence, D. M., Keppel‐Aleks, G., Koven, C. D., Riley, W. J., et al. (2018). The international land model benchmarking (ILAMB) System: Design, theory, and implementation. Journal of Advances in Modeling Earth Systems, 10(11), 2731–2754. https://doi.org/10.1029/2018MS001354

Compo, G. P., Whitaker, J. S., Sardeshmukh, P. D., Matsui, N., Allan, R. J., Yin, X., et al. (2011). The twentieth century reanalysis project. Quarterly Journal of the Royal Meteorological Society, 137(654), 1–28. https://doi.org/10.1002/qj.776

Covington, W. W. (1981). Changes in Forest floor organic matter and nutrient content following clear cutting in Northern hardwoods. Ecology, 62(1), 41–48. https://doi.org/10.2307/1936666

Curasi, S. R., Melton, J. R., Humphreys, E. R., Hermosilla, T., & Wulder, M. A. (2024). Implementing a dynamic representation of fire and harvest including subgrid‐scale heterogeneity in the tile‐based land surface model CLASSIC v1.45. Geoscientific Model Development, 17(7), 2683–2704. https://doi.org/10.5194/gmd‐17‐2683‐2024

De Hertog, S. J., Havermann, F., Vanderkelen, I., Guo, S., Luo, F., Manola, I., et al. (2023). The biogeophysical effects of idealized land cover and land management changes in Earth system models. Earth System Dynamics, 14(3), 629–667. https://doi.org/10.5194/esd‐14‐629‐2023

Dickinson, R. E., & Kennedy, P. (1992). Impacts on regional climate of Amazon deforestation. Geophysical Research Letters, 19(19), 1947–1950. https://doi.org/10.1029/92gl01905

Di Vittorio, A. V., Vernon, C. R., & Shu, S. (2020). Moirai version 3: A data processing System to generate recent historical land inputs for global modeling applications at various scales. Journal of Open Research Software, 8(1), 1. https://doi.org/10.5334/jors.266

Duveiller, G., Hooker, J., & Cescatti, A. (2018). The mark of vegetation change on Earth’s surface energy balance. Nature Communications, 9(1), 679. https://doi.org/10.1038/s41467‐017‐02810‐8

FATES Development Team. (2024). The Functionally Assembled Terrestrial Ecosystem Simulator (FATES). https://doi.org/10.5281/zenodo.3825474

FATES Development Team. (2025). The Functionally Assembled Terrestrial Ecosystem Simulator (FATES) (archive/ShuS_JAMES_2025_ELMFATES) [Computer software]. Zenodo. https://doi.org/10.5281/ZENODO.17232785

Fisher, R. A., Koven, C. D., Anderegg, W. R. L., Christoffersen, B. O., Dietze, M. C., Farrior, C. E., et al. (2018). Vegetation demographics in Earth System Models: A review of progress and priorities. Global Change Biology, 24(1), 35–54. https://doi.org/10.1111/gcb.13910

Fisher, R. A., Muszala, S., Verteinstein, M., Lawrence, P., Xu, C., McDowell, N. G., et al. (2015). Taking off the training wheels: The properties of a dynamic vegetation model without climate envelopes, CLM4.5. Geoscientific Model Development, 8(11), 3593–3619. https://doi.org/10.5194/gmd‐8‐3593‐2015

Foucar, J., Edwards, J., Petersen, M., Bertagna, L., Hoffman, M., Jacobsen, D. (2025). sshu88/E3SM: ELMFATES_global_wood_harvest (archive/ShuS_JAMES_2025_ELMFATES) [Computer software]. Zenodo. https://doi.org/10.5281/ZENODO.17232744

Fourrier, A., Pothier, D., & Bouchard, M. (2013). A comparative study of long‐term stand growth in eastern Canadian boreal forest: Fire versus clear‐cut. Forest Ecology and Management, 310, 10–18. https://doi.org/10.1016/j.foreco.2013.08.011

Friedlingstein, P., O’Sullivan, M., Jones, M. W., Andrew, R. M., Gregor, L., Hauck, J., et al. (2022). Global Carbon Budget 2022. Earth System Science Data, 14(11), 4811–4900. https://doi.org/10.5194/essd‐14‐4811‐2022

Golaz, J.‐C., Van Roekel, L. P., Zheng, X., Roberts, A. F., Wolfe, J. D., Lin, W., et al. (2022). The DOE E3SM Model version 2: Overview of the physical model and initial model evaluation. Journal of Advances in Modeling Earth Systems, 14(12), e2022MS003156. https://doi.org/10.1029/2022MS003156

Halim, M. A., Chen, H. Y. H., & Thomas, S. C. (2019). Stand age and species composition effects on surface albedo in a mixedwood boreal forest. Biogeosciences, 16(22), 4357–4375. https://doi.org/10.5194/bg‐16‐4357‐2019

He, T., Liang, S., & Song, D.‐X. (2014). Analysis of global land surface albedo climatology and spatial‐temporal variation during 1981–2010 from multiple satellite products. Journal of Geophysical Research: Atmospheres, 119(17), 10281–10298. https://doi.org/10.1002/2014JD021667

Hobeichi, S., Abramowitz, G., & Evans, J. (2020). Conserving land–atmosphere synthesis suite (CLASS). Journal of Climate, 33(5), 1821–1844. https://doi.org/10.1175/JCLI‐D‐19‐0036.1

Hoffmann, W. A., & Jackson, R. B. (2000). Vegetation–Climate feedbacks in the conversion of tropical savanna to grassland. Journal of Climate, 13(9), 1593–1602. https://doi.org/10.1175/1520‐0442(2000)013<1593:vcfitc>2.0.co;2

Holm, J. A., Knox, R. G., Zhu, Q., Fisher, R. A., Koven, C. D., Nogueira Lima, A. J., et al. (2020). The central Amazon biomass sink under Current and future atmospheric CO2: Predictions from big‐leaf and demographic vegetation models. Journal of Geophysical Research: Biogeosciences, 125(3), e2019JG005500. https://doi.org/10.1029/2019JG005500

Houghton, R. A., & Nassikas, A. A. (2018). Negative emissions from stopping deforestation and forest degradation, globally. Global Change Biology, 24(1), 350–359. https://doi.org/10.1111/gcb.13876

Huang, M., & Asner, G. P. (2010). Long‐term carbon loss and recovery following selective logging in Amazon forests. Global Biogeochemical Cycles, 24(3), Gb3028. https://doi.org/10.1029/2009GB003727

Huang, M., Xu, Y., Longo, M., Keller, M., Knox, R. G., Koven, C. D., & Fisher, R. A. (2020). Assessing impacts of selective logging on water, energy, and carbon budgets and ecosystem dynamics in Amazon forests using the functionally assembled terrestrial ecosystem simulator. Biogeosciences, 17(20), 4999–5023. https://doi.org/10.5194/bg‐17‐4999‐2020

Hurtt, G. C., Chini, L., Sahajpal, R., Frolking, S., Bodirsky, B. L., Calvin, K., et al. (2020). Harmonization of global land use change and management for the period 850–2100 (LUH2) for CMIP6. Geoscientific Model Development, 13(11), 5425–5464. https://doi.org/10.5194/gmd‐13‐5425‐2020

Hurtt, G. C., Chini, L. P., Frolking, S., Betts, R. A., Feddema, J., Fischer, G., et al. (2011). Harmonization of land‐use scenarios for the period 1500–2100: 600 years of global gridded annual land‐use transitions, wood harvest, and resulting secondary lands. Climatic Change, 109(1), 117–161. https://doi.org/10.1007/s10584‐011‐0153‐2

Jung, M., Reichstein, M., Schwalm, C. R., Huntingford, C., Sitch, S., Ahlström, A., et al. (2017). Compensatory water effects link yearly global land CO2 sink changes to temperature. Nature, 541(7638), 516–520. https://doi.org/10.1038/nature20780

Kalliokoski, T., Bäck, J., Boy, M., Kulmala, M., Kuusinen, N., Mäkelä, A., et al. (2020). Mitigation impact of different harvest scenarios of finnish forests that account for Albedo, aerosols, and trade‐offs of carbon sequestration and avoided emissions. Frontiers in Forests and Global Change, 3, 562044. https://doi.org/10.3389/ffgc.2020.562044

Kellomäki, S., Väisänen, H., Kirschbaum, M. U. F., Kirsikka‐Aho, S., & Peltola, H. (2021). Effects of different management options of Norway spruce on radiative forcing through changes in carbon stocks and albedo. Forestry: International Journal of Financial Research, 94(4), 588–597. https://doi.org/10.1093/forestry/cpab010

Kim, H. (2017). Global soil wetness project phase 3 atmospheric boundary conditions (experiment 1) [Dataset]. Data Integration and Analysis System (DIAS). https://doi.org/10.20783/DIAS.501

Knox, R. G., Koven, C. D., Riley, W. J., Walker, A. P., Wright, S. J., Holm, J. A., et al. (2024). Nutrient dynamics in a coupled terrestrial biosphere and land model (ELM‐FATES‐CNP). Journal of Advances in Modeling Earth Systems, 16(3), e2023MS003689. https://doi.org/10.1029/2023MS003689

Korkiakoski, M., Tuovinen, J.‐P., Penttilä, T., Sarkkola, S., Ojanen, P., Minkkinen, K., et al. (2019). Greenhouse gas and energy fluxes in a boreal peatland forest after clear‐cutting. Biogeosciences, 16(19), 3703–3723. https://doi.org/10.5194/bg‐16‐3703‐2019

Koven, C. D., Knox, R. G., Fisher, R. A., Chambers, J. Q., Christoffersen, B. O., Davies, S. J., et al. (2020). Benchmarking and parameter sensitivity of physiological and vegetation dynamics using the Functionally Assembled Terrestrial Ecosystem Simulator (FATES) at Barro Colorado Island, Panama. Biogeosciences, 17(11), 3017–3044. https://doi.org/10.5194/bg‐17‐3017‐2020

Kuusinen, N., Tomppo, E., Shuai, Y., & Berninger, F. (2014). Effects of forest age on albedo in boreal forests estimated from MODIS and Landsat albedo retrievals. Remote Sensing of Environment, 145, 145–153. https://doi.org/10.1016/j.rse.2014.02.005

Laguë, M. M., Bonan, G. B., & Swann, A. L. S. (2019). Separating the impact of individual land surface properties on the terrestrial surface energy budget in both the coupled and uncoupled land–atmosphere System. Journal of Climate, 32(18), 5725–5744. https://doi.org/10.1175/JCLI‐D‐18‐0812.1

Lambert, M. S. A., Tang, H., Aas, K. S., Stordal, F., Fisher, R. A., Bjerke, J. W., et al. (2023). Integration of a frost mortality scheme into the demographic vegetation model FATES. Journal of Advances in Modeling Earth Systems, 15(7), e2022MS003333. https://doi.org/10.1029/2022MS003333

Law, B. E., Hudiburg, T. W., Berner, L. T., Kent, J. J., Buotte, P. C., & Harmon, M. E. (2018). Land use strategies to mitigate climate change in carbon dense temperate forests. Proceedings of the National Academy of Sciences, 115(14), 3663–3668. https://doi.org/10.1073/pnas.1720064115

Lawrence, P. J., & Chase, T. N. (2007). Representing a new MODIS consistent land surface in the Community Land Model (CLM 3.0). Journal of Geophysical Research, 112(G1), G01023. https://doi.org/10.1029/2006JG000168

Lawrence, P. J., Feddema, J. J., Bonan, G. B., Meehl, G. A., O’Neill, B. C., Oleson, K. W., et al. (2012). Simulating the biogeochemical and biogeophysical impacts of transient land cover change and wood harvest in the Community Climate System Model (CCSM4) from 1850 to 2100. Journal of Climate, 25(9), 3071–3095. https://doi.org/10.1175/JCLI‐D‐11‐00256.1

Lee, J., Morrison, I. K., Leblanc, J.‐D., Dumas, M. T., & Cameron, D. A. (2002). Carbon sequestration in trees and regrowth vegetation as affected by clearcut and partial cut harvesting in a second‐growth boreal mixedwood. Forest Ecology and Management, 169(1), 83–101. https://doi.org/10.1016/S0378‐1127(02)00300‐6

Leppä, K., Korkiakoski, M., Nieminen, M., Laiho, R., Hotanen, J.‐P., Kieloaho, A.‐J., et al. (2020). Vegetation controls of water and energy balance of a drained peatland forest: Responses to alternative harvesting practices. Agricultural and Forest Meteorology, 295, 108198. https://doi.org/10.1016/j.agrformet.2020.108198

Li, L., Fang, Y., Zheng, Z., Shi, M., Longo, M., Koven, C. D., et al. (2023). A machine learning approach targeting parameter estimation for plant functional type coexistence modeling using ELM‐FATES (v2.0). Geoscientific Model Development, 16(14), 4017–4040. https://doi.org/10.5194/gmd‐16‐4017‐2023

Liao, W., Liu, X., Burakowski, E., Wang, D., Wang, L., & Li, D. (2020). Sensitivities and responses of land surface temperature to deforestation‐induced biophysical changes in two global Earth System models. Journal of Climate, 33(19), 8381–8399. https://doi.org/10.1175/JCLI‐D‐19‐0725.1

Lindeskog, M., Smith, B., Lagergren, F., Sycheva, E., Ficko, A., Pretzsch, H., & Rammig, A. (2021). Accounting for forest management in the estimation of forest carbon balance using the dynamic vegetation model LPJ‐GUESS (v4.0, r9710): Implementation and evaluation of simulations for Europe. Geoscientific Model Development, 14(10), 6071–6112. https://doi.org/10.5194/gmd‐14‐6071‐2021

Littleton, E. W., Harper, A. B., Vaughan, N. E., Oliver, R. J., Duran‐Rojas, M. C., & Lenton, T. M. (2020). JULES‐BE: Representation of bioenergy crops and harvesting in the Joint UK Land Environment Simulator vn5.1. Geoscientific Model Development, 13(3), 1123–1136. https://doi.org/10.5194/gmd‐13‐1123‐2020

Liu, S., Bond‐Lamberty, B., Hicke, J. A., Vargas, R., Zhao, S., Chen, J., et al. (2011). Simulating the impacts of disturbances on forest carbon cycling in North America: Processes, data, models, and challenges. Journal of Geophysical Research, 116(G4), G00K08. https://doi.org/10.1029/2010JG001585

Luyssaert, S., Jammet, M., Stoy, P. C., Estel, S., Pongratz, J., Ceschia, E., et al. (2014). Land management and land‐cover change have impacts of similar magnitude on surface temperature. Nature Climate Change, 4(5), 389–393. https://doi.org/10.1038/nclimate2196

Ma, L., Hurtt, G., Ott, L., Sahajpal, R., Fisk, J., Lamb, R., et al. (2022). Global evaluation of the Ecosystem Demography model (ED v3.0). Geoscientific Model Development, 15(5), 1971–1994. https://doi.org/10.5194/gmd‐15‐1971‐2022

Mamkin, V., Kurbatova, J., Avilov, V., Mukhartova, Y., Krupenko, A., Ivanov, D., et al. (2016). Changes in net ecosystem exchange of CO2, latent and sensible heat fluxes in a recently clear‐cut spruce forest in western Russia: Results from an experimental and modeling analysis. Environmental Research Letters, 11(12), 125012. https://doi.org/10.1088/1748‐9326/aa5189

Martens, B., Miralles, D. G., Lievens, H., van der Schalie, R., de Jeu, R. A. M., Fernández‐Prieto, D., et al. (2017). GLEAM v3: Satellite‐based land evaporation and root‐zone soil moisture. Geoscientific Model Development, 10(5), 1903–1925. https://doi.org/10.5194/gmd‐10‐1903‐2017

Maurer, K. D., Hardiman, B. S., Vogel, C. S., & Bohrer, G. (2013). Canopy‐structure effects on surface roughness parameters: Observations in a Great Lakes mixed‐deciduous forest. Agricultural and Forest Meteorology, 177, 24–34. https://doi.org/10.1016/j.agrformet.2013.04.002

Meiyappan, P., Jain, A. K., & House, J. I. (2015). Increased influence of nitrogen limitation on CO2 emissions from future land use and land use change. Global Biogeochemical Cycles, 29(9), 1524–1548. https://doi.org/10.1002/2015GB005086

Miller, S. D., Goulden, M. L., Hutyra, L. R., Keller, M., Saleska, S. R., Wofsy, S. C., et al. (2011). Reduced impact logging minimally alters tropical rainforest carbon and energy exchange. Proceedings of the National Academy of Sciences, 108(48), 19431–19435. https://doi.org/10.1073/pnas.1105068108

Moorcroft, P. R., Hurtt, G. C., & Pacala, S. W. (2001). A method for scaling vegetation dynamics: The ecosystem demography model (ed). Ecological Monographs, 71(4), 557–586. https://doi.org/10.1890/0012‐9615(2001)071[0557:AMFSVD]2.0.CO;2

Noblet‐Ducoudré, N. D., Boisier, J.‐P., Pitman, A., Bonan, G. B., Brovkin, V., Cruz, F., et al. (2012). Determining robust impacts of land‐use‐induced land cover changes on surface climate over North America and Eurasia: Results from the first set of LUCID experiments. Journal of Climate, 25(9), 3261–3281. https://doi.org/10.1175/JCLI‐D‐11‐00338.1

Nabel, J. E. M. S., Naudts, K., & Pongratz, J. (2020). Accounting for forest age in the tile‐based dynamic global vegetation model JSBACH4 (4.20p7; git feature/forests)—A land surface model for the ICON‐ESM. Geoscientific Model Development, 13(1), 185–200. https://doi.org/10.5194/gmd‐13‐185‐2020

Nakai, T., Sumida, A., Daikoku, K., Matsumoto, K., van der Molen, M. K., Kodama, Y., et al. (2008). Parameterisation of aerodynamic roughness over boreal, cool‐ and warm‐temperate forests. Agricultural and Forest Meteorology, 148(12), 1916–1925. https://doi.org/10.1016/j.agrformet.2008.03.009

Needham, J. F., Dey, S., Koven, C. D., Fisher, R. A., Knox, R. G., Lamour, J., et al. (2025). Vertical canopy gradients of respiration drive plant carbon budgets and leaf area index. New Phytologist, 246(1), 144–157. https://doi.org/10.1111/nph.20423

Nolet, P., Kneeshaw, D., Messier, C., & Béland, M. (2018). Comparing the effects of even‐ and uneven‐aged silviculture on ecological diversity and processes: A review. Ecology and Evolution, 8(2), 1217–1226. https://doi.org/10.1002/ece3.3737

Norman, J. M. (1979). Modeling the complete crop canopy. In Modification of the aerial environment of crops (pp. 249–280).

Ohkubo, S., Hirano, T., & Kusin, K. (2021). Influence of disturbances and environmental changes on albedo in tropical peat ecosystems. Agricultural and Forest Meteorology, 301–302, 108348. https://doi.org/10.1016/j.agrformet.2021.108348

Otto, J., Berveiller, D., Bréon, F.‐M., Delpierre, N., Geppert, G., Granier, A., et al. (2014). Forest summer albedo is sensitive to species and thinning: How should we account for this in Earth system models? Biogeosciences, 11(8), 2411–2427. https://doi.org/10.5194/bg‐11‐2411‐2014

Pearson, T. R. H., Brown, S., & Casarim, F. M. (2014). Carbon emissions from tropical forest degradation caused by logging. Environmental Research Letters, 9(3), 034017. https://doi.org/10.1088/1748‐9326/9/3/034017

Peichl, M., Martínez‐García, E., Fransson, J. E. S., Wallerman, J., Laudon, H., Lundmark, T., & Nilsson, M. B. (2023). On the uncertainty in estimates of the carbon balance recovery time after forest clear‐cutting. Global Change Biology, 29(15), e1–e3. https://doi.org/10.1111/gcb.16772

Peng, L., Searchinger, T. D., Zionts, J., & Waite, R. (2023). The carbon costs of global wood harvests. Nature, 620(7972), 110–115. https://doi.org/10.1038/s41586‐023‐06187‐1

Pereira, R., Jr., Zweede, J., Asner, G. P., & Keller, M. (2002). Forest canopy damage and recovery in reduced‐impact and conventional selective logging in eastern Para, Brazil. Forest Ecology and Management, 168(1–3), 77–89. https://doi.org/10.1016/S0378‐1127(01)00732‐0

Pongratz, J., Schwingshackl, C., Bultan, S., Obermeier, W., Havermann, F., & Guo, S. (2021). Land use effects on climate: Current State, recent progress, and emerging topics. Current Climate Change Reports, 7(4), 99–120. https://doi.org/10.1007/s40641‐021‐00178‐y

Purves, D. W., Lichstein, J. W., Strigul, N., & Pacala, S. W. (2008). Predicting and understanding forest dynamics using a simple tractable model. Proceedings of the National Academy of Sciences, 105(44), 17018–17022. https://doi.org/10.1073/pnas.0807754105

Rutishauser, E., Hérault, B., Baraloto, C., Blanc, L., Descroix, L., Sotta, E. D., et al. (2015). Rapid tree carbon stock recovery in managed Amazonian forests. Current Biology, 25(18), R787–R788. https://doi.org/10.1016/j.cub.2015.07.034

Saatchi, S. S., Harris, N. L., Brown, S., Lefsky, M., Mitchard, E. T. A., Salas, W., et al. (2011). Benchmark map of forest carbon stocks in tropical regions across three continents. Proceedings of the National Academy of Sciences, 108(24), 9899–9904. https://doi.org/10.1073/pnas.1019576108

Santoro, M., & Cartus, O. (2021). ESA Biomass Climate Change Initiative (Biomass_cci): Global datasets of forest above‐ground biomass for the years 2010, 2017 and 2018, v3. NERC EDS Centre for Environmental Data Analysis. https://doi.org/10.5285/5f331c418e9f4935b8eb1b836f8a91b8

Seedre, M., & Chen, H. Y. H. (2010). Carbon dynamics of aboveground live vegetation of boreal mixed woods after wildfire and clear‐cutting. Canadian Journal of Forest Research, 40(9), 1862–1869. https://doi.org/10.1139/X10‐120

Shevliakova, E., Stouffer, R. J., Malyshev, S., Krasting, J. P., Hurtt, G. C., & Pacala, S. W. (2013). Historical warming reduced due to enhanced land carbon uptake. Proceedings of the National Academy of Sciences, 110(42), 16730–16735. https://doi.org/10.1073/pnas.1314047110

Shorohova, E., Sinkevich, S., Kryshen, A., & Vanha‐Majamaa, I. (2019). Variable retention forestry in European boreal forests in Russia. Ecological Processes, 8(1), 34. https://doi.org/10.1186/s13717‐019‐0183‐7

Shu, S. (2024). ELM‐FATES global wood harvest simulations. Release (version 1.1.0) [Dataset]. Zenodo. https://zenodo.org/records/13853938

Stocker, B. D., Feissli, F., Strassmann, K. M., Spahni, R., & Joos, F. (2014). Past and future carbon fluxes from land use change, shifting cultivation and wood harvest. Tellus B: Chemical and Physical Meteorology, 66(1), 23188. https://doi.org/10.3402/tellusb.v66.23188

Strukelj, M., Brais, S., & Paré, D. (2015). Nine‐year changes in carbon dynamics following different intensities of harvesting in boreal aspen stands. European Journal of Forest Research, 134(5), 737–754. https://doi.org/10.1007/s10342‐015‐0880‐4

Su, Y., Liu, L., Wu, J., Chen, X., Shang, J., Ciais, P., et al. (2019). Quantifying the biophysical effects of forests on local air temperature using a novel three‐layered land surface energy balance model. Environment International, 132, 105080. https://doi.org/10.1016/j.envint.2019.105080

Tang, J., Luyssaert, S., Richardson, A. D., Kutsch, W., & Janssens, I. A. (2014). Steeper declines in forest photosynthesis than respiration explain age‐driven decreases in forest growth. Proceedings of the National Academy of Sciences, 111(24), 8856–8860. https://doi.org/10.1073/pnas.1320761111

Thurner, M., Beer, C., Santoro, M., Carvalhais, N., Wutzler, T., Schepaschenko, D., et al. (2014). Carbon stock and density of boreal and temperate forests. Global Ecology and Biogeography, 23(3), 297–310. https://doi.org/10.1111/geb.12125

Wagner, S., Fischer, H., & Huth, F. (2011). Canopy effects on vegetation caused by harvesting and regeneration treatments. European Journal of Forest Research, 130(1), 17–40. https://doi.org/10.1007/s10342‐010‐0378‐z

Williams, C. A., Vanderhoof, M. K., Khomik, M., & Ghimire, B. (2014). Post‐clearcut dynamics of carbon, water and energy exchanges in a midlatitude temperate, deciduous broadleaf forest environment. Global Change Biology, 20(3), 992–1007. https://doi.org/10.1111/gcb.12388

Winckler, J., Reick, C. H., & Pongratz, J. (2017). Robust identification of local biogeophysical effects of land‐cover change in a global climate model. Journal of Climate, 30(3), 1159–1176. https://doi.org/10.1175/JCLI‐D‐16‐0067.1

Xing, D., Nielsen, S. E., Macdonald, S. E., Spence, J. R., & He, F. (2018). Survival and growth of residual trees in a variable retention harvest experiment in a boreal mixed wood forest. Forest Ecology and Management, 411, 187–194. https://doi.org/10.1016/j.foreco.2018.01.026

Xu, L., Saatchi, S. S., Yang, Y., Yu, Y., Pongratz, J., Bloom, A. A., et al. (2021). Changes in global terrestrial live biomass over the 21st century. Science Advances, 7(27), eabe9829. https://doi.org/10.1126/sciadv.abe9829

Yuan, X., Hamdi, R., Ochege, F. U., De Maeyer, P., Kurban, A., & Chen, X. (2021). Assessment of surface roughness and fractional vegetation coverage in the CoLM for modeling regional land surface temperature. Agricultural and Forest Meteorology, 303, 108390. https://doi.org/10.1016/j.agrformet.2021.108390

Yue, C., Ciais, P., Luyssaert, S., Li, W., McGrath, M. J., Chang, J., & Peng, S. (2018). Representing anthropogenic gross land use change, wood harvest, and forest age dynamics in a global vegetation model ORCHIDEE‐MICT v8.4.2. Geoscientific Model Development, 11(1), 409–428. https://doi.org/10.5194/gmd‐11‐409‐2018

Zeng, N., King, A. W., Zaitchik, B., Wullschleger, S. D., Gregg, J., Wang, S., & Kirk‐Davidoff, D. (2013). Carbon sequestration via wood harvest and storage: An assessment of its harvest potential. Climatic Change, 118(2), 245–257. https://doi.org/10.1007/s10584‐012‐0624‐0