Plant Traits Help Explain the Tight Relationship between Vegetation Indices and Gross Primary Production

: Remotely-sensed Vegetation Indices (VIs) are often tightly correlated with terrestrial ecosystem CO 2 uptake (Gross Primary Production or GPP). These correlations have been exploited to infer GPP at local to global scales and over half-hour to decadal periods, though the underlying mechanisms remain incompletely understood. We used satellite remote sensing and eddy covariance observations at 10 sites across a California climate gradient to explore the relationships between GPP, the Enhanced Vegetation Index (EVI), the Normalized Di ﬀ erence Vegetation Index (NDVI), and the Near InfraRed Vegetation (NIR v ) index. EVI and NIR v were linearly correlated with GPP across both space and time, whereas the relationship between NDVI and GPP was less general. We explored these interactions using radiative transfer and GPP models forced with in-situ plant trait and soil reﬂectance observations. GPP ultimately reﬂects the product of Leaf Area Index (LAI) and leaf level CO 2 uptake (A leaf ); a VI that is sensitive mainly to LAI will lack generality across ecosystems that di ﬀ er in A leaf . EVI and NIR v showed a strong, multiplicative sensitivity to LAI and Leaf Mass per Area (LMA). LMA was correlated with A leaf , and EVI and NIR v consequently mimic GPP’s multiplicative sensitivity to LAI and A leaf , as mediated by LMA. NDVI was most sensitive to LAI, and was relatively insensitive to leaf properties over realistic conditions; NDVI lacked EVI and NIR v ’s sensitivity to both LAI and A leaf . These ﬁndings carry implications for understanding the limitations of current VIs for predicting GPP, and also for devising strategies to improve predictions of GPP.


Introduction
Accurate estimates of Gross Primary Production (GPP, i.e., an ecosystem's ground-area based photosynthetic rate) are needed to better assess ecosystem function and stress, as well as the role of terrestrial ecosystems in the global carbon cycle [1]. The eddy covariance technique provides a relatively direct measure of GPP at scales of a few hundred meters, but the cost and time required to operate eddy covariance sites precludes the deployment of the large network of homogenously-instrumented towers needed to accurately sample global GPP.
Remote sensing provides a cost-effective strategy to extrapolate the GPP observations from eddy covariance towers to larger spatial and longer temporal scales. Remote sensing approaches for GPP extrapolation range from simple empirical relationships to complex models [2,3]. Recent developments using Solar-Induced Fluorescence (SIF) provide a promising approach for extrapolating GPP [4], but the satellite records required for SIF retrieval are brief and/or at coarse spatial resolution (a few years at 10 km or larger; longer at coarser resolution). Vegetation Indices (VIs) provide an alternative, simple where Blue, Red, and NIR refer to the reflectance in the blue, red, and Near InfraRed (NIR) wavelengths. Currently, there is no consensus on which of these VIs provides the most universal predictions of GPP, with some studies showing stronger relationships with EVI and NIR v [5][6][7]18] while others suggest that no single vegetation index performs best in all cases and that relationships are site-specific [8,10]. Likewise, understanding of the causal links between these VIs and GPP remains incomplete, and tends to emphasize the importance of LAI, with less attention to a VI's sensitivity to biophysical attributes related to leaf level metabolism. NDVI was one of the firsts VIs developed and remains widely used. Previous studies have identified the main plant traits driving VIs using observational data or Global Sensitivity Analyses (GSA) of radiative transfer models. Both data-based and GSA studies have shown that NDVI is particularly sensitive to LAI and chlorophyll density; NDVI saturates at high LAI and is also influenced by soil reflectance [19,20]. Given these sensitivities, some degree of correlation between NDVI and GPP is expected because LAI and leaf chlorophyll largely determine the absorption of photosynthetically active radiation (PAR) radiation in canopies, and also because leaf chlorophyll is a driver of leaf-level photosynthesis [21]. Despite these potentially strong physiological linkages to GPP, NDVI has proven to be a poor predictor of GPP in some cases [10,11]. While some studies assume this issue may be related to NDVI's saturation to LAI or undesired influence of environmental conditions [18], it is unclear whether this could be related to traits not captured by NDVI.
EVI and NIR v were subsequently developed with the goals of decreasing saturation in dense canopies and reducing cross sensitivity to soil and atmospheric conditions [6,12,22]. EVI and NIR v share NDVI's sensitivity to LAI, while adding sensitivity to differences in reflectance in the NIR [6,12,20,23]. Prior GSA studies have produced ambiguous relationships between plant traits, EVI and NIR v , suggesting a strong sensitivity to either mean canopy leaf angle or to chlorophyll and leaf mass per area [20,24]. NIR reflectance is related to scattering at air-water interfaces within leaves [25], and studies Remote Sens. 2020, 12, 1405 3 of 18 have shown that NIR reflectance is correlated with LAI and a suite of leaf traits including leaf mass per area, nitrogen content, water content, and canopy traits such as mean leaf inclination angle and canopy shape [26][27][28][29][30]. Some of these leaf traits are correlated with each other, as well as leaf photosynthetic rate [27,29,31,32]. However, we are unaware of previous studies that have quantitatively assessed which of these plant traits are the main drivers of EVI and NIR v across a broader set of biomes, and how these traits mediate or affect the relationship between these VIs and GPP.
We explored the temporal and spatial relationships between NDVI, EVI, NIR v , and GPP across 10 homogenously-instrumented eddy covariance sites in California that spanned a range of climatic conditions. These sites varied markedly in vegetation density and type, and sampled the local analogues of six major terrestrial biomes (grassland, subtropical desert, shrubland, woodland/savanna, temperate forest, and boreal/subalpine forest). We then explored these relationships using radiative transfer and GPP models forced with in-situ plant trait and soil reflectance observations. We focused on two questions: Which Vegetation Index or Indices are most generally correlated with GPP across temporal and spatial scales? What biophysical property or properties mediate the relationships between VIs and GPP?

