The importance of forest structure for carbon fluxes of the Amazon rainforest

The study region covers 7.8 million km² of forests in South America that are categorized as rainforest or moist deciduous rainforest (fewer than 5 dry months in which precipitation (mm) ≤2 times mean temperature (°C) (FAO 2001)), have an annual mean temperature above 18 °C, are located at an elevation below 1000 m and have an AGB >20 t ha⁻¹ (Rödig et al 2017a).

Our analyses are based on a regionalized Amazon-wide version (Rödig et al 2017a) of the forest gap model FORMIND, which has already been applied at various locations in the tropics (Fischer et al 2016). The forest gap model is driven by mean photosynthetic photon flux density (PPFD), mean precipitation (both 0.5°, mean over 2003–2012, Weedon et al 2014), and clay fraction of soil (8 km, Wieder et al 2014). It simulates forest dynamics at the individual tree level considering the following main processes at a yearly time step: tree growth, competition, establishment, and mortality. Growth of an individual tree depends on its location within the forest community, where trees compete for light and space. A gain in tree biomass results from the difference between photosynthesis and respiration losses. A seedling can establish if light intensity on the forest floor is sufficient.

Zoom In Zoom Out Reset image size

In our model, tree mortality is an important driver for forest dynamics. Mortality increases when tree crowns are limited in space (crowding). Falling of large trees can damage surrounding trees (gap building) which causes conditional mortality events. In addition, every tree underlies a basic mortality rate which is determined stochastically. The characteristic of the Amazon-wide forest model is that this basic mortality rate varies spatially depending on temporal mean precipitation and the fraction of clay content in soil (Rödig et al 2017a). That means that in the model, mortality rates vary in space, but not in time. Hence, the forest model accounts for regional forest dynamics under mean climate conditions but does not consider inter-annual climate anomalies. Dead biomass is transferred to a dead wood carbon pool from which carbon is constantly transferred to a soil carbon pool (decomposition) or respired to the atmosphere.

The exchange of carbon between the forest and the atmosphere (net ecosystem productivity NEP [tC ha⁻¹ a⁻¹]), which is sometimes also referred to as net biome productivity (e.g. Jung et al 2011), is described as follows (Paulick et al 2017, Sato et al 2007):

$$\begin{equation} \begin{array}{l@{}l@{}l@{}} {\rm{NEP}} &=& \\ & & \mathop \sum \limits_{{\rm{trees}}} \left( {{\rm{GP}}{{\rm{P}}_{{\rm{tree}}}} - {{{R}}_{{\rm{tree}}}}} \right) + {{{t}}_{{\rm{DA}}}}{{{S}}_{{\rm{dead}}}} + {{{t}}_{{\rm{SA}}}}{{{S}}_{{\rm{slow}}}}+ {{{t}}_{{\rm{FA}}}}{{{S}}_{{\rm{fast}}}} \end{array} \end{equation}$$

The sum of gross primary productivity (GPP_tree [tC ha⁻¹ a⁻¹]) minus autotrophic respiration (R_tree [tC ha⁻¹ a⁻¹]) over all tress equals the woody above-ground NPP (wANPP) of the forest site ([tC ha⁻¹ a⁻¹]). Autotrophic respiration R_tree is calculated as the sum of maintenance and growth respiration which also includes root respiration (figure S1 available at stacks.iop.org/ERL/13/054013/mmedia). Its annual rate is calculated in order to fit observed above-ground biomass growth of a tree. We assume that wANPP is a constant fraction of NPP (NPP = 2.72wANPP, derived for mature forests from Anderson-Teixeira et al 2016). S_dead is the dead wood pool, S_slow the slow decomposing soil carbon pool and S_fast the fast decomposing soil carbon pool (all [tC ha⁻¹]) with its respiration rates to the atmosphere (t_DA S_dead: to atmosphere, t_SA: S_slow to atmosphere, t_FA: S_fast to atmosphere, table S1).

