Controls over aboveground forest carbon density on Barro Colorado Island, Panama

Despite the importance of tropical forests to the global carbon cycle, ecological controls over landscape-level variation in live aboveground carbon density (ACD) in tropical forests are poorly understood. Here, we conducted a spatially comprehensive analysis of ACD variation for a continental tropical forest – Barro Colorado Island, Panama (BCI) – and tested site factors that may control such variation. We mapped ACD over 1256 ha of BCI using airborne Light Detection and Ranging (LiDAR), which was well-correlated with ground-based measurements of ACD in Panamanian forests of various ages ( r2 = 0.84, RMSE = 17 Mg C ha −1, P < 0.0001). We used multiple regression to examine controls over LiDAR-derived ACD, including slope angle, forest age, bedrock, and soil texture. Collectively, these variables explained 14 % of the variation in ACD at 30-m resolution, and explained 33 % at 100-m resolution. At all resolutions, slope (linked to underlying bedrock variation) was the strongest driving factor; standing carbon stocks were generally higher on steeper slopes. This result suggests that physiography may be more important in controlling ACD variation in Neotropical forests than currently thought. Although BCI has been largely undisturbed by humans for a century, past land-use over approximately half of the island still influences ACD variation, with younger forests (80–130 years old) averaging∼15 % less carbon storage than old-growth forests ( >400 years old). If other regions of relatively old tropical secondary forests also store less carbon aboveground than primary forests, the effects on the global carbon cycle could be substantial and difficult to detect with traditional satellite monitoring. Correspondence to: J. Mascaro (jmascaro@stanford.edu)


Introduction
Tropical forests play a critical role in the global carbon cycle, storing an estimated 350 Pg of carbon in their aboveground biomass -more than any other biome (Fischlin et al., 2007).Tropical forest carbon stocks are currently declining, with losses from deforestation and degradation exceeding increases from secondary regrowth in most (but not all) regions (Asner et al., 2010b).In addition, carbon stocks in intact old-growth forest may be increasingly fluctuating due to global environmental changes, including increasing temperatures and changes to seasonality (Clark et al., 2003), increasing CO 2 concentrations (Phillips et al., 1998), and an increasing ratio of forest edge to interior area (Laurance et al., 2004).Given the pressing importance of quantifying tropical forest carbon balance, ecologists have recently improved carbon monitoring techniques in several areas, including development of generalized tree allometry theory (Chave et al., 2005;Niklas, 2006), assembly of global wood density databases (Chave et al., 2006;Swenson and Enquist, 2007), and improved methods to quantify carbon stored in poorly sampled pools such as lianas, small trees, and necromass (Hughes et al., 1999;Schnitzer et al., 2006;Chao et al., 2009).As a result, our ability to estimate the aboveground carbon storage of a given tropical forest inventory plot and quantify uncertainty -while still laborious -has improved considerably, and published estimates for many inventory plots now exist (Malhi et al., 2006;Chave et al., 2008;Lewis et al., 2009).
Lagging far behind improved field methodology, however, is our understanding of fundamental ecological controls over aboveground carbon density (ACD) at landscape scales in tropical forests.Without an improved understanding of these controls, changes in ACD will be difficult to interpret -much Published by Copernicus Publications on behalf of the European Geosciences Union.less predict -and the driving variables behind such changes will remain unclear.Within a given climate regime, physiography (the form of the land), soil properties, species composition, and disturbance history are all known to influence ACD in terrestrial ecosystems, but comparative studies of the relative importance of these driving factors have been few (but see, e.g., de Castilho et al., 2006).Further, the most comprehensive studies of variables controlling ACD variation come from oceanic islands, where steep gradients in site factors such as soil nutrients and slope angle are common (Aplet et al., 1998;Porder et al., 2005).However, these results may not be generalizable to mainland tropical forests, where the vast majority of tropical forest carbon stocks are found.In some mainland forests, landscape-scale variation in physiography and soil properties is quite limited, and explains little if any variation (e.g., Clark and Clark, 2000).Other mainland sites do exhibit ACD variation related to slope position (Valencia et al., 2009), and at broader scales, soil nutrient and textural properties (Malhi et al., 2006), and these patterns warrant further study.
One factor limiting our understanding of landscape-scale carbon storage patterns is that most work thus far has utilized ground-based field inventory plots that are limited in number and spatial extent.Such plots are invaluable for monitoring long-term community dynamics (Hubbell et al., 1999;Carson and Schnitzer, 2008), but their use for monitoring landscape-scale carbon stocks (and fluxes) has a mixed track record (see, e.g., Phillips et al., 1998;Clark et al., 2001).A single, 50-ha forest dynamics plot may contain some 250 000 trees, which provides a powerful sample size for addressing community-level questions, but the same plot lacks the spatial extent that may be needed to address landscape-level ACD patterns with sufficient statistical power.Inventory plots may also be biased toward terra firme forests on gentle slopes that have little history of recent disturbance (de Castilho et al., 2006).These selection criteria limit the utility of such plots in assessing ACD patterns at landscape scales, where habitat type, slope, and land-use history are potentially important driving factors.
Larger scale approaches have great potential to provide insights into landscape-level controls over ACD.Recently, both airborne and spaceborne Light Detection and Ranging (LiDAR) methods have turned the corner from methodological demonstrations to ecological assessment tools for ACD.Spaceborne LiDAR (linked with forest cover data) is now providing global, coarse-scale maps of forest height -a close correlate of carbon density (Lefsky, 2010), while airborne Li-DAR has provided high-resolution maps of ACD in various tropical and subtropical ecosystems types with uncertainty levels as good or better than field campaigns alone (Asner et al., 2008;Asner et al., 2010a;Asner et al., 2010b;Asner et al., 2011).Airborne LiDAR systems can be flown over thousands of hectares per day, greatly increasing the area of study for large-scale analyses and reducing habitat biases.
Here, we used high-resolution airborne LiDAR to investigate ACD over a 1256 ha portion Barro Colorado Island (BCI).We explored how ACD was related to fundamental site factors, including variation in slope angle, parent material, and soil texture.We also considered land-use history; today, BCI is almost entirely forested, but approximately half of the island was deforested in the past 100 to 200 years (Foster and Brokaw, 1996).The degree of recovery of carbon stocks over the subsequent century in these forests is of considerable interest, given past, current, and projected levels of deforestation and degradation in the tropics.To aid our interpretation of ACD patterns on BCI, we also evaluated correlates of mid-and low-canopy density as measured by LiDAR.

