Terrestrial ecosystem production: A process model based on global satellite and surface data

. This paper presents a modeling approach aimed at seasonal resolution of global climatic and edaphic controls on patterns of terrestrial ecosystem production and soil microbial respiration. We use satellite imagery (Advanced Very High Resolution Radiometer and International Satellite Cloud Climatology Project solar radiation), along with historical climate (monthly temperature and precipitation) and soil attributes (texture, C and N contents) from global (1 o) data sets as model inputs. The Carnegie-Ames-Stanford approach (CASA) Biosphere model runs on a monthly time interval to simulate seasonal patterns in net plant carbon fixation, biomass and nutrient allocation, litterfall, soil nitrogen mineralization, and microbial CO2 production. The model estimate of global terrestrial net primary production is 48 Pg C yr -(cid:127) with a maximum light use efficiency of 0.39 g C MJ -(cid:127) PAR. Over 70% of terrestrial net production takes


INTRODUCTION
The amount of carbon fixed annually via terrestrial net primary productivity (NPP) or released by soil microbial respiration (Rs) is about an order of magnitude greater than the annual increase in atmospheric carbon dioxide (CO2) levels due to fossil fuel combustion [Ajtay, 1979;Houghton et al., 1992].Seasonal changes in the balance between photosynthetic carbon fixation by land plants and microbial respiration are of a size sufficient to drive the intra-annual oscillation of atmospheric CO2 concentration [Bacastow et al., 1985;Houghton, 1987].Either carbon fixation or respiration could be affected substantially by components of global change (e.g., warming or elevated CO2 concentrations), which raises the possibility of longterm modifications in the carbon balance of terrestrial ecosystems [Rastetter et al., 1991;Mooney et al., 1991] and feedbacks to global biogeochemistry and radiative forcing.
All global data sets used as inputs to the model (Table 1) were resampled (if necessary) to a 10x 1 o spatial resolution.Specific data sources are described in the following paragraphs.

Vegetation Index
We used monthly NDVI-FASIR (defined below) data sets for 1987 processed by Los et al. [1993]   texture classes and their associated particle size distributions (Table 3).Data used to construct the SMW has been assembled from actual soil surveys (21% of global coverage); field reconnaissance of topography, geology, vegetation and climatic data (40% of global coverage); and general information from the local literature (39% of global coverage) [Gardiner, 1982].Substantial disparities in the reliability of soil type classification have been identified over areas of tropical Central and South America and Africa [Gardiner, 1982;Richter and Babbar, 1991 ].

Soil Carbon and Nitrogen Contents
We created global gridded data sets by mapping Holdridge [ 1967] life zone soil C and N content (g m -3) averages reported by Post et al. [1985] to their corresponding life zone categories produced by Leemans [1990].Spatial interpolation of the resulting data sets from 0.5 ø to 1.0 ø latitude/longitude was accomplished using bidirectional splining, preceded by nearest-neighbor fill to conserve land-water boundary elements.Soil profiles used in this data set were all from natural vegetation and excluded wetlands.Post et al. [1985] considered coverage of the original data to be best for tropical and cool temperate forests; coverage is poorer over extremely wet areas, dry tundra, dry boreal and warm desert life zones.

Other Variable Definitions
Ranges of certain other model variables were estimated from the literature, as discussed below.
These include leaf:root:wood C and N allocation ratios, litter and soil organic matter decomposition rates, and C assimilation efficiency of microbes.To simplify interpretations in this version of the model, we set spatially uniform values for most of these variables.As part of the modeling process, uniform rate constants related to photosynthesis and microbial respiration fluxes are adjusted for temporally and spatially resolved stress effects.

Net Primary Productivity
New production of plant biomass (NPP) at a grid cell (x) in month t is a product of intercepted photosynthetically active radiation (IPAR) and a light utilization efficiency (e) that is modified by temperature and soil moisture (equation ( 1)).Neither IPAR nor e is dependent on ecosystem type.