The advantage of simulating each tree individually is that changes in forest structure are captured throughout all different successional states with different species compositions: from bare ground to climax stage including natural tree death. Differences in species composition are represented by three plant functional types: early successional, mid successional and late successional tree types that differ mainly in productivity, light needed at establishment, and mortality rates (table S1). Forest dynamics differ spatially as spatial variable mortality rates cause changes in competition between PFTs. This study focuses on the impact of forest structure on carbon fluxes under mean climatic conditions. This means that tree growth is limited to light and space, but trees grow (in the mean) under mean annual temperature and water conditions.

Technically, the individual-based model could simulate tree growth for every tree in the Amazon rainforest. However, the computational effort can be reduced for areas with similar environmental conditions at 1 km² resolution. That means that for every area with similar environmental conditions (in total 1280 areas with similar mean precipitation, mean photosynthetic photon flux density (PPFD), and clay content, Rödig et al 2017a), only one model realization on 1 km² (100 ha) is performed representatively (figure S2). Our simulation represents 7.8 million km² of forest with 410 billion individual trees (>10 cm). Forest succession is simulated over 1000 years from bare ground to climax state (figure S3) on a high-performance Unix cluster.

The successional state of each forest site within the Amazon is identified via a canopy height map (Simard et al 2011) as in Rödig et al (2017a). The canopy height map is based on space-borne LIDAR measurements taken aboard ICESat in May and June 2005. Consequently, the identified successional state of the forests represents the state in 2005. For each location in the Amazon rainforest (1 km²), we selected the simulated successional states to which the simulated canopy height equals the height given by the canopy height map. A successional state of the forest is here described as the basal area fraction of late successional trees (figure 1). We could then identify the forest's current carbon stock and its associated carbon fluxes (figure S2).

Zoom In Zoom Out Reset image size

The canopy height map comes with a root mean square error of 6.1 m at 1° resolution (Simard et al 2011). We performed an error analysis on how this uncertainty influences simulated GPP, wANPP, and NEP (see supporting information S1). This also includes those uncertainties that result from simulated canopy height and their different identified states. In addition, we tested how fluxes change when keeping the mortality rate constant throughout the Amazon (table S3).

We compared our simulation results against observed NPP and GPP values from ten inventory sites in the lowland Amazon rainforest (Malhi et al 2015), measured woody above-ground net primary productivity (wANPP) from 193 sites (Brienen et al 2015a, 2015b), and GPP and NEP from two eddy covariance sites (<1 km² footprint, annual means between 2004–2014, GF-Guy, Bonal et al 2008, annual means between 2000–2004, BR-Sa3, Goulden et al 2004). Please, note that the eddy covariance method observes the entire net ecosystem exchange (NEE) and we assumed here –NEE = NEP (Chapin et al 2006). Global mean GPP and NPP estimates from MODIS at 1 km² resolution for the years 2000–2010 (Zhao and Running 2010, Running et al 2004) and GPP estimates from FLUXCOM at 0.5° resolution for the years 2000–2010 (mean over all realizations of different machine learning methods, Jung et al 2017, Tramontana et al 2016) were remapped (nearest neighbor, CDO 2015) to the resolution of our map to be compared against our simulation results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

By linking vegetation modelling with canopy height map information from remote sensing (Simard et al 2011), we were able to derive the successional state of forests in the Amazon rainforest (figure 1), here as basal area fraction of late successional trees. AGB increases throughout the successional state of a forest as expected (figure 2, figure S5 for additional attributes) while GPP peaks at early-mid successional state and is significantly higher than in mid-late and late successional state. wANPP decreases with successional state of the forest. NEP is significantly higher in early successional than in late successional state. NEP is around 0 in mid-late to late successional state. All fluxes in all states show significant differences (except NEP in mid-late and late successional state).

Across the Amazon basin, we obtain a mean GPP of forests of 25.1 tC ha⁻¹ a⁻¹ (figure 3(a), table 1). GPP values are higher along the Amazon river and in the southern Amazon rainforest (figure 3(a)). The pattern of woody aboveground productivity (wANPP, figure 3(b)) resembles the one of GPP with higher values along the rivers and in the southern Amazon. Mean wANPP is 4.2 tC ha⁻¹ a⁻¹. In some parts of the Amazon (e.g. the Guiana Shield), GPP values are higher than the mean while wANPP values are lower than the mean. NEP varies around zero (figure 3(c)) with a mean NEP value of 0.7 tC ha⁻¹ a⁻¹. NEP values of up to 5 tC ha⁻¹ a⁻¹ are reached only along the south-east and north-west.

