Abstract

Remote sensing can provide long-term and large-scale products helpful for ecosystem model evaluation. The authors compare monthly gross primary production (GPP) simulated by the Community Land Model, version 4 (CLM4) at a half-degree resolution with satellite estimates of GPP from the Moderate Resolution Imaging Spectroradiometer (MODIS) GPP product (MOD17) for the 10-yr period January 2000–December 2009. The assessment is presented in terms of long-term mean carbon assimilation, seasonal mean distributions, amplitude and phase of the annual cycle, and intraannual and interannual GPP variability and their responses to climate variables. For the long-term annual and seasonal means, major GPP patterns are clearly demonstrated by both products. Compared to the MODIS product, CLM4 overestimates the magnitude of GPP for tropical evergreen forests. CLM4 has a longer carbon uptake period than MODIS for most plant functional types (PFTs) with an earlier onset of GPP in spring and a later decline of GPP in autumn. Empirical orthogonal function analysis of the monthly GPP changes indicates that, on the intraannual scale, both CLM4 and MODIS display similar spatial representations and temporal patterns for most terrestrial ecosystems except in northeast Russia and in the very dry region of central Australia. For 2000–09, CLM4 simulated increases in annual averaged GPP over both hemispheres; however, estimates from MODIS suggest a reduction in the Southern Hemisphere (−0.2173 PgC yr−1), balancing the significant increase over the Northern Hemisphere (0.2157 PgC yr−1). The evaluations highlight strengths and weaknesses of the CLM4 primary production and illuminate potential improvements and developments.

1. Introduction

To understand historical changes and predict future trends in ecosystems we must improve our understanding of individual ecosystem processes and their interactions with external environmental factors (Fung et al. 2005; Friedlingstein et al. 2006; Thornton et al. 2009). Gross primary production (GPP) is the amount of carbon assimilated via photosynthesis and constitutes an important link in the terrestrial carbon cycle (Ciais et al. 1997). Because of its importance to human society and welfare, observations and simulations of GPP have received considerable attention from both academic communities and government agencies (Solomon et al. 2007; Beer et al. 2010; Bonan et al. 2011).

We can enhance our ability to interpret the advantages and disadvantages of GPP among site level observation, process-based model simulation, statistical calculation, and remote sensing estimates through intercomparison of each method (Randerson et al. 2009; Beer et al. 2010). Mechanistic ecosystem models are generally based on observed ecophysiological, biophysical, and biogeochemical processes, and constrained by changes in environmental conditions (Cramer et al. 1999; Sitch et al. 2008). The online or offline simulation of such computer models at different temporal and spatial scales, therefore, is essential for both research and management to predict changes in GPP with varying climate (Cox et al. 2000; Jones et al. 2009). Previous studies have shown that there is substantial uncertainty in process-based models regarding model structure, parameterization, and driver data (Zaehle et al. 2005; Jung et al. 2007; Baker et al. 2010; Wang et al. 2011; Bonan et al. 2011). Thus, comprehensive evaluations using multiple data sources are essential for ecosystem model improvement.

Evaluating GPP from global ecosystem models is hindered by the lack of extensive observations at continental and global scales. Statistical methods have been used to upscale site-derived GPP to regions or the globe (Xiao et al. 2008; Jung et al. 2009); however, such diagnostic models lack the capacity to explicitly simulate the future behavior of ecosystems under changing environments. Moreover, they rely largely on the quality of eddy covariance flux measurement data, global fraction of absorbed photosynthetically active radiation (FPAR) retrievals, and the statistical techniques employed (Beer et al. 2010).

Remote sensing techniques provide near-real-time metrics related to global vegetation growth and primary production, such as the normalized difference vegetation index (NDVI), at fine spatiotemporal scales (Justice et al. 2002; Nemani et al. 2003; Running et al. 2004). The Moderate Resolution Imaging Spectroradiometer (MODIS) GPP product (MOD17) represents the latest of a series of efforts to characterize global GPP from space (Zhao and Running 2010, 2011). MODIS GPP products have been extensively evaluated (Running et al. 1999; Heinsch et al. 2006; Turner et al. 2006a,b; Kanniah et al. 2009). Although a wide range of site-specific fidelity between the MODIS GPP and flux tower estimates has been found, the satellite-derived GPP generally provides realistic information on widely varying climates, land use, and vegetation physiognomy over the globe, and can provide valuable global-scale constraints to process models (Zhao et al. 2005; Turner et al. 2006a,b).

With the remote sensing GPP, we have evaluated predictions of global GPP from the Community Land Model version 4 (CLM4), the land component of the Community Earth System Model (CESM) (Oleson et al. 2010; Lawrence et al. 2011). Providing the correct simulation and prediction of terrestrial GPP for CLM4 is an essential step in reducing uncertainty of the future trends of carbon–climate feedbacks. Although many studies have used satellite products to evaluate earlier versions of CLM (Oleson et al. 2003; Tian et al. 2004; Wang et al. 2004; Lawrence and Chase 2007; Stöckli et al. 2008), so far no evaluation of model results against remote-sensing-based GPP has been performed.

We have compared CLM4 and MODIS-based vegetation primary production at various temporal and spatial scales for the period 2000–09. Our objectives are to 1) systematically evaluate the CLM4 capability to represent contemporary global patterns of GPP; 2) quantify the similarities and the differences between the two products; 3) investigate GPP responses to changes in climatic forcings such as temperature, precipitation, and radiation; and 4) explain differences in and address the future development of biogeochemical dynamics in CLM4. The data and methods are described in section 2. In section 3 we show the statistical comparisons between the two estimates. Specific findings and interpretations are discussed in section 4, with conclusions given in section 5.

2. Methodology

a. CLM4

CLM4 represents fundamental biogeochemical and biophysical mechanisms of terrestrial ecosystems. It is a coupled model from the biophysics of CLM and the carbon–nitrogen biogeochemistry of the Biome-BGC4.1.2 (Thornton and Rosenbloom 2005). The resulting model is fully prognostic with respect to all water, energy, carbon, and nitrogen variables in the terrestrial ecosystem. Gross primary production in CLM4 is the sum of carbon assimilation from both the sunlit and shaded portions of the canopy. It is affected by temperature through enzyme activity and stomatal conductance, controlled by soil moisture and root distribution, and limited by the availability of soil and leaf nitrogen. A detailed description of the biophysical structures and biogeochemical components that influence the CLM4 primary production is given in Oleson et al. (2010).

