Remote Sensing Based Simple Models of Gpp in Both Disturbed and Undisturbed Piñon-juniper Woodlands in the Southwestern U.s

Remote sensing is a key technology that enables us to scale up our empirical, in situ measurements of carbon uptake made at the site level. In low leaf area index ecosystems typical of semi-arid regions however, many assumptions of these remote sensing approaches fall short, given the complexities of the heterogeneous landscape and frequent disturbance. Here, we investigated the utility of remote sensing data for predicting gross primary production (GPP) in piñon-juniper woodlands in New Mexico (USA). We developed a simple model hierarchy using climate drivers and satellite vegetation indices (VIs) to predict GPP, which we validated against in situ estimates of GPP from eddy-covariance. We tested the influence of pixel size on model fit by comparing model performance when using VIs from RapidEye (5 m) and the VIs from Landsat ETM+ (30 m). We also tested the ability of the normalized difference wetness index (NDWI) and normalized difference red edge (NDRE) to improve model fits. The best predictor of GPP at the undisturbed PJ woodland was Landsat ETM+ derived NDVI (normalized difference vegetation index), whereas at the disturbed site, the red-edge VI performed best (R 2 adj of 0.92 and 0.90 respectively). The RapidEye data did improve model performance, but only after we controlled for the variability in sensor view angle, which had a significant impact on the apparent cover of vegetation in our low fractional cover experimental woodland. At both sites, model performance was best either during non-stressful growth conditions, where NDVI performed best, or during severe ecosystem stress conditions (e.g., during the girdling process), where NDRE and NDWI improved model fit, suggesting the inclusion of red-edge leveraging and moisture sensitive VI in simple, data driven models can constrain GPP estimate uncertainty during periods of high ecosystem stress or disturbance.


