Hurricane Katrina’s Carbon Footprint on U.S. Gulf Coast Forests Jeffrey Q. Chambers,1* Jeremy I. Fisher,1,2 Hongcheng Zeng,1 Elise L. Chapman,1 David B. Baker,1 George C. Hurtt2 orests recovering from land use, the expansion of woody vegetation, and other ecological processes produce a net terrestrial CO2 sink of ~1 to 2 Pg C year−1 (1). The United States contributes an estimated 0.30 to 0.58 Pg C year−1 to this global sink, with 26 to 33% being actively sequestered in forest trees (2). Changes in the strength and sign of this sink over the coming decades are difficult to predict, but as secondary forests mature the sink strength is likely to diminish (3). Another process that can diminish the terrestrial carbon sink is an increase in disturbance frequency and intensity (4), which transfers biomass from live to dead respiring pools and shifts the stem distribution toward smaller average tree size and lower biomass stocks (5). Here, we quantify hurricane Katrina’s carbon impact on Gulf Coast forests by using a synthetic approach combining detailed field investigations, remote sensing image analyses, and empirically based models for regional scaling. To develop spatially explicit maps of hurricane forest impacts, we used spectral mixture analysis (SMA) (6, 7) on Landsat imagery to quantify perpixel fractional abundance of green vegetation (GV), nonphotosynthetic vegetation (NPV: wood,

F

A

dead vegetation, and surface litter), soil, and shade for seasonally matched Landsat 5 images captured before and after the storm. The fractional change in NPV (∆NPV) from 2003 to 2006 provided a quantitative measure of the change in dead vegetation associated with Katrina. A subset for the Pearl River basin was stratified by DNPV to generate disturbance classes, and forest inventory plots were randomly established across the entire DNPV disturbance gradient (Fig. 1A). In each plot, tree mortality and damage, species composition, and biomass loss were quantified. A strong correlation between Landsat-derived DNPV and field-measured tree mortality and damage (fig. S1) enabled development of tree mortality and damage maps from the Landsat imagery. Next, a second scaling function was generated by comparing Landsat- and MODIS-derived DNPV. With the high temporal frequency and large spatial dimension of MODIS imagery, the Landsat-MODIS scaling provided an assessment of hurricane disturbance across the entire impact region (Fig. 1B). To carry out this scaling, we generated distribution functions for stem density and tree biomass from our forest inventory plots and additional U.S. Forest Service data. A Monte Carlo model was devel-

B

oped to estimate stem density, biomass distribution, mortality, and damage for all forested pixels in the MODIS scene (fig. S2) affected by Katrina. Statistically evaluated minimum and maximum values for key model parameters were used to estimate prediction error intervals (table S1). Nominal runs of the model predicted mortality and severe structural damage to 320 million large [>10-cm diameter at breast height (DBH)] trees (range from 290 to 350) with a total biomass loss of 105 Tg C (1 Pg = 1000 Tg) (range from 92 to 112), an amount equivalent to 50 to 140% of the net annual U.S. carbon sink in forest trees (2). Methods for calculating the contribution of forest trees to the terrestrial carbon sink include summing tree recruitment and growth and subtracting mortality (2). Although carbon in coarse woody debris (CWD) from tree mortality and damage is not immediately respired to the atmosphere, this CWD pulse largely represents committed future CO2 emissions (5). Although subsequent forest recovery from disturbance can offset CO2 emissions from decomposing CWD, a sustained increase in disturbance intensity or frequency (or both) will reduce forest tree carbon stocks and ultimately cause ecosystems to act as a net CO2 source (8). If a warming climate causes more extreme events and greater storm intensity (4), elevated forest tree mortality will increase CWD production, resulting in higher ecosystem respiration and a potentially important positive feedback with elevated atmospheric CO2. References and Notes 1. P. Bousquet et al., Science 290, 1342 (2000). 2. S. W. Pacala et al., Science 292, 2316 (2001). 3. G. C. Hurtt et al., Proc. Natl. Acad. Sci. U.S.A. 99, 1389 (2002). 4. Intergovernmental Panel on Climate Change (IPCC), Climate Change 2007: The Physical Science Basis—Summary for Policymakers (Cambridge Univ. Press, Cambridge, 2007). 5. J. Q. Chambers et al., Oecologia 141, 596 (2004). 6. J. Q. Chambers et al., Trends Ecol. Evol. 22, 414 (2007). 7. J. I. Fisher, J. F. Mustard, Remote Sensing Environ. 109, 261 (2007). 8. S. G. McNulty, Environ. Pollut. 116, S17 (2002). 9. M. D. Powell, S. H. Houston, L. R. Amat, N. Morisseau-Leroy, J. Wind Eng. Ind. Aerodyn. 77-78, 53 (1998). 10. Research supported by U.S. Department of Energy’s National Institute for Climatic Change Research (NICCR), NASA’s Large-Scale Biosphere-Atmosphere Experiment in Amazonia (LBA) project (CD-34), and the Long-Term Estuary Assessment Group (LEAG).