b. MODIS GPP

We used the improved MOD17 gross primary production (GPP) (Zhao et al. 2006; Zhao and Running 2010, 2011) and the collection 5 (C5) leaf area index (LAI) (MOD15A2) (Myneni et al. 2002) for years 2000–09 for our comparisons. The C5 MODIS LAI has been cleaned by temporally filling the cloud-contaminated periods based on accompanied quality assessment fields (Zhao et al. 2005). To compare with CLM4 model results, the MODIS monthly GPP at 0.5° resolution was aggregated from the raw 1-km resolution product and calculated daily by the light use efficiency algorithm:

 
formula

where is the maximum light use efficiency; FSDS is surface downward solar radiation, of which 45% is photosynthetically active radiation (PAR); and FPAR is the fraction of PAR being absorbed by the plants. Along with LAI data, it was remotely sensed from the MODIS sensor (Myneni et al. 2002); and are the daily reduction scalar from water stresses and low temperature, respectively. The long-term MODIS GPP/NPP dataset was obtained from the Numerical Terradynamic Simulation Group, School of Forestry, University of Montana (http://www.ntsg.umt.edu/).

c. Meteorological data

Surface meteorological data are critical drivers for both ecosystem model simulation and remotely sensed GPP estimates. The MODIS GPP algorithm used the daily National Centers for Environmental Prediction/Department of Energy Global Reanalysis 2 (NCEP-2) data, including downward solar radiation and minimum temperature as meteorological drivers (Kanamitsu et al. 2002). The half-degree Climate Research Unit–NCEP (CRUNCEP) dataset was applied to force the CLM4 model. It is a combination of two existing datasets: the CRU Time Series, version 2.1 (CRU TS 2.1) (Mitchell and Jones 2005) 0.5° monthly climatology covering the period 1901–2002 and the 2.5° NCEP-2 reanalysis data beginning in 1948 and available in near-real time (more details at http://dods.extra.cea.fr/data/p529viov/cruncep/).

d. Dominant PFTs distribution

GPP comparisons in this study are made at three spatial scales: global, hemispheric, and aggregated by dominant PFT. CLM4 represents a nested subgrid hierarchy in which grid cells comprise a land unit, snow/soil column, and up to 16 plant functional types (PFTs) (Oleson et al. 2010). For the MODIS GPP algorithm, land cover is from the University of Maryland Land Cover Classification System, based on the collection 4 MODIS 1-km land cover (MOD12Q1) (land cover type 2 in the MOD12Q1 dataset). Land cover in every half-degree grid cell is the dominant biome type based on the 1-km collection 4 MODIS products (Fig. 1a). Since CLM4 and MODIS adopt different vegetation classification systems, the dominant PFT distributions from CLM4 were used to analyze both GPP maps. The CLM4 PFT with more than 50% fractional area in a one-half-degree grid point is considered as a dominant PFT. Dominant PFTs, ignoring barren land cover as defined by MOD12Q1, are shown in Fig. 1b. Globally, broadleaf evergreen tropical trees have the highest coverage, while broadleaf deciduous boreal shrubs and needleleaf evergreen boreal trees are the second and third largest areas.

Fig. 1.

(a) Global vegetation map of the University of Maryland Land Cover Classification System. Water body (Water), evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), deciduous needleleaf forest (DNF), deciduous broadleaf forest (DBF), mixed forests (MF), closed shrublands (CShrub), open shrublands (OShrub), woody savannas (WSavanna), savannas (Savanna), grassland (grass), croplands (Crop), urban (Urb), and bare ground (Bare). Note that in this figure, there are no vegetation classifications for types 11, 14, and 15 (Noveg). (b) Dominant plant functional types of CLM4. Area without vegetation (Noveg), needleleaf evergreen temperate tree (NETem Tree), needleleaf evergreen boreal tree (NEBor Tree), needleleaf deciduous boreal tree (NDBor Tree), broadleaf evergreen tropical tree (BETro Tree), broadleaf evergreen temperate tree (BETem Tree), broadleaf deciduous tropical tree (BDTro Tree), broadleaf deciduous temperate tree (BDTem Tree), broadleaf deciduous boreal tree (BDBor Tree), broadleaf evergreen shrub (BE Shrub), broadleaf deciduous temperate shrub (BDTem Shrub), broadleaf deciduous boreal shrub (BDBor Shrub), C3 arctic grass (C3A Grass), C3 nonarctic grass (C3NA Grass), C4 grass, corn, and wheat.

Fig. 1.

(a) Global vegetation map of the University of Maryland Land Cover Classification System. Water body (Water), evergreen needleleaf forest (ENF), evergreen broadleaf forest (EBF), deciduous needleleaf forest (DNF), deciduous broadleaf forest (DBF), mixed forests (MF), closed shrublands (CShrub), open shrublands (OShrub), woody savannas (WSavanna), savannas (Savanna), grassland (grass), croplands (Crop), urban (Urb), and bare ground (Bare). Note that in this figure, there are no vegetation classifications for types 11, 14, and 15 (Noveg). (b) Dominant plant functional types of CLM4. Area without vegetation (Noveg), needleleaf evergreen temperate tree (NETem Tree), needleleaf evergreen boreal tree (NEBor Tree), needleleaf deciduous boreal tree (NDBor Tree), broadleaf evergreen tropical tree (BETro Tree), broadleaf evergreen temperate tree (BETem Tree), broadleaf deciduous tropical tree (BDTro Tree), broadleaf deciduous temperate tree (BDTem Tree), broadleaf deciduous boreal tree (BDBor Tree), broadleaf evergreen shrub (BE Shrub), broadleaf deciduous temperate shrub (BDTem Shrub), broadleaf deciduous boreal shrub (BDBor Shrub), C3 arctic grass (C3A Grass), C3 nonarctic grass (C3NA Grass), C4 grass, corn, and wheat.

e. Experimental design

Surface meteorological data from the CRUNCEP dataset for the period 1901–2009 was applied to drive the CLM4 with fully coupled carbon and nitrogen cycles in offline mode. Other environmental forcings such as historical land use and land cover, atmospheric CO2, and atmospheric nitrogen deposition were prescribed to follow historical trends as in Bonan and Levis (2010), and Shi et al. (2011). Based on a subset (1901–20) of transient climates, CLM4 was spun up using atmospheric CO2, nitrogen deposition, and land cover on year 1850. It was then integrated to 1901, driven by repeating the 20-yr meteorological variables with the transient CO2 concentration, nitrogen deposition, and historical land use data between 1850 and 1901. Beginning with the model state in 1901, CLM4 was finally run to 2009 with the previously mentioned historical forcings. The half-degree monthly GPP output during the last 10 years (2000–09) was selected for direct comparison with MODIS GPP. We compared the CLM4 GPP against several aspects of the satellite data, including the climatological means, and the intraannual and interannual variabilities. For convenience we refer to the CLM4 GPP and LAI as GPP1 and LAI1, respectively, while MODIS GPP, LAI, and FPAR are referred to as GPP2, LAI2, and FPAR2, respectively. CRUNCEP climate drivers (input to CLM4) are referenced as surface temperature (TEMP1), precipitation (PREC1), and surface downward solar radiation (FSDS1).

3. Results

a. Spatial distribution of GPP climatology

1) GPP magnitude

