North American gross primary productivity: regional characterization and interannual variability

Seasonality and interannual variability in North American photosynthetic activity reflect potential patterns of climate variability. We simulate 24 yr (1983–2006) and evaluate regional and seasonal contribution to annual mean gross primary productivity (GPP) as well as its interannual variability. The highest productivity occurs in Mexico, the southeast United States and the Pacific Northwest. Annual variability is largest in tropical Mexico, the desert Southwest and the Midwestern corridor. We find that no single region or season consistently determines continental annual GPP anomaly. GPP variability is dependent upon soil moisture availability in low- and mid-latitudes, and temperature in the north. Soilmoisture is a better predictor than precipitation as it integrates precipitation events temporally. The springtime anomaly is the most frequent seasonal contributor to the annual GPP variability. No climate mode (i.e. ENSO, NAM) can be associated with annual or seasonal variability over the entire continent. We define a region extending from the Northeast United States through the midwest and into the southwestern United States and northernMexico that explains a significant fraction of the variability in springtime GPP. We cannot correlate this region to a single mechanism (i.e. temperature, precipitation or soil moisture) or mode of climate variability.


Introduction
Global atmospheric CO 2 concentrations have been increasing over the past 250 yr in response to anthropogenic sources in the form of human burning of fossil fuel and land cover change (Keeling et al., 1995). The net increase in CO 2 concentration represents the residual CO 2 from the anthropogenic contribution and large exchange of CO 2 between the atmosphere and surface (oceans and land). It has been shown that on an annual basis, only about 50% of the CO 2 emitted by fossil fuel burning and land cover change resides in the atmosphere; about half is taken up by the oceans and terrestrial biosphere (Oeschger et al., 1975;Tans et al., 1990;IPCC, 2007).
The annual atmospheric CO 2 increase is not linear, but has variability determined by seasonality and the spatial configuration of the continents and oceans (Tans et al., 1990). Northern hemisphere CO 2 concentration is lowest during the Boreal summer, when biospheric uptake is large. In winter when veg-etation is dormant, northern hemisphere CO 2 concentration is higher. A similar signal, with smaller magnitude due to smaller land coverage, is seen in the southern hemisphere. There is also interannual variability in the magnitude of the CO 2 increase. For example, volcanic eruptions have been shown to attenuate the rate of atmospheric CO 2 increase (Roderick et al., 2001), presumably a biospheric response to increased diffuse light (Gu et al., 2002;Niyogi et al., 2004) due to increased aerosol loading, or alternatively, a decrease in ecosystem respiration due to decreased temperature.
Increased greenhouse gas concentration in the atmosphere is predicted to modify the radiative forcing of the atmosphere (Cox et al., 2000;IPCC, 2007) which will result in changes to the meteorological forcing at the surface. Changes in surface behaviour will modify the partitioning of energy flux returned to the atmosphere, as well as imposing further changes to the radiative forcing due to changes in surface carbon flux. The atmosphere and terrestrial biosphere are tightly coupled with respect to exchange of energy, mass and momentum; predictions of future climate are critically dependent upon our understanding of processes operating under present conditions.
It is tempting to correlate atmosphere-biosphere carbon flux with climate indexes such as the Northern Annular Mode [NAM; here considered synonymous with the North Atlantic Oscillation (NAO)] or El-Niño/Southern Oscillation (ENSO) as a way to link regional carbon flux to large-scale climate dynamics. Demonstrable causal links between modes of climate variability and ecophysiological behaviour may provide insight into future climate if these modes exhibit secular trends. Russell and Wallace (2004) showed a correlation between NAM and rate of increase of global CO 2 concentration. In this case, the decreased rate of CO 2 growth was tied to greenness as measured by normalized difference vegetation index (NDVI) over the Eurasian continent, in the form in increased growing season length due to stronger warm onshore flow during the positive phase of the NAM. This finding was consistent with that of Schaefer et al. (2002), who also found that interannual variability in tropical net ecosystem exchange (NEE) of carbon was correlated to ENSO, as a response to changes in precipitation patterns. Surface carbon flux feeds back into atmospheric circulation and climate. Identifying strong mechanistic feedbacks, while not completely 'closing the loop' between surface and atmosphere, will strengthen scientific underpinnings for predictions of future climate.
To date, most studies have focused on the Northern Hemisphere, due to the large land fraction. However, with the exception of the aforementioned link between Eurasian carbon flux and the NAM, strong coupling between the terrestrial biosphere and atmosphere has defied explanation. In particular, studies have not shown consistent results in North America (NA). Zhou et al. (2001) looked at 19 yr of NDVI data and found a consistent response in Eurasia, both in greenness and length of growing season, but a more 'fragmented' situation in NA. Buermann et al. (2003) correlated greenness, temperature/precipitation and NAM/ENSO using meteorological and NDVI data for years [1982][1983][1984][1985][1986][1987][1988][1989][1990][1991][1992][1993][1994][1995][1996][1997][1998], and corroborated the strong link between vegetation and warm onshore flow in Eurasia during the positive phase of the NAM. Again, the results for NA were not as coherent, suggesting more complex interactions between meteorological forcing and surface processes.
A strong causal link between the NAM and carbon uptake has been found in Eurasia, while determination of North American ecophysiology and the forcing mechanisms that determine it have been more ambiguous. In this paper, we will focus on gross primary productivity (GPP) over the North American continent, by evaluating 24 yr of model simulations. By taking a bottom-up approach, we hope to gain insight into one aspect of ecosystem behaviour that may inform our understanding of North American biophysics and carbon dynamics. The goal of this study is multiple: Is it possible to partition the continent into regions that are dominant in terms of explaining large-scale interannual GPP variability? Are there patterns in temporal behaviour that we can identify? Can we describe regions in NA that have identifiable reliance on particular atmospheric drivers of GPP? For example, can we identify regions where an early spring is indicative of an annual increase in GPP, or regions where anomalously high midsummer precipitation results in a large positive excursion in annual GPP? Ultimately, we will attempt to correlate North American GPP variability to modes of climate variability such as NAM and ENSO. Although multiple modes of teleconnection have been identified, it has been proposed (Quadrelli and Wallace, 2004) that most, if not all of these modes represent a linear combination of these two most dominant modes. Rather than compare ecosystem behaviour with many climate modes, we prefer to start with a conservative analysis. With a model we can isolate biospheric uptake of carbon from efflux (respiration or anthropogenic sources) and atmospheric transport, and provide a basic description of biospheric response to climate dynamic processes with spatial and temporal resolution higher than the continental/annual scale often reported in inversion studies (i.e. Gurney et al., 2002Gurney et al., , 2008Rödenbeck et al., 2003;Baker et al., 2006). This paper is organized as follows: Section 2 gives a model synopsis and reviews previous results, as well as giving an overview of the statistical techniques used in the analysis. Analysis of NA GPP is contained in Section 3. This encompasses mean behaviour, regional decomposition and statistical evaluation. Summary and conclusions are given in Section 4.

