Exploring the Potential of Satellite Solar-Induced Fluorescence to Constrain Global Transpiration Estimates

1 Laboratory of Hydrology and Water Management, Ghent University, Coupure Links 653, 9000 Ghent, Belgium; Wh.Maes@ugent.be(W.M.); Brecht.Martens@UGent.be (B.M.); Diego.Miralles@UGent.be(D.G.M.) 2 Department of Earth and Environmental Engineering, Columbia University, 842D S.W. Mudd, 500 West 120th St., New York, NY 10027, USA; pg2328@columbia.edu * Correspondence: Brianna.Pagan@ugent.be; Tel.: +32-9-264-6140 Received: 31 December 2018; Accepted: 9 February 2019; Published: 16 February 2019


Introduction
Plant transpiration dominates the global flux of terrestrial evaporation and is a central element of the hydrologic cycle over land [1,2]. The transpiration flux depends on plant phenology and is mediated by stomatal conductance, by which plants tend to optimize to maximize carbon uptake while reducing water losses [3][4][5][6]. As such, the biophysical and biochemical feedbacks on climate depend on vegetation phenology and plant physiological response to environmental conditions and nutrient cycling [7][8][9][10]. These responses vary among plant species and make the modelling of transpiration particularly challenging [11][12][13]. More specifically, the physiological response to water stress remains a longstanding modelling challenge [14], which may ultimately be addressed by incorporating plant hydraulics into land surface models (LSMs) [15][16][17].
T ∝ SIF × WUE (1) where WUE reflects the ratio of GPP to T. While in nature WUE depends on a plethora of factors, it is widely acknowledged that the atmospheric demand for water explains a large fraction of its variability at ecosystem scale [44,61]: under conditions of high potential evaporation (E p ), the water loss per unit gain of carbon increases. As such, here we explore whether the fraction of E p that is transpired (hereafter referred to as 'transpiration efficiency', τ) holds a proportional relationship to the fraction of photosynthetically active radiation (PAR) that is emitted as SIF: The normalization of SIF (mW/m 2 /sr/nm) by PAR (W/m 2 ) allows the correction of incoming radiation changes-independent of phenological variability or physiological stress [62]-in an analogous way that the normalization of T by E p buffers the effect of short-term changes in available energy. Moreover, normalizing by PAR intentionally retains the variability in SIF associated with changes in the fraction of absorbed PAR (i.e., fPAR), which allows us to employ an analogous E p calculation that is merely atmospheric, so that τ reflects the limitations of the actual ecosystem in response to the atmospheric demand for water (whether that ecosystem limitation has a structural or physiological origin). Following the findings in Maes et al. [63], in which multiple E p equations were evaluated, E p (W/m 2 ) is calculated here based on the simple radiation-based formulation by Milly and Dunne [64]: where α is a constant unitless parameter (similar to α in the Priestley and Taylor equation), R n (W/m 2 ) is the surface net radiation and G (W/m 2 ) is the ground heat flux. Maes et al. [63] recommended using a biome-specific parameterization of α rather than the universal value of α = 0.8 proposed by Milly and Dunne [64]-see Table S1 for specific values of α per biome type. The rationale behind the hypothesis confined in Equation (2) is illustrated in Figure 1. Given a certain amount of R n and PAR, under low stress conditions (and/or high leaf area), we expect higher levels of transpiration (i.e., λE), GPP and SIF than under high stress conditions (and/or low leaf area) given the same R n and PAR. For the latter, we expect more sensible heat flux (H), non-radiative decay and non-photochemical quenching (Q). This hypothesis is further sustained by the results from Maes et al. [65], which, by employing an advanced radiative transfer model, showed that most of the variance in T that is left unexplained by SIF variability can be explained by changes in the atmospheric potential to evaporate water.