Downloaded from www.sciencemag.org on November 15, 2007

BREVIA

Supporting Online Material

Kilometers 0

2

Kilometers

4

0 -0.2

+0.1

50

-0.2

100

+0.1

+0.4

+0.4

Fig. 1. (A) Pre- to posthurricane change in the NPV fraction (DNPV) on a Landsat 5 subset for the Pearl River basin (Louisiana-Mississippi state line) provided a quantitative measure of disturbance intensity. By using this map, we established forest inventory plots (white markers) across the disturbance gradient. Open black markers represent (top) moderately resistant, infrequently flooded, bottomland hardwood forest; (middle) minimally resistant, frequently flooded, bottomland hardwood forest; and (bottom) highly resistant, flooded, cypress-tupelo swamp forest. (B) MODIS-derived DNPV from 2005–2006 provided regional estimates of tree mortality and biomass loss across the entire impact region. Isotachs (white lines) represent tropical storm (TS), category 1 (H1), and category 2 (H2) wind fields (9). www.sciencemag.org

SCIENCE

VOL 318

www.sciencemag.org/cgi/content/full/318/5853/1107/DC1 Materials and Methods Fig. S1 Table S1 References 8 August 2007; accepted 2 October 2007 10.1126/science.1148913 1 Ecology and Evolutionary Biology, Tulane University, 400 Lindy Boggs, New Orleans, LA 70118, USA. 2Institute for the Study of Earth, Oceans, and Space, University of New Hampshire, Durham, NH 03824, USA.

*To whom correspondence should be addressed. E-mail: [email protected]

16 NOVEMBER 2007

1107

www.sciencemag.org/cgi/content/full/318/5853/1107/DC1

Supporting Online Material for Hurricane Katrina’s Carbon Footprint on U.S. Gulf Coast Forests Jeffrey Q. Chambers,* Jeremy I. Fisher, Hongcheng Zeng, Elise L. Chapman, David B. Baker, George C. Hurtt *To whom correspondence should be addressed. E-mail: [email protected] Published 16 November 2007, Science 318, 1107 (2007) DOI: 10.1126/science.1148913

This PDF file includes: Materials and Methods Fig. S1 Table S1 References