Study area
We conducted this study on Barro Colorado Island, Panama (BCI; 9 • 9 N, 79 • 50 W), which is part of the Barro Colorado Nature Monument administered by the Smithsonian Tropical Research Institute (STRI), with additional sampling in STRI's "Agua Salud" holdings (9 • 12 N, 79 • 49 W) adjacent to Soberanía National Park.BCI is approximately 15 km 2 , and was separated from adjacent forest and agricultural lands by the creation of Gatun Lake between 1910 and 1914 as part of the construction of the Panama Canal (Leigh, 1999).Forests on BCI are tropical moist (Holdridge, 1947), receiving approximately 2600 mm of rainfall annually with a pronounced dry season between January and April (<100 mm month −1 ).Mean annual temperature is 26 • C and about 12 % of tree species (>10 m in max height) are deciduous (Condit et al., 2001).BCI and its forests are described in detail by Leigh et al. (1996) and Leigh (1999).To ensure that our LiDAR data were calibrated to forests of all possible carbon densities on and beyond BCI (particularly at low levels not found in inventory plots on the island), we conducted additional LIDAR measurements over mapped plots in nearby Agua Salud (Neumann-Cosel et al., 2010), which contains a matrix of grasslands, active pasture, and young to moderateaged forests that have recently reverted from pasture.

