Integrating proximal broad-band vegetation indices and carbon fluxes to model gross primary productivity in a tropical dry forest

The measurement of carbon exchange between vegetation and the atmosphere is vital to quantify the impact of environmental variables on the carbon sequestration capacity of forests, and to predict how they will respond to future climate. In this study we use proximal remote sensing, defined as observations made from non-contact radiometric or imaging sensors in close proximity to the forest canopy (10–20 m), as an intermediate upscaling tool between direct measurements of carbon fluxes and satellite-derived estimations of primary productivity in a tropical dry forest (TDF) in Jalisco, Mexico. Two broad-band vegetation indices (VIs), the normalized difference VI and the enhanced vegetation index 2 (EVI2), were calculated from proximally sensed canopy properties, validated with field estimates of the fraction of absorbed photosynthetically active radiation by photosynthetic tissue (fAPARgreen), and compared to estimates of gross primary productivity (GPP) and net ecosystem exchange of CO2, measured from a flux tower. The VIs captured the phenology of the TDF, both under typical summer rainfall and during an atypically-dry wet season in El Niño of 2009. The VIs also tracked a secondary leaf-flushing in the dry season of 2010. Our study suggests that (1) VIs are the best predictors of gross carbon uptake, able to explain up to 86% of variations in GPP; (2) VIs are accurate predictors of the photosynthetic capacity of green tissue, able to explain up to 99% of fAPARgreen variation; and (3) VIs and soil water content can be used to develop an empirical model that captures the seasonal trajectory of GPP from high respiration after the rain pulses, to rapid leaf development, and finally to slow senescence as the soil dries out. Proximal remote sensing constitutes a useful tool to link field-base measurements of carbon fluxes to satellite- or airborne-derived estimates of carbon exchange.