Model description
The simple biosphere model (SiB) is a land-surface parameterization scheme originally used to simulate biophysical processes in climate models , but later adapted to include ecosystem metabolism Sellers et al., 1996a). SiB is a model that is useful to meteorologists for its ability to simulate exchanges of mass, energy and momentum between the atmosphere and terrestrial biosphere, and useful to ecologists for its ability to do so in a process-based framework that allows for simulation of explicit biophysical mechanisms. The parameterization of photosynthetic carbon assimilation is based on enzyme kinetics originally developed by Farquhar et al. (1980), that are linked to stomatal conductance and thence to the surface energy budget and atmospheric climate (Collatz et al., 1991Sellers et al., 1996a;Randall et al., 1996).
The soil representation is similar to that of CLM (Dai et al., 2003), with 10 soil layers and soil column depth of 10 meters. Root distribution follows Jackson et al. (1996). SiB has been updated to include prognostic calculation of temperature, moisture and trace gases in the canopy air space Vidale and Stöckli, 2003). We refer to this version of the code as SiB3.
Model photosynthesis rate is tightly coupled to total latent heat flux through transpirational losses of moisture through stomates. Photosynthesis is calculated as the minimum of ratelimitation by light, enzyme kinetics, and electron transport (Collatz et al., 1991;Sellers et al., 1992). Photosynthesis is further scaled downward from an optimum rate by limitation imposed by temperature, relative humidity and moisture availability in the soil (Sellers et al., 1992). Leaf-level temperature, humidity and internal CO 2 concentration are coupled via the Ball-Berry process (Ball et al., 1987) and solved simultaneously. Moisture availability, defined by the combination of root density and soil water concentration in individual model soil layers, imposes a fundamental constraint on photosynthesis and hence evaporation. Commonly, models have defined water availability in terms of soil depth or root density alone, which is unrealistic when compared to actual plant behaviour. We find that coupling root and reservoir characteristics  provides a more realistic simulation framework.

Model runs
For this analysis, we ran a global simulation of SiB on a 1 × 1 degree cartesian grid. Vegetation type is determined by maps as described in DeFries and Townshend (1994), and soil character is specified by IGBP (Global Soil Data Task Group, 2000). In this simulation, SiB3 uses a 10-min timestep forced with 6hourly regridded meteorological analysis products from the National Centers for Environmental Prediction (NCEP Reanalysis-2; Kalnay et al., 1996;Kanamitsu et al., 2002) interpolated to the model timestep for the years 1983-2006. SiB3 ingests temperature, pressure, precipitation, wind and radiation as forcing variables. Vegetation phenology is provided by the GIMMSg NDVI product (Brown et al., 2004;Tucker et al., 2005;Pinzon et al., 2006) which is used to calculate Leaf Area Index (LAI) and fraction of Photosynthetically Active Radiation absorbed (fPAR) following Sellers et al. (1996b).
Reanalysis products such as NCEP have known biases in precipitation (i.e. Costa and Foley, 1998) and other variables (Zhao and Running, 2006;Zhang et al., 2007a). As precipitation can be expected to have considerable influence on photosynthetic processes, we scale the NCEP precipitation to a data set that incorporates satellite and surface observations, in this case the Global Precipitation Climatology Project (GPCP; Adler et al., 2003). Using monthly precipitation values from GPCP, we scale the NCEP precipitation for consistency. We do not create precipitation events, although we may remove precipitation if the GPCP product indicates no precipitation at a location for a given month.
The model was initialized with saturated soil, and the entire 24 yr period  was simulated twice as a spinup, with the model re-initialized with ending (31 December 2006) model state at 1 January 1983. Ecosystem respiration was scaled to obtain annual carbon balance following Denning et al. (1996). Model diagnostics were output as monthly averages on the global grid, and subsampled for NA.