LiDAR data collection, calibration, and mapping
In September 2009, we used the Carnegie Airborne Observatory (CAO; Asner et al., 2007) to collect LiDAR data over central Panama, including 98 % of BCI (i.e., less small clouds) and 100 % of Agua Salud.The LiDAR data were collected 2000 m above ground level, with 1.12-m spot spacing, 30-degree field of view, and 50-kHz pulse repetition frequency, for which the aircraft maintained a ground speed ≤157 kph.From these data, we created ground and surface digital elevation models, and measured vegetation height as the difference between the two (follows Asner et al., 2009).
We previously quantified LiDAR spatial errors at less than 0.15 m vertically and 0.36 m horizontally (root mean square error, Asner et al., 2007;Asner et al., 2010b).We also analyzed the vertical profile of vegetation by binning LiDAR returns into volumetric pixels (voxels) of 5-m spatial resolution and 1-m vertical resolution.By combining all voxels in each 5-x-5-m spatial cell, we created vertical histograms representing the full vertical spread of LiDAR returns.These data were broken down to mean canopy profile height (i.e.MCH, the vertical "center" of the canopy; see Lefsky et al., 1999), a simple metric to be linked to ground data.MCH was slightly but consistently better than other predictors of ACD, including top-of-canopy height (also known as RH100, Lefsky, 2010) and quadratic mean canopy height, confirming previous results from other tropical forests (Asner et al., 2010b;Asner et al., 2011).We also summed relativized LiDAR returns over five, 5-x-5-m spatial cells between 12 and 17 m in height to examine structural variation in a mid-canopy band, and five, 5-x-5-m spatial cells between 2 and 7 m in height to examine structural variation in a low-canopy band.We refer to these as "mid-canopy" and "low-canopy" density, terms that are relative to island-wide vegetation (with a mean topof-canopy height of ∼22 m) rather than the local context.We used standardized heights rather than percentiles of the local canopy because our aim was to capture differences in the density of vegetation at particular absolute heights.
To calibrate the LiDAR MCH data, we used vegetation data from the 50-ha forest dynamics plot on BCI (Condit, 1998;Hubbell et al., 1999;Hubbell et al., 2010), as well as twenty-nine, 20-x-50-m (1000 m 2 ) belt transects in Agua Salud.For each site, the available data included diameter measurements and identification to species for all live trees (i.e., lianas were excluded) ≥ 1 cm diameter at breast height (1.3 m) or above buttresses, if present (hereafter dbh).The 50-ha plot data were collected in 2010, while the Agua Salud data were collected in 2009.For Agua Salud, where we targeted very young to medium-aged secondary forests, we used a single, locally generated allometric model to estimate aboveground biomass for all stems based on dbh and wood density (Table 1).For the 50-ha plot, we used a generalized allometric model for tropical moist forests based on dbh, height, and wood density provided by Chave et al. (i.e., "Model I", 2005).Height measurements were available for 1218 trees and these data were used to create a diameter-toheight allometric model to estimate tree height (Wright and Muller-Landau unpublished data).We did not assign actual height measurements to particular stems in the 50 ha plot because considerable time has elapsed since the measurements were taken.Wood density data were a combination of locally sampled estimates and literature values (Chave et al., 2006, Wright et al. unpublished data).Carbon content was calculated as 48 % of dry aboveground biomass (IPCC, 2006).ACD was calculated per unit horizontal area (not per unit surface area); this is the standard means of comparing carbon and biomass densities in space (Dallmeier, 1992;Con-dit, 1998;Clark and Clark, 2000;de Castilho et al., 2006), including in long-term forest inventory plots such as the 50ha forest dynamics plot (Chave et al., 2003).For the Agua Salud data, plots were delineated along the ground and thus carbon stocks were corrected using slope angle to units of ACD in horizontal area.
In August 2010, we collected GPS points at the corners of the 50-ha forest dynamics plot and endpoints of each Agua Salud transect using a Leica GS50 receiver (Leica Geosystems, Inc.; St. Fallen, Switzerland).All points were differentially corrected to known base stations, in most cases achieving <1-m positional error.Using a geographic information system (GIS; ArcView 9.3, ESRI, Redlands, CA), we delineated the plot boundaries.For young forests at Agua Salud, we treated each 20-x-50-m belt transect as a separate plot, while on BCI, we delineated 128, 60-x-60-m plots by combining nine adjacent 20-x-20-m quadrats as mapped in the forest dynamics plot.Differing plot sizes between young and old forests were used to limit heteroscedasticity in the LiDAR-to-carbon relationship; in previous work, such relationships exhibited increasing ACD variance at increasing levels of ACD for constant size plots (i.e., in older forests, Asner et al., 2010b;Asner et al., 2011).In all, 157 plots comprising 47 ha were used.We regressed ground-based ACD against LiDAR MCH (see previous section) in the form of a power model: where x is MCH, y is ACD, and a, b, and k are model parameters.A non-arithmetic error term (i.e., x k * ε) was used to account for heteroscedasticity that typifies MCH-to-ACD relationships (e.g., Asner et al., 2010b).This approach is analogous to fitting a linear model to log-transformed variables, but avoids the need for back-transformation (Baskerville, 1972;Mascaro et al., 2011).The model was fit using maximum likelihood in the R programming language (v.9.2, R Development Core Team 2009).Finally, we applied the model to our MCH data for BCI to produce maps of ACD at spatial resolutions from 30 to 100 m (see below).

Spatial analysis of island-wide variation in carbon density
We considered the influence of slope, forest age, bedrock type, and soil texture on ACD on BCI (Fig. 1) through backwards stepwise multiple regression, as detailed below.