Satellite Observations
SIF observations are taken from the Global Ozone Monitoring Experiment-2 (GOME-2) on-board EUMETSAT's Meteorological Operational Satellite-A (MetOp-A) for the time period 2007-2014. These retrievals are based on the method by Köhler et al. [47], and are obtained from the band 4 using a spectral window of 720-758 nm. Data is available globally at 0.5 • spatial resolution and up to 1.5-day temporal resolution in the absence of clouds. SIF data is taken from GOME-2 rather than newer sensors, such as the Orbiting Carbon Observatory-2 (OCO-2) or the TROPOspheric Monitoring Instrument (TROPOMI), as GOME-2 has global coverage and the longest record length. PAR data from the Clouds and Earth's Radiant Energy System (CERES) [66] are selected based on the closest overpass time to GOME-2 (near 9:30AM), and resampled from 1 • to 0.5 • . SIF and PAR retrievals are cross-masked and temporally gap-filled prior to the calculation of their ratio using a 15-day moving window; when Remote Sens. 2019, 11, 413 4 of 15 necessary, estimates are subsequently spatially gap-filled by weighting the eight neighboring pixels based on distance. Time steps with PAR < 40 W/m 2 are masked to exclude potential nighttime conditions in high latitudes. Finally, GLOBSNOW snow water equivalent data [67] is used to mask for snow conditions.
(W/m 2 ) is the surface net radiation and G (W/m 2 ) is the ground heat flux. Maes et al. [63] recommended using a biome-specific parameterization of α rather than the universal value of α = 0.8 proposed by Milly and Dunne [64]-see Table S1 for specific values of α per biome type.
The rationale behind the hypothesis confined in Equation (2) is illustrated in Figure 1. Given a certain amount of Rn and PAR, under low stress conditions (and/or high leaf area), we expect higher levels of transpiration (i.e. λE), GPP and SIF than under high stress conditions (and/or low leaf area) given the same Rn and PAR. For the latter, we expect more sensible heat flux (H), non-radiative decay and non-photochemical quenching (Q). This hypothesis is further sustained by the results from Maes et al. [65], which, by employing an advanced radiative transfer model, showed that most of the variance in T that is left unexplained by SIF variability can be explained by changes in the atmospheric potential to evaporate water.