Hurricane Katrina’s Carbon Footprint of Gulf Coast Forests Supporting Online Material MATERIALS AND METHODS Hurricane Katrina caused extensive forest damage ranging from downed and dead trees, snapped boles, fractured crowns, and stripped leaves. Here we briefly describe the methodology for quantifying and regionalizing these impacts including: (i) building a Landsat-resolution (30 m) ∆NPV disturbance map, (ii) sampling tree mortality and damage across the disturbance gradient using the Landsat ∆NPV map, (iii) constructing scaling relationships among field, Landsat and MODIS (500 m) data, and (iv) developing a Monte-Carlo modeling approach for regionalizing tree damage and mortality estimates. Two Landsat 5 scenes were acquired from the USGS for May 29th, 2003 and June 6th, 2006. Scenes were atmospherically corrected with ACORN software (ImSpec LLC), and geometrically corrected with state-provided road data to pixel-level accuracy. Forested regions were identified using supervised classification of the pre-Katrina scene. Landsat scenes were decomposed using an unconstrained four-endmember spectral mixture analysis (SMA) following Fisher et al. (S1), using scene-based endmembers of green vegetation (GV), non-photosynthetic vegetation (NPV), soil, and shade. ∆NPV was calculated as the difference in the NPV fraction of 2006 and 2003. Regions within the hurricane windfield which had been subjected to land-use change between 2003 and 2006 were excluded. The Landsat ∆NPV map was used to conduct stratified random field sampling across the disturbance gradient. Forest inventory plots (400 m2) were established in five randomly chosen pixels from each of five intensity classes (n = 25) and located in the field using a handheld GPS receiver accurate to ~5 m. Plot measurements included species identification, diameter at breast

1

height (DBH at 1.3 m), tree condition (dead, snapped but sprouting, damaged, leaning, etc.), height of the remaining portion of snapped boles (using a Laser rangefinder), and percentage of crown loss for damaged trees distributed among five damage classes (S2). Height of snapped trees that were sprouting (i.e. not dead) was used to estimate partial tree volume and biomass loss using methods similar to Chambers et al. (S3). Mortality and damage quantified in the 25 plots was highly correlated with ∆NPV (Fig. S1), enabling development of quantitative landscape disturbance maps. MODIS 500m 7-band reflectance (MOD09A) data were obtained from the NASA DAAC in eight day composites from January 1st, 2004 to December 31st, 2006 (138 scenes). Endmembers for GV, NPV, and soil were selected using published methods (S4), and SMA fractions were generated. To compare Landsat and MODIS ∆NPV estimates, Landsat data were aggregated to MODIS cell size, and only those MODIS pixels comprising > 75% forested Landsat pixels in 2003 were selected. These data were used to quantify the relationship between MODIS and Landsat ∆NPV for scaling field-derived mortality and damage to the entire impacted region. A novel phenological tool was employed to determine the probability that any given pixel was forested. The GV of a deciduous pixel has a high value in the summer and a low value after senescence, while evergreens maintain relatively level GV values throughout the year, and urban and non-vegetated regions have constantly low GV. When graphed on a scatterplot in fraction space as GV in summer (GVmax) vs. GV in winter (GVmin), most regions fall into a coherent ternary region, with corners represented by pure deciduous, evergreen (primarily pine), and nonvegetated behavior (Fig. S2). Deciduous and evergreen points were selected to represent pure deciduous and evergreen behavior, respectively, while areas between these two points

2

represented mixed forests. Absolute perpendicular distance from the forest line was taken as a measure of the proximity to forest, and modeled probabilistically. Data from the (i) processed MODIS imagery, (ii) additional field plots established in the Pearl River basin before the hurricane (38 plots, 1870 trees), and (iii) the USDA Forest Service’s Forest Inventory Analysis (FIA) plots (fia.fs.fed.us) were employed to simulate tree mortality, damage and biomass loss for the entire area impacted by hurricane Katrina using a Monte-Carlo approach. First, each forested MODIS pixel was assigned a number of trees based on a random draw from stem density datasets for both evergreen and deciduous forest. Next, biomass loss from tree mortality and damage was estimated using published allometric equations predicting species-specific whole tree dry-weight biomass, including estimates of partial biomass loss from snapped trees and damaged crowns (S2). Iterating through each forested MODIS pixel, the total number of dead and damaged trees, and total biomass loss were estimated for the entire impacted region. A sensitivity analysis of minimum and maximum values for 9 key parameters (Table S1) enabled constraining mortality and biomass loss values to a range of plausible results.

3

mortality and damage (fraction of stems)