Model Equation Reference
Agua center pixel.At coarser resolutions of 30 to 100 m, slope was assigned to each pixel as the average slope of all 1.12-m pixels.
Forest age was assigned based on a map of forest types drawn by Enders (1935), interpreted through both field survey notes and observations of a 1927 US Army Air Service aerial photograph.The map contains five classes of forest as they appeared to Enders in the early 1930s: (1) clearings, (2) recently reverted forest no more than 30 years old, (3) forest 40-50 years old, (4) primeval forest, and (5) primeval forest with palms, with classes 2 and 4 accounting for most of the forest on the island.Based on these classes, we delineated three basic age categories as of 2009: (1) secondary forest between 80 and 110 years old, which included recently reverted forest plus clearings that were allowed to revert shortly after they were visited by Enders (such as the radio tower clearing near the center of the island); (2) secondary forest between 120 and 130 years old, which included all of the original 40-to-50-year-old class; and (3) old-growth forest, both with and without palms, which is generally agreed to be >400 years old, although estimates range from 200 to >1500 years old (Foster and Brokaw, 1996;Leigh, 1999).Two clearings have been mostly maintained (the laboratory and a canal light on the northern peninsula), and these were Poorly drained excluded from our analysis.Forest age classifications were originally generated in vector format and were assigned to a given pixel at spatial resolutions from 30 to 100 m based on the value in the center of that pixel.Soil type and bedrock were assigned based on a local soil survey commissioned for BCI (Baillie et al., 2006), which provided adjustments to bedrock geology maps by Woodring (1958) and Johnsson and Stallard (1989), as well as a map of 13 local soil types.Bedrock within the Isthmus of Panama is the product of a volcanic arc produced by the subduction of the Nazca plate beneath the Caribbean plate, which began 13 Ma ago and is ongoing.On BCI, there are three primary formations: the Bohio (early Oligocene) is conglomerate with basaltic pebbles, cobbles, and boulders in a matrix of finer basalts.Above the Bohio lies the Caimito formation (late Oligocene), which is delineated into marine and volcanic facies.The marine facies, primarily found on the western side of BCI, are sorted fossiliferous sandstone of various grain sizes interbedded with foraminiferal limestone.The volcanic facies of the Caimito formation are found on the eastern side of the island and are basaltic conglomerate including boulders.Finally, a large, extrusive andesite "cap" (Oligocene and early Miocene) is found near the center of the island, where it forms a plateau that underlies much of the 50-ha forest dynamics plot.In place of the 13 local soil types delineated by Baillie et al. (2006), we used four textural classes: brown fine loam, red light clay, pale swelling clay, and heavy clay (Table 2).Finer discriminations into the 13 local types were based on the intersection of each textural type with an underlying bedrock feature, and were thus redundant with included bedrock variables.We also excluded two textural classes that were not mapped by the survey due to their small spatial extent.As with forest age, soil and bedrock classes were originally generated in vector format and were assigned to a given pixel at spatial resolutions from 30 to 100 m based on the value in the center of that pixel.
We used backwards stepwise multiple regression to assess carbon density controls on 1256 ha of BCI; we excluded areas covered by small clouds and areas within 50 m of a shoreline.The individual data points represented separate square pixels that covered the available area, including only complete pixels.Four possible variables were included in our model: slope angle (as a continuous variable), forest age classes (as a categorical variable with three classes), bedrock (four classes), and soil textural (four classes).For each categorical variable, the class with the largest spatial extent was considered the base class (i.e., old-growth, Bohio, and brown fine loam, respectively); thus, regression coefficients for other classes denote differences from the base class (sensu Suits, 1957).The main model began by fitting the multiple linear regression using all variables, and then re-fitting with each variable removed to determine which removed variable lead to the greatest decline in model fit (as assessed by the Akaike Information Criterion or AIC, Akaike, 1974).This process was repeated until the AIC stabilized.We repeated the procedure to examine controls on mid-canopy and lowcanopy density.
We considered the influence of spatial scale and spatial variation on our model results in a number of ways.First, to assess the influence of spatial resolution (i.e., pixel size) we compared model results at multiple resolutions (from 30 to 100 m as discussed above).Second, to reduce the potential influence of lack of spatial independence among our data points (all pixels in our study area), we conducted a rarefaction analysis using just 10 % of the data and repeated 1000 times at the highest and lowest spatial resolutions www.biogeosciences.net/8/1615/2011/Biogeosciences, 8, 1615-1629, 2011 considered (30 and 100 m).Finally, we mapped relativized model residuals to allow for visual examination of spatial patterns in the residuals (at 30-and 100-m resolution).All analyses were conducted in the R programming language (v.9.2, R Development Core Team 2009).

