Terrestrial Ecosystem Carbon Fluxes Predicted from MODIS Satellite Data and Large-Scale Disturbance Modeling

The CASA (Carnegie-Ames-Stanford) ecosystem model based on satellite greenness observations has been used to estimate monthly carbon fluxes in terrestrial ecosystems from 2000 to 2009. The CASA model was driven by NASA Moderate Resolution Imaging Spectroradiometer (MODIS) vegetation cover properties and large-scale (1-km resolution) disturbance events detected in biweekly time series data. This modeling framework has been implemented to estimate historical as well as current monthly patterns in plant carbon fixation, living biomass increments, and long-term decay of woody (slash) pools before, during, and after land cover disturbance events. Results showed that CASA model predictions closely followed the seasonal timing of Ameriflux tower measurements. At a global level, predicting net ecosystem production (NEP) flux for atmospheric CO2 from 2000 through 2005 showed a roughly balanced terrestrial biosphere carbon cycle. Beginning in 2006, global NEP fluxes became increasingly imbalanced, starting from -0.9 Pg C yr-1 to the largest negative (total net terrestrial source) flux of -2.2 Pg C yr-1 in 2009. In addition, the global sum of CO2 emissions from forest disturbance and biomass burning for 2009 was predicted at 0.51 Pg C yr-1. These results demonstrate the potential to monitor and validate terrestrial carbon fluxes using NASA satellite data as inputs to ecosystem models.


Introduction
The emission of CO 2 from deforestation and other land cover changes is among the most uncertain components of the global carbon cycle.Inconsistent and unverified information about global deforestation patterns has significant implications for balancing the present-day carbon budget and predicting the future evolution of climate change.A number of studies have estimated carbon emissions from tropical deforestation [1][2][3][4][5][6][7] but the estimates vary greatly and are difficult to be compared due to differences in (land cover) data sources, estimated regional extents, and carbon computation methodologies.
A recent review of previous work on estimating carbon emissions from vegetation change by Ramankutty et al. [8] pointed to the importance of considering ecosystem dynamics following land cover conversion, including the fluxes from the decay of products and slash pools, and the fluxes from either newly established agricultural lands or regrowing forest.This review also suggested that accurate carbon-flux estimates should consider historical land-cover changes over at least the previous 20 years.Such results can be highly sensitive to estimates of the partitioning of cleared carbon into instantaneous burning vs long-time scale dead woody pools.Accordingly, the main objective of the present study was to quantify the major controls on carbon flux patterns and processes terrestrial ecosystems worldwide, using NASA satellite data products to drive models of net ecosystem production (NEP) and detect large-scale ecosystem disturbance, leading to detailed estimates of net biome production (NBP).
An updated version of the CASA (Carnegie Ames Stanford Approach) model [9] was used in this study to predict terrestrial ecosystem fluxes using MODIS collection 5 of the Enhanced Vegetation Index (EVI; [10,11]) as inputs at a geographic resolution of 0.5˚ latitude/longitude.This model version was developed at NASA Ames Research Center [2,9,12,13] to estimate monthly

Satellite Data Inputs
The launch of NASA's Terra satellite platform in 1999 with the moderate resolution imaging spectroradiometer (MODIS) instrument on-board initiated a new era in remote sensing of the Earth system with promising implications for carbon cycle research.Direct input of satellite vegetation index "greenness" data from the MODIS sensor into ecosystem simulation models can be used to estimate spatial variability in monthly net primary production (NPP), biomass accumulation, and litter fall inputs to soil carbon pools.Global NPP of vegetation can be predicted using the relationship between leaf reflectance properties and the fraction of absorption of photosynthetically active radiation (FPAR), assuming that net conversion efficiencies of PAR to plant carbon can be approximated for different ecosystems or are nearly constant across all ecosystems.
Whereas previous versions of the NASA-CASA model [9] used a normalized difference vegetation index (NDVI) to estimate FPAR, the current model version instead relies the EVI time series, which has a higher dynamic response across the full range of vegetated cover, does not saturate in medium-to-high biomass areas, and is less susceptible to atmospheric interference [10,14,15].The lower level of saturation of low-to-medium range plant production estimated from CASA modeling with MODIS EVI inputs (compared to MODIS NDVI inputs) will result in lower annual NPP in any zones where primary forest has become increasing mixed with degraded forest and converted agricultural land uses.
The global 0.5˚ (latitude/longitude) resolution MODIS vegetation index (VI) data sets used as inputs to CASA were generated by aggregating monthly 0.05˚ (~6 km) data (MOD13C2 version 005) from the USGS LP DAAC.The VI layer was selected from each MOD13C2 spatial composite file and surface water values are converted to "NoData".To aggregate from a 0.05˚ cell size to 0.5˚ resolution, the VI values for each pixel block were averaged.Each monthly layer was then multiplied by 0.0001 to scale the data to the standard MODIS VI value range.This aggregation procedure provided the greatest assurance of high-quality, cloud-free VI inputs to the carbon cycle model.