The variation of GPP, wANPP, and NEP values at small scales (0.16 ha resolution) is higher than the variation at 1 km² resolution (frequency distributions in figure 3). For example, the analysis of NEP at 0.16 ha resolution shows that forests can release up to 20 tC ha⁻¹ a⁻¹ to the atmosphere (figure 2(d)), while its maximum release at 1 km² is 0.23 tC ha⁻¹ a⁻¹. Within an area of 7.8 Mio km², the Amazon rainforest takes up 0.56 GtC per year (emissions due to fire and outtake of wood are not included).

We compare our obtained values at 1 km² resolution (figure 4) with estimates derived from MODIS (1 km² resolution, Running et al 2004, Zhao and Running 2010), FLUXCOM (0.5° resolution, Jung et al 2017, Tramontana et al 2016), inventory data (ca. 1 ha plots, Brienen et al 2015a, 2015b, Malhi et al 2015), and eddy-covariance measurements (<1 km² footprint, FLUXNET stations GF-Guy and BR-Sa3).

Gross primary production: The overall frequency distribution of simulated GPP resembles the one derived from MODIS with a mean of 25 tC ha⁻¹ a⁻¹ (figure 4(a)). We observe slightly higher GPP values with MODIS than with our approach for forests in late succession while it seems to be inverse for forests in early succession (figure 5). FLUXCOM GPP is higher with a mean of 28.8 tC ha⁻¹ a⁻¹. GPP derived from inventory data (24–42 tC ha⁻¹ a⁻¹, Malhi et al 2015) and from eddy flux measurements at GF-Guy (32–40 tC ha⁻¹ a⁻¹, Bonal et al 2008) and BR-Sa3 eddy flux station (29–42 tC ha⁻¹ a⁻¹, Goulden et al 2004) fall into the upper ranges of our simulated GPP. A one-to-one comparison shows that eddy flux measurements are in rather good agreement with simulation results at the nearest location (figure S8).

Net primary production: MODIS NPP shows a similar pattern as the related GPP distribution and has a clearly defined peak at 9–10 tC ha⁻¹ a⁻¹ with a mean of 10 tC ha⁻¹ a⁻¹ (figure 4(b)). NPP simulated with FORMIND, on the other hand, is broadly distributed between of 8–16 tC ha⁻¹ a⁻¹. In contrast, if the entire region was assumed to be in climax state, there is a single well defined peak in the NPP distribution (figure S9). NPP values derived from forest inventory range from 9–16 tC ha⁻¹ a⁻¹ (Malhi et al 2015). Estimated wANPP from field inventories (Brienen et al 2015b) range from 1.1–4.7 tC ha⁻¹ a⁻¹, while FORMIND reaches high wANPP values of up to 6 tC ha⁻¹ a⁻¹ (figure S10). Note that the field inventories are limited to measurements (193 sites) mainly taken in old-grown forests.

Net ecosystem productivity: Simulated NEP values under mean climate conditions fall into the range of recordings at the GF-Guy eddy flux station (1.57 tC ha⁻¹ a⁻¹) and the BR-Sa3 station (1.65 tC ha⁻¹ a⁻¹).

Figure 6 shows the obtained relations between carbon fluxes and stocks according to the successional state of a forest plot. Forests in early successional state reach AGB values of 100–150 t ha⁻¹. Forests in early-mid and mid-late successional state reach values up to 300 t ha⁻¹. Forests in late successional state have AGB values greater than 250 t ha⁻¹.

wANPP and AGB stock show a bell-shaped relation (figure 6(a)). GPP values and AGB values form a triangle (figure 6(b)). Forests in early successional states reach higher GPP values than forests in late successional state (with the same AGB). Comparing both figures (a and b) in late successional states, it is evident that wANPP decreases with increasing AGB while GPP increases with AGB.