Model validation
Before evaluating model results for detailed analysis of ecosystem response to meteorological variability for NA, we wish to demonstrate a baseline of model performance. Continental-scale observations of surface flux do not exist so we compare SiB results with a variety of other sources as a means of establishing confidence in model output. We can appraise SiB behaviour when compared against leaf-level measurements, against observations of latent and sensible heat and carbon flux as observed by eddy covariance observation towers and against fluxes inverted from flask measurements and atmospheric transport to demonstrate model competence. These comparisons are not the focus of this paper, but can provide a foundation for trust in model simulation of less observable quantities such as regional ecophysiological variability. If we can demonstrate model skill at multiple scales, then inferences made about regional-scale ecophysiology will be more robust. Although the emphasis in this paper is on GPP, there are no direct measurements of fundamental photosynthetic uptake at either the canopy or regional scale. However, SiB is a 'third generation' surface scheme (Sellers et al., 1997) that contains self-consistency in the equation set that links radiative transfer, evaporation and ecophysiology. Internal consistency and multiple constraints make it possible to extend a degree of confidence to modelled GPP by comparing model output to observed quantities. In this regard, SiB has a demonstrable ability to simulate landsurface processes.
SiB was developed with the intended use as a lower boundary for Atmospheric General Circulation Models (AGCMs). Sato et al. (1989) describe implementation of SiB into an AGCM, and show that surface behaviour over diurnal to seasonal scales are realistic in the fully coupled model. Randall et al. (1996) describe the 'greening' of the Colorado State University AGCM when SiB phenological behaviour is upgraded from tabular values to those derived from satellite NDVI observations. Noncoupled [offline, driven by European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis products] global SiB GPP is shown to be consistent with accepted values in Zhang et al. (1996). SiB has been utilized as a lower boundary for the Regional Atmospheric Modeling System (RAMS; Cotton et al., 2003;Pielke et al., 1992) as well Nicholls et al., 2004;Wang et al., 2007;Corbin et al., 2008). Denning et al. (2003) and Nicholls et al. (2004) describe how biological processes couple with meteorological process to define the regional carbon budget over Wisconsin, USA. Wang et al. (2007) and Corbin et al. (2008) extend the analysis to the central North American continent.
From its inception in 1986, SiB (versions 1, 2 and 3) has also been evaluated by confronting the model with local observations. Sellers and Dorman (1986) compare modelled energy flux, stomatal resistance, albedo and leaf water potential to observations over periods from several days to 1 month at multiple sites. Colello et al. (1998) confirm that SiB can reproduce diurnal cycles of energy, moisture and carbon flux at a grassland over periods of several days using FIFE data, and Hanan et al. (2005) show that incorporation of heterogeneous C3/C4 physiology maps into SiB improve modelled annual cycles and interannual variability at a grassland site in Oklahoma. Effective radiative temperature, soil wetness and energy fluxes at a Tibetan prairie are successfully reproduced during growing season months in the study of Gao et al. (2004). Baker et al. (2003) and Schaefer et al. (2008) compare SiB energy and carbon fluxes to eddy covariance observations taken at mid-latitude forest sites: diurnal, monthly and interannual variability are captured by the model, although SiB is biased towards slightly high Bowen ratio when canopy cover is very dense. The mechanisms that control seasonal variability at a site in tropical Amazonia are incorporated into SiB by Baker et al. (2008), with the result that previously out-of-phase simulated annual flux cycles at this site are brought into agreement with observations. Regional carbon exchange is not a directly measurable quantity. Flask observations of CO 2 are inverted with transport data by Gurney et al. (2008) to obtain estimates of carbon flux for the globe when partitioned into 22 oceanic and terrestrial regions. Figure 1 shows SiB NEE for years 1984-2004 (black line, symbols) compared to the results of Gurney et al. (2008) for Boreal and Temperate NA. There are eight global 'networks' for the inversion data, as changing numbers of flask locations makes a single estimate unreasonable. It is important to note that there is frequently disagreement between the inverted fluxes from different networks (i.e. 19971988-1989 in Temperate NA), underscoring the fact that absolute measurement of regional carbon flux does not exist. When our modelled flux is compared against the inversion results, SiB captures most, but not all, of the major variability seen in the inversion results. In Boreal NA, the efflux (positive values) event of 1990-1992 is well represented by SiB, as is the uptake event of 1997-1998, although that event is only captured by two of the inversion networks. The relative maximum of 2001-2002 is consistent between SiB and the inversion results, although the SiB peak is higher in magnitude. In temperate NA, general trends are consistent between SiB and inversion results. Local minima in 1991Local minima in -1992Local minima in and 1996Local minima in -1998Local minima in are consistent, as are maxima in 1994Local minima in , 2000Local minima in and 2002 It is important to reiterate that inversion results are not conclusive, as represented by the disagreement between different inversion networks. Flask coverage is sparse (especially in Boreal NA), and transport was interannually uniform in the inversion exercise, which can be a cause for ambiguity (Gurney et al., 2008). It is crucial to acknowledge that there is uncertainty in regional flux estimates from either models (such as SiB) or inversions. However, we believe that SiB fluxes, when compared to inversion products, are reasonable and demonstrate an ability to capture the larger features of ecosystem variability across broad spatial domains. This ability is critical to the results we present later in this paper.

Statistical tools
A model is merely a tool to assist in our understanding of a particular system, in this case the spatiotemporal behaviour of the terrestrial biosphere. We acknowledge that uncertainty exists in model output due to parameterizations, subgrid scale heterogeneity and reanalyzing meteorological observations into forcing data sets. For these reasons, we keep our statistical toolbox small. With the inherent confidence bounds in the model system, we believe that relationships that cannot be seen with simple tools may simply be an artefact of the mathematical manipulations being applied. Therefore, we attempt to limit statistical tools to linear regression, explained variance and correlation coefficient. To test for significance, we use a two-ended Student's t-test at the 90% level. We take a conservative approach and assume that each year comprises an independent sample, giving us 24 samples in the modelling period (for a two-tailed test we subtract two, giving 22 degrees of freedom). At this level, the probability that a particular value is the result of random chance requires a t-statistic value of t < 1.72.
We also use singular value decomposition to obtain the principle component (PC) time series, which we can regress GPP and meteorological anomalies onto to determine if there are patterns that explain large fractions of continental-scale variability. Eigenvalues are tested for significance following North et al. (1982). We found that the significance of the eigenvalues was critically dependent upon the degrees of freedom specified, which we calculated (following Bretherton et al., 1998) as where N * is the degrees of freedom, r 2 is the lag-one autocorrelation and t is the timestep of the data. Intuitively, we can think of the spatiotemporal persistence of anomalies in GPP as the dependence of an anomaly upon the previous month, season or year's environmental conditions. Plots of lag-one autocorrelation of monthly GPP (not shown) show a spatial variability in this dependence: In the northeast, in an area bounded approximately by 50 • north latitude and 90 • west longitude, there is very little variability explained by the lag-one autocorrelation, less than 20% on even a monthly basis. In the more arid regions of the continent (desert southwest, arctic), the lag-one autocorrelation on a monthly basis is much higher. When annual anomalies are used, the region of spatial independence becomes larger still, with only the desert southwest suggesting that GPP in a given year is dependent upon conditions in the previous year. For this reason, we continue to use 22 as our degrees of freedom, which is not only smaller than the area-mean value found using eq. (1), but consistent with our theme of using conservative statistics throughout the analysis.