Tower Data
Data from the FLUXNET2015 archive of eddy-covariance tower measurements (http://fluxnet. fluxdata.org/data/fluxnet2015-dataset/) are used to create an in situ τ validation dataset by computing the ratio between transpiration and potential evaporation at each daily time step for each tower using daytime values. The measured latent heat at each tower is partitioned into transpiration and soil evaporation using the simple approach introduced by Wei et al. [1], where the ratio of T to total evaporation is a function of land cover type and leaf area index (LAI)-the latter is obtained from the Global Inventory Modeling and Mapping Studies (GIMMS) [68,69]. Using the empirical equations provided by Wei et al. [1], International Geosphere-Biosphere Programme (IGBP) classifications are regrouped into six categories, and the ratio of T to total evaporation is calculated for each time step using empirical LAI regression functions for each land cover type. This ratio is then multiplied by the latent heat at the tower to obtain the flux of transpiration. The E p from the towers is calculated using the simple available energy formulation by Milly and Dunne [64] using tower R n and G, with an α value calibrated per individual tower using a subset of unstressed days, as explained by Maes et al. [63]. A series of quality checks are applied to eliminate towers in which representativeness issues may arise (i.e., lack of observations, large fraction of open water, etc.); a full description of these checks is available in the supplemental information. After all quality checks, a total of 25 FLUXNET towers remain for comparisons, as illustrated in Figure S1 and listed in Table S2.

Land Surface Models (LSMs)
To understand the potential of SIF data to benchmark and constrain transpiration in state-of-the-art LSMs and hydrological models, the τ resulting from a range of offline model runs is also compared against in situ τ estimates. These runs come from the eartH2Observe Tier-2 database, which provides water resource reanalysis datasets at a global scale [57]. Model data comes at 0.25 • spatial resolution forced by the WATCH Forcing Data methodology applied to ERA-Interim reanalysis (WFDEI) and the Multi-Source Weighted-Ensemble Precipitation (MSWEP) [70]. Again, τ for each model is approximated by the ratio of T to E p (Equation (2)). All eartH2Observe models provide offline estimates of T, and while some provide estimates of E p , inconsistencies in the definition of E p across models makes their inter-comparison difficult. For a fair comparison between models and in situ measurements, we calculate potential evaporation using the same approach as for the towers.
In total, six models from the eartH2Observe database were used: the Global Land Evaporation Amsterdam Model (GLEAM) v3a [26,71], the WorldWide Water Resources Assessment (W3RA) based on the hydrology component of the Australian Water Resource Assessment Landscape model (AWRA-L) version 1 [72,73], the Hydrology Tiled European Centre for Medium-Range Weather Forecasts (ECMWF) Scheme for Surface Exchanges over Land-Catchment-based Macro-scale Floodplain model (HTESSEL-CaMa) [74,75], the Joint United Kingdom Land Environment Simulator (JULES) [76], the Organising Carbon and Hydrology in Dynamic Ecosystems (ORCHIDEE) [77], and the surface externalisée with TRIP river routing (SURFEX-TRIP) [78,79]. The latter four are LSMs, W3RA is a simplified LSM that is often considered to be a hydrologic model, and GLEAM is a remote sensing-based evaporation retrieval model. Once again, GLOBSNOW snow water equivalent data are used to mask for snow, and MSWEP is used to mask out rain days to avoid the effects of interception.

Validation against In Situ Data
Average Pearson correlation coefficients (R) [80] are calculated based on biweekly aggregates of τ for the period 2007-2014 in comparison to in situ data. Overall, SIF/PAR exhibits an average R = 0.59 for the selection of flux towers, showing individual correlations ranging from R = 0.02 (US-Me2) to R = 0.92 (US-MMS), being statistically significant (p < 0.01) at 22 of the 25 sites. In comparison to individual models, SURFEX (R = 0.56) shows the highest correlations against in situ data; this performance is followed by W3RA (R = 0.50). The remaining models display lower average correlations, ranging from R = 0.09-0.34. The consistently lower performance of some models is further evaluated by analyzing the time series of T and E p separately against the towers ( Figure S2). Results demonstrate that E p estimates from all models are similar, as expected from their common atmospheric forcing, and that they all compare well to the E p retrieved at the towers. As such, most of the model-to-tower disagreements in τ result from uncertainties in T: models with higher correlations of T also display a better performance for τ. Figure 2a represents the results of the site validations separately for each individual tower. The overall higher correlations of SIF/PAR against the tower retrievals of τ is evident: for 17 of the 25 sites, correlations are greater than the model ensemble average. Consistently, when correlations are grouped by IGBP biome type at the flux tower scale, the highest correlations between SIF/PAR and the in situ τ data are found at densely forested sites ( Figure S3)-some exceptions are further explored in Section 4, including the US-Me2 evergreen needleleaf forest site.
Correlations against in situ data are also explored for months of vegetation growth and decay separately; these periods are defined based on the GPP climatology at each tower. Periods where GPP increases from one month to the next are designated as growing, and those where GPP decreases, designated as decaying. Overall, both SIF/PAR and the models' τ capture the dynamics at the towers more realistically during the growing season ( Figure 2b). Nonetheless, rankings of model performance across both phases are similar to the results over the entire timeframe shown in Figure 2a. SIF/PAR exhibits the highest average R during the growing period (0.64), and compares better to the in situ τ than any of the models. However, during the decaying period, the correlations of SIF/PAR and SURFEX against in situ are comparable (R = 0.53 and R = 0.54, respectively). The remaining models average correlations between 0.30-0.55 during growth and 0.39-0.49 during decay periods.
The ability to capture the timing of transpiration throughout the year is critical for LSMs to be able to simulate the effects of land conditions on temperature, precipitation and radiation adequately [22]. This skill is explored in Figure 2c

Spatial and Temporal Variability in Solar-Induced Fluorescence/Photosynthetically Active Radiation (SIF/PAR)
Global seasonal maps of SIF/PAR follow the expected patterns of general vegetation activity in response to different phenological and physiological constraints (Figure 4a). In the Northern Hemisphere, SIF/PAR is lowest during winter (December, January February, (DJF)), suggesting a low transpiration efficiency. In wintertime LAI is lower, with sub-optimal radiation and temperatures, resulting in a reduction of stomatal conductance and carboxylation rates. The highest values (suggesting high transpiration efficiency) occur during the summer (June, July, August (JJA)), when LAI is expected to be high and the vast majority of plant types are active-granted an absence of extreme heat or water stress. In Amazonia, despite the low seasonal variability, SIF/PAR peaks in September, October November (SON). This coincides with the end of the dry season when light availability through the canopy increases, while water is typically not limiting due to the deep rooting structure of emergent species [81,82]. Across the Sahel, SIF/PAR is lowest during DJF and highest during SON, mimicking the timing of the monsoon season. Likewise, in northern Australia, SIF/PAR is the lowest during March, April, May (MAM) and the highest during SON, in phase with the monsoonal rainfall. These spatial and seasonal patterns agree well with the exploration of SIF/PAR by Madani et al. [62]. Finally, Figure 4b displays the seasonal average τ for 2007-2014 for the ensemble of models. While the model ensemble shows markedly lower estimates in water-limited regions, overall spatial patterns agree with those of SIF/PAR, with these similarities being larger for SURFEX and ORCHIDEE-see Figure S4.

Spatial and Temporal Variability in Solar-Induced Fluorescence/Photosynthetically Active Radiation (SIF/PAR)
Global seasonal maps of SIF/PAR follow the expected patterns of general vegetation activity in response to different phenological and physiological constraints (Figure 4a). In the Northern Hemisphere, SIF/PAR is lowest during winter (December, January February, (DJF)), suggesting a low transpiration efficiency. In wintertime LAI is lower, with sub-optimal radiation and temperatures, resulting in a reduction of stomatal conductance and carboxylation rates. The highest values (suggesting high transpiration efficiency) occur during the summer (June, July, August (JJA)), when LAI is expected to be high and the vast majority of plant types are active-granted an absence of extreme heat or water stress. In Amazonia, despite the low seasonal variability, SIF/PAR peaks in September, October November (SON). This coincides with the end of the dry season when light availability through the canopy increases, while water is typically not limiting due to the deep rooting structure of emergent species [81,82]. Across the Sahel, SIF/PAR is lowest during DJF and highest during SON, mimicking the timing of the monsoon season. Likewise, in northern Australia, SIF/PAR is the lowest during March, April, May (MAM) and the highest during SON, in phase with the monsoonal rainfall. These spatial and seasonal patterns agree well with the exploration of SIF/PAR by Madani et al. [62]. Finally, Figure 4b displays the seasonal average τ for 2007-2014 for the ensemble of models. While the model ensemble shows markedly lower estimates in water-limited regions, overall spatial patterns agree with those of SIF/PAR, with these similarities being larger for SURFEX and ORCHIDEE-see Figure S4.

Discussion
Overall, the results of the validation of SIF/PAR against in situ eddy-covariance estimates reveal the potential of SIF/PAR as a diagnostic of transpiration efficiency. This transpiration efficiency is understood here as the ability of the land cover to meet the atmospheric demand for water, thus may relate to purely phenological variability (e.g. LAI dynamics) as well as to physiological stress related to e.g. extreme climate. SIF/PAR exhibits higher correlations against the in situ data compared to any of the eartH2Observe models, despite its coarser resolution (0.5° versus 0.25°). In a few sites, some models exhibit equivalent or higher performance than SIF/PAR. It is important to note that despite smoothing daily values and using biweekly averages, sensor noise may still be prevalent in SIF/PAR calculations. Nonetheless, the Köhler et al. [47] GOME-2 retrievals tend to exhibit less noise than those from other methods (see ibid.). Future satellites like the European Space Agency's Fluorescence Explorer (FLEX) or the recently-launched Sentinel 5P (mounting the TROPOMI sensor) are expected to improve these insufficiencies.

Discussion
Overall, the results of the validation of SIF/PAR against in situ eddy-covariance estimates reveal the potential of SIF/PAR as a diagnostic of transpiration efficiency. This transpiration efficiency is understood here as the ability of the land cover to meet the atmospheric demand for water, thus may relate to purely phenological variability (e.g., LAI dynamics) as well as to physiological stress related to e.g., extreme climate. SIF/PAR exhibits higher correlations against the in situ data compared to any of the eartH2Observe models, despite its coarser resolution (0.5 • versus 0.25 • ). In a few sites, some models exhibit equivalent or higher performance than SIF/PAR. It is important to note that despite smoothing daily values and using biweekly averages, sensor noise may still be prevalent in SIF/PAR calculations. Nonetheless, the Köhler et al. [47] GOME-2 retrievals tend to exhibit less noise than those from other methods (see ibid.). Future satellites like the European Space Agency's Fluorescence Remote Sens. 2019, 11, 413 9 of 15 Explorer (FLEX) or the recently-launched Sentinel 5P (mounting the TROPOMI sensor) are expected to improve these insufficiencies.
When analyzing the validation results grouped by IGBP biome class, SIF/PAR is particularly correlated with the in situ data over densely forested sites, along with grasslands and croplands. SIF/PAR exhibits marked improvements in correlations (0.67) compared to the model average (0.47) for these particular IGBP classes. In contrast, evergreen broadleaf and needleleaf forest sites, like US-Me2, show lower correlations; Madani et al. [62] found similar results when relating SIF with tower GPP (see Figure S3). Further validations indicate that SIF/PAR captures the timing of transpiration efficiency accurately (Figure 2c), and may thus be used to evaluate drought periods. While other optical vegetation indices can be considered, SIF is the only physically-representative remote sensing variable that is expected to capture short-term dynamics in vegetation activity, thus reflect the early stages of transpiration stress. Consequently, SIF has been shown to exhibit the shortest lags with respect to the in situ estimates of latent heat, when compared to traditional optical indices such as LAI, normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI) (see Maes et al. [44,65]).
A potential limitation of our analysis is the lack of validation towers in the Southern Hemisphere, and the fact that not every IGBP class is represented by the 25 towers used for validation. Therefore, the overwhelming majority of towers used in this analysis are clustered in northern mid to high latitudes. Moreover, these sites experience seasonal snow in winter, corresponding to times of declining GPP. Subsequent precipitation and snow masking limits the number of observations, which may explain why the correlations against the in situ observations are slightly lower during the vegetation decay period (Figure 2c). Furthermore, SIF is not corrected for angular anisotropy, which can negatively impact the quality of the observations [83]. Finally, 'instantaneous' SIF/PAR retrievals are here compared to daily-average τ from models and towers, as required by the daily resolution of the eartH2Observe data. Consequently, the early overpass timing of GOME-2 (9:30AM), may penalize the SIF/PAR results, as vegetation activity is not expected to be highest in the early morning.
Overall, seasonal analyses and validations demonstrate the potential of SIF/PAR to capture the integrated effect of multiple environmental stress factors on transpiration, such as heat and water stress in dry times and regions or sub-optimal radiation and temperature in wet and cold regions. For instance, at high latitudes, water is not typically a limiting factor, yet sub-optimal temperature and radiation levels may reduce photosynthesis and thus transpiration. This is clearly captured by SIF/PAR (Figure 4a). Especially for deciduous forests, transpiration constraints are not obviously captured by most models, with the exception of SURFEX and ORCHIDEE (see e.g., Figure S4). In the Amazon rainforest, the highest transpiration efficiency occurs at the end of the dry season (SON), which is consistent with the expectations from non-drought years as reported by previous studies [62,84,85]. The relatively lower SIF/PAR values over the Amazon may directly be linked to the higher noise levels in the SIF retrievals over South America due to the South Atlantic magnetic anomaly [86]. Nonetheless, especially in high latitude and tropical biomes, non-climatic drivers can explain a large fraction of the variability in vegetation dynamics [87]; i.e., factors that are poorly represented in models, such as nutrient availability and short-term natural and anthropogenic disturbances (e.g., fires, logging, harvesting, insect epidemics) can be crucial in these environments [88][89][90][91]. The effective assimilation of SIF in global transpiration formulations would directly reflect these dynamics without the need to prescribe species-dependent sensitivities to these stress factors, as in the case of process-based models.