The relation between GPP and wANPP displays the woody above-ground carbon use efficiency which varies for different successional states (figure 6(c)). Low GPP values are reached for all successional states with a broad range of wANPP values. NEP is always positive for early successional forests, but NEP values of late successional forests show a large variation (figure 6(d)).

GPP values in our map come, on average, with a standard error of 4.15 tC ha⁻¹ a⁻¹ caused by the uncertainties in the canopy height map (figure S6(b), table S2). Errors in the GPP map are slightly higher for the Brazilian Shield (including the 'Arc of Deforestation') than for the other regions. wANPP has a mean standard error of 0.90 tC ha⁻¹ a⁻¹. Errors in wANPP are highest in East Central Amazon and West Amazon. NEP has an average standard error of 0.65 tC ha⁻¹ a⁻¹ and is highest for the Brazilian Shield where the error is 23% higher than the overall mean. All fluxes show lower uncertainties for the Guiana Shield.

When keeping mortality rates constant throughout the Amazon, the standard deviations across the Amazon stay more or less constant (table S3).

Our analysis shows how the successional states of a forest (figure 1) influence carbon stocks and fluxes. Noticeably, wANPP, GPP, and NEP are highest for forests in early to mid-successional states (figures 2(c) and (d)). Such forests of high productivity can be found (figure 3), for example, along the 'Arc of deforestation' in the south-east (Nogueira et al 2008) where forests are degraded by human activities.

The successional state of a forest is not only influenced by large-scale disturbances, but also by individual tree mortality. This means that even undisturbed, large forests in mature state consist of different successional states at the small scale (e.g. 400 m², figure S2) and can hence show fluctuations in its carbon dynamics. Chambers et al (2013) found, for example, that biomass of a mature forest is stable over time only for sample plots greater than 10 ha. For that reason, derived flux values at small scales show larger variability than at coarser resolution (frequency distributions at 400 m² = 0.16 ha vs. 1 km² resolution in figure 3, and figure 6 vs. Figure S4). This may be related to the finding of Espírito-Santo et al (2014) that individual tree mortality over the entire Amazon cause greater biomass losses than large-scale disturbances. Note here that our study includes the consequences of disturbances (regrowth from earlier successional state), but does not account for emissions due to fire or outtake of wood (logging).The relations derived here between simulated stocks and fluxes (figure 6) show that the amount of standing biomass (AGB) is not the only driver for productivity. Our derived relation between wANPP and AGB (figure 6(a)) is in agreement with a global field study which reports that highly productive forests may limit AGB values due to a dominance of species with low wood density (Keeling and Phillips 2007). Our relation therefore differs from simulation results of 4 DGVMs for the Amazon region presented in Johnson et al (2016). In their study, three DGVMs simulate increasing NPP with increasing AGB. Only one DGVM shows a slight reduction of NPP at high AGB values. They conclude that stem mortality rates need to be included in order to capture the relation between woody productivity and biomass as observed in the field (Johnson et al 2016). In our approach, we explicitly simulate stem mortality which is driven by a spatial variation of precipitation and clay content. As a consequence, spatial differences in dynamics and mean wood density are considered (Rödig et al 2017a). This is particularly important for the Amazon where carbon stocks and dynamics are driven by a spatial variation of mortality (Malhi et al 2015, Galbraith et al 2013, Quesada et al 2017, Baker et al 2004, Phillips et al 2004).

The approach brings along structural uncertainties such as the simplified assumption that NPP is a constant fraction of wANPP, parameter uncertainties, and limitations due to the resolution of input data, here climatological data and the canopy height map (Simard et al 2011). The canopy height map is based on discrete LIDAR shots and provides estimates for canopy height at a resolution of 1 km². Missing information between shots is a source for uncertainties in the map (Simard et al 2011). The effect these uncertainties on carbon fluxes is discussed below.