Mean North American GPP
For this analysis, we define NA as the region northward of 15 • north latitude between 50 • and 170 • west longitude. The spatial distribution of simulated annual mean GPP is shown in Fig. 2 (panel a). Maximum uptake of CO 2 is in the tropical forest of southern Mexico, over 3.5 kg m −2 yr −1 in some areas. The southeastern United States and Pacific Northwest are the next largest in terms of annual assimilation of carbon, at approximately 2 kg m −2 yr −1 . The desert southwest and arctic tundra regions show smallest annual GPP, which is not unexpected given the harshness of the climate and lack of biomass in these regions.
The standard deviation of annual GPP (Fig. 2, panel b) shows that the tropical forest in southern Mexico has large interannual variability, but the productive regions in the southeast and Pacific Northwest do not. GPP in the Monsoon region in Mexico is highly variable, as is the California coast south of San Francisco Bay. There is also a belt of relatively large variability centred on 100 • west longitude, from Southern Texas to the prairie provinces of Canada. This midwestern band of large standard deviation deserves attention because it exists as a natural boundary between the productive east and the relatively dry intermountain west. This region encompasses the 'dry line' that often focuses severe weather, as well as the region effected by the cold, fast-moving airmasses ('Alberta Clippers') that intrude into the mid-latitude United States from Canada. It is intuitive that this region will not only experience large variability in ecosystem function, but that overall variability here has the potential to impact continental-scale carbon characteristics.
The coefficient of variation (defined as the standard deviation divided by the mean; not shown) for GPP is largest in the desert southwest, and smallest in the local maxima GPP regions in the southeast United States and Pacific Northwest, as well as in the Boreal Forest of Canada. In general, large coefficient of variation is found where mean annual GPP is low.
Any analysis of North American GPP must include consideration of seasonality. In the extreme north, cold winters and brief warm summers ensure that annual GPP is completed in only a few months. For example, in Barrow AK there are only 3 months (June, July and August) that have mean monthly temperature above 273 K, with mean annual temperature amplitude of almost 50 K. As one moves south, temperature seasonality is damped and mean annual temperature is larger. The mean monthly temperature for locations such as Miami FL and Mexico City is above freezing in all months, and the amplitude of the annual cycle is 10 K or less. There are periodic cold air intrusions southwards (more so in the SE United States than in Mexico or Central America), but mean conditions are suitable for photosynthetic activity throughout the year. Precipitation seasonality plays a role as well, especially in the southern regions where temperature is not as variable. Figure 3 (panel a) shows the month where mean maximum GPP occurs, and Fig. 3 (panel b) shows the fraction of annual GPP that occurs during the month of greatest activity. In the arctic, maximum monthly GPP happens in July or August, and 1 month can comprise upwards of 40% of the annual total. By contrast, in the southern United States (with the exception of the southern Rockies), maximum photosynthetic activity occurs in April or May, and no individual month contributes more than 15% to the annual total. In the monsoon region of Mexico, maximum GPP occurs in August/ September, after seasonal rains have replenished moisture in the soil.

North American GPP variability
Continents are commonly partitioned by vegetation type (i.e. Peters et al., 2007) as a means to determine more detailed relationships. In this study, we prefer a regional discretization, based roughly on spatial coherence in the mean, standard deviation and coefficient of variability of GPP. These regions are shown in Fig. 4, and while we did not explicitly establish criteria to define continental subregions, these regions implicitly incorporate natural delineations of mean annual GPP, standard deviation and coefficient of variability, vegetation type, topography and climatological variables such as annual mean precipitation and temperature. These regions are listed in Table 2, along with the fraction of NA land area, fraction of mean annual NA GPP occurring in the region and the ratio of GPP fraction to land fraction. Several features of Table 1 are worth noting.
• The SouthWest and SouthEast regions have identical area, but the SouthEast region has over four times the fraction of mean annual NA GPP that the SouthWest region does.
• Over half (54%) of the mean annual NA GPP occurs in three regions: SouthEast, NorthEast and Mexico/CA. However, the variability in the SouthEast and NorthEast regions is small (Fig. 2, panel b), which may reduce the impact these regions impose on continental-scale GPP anomaly.
• The MidWest region contributes around 12% to mean annual NA GPP, but the relatively large standard deviation in annual GPP for the region suggests that this region may play a larger role in continental-scale behaviour.  The annual anomaly in modelled GPP for 1983-2006 for the entire continent is shown in Fig. 5, with the regional anomalies superimposed. In this figure, the continental-scale anomaly is plotted as the solid line and the regional anomalies are sorted by magnitude, with smallest variation from the regional mean plotted closest to the zero line; larger anomalies are plotted further from the origin, demonstrating which regions contribute the most to that year's anomaly on the continental scale. In no year do all regional anomalies have identical sign; both positive and negative anomalies exist on a regional basis in all years. In 4 years (1996,2000,2002,2005), all regions save one have similar sign. Interestingly, the southernmost region is the outlier for each of these years. For years with small overall anomaly (1990,2001,2004), there are relatively large, yet compensating, regional excursions from the mean.
Variability about the mean NA GPP is not consistently dependent on any single region. Of the nine regions, only eastern Boreal Canada (Boreal-East) and the Pacific Northwest (NorthWest) are never the largest anomaly for a given year, where all other regions contribute the largest anomaly to the annual variability at least once. Several regions stand out due to their frequent large excursions from the mean; Mexico/Central America, MidWest and Central/Western Boreal regions. The Boreal regions commonly exhibit similar sign (i.e. 1983, 1985, 1992) to the overall anomaly but not exclusively (2001,2004) suggesting that these adjoining regions may respond to different forcing mechanisms that determine annual GPP variability.
The SouthWest is the largest anomaly three times; even though the mean GPP is small there, the coefficient of variability is the largest on the continent. The SouthWest can be an influence for either positive or negative anomalies-it is not just that it can contribute if it rains, although the largest contribution is during the large El Niño event of the early 1980s. This region is important for continental-scale carbon flux.
The continental-scale annual GPP anomaly is partitioned by seasonal contribution in Fig. 6. On a seasonal basis, all seasons except winter have at least 6 yr where that seasonal anomaly is the same sign as the annual anomaly and is largest during the year. Spring anomalies dominate the year the most (11 times) with Summer contributing the most to the annual anomaly 7 yr, and Fall 6. The two largest anomalies during the year are usually adjoining seasons (i.e. spring and summer, or summer and fall), but not exclusively. It is not uncommon for the two largest seasonal anomalies with the same sign to be opposing seasons (i.e. 1985, 2005), nor is it uncommon for adjoining seasons to have opposing signs (i.e. 1991, 2004) although in this Fig. 5. Regional contribution to annual North America GPP anomaly, in GT carbon, for years 1983-2006. Continental-scale variability is shown as the solid black line, and regional contribution to the total is given by the coloured boxes. Smallest regional anomalies are plotted nearest the zero line, and largest anomalies further from the x-axis. Fig. 6. Seasonal contribution to annual North America GPP anomaly, in GT carbon, for years 1983-2006. Continental-scale variability is shown as the solid black line, and seasonal contribution to the total is given by the coloured boxes. Smallest regional anomalies are plotted nearest the zero line, and largest anomalies further from the x-axis.
case the magnitude of the anomaly is generally small and winter is usually one of the seasons. For each season, the distribution between positive and negative is fairly evenly distributed; both positive and negative excursions from the mean are found in the simulation record.