The NPP formulation allows for regulation in either of the terms on the right side of equation (1).
Several lines of evidence indicate that most of the regulation should be in IPAR, with less in e.One line of evidence comes from surveys which indicate that NPP of many ecosystem types is highly correlated with the annual integral of NDVI [Goward et al., 1985].Another is the constancy of e from many experimental studies of unstressed plants, plus results from several studies indicating that nutrient stress [Garcia et al., 1988] and water stress [Squire et al., 1986] have much larger effects on IPAR than on •.Field [1991]  For mean monthly temperatures of-10 ø C and below, Tel is set equal to zero.The basic motivation for including Tel is to reflect the empirical observation that plants in very cold habitats typically have low maximum growth rates [Chapin, 1980;Grime, 1979] and high root biomass [Sala et al., 1993], potentially imposing large respiratory costs.Plants in very hot environments may have high growth rates [Schulze and Chapin, 1987], but the efficiency of light utilization should be impacted by high rates of respiration [Amthor, 1989;Ryan, 1991 including Te2 is to capture some of the intrinsic limitations on the flexibility of temperature  2 and listed in Table 4).
The motivation for this single step calibration is that we lack the understanding to estimate e* from first principles.The maximum photon yield of photosynthesis sets an absolute upper bound on •* of approximately 2.88 g C MJ -1 PAR [Ehleringer and Bj6rkman, 1977], but this value will always be reduced by saturation in the light response of photosynthesis.At the Harvard Forest, Wofsy et al.
[1993] measured a gross primary productivity light use efficiency during summer months of around 1.10 g C MJ -1 PAR.A 2:1 ratio of gross primary productivity to net primary productivity would argue for a maximum •* of around 0.55 g C MJ -1 PAR at Harvard Forest, fairly close to the •* obtained through the above calibration.

Soil Moisture Submodel
Soil moisture content was calculated at each grid cell using monthly temperature and precipitation in combination with soil texture data and moisture-holding capacity [Saxton et al., 1986].This submodel is a one-layer "bucket" formulation that builds on previous simulation studies of regional and global surface hydrology [Mintz and Serafini, 1981;VOr6smarty et al., 1989;Bouwman et al., 1993].
Monthly soil moisture storage is calculated for each grid cell (x) as a state variable, SOILM, with the potential to accumulate moisture over several months.

(S OILM(x,t-1)).
For months when temperature is less than or equal to 0 ø C, PET and PPT are set equal to zero and there is no net change in SOILM.During these months, precipitation accumulates as snow in a state variable PACK.
PACK is added to PPT in the first month that monthly average air temperature (T) > 0 o C.This function initiates spring snow melt.
PET is calculated with the method of Thomthwaite [1948].Lower and upper limits for SOILM were set at wilting point (WPT) and field capacity (FC), respectively (Table 3).These values were derived from soil texture relationships described by Saxton et     3).5).To avoid the potential for complete depletion of N pools, we fixed the priority for potentially competing N transfers between litter, microbial, soil organic matter, and mineral pools.