Introduction
The balance between carbon sequestration in forest biomass and carbon release through respiration, i.e. the net primary production of forests, is both a provisioning ecosystem service that drives the supply of food and wood for human consumption, and a regulatory ecosystem service as it reduces the amount of CO 2 in the atmosphere, and thus mitigates climate change. The quantification of the temporal and spatial variations of carbon stocks and carbon exchange between forests and the atmosphere is therefore essential to our understanding of the global carbon cycle, and to our ability to predict future climate. At the national and regional levels, carbon monitoring is usually based on a combination of field-based forest inventories and remote sensing of structural attributes of forests, such as aboveground biomass, canopy cover or changes in amount of forested land (e.g. Urbazaev et al 2016). This approach, combined with process or empirical models, may be sufficient to develop nationwide greenhouse gas inventories or carbon stock estimations. However, research-intensive sites that combine field plots, flux towers and remote sensing techniques can be deployed to produce detailed information on carbon dynamics at local and regional levels (Birdsey et al 2014). Such information may be used to generate and test physiological parameters that inform model parameterizations, or to generate emission factors, which in turn may be applied to assess carbon gains and losses in similar forests climate conditions, where implementing a detailed carbon monitoring program is unfeasible.
The majority of flux monitoring sites are located in temperate and boreal ecosystems, with a smaller number of sites situated in evergreen forests of the humid tropics and in tropical savannas (ORNL DAAC 2013). At the global scale, recent work indicates that the drought-sensitive vegetation of mesic biomes (e.g. tropical and temperate rainforests) determines the strength of the land carbon sink. However, there is increasing evidence that drought-adapted vegetation of low latitudes has a dominant effect on both the year-to-year variations of carbon exchange and land greening trends (Ahlström et al 2015, Haverd et al 2017, Poulter et al 2014, and that the response of tropical biomes to water availability underlies their sensitivity to temperature and may control the variation of the terrestrial carbon balance (Wang et al 2014, Jung et al 2017. Among ecosystems where drought-adapted vegetation dominates, tropical dry forests (TDFs) have received very little attention from the flux community and their contribution to the global carbon balance is still poorly understood. As of today only five of 517 FLUXNET sites are reportedly situated in this type of vegetation (ORNL DAAC 2013b). This contrasts directly to the importance, extent, and vulnerability of TDFs to disturbance (Portillo-Quintero and Sanchez-Azofeifa 2010, Portillo et al 2015, Calvo-Rodriguez et al 2017. In the Americas alone, their current extent is approximately 500 000 km 2 , but compared to their potential extent, 66% of these forests have already been converted to other land uses (Portillo-Quintero and Sanchez-Azofeifa 2010). Moreover, estimates of aboveground biomass and CO 2 emissions to the atmosphere during conversion to other land uses suggest these forests might be sites of substantial CO 2 exchange (Kauffman et al 2003, Martínez-Yrizar et al 1996. At larger spatial scales, neither the magnitude of primary production of TDFs nor their sensitivity to environmental drivers is adequately characterized via remote sensing. The widely used MODIS product for the gross and net primary production (GPP/NPP) of Earth's vegetated land surface (Friedl et al 2010, Running et al 2004 relies on the light use efficiency model first proposed by Monteith (1972). In the MODIS algorithm, a biome-specific parameter is extracted from the Biome Property Lookup Table for different vegetation types derived from five land cover classification schemes (Friedl et al 2010), none of which recognize TDF as a distinct vegetation type, since areas classified by Miles et al (2006) as TDF appear either as broadleaf deciduous forests or as savannas in these schemes, and therefore the ecophysiological particularities of tree drought-deciduousness in the tropics are disregarded.
There is growing interest in using carbon flux data to validate and refine the parameters used in satellitebased estimates of primary production (e.g. Canadell et al 2000, Heinsch et al 2006, Yuan et al 2007. A key challenge is the inherent mismatch of the spatial and temporal scales of optical and flux measurements. Fluxes are sampled at high temporal frequency but the size of the source area for a given eddy covariance (EC) system typically ranges from hundreds of square meters to 1-2 km 2 , and varies with atmospheric stability and wind speed (Göckede et al 2004, Schmid 1994. Meanwhile, satellite-borne sensors view a much larger scene at any given time but have limited spatial and temporal resolution.
Recently, many flux sites have deployed optical sampling devices at ground level, above the canopy, as a way to overcome the mismatch between Earth observation platforms and flux measurements (Huemmrich et al 1999, Eklundh et al 2011, Gamon 2015. This proximal remote sensing of the optical properties of the canopy (also referred to as in situ, ground-based, or near-surface remote sensing) provides a closer spatial and temporal match to fluxes while remaining directly scalable to optical satellite products. It is also less affected by atmospheric conditions and the cloud contamination that pervades satellite and aerial observations in tropical ecosystems (Sanchez-Azofeifa et al 2017). Huemmrich et al (1999) first demonstrated the use of above-canopy reflectances in the photosynthetically active radiation (PAR) and near-infrared regions to derive a daily broad-band normalized difference vegetation index (NDVI) to track phenology in boreal sites. Since then, proximal broad-band indices have been used to monitor the seasonality of vegetation in a wide range of ecosystems (e.g. Soudani et al 2012, Wilson andMeyers 2007), compared to satellite-derived estimates of productivity (Fensholt et al 2004, Tittebrand et al 2009, Wilson and Meyers 2007, and used as proxy data to fill missing periods in flux records (Nestola et al 2016, Wohlfahrt et al 2010. As of today, no similar studies have been implemented for a TDF. In this study we concurrently measured the carbon exchange of a pristine TDF and the phenological changes in its photosynthetic capacity using vegetation indices (VIs) derived from proximal remote sensing.
We aimed to answer the following questions: (1) can broad-band VIs be used to track the dynamics of carbon fluxes in a TDF over seasonal to annual time scales?, (2) how well do broad-band VIs correlate with a groundbased physiological parameter like f APARgreen ?, and (3) can broad-band VIs be used as predictors of carbon fluxes and aid in the gap-filling process of flux time series, and ultimately in the parameterization of models that predict the GPP in similar types of forests and climate conditions?

Site description and instrumentation
A micrometeorological tower equipped with an EC system for flux monitoring, and proximal remote sensing instruments were installed at Estación de Biología Chamela (EBCh) (19.5095 • N, 105.0401 • W) of Universidad Nacional Autónoma de México, 3 km inland from the Pacific coastline, in Jalisco, Mexico (figure 1), and operated between 2007 and 2015. The climate at EBCh is warm and sub-humid, with summer rains and an isothermal regime (Garcia 1988). Annual temperature averages at 25.8 • C . Mean annual precipitation is 844 mm , 89% of which falls from May to October. The tower is surrounded by low convex hills covered by old-growth, highly diverse TDF (Bullock and Solís-Magallanes 1990). Canopy height at the tower location is ∼10.5 m, with a dense understory.
The EC system was installed on a steel tower, 12.5 m above ground level. The system consisted of a 3D sonic anemometer and an open-path CO 2 /H 2 O infrared gas analyzer (IRGA) (table 1). Wind velocity components, sonic temperature and gas concentrations were measured at a frequency of 10 Hz. We also recorded 30 min averages of ancillary meteorological variables (table 1). Daily precipitation was measured with a standard rain gauge at the meteorological station of EBCh, 1.5 km away from the tower site.
A detailed explanation of the methods used to estimate carbon fluxes can be found in Gonzalez del Castillo et al (2018). Additional information on flux partitioning, gap filling, proximal broad-band indices, fractions of absorbed PAR and GPP modeling can be found as supplementary material available at stacks.iop.org/ERL/13/065017/mmedia.

Environmental drivers of the temporal patterns of carbon flux
In spite of the high variability of carbon exchange in this TDF, a consistent temporal correlation between rainfall, greenness, and carbon fluxes can be discerned from the time series in figures 2 and 3. Precipitation during the rainy season of 2008 (May-October) was 35% above the 30 yr average (1981-2010, 805 mm) for our site. Due to a moderate warm phase of El Niño-Southern Oscillation, which prevailed from June 2009-May 2010, cumulative rainfall for the rainy season of 2009 was 23% below the long-term mean (720 mm), but 3.2 times the average (86 mm) for the dry season (November 2009-April 2010). Accordingly, the midday average of the net ecosystem exchange of CO 2 (NEE md ) and the daily NEE were respectively 16% and 8% lower during the summer of 2009 compared to those at the beginning of the rainy season in 2008, and 44% and 30% lower than the respective rates during the summer of 2010. The rains in February 2010 triggered a flush of new leaves in the middle of the dry season that resulted in sustained C gains until the end of April. The average daily NEE during this brief growing period (−2.4 ± 1.2 g C m −2 d −1 ) was comparable to the average daily NEE for the , daily gross primary productivity (GPP, dark blue) and ecosystem respiration (R eco , red), daily average of air temperature (orange) and vapor pressure deficit (light blue), soil temperature (green) and soil water content (brown), and total daily precipitation in summer and fall (dark blue) and winter and spring (light blue). Bars indicate the daily range of each variable.  Verduzco et al (2015) found that the respiration-dominated flux accounted for 36%-148% of the total C fixed during the wet season in a TDF, and determined whether the system acted as a source or a sink of C in a given year. Several non-exclusive mechanisms have been proposed for these transient CO 2 outbursts, among them physical displacement of air from soil pores by water (Unger et al 2010), and rapid decomposition and mineralization of organic C pools by soil heterotrophs (Tang and Baldocchi 2005, Unger et al 2010, Xu et al 2004. EC data cannot be used to unequivocally ascertain the origin of the CO 2 efflux, but the signal of an immediate physical expulsion of air from soil is likely missed because such flux would be small and transient (Unger et al 2010), and data gathered during the rain event itself and until the open-path IRGA dries out are discarded. NEE partition reveals that after the initial CO 2 outburst, GPP grew at a slower pace, and its peak lagged behind R eco by 8 days on average in summer 2008, summer 2009 and winter 2010 (figure 2). This suggests that R eco was at first dominated by the activation of heterotrophic respiration of labile C substrates accumulated on the forest floor during the long dry season. Indeed in this same TDF, Anaya et al (2012) found that rain events >10 mm were positively correlated to the decomposition rates of litterfall. Tissue construction and transport of recent photosynthate to the soil take longer and may not start until deeper roots have access to water (Carbone et al 2011, Huxman et al 2004, so the contribution from autotrophic respiration may increase over the days following a rain event. CO 2 uptake also lags behind R eco since it depends on leaf expansion and maturation. We found that the increase in GPP was closely matched by the rise in the proximal greenness indices (EVI2 md and NDVI md ) from their dry season minima to their peak in ∼13 days on average, as leaves expanded and the canopy closed ( figure 3).
The amplitude of the annual thermal oscillation was small; average daily temperatures dropped to 20 ± 3 • C from November through mid-April and went up to 27 ± 3 • C from June to late October (figure 2). Half-hourly minimum temperatures during the study (10 • C-12 • C) were registered in early March 2008 due to the passage of a cold front, while maximum temperatures of 33 • C-35 • C were recorded for several days from July to September 2009. Due to the constant supply of water vapor from the nearby ocean, relative humidity was high throughout the year; only 12% of half-hourly values fell below 0.60. As a consequence, the vapor pressure deficit (VPD) generally remained low (0.2-1.7 kPa) (figure 2).
Every year soil temperature rose as the dry season advanced, reaching daily maxima of 35 • C-39 • C right before the onset of the rains, during periods when almost negligible NEE was observed. Then, in June 2008, June 2009 and late July 2010, the average T soil cooled overnight by 3 • C-5 • C after the first substantial rain fell (figure 2). The daily T soil range also diminished from an average of 4 • C prior to the rains to less than 1.5 • C during the rainy season. After the drop, the average T soil slowly diminished again toward its annual winter minimum. The negative correlation of T soil with the volumetric soil water content (SWC) and its limited variability made it an unsuitable predictor of nighttime R eco , as reported for other tropical locations (Davidson et al 2000, Hutyra et al 2007; thus, SWC was preferred for NEE partitioning. The seemingly small effect of T soil on C fluxes was nevertheless incorporated into R eco modeling since a separate equation for the relation SWC-R eco was fitted for different T soil bins (table 2). After the rains in the dry season of 2010, SWC reached the same level of wet season rain pulses, but lower estimates for R eco at the lowest T soil class resulted in lower daily R eco rates, since it coincided with the annual minimum T soil . During this period, canopy greenness peaked at an EVI2 md of only 0.45, which is 30% lower than the average summer-fall maximum (0.65) (figure 3), but net C fixation rates were similar to those of the drier-thanaverage summer and fall of 2009, because the smaller respiratory losses at cooler T soil were soon exceeded by the daily GPP. As could be expected from the described asynchrony between these factors and the C exchange at Chamela, T a , VPD, PAR and T soil exhibit low (some non-significant) correlations with C flux (table 3). In contrast, SWC is able to explain 79% and 58% of R eco and GPP variation, respectively; the strong correlation between soil humidity and C fluxes is expected since they were partially modeled based on SWC. Less predictably, SWC alone is also able to explain almost 40% of NEE md variability.

Proximal VIs as predictors of carbon fluxes
Proximal, broad-band VIs were able to capture the vegetative phenophases of the deciduous canopy at the Chamela Biological Station (figure 3). Both NDVI md and EVI2 md rose from their dry-season minima (around 0.44 and 0.25, respectively) to their peaks (0.78 and 0.65) during quick green-up periods of 14, 11 and 14 days after the first rains in June 2008, June 2009, and July 2010, respectively. In the wet season of 2008, the fully developed, mature canopy, maintained a level of greenness close to peak values until late October. Leaves then senesced and fell, and the VIs gradually descended and reached their annual minima in April 2009. After the dry-season rains in February 2010, NDVI md and EVI2 md increased from higher minimal values (0.51 and 0.27) since the canopy was not completely leafless, and reached lower maxima of 0.68 and 0.46, respectively. Both indices captured interannual variation in canopy development; the peak NDVI md and EVI2 md values for the drier wet season of 2009 were lower than those in 2008 and 2010, but EVI2 md showed more intra-seasonal variability at the fully developed canopy stage than NDVI. At high values, EVI2 md also exhibited more day-to-day noise than NDVI md . Proximal EVI2 md and NDVI md were the most important predictors of C exchange, statistically explaining, respectively, 77% and 71% of NEE md, and 86% and 81% of daily GPP, but only 59% and 53% of R eco variations (table 3, figure 4). Over seasonal and annual time scales, a strong influence of the photosynthetic capacity of the canopy on R eco is to be expected, as root respiration depends on the allocation of photosynthates recently produced by the aerial parts of the plant (Högberg et al 2005, Janssens et al 2001). However, controls exerted on the heterotrophic component of R eco may dominate over shorter periods, as discussed in section 3.1. An inspection of the high R eco data points that depart from the bulk of the dispersion at EVI2 md <0.45 and NDVI md <0.65 reveals that they correspond to periods of steep increase in VIs after the first rains, before peak photosynthetic capacity is reached. It is likely that in these conditions high R eco is produced by elevated metabolic rates during the active construction of new plant tissues combined with large heterotrophic CO 2 efflux driven by high SWC. Not surprising either is the low correlation of daily NEE to VIs, since NEE carries the signal of two opposing processes (GPP and R eco ) that summed over diurnal cycles are of almost equal magnitude (Reichstein et al 2005).
We also explored the possibility of multiplicative effects of environmental variables over the single best predictor (EVI2 md ) of GPP; none of the direct products of EVI2 md with PAR, T soil , or SWC improved the obtained correlation to EVI2 md alone (table 3).
The advantage of using EVI2 over the most widely used NDVI is illustrated in figure 4. Some saturation of EVI2 md at high levels of GPP was evident; variations in the photosynthetic capacity of the forest at the peak of the growing season were not matched by variations in the optical properties of the canopy. However, the relationships of the main C fluxes to EVI md were mostly linear. In contrast, NEE md , R eco and GPP exhibit good sensitivity to changes in NDVI md < 0.7, but above this value the VI is almost invariant while the dynamic range of C fluxes may span 6-8 g m −2 d −1 for GPP and R eco respectively, and up to 17 mol m −2 s −1 for NEE md .

Ground-based measurements of radiation absorption
Field estimates of f APAR and f APARgreen identified senescence in late 2009 and winter growth, green-up, and maturity peak phenophases in the wet season of 2010 ( figure 5). However, f APARgreen dropped more significantly in the fall and winter of 2009, better mirroring the behavior of the VIs.
There were large temporal variations in the spatial repeatability of f APAR and f APARgreen ; the SDs of both estimates decreased five to six times from the dry to the wet season ( figure 5). This points to a mismatch between the footprint of the sensors at the tower and the ceptometer on the ground, which    is difficult to circumvent. Down-looking instruments above the canopy sense a mixture of the radiation reflected by numerous individual elements of the vegetation and the ground within their field of view. At ground level this mixture is sampled at discrete locations, so a topographical feature, a large tree, or a bare patch of soil may disproportionately influence the small area seen by the ceptometer at any given sampling point. The influence of such point-to-point variation is larger during the dry season and under thinner canopies compared to the uniform shade provided by full foliage. Weighted regressions ( figure 5, table 4) show that there is an apparent non-linearity of the EVI2 mdf APAR relationship that disappears for f APARgreen , which confirms the greater ability of f APARgreen to capture the canopy photosynthetic capacity, without the influence of PAR absorption by non-photosynthetic tissues Gamon 2015, Nestola et al 2016). Narrow-band NDVI and EVI2 exhibit an asymptotic relationship with f APARgreen in crops, becoming insensitive to f APARgreen changes for NDVI > 0.7 and EVI2 > 0.6, respectively (Goward andHuemmrich 1992, Viña andGitelson 2005). The dynamic range of the broad-band VIs at Chamela fell almost entirely within the linear portion of the VI-f APARgreen relationship, so both indices were accurate predictors of f APARgreen (figure 5).

Potential of proximal VIs for GPP modeling
We found moderate correlations between the parameters of the light response curve and EVI md (figure 6). The relation to the apparent quantum yield ( ) was linear, but in contrast to Wohlfahrt et al (2010), both the photosynthetic rate at the saturating PAR in ( ) and daytime R eco (R eco−dt ) showed an asymptotic response to EVI md (table 5). Correlation coefficients (−0.74, −0.85 and 0.71 for , and R eco−dt , respectively) are within the range found by Wohlfahrt et al (2010) for two mountain grassland sites.
From June 2009 to February 2010, absorbed radiation (APAR) derived from EVI2 md was able to explain less GPP variability than EVI2 md directly (57% and 72%, respectively), although the relation APAR-GPP was more linear (figure 7). The addition of the second best predictor (SWC) reveals a U-shaped trajectory response: starting from the right, upper, back side of the cube in figure 7, the cubical points (dark green) that depart clearly from the rest of the data correspond to the initial CO 2 released immediately after the first rains but before canopy closure. The trajectory then depicts an increase toward higher photosynthetic activity (high APAR, lower GPP) as the canopy developed rapidly (dark to medium green dots) under high SWC. Later, as soil dried and photosynthetic capacity diminished, C fixation was reduced during the senescence phase (light green to orange to pink dots).    Orange cubes represent respiration-dominated days when rain pulses interrupted intra-seasonal dry spells. This trajectory was consistent among all seasons (not shown). When the relationship of GPP with APAR derived from EVI2 md is examined with the exclusion of SWC, the cubical points are not distinguishable from the rest of the points in the scatterplot; the fact that they clearly stand out when the SWC is considered, reveals that different mechanisms contribute to low uptake at the beginning of the rainy season and at the end of the senescing phase, and adds complexity to our understanding of C dynamics' dependence on water availability. The best-fit surface for this relation is a Gaussian model with slightly higher explanatory power (77%) than EVI2 md alone (table 6). This model, although purely statistical and therefore site-specific, produced GPP estimates that closely compared to GPP series gap-filled using look-up tables of observed fluxes (figure 8, left), while the procedures that modeled the light response curve based on EVI2 md (figure 8, middle) or directly based on PAR (figure 8, right) tended to overestimate high GPP and underestimate low GPP (figure 8, table 7). Because large data gaps remained in the record (Gonzalez del Castillo, unpublished results), we did not attempt to obtain seasonal or annual sums of C fluxes, but if the tendencies depicted in figure 8 prevailed for the whole data sets, models relying on the light response curve would likely overestimate C gained during the summer but would underestimate late-fall fixation rates.

Conclusions
Our study provided two important conclusions for TDFs in this region. The first conclusion is mechanistic in nature. As expected for these ecosystems-but not fully documented or linked to carbon fluxesthere is a clear controlling influence of precipitation and soil moisture on carbon release and uptake. Water availability drives both the respiratory release of Table 7. Linear regression coefficients and statistics for the relationship between daily GPP that was gap-filled using different procedures and daily GPP gap-filled using standard look-up tables of half-hourly PAR and SWC. carbon from soils at the onset of the rains and the photosynthetic capacity of the canopy throughout the growing period. The response of the TDF to the atypical rains at the peak of the dry season is a clear example of this dependence. Our second conclusion is related to the role that proximal remote sensing plays as a tool to estimate carbon fluxes in a TDF environment. The integration of proximal remote sensing indices with the flux and meteorological data showed two main advantages: (1) it yielded a more complex picture of the way water availability controls carbon flux in this ecosystem, and (2) it allowed the prediction of GPP values comparable to those obtained using past relations between meteorological drivers and observed fluxes, and therefore may be used as a valuable tool in GPP gap-filling and modeling. It is fundamentally clear that, once spectral vegetation indices have been properly validated against field estimates of light absorption and conversion to biomass in a particular forest, they show good promise as aids in thorough carbon monitoring programs at the local scale.