Correlation to physical mechanism
3.3.1. Annual GPP variability. Ultimately, we wish to link ecophysiological behaviour to modes of climate variability, but doing so directly requires a causal jump across the meteoro-logical mechanisms that influence ecosystem behaviour. In this section we make simple statistical comparisons of GPP to meteorological forcing such as temperature, radiation, and precipitation. We also include a comparison to soil moisture availability, as precipitation alone may neglect features of precipitation distribution that may be misleading. For example, a large precipitation event will suggest a positive anomaly, but if a large fraction of that precipitation is lost as runoff then the potential benefit to vegetation will be lost. Furthermore, the effect of wintertime precipitation anomalies may be attenuated or exacerbated by the nature of the spring warmup. Using soil moisture as a diagnostic tool takes advantage of its integrative nature, which may not be possible if precipitation alone is used.
We regress variability in annual GPP against variability in annual precipitation, soil moisture availability, temperature, and radiation. The fraction of GPP variability explained by the various mechanisms is calculated along with a Pearson's correlation coefficient to help determine the nature of the relationship. Significance is determined using a Student's t-statistic, calculated for 90% significance assuming 22 degrees of freedom (assuming each year as an independent sample).
For the North American continent, we aggregate the results and show the mechanism that explains the largest fraction in annual GPP in Fig. 7 (panel a). Figure 7 (panel b) shows the amount of variance explained by the mechanism that explains the most variance. We plot the mechanisms as grid-boxes (not contoured) because the classification used is discrete. Blank spots reflect gridcells where no single mechanism is able to explain interannual variability in GPP at the 90% significance level. Soil moisture availability explains the most GPP variability for a significant fraction of the continent, from Mexico and Central America through the southern tier of the United States. The midwestern United States and the southern boundary of the Canadian Prairie provinces are responsive to soil moisture availability, as is most of Alaska and parts of Arctic Canada. The fraction of variance explained by soil moisture availability is largest in the Yucatan Peninsula, northwest Mexico/Southwest United States and the Colorado Plateau. In this region, variance explained can exceed 95%. In the southern plains, variance explained is in the 20-40% range, while there is a local maximum in soil moisture influence over the Dakotas and southern Manitoba of around 75%. In the arctic, the variance explained is lower, generally between 20% and 40% (with small local maxima).
Temperature explains most of the GPP variability in the eastcentral United States and over most of Canada. Temperature influence in Canada is intuitive, as an early spring can be expected to lead to anomalously large GPP. However, the fraction of variability explained by temperature is small everywhere where temperature is the dominant mechanism, generally between 10% and 40%. Radiation variability explains the most variance over small pockets in the Pacific Northwest and Ohio River Valley. The fraction of GPP variability explained by radiation is generally small.
3.3.2. Spring GPP variability. We identified spring as the season that most frequently makes the largest contribution to annual GPP variability (Fig. 6). We can repeat the analysis that we performed on annual anomalies for springtime GPP, and in addition to regressing springtime GPP onto springtime meteorology, we can regress springtime GPP onto winter meteorology to look for a lagged response.
When focus is isolated on spring, the intuitive result is borne out from the statistical analysis: Spring GPP responds positively to anomalously warm temperature across a large fraction of NA. The largest response, in terms of fraction of spring GPP variability explained by temperature variability, extends along a band from New England, along the United States-Canada border into the Pacific Northwest. There is a narrow strip through central Mexico where 20-50% of the springtime GPP variability is explained by temperature, yet in this region the relationship is inverted-cooler spring enhances productivity. The annual pattern shown previously (Fig. 5) in NW Mexico/SW United States holds in the spring as well-increased soil moisture availability is strongly correlated with enhanced GPP.
When we regress springtime GPP against winter meteorological variability, the same relationships seen during the same season (spring/spring) regression hold, albeit with slightly different spatial structure. This is likely due to the fact that anomalously large winter precipitation has an opportunity to be integrated into the soil prior to springtime GPP onset for large areas of the continent. When the seasonal lag is included into the analysis, anomalous precipitation in winter or spring is similar for the purpose of statistical regression and the patterns we retrieve. The positive temperature relationship to GPP is significant over a much smaller area, centred on Northwestern Quebec.