Both CLM4 and MODIS annual mean GPP show high spatial heterogeneity (Figs. 2a,b). Both show the highest carbon assimilation in tropical ecosystems, followed by temperate and boreal forests. High-latitude regions with short growing seasons and dry areas have the lowest GPP. Globally, mean total GPP from CLM4 is 146.30 PgC yr−1 and from MODIS it is 111.58 PgC yr−1, ignoring barren land cover defined using MOD12Q1 (Table 1). The spatial correlation coefficient between the two maps is 0.90, indicating that the CLM4 and the MODIS are in good agreement in describing the large-scale GPP distribution. CLM4 GPP is consistently higher than MODIS GPP across the tropics and subtropics, and lower than MODIS GPP at high northern latitudes and in southern Africa (Fig. 2c). Tropical GPP is higher for CLM4 than for MODIS throughout the year (Figs. 3a–d). In Arctic tundra regions the low bias in CLM4 GPP relative to MODIS is most severe during the summer (Fig. 3b).

Fig. 2.

Annual mean (ANN) GPP (gC m−2 month−1) between 2000 and 2009 for (a) CLM4, (b) MODIS, and (c) the difference (gC m−2 month−1) between CLM4 and MODIS.

Fig. 2.

Annual mean (ANN) GPP (gC m−2 month−1) between 2000 and 2009 for (a) CLM4, (b) MODIS, and (c) the difference (gC m−2 month−1) between CLM4 and MODIS.

Table 1.

Mean values, standard deviations, and trends of spatially averaged annual GPP and climate variables for the globe, the Northern Hemisphere, and the Southern Hemisphere from 2000 to 2009. Bold values represent trends with significance (P < 0.05) during the study period.

Mean values, standard deviations, and trends of spatially averaged annual GPP and climate variables for the globe, the Northern Hemisphere, and the Southern Hemisphere from 2000 to 2009. Bold values represent trends with significance (P < 0.05) during the study period.
Mean values, standard deviations, and trends of spatially averaged annual GPP and climate variables for the globe, the Northern Hemisphere, and the Southern Hemisphere from 2000 to 2009. Bold values represent trends with significance (P < 0.05) during the study period.
Fig. 3.

The spatial difference pattern of the climatological mean GPP (gC m−2 month−1) during 2000–09 between CLM4 and MODIS in (a) March–May (MAM), (b) June–August (JJA), (c) September–November (SON), and (d) December–February (DJF).

Fig. 3.

The spatial difference pattern of the climatological mean GPP (gC m−2 month−1) during 2000–09 between CLM4 and MODIS in (a) March–May (MAM), (b) June–August (JJA), (c) September–November (SON), and (d) December–February (DJF).

2) GPP phase

To evaluate differences in seasonal timing of GPP variation independent of biases in GPP magnitude, we normalized the long-term mean seasonal cycles of CLM4 and MODIS GPP by their climatological maximum values in each grid cell and compared the seasonal differences in these normalized values between CLM4 and MODIS (Fig. 4). Normalized GPP from CLM4 is higher than MODIS across the boreal forest zone in spring (Fig. 4a) and across much of the Northern Hemisphere midlatitudes in fall (Fig. 4c). Seasonality of GPP is also notably different (sometimes higher, sometimes lower) between CLM4 and MODIS in some dry regions such as interior Australia, interior western North America, and the Kirghiz Steppe north of the Caspian and Aral Seas. These dry regions are classified as grasslands or open shrub lands in the MODIS products, while the dominant cover type is bare ground in the CLM4 subgrid classification (Fig. 1).

Fig. 4.

The spatial difference pattern of the normalized climatological mean GPP during 2000–09 between CLM4 and MODIS in (a) MAM, (b) JJA, (c) SON, and (d) DJF.

Fig. 4.

The spatial difference pattern of the normalized climatological mean GPP during 2000–09 between CLM4 and MODIS in (a) MAM, (b) JJA, (c) SON, and (d) DJF.

b. Seasonal cycle of GPP summarized by vegetation zones

Both the half-degree CLM4 and MODIS GPP are derived from subgrid-scale GPP based on different single-type definitions. CLM4 includes multiple plant functional types at the subgrid scale, while the half-degree MODIS product is aggregated from the raw dataset at 1-km resolution. To minimize ambiguity in our comparisons, we summarized seasonality of GPP and LAI from both CLM4 and MODIS over regions dominated by individual plant function types in the CLM4 classification (Figs. 1, 5, and 6).

Fig. 5.

The 2000–09 climatological annual cycle of CLM4 (red line) and MODIS (black line) GPP (gC m−2 month−1) on selected dominant PFTs defined in Fig. 1b: (a)–(k) over the globe, (l),(n) the Southern Hemisphere, and (m),(o) the Northern Hemisphere.

Fig. 5.

The 2000–09 climatological annual cycle of CLM4 (red line) and MODIS (black line) GPP (gC m−2 month−1) on selected dominant PFTs defined in Fig. 1b: (a)–(k) over the globe, (l),(n) the Southern Hemisphere, and (m),(o) the Northern Hemisphere.

Fig. 6.

The 2000–09 normalized climatological annual cycle of GPP (solid line) and LAI (dashed line) for CLM4 (red line) and MODIS (black line) on selected dominant PFTs defined in Fig. 1b: (a)–(k) over the globe, (l),(n) the Southern Hemisphere, and (m),(o) the Northern Hemisphere.