Results
Airborne LiDAR revealed wide variation in forest canopy height across the island (Fig. 2a).Taller canopies (red and orange) were found in three primary areas: (1) a series of folded ridges and valleys along the north side of the island (stretching from the Standley to Fairchild peninsulas), (2) a band along the fault line (see Fig. 1b) beginning near the laboratory and extending to the southern peninsula, and (3) the southern portion of the andesite cap at the center of the island.These are all areas of older forest (see Fig. 1d).By contrast, shorter canopies (blue and blue-green) are concentrated on both the east and west sides of the island, near the shoreline, and on many of the peninsulas -mostly areas of younger forest.Throughout these shorter canopy regions, clumps of emergent trees are visible (isolated red dots encircled by orange/yellow).In general, the canopy is taller on steeper slopes and in older forests.Mid-canopy (12-17 m) density is generally highest in areas of low canopy height (compare Fig. 2a, b), with high mid-canopy density areas on the eastern side of the island, a hill just northeast of center, the northern peninsulas, and the western edge.Low-canopy (2-7 m) density is largely concentrated in small, very dense patches that track areas of low canopy height (compare Fig. 2a, c).These dense low-canopy patches are found mostly on the western side of the island, especially the western shoreline, with more diffuse scattering across the central hill and plateau, as well as the western slope.Other concentrations include the laboratory clearing and a feature on the spine of the Miller peninsula.
LiDAR-derived MCH strongly predicted ACD as measured in field inventory plots (r 2 = 0.84, RMSE = 17 Mg C ha −1 , P < 0.0001; Fig. 3), and the relationship was consistent across sparse pastures and very young secondary forests in Agua Salud, as well as mature forest or patches in various stages of gap recovery within the 50-ha forest dynamics plot (Fig. 3).The data were heteroscedastic; errors when considering the BCI 50-ha plot data alone were 19 Mg C ha −1 (RMSE) compared to 8 Mg C ha −1 for the Agua Salud data.
ACD was higher on steeper slopes and in older forests (Table 3), and was significantly related to slope, forest age, bedrock, and soil texture (Table 4).Comparisons between observed ACD, predicted ACD, and residual ACD are shown in Fig. 4. At 30-m resolution, the multiple regression model explained 14.7 % of the variation in ACD on BCI (F 9,13948 = 267.2,P < 0.0001; Table 4, Fig. 5a).The percentage of variation explained by the model increased as spatial resolution decreased (i.e., with increasing pixel size); at 100-m resolution, 33.1 % of the variation was explained (F 9,1038 = 58.8,P < 0.0001).At all resolutions examined, slope was the most significant variable in the model, explaining 7.5 % of the variation at 30-m resolution and 19.1 % at 100-m resolution (Fig. 5a).Forest age was generally the next most significant variable, ranging between 3.1 and 6.8 % of variation explained.Throughout BCI, secondary forests (80-130 years old) store 10-20 % less carbon than do oldgrowth forests (>400 years old), depending on slope angle (Table 3).Bedrock and soil texture had more limited overall influence (Fig. 5a), but both Caimito bedrock types as well as pale swelling clay were associated with lower ACD (Table 4).The fitted model reproduced the broad island-wide patterns in ACD (compare Fig. 4a v. 4c,4b v. 4d); some spatial autocorrelation was evident in the residuals (Fig. 4e, f).
The factors controlling ACD variation remained consistent in both sign and magnitude following rarefaction of the dataset, while declining across the board in significance level due to declining sample size (Table 4).The residuals of the spatial model highlighted areas in which ACD departed from the model predictions (Fig. 4e, f).The model underpredicted ACD in forests inland from Shannon cove, and a strip of forest just inland from the Standley peninsula (see landmark names in Fig. 2).In contrast, the model overpredicted ACD in a line of forest on the spine of the Drayton peninsula, the area north of the laboratory, the hill near the island's center, and several patches on the western slope.Correlations were common among the independent variables examined, but only in two cases were variables outside of mutually exclusive categories (e.g., forest age) related by a Pearson's R>0.5 (Table 5).The influence of model parameters on mid-canopy density was of comparable strength, and essentially inverted in sign relative to the controls over ACD.At 30-m resolution, the model explained 16.6 % of the variation in mid-canopy density (F 9,13948 = 310.6,P < 0.0001), and this increased to 39.4 % at 100-m resolution (F 9,1038 = 76.6,P < 0.0001) (Fig. 5b).Mid-canopy density was lower on steeper slopes, and was higher in younger forests and on Caimito (volcanic) bedrock (Table 4).Controls over low canopy density were markedly different.Only 8.8 % of variation was explained at 30-m resolution (F 9,13948 = 150.5,P < 0.0001), and this rose to 23.8 % at 100-m resolution (F 9,1038 = 47.7,P < 0.0001; Fig. 5c).Bedrock type and soil texture were the strongest drivers of low-canopy density, while the influence of forest age was very weak and was not retained in the model at 100-m resolution (Fig. 5c).Pale swelling clay in particular was associated with 37 % higher low-canopy density than for forests on brown fine loam (Table 4).However, bedrock and soil texture (with the exception of pale swelling clay) were no longer significantly associated with low-canopy density after rarefaction.

Influence of physiography on carbon density
Our model indicated that slope angle was the single most important parameter controlling aboveground carbon density on the island.Moderate and steep slopes support higher ACD www.biogeosciences.net/8/1615/2011/Biogeosciences, 8, 1615-1629, 2011 than gentle slopes -a pattern that was consistent for both young and old forest (Table 3).In part, this result can be explained by greater access to space for trees within sloped areas; these areas have greater surface area and thus greater soil and light resources over a given spatial extent (though this effect is nonlinearly related to slope).However, even after adjusting for surface area, slope remained a significant predictor of ACD (i.e., it accounted for 11 % -down from 19 % -of all surface-area-adjusted ACD variation at 100 m resolution, model-wide F 9,1038 = 42.0,P < 0.0001; Supplement).Taken together, the results suggest that physiography (the form of the land) is the primary driving factor behind ACD variation on BCI.We propose three non-mutually exclusive hypotheses to explain the greater carbon density on steeper slopes: nutrient availability, water availability, and competition for light.
Nutrient and water availability are driven by linkages between slope, bedrock, and soil texture and these warrant careful examination.Inclined slopes generally lead to greater rates of sediment transport than weathering, whereas gentle slopes accumulate highly weathered material -and these patterns may have consequences for plant growth.On BCI, Johnsson and Stallard (1989) found that slope angle, not the mineralogy of the underlying bedrock, was the principal Table 4. Fitted coefficients of multiple regression analysis of variables controlling carbon density variation on BCI, Panama.With the exception of slope, all variables are categorical and thus coefficients reflect predicted deviation from the base classes (these are "Old-growth" in forest type, "Bohio" in bedrock type, and "Brown Fine Loam" in soil type).The primary model utilized all areas inside a 50-m shoreline mask on BCI, shown here at two resolutions: 13 959 0.09-ha pixels (30 m), and 1048 1.0-ha pixels (100 m).The rarefaction analysis repeated the model 1000 times using 10 % of all pixels at each resolution; for this analysis, median P-values from all runs are presented.