CASA Ecosystem Carbon Fluxes
As documented in Potter [2] the monthly NPP flux, defined as net fixation of CO 2 by vegetation, is computed in CASA on the basis of light-use efficiency [16].Monthly production of plant biomass is estimated as a product of time-varying surface solar irradiance, Sr, and EVI from the MODIS satellite, plus a constant light utilization efficiency term (e max ) that is modified by time-varying stress scalar terms for temperature (T) and moisture (W) effects (Equation (1)).The e max term is set uniformly at 0.55 g C MJ −1 PAR, a value that derives from calibration of predicted annual NPP to previous field estimates [12].This model calibration has been validated globally by comparing predicted annual NPP to more than 1900 field measurements of NPP [12,13,17].Climate drivers for the CASA model were from the National Center for Environmental Prediction (NCEP) reanalysis data set (version NCEP/DOE II and may eventually leave the system as seepage and runoff.Freeze-thaw dynamics with soil depth operate according to the empirical degree-day accumulation method [24], as described by Bonan [25].The T stress scalar is computed with reference to derivation of optimal temperatures (T opt ) for plant production.The T opt setting will vary by latitude and longitude, ranging from near 0˚C in the Arctic to the middle thirties in low latitude deserts.The W stress scalar is estimated from monthly water deficits, based on a comparison of moisture supply (precipitation and stored soil water) to potential evapotranspiration (PET) demand using the method of Priestly and Taylor [20].
Based on plant production as the primary carbon and nitrogen cycling source, the NASA-CASA model is designed to couple daily and seasonal patterns in soil nutrient mineralization and soil heterotrophic respiration (R h ) of CO 2 from soils worldwide.Net ecosystem production (NEP) can be computed as NPP minus R h fluxes, excluding the effects of small-scale fires and other localized disturbances or vegetation regrowth patterns on carbon fluxes [26].The NASA-CASA soil model uses a set of compartmental difference equations.First-order decay equations simulate exchanges of decomposing plant residue (metabolic and structural fractions) at the soil surface.The model also simulates surface soil organic matter (SOM) fractions that presumably vary in age and chemical composition.Turnover of active (microbial biomass and labile substrates), slow (chemically protected), and passive (physically protected) fractions of the SOM are represented.
Evapotranspiration is connected to water content in the soil profile layers (Figure 1), as estimated using the NASA-CASA algorithms described by Potter [2].The soil model design includes three-layer (M1-M3) heat and moisture content computations: surface organic matter, topsoil (0.3 m), and subsoil to rooting depth (1 to 10 m).These layers can differ in soil texture, moisture holding capacity, and carbon-nitrogen dynamics.The setting for deep rooting depths (up to 10 meters) in tropical forest biomes follows the findings from studies of primary production seasonality in these regions [21][22][23].Water balance in the soil is modeled as the difference between precipitation or volumetric percolation inputs, monthly estimates of PET, and the drainage output for each layer.Inputs from rainfall can recharge the soil layers to field capacity.Excess water percolates through to lower layers

Deforestation Carbon Fluxes
The general method used in this study to compute biomass burning gas emissions is based on the approach described by Potter et al. [27].To estimate regional trace gas emissions from vegetation fires, we apply the following equation: where E t (Pg; 1 Pg = 10 15 g) is the regional emissions total at time t (d), B is the biomass subjected to burn at location x, C F is the biomass combustion fraction associated with a particular plant tissue fraction (leaf versus wood), e f is the emission factor (flaming and/or smoldering) associated with a particular trace gas, A is the area burned (km 2 ) at location x and time t.
To estimate the B term in Equation ( 2), maps of vegetation biomass can be derived by one of two general methods.The first is by spatial interpolation, using what is normally a small number (<100) of intensive field site measurements of aboveground plant mass [28].A weakness of any interpolation method is that a small number of measurements may not adequately represent the variability of biomass growth patterns.The second general method, and the one used in this study, is developed through our combination of satellite remote sensing and ecosystem carbon flux modeling.Satellite imagery can be transformed using plant production models to provide relatively high spatial resolution maps of above-ground biomass over a regional area of interest [2].
We adopted default C F values largely from the FLAMBE modeling system [29], with several modifications.We used a globally uniform C F value of 0.95 for leaf material [27].In tropical rainforest zones, we adopted a C F value of 0.45 for wood material, derived from Amazon forest slash burning studies [30,31].Following the approach of Potter et al. [27], we based these estimates for typical C F values on studies that were conducted in the Amazon on small-holder properties, where all decisions as to which vegetation to burn, size of area slashed, location, and the timing of the slash and burn process (how to slash, how long to dry, when to burn) were entirely left to property owners.The estimated C F values used in our analysis are typical of those reported in several other studies of tropical biomass burning [32].
The e f term in Equation ( 2) is defined as the amount of a compound released per amount of fuel consumed (g dry matter).Calculation of this parameter requires knowledge of the carbon content of the biomass burned and the carbon budget of the fire usually expressed as the C F term [33].Where fuel and residue data at the ground level are not available, an overall fuel carbon content of 45% -50% is commonly assumed [34].
The e f vales we used in Equation (2) were estimated by Scholes et al. [35] based on a review of some 70 publications, a large fraction of which were produced as a result of International Geosphere-Biosphere Programme (IGBP) Biomass Burning Experiment (BIBEX) campaigns.It appears from this compilation of published e f values that adequate data exist for CO 2 emissions for savanna and tropical forest fires.Post-disturbance decomposition of residual biomass carbon pools in wood and soils followed the methods of Potter et al. [7].

Ecosystem Disturbance Events
In this section, we describe algorithms implemented for mapping global land cover change and wildfires based on satellite observations from MODIS data by Mithal et al. [36].The new forest cover change algorithms were unsupervised in nature and exploit both the temporal and spatial structure in the MODIS data.Using independent wildfire perimeter data sets, we have comprehensively evaluated these algorithms, as well as those from alternate methods, across different forest climatic zones.The framework depends upon the Enhanced Vegetation Index (EVI) from MODIS 16-day Level 3 1 km Vegetation Indices (MOD13A2) products.
The following is a brief description of the algorithms used in this framework for mapping of global forest cover change and wildfires, and subsequently linked to CASA forest biomass pools for carbon emission fluxes.In this stratified framework for mapping forest disturbances, we have employed multiple, complementary scoring mechanisms using 1 km MODIS EVI time series products.The main assumption behind these algorithms is that in mature forests, EVI values for future time steps tend to be similar to previous years when accounting for seasonal variation.On the other hand, changes like wildfires and deforestation are characterized by an abrupt decrease in EVI after the event.The algorithms build a model that is used for predicting the expected EVI values for the future years.Deviation of the future observations from the predicted value indicates a change.A measure that quantities the deviation of future observations from the model prediction is used to assign the change score.
The change score should reflect the significance of deviations with respect to the natural variation of vegetation response for a given forest location.The VID (Vegetation-Independent Yearly Delta) algorithm developed by Mithal et al. [36] addressed this requirement by including the standard deviation of the variability in the change score.It assumed that the random fluctuations in mean annual EVI for a particular vegetation type are normally distributed for a location and estimates the mean (µvar) and standard deviation (∂var) of the variability score distribution as the maximum likelihood estimates for the distribution.
The new VID score was computed as Equation (3):

 
Yearly Delta score var var Mithal et al. [36] examined the performance of the new VID forest disturbance algorithm in several regions around the world, including the states of California (US) and the Yukon (Canada).Results demonstrated high accuracy levels for all major wildfires mapped from aerial surveys in these diverse forest regions.Mithal et al. [36] quantitatively showed that the VID forest change algorithms perform better than the well-known MODIS burned area (BA) framework [37] in the state of California (US) and were comparable in results to the BA framework for wildfire disturbances in the Yukon (Canada).

CASA Validation with Tower Flux Measurements
Flux estimates from eddy-correlation analysis were obtained from AmeriFlux tower flux sites that could meet certain criteria for CASA model comparisons.First, at least three complete years of site flux measurements were required to evaluate model predictions of interannual variations in CASA NEP fluxes.Second, winter (and/or dormant/dry) season NEP fluxes were required from a site to evaluate model predictions of soil CO 2 emissions on a year-round basis.Third, tower sites were required to be representative of the predominant vegetation class setting in the global land cover data used as input to the CASA model.
For sites meeting all of these criteria, AmeriFlux data sets were obtained from the central data repository located at the Carbon Dioxide Information Analysis Center (CDIAC; public.ornl.gov/ameriflux/dataproducts.shtml).Level 4 AmeriFlux records contained gap-filled and µstar filtered records, complete with calculated gross productivity and total ecosystem respiration terms on varying time intervals including hourly, daily, weekly, and monthly with flags for the quality of the original and gapfilled data.
CASA monthly NEP predictions from the MOD13C2 EVI data values closest to the tower location were compared to AmeriFlux eddy-correlation monthly estimates of the corresponding NEP fluxes.We note that the monthly MODIS EVI values in practically every grid cell of the global CASA model will be influenced by periodic land cover disturbances and (some naturally occurring) areas of sparse vegetation cover, including development, roads, water bodies.It was expected, therefore, that CASA model NEP flux predictions would be systematically lower than tower measurements of these carbon fluxes, since tower footprints tend to be far less affected by wildfire and other disturbances (such as logging and forest thinning), compared for instance to the surrounding MODIS grid cell area in which they are located.
A total of four AmeriFlux tower sites, together reporting 196 monthly measurements, were found to meet the criteria cited above for comparison to CASA model NEP predictions (Figure 2).CASA model predictions closely followed the seasonal timing of AmeriFlux tower measurements at each site.The linear regression correlation coefficient between CASA NEP predictions and tower fluxes was estimated at R 2 = 0.41 (p < 0.05) for all sites combined.Most of the unexplained variance in this model-to-tower flux comparison resulted from the Blodgett evergreen needle-leaf forest tower site, which showed a sudden shift upwards in NEP (greatly enhanced ecosystem carbon sink) after 2001 in the AmeriFlux data sets, presumably due to stand thinning in 2000 [38].We note that Tang et al. [39] reported a similar correlation value (R 2 = 0.45) in model-to-tower NEP flux comparesons at the Blodgett forest location, compared to higher correlations (all with R 2 > 0.50) at seven other AmeriFlux forest sites, which implied that small-scale forest management at the Blodgett site is an unaccounted source of uncertainty in our model flux comparisons.

Global Net Primary Production
As previously reported by Potter et al. [13], predicted terrestrial NPP for the globe in 2009 was 50.1 Pg C, a total carbon flux in the middle of the range of previous vegetation NPP predictions of between 44 to 66 Pg C per year for the period 1982-1998 [40,41].We estimate that global terrestrial NPP increased by +0.14 Pg C over the time period of 2000 to 2009, due almost entirely to a strong upward trend in the Northern Hemisphere (Figure 3).Annual NPP was predicted to have increased between the years 2000 and 2007 in the regions of high-latitude (>50˚N) North America and Eurasia, and also in South Asia, West and Central Africa, and the western Amazon.This upward trend in high-latitude NPP was controlled by a combination of rapidly warming temperatures from 2004 to 2005 and by elevated MODIS EVI patterns, which in turn were closely correlated with precipitation amounts [13].Periodic declines in regional NPP levels were predicted for the southern Untied States, the southern Amazon, western Europe, southern and eastern Africa, and Australia; the timing of negative NPP anomalies in each of these regions was associated with severe droughts and, in some cases, extreme heat waves [42].was significantly associated with positive NEP (net sink) fluxes in high northern latitude tundra, grasslands, and boreal forest areas, whereas from 2005-2009, major drought events were associated with negative NEP (net source fluxes) in tropical evergreen forests, temperate deciduous forests, croplands, grasslands, and savannas worldwide (Figure 5).

Ecosystem Disturbance Emissions
Comparison of areas of disturbed forest land between CASA model inputs and national reports from the Food and Agriculture Organization [43] of the United Nations showed a close match among the 30 leading counties in terms of hectares of forest converted (Table 1).Six of the top ten counties ranked by the FAO in terms of annual forest conversion rates (2005)(2006)(2007)(2008)(2009) were also among the top ten counties mapped for forest lands converted in our CASA model inputs from 2005-2009 MODIS data.The notable exceptions in the global comparison shown in Table 1 were that of India and Papua New Guinea (FAO ranks 3 and 6, respectively).These two Asian countries were still ranked among the 30 leading counties for forest area converted in CASA model inputs.According to the FAO [43], India recorded over 25 million hectares of forests as being affected by grazing by domestic animals, a cover change that involves the type of gradual forest degradation processes that are difficult to detect by satellite remote sensing from MODIS.
CASA model predictions of CO 2 emissions from forest disturbance and biomass burning (Figure 6) were highest on a unit area (g C m −2 yr −1 ) basis in the region of Southeast Asia (specifically in the countries of Myanmar, Malaysia, Cambodia, Vietnam, and Indonesia).Although the CO 2 emissions on a unit area basis from forest disturbance in the United States and Brazil were estimated at less than half of those estimated from Southeast Asian countries (Table 2), the total areas of forest disturbed annually gave the United States and Brazil the highest national totals of carbon lost from biomass burning in 2009.Forested regions of the Pacific Northwest, the southeastern US and Gulf Coast, and the Amazon rainforest were consequently major contributors to global biomass emissions to the atmosphere (Figure 6).The global sum of CO 2 emissions from forest disturbance and biomass burning for 2009 alone was predicted at 0.51 Pg C yr −1 .Decomposition emissions of residual (dead) forest biomass in the CASA model from three years (2007 to 2009) of deforestation globally added 0.15 Pg C yr −1 to the atmosphere.

Discussion
A recent study by Zhao and Running [19] reported a decreasing trend in global terrestrial NPP from 2000 to 2009, using MODIS satellite inputs and the same NCEP reanalysis data set as in the CASA model for climate inputs.On a regional basis, our CASA model results differed from Zhao and Running [19] which reported that NPP in the tropical zones (23.5˚S to 23.5˚N) explained 93% of variations in the global NPP.In contrast, we found that NPP in the tropical zones explained only 50% -60% of variations in the global NPP, whereas NPP in the latitude zone between 30˚N and 60˚N could explain between 40% and 50% of variations in the global NPP [13].Notwithstanding the difference in the global trend of NPP between CASA and the predictions from Zhao and Running [19]  The more complete CASA model NEP results reported in this paper suggest that surface temperature warming together with regional droughts (e.g., from 2006 to 2009) can drive ecosystem carbon losses to the atmosphere of more than 1 Pg C yr −1 in excess of long-term average terrestrial NEP fluxes.These predicted changes in NEP fluxes over a decade of model results would have made a larger impact on atmospheric CO 2 concentrations than NPP trends alone, and our results highlight the importance of including the annual variations in soil R h decomposition fluxes from downed and burned plant biomass in global carbon cycle.Variations in plant production alone can account for less than one-third of terrestrial ecosystem fluxes in most years.
The forest fire and land cover change mapping framework presented in this paper has limitations under several scenarios.These include situations where 1) the vegetation rapidly recovers after a fire or if there are multiple fires in rapid succession; 2) the loss in vegetation green cover associated with land cover conversion is insignificant, such as in crop fallow and rotation practices; 3) the vegetation cover has high natural variability in seasonal greenness, which is common in mixed woodland-grassland ecosystems.Each of these scenarios poses distinct challenges for our current land cover change detection framework that are being addressed in future validation studies with extensive ground-truth data sets.
Nevertheless, we have identified numerous relatively small-scale patterns throughout the world where terrestrial carbon fluxes may vary between net annual sources and sinks from one year to the next.We conclude that accurately monitoring of NEP for these areas of high interannual variability will require further validation of carbon model estimates, with a focus on both flux algorithm mechanisms and potential scaling errors to the regional level.