Fig. 6.

The 2000–09 normalized climatological annual cycle of GPP (solid line) and LAI (dashed line) for CLM4 (red line) and MODIS (black line) on selected dominant PFTs defined in Fig. 1b: (a)–(k) over the globe, (l),(n) the Southern Hemisphere, and (m),(o) the Northern Hemisphere.

In the boreal zone, CLM4 evergreen GPP has a midsummer peak similar to MODIS, with an increase earlier in the spring and a later decline in the fall (Fig. 5b). For boreal deciduous needleleaf trees, CLM4 has a lower midsummer GPP than MODIS (Fig. 5c). For boreal deciduous broadleaf trees, however, CLM4 GPP has a higher midsummer peak than MODIS and a later decline in fall (Fig. 5h). Similar differences are seen in the timing of CLM4 and MODIS LAI for these types (Fig. 6). CLM4 has relatively little seasonal variation in boreal evergreen LAI, while MODIS shows a sharp decline in LAI over the winter (Fig. 6b). CLM4 LAI for boreal deciduous trees is maintained at midseason peak values for an extra month in the fall compared to MODIS (Figs. 6c,h).

In northern temperate regions the timing of spring and fall changes in CLM4 GPP are in good agreement with MODIS for both evergreen and deciduous types, but CLM4 has an overall higher GPP for evergreens (Figs. 5a,g; 6a,g). CLM4 temperate evergreen forest LAI is relatively constant over the year, while MODIS shows variation of approximately a factor of 2 from winter to summer (Fig. 6a). Similar to the boreal zone, CLM4 temperate deciduous forest LAI is maintained at peak values for about a month longer in the fall than MODIS LAI (Fig. 6g).

Neither CLM4 nor MODIS shows significant seasonality for either GPP or LAI in tropical deciduous or evergreen forests (Figs. 6d,f). CLM4 GPP is higher than MODIS GPP in all months for tropical evergreen forests (Fig. 5d), but CLM4 and MODIS GPP are of similar magnitude in tropical deciduous forests (Figs. 5f).

Arctic shrub GPP is lower in CLM4 than in MODIS (Fig. 5j) and, while the seasonal timing of GPP between the two products is similar, CLM4 LAI remains near its midsummer peak almost two months longer in the fall than with MODIS LAI (Fig. 6j). Arctic grass has a similar peak GPP in CLM4 and MODIS, but here again the CLM4 LAI remains high for about two months longer than with the MODIS LAI (Figs. 5k and 6k), corresponding to higher autumn GPP in CLM4 compared to MODIS for this vegetation type (Fig. 5k).

Because C3 and C4 grasses have significant distribution in both Northern and Southern Hemispheres, with distinctive seasonality, we analyzed the two hemispheres separately for these types. In both hemispheres the seasonality and magnitude of C3 grass GPP are similar for CLM4 and MODIS (Figs. 5l,m). For C4 grass, CLM4 and MODIS GPP have similar timing in both hemispheres, but the Northern Hemisphere CLM4 has a higher magnitude throughout the year (Figs. 5n,o).

c. Intraannual variability

We used an empirical orthogonal function (EOF) analysis to characterize the dominant modes of monthly variability in both CLM4 and MODIS GPP (Fig. 7). Following North et al. (1982), we calculated the percentage variance accounted for by the first four modes of the two monthly GPPs (Fig. 7a). During each month, the first and second modes of the MODIS and CLM4 GPP are well distinguished from the rest of the monthly EOFs in terms of the one standard deviation (SD) of the sampling errors. The first two eigenfunctions (EOF1 and EOF2) explain, respectively, 48.5% and 12.7% of the monthly GPP variability in CLM4 and 51.8% and 16.2% of the monthly GPP variability in MODIS. The spatial and temporal patterns of EOF1 are similar for CLM4 and MODIS, and capture the annual cycle of peak GPP as it shifts between the Northern and Southern Hemispheres under the fundamental influence of incident radiation and its effects on temperature. Some exceptions to this overall pattern are seen in central Australia and the northern Amazon basin. Spatial and temporal patterns of EOF2 also are remarkably similar for CLM4 and MODIS GPP, capturing a weaker seasonality in the tropics and are also apparently influenced by monsoon regions in Asia, Australia, and southwestern North America. A notable exception for this EOF is over northeast Russia where CLM4 and MODIS have opposite signs.

Fig. 7.

The 2000–09 variance (%) explained by the (a) first four EOF modes of monthly GPP of CLM4 and MODIS, spatial patterns of the first two EOFs from the monthly GPP of (b),(c) CLM4 and (d),(e) MODIS, and time evolutions of the first two EOFs from the monthly GPP of (f),(g) CLM4 and (h),(i) MODIS.

Fig. 7.

The 2000–09 variance (%) explained by the (a) first four EOF modes of monthly GPP of CLM4 and MODIS, spatial patterns of the first two EOFs from the monthly GPP of (b),(c) CLM4 and (d),(e) MODIS, and time evolutions of the first two EOFs from the monthly GPP of (f),(g) CLM4 and (h),(i) MODIS.

d. Interannual variation

Interannual variability in global GPP is higher for CLM4 than for MODIS, with interannual standard deviations (SDs) of 3.09 and 0.85 PgC yr−1, respectively, for the period 2000–09 (Table 1 and Fig. 8). The spatial distributions of the amplitude of interannual variation of GPP for CLM4 and MODIS are shown in Fig. S1. Globally, main features of the 10-yr standard deviations of the two GPP products are generally similar, with larger interannual variability over the low and midlatitudes, such as tropical Africa and southeastern America, than in high latitudes, such as the tundra of northeastern Russia (Figs. S1a,b). The amplitude of the interannual variation of CLM4 GPP is stronger than the satellite product over most ecosystems, especially in the tropical trees, C4 grasses, and low and midlatitude crops (Fig. S1c and Fig. 1). Major underestimations exist in the tropical Amazon, in southern Africa, and in the boreal shrub areas.

Fig. 8.

Spatially averaged time series of standardized annual CLM4 and MODIS GPP and LAI, MODIS FPAR, CRUNCEP temperature, precipitation, and surface downward solar radiation over (a) the globe, (b) the Northern Hemisphere and (c) the Southern Hemisphere from 2000 to 2009.

Fig. 8.