Conclusions
We present a novel diagnostic of transpiration efficiency, understood here as the ability of the land surface to meet the atmospheric demand for water vapor. This diagnostic is based on the ratio of SIF over PAR. Validation suggests that SIF satellite observations can potentially be used to constrain transpiration estimates in land surface models. SIF/PAR combines phenological and physiological constraints on the atmospheric demand for water, and adequately reflects the timing of stress. In congruency with the seasonal results, SIF/PAR appears to be capturing more than just water stress and leaf area changes, by demonstrating sensitivity to other factors such as sub-optimal radiation and temperature levels. However, further analysis is needed to determine the capability of SIF/PAR to capture the impact of drought events or heatwaves on transpiration. Estimates of SIF/PAR can be incorporated into transpiration retrieval models (such as GLEAM) or offline land surface models via data assimilation. Finally, our results support the potential of SIF observations to benchmark climate models. This applies not just to the representation of ecosystem carbon exchange, but also to the estimates of surface water and energy fluxes, at a time in which ecosystem stress is still incorporated in climate models via empirical parameterizations.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/4/413/s1, Figure S1: Global map of FLUXNET towers used for comparisons, Figure S2: Taylor diagrams comparing transpiration and potential evaporation from each eartH2Observe model with FLUXNET towers, Figure S3: SIF/PAR correlations against tower τ grouped by IGBP classifications, Figure S4: Average values for τ from each eartH2Observe model globally. Table S1: List of adjusted alpha values implemented per biome type, Table S2: List of FLUXNET sites used for this analysis.