Component
All  factor controlling the supply of cation-rich weathering products.They also note that bedrock composition was ultimately behind the present-day topographic relief on the island.Once erosion penetrated the plateau, the underlying Bohio formation weathered and eroded rapidly, a process that continues to yield smectitic clays with abundant cation exchange complexes.On shallow slopes, weathering proceeds in place, resulting in deep (>2 m), oxic soils with high Al content.
Several studies have demonstrated that slope angle can affect primary mineral nutrient release, and that such effects may be measureable in the vegetation (Tanner, 1977;Scatena and Lugo, 1995;Chen et al., 1997;Porder et al., 2005), although most of these studies have been conducted on islands where tree diversity is somewhat (or much) lower than mainland forests.One Amazonian study found no effect of slope on total carbon stocks ( de Castilho et al., 2006), although they found that emergent trees were less abundant -and smaller trees more abundant -on steeper terrain, suggesting that forests in sloped areas were more exposed to periodic mass wasting and canopy damage.In addition to nutrient release, water availability appears to be higher on sloped terrain on BCI relative to the hilltop and plateau, particularly during the dry season.Hubbell and Foster ( 1983) found higher soil water potential on sloped terrain on BCI, as well as higher leaf water potentials.Daws et al. (2002) also found that sloped areas consistently had higher water availability than flat areas.They further found that sloped areas experienced less extreme seasonality than flat terrain because the increased water availability mitigated the effects of the dry season.All sampled sloped areas experienced <50 drought days per year (defined as those with a soil water potential of less than −1500 kPa), while more than half of sampled plateau areas experienced >50 drought days per year, and several experienced >80 drought days.
Finally, competition for light may also drive trees to attain greater heights on steeply sloped terrain -particularly along gully walls or bottoms, where neighboring tree canopies on ridges have easier access to light (Osada et al., 2004;Lang et al., 2010).
Our study provides the strongest evidence to date of physiographic control over carbon density variation for a Neotropical forest site.Some studies have suggested that physiographic variation may play little role in ACD variation in mainland forests primarily because such variation tends to be minor (Clark and Clark, 2000).While BCI is more finely dissected than some tropical forest sites, the level of physiographic variation is by no means extreme, and yet we found that this variation was a primary controlling variable on ACD.
We are confident that our results on the effects of slope do not reflect sampling artifacts related to LiDAR.LiDAR can produce spurious estimates of vegetation height on cliff faces and in slumping terrain, where the distance between a portion of a tree's canopy and the ground DEM -while real -does not reflect the height of that vegetation allometrically.This effect can manifest on steeply sloped terrain as well, because downslope portions of tree canopies add to top-of-canopy height, while upslope potions are buried below other canopies and thus do not subtract from top-ofcanopy height.Importantly, however, we used MCH rather than top-of-canopy height to estimate ACD, and this variable acts quite differently.In the case of MCH, upslope portions of tree canopies that are lower to the ground cause MCH to decrease and balance the increase produced by downslope canopies.The slope effect reported here matches the magnitude of that reported by Chave et al. (2003) based on steep versus level terrain in the 50-ha plot.Finally, the influence of slope extended to pixels with a slope of just 10 degrees (Fig. 6).