Figure 1 .
Figure 1.Schematic representation of components in the NASA-CASA model.The soil profile component (a) is layered with depth into a surface ponded layer (M 0 ), a surface organic layer (M 1 ), a surface organic-mineral layer (M 2 ), and a subsurface mineral layer (M 3 ), showing typical levels of soil water content (shaded) in three general vegetation types.The production and decomposition component (b) shows separate pools for carbon cycling among pools of leaf litter, root litter, woody detrius, microbes, and soil organic matter, with dependence on litter quality (Lit q).t Subtracting monthly R h fluxes from monthly NPP fluxes yields NEP flux estimates.Predicted global NEP fluxes from 2000 through 2005 showed a roughly balanced terrestrial biosphere carbon cycle, with variations less than ±0.5 Pg C yr −1 (Figure 4).Nonetheless, beginning in 2006, global NEP fluxes became increasingly imbalanced, starting from −0.9 Pg C yr −1 to the largest negative (net terrestrial source) flux of −2.2 Pg C yr −1 in 2009.Notable surface temperature warming from 2000-2005

Figure 4 .
Figure 4. Interannual variations from 2000 through 2009 in anomalies of annual total NEP for the CASA model for the Northern Hemisphere (NH-green circles) and the Southern Hemisphere (SH-red circles).Units are Pg C yr −1 (1 Pg = 10 15 g), with positive values indicating net ecosystem sink fluxes and negative values indicating net ecosystem source fluxes.