Spatially averaged time series of standardized annual CLM4 and MODIS GPP and LAI, MODIS FPAR, CRUNCEP temperature, precipitation, and surface downward solar radiation over (a) the globe, (b) the Northern Hemisphere and (c) the Southern Hemisphere from 2000 to 2009.

At the global scale, CLM4 GPP shows a significant positive trend over this period (0.8260 PgC yr−1), while MODIS GPP shows no significant trend (Table 1). Both CLM4 and MODIS GPP show significant positive trends over the Northern Hemisphere but not over the Southern Hemisphere (Table 1). Both CLM4 and MODIS GPP estimates are anomalously low in 2005 over the Southern Hemisphere, associated with low precipitation and high temperatures (Fig. 8c). Pearson correlation analysis was used to explore relationships among GPP, LAI, FPAR, and climate forcings at global and hemispheric scales based on interannual variability (Table 2). CLM4 GPP has a significant positive correlation with CLM4 LAI on a global basis and in each hemisphere, while MODIS GPP is not significantly correlated with MODIS LAI. MODIS GPP has a significant positive correlation with its FPAR in the Southern Hemisphere. CLM4 GPP and LAI show no significant correlation with temperature, while MODIS GPP has a significant positive correlation with temperature in the Northern Hemisphere and a significant negative correlation in the Southern Hemisphere. CLM4 GPP and LAI show significant positive correlations with PREC1 for the globe and in both hemispheres, while MODIS GPP shows a significant positive correlation with PREC1 for only the Northern Hemisphere. Both CLM4 GPP and LAI show significant negative correlation with incident shortwave radiation for the globe and in both hemispheres, while MODIS GPP has a significant negative correlation with FSDS1 in the Northern Hemisphere only.

Table 2.

Pearson correlations among the spatially averaged annual GPP, LAI, FPAR, and climate variables for the globe and two hemispheres from 2000 to 2009. Bold values mean correlations with significance (P < 0.05) during the 10-yr period.

Pearson correlations among the spatially averaged annual GPP, LAI, FPAR, and climate variables for the globe and two hemispheres from 2000 to 2009. Bold values mean correlations with significance (P < 0.05) during the 10-yr period.
Pearson correlations among the spatially averaged annual GPP, LAI, FPAR, and climate variables for the globe and two hemispheres from 2000 to 2009. Bold values mean correlations with significance (P < 0.05) during the 10-yr period.

4. Discussion

a. Uncertainties of the remote sensing GPP

Our analysis examined differences in the GPP using two different products at different temporal and spatial scales. The comparisons, however, are complicated since there are a number of possible reasons for the differences in CLM4 and MODIS GPP estimates due to different parameter values, calculation algorithms, and environmental inputs. Strictly speaking, the MODIS GPP is also a modeled product. The performance of the MODIS GPP can be largely influenced by uncertainties from inputs such as land cover, FPAR/LAI, reanalyzed meteorological observations, and the algorithm itself (Zhao et al. 2006). Also, the scaling of high resolution (1 km) raw data to the coarse (0.5°) product can introduce considerable uncertainties of the remote sensing GPP. As a result, it may contain systematic errors in some regions. More detailed discussions of these aspects can be found elsewhere (Coops et al. 2009; Heinsch et al. 2006; Zhao et al. 2005).

Nonetheless, among a variety of measurements and simulations, the MODIS satellite remains a unique tool for monitoring terrestrial primary productivity and is a benchmark for gauging ecosystem model improvement because of its globally high spatiotemporal resolution and continuity. More importantly, at the global scale the MODIS GPP estimates are likely very reasonable and some general observations may be warranted (Zhao and Running 2010, 2011; Frankenberg et al. 2011). Hence, in the following we mainly focus on discussions of possible causes from the perspective of CLM4 for the evident discrepancies of GPP as shown in the results section.

b. Evaluation of tropical GPP

Spatial correlation [ section 3a(1)] and EOF analysis (Fig. 7) suggest that CLM4 and MODIS algorithms, in spite fundamental differences in structure and inputs, generate very similar information regarding the climatological mean state of the global GPP distribution at annual and seasonal scales. First-order differences between climatological means for the two approaches include high estimates of tropical evergreen forest GPP from CLM4, longer high-latitude growing season in CLM4, and stronger seasonality of needleleaf evergreen leaf area from MODIS.

Mean GPP for tropical evergreen forest is 360 gC m−2 month−1 in CLM4, compared to 214 gC m−2 month−1 for MODIS and 194 gC m−2 month−1 from a third independent estimate (Beer et al. 2010). A possible explanation for this difference is related to the current treatment of nutrient dynamics and nutrient limitations on GPP in CLM4. The availability of mineral nitrogen provides an important constraint on GPP in CLM4, but other nutrient cycles have so far been ignored. Of particular concern is the phosphorus cycle since phosphorus is generally believed to be more limiting than nitrogen to tree growth in lowland tropical moist forests (Townsend et al. 2011). While the current model predicts both instantaneous and long-term constraints on GPP from nitrogen limitation in this region (Fig. 8 in Thornton and Zimmermann 2007), we expect that the real limitation from phosphorus availability in this region is even stronger. We are currently developing a version of CLM4 that includes a prognostic phosphorus cycle as a partial test of this hypothesis. Additional evaluation of this hypothesis against site-level observations and experimental manipulations is also warranted.

Another possible explanation for high CLM4 GPP in the tropics has to do with the representation of the vertical canopy gradient in specific leaf areas. The linear relationship proposed by Thornton and Zimmerman (2007) is a good empirical fit to their observations, but it permits unrealistically high leaf area at high canopy leaf carbon content. Truncation of the relationship due to physical limits on (low) leaf thickness would reduce high bias in LAI. However, because the model assumes, in correspondence with observations, a constant leaf C:N ratio regardless of canopy position, this modification is unlikely to have a strong influence on canopy-scale GPP. Preliminary studies suggest that CLM4 GPP and biomass in this region are sensitive to the specification of tree mortality rates. As we have specified a single nonfire mortality rate for all vegetation types worldwide (0.5% yr−1), it is possible that refinements by climate zone or PFT could influence the tropical forest GPP differences.

c. Evaluation of high-latitude GPP