Study Sites
The study focused on existing and new observations at 10 eddy covariance sites across a climate gradient in California ( Figure 1). All sites use matched instrumentation and data processing and had existing 5-to-10 year-long eddy covariance records. The sites were located along a broad topographic and climatic gradient ranging from 275 m to 2700 m above sea level, a mean annual precipitation from 129 mm yr −1 to 1078 mm yr −1 and a mean annual temperature from 4.2 • C to 20.9 • C [33]. The sites sample 6 of the 10 major biomes of the world (according to [34]): grassland, subtropical desert, shrubland, woodland/savanna, temperate forest, and boreal/subalpine forest.

Eddy Covariance and Satellite Observations
Eddy covariance data from these sites has been used to explore the patterns and trends in carbon dioxide and evapotranspiration in California and the Southwestern US [33,35]. The ten sites share the same instrumentation, data processing, and maintenance [36]. The eddy covariance systems consist of a closed-path infrared gas analyzer (LI-7000, LiCor Biosciences) and a 3D sonic anemometer (CSAT-3, Campbell Scientific). Data processing includes corrections for sensor lag, mean wind rotation, and energy budget closure. The net CO 2 fluxes were partitioned into respiration and gross uptake (half-hour GPP) by extrapolating light-response curves to darkness. Our analysis focused on the gross uptake during sunny, midday (11:00 to 13:00 local time) periods and we removed observations during cloudy periods or with high vapor pressure deficits (>2 KPa). We focused on midday observations to avoid the confounding effects of seasonal changes in sunlight duration, and the effects of meteorological constraints on photosynthesis, such as low light or high evaporative demand.
We used surface reflectance data for the sites from the MODerate resolution Imaging Spectrometer (MODIS) MCD43A4.006 product (downloaded from AppEEARS; https://lpdaacsvc.cr.usgs.gov/). This product provides daily coverage of surface reflectance (atmospherically corrected) at 500 m spatial resolution adjusted to a nadir view angle and local noon. For most sites, we extracted the time series of the pixel containing the tower location and VIs were estimated from bands 1 (Red), 2 (NIR), and 3 (Blue) according to Equations (1)- (3). Some sites in Southern California have relatively small areas of homogeneous vegetation around the towers relative to the MODIS pixel, and for these sites we used observations for nearby larger patches of well-matched vegetation [33] (Supplementary Material 1). We selected the nearby pixels as the closest ones showing the most similar VIs values to the area around the towers using Landsat imagery, and the similarity in vegetation type and composition at the pixel area was confirmed by visual inspection. MODIS data were filtered for clouds, high aerosols, Remote Sens. 2020, 12, 1405 4 of 18 and conditions leading to poor corrections of sensor view and sun angles according to the quality bands in the MCD43A4.006 product, and for patchy snow based on the Normalized Difference Snow Index (NDSI).

Study Sites
The study focused on existing and new observations at 10 eddy covariance sites across a climate gradient in California ( Figure 1). All sites use matched instrumentation and data processing and had existing 5-to-10 year-long eddy covariance records. The sites were located along a broad topographic and climatic gradient ranging from 275 m to 2700 m above sea level, a mean annual precipitation from 129 mm yr -1 to 1078 mm yr -1 and a mean annual temperature from 4.2 °C to 20.9 °C [33]. The sites sample 6 of the 10 major biomes of the world (according to [34]): grassland, subtropical desert, shrubland, woodland/savanna, temperate forest, and boreal/subalpine forest. The time spans of eddy covariance and satellite observations used are shown on Figure 1, and ranges from 5 to 10 years of observations for each site, including all seasons of the year and all stages of canopy development.

Plant Traits and Soil Reflectance at the Study Sites
We sampled plant species cover, LAI, leaf traits, and soil reflectance at five of the study sites around their peak growing season between April and July 2019 ( Figure 1, see Supplementary Material 2 for specific sampling dates at each site). In-situ sampling focused on a 90 m by 90 m block in the immediate mean upwind direction from each tower. Plant species cover was measured by intercept along three 90m line-transects in the block. LAI was measured at 5 m intervals along each transect with a LI-COR LAI-2000 (LI-COR Biosciences, Lincoln, NE, USA). LAI measurements were made with a 90 o field of view restrictor and around sunrise or sunset to avoid direct sunlight.
Leaf samples were collected and leaf gas exchange measurements were made on five or six of the most abundant species at each site (according to species cover data). Five mature plants were sampled for each species, and five or six leaves (or approximately 0.2 g of leaves for small leaved species) were collected. Fresh leaf area and weight and dry weight were measured and used to calculate the leaf mass per area (LMA) and leaf water content (W c ) [37]. The light-saturated leaf photosynthesis rate per leaf area (A leaf ) was measured on one sun-exposed leaf of each plant using a LI-COR LI-6400 (LI-COR Biosciences) with a blue-red LED light. Gas exchange was measured at near ambient temperature and humidity, 2000 µmol m −2 s −2 photosynthetically active radiation and 400 µmol mol −1 CO 2 . The community weighted trait means were calculated for each site by weighting traits by species cover [38]. The bare soil spectral reflectance was measured at five locations using an ASD FieldSpec HandHeld 2 (ASD Inc, Boulder CO) at midday. Each measurement was the average of five readings taken with a 34 ms integration time, corrected by dark current and calibrated with a Spectralon panel. The median of the five measurements per site was used for further analysis.

Radiative Transfer and Gross Primary Production (GPP) Models
A simple canopy photosynthesis model [50] (Thornley's model) was used to explore the links between plant traits and GPP. This model integrates a non-rectangular hyperbola leaf photosynthetic light response curve through the canopy. Its main inputs are incident photosynthetically active radiation (PAR), leaf light response parameters, leaf transmittance, light saturated leaf photosynthetic rate (A leaf ), leaf area index (LAI), and canopy light extinction coefficient.
The PROSAIL radiative transfer model was used to explore the links between plant traits and spectral reflectance. PROSAIL couples leaf (PROSPECT 5) and canopy (4SAIL) radiative transfer models to simulate the effects of plant traits, soil reflectance, and view and sun geometry on land surface reflectance [51]. PROSAIL inputs include leaf chlorophyll per area (C ab ), Leaf Mass per Area (LMA), leaf structure (number of air-cell wall interfaces), carotenoid and brown pigment content, water content per leaf area (W c ), and Leaf Area Index (LAI).
We applied variance-based global sensitivity analyses [52,53] to the Thornley and PROSAIL models to identify the plant traits and other inputs that drive GPP and the VIs. The analysis of the Thornley model was run using the R sensitivity package v1.18.0. The analysis on PROSAIL was run using the Global Sensitivity Analysis module v.1.08 of the ARTMO toolbox v.3.25 [54]. Previous sensitivity analysis on VIs have emphasized crop species, and parameter ranges has been either tailored for such ecosystems or arbitrarily selected [20,24]. The results of global sensitivity analyses depend heavily on the range of inputs specified, and we selected ranges that represent the variation across our study sites or from literature datasets that covered a diverse array of species and biomes (see Appendix A for details). The global sensitivity analyses of both models used 2000 samples per input parameter following a Sobol sampling scheme.
Once the drivers and factors driving GPP and VIs were identified, we compared the results of Thornley and PROSAIL models forced with in-situ information from these main drivers to further explore how plant traits link VIs to GPP (see Appendix B for details).

Correlations between Vegetation Indices (Vis) and GPP
We analyzed the relationships between satellite VIs and midday GPP within individual sites to assess the effects of temporal variation within sites (separate regression lines for each site in Figure 2), and across all sites to assess the between-ecosystem type generality (dashed lines in Figure 2). The within site relationships between GPP and both EVI and NIR v were tighter than those for NDVI (Supplementary Material 3). The within site relationships were weakest in the three sites that were dominated by evergreens (oak-pine, sierra mixed conifer, and subalpine forests), especially for NDVI. Comparatively weak temporal relationships at evergreen sites could be expected given the low seasonal variation in canopy biophysical properties [11]. The temporal relationships between GPP and both EVI and NIR v had more consistent intercepts and slopes between sites than did those with NDVI ( Figure 2, Table S3). This resulted in stronger across-site relationships for EVI and NIR v than for NDVI. Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 18 NDVI ( Figure 2, Table S3). This resulted in stronger across-site relationships for EVI and NIRv than for NDVI.

Plant Traits Influencing GPP
Our global sensitivity analysis of Thornley's GPP model showed that LAI (m 2 leaves m -2 ground) and leaf CO2 uptake (Aleaf; the CO2 uptake m -2 leaf) are the main determinants of short-term GPP (CO2 uptake m -2 ground) during well-illuminated, unstressed conditions ( Figure 3a). Both aspects have a similar contribution to GPP (LAI accounts for 48%, and Aleaf accounts for 42%, of overall GPP variation). All of the other attributes, including mean leaf angle, have a relatively minor contribution. GPP modeled with these two traits closely matched the observed variation of GPP among sites (Supplementary Material 4).
We further explored how LAI and Aleaf influence GPP (Figure 3b). GPP increases with increasing LAI and the fraction of light intercepted by leaves before saturating and reaching 80% of its maximum

Plant Traits Influencing GPP
Our global sensitivity analysis of Thornley's GPP model showed that LAI (m 2 leaves m −2 ground) and leaf CO 2 uptake (A leaf ; the CO 2 uptake m −2 leaf) are the main determinants of short-term GPP (CO 2 uptake m −2 ground) during well-illuminated, unstressed conditions ( Figure 3a). Both aspects have a similar contribution to GPP (LAI accounts for 48%, and A leaf accounts for 42%, of overall GPP variation). All of the other attributes, including mean leaf angle, have a relatively minor contribution. GPP modeled with these two traits closely matched the observed variation of GPP among sites (Supplementary Material 4).
We further explored how LAI and A leaf influence GPP (Figure 3b). GPP increases with increasing LAI and the fraction of light intercepted by leaves before saturating and reaching 80% of its maximum response around 2.5-3 m 2 m −2 . Conversely, increasing A leaf caused a comparatively steady increase  Our analysis of leaf traits at the study sites showed that light-saturated leaf photosynthesis (Aleaf) was negatively correlated to LMA (log-log relationship, r 2 = 0.32, p = 0.002), with relatively thick, high LMA leaves showing a lower Aleaf. A similar pattern occurred across sites, where community weighted Aleaf decreased with increasing LMA (Figure 4). Previous analyses have also reported a correlation between LMA and Aleaf [56,57]. The relationship between LMA and leaf photosynthesis reflects a series of tradeoffs and convergent evolutionary strategies along the leaf economics spectrum: thicker leaves are longer lived, with lower nutrient content and slower rates of photosynthesis [31,32].  Our analysis of leaf traits at the study sites showed that light-saturated leaf photosynthesis (A leaf ) was negatively correlated to LMA (log-log relationship, r 2 = 0.32, p = 0.002), with relatively thick, high LMA leaves showing a lower A leaf . A similar pattern occurred across sites, where community weighted A leaf decreased with increasing LMA (Figure 4). Previous analyses have also reported a correlation between LMA and A leaf [56,57]. The relationship between LMA and leaf photosynthesis reflects a series of tradeoffs and convergent evolutionary strategies along the leaf economics spectrum: thicker leaves are longer lived, with lower nutrient content and slower rates of photosynthesis [31,32].
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 18 response around 2.5-3 m 2 m -2 . Conversely, increasing Aleaf caused a comparatively steady increase in GPP, with only minor saturation. LAI and Aleaf have a multiplicative effect over GPP, such that the sensitivity of GPP to Aleaf increases with LAI. Our analysis of leaf traits at the study sites showed that light-saturated leaf photosynthesis (Aleaf) was negatively correlated to LMA (log-log relationship, r 2 = 0.32, p = 0.002), with relatively thick, high LMA leaves showing a lower Aleaf. A similar pattern occurred across sites, where community weighted Aleaf decreased with increasing LMA (Figure 4). Previous analyses have also reported a correlation between LMA and Aleaf [56,57]. The relationship between LMA and leaf photosynthesis reflects a series of tradeoffs and convergent evolutionary strategies along the leaf economics spectrum: thicker leaves are longer lived, with lower nutrient content and slower rates of photosynthesis [31,32].

Plant Traits Influencing VIs and Their Linkage to GPP
Global sensitivity analyses with PROSAIL were used to investigate the importance of ecosystem properties on spectral reflectance (Figure 5a) and Vegetation Indices (Figure 5b). Red reflectance is strongly driven by LAI and chlorophyll concentration (C ab ), while NIR reflectance is influenced by LMA and to a lesser extent LAI. NDVI is particularly sensitive to LAI and C ab (driving 65 and 23% of overall NDVI variation, respectively), and mirroring the importance of these properties in driving red reflectance. NDVI is largely insensitive to LMA, apparently because the calculation of NDVI under typical conditions deemphasizes changes in NIR reflectance. EVI and NIR v were sensitive to LAI and LMA, with both properties showing a similar importance. These two VIs showed increased sensitivity to NIR reflectance relative to NDVI, a pattern that has been discussed previously [6,12]. Other parameters, including sun and sensor view angles and leaf angle, have a smaller effect on these VIs. The VIs predicted using PROSAIL driven with site-specific data captured the observed variation of MODIS VI at the study sites (Supplementary Material 4).

Plant Traits Influencing VIs and Their Linkage to GPP
Global sensitivity analyses with PROSAIL were used to investigate the importance of ecosystem properties on spectral reflectance (Figure 5a) and Vegetation Indices (Figure 5b). Red reflectance is strongly driven by LAI and chlorophyll concentration (Cab), while NIR reflectance is influenced by LMA and to a lesser extent LAI. NDVI is particularly sensitive to LAI and Cab (driving 65 and 23% of overall NDVI variation, respectively), and mirroring the importance of these properties in driving red reflectance. NDVI is largely insensitive to LMA, apparently because the calculation of NDVI under typical conditions deemphasizes changes in NIR reflectance. EVI and NIRv were sensitive to LAI and LMA, with both properties showing a similar importance. These two VIs showed increased sensitivity to NIR reflectance relative to NDVI, a pattern that has been discussed previously [6,12]. Other parameters, including sun and sensor view angles and leaf angle, have a smaller effect on these VIs. The VIs predicted using PROSAIL driven with site-specific data captured the observed variation of MODIS VI at the study sites (Supplementary Material 4). We further explored how LAI, Cab, and LMA and their interactions control the value of VIs ( Figure 6). NDVI showed the sharpest saturation to LAI, reaching 80% of maximum LAI response at ~1.7 m 2 m -2 . NDVI also showed a sharp saturation to Cab, reaching 80% of maximum response at a Cab of ~17 µg cm -2 , which is within the lower range of values observed at the sites (Supplementary material 5). NDVI's sharp saturation at low Cab and its low sensitivity to LMA prevents NDVI from fully capturing an optical signal of a leaf trait that is related to Aleaf (Figures 6a-b).
EVI and NIRv had a comparatively less sharp saturation to LAI than NDVI, with its saturation to LAI varying with LMA up to 3 m 2 m -2 . Furthermore, the value of EVI and NIRv saturated modestly with increasing LMA across the range of values observed at the sites. This implies that decreasing LMA from the lowest value (similar to that of the Pinyon-juniper woodland) to the highest value (similar to that of the grassland), causes a relatively large increase in EVI and NIRv. EVI and NIRv's balanced sensitivity to LAI and LMA, low saturation to the signal of a leaf trait related to Aleaf, and multiplicative sensitivity to LMA and LAI, allowed both VIs to largely mimic the response of GPP to LAI and Aleaf (Figure 3b and Figures 6c-d). We further explored how LAI, C ab , and LMA and their interactions control the value of VIs ( Figure 6). NDVI showed the sharpest saturation to LAI, reaching 80% of maximum LAI response at~1.7 m 2 m −2 . NDVI also showed a sharp saturation to C ab , reaching 80% of maximum response at a C ab of~17 µg cm −2 , which is within the lower range of values observed at the sites (Supplementary Material 5). NDVI's sharp saturation at low C ab and its low sensitivity to LMA prevents NDVI from fully capturing an optical signal of a leaf trait that is related to A leaf (Figure 6a,b). EVI and NIR v had a comparatively less sharp saturation to LAI than NDVI, with its saturation to LAI varying with LMA up to 3 m 2 m −2 . Furthermore, the value of EVI and NIR v saturated modestly with increasing LMA across the range of values observed at the sites. This implies that decreasing LMA from the lowest value (similar to that of the Pinyon-juniper woodland) to the highest value (similar to that of the grassland), causes a relatively large increase in EVI and NIR v . EVI and NIR v 's balanced sensitivity to LAI and LMA, low saturation to the signal of a leaf trait related to A leaf , and multiplicative sensitivity to LMA and LAI, allowed both VIs to largely mimic the response of GPP to LAI and A leaf (Figures 3b and 6c,d).

Relationships among VIs
EVI and NIRv showed similar patterns throughout the analyses (Figures 2b-c, Figure 5b, Figures  6c-d), leading us to wonder how tightly correlated these VIs are across typical conditions and whether they provide complementary or redundant information. A further exploration of the MODIS observations revealed a tight correlation between EVI and NIRv, and relatively weak relationships between both of these VIs and NDVI (Figure 7). We found similar relationships in the PROSAIL simulations, suggesting that a tight correlation is a basic characteristic of EVI and NIRv (Supplementary Material 6).

Relationships among VIs
EVI and NIR v showed similar patterns throughout the analyses (Figure 2b,c, Figure 5b, Figure 6c,d), leading us to wonder how tightly correlated these VIs are across typical conditions and whether they provide complementary or redundant information. A further exploration of the MODIS observations revealed a tight correlation between EVI and NIR v , and relatively weak relationships between both of these VIs and NDVI (Figure 7). We found similar relationships in the PROSAIL simulations, suggesting that a tight correlation is a basic characteristic of EVI and NIR v (Supplementary Material 6).

Vegetation Indices as General Predictors of GPP
Previous studies have found that no single VI provides a universal predictor of GPP and that the relationships between GPP and VIs are often site-or biome-specific [8,10,18]. We found that the within-site relationships between NDVI and GPP varied between ecosystem types and differed from that observed in the across-site comparison (Figure 2). EVI and NIR v , on the other hand, provided more general predictions of GPP, with temporal relationships that were comparatively consistent across sites. Even at evergreen forests, where all VIs were weaker predictors in our study and others [7,8,11,18], the GPP relationships with EVI and NIR v fall within the range of other sites/biomes. These more universal relationships underscore the utility of EVI and NIR v for extrapolating GPP in the absence of further site-or biome-specific information.
Our analyses revealed strong similarities between EVI and NIR v , and subsequent comparisons showed a tight correlation between EVI and NIR v (Figure 7 and Supplementary Material 6). Previous studies have also pointed to a similar performance between EVI and NIR v [10]. This correlation was unexpected given the mathematical differences between the two VIs. EVI was designed to reduce cross sensitivity to soil and atmospheric conditions, while reducing saturation in dense canopies [22]; previous studies have shown EVI has enhanced sensitivity to variation in NIR reflectance [23]. In a similar way, the design and derivation of NIR v arose from an attempt to isolate the NIR reflectance of the vegetation from that of the soil [6,58]. Our analyses confirm the strong similarity between EVI and NIR v is a consequence of the shared sensitivity to traits influencing NIR, and the reduced influence from Red reflectance and soil reflectance. We conclude that the information provided by EVI and NIR v is largely redundant, and that both are similarly useful for inferring GPP, phenology, and other ecosystem properties [5,6,58,59].

Explaining the Differential Performance of Vegetation Indices
The PROSAIL analyses showed that NDVI is particularly sensitive to LAI and C ab ( Figure 5; [11,19]), while EVI and NIR v are sensitive to LAI and LMA. EVI and NIR v 's increased sensitivity to LMA is the main characteristic that separates these VIs from NDVI, and that we believe leads to their more general relationships with GPP. The relationships between LMA and EVI and NIR v are apparently driven by a decrease in NIR reflectance with progressively thicker leaves. LMA is not a dominant controller of leaf-level NIR reflectance [60], but becomes important at the canopy level due to the absorption of light by lignin and cellulose [61]. Thicker, higher LMA leaves with more lignin and cellulose absorb more NIR for a given LAI, leading to increased light absorption, and reduced canopy-level NIR reflectance, and consequently reduced EVI and NIR v [51].
LMA is one of several leaf traits that co-vary along the leaf economic spectrum [32]. Thicker leaves with a higher LMA are globally correlated with longer leaf life spans, lower leaf nitrogen concentrations, and lower rates of leaf photosynthesis. The leaf economic relationships are strongest when compared on a mass vs mass basis, whereas our LMA vs A leaf relationship is presumably mediated by leaf biochemistry and should be thought of as a mass vs area comparison [31,32,56,57]. Nonetheless, previous studies have found that leaf nitrogen content per mass combined with LAI largely explains GPP variation across sites and biomes [15,28,62].
LMA may be thought of as a proxy for leaf nitrogen content and photosynthetic rate, and VIs that are influenced by both LAI and LMA will provide more general predictions of GPP than ones that are influenced mainly by LAI. Hence, seasonal or interannual GPP shifts associated with the initial flushing of young, lower LMA, rapidly photosynthesizing leaves, and the gradual aging to thicker leaves with lower photosynthetic rates [39,[63][64][65][66] may be effectively explored with VIs such as EVI and NIR v . Likewise, the GPP effects of seasonal, interannual, or successional shifts from herbaceous to evergreen leaf importance may be detected using VIs such as EVI and NIR v .
We have argued that EVI and NIR v are particularly well correlated with GPP because they are influenced by both LAI and LMA. However, we have also shown that NDVI is sensitive to C ab content (Figure 5b), which is also correlated with leaf level CO 2 uptake [21], raising the question of the apparent insensitivity of NDVI to the leaf-level component of GPP. Previous studies have assumed that a weaker relationship between NDVI and GPP is caused by saturation at high LAI or environmental influences that are not related to plant physiology [18]. Here we have shown that NDVI saturates at relatively low values of C ab (Figure 6a) due to strong chlorophyll absorption at red wavelengths ( Figure 5a). Hence, NDVI is relatively insensitive to C ab over the typical conditions observed for upper-canopy leaves, and NDVI lacks strong sensitivity to a leaf trait that is well correlated with A leaf . In contrast, the effect of LMA on EVI and NIR v shows minor saturation across a realistic range of leaf thicknesses (Figure 5c,d).

Implications and Opportunities
A site's GPP ultimately reflects the combination of LAI and A leaf (Figure 3). The ideal VI should show sensitivity to LAI and a similarly strong, multiplicative sensitivity to some aspect of the land surface that is closely related to A leaf (Figure 4, [50]). EVI and NIR v do a better job of mimicking this sensitivity than NDVI does (Figure 6c,d vs Figure 3b), and as a result provide tighter relationships with GPP ( Figure 2). Nonetheless, there is no reason to believe that EVI or NIR v represent the optimal solution to this problem and alternative VIs may do even better jobs of predicting GPP.
Recent advances in remote sensing and plant physiology may further improve predictions of GPP. Studies are exploring the relationships between plant traits and GPP across ecosystem types [14,15,28,62]. Plant trait complexes have relatively strong spectral signals over visible and near infrared broad-bands [67], and advances in simulating canopy radiation transfer should allow further explorations of the links between GPP and surface reflectance [68]. Alternative VIs that better mimic the relative importance of LAI and A leaf in controlling GPP, or that more correctly account for the rate that the LAI and A leaf controls on GPP saturate (Figure 3, Figure 6), may provide even more general predictions of GPP. Likewise, alternative VIs or retrieval approaches that combine sensitivity to LAI with a proxy that is better correlated with A leaf than is LMA (Figure 4; [21,31,32]) may provide improved predictions of GPP.

Conclusions
We explored the relationships and underlying linkages between Gross Primary Production (GPP) and three commonly used Vegetation Indices: NDVI, EVI, and NIR v . We analyzed data from 10 eddy covariance sites located in the California analogues of six major terrestrial biomes, and found that GPP, EVI and NIR v have strong and convergent within-and between-site relationships, whereas NDVI and GPP showed weaker and more site-specific relationships.
Our analyses showed that GPP is driven by the multiplicative interaction between leaf area index (LAI) and leaf-level CO 2 uptake (A leaf ), with the importance of A leaf increasing with LAI. This analysis led us to conclude that a VI that is sensitive to the interacting effect of LAI and A leaf should provide better predictions of GPP than a VI that is sensitive to just one of these properties. Moreover, we found that leaf mass per area (LMA) provides a proxy for A leaf , with thicker leaves showing lower rates of A leaf. Finally, we found that EVI and NIR v are strongly sensitive to LMA and LAI, whereas NDVI is most strongly controlled by LAI over a range of realistic conditions. Hence, we conclude that the comparatively tight relationships between GPP, EVI, and NIR v are driven by the multiplicative sensitivity of these two VIs to LAI and LMA, which largely mimics the response of GPP to LAI and A leaf . NDVI is comparatively less sensitive to a leaf trait that is tied to leaf photosynthesis, which explains its weaker correlation with GPP.
Our analysis provides a roadmap for more mechanistically linking VIs, ecosystem biophysical properties, and GPP, and implies that VIs that better mimic GPP's sensitivity to LAI and A leaf may provide more general predictions of canopy photosynthesis.  where I 0 is the beam irradiance and θ sz is the sun zenith angle. Following Ross [70] and Fuchs et al. [71], the sun zenith angle and the canopy mean leaf angle θ l influence the canopy extinction coefficient K according to Our sensitivity analysis of GPP did not explicitly incorporate the effect of other potentially important variables such as air or leaves temperature. However, previous studies have shown that the important controlling effect of plant traits over GPP holds regardless of variation in these meteorological conditions [15,62].
We ran the global sensitivity analysis of the PROSAIL model over the following ranges: 9-58 µg cm −2 for leaf chlorophyll content (see Supplementary Material 5, [72]); 1-4 for the leaf structure parameter [73]; 0.0064-0.0526 g cm −2 for LMA; 0.0146-0.0771 g cm −2 for W c , and 0.1-5 m 2 m −2 LAI. The ranges of LMA and W c were determined by the variation in community weighted mean traits observed at the sites. The range of LAI exceeded the actual ranges observed at the sites (observed range 0.2-3.0 m 2 m −2 at 5 sites). Ranges for other parameters were 25-63 • for mean leaf angle (based on a review of published values [74], see Supplementary Material 7), 0.0001-0.01 for hotspot (approximate range for the sites according to leaf width divided by canopy height), 0-1 for the soil parameter, 0-1 for proportion of diffuse to direct radiation, 14-65 o for solar zenith angle, 0-60 • for view zenith angle of the sensor, and 0-180 o for the relative azimuth angle between the sun and sensor. The ranges for solar and view angles are based on exploratory analyses of MODIS imagery and MODIS specifications. We ignored variation in carotenoid and brown pigments, given the paucity of field observations and lack of a major effect on reflectance [75].
Previous global sensitivity analyses have used the maximum possible range for mean canopy leaf angles (i.e., 0-90 • ), and have shown that leaf inclination exerts a strong control over NIR and vegetation indices such as EVI [20,54]. Our review of published measurements of mean leaf angles shows a comparatively narrower range of leaf inclination angles (Supplementary Material 7), and suggests that some previous analyses may have overestimated the importance of leaf inclination in NIR reflectance.

Appendix B. Analysis with Thornley (2002) and PROSAIL Models
We used the model of Thornley [50] to predict gross primary production (GPP) in Figure 3b. This model integrates a non-rectangular hyperbola leaf light response curve through the canopy, accounting for leaf acclimation to sun and shade conditions. It models GPP from incident photosynthetically active radiation (PAR), the leaf light response parameters of light saturated photosynthetic rate sunlit leaves (A leaf ), maximum quantum yield, and a light response curvature parameter. Canopy inputs are leaf area index (LAI), light extinction coefficient, and leaf transmission coefficient. For Figure 3b, we set the quantum yield to 9.118121 × 10 −9 Kg CO 2 J −1 (the average value for C 3 species from [69]), a curvature parameter of 0.9 [50], a proportion of direct radiation of 0.9, an extinction coefficient of 0.5 (spherical leaf angle distributions, [76]), leaf transmission coefficient of 0.1 (see Appendix A, [50]), and a PAR of 435 W PAR m −2 (equivalent to about 2000 µmol m −2 s −1 ). We varied LAI from 0-5 m 2 m −2 and A leaf 5-20 µmol m −2 s −1 , approximating the range of community weighted means found at our sites (5.8-19.4 µmol m −2 s −1 ).
The PROSAIL model predictions in Figure 6 were calculated with the following inputs: 1.5 for leaf structure parameter (typical value for healthy leaves, [73]), 0 for leaf carotenoid and brown pigment content per area (see Appendix A), the average leaf water observed at the field sites (0.0384 g cm −2 ), a spherical leaf angle distribution, hot spot value of 0.01, 10% fraction of diffuse to direct radiation, and a fixed soil reflectance; 0 • for view angle (from nadir), 0 • view-sun angle azimuth difference, and sun zenith angle of 22 • (average value at midday during spring at our sites). LAI was varied from 0 to 5 m 2 m −2 , and leaf chlorophyll content (C ab ) or leaf mass per area (LMA) were varied or fixed as specified in the figure panel ( Figure 6); fixed values were 30 µg cm −2 for C ab (see Supplementary Material 5) and 0.024 g cm −2 for LMA (the average at our sites). The levels of C ab or LMA used in Figure 6 correspond to the leaf photosynthetic rates shown in Figure 3b (see Figure 6 legend), and approximate the ranges at the field sites ( Figure 4, Supplementary Material 5). The quantitative crosswalk between LMA in Figure 6 and leaf photosynthetic rate ( Figure 5 legend) was based on Figure 4. The crosswalk between C ab and leaf photosynthetic rate (Figure 6 legend) was based on [21], which related C ab to the leaf maximum carboxylation and electron transport capacities, and then extrapolated to A leaf following [77] implemented in plantecophys R package [78]. The leaf photosynthesis model assumed a leaf temperature of 25 • C and an intercellular CO 2 concentration of 256 µmol mol −1 (the average values in our field measurements).