Modes of climate variability
Ultimately, we wish to link modes of large-scale climate variability (MCV) such as ENSO or NAM to continental-scale photosynthetic behaviour. However, partitioning continental GPP into regional or temporal components (Figs 5 and 6) does not reveal consistent or coherent patterns, suggesting that ecosystem linkage to large-scale climatic modes may be a difficult prospect.
However, the lack of continental-scale coherence does not preclude regional relationships between modes of climate variability and ecosystem function. To test this idea, we regress GPP anomalies against the two main indices that measure modes of climate variability: the Multivariate ENSO Index (MEI: Wolter and Timlin, 1993;Wolter and Timlin, 1998), and the Northern Annular Mode (NAM; Thompson and Wallace, 2000). We regress anomalies of GPP, temperature and precipitation against these indices as a way to draw out mechanisms that influence vegetation behaviour. As we did for meteorological mechanisms, we calculate a Pearson's correlation coefficient, and test for significance using a two-tail Student's t-statistic at the 90% level.
The data may be dissected in several ways. Modes of climate variability show power over multiple timescales (Rasmussen and Carpenter, 1982;Barnston and Livezy, 1987;Enfield and Birkes, 1993;Hurrell, 1995), and many index values are only evaluated during winter/spring months. We might expect to see a coupling between annual GPP anomaly and annual MCV index in some cases, and a higher-frequency response in others. We might also expect to see a lag between the mechanism invoked by variability of a climate index and environmental response. Furthermore, we focus seasonal attention on the spring (MAM) and summer (JJA) seasons, as these are the seasons with largest GPP and variability, and thus the seasons most likely to influence annual carbon uptake anomaly. Therefore, we conducted the following regressions: • Annual GPP anomaly onto annual index: Can we relate annual GPP anomaly to low-frequency (annual) variability in a climate index?
• Annual GPP anomaly onto winter, spring, or summer index: Does seasonal-scale variability in a climate index translate to an annual anomaly in GPP?
• Spring GPP anomaly onto winter index, or summer GPP onto spring index (lag-one seasonal comparison): Is there a delay in atmospheric or ecosystem response to climatic forcing?
• Spring GPP anomaly onto spring index or summer onto summer (lag-zero): Immediate ecosystem response to changes in forcing invoked by index variability.
Results are summarized later; no obvious relationship between a single MCV and GPP variability emerged on the largescale, although there are multiple responses on the regional scale between climate indices, GPP, and the mechanisms that influence GPP that are intuitive and consistent with current knowledge of manifestation of climate variability on meteorology.
3.4.1. Multivariate ENSO Index: MEI. For example, we have shown a relationship between soil moisture availability and GPP in the desert southwest (i.e. Fig. 7, panel a). There is a wellknown link between positive ENSO index and high winter precipitation in this region (Sheppard et al., 2002;Mauget, 2003), although an inverse relationship is seen during the monsoon season of the late summer (Higgins et al., 1999). Therefore, we might expect a correlation between GPP and ENSO index here. We do not see large areas of correlation between annual GPP and annual MEI, but we do see a significant relationship in the desert southwest and Pacific Northwest from winter through spring. When spring (MAM) GPP is regressed against the winter MEI, we see between 20% and 30% of the variability in the desert southwest explained, and up to 45% of the GPP variability in Washington, British Columbia and western Alberta explained. This relationship is reinforced when spring GPP is regressed against spring MEI. By summer, the influence in the southwest has moved slightly northeast into the Colorado Plateau, and the correlation in the Pacific Northwest is gone. In the desert southwest, the ENSO influence is in the form of enhanced precipitation, and the positive correlation between MEI and GPP to the north is due to warmer temperatures resulting in extended growing season.
3.4.2. Northern annular mode. The characteristic frequency of the NAM is much shorter than that of ENSO, so we may expect that correlation between annual mean index and annual mean GPP are non-existent. However, we do see a suppression of GPP along the Canadian arctic coast between 100 and 120 west longitude, although it explains less than 25% of the variability. There is no statistical significance between NAM and temperature or precipitation in this region, so the exact mechanism by which NAM influences GPP is not clear.
We see an influence on annual temperature anomaly by seasonal NAM index. A high springtime NAM is associated with positive annual temperature anomaly in a region from central Alaska through the Yukon-NW Territories border, where the NAM explains up to 30% of the temperature variability. Similarly, a high summertime NAM can explain up to 40% of the variability in annual temperature over a sizable fraction of Wyoming and eastern Montana. However, in neither of these regions do we see a correlation between GPP and either spring or summer NAM. This suggests that other mechanisms (water availability, relative humidity) play a larger role in regulating GPP in these regions than temperature does.
One-season lagged comparison shows no significance when spring GPP is compared to winter NAM. This is likely due to low ecosystem activity in the region of influence (generally far north latitudes) until later in the year. When we correlate summer ecosystem activity to spring NAM index, we see that high NAM is associated with suppressed GPP along the Canadian arctic coast noted earlier.
When NAM is compared to same-season quantities, we see an arc of anomalously high GPP along an arc from central Alaska along the eastern slope of the Canadian Rockies to the Canada-USA border. This GPP anomaly is associated with a positive temperature anomaly, implying early warming and extension of growing season.
The net result is that no continental-scale relationship emerges, with respect to temperature, precipitation or GPP. This is consistent with the previously cited studies such as Zhou et al. (2001) or Buermann et al. (2003). The extensive continental dependence on soil moisture suggests that subtle interactions between precipitation, temperature and radiation are responsible for large-scale GPP variability, and these interactions are heterogeneous in space and time. We are not, at this time, able to make predictions about NA ecophysiological behaviour based on ENSO or NAM index.