Several factors make interpreting high-latitude differences in the predicted length of growing season and seasonal cycles of GPP and LAI difficult. For example, snow cover and low light obscure remote sensing signals in early and late winter seasons at high latitudes. Also, the phenology algorithms in CLM4 have been derived mainly from temperate zone observations, so we have a low a priori confidence in model performance in the boreal forest. We find the strong seasonality in boreal evergreen forest LAI from MODIS suspicious, given the multiyear leaf longevity observed for these forests and the dominance of forest cover imposed by our analysis (section 2d). This suggests that some combination of snow cover and low light is producing a low bias in the MODIS LAI product for this PFT and that snow cover may also play a role in the strong seasonality of MODIS LAI in temperate evergreen needleleaf forest (Figs. 6a,b). It is also possible that the model controls on GPP at low temperatures are not stringent enough, leading to higher CLM GPP in the shoulder seasons in the boreal zone. Evaluations of model predictions against observations at individual sites, including boreal deciduous and evergreen sites, are under way and should help to clarify this issue.

d. Evaluation of interannual GPP variability in long-term trends

CLM and MODIS in GPP variability on seasonal time scales (Fig. 7) are in broad agreement; however, their interannual variability, long-term trends, and correlations of GPP with potential forcing factors and related ecosystem states are quite different. We note with particular interest that interannual variations in GPP and LAI as predicted by CLM are highly correlated in both hemispheres, while GPP and LAI derived from MODIS observations have no correlation in the Northern Hemisphere and only a weak (nonsignificant) correlation in the Southern Hemisphere. CLM GPP and LAI also show strong positive correlation with precipitation and strong negative correlation with incident shortwave radiation in both hemispheres. MODIS GPP shows weaker, but still significant, correlations with these variables in the Northern Hemisphere and nonsignificant correlations in the Southern Hemisphere. MODIS LAI has no significant correlations with precipitation or incident shortwave radiation in either hemisphere. Both MODIS and CLM GPP show positive (negative) correlation with temperature in the Northern (Southern) Hemisphere, but these correlations are strong and significant for MODIS and nonsignificant for CLM.

It is not clear if these differences reflect biases in CLM, MODIS products, or both. Our analysis does, however, suggest a testable hypothesis, which would help to resolve these differences. We hypothesize that interannual variations in precipitation have a direct and positive influence on interannual variation in GPP at global and hemispheric scales and that increased cloudiness associated with higher precipitation leads to a negative correlation between GPP and incident radiation. An additional factor in this causal framework is the direct influence of increased fraction of diffuse radiation under cloudier conditions, which could also lead to enhanced GPP for limited ranges of variation in total incident shortwave radiation. We further hypothesize that the same causal framework has a direct effect on LAI, as LAI is determined at least in part by current year GPP. Finally, we hypothesize a reinforcing effect of the variation in LAI on variation in GPP through the variation in fractional absorption of photosynthetically active radiation at higher LAI. Evaluations against observed GPP at individual flux measurement sites in both hemispheres are underway and should provide at least a partial test of this set of hypotheses, helping to quantify the prediction error associated with the CLM and MODIS GPP and LAI products.

e. Next steps for hypothesis testing

Our extension of the evaluation process using satellite GPP provides alternatives for understanding the strengths and weaknesses of process-oriented models. It also identifies several areas that should have priority in further evaluating and improving GPP in the satellite observation. From the perspective of remote sensing products, the 10-yr MODIS GPP is relatively short and restricts further intercomparison on the interannual and even decadal scales. Therefore, we speculate that more long-term datasets and the range of uncertainty associated with each product still need to be created. In view of the discrepancies of CLM4 performance on GPP, we would thus need to test the various factors controlling photosynthetic production and phenological parameterizations in the CLM4 to quantify potential errors arising from the representation of ecosystem mechanisms themselves. Future work includes tests such as exploring how uncertainties from external inputs, internal parameters, and model structure propagate to GPP. Furthermore, because no single measure is adequate to evaluate the performance of CLM4 GPP, we strongly propose more comprehensive estimations such as tower flux measurements, data-based diagnostic approaches, and multimodel intercomparisons to systematically investigate and improve the predictions of GPP, which would also improve simulations of the entire biogeochemical cycle in CESM.

5. Conclusions

Global intercomparisons and multistatistics have been analyzed to assess the GPP biases between the CLM4 and MOD17 GPP for a full 10-yr period between 2000 and 2009. We conclude that the CLM4 is in rough agreement with the remotely sensed primary production in describing the large-scale seasonal–spatial patterns and temporal variations. However, some substantial GPP differences still remain. CLM4 systematically overestimated tropical GPP, underestimated the boreal deciduous needleleaf tree and boreal deciduous shrub primary production, predicted a stronger and longer GPP seasonal cycle over some plant function types, and showed stronger sensitivity to interannual variation in precipitation. We explored the possible reasons behind those GPP differences between the two products. Several interpretations and testable hypotheses about existing GPP errors were addressed. These understandings, combined with our intercomparison results, highlight key steps necessary for the improvement of CLM4 diagnostics of biogeochemical cycles at multiple temporal and spatial scales.

Acknowledgments

This research is supported in part by the U.S. Department of Energy (DOE), Office of Science, Biological and Environmental Research. Oak Ridge National Laboratory is managed by UT-BATTELLE for the DOE under Contract DE-AC05-00OR22725. Special thanks are given to Dr. Sam Levis at NCAR for his help with CLM4 simulations and to Terry Copeland Pfeiffer at ORNL for her text editing.

REFERENCES