Land use history and recent disturbance
Although most secondary forests on BCI are relatively old (>80 yr), we found that forest age remains an important controlling variable over ACD patterns.On gentle slopes (0-6 degrees), mean ACD in secondary forests 80 to 130 years old was 83 % of that in old-growth forests (79 v. 96 Mg C ha −1 ), and 86 % on moderate slopes (6-15 degrees; 85 v. 98 Mg C ha −1 ; Table 3).A similar difference was reported by Chave et al. (2003) within the 50-ha forest dynamics plot, where secondary forest ACD averaged 79 % of that in old-growth forests (103 v. 131 Mg C ha −1 ). Brown and Lugo (1990) found that 80 years was generally not sufficient time to recover ACD equivalent to primary forests in tropical moist forests, and our analysis suggests that such a disparity has persisted even longer on BCI.Most secondary forests on BCI are closer to 130 years old according to assessments by Enders (1935) and other naturalists (citiations in Foster and Brokaw, 1996).A matrix of many ages is most likely, and in this case Foster and Brokaw (1996) suggest that the oldest patches of secondary forest may be 200 years old.One region that can be dated with greater certainty is the small hill just north of the center of the island, which was the location of a French distillery that was abandoned in 1910.One hundred years later, this area is clearly visible as a low carbon density section, with dense mid-and low-canopy vegetation (Figs. 2, 4).
There are relatively few studies of carbon storage in secondary forests ≥80 years old for comparison with our estimates for BCI, and comparisons should be considered in the context of climate and edaphic factors, proximity to primary forest, land-use history, and even differing allometries.Still, our results suggest that the recovery has been slower than what might be predicted for Central America, where soils (including those on BCI) are relatively fertile compared to the global tropics.In the wetter (>4000 mm) Los Tuxtlas region of Mexico, for example, Hughes et al. (1999) projected that secondary forests would recover to 80 % of mature forest carbon storage after 25 to 60 years (mean 55 years), depending on the duration of previous land-use activity.Secondary montane wet forests (∼2000 mm) in Puerto Rico grew even faster -far exceeding primary forest carbon storage after 80 years (125 v. 77 Mg C ha −1 ) (Marín-Spiotta et al., 2007).In contrast, 80-year old secondary moist forests (∼3500 mm) on very poor soils in the Northwestern Amazon recovered to only 70 % of primary forest carbon stocks (86 vs. 122 Mg C ha −1 ) (Saldarriaga et al., 1988).However, other poor sites have exhibited faster recovery of carbon stocks, such as moist forests in Rondonia and Cameroon which respectively recovered to ∼60 % and ∼70 % of primary forests values after <20 years (Alves et al., 1997;Kotto-Same et al., 1997).In combination, these studies highlight the variability in recovery rates among sites, as well as the strong deceleration in recovery with time; the first 50 % is regained fairly quickly, but full recovery may take much longer.
Although not included in our explanatory model, our results support a pronounced role of forest disturbance interval affecting ACD as well as mid-and low-canopy structure.For instance, we found that low-canopy density was lower on steep slopes (Table 4, Fig. 2c), even though steeply sloped areas are often thought of as having higher rates of gap formation (see, e.g., de Castilho et al., 2006).With only one LiDAR flight, we cannot test controls over the rate of gap formation, but our static sample suggests that lowcanopy density is highest on the gentle, west slope of the island where pale swelling clays are abundant.Foster and Brokaw (1996) suggested that treefalls may be more common on the pale swelling clays because shrink-swell activity (due to high montmorillonite content) would destabilize soils.In fact, we found that the presence of pale swelling clays was strongly positively associated with low-canopy vegetation density, and, aside from slope, the only parameter to remain significant at coarse scale after rarefaction of our dataset (Table 4).An alternative, but not mutually exclusive explanation, is that the west side of BCI is more exposed to damaging wind events.Foster and Brokaw (1996) cite an unpublished report suggesting that the creation of Gatun Lake (1910)(1911)(1912)(1913)(1914) may have altered wind patterns and increased wind speeds, and suggested that "the periodic blasting of the half-grown regeneration on the west side of the island appears to be a new situation."Since Ender's (1935) survey of forest age, several large blowdowns have indeed impacted BCI's western forests, and some of these are visible in the LiDAR-derived aboveground carbon map.Notably, one large blowdown affected an area on the west side of the island (south of the Standley peninsula) in 1957 (Foster and Brokaw, 1996).Another large area with ACD lower than predicted is along the ridge of the Drayton trail just north of the southern peninsula and clearly tracks the path of a large blowdown that occurred there in 1989 (S. J. Wright, personal communication, 2010, Fig. 4a).