EOF/PC analysis
To determine if coherent patterns exist for GPP variability on annual and seasonal scales independent of reliance on physical mechanism, we performed singular value decomposition (SVD) analysis to obtain the PC time series for GPP anomalies. We tested eigenvalues for significance following North et al. (1982).
On an annual basis, we found no significant pattern, nor did one emerge for summertime GPP. However, the first eigenvalue for springtime showed separation, suggesting that a spatially consistent pattern of GPP variability exists for this season. The GPP anomaly regressed upon the first PC time series is shown in Fig. 8, and shows a spatially coherent region extending from around 90 • to 105 • west longitude, bounded by the Rio Grande in the south and extending into the Prairie Provinces of Canada to the north. This region branches northeastward up the Ohio River valley towards New England. There is also an associated region of coherence in the Monsoon region of Mexico and the southwestern United States. This EOF implies that there is a large region that behaves consistently, and explains most of   . 9. Regional contribution to springtime North America GPP anomaly, in GT carbon for years . Continental-scale variability is shown as the solid black line, and regional contribution to the total is given by the coloured boxes. Smallest regional anomalies are plotted nearest the zero line, and largest anomalies further from the x-axis.
the variability in springtime GPP. If we look at the subregions that describe this region (Fig. 4), it falls mainly in the Mid-Western, NorthEast and Boreal-Central regions. The spring GPP anomalies, broken out by subregion contribution, are shown in Fig. 9, and it is easily seen that the three aforementioned regions play a large role in determining the seasonal variability. The largest spring anomaly (of the same sign as the continental-scale seasonal anomaly) is from the MidWest, NorthEast or Boreal-Central region in 16 of 24 yr, one of the top two in 20 of 24 yr, and at least one of these regions is in the top three anomalies in 23 of 24 yr. In 4 yr (1983,1987,1997,2002), these three regions all have the same sign on their anomaly, and are the top three anomalies in terms of the magnitude of GPP variability that determines continental-scale GPP variability in spring.
Temperature is the driving mechanism for a large portion of the springtime variability in this region. Generally, an early spring results in increased GPP. However, there are also portions of Texas and the Dakotas where, even in spring, the available soil moisture explains most of the variability when compared with other mechanisms.
Springtime GPP variability, when correlated with modes of climate variability, show no strong relationships over the region shown in Fig. 8 when compared concurrently or with a 1-season lag. ENSO index is positively correlated with warmer temperatures along the southern edge of the Great Lakes and in Alaska/Yukon when springtime temperature is compared with winter index. The northern area moves eastward and expands in size when spring temperature is correlated to spring ENSO index. The correlation to GPP is smaller, partially due to low spring GPP in the north. There is virtually no correlation between the NAM and either temperature or GPP at a lag of one season (spring GPP/temperature correlated with winter NAM index). When concurrent NAM index and temperature/GPP are compared in spring, there is a large band extending from Alaska through the Prairie Provinces where temperature is positively correlated to the NAM. The GPP correlation is strongest for a region extending from Oklahoma through Louisiana. Therefore, although there are portions of the pattern shown in Fig. 8 where GPP anomaly can be tied a particular climate index, no one index or mechanism explains the region-wide behaviour.

Comparison with observations
Rigorous comparison of our results with observational data is not possible, but we can compare our results from Fig. 7 with studies from the literature. In general, we find that grasslands are most sensitive to soil moisture as a mechanism that influences GPP. This is consistent with the results of Flanagan et al. (2002), who looked at 3 yr of meteorological and eddy covariance flux data from Lethbridge (49 • N/112 • W). Meyers (2001) also show strong dependence between soil moisture and GPP in an Oklahoma grassland (35 • N/98 • W) when evaluating multiple years of observations. Kjelgaard et al. (2008) report a strong correlation between precipitation and GPP in the Texas 'Hill Country' of the southern plains (30-32 • N/100-105 • W), which is inconsistent with our findings. However, in this case the locally shallow soil (<50 cm) may reduce the importance of soil moisture storage. SiB uses globally uniform soil depth, and so cannot reproduce this result at this time. Zhang et al. (2007b) compared a piecewise regression model of GPP to flux tower data and to GPP modelled using MODIS over the high plains region of the United States. Although soil moisture was not a component of the model, the authors allude to the importance of moisture holding capacity to GPP in the region.
We find that GPP variability in the northeastern quadrant of NA is principally dependent on temperature. Conceptually, this can be thought of as a lengthening of the growing season correlating to high seasonal or annual GPP. This result is consistent with that of Richardson et al. (2009), who found growing season length was a significant driver of GPP variability at two flux tower sites in the northeastern United States. Urbanski et al. (2007) found temperature to be a driver of GPP variability at Harvard Forest in Massachusetts. During the 13 yr of their study, the largest annual GPP occurred during a year with extremely early canopy development. The year with the smallest annual GPP was also a year with an early canopy. However, in this case low temperature and insolation in early summer acted to suppress GPP. They found that soil moisture was unrelated to seasonal or annual GPP variability, except during the late summer. This suggests that temperature, while an important driver of GPP variability, is not the sole influencing factor. This is supported by Fig. 7, which shows that while temperature is the dominant driver of GPP variability the fraction of variability explained is relatively low.
Three different-aged Douglas-fir stands on Vancouver Island (British Columbia, Canada) were studied by Jassal et al. (2009). They found that summertime GPP was water limited. However, they also noted that winter season GPP is energy (or light) limited, and the largest annual GPP occurred during the warmest year of their observational record. This variability in the Pacific Northwest is not wholly inconsistent with our results. We find a heterogeneous situation in the region (Fig. 7, panel a), with pixels showing highest GPP dependence on soil moisture, temperature and radiation in close proximity to each other. These results, although not conclusive, indicate a general correspondence between several observational studies and our results as indicated in Fig. 7. We have not found observational studies that directly contradict our analysis of the mechanisms that drive variability in GPP.