REFERENCES
Baker
,
I. T.
,
A. S.
Denning
, and
R.
Stöckli
,
2010
:
North American gross primary productivity: Regional characterization and interannual variability
.
Tellus
,
62B
,
533
549
.
Beer
,
C.
, and
Coauthors
,
2010
:
Terrestrial gross carbon dioxide uptake: Global distribution and covariation with climate
.
Science
,
329
,
834
838
.
Bonan
,
G. B.
, and
S.
Levis
,
2010
:
Quantifying carbon-nitrogen feedbacks in the Community Land Model (CLM4)
.
Geophys. Res. Lett.
,
37
,
L07401
,
doi:10.1029/2010GL042430
.
Bonan
,
G. B.
,
P. J.
Lawrence
,
K. W.
Oleson
,
S.
Levis
,
M.
Jung
,
M.
Reichstein
,
D. M.
Lawrence
, and
S. C.
Swenson
,
2011
:
Improving canopy processes in the Community Land Model version 4 (CLM4) using global flux fields empirically inferred from FLUXNET data
.
J. Geophys. Res.
,
116
,
G02014
,
doi:10.1029/2010JG001593
.
Ciais
,
P.
, and
Coauthors
,
1997
:
A three-dimensional synthesis study of δ18O in atmospheric CO2 1. Surface fluxes
.
J. Geophys. Res.
,
102
(
D5
),
5857
5872
.
Coops
,
N. C.
,
C. J.
Ferster
,
R. H.
Waring
, and
J.
Nightingale
,
2009
:
Comparison of three models for predicting gross primary production across and within forested ecoregions in the contiguous United States
.
Remote Sens. Environ.
,
113
,
680
690
.
Cox
,
P. M.
,
R. A.
Betts
,
C. D.
Jones
,
S. A.
Spall
, and
I. J.
Totterdell
,
2000
:
Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model
.
Nature
,
408
,
184
187
.
Cramer
,
W.
, and
Coauthors
,
1999
:
Comparing global models of terrestrial net primary productivity (NPP): Overview and key results
.
Global Change Biol.
,
5
,
1
15
.
Frankenberg
,
C.
, and
Coauthors
,
2011
:
New global observations of the terrestrial carbon cycle from GOSAT: Patterns of plant fluorescence with gross primary productivity
.
Geophys. Res. Lett.
,
38
,
L17706
,
doi:10.1029/2011GL048738
.
Friedlingstein
,
P.
, and
Coauthors
,
2006
:
Climate–carbon cycle feedback analysis: Results from the C4MIP model intercomparison
.
J. Climate
,
19
,
3337
3353
.
Fung
,
I. Y.
,
S. C.
Doney
,
K.
Lindsay
, and
J.
John
,
2005
:
Evolution of carbon sinks in a changing climate
.
Proc. Natl. Acad. Sci. USA
,
102
,
11 201
11 206
.
Heinsch
,
F. A.
, and
Coauthors
,
2006
:
Evaluation of remote sensing based terrestrial productivity from MODIS using regional tower eddy flux network observations
.
IEEE Trans. Geosci. Remote Sens.
,
44
,
1908
1925
.
Jones
,
C.
,
J.
Lowe
,
S.
Liddicoat
, and
R.
Betts
,
2009
:
Committed terrestrial ecosystem changes due to climate change
.
Nat. Geosci.
,
2
,
484
487
.
Jung
,
M.
, and
Coauthors
,
2007
:
Uncertainties of modeling gross primary productivity over Europe: A systematic study on the effects of using different drivers and terrestrial biosphere models
.
Global Biogeochem. Cycles
,
21
,
GB4021
,
doi:10.1029/2006GB002915
.
Jung
,
M.
,
M.
Reichstein
, and
A.
Bondeau
,
2009
:
Towards global empirical upscaling of FLUXNET eddy covariance observations: Validation of a model tree ensemble approach using a biosphere model
.
Biogeosciences
,
6
,
5271
5304
.
Justice
,
C. O.
,
J. R. G.
Townshend
,
E. F.
Vermote
,
E.
Masuoka
,
R. E.
Wolfe
,
N.
Saleous
,
D. P.
Roy
, and
J. T.
Morisette
,
2002
:
An overview of MODIS land data processing and product status
.
Remote Sens. Environ.
,
83
,
3
15
.
Kanamitsu
,
M.
,
W.
Ebisuzaki
,
J.
Woolen
,
S.-K.
Yang
,
J. J.
Hnilo
,
M.
Fiorino
, and
G. L.
Potter
,
2002
:
NCEP–DOE AMIP-II Reanalysis (R-2)
.
Bull. Amer. Meteor. Soc.
,
83
,
1631
1643
.
Kanniah
,
K. D.
,
J.
Beringer
,
L. B.
Hutley
,
N. J.
Tapper
, and
X.
Zhu
,
2009
:
Evaluation of collections 4 and 5 of the MODIS gross primary productivity product and algorithm improvement at a tropical savanna site in northern Australia
.
Remote Sens. Environ.
,
113
,
1808
1822
.
Lawrence
,
D. M.
, and
Coauthors
,
2011
:
Parameterization improvements and functional and structural advances in version 4 of the Community Land Model
.
J. Adv. Model. Earth Syst.
,
3
,
M03001
,
doi:10.1029/2011MS000045
.
Lawrence
,
P. J.
, and
T. N.
Chase
,
2007
:
Representing a new MODIS consistent land surface in the Community Land Model (CLM 3.0)
.
J. Geophys. Res.
,
112
,
G01023
,
doi:10.1029/2006JG000168
.
Mitchell
,
T. D.
, and
P. D.
Jones
,
2005
:
An improved method of constructing a database of monthly climate observations and associated high-resolution grids
.
Int. J. Climatol.
,
25
,
693
712
.
Myneni
,
R. B.
, and
Coauthors
,
2002
:
Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data
.
Remote Sens. Environ.
,
83
,
214
231
.
Nemani
,
R. R.
,
C. D.
Keeling
,
H.
Hashimoto
,
W. M.
Jolly
,
S. C.
Piper
,
C. J.
Tucker
,
R. B.
Myneni
, and
S. W.
Running
,
2003
:
Climate-driven increases in global terrestrial net primary production from 1982 to 1999
.
Science
,
300
,
1560
1563
.
North
,
G. R.
,
T. L.
Bell
,
R. F.
Cahalan
, and
F. J.
Moeng
,
1982
:
Sampling errors in the estimation of empirical orthogonal functions
.
Mon. Wea. Rev.
,
110
,
699
706
.
Oleson
,
K. W.
,
G. B.
Bonan
,
C.
Schaaf
,
F.
Gao
,
Y. F.
Jin
, and
A.
Strahler
,
2003
:
Assessment of global climate model land surface albedo using MODIS data
.
Geophys. Res. Lett.
,
30
,
1443
,
doi:10.1029/2002GL016749
.
Oleson
,
K. W.
, and
Coauthors
,
2010
:
Technical description of version 4.0 of the Community Land Model (CLM). NCAR Tech. Note NCAR/TN-478+STR, 257 pp
.
Randerson
,
J. T.
, and
Coauthors
,
2009
:
Systematic assessment of terrestrial biogeochemistry in coupled climate–carbon models
.
Global Change Biol.
,
15
,
2462
2484
.
Running
,
S. W.
,
D. D.
Baldocchi
,
D. P.
Turner
,
S. T.
Gower
,
P. S.
Bakwin
, and
K. A.
Hibbard
,
1999
:
A global terrestrial monitoring network integrating tower fluxes, flask sampling, ecosystem modeling and EOS satellite data
.
Remote Sens. Environ.
,
70
,
108
127
.
Running
,
S. W.
,
R. R.
Nemani
,
F. A.
Heinsch
,
M. S.
Zhao
,
M.
Reeves
, and
H.
Hashimoto
,
2004
:
A continuous satellite-derived measure of global terrestrial primary production
.
Bioscience
,
54
,
547
560
.
Shi
,
X.
,
J.
Mao
,
P. E.
Thornton
,
M. H.
Forrest
, and
M. P.
Wilfred
,
2011
:
The impact of climate, CO2, nitrogen deposition and land use change on simulated contemporary global river flow
.
Geophys. Res. Lett.
,
38
,
L08704
,
doi:10.1029/2011GL046773
.
Sitch
,
S.
, and
Coauthors
,
2008
:
Evaluation of the terrestrial carbon cycle, future plant geography and climate-carbon cycle feedbacks using five Dynamic Global Vegetation Models (DGVMs)
.
Global Change Biol.
,
14
,
2015
2039
.
Solomon
,
S.
,
D.
Qin
,
M.
Manning
,
M.
Marquis
,
K.
Averyt
,
M. M. B.
Tignor
,
H. L.
Miller
Jr.
, and
Z.
Chen
, Eds.,
2007
:
Climate Change 2007: The Physical Science Basis. Cambridge University Press, 996 pp
.
Stöckli
,
R.
,
T.
Rutishauser
,
D.
Dragoni
,
J.
O’Keefe
,
P. E.
Thornton
,
M.
Jolly
,
L.
Lu
, and
A. S.
Denning
,
2008
:
Remote sensing data assimilation for a prognostic phenology model
.
J. Geophys. Res.
,
113
,
G04021
,
doi:10.1029/2008JG000781
.
Thornton
,
P. E.
, and
N. A.
Rosenbloom
,
2005
:
Ecosystem model spin-up: Estimating steady state conditions in a coupled terrestrial carbon and nitrogen cycle model
.
Ecol. Modell.
,
189
,
25
48
.
Thornton
,
P. E.
, and
N. E.
Zimmermann
,
2007
:
An improved canopy integration scheme for a land surface model with prognostic canopy structure
.
J. Climate
,
20
,
3902
3923
.
Thornton
,
P. E.
, and
Coauthors
,
2009
:
Carbon-nitrogen interactions regulate climate-carbon cycle feedbacks: Results from an atmosphere-ocean general circulation model
.
Biogeosciences
,
6
,
2099
2120
.
Tian
,
Y.
, and
Coauthors
,
2004
:
Comparison of seasonal and spatial variations of leaf area index and fraction of absorbed photosynthetically active radiation from Moderate Resolution Imaging Spectroradiometer (MODIS) and Common Land Model
.
J. Geophys. Res.
,
109
,
D01103
,
doi:10.1029/2003JD003777
.
Townsend
,
A. R.
,
C. C.
Cleveland
,
B. Z.
Houlton
,
C. B.
Alden
, and
J. W. C.
White
,
2011
:
Multi-element regulation of the tropical forest carbon cycle
.
Front. Ecol. Environ.
,
9
,
9
17
,
doi:10.1890/100047
.
Turner
,
D. P.
,
W. D.
Ritts
,
M. S.
Zhao
,
S. A.
Kurc
,
A. L.
Dunn
,
S. C.
Wofsy
,
E. E.
Small
, and
S. W.
Running
,
2006a
:
Assessing interannual variation in MODIS-based estimates of gross primary production
.
IEEE Trans. Geosci. Remote Sens.
,
44
,
1899
1907
.
Turner
,
D. P.
, and
Coauthors
,
2006b
:
Evaluation of MODIS NPP and GPP products across multiple biomes
.
Remote Sens. Environ.
,
102
,
282
292
.
Wang
,
W.
,
J.
Dungan
,
H.
Hashimoto
,
A. R.
Michaelis
,
C.
Milesi
,
K.
Ichil
, and
R. R.
Nemani
,
2011
:
Diagnosing and assessing uncertainties of terrestrial ecosystem models in a multimodel ensemble experiment: 1. Primary production
.
Global Change Biol.
,
17
,
1350
1366
.
Wang
,
Z.
,
X.
Zeng
,
M.
Barlage
,
R. E.
Dickinson
,
F.
Gao
, and
C. B.
Schaaf
,
2004
:
Using MODIS BRDF and albedo data to evaluate global model land surface albedo
.
J. Hydrometeor.
,
5
,
3
14
.
Xiao
,
J.
, and
Coauthors
,
2008
:
Estimation of net ecosystem carbon exchange for the conterminous United States by combining MODIS and AmeriFlux data
.
Agric. For. Meteor.
,
148
,
1827
1847
.
Zaehle
,
S.
,
S.
Sitch
,
B.
Smith
, and
F.
Hatterman
,
2005
:
Effects of parameter uncertainties on the modeling of terrestrial biosphere dynamics
.
Global Biogeochem. Cycles
,
19
,
GB3020
,
doi:10.1029/2004GB002395
.
Zhao
,
M.
, and
S. W.
Running
,
2010
:
Drought-induced reduction in global terrestrial net primary production from 2000 through 2009
.
Science
,
329
,
940
943
.
Zhao
,
M.
, and
S. W.
Running
,
2011
:
Response to comments on “Drought-induced reduction in global terrestrial net primary production from 2000 through 2009.”
Science
,
333
,
1093
.
Zhao
,
M.
,
F. A.
Heinsch
,
R. R.
Nemani
, and
S. W.
Running
,
2005
:
Improvements of the MODIS terrestrial gross and net primary production global data set
.
Remote Sens. Environ.
,
95
,
164
176
.
Zhao
,
M.
,
S. W.
Running
, and
R. R.
Nemani
,
2006
:
Sensitivity of Moderate Resolution Imaging Spectroradiometer (MODIS) terrestrial primary production to the accuracy of meteorological reanalyses
.
J. Geophys. Res.
,
111
,
G01002
,
doi:10.1029/2004JG000004
.

Footnotes

*

Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-11-00401.s1.

Supplemental Material