LiDAR estimation of carbon density
We were able to predict ACD with high confidence using LiDAR-derived MCH in forests spanning many different ages and stages of recovery from treefall gaps (Fig. 3).This result, along with other studies (e.g., Asner et al., 2008), demonstrates the ability of locally-calibrated LiDAR metrics to estimate ACD in tropical forests.The overall LiDARto-carbon model error (∼ 18 %) is comparable to allometric uncertainties in plot-based approaches (Chave et al., 2004).When the model is applied over a large region, relative error declines rapidly with increasing sample area, because total errors increase more slowly than total carbon (Asner et al., 2010b).Thus, we can estimate total carbon stores on BCI with much higher confidence than for a single plot; overall, we estimate that interior BCI forests (>50 m from the shoreline) store 121.2 ± 0.6 Gg C in live aboveground tree biomass.
The signal strength of the variables driving ACD variation depended to a high degree on the spatial resolution of our analysis, as we would expect.At 30-m resolution, the random inclusion or exclusion of a single large tree (in whole or in part) may change the true ACD for that pixel by an order of magnitude.At lower spatial resolution (i.e., larger pixel size), this effect declines rapidly because small-scale variation in forest structure caused by gaps, tree spacing, and standing dead trees balances, and thus there is lower overall variation to be explained by the model.
We caution that LiDAR estimates of ACD are inherently dependent on the underlying allometric equations used to estimate aboveground carbon in field inventory plots.We used a combination of local and general models that we believe are appropriate for the study region, but other models would give different results, and any differences would be additive at a regional scale.As new allometric models are developed, or wood density values for BCI tree species are refined (e.g., Williamson and Wiemann, 2010), the gross estimates presented here may require adjustment.We note, however, that any adjustments would likely have little effect on our conclusions concerning drivers of ACD, as changes would mostly increase or decrease all points in concert, thus changing the slope of the LiDAR-to-carbon relationship but not the relative position of points.An additional consideration is the possibility that tree diameter-to-height relationships may differ in steep versus flat areas, which could cause comparisons between LiDAR-derived carbon stocks on steep versus flat areas to be biased.Available tropical tree allometries are not sufficient to address this hypothesis at present (see Feldpausch et al., 2010); on BCI, however, generalized diameterto-height relationships from flat areas in the 50-ha plot and ridges near six radio towers did not differ (Larjavaara and Muller-Landau unpublished data).

Implications
Physiography may be more important in controlling ACD variation in mainland tropical forests than currently thought.Present discussion of controls over ACD variation is focused largely on species composition, namely the role of wood density and growth form variation (Laurance et al., 2004;Bunker et al., 2005;Schnitzer and Carson, 2010;Schnitzer and Bongers, 2011); this variation (and the prospect of community composition changing with increased environmental change) is clearly important, but we suggest that future research should also consider physiographic variation as a proximate control over baseline ACD variation upon which future trends will act.Additionally, the location of field inventory plots, forest reserves, and cultivation activity are all generally biased toward gentle slopes (Malhi et al., 2006), and these biases may have implications for carbon stocks measured, protected, and released to the atmosphere.
Our results also suggest that increased attention is needed on the carbon stocks of relatively old secondary forests (> 80 yr).If large regions of these forests store significantly less aboveground carbon than primary forests (as is the case on BCI), the effects on the global carbon cycle could be substantial.Secondary forests are increasing in abundance in tropical regions (as a portion of all forests), and are expected to increase further in countries that undergo increased urbanization (Wright and Muller Landau, 2006).When these forests are very young, their effects on carbon storage at the regional scale should be relatively accessible through satellite monitoring.But old secondary forests on BCI, where tree cover is similar to that in primary forests, appear to be in a long, slow climb toward primary forest carbon densities.Quantifying carbon variation caused by such old secondary forests in other (and larger) regions should be a high priority.

Fig. 1 .
Fig. 1.Spatial variation in site factors on Barro Colorado Island, Panama (BCI): (a) slope was calculated at 1.12-m resolution from the LiDAR ground DEM, and is shown here at 5-m resolution (it was coarsened by averaging to 30-m resolution and above for input into the multiple regression model); (b) bedrock geology, mapped by Woodring (1958) with updates by Johnsson and Stallard (1989) and Baillie et al. (2006); note strike-slip fault separating the Bohio and Caimito volcanic bedrock features; (c) soil textural types mapped by Baillie et al. (2006); and (d) forest age, based on land-use history mapped by Enders in 1935.

Fig. 4 .
Fig. 4. Observed LiDAR-derived carbon density at (a) 30-and (b) 100-m resolutions, (c, d) carbon density as predicted by multiple linear regression analysis, (e, f) relativized model residuals.The color ramp for carbon density reflects a minimum and maximum of one standard deviation from the mean observed carbon density on BCI (96.5 Mg C ha −1 ).

Fig. 6 .
Fig. 6.Change in relative amount of carbon density variation explained by slope angle as influenced by the maximum slope angle considered in the multiple regression model at 30-m pixel resolution.

Table 2 .
Soil textural types as delineated during a local soil survey byBallie et al. (2006).Textural types are also mapped in Fig.1c.Stony clay and fine loam over weathered saprolite.Clay fraction is dominated by kaolinite with substantial montmorillonite.These soils are immature with abundant exchangeable Ca and Mg but low K and little Al.Shallow mottled clay Less weathered analogues of the pale clays, with high montmorillonite content.Very high exchangeable Ca and Mg and low Al, plus higher K content then other BCI soils.

Table 3 .
Mean and standard deviation (SD) aboveground carbon density (Mg ha −1 ) values for different slope and forest age classifications on BCI, Panama.

Table 5 .
Correlation coefficients (Pearson's R) for model variables examined for controlling carbon density on BCI, Panama.Grey areas indicate mutually exclusive categorical variables (forest age, bedrock, soil texture).Only correlations significant at P < 0.0007 are shown (after Bonferroni correction for 66 comparisons).Correlations of R > 0.5 are highlighted in bold.