Summary and Conclusions
As global atmospheric CO 2 levels rise, the exchange between the atmosphere and terrestrial biosphere is important for two reasons: First, around half of the anthropogenic CO 2 currently emitted is taken up by the oceans and land (the 'missing sink'), and secondly, vegetation behaviour plays an important role in determining the exchange of energy, mass and momentum between the atmosphere and land surface. Understanding of present-day ecophysiological behaviour is critical to predictions of future climate.
There has been an observed increase in the positive phase of the NAM (or NAO; Hurrell, 1995;Hurrell et al., 2001) in recent years, which has been positively correlated to increased carbon uptake in Eurasia (Schaefer et al., 2002;Buermann et al., 2003;Russell and Wallace, 2004). The causal links are simple to follow. A tightening of the polar vortex results in anomalous intrusions of warm maritime air onto the continent, resulting in an extension of growing season due to warmer spring and/or fall. If we have a scientific basis for predicting a continued persistence in the positive phase of the NAM, we can make predictions of ecophysiological response with considerable confidence.
In NA, the situation is much less clear. We have not been able to find a consistent, continent-wide vegetation response to either ENSO or NAM, and this result is consistent with previous work (i.e. Zhou et al., 2001;Schaefer et al., 2002;Russell and Wallace, 2004). Part of the problem is the large latitudinal extent of NA, further complicated by the presence of large mountain chains running along almost the entire western continental boundary from north to south. These mountains influence and disrupt circulation, and have a large impact on weather and climate. Secondly, while the situation in Eurasia can be easily attributed to temperature, in NA temperature and precipitation are secondary in importance to soil moisture availability. This can be seen in Fig. 10, which shows the lag-correlated coefficient of determination (R 2 ), calculated as where R 2 is the coefficient of determination, X is the meteorogical forcing anomaly and G is the GPP anomaly. In eq. (2), the overbar represents the mean of all gridcells on the continent (weighted by cosine of latitude), and covariance is calculated for multiple lags (zero, one, two, etc. months to the end of the year). With these plots, we can determine what fraction of the variability of a particular month's GPP anomaly is explained by the variability in meteorological forcing during the same or a prior month. To read Fig. 10, locate the desired month on the y-axis and move horizontally to the diagonal line. The value or colour at this point represents the lag-zero relationship between GPP and the chosen meteorological variable. Moving to the right, the explained variability in GPP at increasing lag is shown. The vertical axis shows the month of the forcing anomaly, and the x-axis shows the month of GPP response.
There are several interesting features of Fig. 10. First, the lack of variability explained solely by temperature is dramatic. A slight signal in the spring and fall can be seen, with spring having a longer influence than fall. There is a slight negative influence of temperature in midsummer, which is intuitive. A hot summer can impose stress on vegetation. Precipitation has more overall influence than temperature, mainly during the summer months when ecosystems can respond quickly to precipitation events. Radiation has almost no influence, on a continental basis.
The idea that soil moisture availability has the most power to explain North American GPP variability is supported by the lower right panel in Fig. 10. The fraction of variability explained is much higher that any other mechanism, and the lag covariance has influence for a much longer period of time. Up to 25% of Fig. 10. Lag covariance (divided by total variance) of GPP to anomalies in temperature (upper left), precipitation (upper right), radiation (lower left) and soil moisture availability (lower right). Axes represent time, shading shows increasing covariance from light to dark. Lag covariance is found by following the position of a month on the vertical axis towards the right along the dashed horizontal line. Numerical value represents fraction of total variance for GPP in a given month explained by variance in forcing mechanism.
August GPP variability can be explained by January/February soil moisture anomaly. This reinforces the idea that we cannot point to a single meteorological driver to explain vegetation response on the continent. Soil moisture availability is defined by unique combinations of precipitation and temperature, and is necessarily integrated through time with respect to both snow and infiltration rate.
When continental GPP is partitioned temporally, spring is the season that most frequently contributes more to the annual anomaly than any other season. However, for a given region, the annual variability is not confined to spring-it is common to see almost any season except winter contain the largest fraction of the annual anomaly for a region. Springtime meteorological variability, generally in the form of anomalously warm (cool) temperature, is correlated with anomalously high (low) GPP over large areas north of the 40th parallel. The desert southwest, not surprisingly, is tightly coupled with precipitation and soil moisture availability on both an annual and seasonal basis.
It is intuitive that spring is the season that most commonly determines the annual GPP anomaly. One would expect that an early or late spring bud-burst and leaf-out would have an impact on annual carbon budget. EOF analysis reveals a coherent region of the continent, extending along 95 • degrees west longitude from the Gulf of Mexico through the prairie provinces of Canada and extending into New England, that exhibits con-sistent springtime variability. This region commonly determines the sign of the continental GPP anomaly for the spring. As spring is the largest seasonal anomaly for almost half of the simulated years, this region is the closest we can find to a bellwether for continental-scale GPP variability. This region is not tightly coupled to any single measurement of climate variability.
At present, we cannot make statements about North American carbon uptake based on the values of climate indices. For example, Hurrell (1995) reports on a persistent elevation of the NAO index during the 1980s, but we see no corresponding response in either regional-or continental-scale GPP. A high positive-phase ENSO may suggest anomalously large GPP in the southwestern United States in the spring, or a tightening of the polar vortex may imply an early spring in northern Canada, but neither effect is large enough to have implications on continental-scale carbon flux on a seasonal or annual basis. The pattern in coherent springtime GPP variability that we find must be explained by another method.

Acknowledgments
This research was sponsored by the National Science Foundation Science and Technology Center for Multi-Scale Modeling of Atmospheric Processes, managed by Colorado State University under cooperative agreement No. ATM-0425247. This research was also funded by Department of Commerce/National Oceanic and Atmospheric Administration contract NA08AR4320893, NASA contracts NNX06AC75G and NNX08AM56G, Department of Energy contract DE-FG02-06ER64317 and NICCR contract MTU050516Z14. The authors would like to acknowledge D.W.J. Thompson and A.H. Butler for assistance with statistical analysis, and M.D. Branson for graphical support.