0.80 Mortality Mortality plus snapped trees

0.70

0.60

0.50

0.40

0.30

0.20

0.10

0.00 0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

remote sensing metric (∆NPV)

Figure S1. Tree damage and mortality rates from the 25 forest inventory plots and Landsatderived ∆NPV. Linear regressions for tree mortality (black line: NM = 0.553·∆NPV, r2adj = 0.85, p < 0.0001) and total tree mortality plus snapped trees (i.e. severe forest structural damage) (grey line: NT = 0.996·∆NPV, r2adj = 0.88, p < 0.0001) enabled estimating tree mortality and damage from the ∆NPV remote sensing metric.

4

Figure S1. Scatter-plot showing the density of all non-water points (including non-forested) in the overall MODIS image. Point C represents maximum evergreen phenological behavior, and D deciduous behavior, with a randomly placed example point with a distance 0.13 from the forest line (well within the "acceptable" forest bounds), and a distance of 0.27 from D and 0.33 from C, identifying a forest with a predicted composition of 45% deciduous and 55% evergreen species.

5

Table S1. Principal parameters for developing the Monte-Carlo model for predicting the number of dead and damaged trees and total biomass loss in the Katrina impact region. Lmax and Lmin are the “distance to forest” from the MODIS data scatter plot (Fig. S2), describing the probability of a MODIS pixel being forested, btot is the total number of dead and snapped trees from the upper regression line of Fig. 1b, ρstem is mean stem density (trees ha-1) for evergreen and deciduous forests, Ln(Btree) is natural-log transformed mean tree biomass, fdead is the fraction of dead stems from the total dead and snapped, fsnap is the mean biomass loss from a snapped live tree (i.e. not considered mortality), fcrown is fractional increase in the number of crown damaged trees from the total number of dead and snapped trees, and a-bNPV are the slope and intercept for scaling MODIS to Landsat ∆NPV. Model runs for nominal (best estimate) values, and max and min values for each parameter, provided an error analyses of model predictions.

Parameters

Nominal Pareamter Values

Low Pareamter Values

Total Dead Snap trees trees trees Biomass Carbon (million) (million) (million) Loss (Tg) (Tg)

Total Dead Snap High Pareamter trees trees trees Biomass Carbon Values (million) (million) (million) Loss (Tg) (Tg)

Lmax

0.126

0.0707

290

100

82

187

93

0.15

332

115

94

214

107

Lmin

0.43

0.324

286

99

81

184

92

0.52

340

118

96

219

109

b tot

0.996

0.847

292

101

83

190

95

1.145

352

122

100

225

112

310

108

88

200

100

deciduous 292 coniferous 395

332

115

94

214

107

deciduous 5.52 coniferous 5.56

321

111

91

215

107

ρstem

deciduous 283 coniferous 384

deciduous 273 coniferous 372

Ln(Btree)

deciduous 5.45 coniferous 5.54

deciduous 5.37 coniferous 5.51

321

111

91

199

100

f dead

0.55

0.454

317

92

108

205

103

0.646

321

131

72

211

105

f snap

0.809

0.793

321

111

91

206

103

0.826

321

111

91

208

105

0.423

321

111

91

201

101

0.752

321

111

91

217

108

329

114

93

213

106

f crown a -b NPV nominal

0.588 intercept 0.180 slope 0.72

intercept 0.176 slope 0.69

314

109

89

204

102

321

111

91

207

104

6

intercept 0.182 slope 0.73

References S1. S2. S3. S4.

J. I. Fisher, J. F. Mustard, M. A. Vadeboncoeur, Remote Sens Environ 100, 265 (Jan 30, 2006). E. L. Chapman, Masters Thesis, Tulane University (2006). J. Q. Chambers, J. Santos, R. J. Ribeiro, N. Higuchi, Forest Ecology and Management 152, 73 (2001). J. I. Fisher, J. F. Mustard, Remote Sens Environ 109, 261 (2007).

7