Figure 3 .
Figure 3. Spatial pattern of terrestrial NPP linear trends from 2000 through 2009.

Figure 5 .
Figure 5. Global predicted NEP fluxes from (a) 2003 (highest total ecosystem sink flux) and (b) 2009 (highest total ecosystem source flux).Units are g C m −2 yr −1 gridded at 0.5 deg spatial resolution.Europe, Canada, Russia, and China experienced summer periods of warming temperatures.From 2004 to 2008, parts of Pakistan, Australia, the United States, Canada, and Europe continued to experience extreme heat waves.On the other hand, extreme cold winter temperatures were reported repeatedly in Russia, India, and China over the period 2000 to 2008.Large parts of Europe experienced unusually cold temperatures in the summer of 2001, as did Australia and South Africa in 2007.There were many instances of severe drought across the globe during the period of 2000 to 2008, mainly affecting regions of the central North America, Africa, Brazil, and China (World Meteorological Organization, [42]).Beginning with major droughts in Brazil, the Horn of Africa, the Middle East, Central and South Asia, and China in 2000 and 2001, these events were followed by most of North America, southern Africa, and Australia experiencing record low precipitation amounts in 2002, 2003, and 2004.Large areas of Europe, southern Africa, Brazil, and Paraguay were affected by severe droughts in 2005.From 2006 though 2008, much of the United States, eastern and southern Africa, China, and Australia experienced continued deficits of precipitation.Strongly negative NEP fluxes were predicted to be associated with droughts reported in South Asia, eastern Africa, northern China, and northern and eastern coastal South America.Strongly positive NEP fluxes predicted by the CASA model were associated with periodically heavy rainfall

Figure 6 .
Figure 6.Global predicted biomass burning fluxes of CO 2 in 2009.Units are g C m −2 yr −1 gridded at 0.5 degree spatial resolution.
the overall patterns of interannual variations in Northern and Southern Hemisphere NPP anomalies were similar between the two model results.NPP anomalies in the Northern Hemisphere were nega-tive from 2000-2003 and then became strongly positive from 2004-2008, closely following the 0.1˚ yr −1 surfacewarming trend in the model input data.NPP anomalies in the Southern Hemisphere were positive from 2000-2003 and then turned negative between 2004-2008, with 2005 being the most strongly negative anomaly year.