We here make simplifications regarding the handling of disturbance histories of the forests (i.e. fire event, logging). Currently, we cannot reconstruct the type of disturbance event from the canopy height map. Consequently, a disturbed forest site is based on the assumption that it either regenerates after deforestation (simulation from bare ground with initialized dead wood and soil pools, table S1), or is influenced by natural tree fall and its resulting mortality (gap-dynamics).

The study presented here focuses on the spatial variability in carbon dynamics resulting from forest structural and environmental differences under current climatic conditions. Consequently, results capture all fluxes for only one point in time (i.e. an average over the period 2003–2012 with structural information of the canopy height map in 2005, Rödig et al 2017a). The study design does not yet consider the effects of inter- or intra-annual variations of temperature or atmospheric CO₂ on forest structure (Sakschewski et al 2015, Haverd et al 2013). It is anticipated to integrate such effects (as it has already been tested for temperate forests in Bohn and Huth 2017 and Rödig et al., 2017) into future work in order to also have the possibility effects of droughts.

Testing the representativeness of our GPP map relies on the comparison with other maps derived from remote sensing and simulation (MODIS, Running et al 2004, Zhao and Running 2010) and upscaled eddy flux measurements (FLUXCOM, Jung et al 2017, Tramontana et al 2016). Direct carbon flux measurements in the Amazon rainforest are rare (10 inventory plots (Malhi et al 2015) and two eddy flux measurements) and only partly suitable for large-scale estimates.

GPP: The frequency distribution of the MODIS GPP (Running et al 2004, Zhao and Running 2010) resembles our simulated GPP distribution (figure 4(a)) although both maps are derived from different remote sensing techniques (NDVI vs. LIDAR) which lends confidence to our simulations. MODIS and FLUXCOM show particularly low GPP values for the Brazilian Shield for which FORMIND produces higher GPP values than for the rest of the Amazon (table 1, figure 5). These differences probably arise from the fact that our approach uses a representation of forest succession which leads to higher GPP values in early-mid successional states (figure 2(b)). This may also be the reason for differences in GPP values towards the Andes (figure 5, figure S7). In addition, discrepancies may arise from the fact that our approach includes spatial variable mortality rates which lead to differences in forest gap-dynamics across the Amazon. Our simulations at 1 ha resolution reach GPP values comparable to field observations (figure S4), whereas simulated values at 1 km² are lower.

NPP: Our approach (figure 4(b)) displays more forests with high NPP values (>12 tC ha⁻¹ a⁻¹) than the MODIS product. Please note that the NPP pattern of the MODIS product correlates with its GPP pattern (table 1) whereas our approach indicates forest structure as an important additional factor (GPP vs. NPP, carbon use efficiency is not constant, figure 6(c), figure S5(a)). We identify two potential explanations for such discrepancy. First, the representation of different PFTs, which allows for simulating forest succession, leads to higher NPP values (figure 2(c)), especially in highly disturbed regions like the Brazilian shield (table 1). The successional state identified with the canopy height map seems to have more influence on NPP than internal gap-dynamics due to variable mortality rates (table S3). Second, the MODIS product, which is often used for global estimates (Jung et al 2017, Bloom et al 2016, Zhao and Running 2010), is limited by the fact that it is derived from Normalized Difference Vegetation Index (NDVI) values. The NDVI tends to saturate in dense mature forests and, thus, has limited capability to identify spatial variations, such as in the Amazon (Hall et al 2011, Myneni et al 2001). The method presented here uses a canopy height map (Simard et al 2011) that is derived from LIDAR measurements. Active remote sensing tracers (here LIDAR) have the advantage that they are highly sensitive to structural variations (Lefsky et al 2002).

NEP: We here estimate that the Amazon forest gains, on average, 0.56 Gt of carbon per year and as such yields a positive NEP. Even when considering a mean standard error of ±0.51 Gt of carbon per year due to uncertainties in the canopy height map, the Amazon is identified as a sink of atmospheric carbon. As this value is the potential carbon uptake due to the successional states of the forests, our estimate compensates approximately the amount of carbon currently released due to deforestation and land-use change in South America (Baccini et al 2012). In particular, forests in earlier successional states contribute to an uptake of atmospheric carbon (Poorter et al 2014, figure 2(d)).