Introduction
Semi-arid regions and drylands together cover more than 40% of the globe.In spite of the low fractional cover of vegetation and minimal annual precipitation, the contribution of semi-arid biomes to the global uptake of carbon dioxide from the atmosphere can be significant on annual time scales [1,2].The inter-annual variability of precipitation in these ecosystems can result in both rapid carbon stock accumulation, and subsequent turnover due to frequent disturbance (e.g., drought, insect outbreak, fire), but how these relationships alter the net storage of carbon is poorly understood (e.g., [1]).The climate in the southwestern USA has recently experienced both increased regional air temperature, and decreased precipitation, both trends which are projected to continue in the coming decades [3][4][5][6][7].The combined effects of these changes in climate have increased drought severity, exacerbated ecosystem stress, and ultimately triggered widespread forest mortality across the region [8][9][10].Given the rapid state transitions of vegetation in these climatically sensitive biomes, and the projected changes in climate and associated mortality for this region [11], accurate monitoring of carbon uptake dynamics in these semi-arid ecosystems is an essential part of constraining uncertainties associated with regional carbon balance.
Monitoring ecosystem carbon uptake over large geographic extents requires the use of remote sensing.Typically, remote sensing driven models of ecosystem function predict Gross Primary Productivity (GPP), the total atmospheric carbon taken up via photosynthesis (e.g., [12,13]).Satellite remote sensing provides the means to predict GPP by employing light use efficiency-based algorithms [14,15] that rely on vegetation indices (VIs), such as the normalized difference vegetation index (NDVI, [16]).This approach has been used at a variety of sites and scales, both by employing empirical (e.g., [17,18]) and light use efficiency/process driven techniques (e.g., [19]), generally using high temporal resolution, coarse spatial resolution (ě250 m) remote sensing data from sensors, such as the Moderate Resolution Imaging Spectrometer (MODIS) (e.g., [12]).To this end, the success of most space borne estimates of GPP hinge on a robust relationship between satellite estimated light interception and photosynthesis, via light use efficiency models, which generally perform well relative to in situ measurements from flux tower networks [20].
However, in semi-arid regions, vegetation function is constrained by water availability for the majority of the year, with periods of drought and heat interrupted briefly by spring snow melt or episodic pulses of precipitation.Consequently, in these biomes, changes in NDVI are constrained primarily by the availability of water, rather than light or temperature [21].This results in NDVI and GPP often being decoupled in semi-arid regions due to low soil moisture, particularly where evergreen plants are present [22].Other characteristics unique to semi-arid ecosystems that challenge the use of NDVI for characterizing changes in GPP in these biomes are highly variable precipitation patterns which trigger high inter-annual variability in GPP due to seasonal water limitations, but low variability in LAI and/or chlorophyll concentration (subsequently low variability in NDVI).Further, the low LAI (<1.5 mean LAI across landscapes, e.g., [23]) and spatially heterogeneous plant canopies typical of these systems can result in further uncertainties due to high reflectance by the soil background, thus confounding spectral signals relating to plant function (e.g., [24][25][26]).
One of the most universal responses to leaf stress is increasing visible reflectance [27], due to a combination of stressors which ultimately reduce chloropyll a + b, and consequently reduce the absorption of incident light [27,28].Chlorophyll a + b strongly absorb in the red portion of the visible spectrum, resulting in saturation of the red band at low Chlorophyll a + b, and reducing its potential to track initial chlorophyll loss at the onset of stress.However, the red-edge has been shown to have a more linear response to a wide range of chlorophyll concentrations, increasing its potential as a stress indicator in vegetated systems over conventional red-NIR combinations [26,[29][30][31][32].Given the increasing availability of the red-edge waveband in commercial (e.g., RapidEye, WorldView-2, WorldView-3) sensors, and freely available data (Sentinel-2, launched June 2015), testing the potential of the red-edge waveband to improve modeled estimates of GPP in semi-arid ecosystems is becoming a feasible task.
Here, we test the ability of spectral VI's other than NDVI to model GPP in semi-arid ecosystems.We focus particularly on piñon-juniper woodlands for several reasons.First, because it is the largest biome in the Southwestern US, covering 18 million ha in New Mexico, Arizona, Colorado and Utah.Second, the changes in climate in this region have triggered a significant amount of mortality in this biome [10,33], and thus, quantifying the extent of this disturbance on productivity throughout the region is crucial to understanding how this extensive mortality has impacted both current and future carbon dynamics.Finally, we take advantage of an existing experimental manipulation in a Piñon-Juniper woodland in central New Mexico that was girdled in 2009 to simulate the widespread piñon mortality observed throughout the region.The advantage of using this experimental manipulation is that since the girdling in 2009, changes in ecosystem productivity triggered by this mortality have been continuously monitored using eddy covariance, and compared to similar measurements in a nearby intact PJ woodland that serves as a control.
Recent research [34] conducted in our experimental manipulation suggested that the red-edge employing normalized difference red-edge index (NDRE) [35] was more sensitive than NDVI to the observed initial decrease in leaf chlorophyll concentration triggered by piñon girdling in this PJ woodland, adding evidence to its potential as a stress monitoring component in GPP modeling efforts.The NDRE was determined to be more sensitive to changes in greenup of the low LAI herbaceous vegetation following mortality in this same system [22].Based on these findings, the goal of this study was to test three specific hypotheses.The first hypothesis being that the observed variability in NDRE will allow more accurate estimation of GPP in both the disturbed and undisturbed PJ woodland in this experimental manipulation.Secondly, due to the inherent dependence of productivity on water availability in this biome, adding the normalized difference wetness index (NDWI), a VI that is sensitive to changes in foliar water content, will significantly improve the model fit by constraining the model error during the dry season.Finally, given the low vegetation cover and highly heterogeneous nature of PJ woodlands, that VIs generated from higher spatial resolution (5 m) data will provide more accurate estimates of GPP relative to traditional, moderate resolution (30 m) remote sensing data.

Site Description
The study site includes two piñon-juniper woodlands separated by 3 km and located south of Mountainair, NM.In September of 2009, 1632 adult piñon (>7 cm diameter at breast height) in a 4 ha plot in one of the sites (hereafter referred as the girdled site; 34 ˝26 1 48.54"N, 106 ˝12 1 48.63"W) were mechanically girdled by severing the sapwood with a chainsaw at breast height (1.4 m) followed by the application of a 5% glyphosate solution directly applied to the cut to ensure rapid loss of conductivity and subsequent mortality of the trees.The total decoupling of the leaves with the roots following the manipulation was designed to replicate the landscape following the mortality that occurred in the response to the 1999-2002 drought [9,27].The second PJ woodland (hereafter referred to as the control site; 34 ˝26 1 18.420"N, 106 ˝14 1 15.698"W) remained intact.We have continually monitored surface fluxes of carbon, water and energy at each site using tower mounted eddy covariance and associated micrometeorological sensors since 8 months prior to the girdling manipulation in February 2009.
The soil at both sites is Turkey Springs stony loam, characterized by an abundance of alluvially deposited limestone (Soil Survey Staff).The climate of the region is semi-arid, with a mean annual precipitation of 372 mm (˘86.8 mm standard deviation, sd) and max, and min monthly temperatures of 19.8 ˝C (˘0.77 ˝C sd) and 2.32 ˝C (˘0.64 ˝C sd) respectively, over the past 20 years (PRISM Climate Group, Oregon State University).Incoming moisture to the site is largely bimodal, broken into winter snow melt from January to March and seasonal monsoon precipitation between August and October with a pronounced dry season occurring from April through July.

Gross Primary Productivity
Data processing and instrumentation is identical at both sites.Eddy covariance (EC) derived surface fluxes of carbon, water and energy were measured at 10 Hz using a 3-axis sonic anemometer (CSAT-3, Campbell Scientific, Logan, UT, USA) and an open-path infrared gas analyzer (Li-7500, LiCor Biosciences; Lincoln, NE, USA).Continuous measures of net radiation, air temperature and relative humidity, soil temperature, photosynthetically active radiation, and soil moisture (volumetric water content) were made using a CNR1 (Kipp and Zonen, Delft, The Netherlands), HMP45C (Vaisala RH probe with aspirated radiation shield, Vantaa, Finland), TCAV (averaging thermocouple probes, 27 per site), up facing quantum sensors (Li-190SB, Licor Biosciences, Lincoln, NE, USA) and ECHO probes (TE 5 cm, Decagon Devices, Pullman, WA, USA, 27 per site) respectively.The fluxes were aggregated to 30 min intervals and were corrected for temperature and moisture variations (WPL, [36]) as well as frequency responses according to Massman [37].Anemometer tilt due to terrain variability was corrected using a planar fit method.We used a friction velocity (u*) filter to reject data obtained when turbulence is low (u* less than a threshold value).Data gaps created by the u* filter, malfunctioning instruments, and rain were filled following Lasslop et al. [38].We partitioned the gapfilled net ecosystem exchange (NEE) into the components of total C uptake through photosynthesis (i.e., GPP) and total carbon leaving the ecosystem through both autotrophic and heterotrophic respiration (Re).We used exponential relationships between nighttime NEE and temperature following the methods of Lasslop et al. [38] to calculate continuous ecosystem respiration during the day.GPP was then calculated as NEE-Re [39].
We used a flux source area model [40] modified for 2 dimensions by Detto et al. [41], to characterize the experimental region at both sites that are measured by the flux towers.The source area model suggested that the four hectare analysis region at each site accounted for approximately 80% recovery at the tower level.In this way we ensured that the remote sensing data and the EC measurements were representing the same patch of vegetation.

Landsat ETM+
We used the web enabled Landsat database (WELD) web service [42] to download time series of Landsat ETM+ data from 2009 to 2011 for 36 pixels within the four hectare area measured by each flux tower.Each pixel time series was cleaned (i.e., filtered to meet the following conditions) using the QA/QC data associated with WELD products [42]; only the pixels classified as containing no saturated bands and not cloudy (determined by the automatic cloud cover assessment algorithm, ACCA) [43] were retained.To remove effects of snow cover on the VI, the normalized difference snow index (NDSI) was used to remove snow covered pixels from the analysis.If more than 25% of the total analysis region for either the control or girdled site was covered by snow, the entire acquisition date was removed from the analysis.Hall et al. [44] suggest mapping snowy pixels as NDSI >0.4.However, we used a lower threshold of 0.2 in this study, which we validated with ground based imagery and field technician reports.The majority (82%) of the Landsat ETM+ data were acquired within ď3 days of a RapidEye acquisition (see Section 2.2.3), with the greatest time separation being 10 days.The Landsat ETM+ pixel time series was then linearly interpolated to daily intervals for comparison with the RapidEye data.In this manner, the same number of Landsat ETM+ VI and RapidEye VI were used in each of the models tested.Geolocation error for the Landsat ETM+ imagery was generally <0.5 pixel according to the WELD documentation.

RapidEye
A time series of 46 images (September 2009 to October 2011) were used with a top of atmosphere, dark object subtraction (TOA-DOS) correction applied by RapidEye, Inc. (now BlackBridge Ltd. of Berlin, Germany).The RapidEye data were reasonably registered on delivery, but to minimize variability in the spectral signature due to pixel shifts between scenes, we manually co-registered each image to a master image (September 2009) using 50 ground control points and a second order polynomial transformation.A total of 15-20 check points per image resulted in each scene being co-registered to a root mean square error (RMSE) <2.5 m (0.5 pixels).Eighteen images were excluded from the time series due to snow or cloud cover.Existing site boundary shapefiles were used to clip the images to the extent of both the control and girdled sites, coincident with the pixels extracted for the Landsat ETM+ analysis (Python Software Foundation.Python Language Reference, version 2.7). Figure 1 illustrates the distribution of RapidEye images throughout the duration of this experiment for the girdle and control site.

Spectral Vegetation Indices
NDVI (Normalized Difference Vegetation Index), NDRE (Normalized Difference Red Edge), and NDWI (Normalized Difference Wetness Index), were calculated to inform our simple models of GPP and to test our hypotheses.NDRE and NDWI are only available using RapidEye and Landsat ETM+, respectively, and indicated by their subscripts in Table 1.The red edge here refers to a waveband on the RapidEye constellation of sensors that spans 690 to 730 nm.

Spectral Vegetation Indices
NDVI (Normalized Difference Vegetation Index), NDRE (Normalized Difference Red Edge), and NDWI (Normalized Difference Wetness Index), were calculated to inform our simple models of GPP and to test our hypotheses.NDRE and NDWI are only available using RapidEye and Landsat ETM+, respectively, and indicated by their subscripts in Table 1.The red edge here refers to a waveband on the RapidEye constellation of sensors that spans 690 to 730 nm.

Statistical Analysis
We fitted six multiple linear regression models in open-source software package R 3.0.1 (R Development Core Team, 2013).Each regression model is structured as follows: where PAR is photosynthetic active radiation, TA is air temperature, and VI depicts one or more vegetation index that are shown in Table 1.GPP, PAR, and TA were derived from on-site instrumentation (for GPP, see Section 2.2.1).The VI were calculated as the mean of all pixels that passed the qa/qc steps mentioned in Sections 2.2.2 and 2.2.3.Thus the spatial resolution of the driving VI was either 5 m (from RapidEye), 30 m (from Landsat ETM+), or a combination in the case of the mixed models.Within each model group, we created a hierarchical structure of models containing all possible combinations of parameters and one first order interaction term.In this way we aimed to compare model groups by comparing the models with the best set of parameters within each group.Model performance was characterized using the root means squared error (RMSE) and the adjusted R 2 (R 2 adj ) that adds an additional penalty for increasing model complexity: where n = the number of observations.To identify the best model, we further calculated the Akaike information criterion (AIC adj ) : where k = the number of parameters .We chose the AIC adj to account for the relatively small number of observations and to avoid model over fitting [45].The best model for both the control and girdle site was chosen based on the largest AIC adj weight, and largest relative log likelihood, a measure of how likely a given model is actually the best.

Models of GPP in PJ Woodlands Had the Lowest Error at the Disturbed Site
Simple linear models driven by photosynthetically active radiation (PAR), air temperature (TA), and remotely sensed VI, summarized in Table 1, were able to explain between 80% and 90% of the variability of gross GPP in both intact and disturbed PJ woodlands (e.g., Table 2).However, the strength of this relationship was contingent on a tight coupling between the variability in VI and canopy physiological processes that affect a given index.For this reason, the model was run on individual years, as well as on all years together.Overall, the best performing model VI at the undisturbed PJ woodland was Landsat ETM+ derived NDVI, whereas at the disturbed site the red-edge VI performed best (R 2 adj of 0.90 and 0.92, respectively).When we drove the model hierarchy with data from all three years in the study, the collection of models performed poorly across the board (Figure 2, Table 3), with the Landsat ETM+ based models performed best (R 2 adj = 0.40 at control, R 2 adj = 0.50 at girdle).We believe this to be due to the coupling strength between changes in canopy VI and canopy function during as a function of interannual variability in climate.As an example, the fit statistics from both the Landsat ETM+ and RapidEye driven models improved significantly when a single year (e.g., 2009) was used to drive the regression (Figure 3), dramatically increasing the fit of the simple Landsat ETM+ NDVI model at the control site (R 2 adj from 0.40 to 0.90) and increased the fit of the RapidEye NDRE model at the girdled site (R 2 adj from 0.14 to 0.92) (Table 2, discussed in following sections).At the control site, Fall 2009 was very productive, with a strong monsoon bringing sufficient moisture to the ecosystem, and consequently the coupling of Landsat ETM+ NDVI, PAR, and air temperature explained almost 90% of the variability seen in GPP.In the case of girdled site during the same period, the manipulation in piñon mortality drove both the rapid decrease in VI, and the sudden decrease in canopy uptake of carbon.Both during this period of significant canopy stress and throughout the entire analysis period, the best model fits were seen at the girdled site, possibly due to the drastic decreases in LAI and subsequent increased coupling between the VI and GPP, previously documented as a potential effect of the manipulation itself [22].

7
of the variability seen in GPP.In the case of girdled site during the same period, the manipulation in piñon mortality drove both the rapid decrease in VI, and the sudden decrease in canopy uptake of carbon.Both during this period of significant canopy stress and throughout the entire analysis period, the best model fits were seen at the girdled site, possibly due to the drastic decreases in LAI and subsequent increased coupling between the VI and GPP, previously documented as a potential effect of the manipulation itself [22].One to one plots of tower derived (measured) and predicted (modeled) gross primary productivity for the control (C, black) and girdled (G, grey) sites for all years of the study.See Table 1 for a description of model names.
Table 3. Model performance statistics for all years of data.See Table 1 for a description of model names.R 2 adj is the adjusted R 2 , ∆AICc is the difference in model AICc relative to the lowest scoring model, RelLL is the relative log likelihood, or how likely it is that each model is actually the best, and Weights is a list of the AICc weights for each model.   1 for a description of model names.  1 for a description of model names.

NDRE and NDWI Reduced Model Error only during Periods of Significant Stress
In our linear model hierarchy, including NDRE did not significantly improve model performance relative to models that used only NDVI (Table 2).Previous work in this region suggests NDRE is more responsive to decreases in canopy chlorophyll concentration relative to NDVI, yet following a lag period (roughly 2 weeks), the NDVI performed equally well in quantifying foliar concentration of chlorophyll [34].Given that our time series spanned multiple years, with rapid decreases in canopy chlorophyll occurring only at the girdled site in 2009, we ran the same model structures on only data from 2009, which captured the girdling process (Figure 3).Over this smaller time period of significant canopy stress, the NDRE dramatically improved model performance at the girdled site (R 2 adj from 0.136 for NDVI to 0.92 for NDRE, Table 2).However, we did not observe significant improvement in NDRE driven model fit at the control site, which is consistent with the idea that NDRE explains variability in canopy function best during early periods of rapid canopy change.
The incorporation of Landsat ETM+-derived NDWI into the model structure did not improve model fit in either Landsat ETM+ or RapidEye based models when the entire time series was used in the regression (Figure 2, Table 3).NDWI did however slightly reduce prediction error in both Landsat ETM+ and RapidEye NDVI based models at the girdled site, in 2009 alone (Figure 3, Table 2).While the inclusion of the moisture index slightly improved model performance, the added complexity was penalized fairly heavily given the small number of observations used to parameterize the model.
Previous work has shown that including Landsat ETM+ moisture indices can improve estimates of productivity in some dryland ecosystems [21,46], although it appears to be less effective in sparse canopy systems [47].Our data suggest NDWI may only improve estimates of productivity in semi-arid coniferous woodlands during rapid changes in canopy physiological status, such as during canopy mortality, yet do not provide sufficient evidence to justify increasing the complexity of the model structure in this case.Moisture sensitive VI are less effective at representing water content during periods of low to moderate drought stress [48,49], suggesting that GPP predictions during periods of the year when water is less limiting may not benefit from the inclusion of a canopy moisture VI as much as during periods of drought stress.This is evident at the girdle site during 2011, a year which was characterized by significant drought across the southwestern US.Model performance improved at the girdle site with the inclusion of Landsat ETM+ NDWI.However, the increased complexity of the model and our small sample size limit the generalizability of the finding that NDWI may constrain model uncertainty in these systems during drought periods (Table 4).Our results from the entire 3 year period do not support the hypothesis that VI generated from higher spatial resolution data should be a better predictor of ecosystem carbon uptake (Table 2).In fact, using higher spatial resolution data from RapidEye surprisingly decreased model performance at the control site in some cases, and did not significantly improve model performance at the girdled site.Some of this can be explained by the inconsistent scene to scene sensor view angle in our RapidEye time series, which influenced the apparent fractional cover of vegetation (Figure 4A).In an effort to test whether or not the variability in view angle was imposing noise into the RapidEye data, we binned the analysis by sensor view angles <7 ˝C off nadir, given the majority of our data were collected at more than 7 degrees off nadir (60%, Figure 4B).Limiting the analysis to scenes with absolute sensor view angle <7 degrees significantly improved model fits in NDVI RapidEye mods (at the girdled site (R 2 adj increased from 0.14 to 0.77), yet decreased model performance at the control site; Table 5, Figure 5), potentially due to the changing apparent fraction of under-story and inter-canopy vegetation as a function of sensor view angle, which has a more pronounced impact at the girdle site relative to the control due to the lower fractional cover of overstory vegetation post manipulation.The inclusion of view angle into the model hierarchy as a parameter resulted in model over-fitting given the small sample size.The high inter-annual variability in system state (e.g., monsoonal precipitation or severe drought) at these sites further confounded our ability to account for view angle effects, due to the coupling between canopy function and the remotely derived VI for some dates being much tighter than others (see Section 2.1).

10
(at the girdled site (R 2 adj increased from 0.14 to 0.77), yet decreased model performance at the control site; Table 5, Figure 5), potentially due to the changing apparent fraction of under-story and intercanopy vegetation as a function of sensor view angle, which has a more pronounced impact at the girdle site relative to the control due to the lower fractional cover of overstory vegetation post manipulation.The inclusion of view angle into the model hierarchy as a parameter resulted in model over-fitting given the small sample size.The high inter-annual variability in system state (e.g., monsoonal precipitation or severe drought) at these sites further confounded our ability to account for view angle effects, due to the coupling between canopy function and the remotely derived VI for some dates being much tighter than others (see Section 2.1).One to one plots of tower derived (measured) and predicted (modeled) gross primary productivity for the control (C, black) and girdled (G, grey) sites for RapidEye acquisition angles < 7 degrees off nadir.See Table 1 for a description of model names.One to one plots of tower derived (measured) and predicted (modeled) gross primary productivity for the control (C, black) and girdled (G, grey) sites for RapidEye acquisition angles <7 degrees off nadir.See Table 1 for a description of model names.
Previous work indicates the effects of sensor view angle on the fractional cover of various feature types (e.g., [50], pertaining to composites derived from the advanced very high radiometric resolution spectrometer (AVHRR)), and generally near-nadir view angles are a selection criteria for multi-image or multi-sensor analyses (e.g., [51]).Given the increasing amount and complexity of remote sensing data at the disposal of terrestrial scientists, multi-resolution and time series data assimilation approaches to developing temporally and spatially resolved remote sensing products will need to address the confounding effects of view angle on apparent vegetation cover, especially in sparse canopy systems.

Implications for Regional Remote Sensing Based Estimations of GPP in PJ Woodlands
The results of our case study suggest that remote sensing driven, simple linear models of GPP have the potential to accurately describe patterns in regional carbon uptake, using locally measured PAR and air temperature as covariates.In this study, we chose to use Landsat ETM+ data, given its ease of access via the WELD portal.The highly heterogeneous composition of PJ woodlands, and the small scale of our manipulation experiment (200 m ˆ200 m), precluded the use of MODIS sized pixels (ě250 m ˆ250 m), and given the heterogeneous patterns of mortality and subsequent recovery typical of disturbed PJ woodlands (e.g., [9,22,27]), Landsat scale or finer resolution (ď30 m) products are of the appropriate spatial scale to resolve the patchy, heterogeneous patterns of mortality typical in disturbed PJ woodlands.In spite of the sensor view angle complications described in Section 3.3, limiting image acquisition angle variability impacted model performance in some cases, suggesting that increased spatial resolution data may have a role in constraining GPP predictions in these heterogeneous woodlands in the future, yet we were unable to fully test the hypothesis due to inadequate sample sizes when considering the variation in sensor view angle.
The inclusion of the red-edge leveraging NDRE from RapidEye only improved model performance during periods of significant stress, such as during the selective mortality of piñon pine that we imposed on the girdled site in Fall 2009, or during periods of severe drought such as in 2011 (see Section 3.2).Given the imminent transition in climate projected for the southwestern US, and the already evident impacts on piñon mortality regionally, GPP prediction efforts in PJ woodlands should benefit from the incorporation of NDRE during mortality.Given recent launch of Sentinel-2, which can leverage the red-edge and infrared wavebands required to compute NDRE and NDWI respectively, the potential to incorporate NDRE into near-future regional models of GPP as a consequence of piñon mortality may significantly constrain the current uncertainty in estimating PJ woodland ecosystem function.Further, while NDWI did not significantly contribute to model performance in this study, the measured productivity of semi-arid plant canopies is strongly a function of available water, which ultimately affects the coupling between canopy VI and GPP (e.g., [22]), and consequently the predictive power of VI driven models of GPP in these regions.Previous work has shown NDWI to correlate reasonably well with foliar water content, and water potential in piñon pine, but not juniper [52], suggesting further that the inclusion of these VI to empirical modeling attempts should be perhaps constrained to piñon dominated pixels, requiring higher resolution imagery.Remotely sensed regional estimates of surface soil water content from the soon to launch Soil Moisture Active Passive sensor (SMAP), may provide an alternative to constrain the inter-annual variability in the VI ~GPP relationship in semi-arid PJ woodlands.

Conclusions
Our results are promising in that we can use simple linear models to estimate GPP in both disturbed and undisturbed PJ woodlands driven by remotely sensed datasets.While structurally sensitive, NDVI is more informative to models of GPP than NDRE except during periods of extreme stress or disturbance.Similarly, we only saw a significant improvement in model performance using NDWI at the girdled site, during the manipulation event that took place in the Fall of 2009.Finally, the use of the RapidEye data did slightly improve estimates of GPP in both the control and girdled sites relative to Landsat ETM+, however this was only true when we reduced the variability in scene to scene sensor view angle in our RapidEye time series.This apparent advantage of the RapidEye data may be due to a combination of factors, including spatial resolution (5 m pixels vs. 30 m pixels) and spectral sensitivity of the sensor.While this may not play a strong role in more homogeneous, closed canopy systems, sensor view angle in this study often imposed more variability on NDVI than natural seasonal variability.Consequently, we recommend that remote sensing efforts to model VI sensitive processes in heterogeneous, low fractional cover systems, place high constraints on acquisition angles for time series, or bin analyses by viewing angle to minimize the potential confounding effects.
We recognize that the temporally resolved RapidEye data set we utilized for this study is not a common commodity and currently carries with it a large cost.Using red-edge data added sensor and illumination geometry complexity, but did improve estimates of GPP during periods of ecosystem stress despite it.Our results suggest high resolution, red-edge employing platforms will potentially be very useful for resolving changes in canopy function during periods of rapid disturbance or recovery where LAI may be changing slowly in relation to chlorophyll content.The recently launched Sentinel-2 satellite mission will allow this to be tested on a broader scale by providing greater spatial and temporal resolution than Landsat, as well as the ability to calculate NDRE and NDWI, and be freely available.Secondly, the upcoming soil moisture active passive sensor (SMAP) may provide either direct measurements, or modeled estimates of soil moisture, providing further predictive power to estimate carbon uptake rates in semi-arid ecosystems.

20 Figure 1 .
Figure 1.Time series of tower-measured gross primary production (A) and precipitation (B) at both the control (grey) and girdled (black) sites.RapidEye acquisition dates are indicated by the rug plot above the top plot.

Figure 1 .
Figure 1.Time series of tower-measured gross primary production (A) and precipitation (B) at both the control (grey) and girdled (black) sites.RapidEye acquisition dates are indicated by the rug plot above the top plot.

Figure 2 .
Figure 2.One to one plots of tower derived (measured) and predicted (modeled) gross primary productivity for the control (C, black) and girdled (G, grey) sites for all years of the study.See Table1for a description of model names.

Figure 3 .
Figure 3.One to one plots of tower derived (measured) and predicted (modeled) gross primary productivity for the control (C, black) and girdled (G, grey) sites for 2009 only.See Table1for a description of model names.

Figure 3 .
Figure 3.One to one plots of tower derived (measured) and predicted (modeled) gross primary productivity for the control (C, black) and girdled (G, grey) sites for 2009 only.See Table1for a description of model names.

Figure 4 .
Figure 4. (A) RapidEye false color composites of the PJ control site (top) for two image dates 3 days apart in 2009, and the corresponding fraction of pixels with NDVI > 0.30 (bottom), with detailed subsets on the far right; (B) Histogram of view angles that comprised the entire RapidEye time series.View angles are represented as absolute values from nadir.

Figure 4 .
Figure 4. (A) RapidEye false color composites of the PJ control site (top) for two image dates 3 days apart in 2009, and the corresponding fraction of pixels with NDVI >0.30 (bottom), with detailed subsets on the far right; (B) Histogram of view angles that comprised the entire RapidEye time series.View angles are represented as absolute values from nadir.

Figure 5 .
Figure 5.One to one plots of tower derived (measured) and predicted (modeled) gross primary productivity for the control (C, black) and girdled (G, grey) sites for RapidEye acquisition angles < 7 degrees off nadir.See Table1for a description of model names.

Figure 5 .
Figure 5.One to one plots of tower derived (measured) and predicted (modeled) gross primary productivity for the control (C, black) and girdled (G, grey) sites for RapidEye acquisition angles <7 degrees off nadir.See Table1for a description of model names.

Table 1 .
Model framework description.Each model is structured as GPP ~ PAR × TA × VI where VI is represented by the corresponding Table entry above, while PAR and TA are derived from on-site instrumentation.Acronym definition: GPP (gross primary productivity), PAR (photosynthetically active radiation), TA (air temperature), VI (vegetation index), NDVI (normalized difference vegetation index), NDRE (normalized difference red-edge index), NDWI (normalized difference

Table 1 .
Model framework description.Each model is structured as GPP ~PAR ˆTA ˆVI where VI is represented by the corresponding Table entry above, while PAR and TA are derived from on-site instrumentation.Acronym definition: GPP (gross primary productivity), PAR (photosynthetically active radiation), TA (air temperature), VI (vegetation index), NDVI (normalized difference vegetation index), NDRE (normalized difference red-edge index), NDWI (normalized difference wetness index).

Table 2 .
Model performance statistics for only year 2009 of the experiment.See Table1for a description of model names.R 2 adj is the adjusted R 2 , ∆AICc is the difference in model AICc (Akaike information criterion) relative to the lowest scoring model, RelLL is the relative log likelihood, or how likely it is that each model is actually the best, and Weights is a list of the AICc weights for each model.

Table 4 .
Model performance statistics for only year 2011 of the experiment, a significant drought year.See Table1for a description of model names.R 2 adj is the adjusted R 2 , ∆AICc is the difference in model AICc relative to the lowest scoring model, RelLL is the relative log likelihood, or how likely it is that each model is actually the best, and Weights is a list of the AICc weights for each model.

Table 5 .
Model performance statistics for RapidEye sensor view angles less than 7 degrees off nadir.See Table1for a description of model names.R 2 adj is the adjusted R 2 , ∆AICc is the difference in model AICc relative to the lowest scoring model, RelLL is the relative log likelihood, or how likely it is that each model is actually the best, and Weights is a list of the AICc weights for each model.