Comparison with Other Models
There are several notable differences between the CASA-Biosphere model (this study) and previous ecosystem models with carbon-nitrogen feedbacks that have been applied at regional and global scales, particularly CENTURY [Parton et al., 1987[Parton et al., , 1988[Parton et al., , and 1989a] ]  The CASA soil submodel is baged on the CENTURY concept, but has been simplified for global applications.For example, the CASA version does not include optional effects of pH, phosphorus, sulfur, and atmospheric-hydrologic fluxes of nitrogen.CENTURY also includes algorithms for competition between trees and grasses for soil N.
Unlike CENTURY, but like TEM, we make no correction for differences between air and soil surface temperature.In the CASA model, the C-to-N ratios of SLOW and OLD soil pools are set according to soil texture [Zobler, 1986], whereas in CENTURY these ratios float over a prescribed range.In both models, N is transferred in proportion to monthly NPP at the amount needed to produce the highest N quality litter possible over a prescribed range.

Net Primary Productivity
The CASA-Biosphere model was applied using global data inputs as described above.On the basis of monthly calculations over 14,713 non-ice terrestrial grid cells, our estimate of total terrestrial net primary production (NPP) is 48 Pg C yr -1.Over 34 Pg C yr -1, or 70% of terrestrial NPP, occurs between 30 ø N and 30 ø S latitude, in contrast to 4% at latitudes higher than 60 ø N and S (Table 7)    Note that all stresses on ecosystem production also show up in IPAR (Table 8), and that the sum of stresses on both terms on the right hand side of equation ( 1

Litterf all Patterns
The seasonal timing of litterfall is latitude dependent (Figure 6  Here S is the total annual amount of PAR at the surface, z

Biome Contributions to Global NEP
For northern hemisphere broadleaf evergreen forests, mean NPP is greater than mean Rs (g m -2 mo-1 over 1 o latitude zones) from September to February (class 1; Figure 7a).In the southern hemisphere, the NEP pattern for the same biome is reversed.Northern midlatitude to high-latitude forests show a positive NEP from June to August (Figure 7b).If these per-unit-area averages are corrected for global area coverage of the biomes, needleleaf evergreen forests (class 4) exchange about six times more carbon than broadleaf deciduous forests (class 2), an observation explained largely by over sevenfold greater area coverage by class 4. Northern grasslands (class 7) and cultivated (class 12) vegetation (Figure 7c) have a low seasonal NEP amplitude relative to temperate forest and northern tundra biome types (class 10; Figure 7d).The lowest seasonal amplitude in NEP is in biomes characterized as bare soil and deserts (class 11), paralleling their small contributions to NPP (Table 7).Litter and soil carbon pools (g C m -2) are greatest in mixed temperate forests (class 3, Table 9).Perennial grasslands (class 7) and deserts (class 11), have the lowest litter pools, both per unit area and integrated over area (Table 9).
Turnover times of four litter-soil carbon pools (expressed as the steady state pool size divided by the sum of annual losses) are fastest for biomes found chiefly in the tropics (broadleaf evergreen forests and savanna vegetation classes) and for cultivated areas (Table 10).Carbon turnover in needleleaf deciduous forests and tundra is slowest.

Correlations with Atmospheric C02
As a consistency check of CASA-Biosphere model predictions, we compared our monthly NEP estimates to average atmospheric CO2 concentrations from three stations in the NOAA/Geophysical Monitoring for Climate Change (GMCC) Flask

Model Result Comparisons of Production
The  9).One reason for this pattern is that overall litter decomposition rates can be underestimated using the fractionation algorithm shown in equation ( 16

Fig. 1 .
Fig. 1.Model integration framework.Global climate data sets are combined with soil texture settings to compute the monthly water balance, which controls NPP and soil microbial activity.
considers ecological factors that should tend to constrain investments in light harvesting (which are manifested as IPAR) in relation to whatever resource or resources are limiting to growth so that all of the IPAR can be used for growth.A strong relationship between NPP and IPAR does not necessarily indicate that light is the primary resource limiting to growth [Te2 account for effects of temperature stress, We accounts for effects of water stress, and e* is the maximum possible efficiency.The two temperature stress terms serve to depress e at very high and very low temperatures (Tel) and to depress e when the temperature is above or below e[ 0'3 (-Topt(x) -10 + T(x,t))] } (7) Te2 falls to half its value at Topt at temperatures 10 ø C above or 13 ø C below Top t.The idea behind

Fig. 3 .Fig. 4 .
Fig. 3. Soil drying curves derived using a mathematical transformation equation (12) of the relationship between soil water potential and volumetric moisture content presented by Saxton et al. [ 1986].

Fig. 5 .[
Fig. 5. Ecosystem carbon-nitrogen model.Carbon pools are outlined in black and labeled with C-to-N ratios, C fluxes in solid arrows, CO2 production in stippled arrows; Nitrogen pools in gray, N fluxes in gray arrows.Levels of litter, microbe (MIC), and soil organic (SLOW and OLD) pools are shown.Structural (S) and metabolic (M) pools are shown for leaf and root litter.
is satisfied first by computing the amount of N needed, relative to the corresponding carbon side flow, to maintain the critical C:N of the SLOW pool.If there is insufficient N available in the structural pool to meet 100% of this demand, no outflows from the structural source pool occurs during that time step.If there is sufficient N to satisfy the Snl flow, Sn2 is then computed as the lesser of (1) the amount needed to maintain microbe pool C:N at 10, or (2) the difference between N available in the source pool and Snl.Last, if there is sufficient N available in the source pool to satisfy both the Snl and Sn2 flows, Sn3 is computed as difference between the potential amount of N lost from the source pool (relative to carbon decomposition) and (Snl + Sn2).Likewise, for the metabolic leaf and root litter pools, there are two possible N outflow paths; transfer to microbial pools is satisfied first by computing the amount of N needed to maintain the critical C:N of microbe pools.The three possible outflow paths from microbial N pools, and the two possible outflow paths from each soil (SLOW and OLD) N pool are computed according to the same rules and priorities.There is no direct flux from the SLOWn to the OLDn pool; transfer is indirect via the soil microbial pool.Following mineralization calculations, immobilization flow from the MIN n pool to each litter pool is computed.These flows are dependent on total immobilization demand, which is determined as the sum of all N flows needed to meet the critical C:N ratio of individual litter pools.If total demand exceeds the current MINn pool level, each litter pool whose C:N exceeds the critical level receives an equal portion of the MINn pool.If total demand is less than the MINn level, demand is satisfied completely for each litter pool whose C:N is higher than the critical level.The remaining MINn is available for uptake by vegetation.Impacts of Cultivation on Microbial RespirationConversion of native ecosystems to cropping systems through cultivation involves some degree of vegetation clearing, soil tillage and fertilization.We make two modifications for cultivated grid cells.First, plant litter lignin concentrations are set at 5% and the minimum C-to-N ratio of litter is adjusted and TEM [Raich et al., 1991; McGuire et al., 1992; Melillo et al., 1993].The major differences (as summarized in Table 6) include sources and spatial resolutions of climate and vegetation data sets used, the degree to which parameter values are calibrated to specific vegetation types, use of remote sensing observations, and detail of the soil moisture and C/N cycling submodels.In the TEM, NPP is computed as gross primary productivity (GPP) minus respiration, using atmospheric CO2 concentrations and solar radiation drivers.It includes unitless scalars for the effects of air temperature, relative N availability, and plant phenology.In this version of the CASA model, we do not compute plant respiration, nor do we consider CO2 concentrations and N limitation explicitly as controllers of productivity.The TEM equation for N uptake is hyperbolic function of mineral N and temperature, carbon availability, and a biome-specific maximum uptake constant.Other biome-specific parameters in TEM that are not biome-specific in CASA include the timing of litterfall and microbial respiration rates.The TEM uses different climate data sources, vegetation type classification, and a finer spatial resolution (0.5 ø versus 1.0o).
. For individual grid cells, annual NPP ranges from values less than 20 g C m -2 in some tundra and desert locations to values greater than 1400 g C m -2 in broadleaf evergreen forests of South America and Central Africa (Plate 2).Approximately one third of all solar radiation impinging on terrestrial surfaces is intercepted by Global map of NPP computed according to the model in equation (1).
) is adjusted in the calibration process (Figure 2).As an example of the co-occurrence of high efficiencies and high IPAR valuesaverage light use efficiency is approximately 54% of •*.These two ecosystems comprise a significant fraction of the total land surface area in which growth is constrained by water and temperature stress on the light use efficiency term.Production estimates for these areas are depressed far below what is allowed by •*.The overall effect of making a sensitive to water and temperature stress is to expand the dynamic range of possible NPP estimates consistent with a given distribution of IPAR.
).The highest monthly percent litterfall occurs during October between 60 ø and 65 ø N. Strong seasonality in litterfall is seen in the northern temperate zone as far south 40 ø N. Distinct seasonality can be detected in the tropics from 5 Fig. 6.Seasonal patterns of litterfall over averaged 1 o latitude zones.Contour labels are percent of total annual litterfall.
Fig. 8. Latitude-based litter and soil carbon pool totals (Pg C).

TEM
those used in TEM for our NPP estimates, the modeling approaches are sufficiently different that continental estimates are only loosely constrained.On a continental basis, our model estimates production for North and South America at 6.1 and 14.4 Pg C yr -1, respectively, compared to than the CASA model prediction.The lower CASA NPP estimates for many biomes may be due in part to the sensitivity of our model to IPAR.Even though the mean climate -NPP regressions developed in the MIAMI model may suggest a potential NPP for a grid cell, there can be areas within the cell where plants do not persist (e.g., surface area covered by rock, lakes or asphalt) or times of the year when human management has altered the land cover.These bare areas will lower the average IPAR for the grid cell and thus reduce the average NPP calculated in CASA for the biome type on a per square meter basis.For example, the CASA predicted NPP for class 11 (bare soil and desertsdepth.Our prediction of nonwoody surface litter (51 Pg C) is close to Schlesinger's [1977] estimate of 55 Pg C. The CASA prediction for standing litter pools in the tropics is somewhat higher than expected, compared to those for temperate forests (Table
Fig. 10.Comparison of model NPP model predictions by vegetation classes from Dorman and Sellers [ 1989].Shaded bars are MIAMI model [Lieth, 1975] predictions; hatched bars are CASA model (this study) predictions.

Organic soils were assigned to the coarse/medium texture class [Bouwman et al., 1993]. aFrom Zobler [ 1986] bField capacity (m) for forested grid cells; FC(x) for other vegetation types are 50% of these values. Computed based on equation (2) in the work by Saxton et al. [ 1986] for soil water tension greater than or equal to 10 kPa. Tension was assumed to be 33 kPa for medium to fine textures and 10 kPa for coarse textures [Papendick and Campbell, 1980]. cWilting point (m) for forested grid cells; WP(x) for other vegetation types are 50% of these values. Computed based on equation (2) in the work by Saxton et al. [1986] for soil water tension equal to 1500 kPa. dFrom Saxton et al
. [1986]; used in calculation of RDR.eOn the basis of weighted average particle-size C:N values reported by Anderson et al. [1981 ], Hinds and Lowe [1980], and Cameron and Posner [1979].

5 (2) where SOL is the total solar radiation incident on grid cell x in month t, from the database of Bishop and Rossow [1991 ], FPAR is the fraction of the incoming PAR intercepted by green vegetation, and the factor of 0.5 accounts for the fact that approximately half of the incoming solar radiation is in the PAR waveband (0.4-0.7 [tm) [McCree, 1981]. FPAR is calculated as a linear function of the
IPAR(x,t) = SOL(x,t) FPAR(x,t) 0.

where PPT is average precipitation at month t, PET is potential evapotranspiration at month t, and RDR is a relative drying rate scalar for potential water extraction as a function of soil moisture
[PET(x,t) -PPT(x,t)] RDR For PPT(x,t) < PET(x,t) SOILM(x,t) = SOILM(x,t-1) + [PPT(x,t) -PET(x,t)] For PPT(x,t) > PET(x,t)

(14c) NPPann is total annual NPP. The LAImin and LAIave terms are the minimum and the average of 12 monthly LAI values, respectively. When LTcon is greater than zero, a fraction of NPP is distributed evenly throughout the year. The LTvar term was computed from annual z/XLAIn, the sum of all monthly time steps for which the difference between month t and month t-1 was less than zero. For completely deciduous grid cell locations, the LTcon term is zero. For completely constant evergreen grid cells, the LTvar term drops to zero for all months. In the model experiments discussed here, it was assumed that no leaves remain in the canopy longer than one year, that the timing of leah root, and wood litterfall is synchronous, and that senescent plant material immediately enters litter pools. Litter C-to-N ratios are continually adjusted over the course of a model run in accordance with computed annual totals of available mineral nitrogen and a prescribed maximum leaf nitrogen concentration (Table 5; based on values reported by Rodin and Bazilevich [1967], Titlyanova and Bazilevich [1979], Coupland and Van Dyne [1979], Cole and Rapp [ 1981 ], Gosz [ 1981 ], Heal et al. [ 1981 ]). Mineralized N (MINn) available to vegetation (in excess of the amount needed to meet total N litter, microbial and soil pool demands; see section below) is transferred to monthly plant production (NPP) up to the amount needed to produce litter of the highest N quality (lowest possible C-to-N ratio). If the size of the MINn pool is less that the amount needed to support NPP at the highest possible litter quality level, the contents of the MINn pool is transferred to NPP, and the C-to-N ratio of NPP floats upward accordingly. Because we assumed that our NPP estimates largely account for effects of N limitation, no minimum N level is set for new production. Total N returned in litterfall is calculated according to the yearly amount of woody and nonwoody biomass produced divided by mineral N uptake over the previous 12 months. This continually updated ratio is used in combination with fixed (biome specific) lignin levels (Table 5) to partition litterfall into structural and metabolic subpools (see following sections). Soil Carbon and Nitrogen Fluxes Major soil carbon-nitrogen pools in the model and transformation processes connecting them at the TABLE 5. Litter Quality Settings for Vegetation Classes Class Minimum CN %Lignin
where LTcon(X) = [LAImin(x)/LAIave(X)]/12 (14b) and LTvar(X,t) = [6LAIn(x,t)/Z 6LAIn] [1 -{ LAImin(x)/LAIave(X) } ]